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