Home > NoiseTools > nt_star2.m

# nt_star2

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

     [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

## 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_detrend [y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_growmask ww=nt_growmask(w,margin) - widen mask
- nt_inpaint function y=nt_inpaint(x,w) - weighted interpolation based on correlation structure
- 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_read_data [p,data]=nt_read_data(fname,flag) - read data from file
- nt_star2 [y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression
- nt_unfold y=nt_fold(x) - unfold 3D to 2D
- nt_unique [C,IA,IC,N] = nt_unique(A, varargin) - unique with counts
- nt_wpwr [y,tweight]=nt_wpwr(x,w) - weighted power

This function is called by:

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

## SUBFUNCTIONS

- function [ww,nOccurrences,iBack]=patternDict(w);
- function [w,x]=coalesceWeights(w,x)

## SOURCE CODE

    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