Home > NoiseTools > nt_pca0.m

nt_pca0

PURPOSE ^

[topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca

SYNOPSIS ^

function [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w)

DESCRIPTION ^

[topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca

  topcs: matrix to convert data to PCs
  pwr: power per PC

  x: data matrix
  shifts: array of shifts to apply
  nkeep: number of PCs to keep
  w: weight (see nt_cov)
  threshold: remove components with normalized eigenvalues smaller than threshold (default: 0)

 mean is NOT removed prior to processing

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w)
0002 %[topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
0003 %
0004 %  topcs: matrix to convert data to PCs
0005 %  pwr: power per PC
0006 %
0007 %  x: data matrix
0008 %  shifts: array of shifts to apply
0009 %  nkeep: number of PCs to keep
0010 %  w: weight (see nt_cov)
0011 %  threshold: remove components with normalized eigenvalues smaller than threshold (default: 0)
0012 %
0013 % mean is NOT removed prior to processing
0014 
0015 
0016 if nargin<1; error('!'); end
0017 if nargin<2||isempty(shifts); shifts=[0]; end
0018 if nargin<3; nkeep=[]; end
0019 if nargin<4||isempty(threshold); threshold=0; end
0020 if nargin<5; w=[]; end
0021 
0022 [m,n,o]=size(x);
0023 
0024 % remove mean
0025 %x=fold(demean(unfold(x)),size(x,1));
0026 
0027 % covariance
0028 if isempty(w);
0029     c=nt_cov(x,shifts);
0030 else
0031     c=nt_cov(x,shifts,w);
0032 end
0033 
0034 % PCA matrix
0035 if ~isempty(nkeep)
0036     [topcs,ev]=nt_pcarot(c,nkeep);
0037 else
0038     [topcs,ev]=nt_pcarot(c);
0039 end
0040 
0041 %if ~isempty(nkeep); topcs=topcs(:,1:nkeep); end
0042 
0043 % power per PC
0044 pwr=diag(topcs'*c*topcs)/(m*o);
0045 idx=find(pwr>=threshold*max(pwr));
0046 pwr=pwr(idx)';
0047 topcs=topcs(:,idx);
0048 
0049 % PCs
0050 if nargout>2
0051     y=nt_mmat(x,topcs);
0052 end
0053 
0054 %% test code
0055 if 0
0056     x=randn(1000,10);
0057     [topcs,pwr,y]=nt_pca0(x);
0058     figure(1); plot(pwr);
0059     figure(2); subplot 121; plot(y); subplot 122; plot(x*topcs);
0060 end
0061 if 0
0062     x=zeros(1000,10);
0063     [topcs,pwr,y]=nt_pca0(x);
0064     figure(1); plot(pwr);
0065     figure(2); subplot 121; plot(y); subplot 122; plot(x*topcs);
0066 end
0067 if 0 
0068     x=sin(2*pi*3*(1:1000)'/1000)*randn(1,10);
0069     x=2*x + randn(size(x));
0070      [topcs,pwr,y]=nt_pca0(x);
0071     figure(1); plot(pwr);
0072     figure(2); subplot 121; plot(x); subplot 122; plot(x*topcs);
0073 end   
0074 
0075 
0076

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005