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