nt_star2

 [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

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 function [ww,nOccurrences,iBack]=patternDict(w);
0148 % ww: dictionary of patterns
0149 % nOccurrences: number of times each pattern occurred
0150 % iBack: index to reconstruct input from dictionary
0151 [ww,~,IC,nOccurrences]=nt_unique(w,'rows');
0152 [nOccurrences,iSort]=sort(nOccurrences, 'descend'); % sort by decreasing number
0153 [~,iReverse]=sort(iSort); %
0154 ww=ww(iSort,:); % same order for patterns, w = ww(iReverse1(IC),:)
0155 iBack=iReverse(IC); % w = ww(iBack,:)
0156 
0157 function [w,x]=coalesceWeights(w,x)
0158 % reduce the number of weight patterns based on pruning & interpolation, to
0159 % save time
0160 nchans=size(w,2);
0161 w1=w(1:end-2,:); 
0162 w2=w(2:end-1,:); 
0163 w3=w(3:end,:); 
0164 x1=x(1:end-2,:); 
0165 x2=x(2:end-1,:); 
0166 x3=x(3:end,:); 
0167 % weights that correspond to a single sample
0168 idx=find(~w1&w2&~w3); % isolated good, set to bad
0169 w2(idx)=0; w(2:end-1,:)=w2;
0170 idx=find(w1&~w2&w3); % isolated bad, interpolate, set to good
0171 w2(idx)=1; w(2:end-1,:)=w2; 
0172 x2(idx)=(x1(idx))+x3(idx)/2; x(2:end-1,:)=x2;
0173 % weight patterns that occur rarely
0174 [wDict,nOccurrences,iBack]=patternDict(w); 
0175 idx=find(nOccurrences<=3);
0176 isolated=iBack(idx);
0177 for k=1:numel(isolated)
0178     iIsolated=isolated(k);
0179     if iIsolated-1>=1 && iIsolated+1<=size(w,1)            
0180         if sum(abs(w(iIsolated-1,:)-w(iIsolated,:)))<sum(abs(w(iIsolated,:)-w(iIsolated+1,:)));
0181             iFix=iIsolated-1;% preceding is closest
0182         else
0183             iFix=iIsolated+1;% following is closest
0184         end
0185         idx=find(w(iIsolated,:)&~w(iFix,:));
0186         w(iIsolated,idx)=0;
0187         idx=find(~w(iIsolated,:)&w(iFix,:));
0188         w(iIsolated,idx)=1;
0189         x(iIsolated,idx)=(x(iIsolated-1,idx)+x(iIsolated+1,idx))/2;
0190     end
0191 end
0192 
0193     
0194 
0195 
0196 %% test code
0197 if 0
0198     %[p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/BQ/BQ_aespa_speech_47.mat');
0199     %[p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/MC/MC_aespa_speech_14.mat');
0200     [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/NM/NM_aespa_speech_12.mat');
0201     x=x'; x=x(:,1:128); x=x(:,:);
0202     x=nt_demean(x);
0203     [x,w]=nt_detrend(x,10); 
0204     iBad=nt_badChannels(x,[],5000); 
0205     if ~isempty(iBad)
0206         disp(['bad channels: ', num2str(iBad)]);
0207         w(:,iBad)=0;
0208     end
0209     
0210     [xx,w]=nt_star2(x,3,[],w); 
0211     figure(1); clf; 
0212     subplot 311; plot(x); title('raw'); subplot 312; plot(xx); title('clean'); subplot 313; plot (x-xx); title('raw-clean');
0213     ch=1; figure(2); clf; subplot 211; plot ([x(:,ch),xx(:,ch)]); title('single channel'); legend('raw','clean'); subplot 212; plot (x(:,ch)-xx(:,ch)); title('raw-clean')
0214 end   
0215 
0216         
0217