Home > NoiseTools > nt_sns4.m

# nt_sns4

## PURPOSE

[y,w]=nt_sns2(x,thresh) - sensor noise suppression, new version

## SYNOPSIS

function [x,w]=nt_sns2(x,thresh)

## DESCRIPTION

     [y,w]=nt_sns2(x,thresh) - sensor noise suppression, new version

      y: denoised data
      w: 0 for parts that needed fixing, 1 elsewhere (time*chans)

      x: data to denoise (time*chans or time*chans*trials)
      thresh: threshold for Mahalanobis distance (default:1);

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_cov2 [c,tw]=nt_cov2(x,w) - weighted covariance
- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1
- nt_pcarot [topcs,eigenvalues]=pcarot(cov,N) - PCA matrix from covariance
- nt_unfold y=nt_fold(x) - unfold 3D to 2D

This function is called by:

## SOURCE CODE

    0001 function [x,w]=nt_sns2(x,thresh)
    0002 % [y,w]=nt_sns2(x,thresh) - sensor noise suppression, new version
    0003 %
    0004 %  y: denoised data
    0005 %  w: 0 for parts that needed fixing, 1 elsewhere (time*chans)
    0006 %
    0007 %  x: data to denoise (time*chans or time*chans*trials)
    0008 %  thresh: threshold for Mahalanobis distance (default:1);
    0009 %
    0010 nt_greetings;
    0011 if nargin<1; error; end
    0012 if nargin<2 || isempty(thresh); thresh=1; end
    0013 PCA_THRESH=10^-15;
    0014 [nsample,nchan,~]=size(x);
    0015
    0016 x=nt_unfold(x);
    0017 mn=mean(x); % save means
    0018 x=nt_demean(x);
    0019 nn=sqrt(mean(x.^2)); % save norm
    0020 x=nt_normcol(x);
    0021
    0022 %{
    0023 For each channel, find sections for which it does not fit the
    0024 subspace spanned by other sensors.  The calculation is repeated
    0025 and the projection matrix is refined at each step.
    0026 %}
    0027
    0028 w=ones(size(x));
    0029
    0030 NITER=3; % iterations to refine c0
    0031 for k=1:NITER
    0032
    0033     % c0: covariance of non-artifact part
    0034     ww=min(w,[],2);
    0035     x=nt_demean(x,ww);
    0036     [c0,tw]=nt_cov2(x,ww);
    0037     c0=c0./tw; c0(isnan(c0))=0;
    0038
    0039     for iChan=1:nchan
    0040
    0041         % regress this channel on all others
    0042         oChan=setdiff(1:nchan,iChan);
    0043
    0044         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
    0045         idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
    0046         topcs=topcs(:,idx);
    0047         b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % matrix to project on other channels
    0048         y(:,iChan)=(x(:,oChan)*topcs)*b'; % projection
    0049
    0050         % fix parts where other channels are bad
    0051         for iChan2=oChan
    0052             oChan2=setdiff(oChan,iChan2);
    0053             iBad=find(w(:,iChan2)==0);
    0054             if ~isempty(iBad)
    0055                 [topcs,eigenvalues]=nt_pcarot(c0(oChan2,oChan2)); % PCA
    0056                 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
    0057                 topcs=topcs(:,idx);
    0058                 b=c0(iChan,oChan2)*topcs/(topcs'*c0(oChan2,oChan2)*topcs); % matrix to project on other channels
    0059                 y(iBad,iChan)=(x(iBad,oChan2)*topcs)*b'; % projection
    0060             end
    0061         end
    0062
    0063
    0064
    0065         dd=y(:,iChan)-x(:,iChan); % difference from projection
    0066         %plot([y(:,iChan),x(:,iChan)]); pause
    0067         d=mahal(dd,dd)/thresh; % excentricity of each sample
    0068
    0069         w(:,iChan)=(d<1);
    0070
    0071     end
    0072     %figure(10); clf; imagescc(w);  pause
    0073     disp(mean(w(:)))
    0074 end
    0075
    0076 x=y;
    0077
    0078 %{
    0079 To do:
    0080 Record all the DC shifts introduced when x mean is removed, so as to
    0081 restore accurately.
    0082 %}
    0083
    0084 x=nt_demean(x);
    0085 x=bsxfun(@times,x,nn);
    0086 x=bsxfun(@plus,x,mn);
    0087
    0088 x=nt_fold(x,nsample);
    0089 w=nt_fold(w,nsample);
    0090
    0091
    0092