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:

- 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_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_greetings nt_greetings - display message the first time the toolbox is used
- 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
    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