Home > NoiseTools > nt_star.m

# nt_star

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

     [y,w,ww]=nt_star(x,thresh,closest,depth) - 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)
      depth: maximum number of channels to fix at each sample (default 1)

     See also: nt_sns, nt_proximity

## 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,nkeep,threshold,N) - PCA matrix from covariance
- nt_star [y,w,ww]=nt_star(x,thresh,closest,depth) - sensor noise suppression
- nt_unfold y=nt_fold(x) - unfold 3D to 2D
- nt_wpwr [y,tweight]=nt_wpwr(x,w) - weighted power

This function is called by:

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

## SOURCE CODE

    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