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

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