nt_tsxcov

[c,tw]=nt_tsxcov(x,y,shifts,w) - cross-covariance of X and time-shifted Y


  c: cross-covariance matrix
  tw: total weight

  x,y: data to cross correlate
  shifts: array of time shifts (must be non-negative)
  w: weights
  
 This function calculates, for each pair of columns (Xi,Yj) of X and Y, the
 scalar products between Xi and time-shifted versions of Yj. 
 Shifts are taken from array SHIFTS. 

 The weights are applied to X.

 X can be 1D, 2D or 3D.  W is 1D (if X is 1D or 2D) or 2D (if X is 3D).
 
 Output is a 2D matrix with dimensions ncols(X)*(ncols(Y)*nshifts).
 
 NoiseTools

0001 function [c,tw]=nt_tsxcov(x,y,shifts,w);
0002 %[c,tw]=nt_tsxcov(x,y,shifts,w) - cross-covariance of X and time-shifted Y
0003 %
0004 %
0005 %  c: cross-covariance matrix
0006 %  tw: total weight
0007 %
0008 %  x,y: data to cross correlate
0009 %  shifts: array of time shifts (must be non-negative)
0010 %  w: weights
0011 %
0012 % This function calculates, for each pair of columns (Xi,Yj) of X and Y, the
0013 % scalar products between Xi and time-shifted versions of Yj.
0014 % Shifts are taken from array SHIFTS.
0015 %
0016 % The weights are applied to X.
0017 %
0018 % X can be 1D, 2D or 3D.  W is 1D (if X is 1D or 2D) or 2D (if X is 3D).
0019 %
0020 % Output is a 2D matrix with dimensions ncols(X)*(ncols(Y)*nshifts).
0021 %
0022 % NoiseTools
0023 
0024 if nargin<4; w=[]; end;
0025 if nargin<3||isempty(shifts); shifts=0; end;
0026 
0027 if ~isempty(w) && size(x,1)~=size(w,1); error('X and W should have same nrows'); end
0028 if size(x,3)~=size(y,3); error('X and Y should have same npages'); end
0029 if ~isempty(w) && size(x,3)~=size(w,3); error('X and W should have same npages'); end
0030 
0031 shifts=shifts(:); 
0032 nshifts=numel(shifts); 
0033 
0034 [mx,nx,ox]=size(x);
0035 [my,ny,oy]=size(y);
0036 c=zeros(nx,ny*nshifts);
0037 
0038 if ~isempty(w)
0039     x=nt_fold(nt_vecmult(nt_unfold(x),nt_unfold(w)),mx);
0040 end
0041 
0042 % cross covariance
0043 for k=1:ox
0044     yy=nt_multishift(y(:,:,k),shifts);
0045     xx=x(1:size(yy,1),:,k);
0046     c=c+xx'*yy;
0047 end
0048 
0049 if isempty(w)
0050     tw=ox*ny*size(yy,1);
0051 else
0052     w=w(1:size(yy,1),:,:);
0053     tw=sum(w(:));
0054 end