nt_whiten

[A,y]=nt_whiten(x,N) - whiten spectrally using pca

  A: whitening matrix (to be applied to time-shifted x)
  y: whitened signal 

  x: signal to whiten
  N: order (number of time shifts)

0001 function [A,y,AA]=nt_whiten(x,N)
0002 %[A,y]=nt_whiten(x,N) - whiten spectrally using pca
0003 %
0004 %  A: whitening matrix (to be applied to time-shifted x)
0005 %  y: whitened signal
0006 %
0007 %  x: signal to whiten
0008 %  N: order (number of time shifts)
0009 %
0010 
0011 if nargin<2; error('!'); end
0012 
0013 % calculate covariance across time shifts, looping over columns to save space
0014 sz=size(x);
0015 xx=reshape(x,sz(1),prod(sz(2:end)));
0016 C=zeros(N); % covariance of time-shifted data
0017 for iCol=1:size(xx,2);
0018     xxx=nt_multishift(xx(:,iCol),0:N-1); 
0019     C=C+xxx'*xxx;
0020 end
0021 C=C/size(x,1);
0022 
0023 % PCA, normalize, inverse
0024 [topcs,evs]=nt_pcarot(C);
0025 tmp=1./evs; tmp(find(evs<=0))=0;
0026 A=topcs*diag(tmp)*pinv(topcs);
0027 B=pinv(A);
0028 AA=A(:,2:end)*B(2:end,:);
0029 
0030 
0031 if nargout>1
0032     % apply whitening matrix, keeping only 1st column
0033     yy=zeros(size(xx));
0034     for iCol=1:size(yy,2)
0035         xxx=nt_multishift(xx(:,iCol),0:N-1);
0036         yy(1:size(xxx,1),iCol)=xxx*A(:,1);
0037     end
0038     y=reshape(yy,[size(yy,1),sz(2:end)]);
0039 end
0040