Home > NoiseTools > nt_sns2.m

# nt_sns2

## 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_cov [c,tw]=nt_cov(x,shifts,w) - time shift 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_pca [z,idx]=nt_pca(x,shifts,nkeep,threshold,w) - time-shift pca
- 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
    0012 PCA_THRESH=10^-15;
    0013
    0014 if nargin<1; error; end
    0015 if nargin<2 || isempty(thresh); thresh=1; end
    0016
    0017 [nsample,nchan,~]=size(x);
    0018 x=nt_unfold(x);
    0019
    0020 mn=mean(x); % save means
    0021 x=nt_demean(x);
    0022 nn=sqrt(mean(x.^2)); % save norm
    0023 x=nt_normcol(x);
    0024
    0025 %{
    0026 For each channel, find sections for which it does not fit the
    0027 subspace spanned by other sensors.  The calculation is repeated
    0028 and the projection matrix is refined at each step.
    0029 %}
    0030
    0031 w=ones(size(x,1),1);
    0032
    0033 NITER=2; % iterations to refine c0
    0034 for k=1:NITER
    0035
    0036     % c0: covariance of non-artifact part
    0037     x=nt_demean(x,w);
    0038     c0=nt_cov(x,[],w);
    0039
    0040     for iChan=1:nchan
    0041
    0042         % regression matrix for this channel on all others
    0043         oChan=setdiff(1:nchan,iChan);
    0044
    0045         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
    0046         idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
    0047         topcs=topcs(:,idx);
    0048         b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % matrix to project on other channels
    0049         y(:,iChan)=(x(:,oChan)*topcs)*b'; % projection
    0050         dd=y(:,iChan)-x(:,iChan); % difference from projection
    0051         %plot([y(:,iChan),x(:,iChan)]); pause
    0052         d=mahal(dd,dd)/thresh; % excentricity of each sample and each channel
    0053         w=min(w,(d<1));
    0054
    0055     end
    0056     disp(mean(w))
    0057 end
    0058
    0059 x=nt_demean(x,w);
    0060 c0=nt_cov(x,[],w);
    0061
    0062
    0063 %{
    0064 For each channel, find the part for which it dominates the artifact.
    0065 Replace that part based on projection on the other channels.
    0066 %}
    0067
    0068 %plot(w); pause
    0069
    0070 iBad=find(~w);
    0071 xx=x(iBad,:);
    0072 z=nt_pca(xx);
    0073 xx=abs(nt_normcol(xx));
    0074
    0075 %figure(1); clf; plot(xx); pause
    0076
    0077 %xx=nt_smooth(xx,20,[],1);
    0078
    0079 for iChan=1:nchan
    0080
    0081     iBad2=find(xx(:,iChan)>max(xx(:,setdiff(1:nchan,iChan)),[],2)); % this channel dominates others
    0082     oChan=setdiff(1:nchan,iChan);
    0083     [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
    0084     idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
    0085     topcs=topcs(:,idx);
    0086     b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % matrix to project on other channels
    0087     x(iBad(iBad2),iChan)=(x(iBad(iBad2),oChan)*topcs)*b'; % projection
    0088 end
    0089
    0090 %{
    0091 To do:
    0092 Record all the DC shifts introduced when x mean is removed, so as to
    0093 restore accurately.
    0094 %}
    0095
    0096
    0097 x=nt_demean(x);
    0098 x=bsxfun(@times,x,nn);
    0099 x=bsxfun(@plus,x,mn);
    0100
    0101 x=nt_fold(x,nsample);
    0102 w=nt_fold(w,nsample);
    0103
    0104
    0105