1 function [CC,NN] = covm(X,Y,Mode,W)
2 % COVM generates covariance matrix
3 % X and Y can contain missing values encoded with NaN.
4 % NaN's are skipped, NaN do not result in a NaN output.
5 % The output gives NaN only if there are insufficient input data
8 % calculates the (auto-)correlation matrix of X
10 % calculates the crosscorrelation between X and Y
12 % weighted crosscorrelation
14 % Mode = 'M' minimum or standard mode [default]
15 % C = X'*X; or X'*Y correlation matrix
17 % Mode = 'E' extended mode
18 % C = [1 X]'*[1 X]; % l is a matching column of 1's
19 % C is additive, i.e. it can be applied to subsequent blocks and summed up afterwards
20 % the mean (or sum) is stored on the 1st row and column of C
22 % Mode = 'D' or 'D0' detrended mode
23 % the mean of X (and Y) is removed. If combined with extended mode (Mode='DE'),
24 % the mean (or sum) is stored in the 1st row and column of C.
25 % The default scaling is factor (N-1).
26 % Mode = 'D1' is the same as 'D' but uses N for scaling.
29 % C is the scaled by N in Mode M and by (N-1) in mode D.
31 % C is not scaled, provides the scaling factor N
32 % C./N gives the scaled version.
34 % see also: DECOVM, XCOVF
36 % $Id: covm.m 9032 2011-11-08 20:25:36Z schloegl $
37 % Copyright (C) 2000-2005,2009 by Alois Schloegl <alois.schloegl@gmail.com>
38 % This function is part of the NaN-toolbox
39 % http://pub.ist.ac.at/~schloegl/matlab/NaN/
41 % This program is free software; you can redistribute it and/or modify
42 % it under the terms of the GNU General Public License as published by
43 % the Free Software Foundation; either version 3 of the License, or
44 % (at your option) any later version.
46 % This program is distributed in the hope that it will be useful,
47 % but WITHOUT ANY WARRANTY; without even the implied warranty of
48 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
49 % GNU General Public License for more details.
51 % You should have received a copy of the GNU General Public License
52 % along with this program; If not, see <http://www.gnu.org/licenses/>.
55 global FLAG_NANS_OCCURED;
70 error('Missing argument(s)');
73 elseif (nargin==3) && isnumeric(Y) && ~isnumeric(Mode);
76 elseif (nargin==3) && ~isnumeric(Y) && isnumeric(Mode);
81 elseif (nargin==4) && ~isnumeric(Mode) && isnumeric(Y);
84 error('invalid input arguments');
93 error('X and Y must have the same number of observations (rows).');
102 mexFLAG2 = exist('covm_mex','file');
105 mexFLAG = exist('sumskipnan_mex','file');
112 error('Error COVM: size of weight vector does not fit number of rows');
114 %w = spdiags(W(:),0,numel(W),numel(W));
122 if mexFLAG2 && mexFLAG && ~isempty(W),
123 %% the mex-functions here are much slower than the m-scripts below
124 %% however, the mex-functions support weighting of samples.
125 if isempty(FLAG_NANS_OCCURED),
126 %% mex-files require that FLAG_NANS_OCCURED is not empty,
127 %% otherwise, the status of NAN occurence can not be returned.
128 FLAG_NANS_OCCURED = logical(0); % default value
131 if any(Mode=='D') || any(Mode=='E'),
132 [S1,N1] = sumskipnan(X,1,W);
134 [S2,N2] = sumskipnan(Y,1,W);
138 if any(Mode=='D'), % detrending mode
139 X = X - ones(r1,1)*(S1./N1);
141 Y = Y - ones(r1,1)*(S2./N2);
146 [CC,NN] = covm_mex(real(X), real(Y), FLAG_NANS_OCCURED, W);
148 if ~isreal(X) && ~isreal(Y)
149 [iCC,inn] = covm_mex(imag(X), imag(Y), FLAG_NANS_OCCURED, W);
152 if isempty(Y) Y = X; end;
154 [iCC,inn] = covm_mex(imag(X), real(Y), FLAG_NANS_OCCURED, W);
158 [iCC,inn] = covm_mex(real(X), imag(Y), FLAG_NANS_OCCURED, W);
162 if any(Mode=='D') && ~any(Mode=='1'), % 'D1'
165 if any(Mode=='E'), % extended mode
166 NN = [nn, N2; N1', NN];
167 CC = [nn, S2; S1', CC];
173 error('Error COVM: weighted COVM requires sumskipnan_mex and covm_mex but it is not available');
175 %% weighted covm without mex-file support
176 %% this part is not working.
179 if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
180 NN = real(X==X)'*real(Y==Y);
181 FLAG_NANS_OCCURED = any(NN(:)<nn);
182 X(X~=X) = 0; % skip NaN's
183 Y(Y~=Y) = 0; % skip NaN's
186 else % if any(Mode=='D') | any(Mode=='E'),
187 [S1,N1] = sumskipnan(X,1);
188 [S2,N2] = sumskipnan(Y,1);
189 NN = real(X==X)'*real(Y==Y);
191 if any(Mode=='D'), % detrending mode
192 X = X - ones(r1,1)*(S1./N1);
193 Y = Y - ones(r1,1)*(S2./N2);
194 if any(Mode=='1'), % 'D1'
200 X(X~=X) = 0; % skip NaN's
201 Y(Y~=Y) = 0; % skip NaN's
204 if any(Mode=='E'), % extended mode
205 NN = [nn, N2; N1', NN];
206 CC = [nn, S2; S1', CC];
211 if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
214 X(X~=X) = 0; % skip NaN's
216 FLAG_NANS_OCCURED = any(NN(:)<nn);
218 else % if any(Mode=='D') | any(Mode=='E'),
219 [S,N] = sumskipnan(X,1);
222 if any(Mode=='D'), % detrending mode
223 X = X - ones(r1,1)*(S./N);
224 if any(Mode=='1'), % 'D1'
231 X(X~=X) = 0; % skip NaN's
233 if any(Mode=='E'), % extended mode
234 NN = [nn, N; N', NN];
235 CC = [nn, S; S', CC];
243 CC = CC./NN; % unbiased
246 %!assert(covm([1;NaN;2],'D'),0.5)
247 %!assert(covm([1;NaN;2],'M'),2.5)
248 %!assert(covm([1;NaN;2],'E'),[1,1.5;1.5,2.5])