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: 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

Generated on Wed 29-Apr-2015 15:09:19 by m2html © 2005