Home > NoiseTools > nt_sns0.m

# nt_sns0

## PURPOSE

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

## SYNOPSIS

function r=nt_sns0(c,nneighbors,skip,wc)

## DESCRIPTION

     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

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_pcarot [topcs,eigenvalues]=pcarot(cov,nkeep,threshold,N) - PCA matrix from covariance
- nt_vecmult y=nt_vecmult(x,v) - multiply all rows or columns of matrix by vector

This function is called by:

- nt_sns y=nt_sns(x,nneigbors,skip,w) - sensor noise suppression

## SOURCE CODE

    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