Home > NoiseTools > nt_star3.m

nt_star3

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

  y: denoised data 
  w: 0 for parts that needed fixing, 1 elsewhere (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)
 
 See also: nt_sns, nt_proximity

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x,w,ww]=nt_star3(x,thresh,closest,w)
0002 % [y,w,ww]=nt_star3(x,thresh,closest,w) - sensor noise suppression
0003 %
0004 %  y: denoised data
0005 %  w: 0 for parts that needed fixing, 1 elsewhere (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 %
0012 % See also: nt_sns, nt_proximity
0013 
0014 %{
0015 26 Jan 2017
0016 Regress also on polynomial to handle trends.
0017 %}
0018 
0019 
0020 
0021 nt_greetings;
0022 
0023 PCA_THRESH=10^-15;  % threshold for discarding weak PCs
0024 NSMOOTH=10;         % samples, smoothing to apply to excentricity
0025 MINPROP=0.3;        % minimum proportion of artifact free at first iteration
0026 NITER=3;            % iterations to refine c0
0027 VERBOSE=1;          % set to 0 to shut up
0028 
0029 if nargin<1; error; end
0030 if nargin<2 || isempty(thresh); thresh=1; end
0031 if nargin<3; closest=[]; end
0032 if ~isempty(closest)&&size(closest,1)~=size(x,2);
0033     error('closest array should have as many rows as channels of x'); 
0034 end
0035 if nargin<4 || isempty(w); w=1; end
0036 
0037 if nargout==0 % plot, don't return result
0038     [y,w,ww]=nt_star3(x,thresh,closest,w);
0039     disp([mean(w(:)), mean(ww(:))])
0040     figure(1); clf;
0041     subplot 311; plot(x);
0042     subplot 312; plot(y);
0043     subplot 313; plot(w, '.-'); ylim([0 1.1]);
0044     clear x w ww
0045     return
0046 end
0047     
0048 
0049 [nsample,nchan,~]=size(x);
0050 x=nt_unfold(x);
0051 
0052 p0=nt_wpwr(x);
0053 mn=mean(x); % save means
0054 x=nt_demean(x);
0055 nn=sqrt(mean(x.^2)); % save norm
0056 x=nt_normcol(x);
0057 p00=nt_wpwr(x);
0058 
0059 % NaN and zero channels are set to rand, which effectively excludes them
0060 iNan=find(all(isnan(x)));
0061 iZero=find(all(x==0));
0062 x(:,iNan)=randn(size(x,1),numel(iNan));
0063 x(:,iZero)=randn(size(x,1),numel(iZero));
0064 
0065 order=1;
0066 r=zeros(size(x,1),numel(order));
0067 lin=linspace(-1,1,size(x,1));
0068 for k=1:order
0069     r(:,k)=lin.^k;
0070 end
0071 
0072 x=nt_demean(x);
0073 c0=nt_cov(x); % initial covariance estimate
0074 
0075 %{
0076 Find time intervals where at least one channel is excentric --> w==0.
0077 %}
0078 
0079 iIter=NITER;
0080 while iIter>0
0081     
0082     
0083     w=ones(size(x,1),1);
0084     for iChan=1:nchan
0085 
0086         % other channels
0087         if ~isempty(closest); 
0088             oChan=closest(iChan,:);
0089         else
0090             oChan=setdiff(1:nchan,iChan);
0091         end
0092         oChan(oChan>nchan)=[];
0093         
0094         % PCA other channels to remove weak dimensions
0095         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
0096         idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
0097         topcs=topcs(:,idx);
0098         
0099         % project this channel on other channels
0100         b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % projection matrix
0101         y=x(:,oChan)*(topcs*b'); % projection
0102         dx=abs(y-x(:,iChan));   % difference from projection
0103         dx=dx+eps;              % avoids error on simulated data
0104         
0105         d=dx/mean(dx(find(w))); % excentricity measure
0106         if NSMOOTH>0; 
0107             d=filtfilt(ones(NSMOOTH,1)/NSMOOTH,1,d);
0108         end
0109         
0110         d=d/thresh;
0111         w=min(w,(d<1)); % w==0 for artifact part
0112         
0113     end    
0114     
0115     prop=mean(w);
0116     if VERBOSE>0; disp(['proportion artifact free: ', num2str(prop)]); end
0117     
0118     if iIter==NITER && prop<MINPROP
0119         thresh=thresh*1.1;
0120         if VERBOSE>0; disp(['Warning: nt_star increasing threshold to ', num2str(thresh)]); end
0121         w=ones(size(w));
0122     else
0123         iIter=iIter-1;
0124     end
0125     
0126     x=nt_demean(x,w);
0127     c0=nt_cov(x,[],w); % restrict covariance estimate to non-artifact part
0128 end
0129 
0130 %{
0131 We now know which part contains channel-specific artifacts (w==0 for artifact part), 
0132 and we have an estimate of the covariance matrix of the artifact-free part.
0133 %}
0134 
0135 %{
0136 Find which channel is most excentric at each time point.
0137 Here we use an excentricity measure based on the absolute value of the signal,
0138 rather than the difference between signal and projection.
0139 %}
0140 
0141 xx=abs(x);
0142 xx=bsxfun(@times,xx, 1 ./ sqrt(mean(xx(find(w),:).^2))); % divide by std over non-artifact part
0143 if NSMOOTH>0; 
0144     xx=filtfilt(ones(NSMOOTH,1)/NSMOOTH,1,xx);
0145 end
0146 [~,rank]=sort(xx','descend'); 
0147 rank=rank';
0148 rank(find(w),:)=0;      % exclude parts that were not judged excentric
0149 
0150 depth=1;
0151 depth=min(depth,nchan-1);
0152 ww=ones(size(x));
0153 for iDepth=1:depth
0154 
0155     %{
0156     Fix each channel by projecting on other channels.
0157     %}
0158     
0159     iFixed=nchan;
0160     nFixed=0;
0161     for iChan=1:nchan
0162 
0163         bad_samples=find(iChan==rank(:,iDepth)); % samples where this channel is the most excentric
0164         if iDepth ~=1; 
0165             bad_samples(find(xx(bad_samples,iChan)<thresh)) =[]; % exclude if not very bad
0166         end
0167         
0168         nFixed=nFixed+numel(bad_samples);
0169         if isempty(bad_samples); 
0170             iFixed=iFixed-1;
0171             continue;
0172         end
0173         ww(bad_samples,iChan)=0;
0174 
0175         % other channels
0176         if ~isempty(closest); 
0177             oChan=closest(iChan,:);
0178         else
0179             oChan=setdiff(1:nchan,iChan);
0180         end
0181         oChan(oChan>nchan)=[]; % in case closest includes channels not in data
0182 
0183         % PCA other channels to remove weak dimensions
0184         [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan)); % PCA
0185         idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH); % discard weak dims
0186         topcs=topcs(:,idx);
0187 
0188         % project this channel on other channels
0189         b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs); % projection matrix
0190         y=x(bad_samples,oChan)*(topcs*b'); % projection
0191 
0192         x(bad_samples,iChan)=y(:); % fix
0193 
0194     end
0195     
0196     if VERBOSE>0; 
0197         disp(['depth: ', num2str(iDepth), ', n fixed channels: ',num2str(iFixed),...
0198             ', n fixed samples: ', num2str(nFixed), ', ratio: ',num2str(nt_wpwr(x)/p00)]);
0199     end
0200 end
0201 
0202 x=nt_demean(x);
0203 x=bsxfun(@times,x,nn);
0204 x=bsxfun(@plus,x,mn);
0205 
0206 x=nt_fold(x,nsample);
0207 w=nt_fold(w,nsample);
0208 ww=nt_fold(ww,nsample);
0209 
0210 
0211 
0212 
0213 x(:,iNan)=nan;
0214 x(:,iZero)=0;
0215 
0216 if VERBOSE>0; disp(['power ratio: ', num2str(nt_wpwr(x)/p0)]); end

Generated on Tue 09-May-2017 13:19:57 by m2html © 2005