Home > NoiseTools > nt_star2.m

nt_star2

PURPOSE ^

[y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression

SYNOPSIS ^

function [x,w,ww]=nt_star2(x,thresh,closest,w)

DESCRIPTION ^

 [y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression

  y: denoised data 
  w: 1 for parts used to calculate covariance (time*1)
  ww: 0 for parts that needed fixing, 1 elsewhere (time*chans)

  x: data to denoise (time*chans or time*chans*trials)
  thresh: threshold for excentricity measure (default:1);
  closest: indices of channels that are closest to each channel (default: all)
  w: initial covariance weight matrix (time*1)
 
 See also: nt_sns, nt_proximity

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [x,w,ww]=nt_star2(x,thresh,closest,w)
0002 % [y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression
0003 %
0004 %  y: denoised data
0005 %  w: 1 for parts used to calculate covariance (time*1)
0006 %  ww: 0 for parts that needed fixing, 1 elsewhere (time*chans)
0007 %
0008 %  x: data to denoise (time*chans or time*chans*trials)
0009 %  thresh: threshold for excentricity measure (default:1);
0010 %  closest: indices of channels that are closest to each channel (default: all)
0011 %  w: initial covariance weight matrix (time*1)
0012 %
0013 % See also: nt_sns, nt_proximity
0014 
0015 nt_greetings;
0016 
0017 PCA_THRESH=10^-15;  % threshold for discarding weak PCs
0018 NSMOOTH=10;         % samples, smoothing to apply to excentricity
0019 MINPROP=0.4;        % minimum proportion of artifact free at first iteration
0020 NITER=4;            % iterations to refine c0
0021 VERBOSE=1;          % set to 0 to shut up
0022 THRESH_RATIO=2;     % ratio of shared-weight and per-channel-weight thresholds
0023 
0024 if nargin<1; error; end
0025 if nargin<2 || isempty(thresh); thresh=1; end
0026 if nargin<3; closest=[]; end
0027 if ~isempty(closest)&&size(closest,1)~=size(x,2);
0028     error('closest array should have as many rows as channels of x'); 
0029 end
0030 if nargin<4; w=[]; end
0031 
0032 if nargout==0 % plot, don't return result
0033     [y,w,ww]=nt_star2(x,thresh,closest,w);
0034     disp([mean(w(:)), mean(ww(:))])
0035     figure(1); clf;
0036     subplot 411; plot(x); mx=max(x(:)); mn=min(x(:)); ylim([mn mx]); title('raw');
0037     subplot 412; plot(y);ylim([mn,mx]); title('clean')
0038     subplot 413; plot(x-y); title('diff');
0039     subplot 414; plot(w, '.-'); ylim([0 1.1]); title('weight');
0040     clear x w ww
0041     return
0042 end
0043 
0044 [nsample,nchan,~]=size(x);
0045 x=nt_unfold(x);
0046 
0047 p0=nt_wpwr(x);
0048 mn=mean(x); % save means
0049 x=nt_demean(x);
0050 nn=sqrt(mean(x.^2)); % save norm
0051 x=nt_normcol(x);
0052 p00=nt_wpwr(x);
0053 
0054 % NaN, zero, and zero-weight channels are set to rand, which effectively excludes them
0055 iNan=find(all(isnan(x)));
0056 iZero=find(all(x==0));
0057 x(:,iNan)=randn(size(x,1),numel(iNan));
0058 x(:,iZero)=randn(size(x,1),numel(iZero));
0059 
0060 x=nt_demean(x);
0061 
0062 if isempty(w); w=ones(size(x,1),1); end
0063 if size(w,2)>1; w=min(w,[],2); end
0064 c0=nt_cov(x,[],w);          % initial covariance estimate
0065 ww=ones(size(x));           % per-channel weights
0066 
0067 %{
0068 Find time intervals where at least one channel is excentric --> w==0.
0069 %}
0070 
0071 iIter=NITER;
0072 while iIter>0
0073       
0074     for iChan=1:nchan
0075 
0076         % other channels
0077         if ~isempty(closest); 
0078             oChan=closest(iChan,:);
0079         else
0080             oChan=setdiff(1:nchan,iChan);
0081         end
0082         oChan(oChan>nchan)=[];
0083         
0084         % PCA other channels to remove weak dimensions
0085         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
0086         idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
0087         topcs=topcs(:,idx);
0088         
0089         % project this channel on other channels
0090         b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % projection matrix
0091         y=x(:,oChan)*(topcs*b'); % projection
0092                 
0093         dx=abs(y-x(:,iChan));   % difference from projection
0094         dx=dx+eps;              % avoids error on simulated data
0095         d=dx/mean(dx(find(w))); % excentricity measure
0096         d(find(isnan(d)))=0;
0097         if NSMOOTH>0; 
0098             d=filtfilt(ones(NSMOOTH,1)/NSMOOTH,1,d);
0099         end
0100        
0101         thresh1=thresh;     % threshold for shared weight, for estimating variance
0102         thresh2=thresh/THRESH_RATIO;   % threshold for channel-specific weight, to define clean subspace
0103 
0104         w=min(w,d<thresh1);     % w==0 for artifact part
0105         ww(:,iChan)=d<thresh2;  % weights specific to each channel
0106         
0107    end    
0108     
0109     prop=mean(w);
0110     if VERBOSE>0; disp(['proportion artifact free: ', num2str(prop)]); end    
0111     if iIter==NITER && prop<MINPROP
0112         thresh=thresh*1.1;
0113         if VERBOSE>0; disp(['Warning: nt_star increasing threshold to ', num2str(thresh)]); end
0114         w=ones(size(w));
0115     else
0116         iIter=iIter-1;
0117     end
0118     
0119     x=nt_demean(x,w);
0120     c0=nt_cov(x,[],w);      % recalculate covariance on non-artifact part
0121     w=ones(size(x,1),1);    % reset
0122 end
0123 ww=nt_growmask(ww,10);
0124 
0125 
0126 %{
0127 We now know which part contains channel-specific artifacts (w==0 for artifact part), 
0128 and we have an estimate of the covariance matrix of the artifact-free part.
0129 %}
0130 
0131 figure(1); clf; plot(mean(ww,2)); 
0132 y=nt_inpaint(x,ww,w);
0133 
0134 x=y;
0135 
0136 
0137 x=nt_demean(x);
0138 x=bsxfun(@times,x,nn);
0139 x=bsxfun(@plus,x,mn);
0140 x=nt_fold(x,nsample);
0141 w=nt_fold(w,nsample);
0142 ww=nt_fold(ww,nsample);
0143 x(:,iNan)=nan;
0144 x(:,iZero)=0;
0145 if VERBOSE>0; disp(['power ratio: ', num2str(nt_wpwr(x)/p0)]); end
0146 
0147 
0148 %% test code
0149 if 0
0150     %[p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/BQ/BQ_aespa_speech_47.mat');
0151     [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/MC/MC_aespa_speech_14.mat');
0152     %[p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/NM/NM_aespa_speech_12.mat');
0153     x=x'; x=x(:,1:128); x=x(:,:);
0154     x=nt_demean(x);
0155     [x,w]=nt_detrend_robust(x,10); 
0156     iBad=nt_badChannels(x,[],5000); 
0157     if ~isempty(iBad)
0158         disp(['bad channels: ', num2str(iBad)]);
0159         w(:,iBad)=0;
0160     end
0161     
0162     [xx,w]=nt_star2(x,3,[],w); 
0163     figure(1); clf; 
0164     subplot 311; plot(x); subplot 312; plot(xx); subplot 313; plot (x-xx);
0165     ch=1; figure(2); clf; subplot 211; plot ([x(:,ch),xx(:,ch)]); subplot 212; plot (x(:,ch)-xx(:,ch));
0166 end   
0167 
0168 function [ww,nOccurrences,iBack]=patternDict(w);
0169 % ww: dictionary of patterns
0170 % nOccurrences: number of times each pattern occurred
0171 % iBack: index to reconstruct input from dictionary
0172 [ww,~,IC,nOccurrences]=nt_unique(w,'rows');
0173 [nOccurrences,iSort]=sort(nOccurrences, 'descend'); % sort by decreasing number
0174 [~,iReverse]=sort(iSort); %
0175 ww=ww(iSort,:); % same order for patterns, w = ww(iReverse1(IC),:)
0176 iBack=iReverse(IC); % w = ww(iBack,:)
0177 
0178 function [w,x]=coalesceWeights(w,x)
0179 % reduce the number of weight patterns based on pruning & interpolation, to
0180 % save time
0181 nchans=size(w,2);
0182 w1=w(1:end-2,:); 
0183 w2=w(2:end-1,:); 
0184 w3=w(3:end,:); 
0185 x1=x(1:end-2,:); 
0186 x2=x(2:end-1,:); 
0187 x3=x(3:end,:); 
0188 % weights that correspond to a single sample
0189 idx=find(~w1&w2&~w3); % isolated good, set to bad
0190 w2(idx)=0; w(2:end-1,:)=w2;
0191 idx=find(w1&~w2&w3); % isolated bad, interpolate, set to good
0192 w2(idx)=1; w(2:end-1,:)=w2; 
0193 x2(idx)=(x1(idx))+x3(idx)/2; x(2:end-1,:)=x2;
0194 % weight patterns that occur rarely
0195 [wDict,nOccurrences,iBack]=patternDict(w); 
0196 idx=find(nOccurrences<=3);
0197 isolated=iBack(idx);
0198 for k=1:numel(isolated)
0199     iIsolated=isolated(k);
0200     if iIsolated-1>=1 && iIsolated+1<=size(w,1)            
0201         if sum(abs(w(iIsolated-1,:)-w(iIsolated,:)))<sum(abs(w(iIsolated,:)-w(iIsolated+1,:)));
0202             iFix=iIsolated-1;% preceding is closest
0203         else
0204             iFix=iIsolated+1;% following is closest
0205         end
0206         idx=find(w(iIsolated,:)&~w(iFix,:));
0207         w(iIsolated,idx)=0;
0208         idx=find(~w(iIsolated,:)&w(iFix,:));
0209         w(iIsolated,idx)=1;
0210         x(iIsolated,idx)=(x(iIsolated-1,idx)+x(iIsolated+1,idx))/2;
0211     end
0212 end
0213 
0214     
0215     
0216         
0217

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005