[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. 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. 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 if nargin<4||isempty(demean_flag); demean_flag=0; end 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 0030 if min(shifts)<0; error('shifts should be non-negative'); end 0031 0032 shifts=shifts(:); % --> column vector 0033 nshifts=numel(shifts); 0034 0035 [m,n,o]=size(x); 0036 c=zeros(n*nshifts); 0037 0038 if isempty(w) 0039 % no weights 0040 0041 for k=1:o 0042 xx=nt_multishift(x(:,:,k),shifts); 0043 if demean_flag 0044 xx=nt_demean(xx); 0045 end 0046 c=c+xx'*xx; 0047 end 0048 tw=size(xx,1)*o; 0049 0050 else 0051 % weights 0052 if size(w,2)>1; error('w should have single column'); end 0053 0054 for k=1:o 0055 if ~all(shifts == [0]); 0056 xx=nt_multishift(x(:,:,k),shifts); 0057 ww=nt_multishift(w(:,:,k),shifts); 0058 ww=min(ww,[],2); 0059 else 0060 xx=x(:,:,k); ww=w(:,:,k); 0061 end 0062 xx=nt_vecmult(xx,ww); 0063 if demean_flag; xx=nt_demean(xx); end 0064 c=c+xx'*xx; 0065 end 0066 tw=sum(w(:)); 0067 end 0068