nt_tsr

[y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh) - time-shift regression (TSPCA)

  y: denoised data 
  idx: x(idx) is aligned with y
  w: weights applied by tsr
 
  x: data to denoise (time * channels * trials)
  ref: reference (time * channels * trials)
  shifts: array of shifts to apply to ref (default: [0])
  wx: weights to apply to x (time * 1 * trials);
  wref: weights to apply to ref (time * 1 * trials);
  keep: number of shifted-ref PCs to retain (default: all)
  thresh: ignore shifted-ref PCs smaller than thresh (default: 10.^-12)

 NoiseTools

0001 function [y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh)
0002 %[y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh) - time-shift regression (TSPCA)
0003 %
0004 %  y: denoised data
0005 %  idx: x(idx) is aligned with y
0006 %  w: weights applied by tsr
0007 %
0008 %  x: data to denoise (time * channels * trials)
0009 %  ref: reference (time * channels * trials)
0010 %  shifts: array of shifts to apply to ref (default: [0])
0011 %  wx: weights to apply to x (time * 1 * trials);
0012 %  wref: weights to apply to ref (time * 1 * trials);
0013 %  keep: number of shifted-ref PCs to retain (default: all)
0014 %  thresh: ignore shifted-ref PCs smaller than thresh (default: 10.^-12)
0015 %
0016 % NoiseTools
0017 nt_greetings;
0018 
0019 
0020 % Copyright 2007, 2008 Alain de Cheveigne
0021 
0022 % See:
0023 % de Cheveign\'e, A. and Simon, J. Z. (2007). "Denoising based on
0024 % Time-Shift PCA." Journal of Neuroscience Methods 165: 297-305.
0025 %
0026 % The basic idea is to project the signal X on a basis formed by the
0027 % orthogonalized time-shifted REF, and remove the projection. Supposing REF
0028 % gives a good observation of the noise that contaminates X, the noise is
0029 % removed. By allowing time shifts, the algorithm finds the optimal FIR filter
0030 % to apply to REF so as to compensate for any convolutional mismatch
0031 % between X and REF.
0032 %
0033 % Implementation issues:
0034 % - Data are often available as an array of epochs. This implementation
0035 % caters for 3D data (time * channnels * trials);
0036 % - It is important to deemphasize high amplitude artifacts and glitches
0037 % so that they do not dominate the solution.  This implementation uses
0038 % weighted covariance and means.
0039 % - Processing assumes zero-means data. Means are calculated with weights.
0040 % - The implementation tries to be efficent and minimize memory requirements
0041 % so as to handle large data sets.
0042 %
0043 % Larger data sets (disk based) could be handled by performing mean and
0044 % covariance calculations block-by-block, in several passes.
0045 %
0046 
0047 if nargin<2; error('too few arguments'); end
0048 if nargin<3 || isempty(shifts); shifts=0; end
0049 if nargin<4; wx=[]; end
0050 if nargin<5; wref=[]; end
0051 if nargin<6 || isempty(keep); keep=[]; end
0052 if nargin<7 || isempty(thresh); thresh=10.^-20; end
0053 
0054 % check argument values for sanity
0055 if size(x,1)~=size(ref,1); error('X and REF should have same nrows'); end
0056 if size(x,3)~=size(ref,3); error('X and REF should have same npages'); end
0057 if ~isempty(wx) && size(x,3)~=size(wx,3); error('X and WX should have same npages'); end
0058 if ~isempty(wx) && size(x,1)~=size(wx,1); error('X and WX should have same nrows'); end
0059 if ~isempty(wref) && size(ref,1)~=size(wref,1); error('REF and WREF should have same nrows'); end
0060 if ~isempty(wref) && size(ref,3)~=size(wref,3); error('REF and WREF should have same npages'); end
0061 if max(shifts)-min(0,min(shifts)) >= size(x,1); error('X has too few samples to support SHIFTS'); end
0062 if ~isempty(wx) && size(wx,2)~=1; error('wx should have ncols=1'); end
0063 if ~isempty(wref) && size(wref,2)~=1; error('wref should have ncols=1'); end
0064 if ~isempty(wx) && sum(wx(:))==0; error('weights on x are all zero!'); end
0065 if ~isempty(wref) && sum(wref(:))==0; error('weights on ref are all zero!'); end
0066 
0067 if numel(shifts)>1000; error(['numel(shifts)=',num2str(numel(shifts)), ' (if OK comment out this line)']); end
0068 
0069 
0070 % We need to adjust x and ref to ensure that shifts are non-negative.
0071 % If some values of shifts are negative, we increment shifts and truncate x.
0072 
0073 % adjust x to make shifts non-negative
0074 offset1=max(0,-min(shifts));
0075 idx=1+offset1:size(x,1);
0076 x=x(idx,:,:);                             % truncate x
0077 if ~isempty(wx); wx=wx(idx,:,:); end
0078 shifts=shifts+offset1;                    % shifts are now positive
0079 
0080 % adjust size of x
0081 offset2=max(0,max(shifts)); 
0082 idx=1: size(ref,1)-offset2; 
0083 x=x(idx,:,:);                           % part of x that overlaps with time-shifted refs
0084 if ~isempty(wx); wx=wx(idx,:,:); end
0085 
0086 [mx,nx,ox]=size(x);
0087 [mref,nref,oref]=size(ref);
0088 
0089 % consolidate weights into single weight matrix
0090 w=zeros([mx,1,oref]);
0091 if isempty(wx) && isempty(wref)
0092     w(1:mx,:,:)=1;
0093 elseif isempty(wref);
0094     w(:,:,:)=wx(:,:,:);
0095 elseif isempty(wx)
0096     for k=1:ox
0097         wr=wref(:,:,k);
0098         wr=ts_multishift(wr,shifts);
0099         wr=min(wr,[],2);
0100         w(:,:,k)=wr;
0101     end;
0102 else
0103     for k=1:ox
0104         wr=wref(:,:,k);
0105         wr=nt_multishift(wr,shifts);
0106         wr=min(wr,[],2);
0107         wr=min(wr,wx(1:size(wr,1),:,k));
0108         w(:,:,k)=wr;
0109     end
0110 end
0111 wx=w;
0112 wref=zeros(mref,1,oref);
0113 wref(idx,:,:)=w;
0114 
0115 % remove weighted means
0116 x0=x;
0117 x=nt_demean(x,wx);
0118 mn1=x-x0;
0119 ref=nt_demean(ref,wref);
0120 
0121 % equalize power of ref chans, then equalize power of ref PCs
0122 ref=nt_normcol(ref,wref);
0123 ref=nt_pca(ref,0,[],10^-6);
0124 ref=nt_normcol(ref,wref);
0125 
0126 % covariances and cross covariance with time-shifted refs
0127 [cref,twcref]=nt_cov(ref,shifts,wref);
0128 [cxref,twcxref]=nt_xcov(x,ref,shifts,wx);
0129 
0130 % regression matrix of x on time-shifted refs
0131 r=nt_regcov(cxref/twcxref,cref/twcref,keep,thresh);
0132 
0133 %r=r*0.765;
0134 
0135 % TSPCA: clean x by removing regression on time-shifted refs
0136 y=zeros(mx,nx,ox);
0137 for k=1:ox
0138     z=nt_multishift(ref(:,:,k),shifts)*r;
0139     %plot([x(1:size(z,1),1,k), z(:,1)]); pause
0140     y(:,:,k)=x(1:size(z,1),:,k)-z;
0141 end
0142 y0=y;
0143 y=nt_demean(y,wx);    % multishift(ref) is not necessarily 0 mean
0144 mn2=y-y0;
0145 
0146 %idx=1+offset1:n0-offset2;
0147 idx=1+offset1:size(y,1)+offset1;
0148 mn=mn1+mn2;
0149 w=wref;
0150 
0151 %y=vecadd(y,mn);