Home > NoiseTools > nt_tsr.m

# nt_tsr

## PURPOSE

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

## SYNOPSIS

function [y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh)

## DESCRIPTION

    [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

## 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_greetings nt_greetings - display message the first time the toolbox is used
- nt_multishift z=nt_multishift(x,shifts,amplitudes) - apply multiple shifts to matrix
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1
- nt_pca [z,idx]=nt_pca(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_regcov r=nt_regcov(cxy,cyy,keep,threshold) - regression matrix from cross covariance
- nt_xcov [c,tw]=nt_xcov(x,y,shifts,w) - cross-covariance of X and time-shifted Y

This function is called by:

- nt_sns1 y=nt_sns1(x,nneigbors,skip,w,threshold) - sensor noise suppression
- nt_zapline [y,yy]=nt_zapline(x,fline,nremove,p,plotflag) - remove power line artifact

## SOURCE CODE

    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);