1function [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT) 2% XCOVF generates cross-covariance function. 3% XCOVF is the same as XCORR except 4% X and Y can contain missing values encoded with NaN. 5% NaN's are skipped, NaN do not result in a NaN output. 6% The output gives NaN only if there are insufficient input data 7% 8% [C,N,LAGS] = xcovf(X,MAXLAG,SCALEOPT); 9% calculates the (auto-)correlation function of X 10% [C,N,LAGS] = xcovf(X,Y,MAXLAG,SCALEOPT); 11% calculates the crosscorrelation function between X and Y 12% 13% SCALEOPT [character string] specifies the type of scaling applied 14% to the correlation vector (or matrix). is one of: 15% 'none' return the unscaled correlation, R, 16% 'biased' return the biased average, R/N, 17% 'unbiased' return the unbiassed average, R(k)/(N-|k|), 18% 'coeff' return the correlation coefficient, R/(rms(x).rms(y)), 19% where "k" is the lag, and "N" is the length of X. 20% If omitted, the default value is "none". 21% If Y is supplied but does not have the ame length as X, 22% scale must be "none". 23% 24% 25% see also: COVM, XCORR 26 27% Copyright (C) 2005,2010,2011 by Alois Schloegl <alois.schloegl@gmail.com> 28% This function is part of the NaN-toolbox 29% http://pub.ist.ac.at/~schloegl/matlab/NaN/ 30 31% This program is free software; you can redistribute it and/or modify 32% it under the terms of the GNU General Public License as published by 33% the Free Software Foundation; either version 3 of the License, or 34% (at your option) any later version. 35% 36% This program is distributed in the hope that it will be useful, 37% but WITHOUT ANY WARRANTY; without even the implied warranty of 38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39% GNU General Public License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with this program; If not, see <http://www.gnu.org/licenses/>. 43 44if nargin<2, 45 Y = []; 46 MAXLAG = []; 47 SCALEOPT = 'none'; 48elseif ischar(Y), 49 MAXLAG = Y; 50 SCALEOPT=MAXLAG; 51 Y=[]; 52elseif all(size(Y)==1), 53 if nargin<3 54 SCALEOPT = 'none'; 55 else 56 SCALEOPT = MAXLAG; 57 end; 58 MAXLAG = Y; 59 Y = []; 60end; 61 62if 0, 63 64elseif isempty(Y) && isempty(MAXLAG) 65 NX = isnan(X); 66 X(NX) = 0; 67 [C,LAGS] = xcorr(X,'none'); 68 [N,LAGS] = xcorr(1-NX,'none'); 69elseif ~isempty(Y) && isempty(MAXLAG) 70 NX = isnan(X); 71 NY = isnan(Y); 72 X(NX) = 0; 73 Y(NY) = 0; 74 [C,LAGS] = xcorr(X,Y,'none'); 75 [N,LAGS] = xcorr(1-NX,1-NY,'none'); 76elseif isempty(Y) && ~isempty(MAXLAG) 77 NX = isnan(X); 78 X(NX) = 0; 79 [C,LAGS] = xcorr(X,MAXLAG,'none'); 80 [N,LAGS] = xcorr(1-NX,MAXLAG,'none'); 81elseif ~isempty(Y) && ~isempty(MAXLAG) 82 NX = isnan(X); 83 NY = isnan(Y); 84 X(NX) = 0; 85 Y(NY) = 0; 86 [C,LAGS] = xcorr(X,Y,MAXLAG,'none'); 87 [N,LAGS] = xcorr(1-NX,1-NY,MAXLAG,'none'); 88end; 89 90if 0, 91 92elseif strcmp(SCALEOPT,'none') 93 % done 94 95elseif strcmp(SCALEOPT,'coeff') 96 ix = find(LAGS==0); 97 if ~any(size(X)==1), %% ~isvector(X) 98 c = C(ix,1:size(X,2)+1:end); %% diagonal elements 99 v = c.^-0.5; % sqrt(1./c(:)); 100 v = v'*v; 101 C = C.*repmat(v(:).',size(C,1),1); 102 elseif isempty(Y) 103 C = C/C(ix); 104 else 105 C = C/sqrt(sumsq(X)*sumsq(Y)); 106 end; 107 108elseif strcmp(SCALEOPT,'biased') 109 C = C./repmat(max(N),size(C,1),1); 110 111elseif strcmp(SCALEOPT,'unbiased') 112 C = C./(repmat(max(N),size(C,1),1)-repmat(LAGS,1,size(C,2))); 113 114else 115 warning('invalid SCALEOPT - not supported'); 116end; 117 118