Home > NoiseTools > nt_sns3_old.m

nt_sns3_old

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 
0029 y=zeros(size(x));
0030 d=ones(size(x,1),nchan); % distance from projection, d>thresh means artifact
0031 
0032 NITER=2; % iterations to refine c0
0033 for k=1:NITER
0034     
0035     % c0: covariance of non-artifact part
0036     dd=max(d,[],2);
0037     c0=nt_cov(x,[],dd<=1);
0038     
0039     d0=d;
0040     for iChan=1:nchan
0041 
0042         oChan=setdiff(1:nchan,iChan); % other channels
0043         
0044         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
0045         PCA_THRESH=10^-6;
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 
0051         y0=y;
0052         
0053         % fix to exclude the most serious competing channels from projection
0054         competitors=zeros(size(d0(:,oChan)));
0055         iCompetitors=d0(:,oChan)> 1 ... % bad
0056             & d0(:,oChan) > repmat(d0(:,iChan),1,nchan-1); % and worse than this one
0057         competitors(find(iCompetitors))=d0(find(iCompetitors));
0058         idx=find(~min(competitors,[],2)); % samples with no competitors
0059         [~,worst]=max(competitors,[],2); % which competitor dominates?
0060         worst(idx)=0; % set to zero if there is none
0061         for k=1:numel(oChan);
0062             
0063             iReplace=find(k==worst);
0064             if iReplace
0065                 ooChan=setdiff(1:nchan,[iChan,oChan(k)]); % channels other than this and worst
0066                 [topcs,eigenvalues]=nt_pcarot(c0(ooChan,ooChan)); % PCA
0067                 PCA_THRESH=10^-6;
0068                 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
0069                 topcs=topcs(:,idx);
0070                 b=c0(iChan,ooChan)*topcs/(topcs'*c0(ooChan,ooChan)*topcs);
0071                 y(iReplace,iChan)=(x(iReplace,ooChan)*topcs)*b'; % projection
0072             end
0073         end
0074         
0075         
0076         ddd=y(:,iChan)-x(:,iChan); % difference from projection
0077         
0078         d(:,iChan)=mahal(ddd,ddd)/thresh; % excentricity of each sample and each channel
0079 
0080         figure(1); clf
0081         FOCUS=861600+(1:500); subplot 311; plot(x(FOCUS,:)); title('raw');
0082         subplot 312; plot([x(FOCUS,iChan), y0(FOCUS,iChan), y(FOCUS,iChan)]); title (num2str(iChan)); legend('raw','proj','fixed');
0083         subplot 313; plot(d(FOCUS,:)); pause
0084 
0085     end    
0086     %disp(mean(w))
0087 end
0088 
0089 dd=max(d,[],2);
0090 %x=nt_demean(x,dd<1);
0091 c0=nt_cov(x,[],dd<1);
0092 
0093 
0094 %{
0095 For each channel, find the part for which it dominates the artifact.
0096 Replace that part based on projection on the other channels. 
0097 %}
0098 
0099 idx=find(dd>1);
0100 xx=x(idx,:);
0101 xx=abs(nt_normcol(xx));
0102 
0103 for iChan=1:nchan
0104     
0105 %     oChan=setdiff(1:nchan,iChan);
0106 %     b=c0(iChan,oChan)/c0(oChan,oChan);    % projection matrix
0107 %     idx2=find(xx(:,iChan)>max(xx(:,oChan),[],2)); % find sections where this channel is largest
0108 %     idx2=find(xx(:,iChan)>max(xx(:,oChan),[],2)); % find sections where this channel is largest
0109 %     x(idx(idx2),iChan)=x(idx(idx2),oChan)*b'; % replace by projection
0110 
0111     oChan=setdiff(1:nchan,iChan);
0112     b=c0(iChan,oChan)/c0(oChan,oChan);    % projection matrix
0113     idx2=find(xx(:,iChan)>max(xx(:,oChan),[],2)); % find sections where this channel is largest
0114     
0115     x(idx(idx2),iChan)=x(idx(idx2),oChan)*b'; % replace by projection
0116     
0117     idx2=find(d(:,iChan)>1 & d(:,iChan)>max(d(:,oChan),[],2));
0118     %x(idx2,iChan)=x(idx2,oChan)*b';
0119 
0120 end
0121 
0122 %{
0123 To do:
0124 - avoid problems with rank-deficient data
0125 - exclude borderline channels from projection
0126 %}
0127 
0128 x=bsxfun(@times,x,nn);
0129 x=bsxfun(@plus,x,mn);
0130 
0131 x=nt_fold(x,nsample);
0132 d=nt_fold(d,nsample);
0133 
0134 
0135

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