nt_cov

[c,tw]=nt_cov(x,shifts,w) - time shift covariance

  c: covariance matrix
  tw: total weight (c/tw is normalized covariance)

  x: data
  shifts: array of time shifts (must be non-negative)
  w: weights
  
 If SHIFTS==[0], this function calculates the covariance matrix of X, 
 weighted by W.  In the general case it calculates the covariance matrix
 of time-shifted X.

 X can be 1D, 2D or 3D.  It can also be a cell array. 
 W can be 1D (if X is 1D or 2D) or 2D (if X is 3D). The same weight is
 applied to each column.
 
 Output is a 2D matrix with dimensions (ncols(X)*numel(SHIFTS))^2.
 It is made up of an ncols(X)*ncols(X) matrix of submatrices, each of 
 dimensions numel(SHIFTS)*numel(SHIFTS).

 NoiseTools

0001 function [c,tw]=nt_cov(x,shifts,w);
0002 %[c,tw]=nt_cov(x,shifts,w) - time shift covariance
0003 %
0004 %  c: covariance matrix
0005 %  tw: total weight (c/tw is normalized covariance)
0006 %
0007 %  x: data
0008 %  shifts: array of time shifts (must be non-negative)
0009 %  w: weights
0010 %
0011 % If SHIFTS==[0], this function calculates the covariance matrix of X,
0012 % weighted by W.  In the general case it calculates the covariance matrix
0013 % of time-shifted X.
0014 %
0015 % X can be 1D, 2D or 3D.  It can also be a cell array.
0016 % W can be 1D (if X is 1D or 2D) or 2D (if X is 3D). The same weight is
0017 % applied to each column.
0018 %
0019 % Output is a 2D matrix with dimensions (ncols(X)*numel(SHIFTS))^2.
0020 % It is made up of an ncols(X)*ncols(X) matrix of submatrices, each of
0021 % dimensions numel(SHIFTS)*numel(SHIFTS).
0022 %
0023 % NoiseTools
0024 
0025 %% arguments
0026 if nargin<3; w=[]; end;
0027 if nargin<2||isempty(shifts); shifts=0; end;
0028 if prod(size(x))==0; error('data empty'); end
0029 shifts=shifts(:);           % --> column vector
0030 nshifts=numel(shifts); 
0031 
0032 if isempty(x); error('data empty'); end
0033 
0034 if isempty(w) 
0035     %% no weights
0036     if isnumeric(x)
0037         %% matrix
0038         [m,n,o]=size(x);
0039         c=zeros(n*nshifts);
0040         for k=1:o
0041             xx=nt_multishift(x(:,:,k),shifts);
0042             c=c+xx'*xx;
0043         end
0044         tw=size(xx,1)*o;
0045     elseif iscell(x)
0046         %% cell array
0047         [m,n]=size(x{1});
0048         o=length(x);
0049         c=zeros(n*nshifts);
0050         for k=1:o;
0051             if size(x{k},2)~=n;  disp([size(x{k})]); disp([n, k]); error('!'); end
0052             xx=nt_multishift(x{k}(:,:),shifts);
0053             c=c+xx'*xx;
0054         end
0055         tw=size(xx,1)*o;
0056     else error('!'); end   
0057 else
0058     %% weights
0059     if isnumeric(x)
0060         %% matrix
0061         [m,n,o]=size(x);
0062         c=zeros(n*nshifts);
0063         for k=1:o
0064             if ~all(shifts == [0]); 
0065                 xx=nt_multishift(x(:,:,k),shifts); 
0066                 ww=nt_multishift(w(:,:,k),shifts);
0067                 ww=min(ww,[],2);
0068             else
0069                 xx=x(:,:,k); ww=w(:,:,k);
0070             end
0071             xx=nt_vecmult(xx,ww);
0072             c=c+xx'*xx;
0073         end
0074         tw=sum(w(:));
0075     else
0076         %% cell array
0077         c=zeros(n*nshifts);
0078         [m,n]=size(x{1});
0079         o=length(x);
0080         for k=1:o
0081             if ~all(shifts == [0]); 
0082                 xx=nt_multishift(x{k}(:,:),shifts); 
0083                 ww=nt_multishift(w{k}(:,:),shifts);
0084                 ww=min(ww,[],2);
0085             else
0086                 xx=x{k}(:,:); ww=w{k}(:,:);
0087             end
0088             xx=nt_vecmult(xx,ww);
0089             c=c+xx'*xx;
0090         end
0091         tw=sum(w(:));
0092     end       
0093 end
0094