nt_sns0

 r=nt_sns0(c,nneigbors,skip,wc) - sensor noise suppression

   r: denoising matrix

   c: full covariance of data to denoise
   nneighbors: number of channels to use in projection 
   skip: number of neighbors to skip [default: 0]
   wc: weighted covariance

0001 function r=nt_sns0(c,nneighbors,skip,wc)
0002 % r=nt_sns0(c,nneigbors,skip,wc) - sensor noise suppression
0003 %
0004 %   r: denoising matrix
0005 %
0006 %   c: full covariance of data to denoise
0007 %   nneighbors: number of channels to use in projection
0008 %   skip: number of neighbors to skip [default: 0]
0009 %   wc: weighted covariance
0010 %
0011 %
0012 
0013 
0014 n=size(c,1);
0015 
0016 if nargin<2 || isempty(nneighbors); error('need to specify nneighbors'); end
0017 if nargin<3 || isempty(skip); skip=0; end
0018 if nargin<4 || isempty(wc); wc=c; end
0019 
0020 nneighbors=min(nneighbors,n-skip-1);
0021 
0022 r=zeros(size(c));
0023 
0024 % normalize
0025 d=sqrt(1./(diag(c)+eps));
0026 c=nt_vecmult(nt_vecmult(c,d),d');
0027 
0028 for iChan=1:n
0029  
0030     c1=c(:,iChan);                      % correlation of channel with all other channels
0031     [c1,idx]=sort(c1.^2,1,'descend');   % sort by correlation
0032     idx=idx(skip+2:skip+1+nneighbors);  % keep best
0033 
0034     % pca neighbors to orthogonalize them
0035     c2=wc(idx,idx);
0036     [topcs,eigenvalues]=nt_pcarot(c2);
0037     topcs=topcs*diag(1./sqrt(eigenvalues));
0038     topcs(find(isinf(topcs)|isnan(topcs)))=0;
0039     
0040     % augment rotation matrix to include this channel
0041     topcs=[1,zeros(1,nneighbors);zeros(nneighbors,1),topcs];
0042     
0043     % correlation matrix for rotated data
0044     c3=topcs'*wc([iChan;idx],[iChan;idx])*topcs;
0045     
0046     % first row defines projection to clean component iChan
0047     c4=c3(1,2:end)*topcs(2:end,2:end)';
0048 
0049     % insert new column into denoising matrix
0050     r(idx,iChan)=c4;
0051 
0052 end