Home > NoiseTools > nt_sns3.m

# nt_sns3

## PURPOSE

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

## SYNOPSIS

function [x,d]=nt_sns3(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_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- 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,d]=nt_sns3(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 if nargin<1; error; end
    0013 if nargin<2 || isempty(thresh); thresh=1; end
    0014
    0015 [nsample,nchan,~]=size(x);
    0016 x=nt_unfold(x);
    0017
    0018 mn=mean(x); % save means
    0019 x=nt_demean(x);
    0020 nn=sqrt(mean(x.^2)); % save norms
    0021 x=nt_normcol(x);
    0022
    0023 %{
    0024 For each channel, find sections for which it does not fit the
    0025 subspace spanned by other sensors.  The calculation is repeated
    0026 and the projection matrix is refined at each step.
    0027 %}
    0028 PCA_THRESH=10^-15;
    0029
    0030 d=zeros(size(x));
    0031
    0032
    0033 % c0: covariance of non-artifact part
    0034 x=nt_demean(x);
    0035 c0=nt_cov(x);
    0036
    0037 for iChan=1:nchan
    0038
    0039     oChan=setdiff(1:nchan,iChan);
    0040
    0041     % PCA to avoid problems with rank-deficient data
    0042     [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
    0043     idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
    0044     topcs=topcs(:,idx);
    0045
    0046     % regression matrix for this channel on all others
    0047     b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs);
    0048     y(:,iChan)=(x(:,oChan)*topcs)*b'; % projection
    0049     dd=y(:,iChan)-x(:,iChan); % difference from projection
    0050
    0051     d(:,iChan)=mahal(dd,dd)/thresh; % excentricity of each sample from distribution
    0052
    0053 end
    0054 %figure(1); clf; plot(d); pause
    0055
    0056 %{
    0057 Find which channel is most excentric at each moment.
    0058 %}
    0059 w=max(d,[],2)<1;
    0060 [~,iWorstChan]=sort(d','descend');
    0061 iWorstChan=iWorstChan';
    0062 iWorstChan(find(w),:)=nan;
    0063
    0064 x=nt_demean(x,w);
    0065 c0=nt_cov(x(find(w),:)); % covariance of non-artifact part
    0066 for iChan=1:nchan
    0067
    0068     idx=find(iWorstChan(:,1)==iChan);
    0069
    0070     ww=zeros(size(w));
    0071     ww(idx)=1;
    0072
    0073     if numel(idx)>0
    0074         c1=nt_cov(x(idx,:));
    0075         todss=nt_dss0(c0,c1);
    0076         z=nt_mmat(x,todss);
    0077         fromdss=pinv(todss);
    0078         x(idx,iChan)=nt_mmat(z(idx,2:end),fromdss(2:end,iChan));
    0079     end
    0080
    0081
    0082 end
    0083
    0084
    0085 x=bsxfun(@times,x,nn);
    0086 x=bsxfun(@plus,x,mn);
    0087
    0088 x=nt_fold(x,nsample);
    0089 d=nt_fold(d,nsample);
    0090