Description of nt

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
