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

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