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:

- 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_star3 [y,w,ww]=nt_star3(x,thresh,closest,w) - 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_star3 [y,w,ww]=nt_star3(x,thresh,closest,w) - sensor noise suppression

## 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