nt_pca0

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

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

  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

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