1function [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 6% 7% COVM(X,Mode); 8% calculates the (auto-)correlation matrix of X 9% COVM(X,Y,Mode); 10% calculates the crosscorrelation between X and Y 11% COVM(...,W); 12% weighted crosscorrelation 13% 14% Mode = 'M' minimum or standard mode [default] 15% C = X'*X; or X'*Y correlation matrix 16% 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 21% 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. 27% 28% C = covm(...); 29% C is the scaled by N in Mode M and by (N-1) in mode D. 30% [C,N] = covm(...); 31% C is not scaled, provides the scaling factor N 32% C./N gives the scaled version. 33% 34% see also: DECOVM, XCOVF 35 36% Copyright (C) 2000-2005,2009 by Alois Schloegl <alois.schloegl@gmail.com> 37% This function is part of the NaN-toolbox 38% http://pub.ist.ac.at/~schloegl/matlab/NaN/ 39 40% This program is free software; you can redistribute it and/or modify 41% it under the terms of the GNU General Public License as published by 42% the Free Software Foundation; either version 3 of the License, or 43% (at your option) any later version. 44% 45% This program is distributed in the hope that it will be useful, 46% but WITHOUT ANY WARRANTY; without even the implied warranty of 47% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 48% GNU General Public License for more details. 49% 50% You should have received a copy of the GNU General Public License 51% along with this program; If not, see <http://www.gnu.org/licenses/>. 52 53 54global FLAG_NANS_OCCURED; 55 56if nargin<3, 57 W = []; 58 if nargin==2, 59 if isnumeric(Y), 60 Mode='M'; 61 else 62 Mode=Y; 63 Y=[]; 64 end; 65 elseif nargin==1, 66 Mode = 'M'; 67 Y = []; 68 elseif nargin==0, 69 error('Missing argument(s)'); 70 end; 71 72elseif (nargin==3) && isnumeric(Y) && ~isnumeric(Mode); 73 W = []; 74 75elseif (nargin==3) && ~isnumeric(Y) && isnumeric(Mode); 76 W = Mode; 77 Mode = Y; 78 Y = []; 79 80elseif (nargin==4) && ~isnumeric(Mode) && isnumeric(Y); 81 ; %% ok 82else 83 error('invalid input arguments'); 84end; 85 86Mode = upper(Mode); 87 88[r1,c1]=size(X); 89if ~isempty(Y) 90 [r2,c2]=size(Y); 91 if r1~=r2, 92 error('X and Y must have the same number of observations (rows).'); 93 end; 94else 95 [r2,c2]=size(X); 96end; 97 98persistent mexFLAG2; 99persistent mexFLAG; 100if isempty(mexFLAG2) 101 mexFLAG2 = exist('covm_mex','file'); 102end; 103if isempty(mexFLAG) 104 mexFLAG = exist('sumskipnan_mex','file'); 105end; 106 107 108if ~isempty(W) 109 W = W(:); 110 if (r1~=numel(W)) 111 error('Error COVM: size of weight vector does not fit number of rows'); 112 end; 113 %w = spdiags(W(:),0,numel(W),numel(W)); 114 %nn = sum(W(:)); 115 nn = sum(W); 116else 117 nn = r1; 118end; 119 120 121if mexFLAG2 && mexFLAG && ~isempty(W), 122 %% the mex-functions here are much slower than the m-scripts below 123 %% however, the mex-functions support weighting of samples. 124 if isempty(FLAG_NANS_OCCURED), 125 %% mex-files require that FLAG_NANS_OCCURED is not empty, 126 %% otherwise, the status of NAN occurence can not be returned. 127 FLAG_NANS_OCCURED = logical(0); % default value 128 end; 129 130 if any(Mode=='D') || any(Mode=='E'), 131 [S1,N1] = sumskipnan(X,1,W); 132 if ~isempty(Y) 133 [S2,N2] = sumskipnan(Y,1,W); 134 else 135 S2 = S1; N2 = N1; 136 end; 137 if any(Mode=='D'), % detrending mode 138 X = X - ones(r1,1)*(S1./N1); 139 if ~isempty(Y) 140 Y = Y - ones(r1,1)*(S2./N2); 141 end; 142 end; 143 end; 144 145 if issparse(X) || issparse(Y), 146 fprintf(2,'sumskipnan: sparse matrix converted to full matrix\n'); 147 X=full(X); 148 Y=full(Y); 149 end; 150 151 [CC,NN] = covm_mex(real(X), real(Y), FLAG_NANS_OCCURED, W); 152 %% complex matrices 153 if ~isreal(X) && ~isreal(Y) 154 [iCC,inn] = covm_mex(imag(X), imag(Y), FLAG_NANS_OCCURED, W); 155 CC = CC + iCC; 156 end; 157 if isempty(Y) Y = X; end; 158 if ~isreal(X) 159 [iCC,inn] = covm_mex(imag(X), real(Y), FLAG_NANS_OCCURED, W); 160 CC = CC - i*iCC; 161 end; 162 if ~isreal(Y) 163 [iCC,inn] = covm_mex(real(X), imag(Y), FLAG_NANS_OCCURED, W); 164 CC = CC + i*iCC; 165 end; 166 167 if any(Mode=='D') && ~any(Mode=='1'), % 'D1' 168 NN = max(NN-1,0); 169 end; 170 if any(Mode=='E'), % extended mode 171 NN = [nn, N2; N1', NN]; 172 CC = [nn, S2; S1', CC]; 173 end; 174 175 176elseif ~isempty(W), 177 178 error('Error COVM: weighted COVM requires sumskipnan_mex and covm_mex but it is not available'); 179 180 %% weighted covm without mex-file support 181 %% this part is not working. 182 183elseif ~isempty(Y), 184 if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M 185 NN = real(X==X)'*real(Y==Y); 186 FLAG_NANS_OCCURED = any(NN(:)<nn); 187 X(X~=X) = 0; % skip NaN's 188 Y(Y~=Y) = 0; % skip NaN's 189 CC = X'*Y; 190 191 else % if any(Mode=='D') | any(Mode=='E'), 192 [S1,N1] = sumskipnan(X,1); 193 [S2,N2] = sumskipnan(Y,1); 194 NN = real(X==X)'*real(Y==Y); 195 196 if any(Mode=='D'), % detrending mode 197 X = X - ones(r1,1)*(S1./N1); 198 Y = Y - ones(r1,1)*(S2./N2); 199 if any(Mode=='1'), % 'D1' 200 NN = NN; 201 else % 'D0' 202 NN = max(NN-1,0); 203 end; 204 end; 205 X(X~=X) = 0; % skip NaN's 206 Y(Y~=Y) = 0; % skip NaN's 207 CC = X'*Y; 208 209 if any(Mode=='E'), % extended mode 210 NN = [nn, N2; N1', NN]; 211 CC = [nn, S2; S1', CC]; 212 end; 213 end; 214 215else 216 if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M 217 tmp = real(X==X); 218 NN = tmp'*tmp; 219 X(X~=X) = 0; % skip NaN's 220 CC = X'*X; 221 FLAG_NANS_OCCURED = any(NN(:)<nn); 222 223 else % if any(Mode=='D') | any(Mode=='E'), 224 [S,N] = sumskipnan(X,1); 225 tmp = real(X==X); 226 NN = tmp'*tmp; 227 if any(Mode=='D'), % detrending mode 228 X = X - ones(r1,1)*(S./N); 229 if any(Mode=='1'), % 'D1' 230 NN = NN; 231 else % 'D0' 232 NN = max(NN-1,0); 233 end; 234 end; 235 236 X(X~=X) = 0; % skip NaN's 237 CC = X'*X; 238 if any(Mode=='E'), % extended mode 239 NN = [nn, N; N', NN]; 240 CC = [nn, S; S', CC]; 241 end; 242 end 243 244end; 245 246 247if nargout<2 248 CC = CC./NN; % unbiased 249end; 250 251%!assert(covm([1;NaN;2],'D'),0.5) 252%!assert(covm([1;NaN;2],'M'),2.5) 253%!assert(covm([1;NaN;2],'E'),[1,1.5;1.5,2.5]) 254