1## Copyright (C) 1999-2001 Paul Kienzle <pkienzle@users.sf.net> 2## Copyright (C) 2004 <asbjorn.sabo@broadpark.no> 3## Copyright (C) 2008, 2010 Peter Lanspeary <peter.lanspeary@.adelaide.edu.au> 4## 5## This program is free software: you can redistribute it and/or modify 6## it under the terms of the GNU General Public License as published by 7## the Free Software Foundation, either version 3 of the License, or 8## (at your option) any later version. 9## 10## This program is distributed in the hope that it will be useful, 11## but WITHOUT ANY WARRANTY; without even the implied warranty of 12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13## GNU General Public License for more details. 14## 15## You should have received a copy of the GNU General Public License 16## along with this program; see the file COPYING. If not, see 17## <https://www.gnu.org/licenses/>. 18 19## -*- texinfo -*- 20## @deftypefn {Function File} {[@var{R}, @var{lag}] =} xcorr ( @var{X} ) 21## @deftypefnx {Function File} {@dots{} =} xcorr ( @var{X}, @var{Y} ) 22## @deftypefnx {Function File} {@dots{} =} xcorr ( @dots{}, @var{maxlag}) 23## @deftypefnx {Function File} {@dots{} =} xcorr ( @dots{}, @var{scale}) 24## Estimates the cross-correlation. 25## 26## Estimate the cross correlation R_xy(k) of vector arguments @var{X} and @var{Y} 27## or, if @var{Y} is omitted, estimate autocorrelation R_xx(k) of vector @var{X}, 28## for a range of lags k specified by argument "maxlag". If @var{X} is a 29## matrix, each column of @var{X} is correlated with itself and every other 30## column. 31## 32## The cross-correlation estimate between vectors "x" and "y" (of 33## length N) for lag "k" is given by 34## 35## @tex 36## $$ R_{xy}(k) = \sum_{i=1}^{N} x_{i+k} \conj(y_i), 37## @end tex 38## @ifnottex 39## @example 40## @group 41## N 42## R_xy(k) = sum x_@{i+k@} conj(y_i), 43## i=1 44## @end group 45## @end example 46## @end ifnottex 47## 48## where data not provided (for example x(-1), y(N+1)) is zero. Note the 49## definition of cross-correlation given above. To compute a 50## cross-correlation consistent with the field of statistics, see @command{xcov}. 51## 52## @strong{ARGUMENTS} 53## @table @var 54## @item X 55## [non-empty; real or complex; vector or matrix] data 56## 57## @item Y 58## [real or complex vector] data 59## 60## If @var{X} is a matrix (not a vector), @var{Y} must be omitted. 61## @var{Y} may be omitted if @var{X} is a vector; in this case xcorr 62## estimates the autocorrelation of @var{X}. 63## 64## @item maxlag 65## [integer scalar] maximum correlation lag 66## If omitted, the default value is N-1, where N is the 67## greater of the lengths of @var{X} and @var{Y} or, if @var{X} is a matrix, 68## the number of rows in @var{X}. 69## 70## @item scale 71## [character string] specifies the type of scaling applied 72## to the correlation vector (or matrix). is one of: 73## @table @samp 74## @item none 75## return the unscaled correlation, R, 76## @item biased 77## return the biased average, R/N, 78## @item unbiased 79## return the unbiased average, R(k)/(N-|k|), 80## @item coeff 81## return the correlation coefficient, R/(rms(x).rms(y)), 82## where "k" is the lag, and "N" is the length of @var{X}. 83## If omitted, the default value is "none". 84## If @var{Y} is supplied but does not have the same length as @var{X}, 85## scale must be "none". 86## @end table 87## @end table 88## 89## @strong{RETURNED VARIABLES} 90## @table @var 91## @item R 92## array of correlation estimates 93## @item lag 94## row vector of correlation lags [-maxlag:maxlag] 95## @end table 96## 97## The array of correlation estimates has one of the following forms: 98## (1) Cross-correlation estimate if @var{X} and @var{Y} are vectors. 99## 100## (2) Autocorrelation estimate if is a vector and @var{Y} is omitted. 101## 102## (3) If @var{X} is a matrix, R is an matrix containing the cross-correlation 103## estimate of each column with every other column. Lag varies with the first 104## index so that R has 2*maxlag+1 rows and P^2 columns where P is the number of 105## columns in @var{X}. 106## 107## If Rij(k) is the correlation between columns i and j of @var{X} 108## 109## @code{R(k+maxlag+1,P*(i-1)+j) == Rij(k)} 110## 111## for lag k in [-maxlag:maxlag], or 112## 113## @code{R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j))}. 114## 115## @code{reshape(R(k,:),P,P)} is the cross-correlation matrix for @code{X(k,:)}. 116## 117## @seealso{xcov} 118## @end deftypefn 119 120## The cross-correlation estimate is calculated by a "spectral" method 121## in which the FFT of the first vector is multiplied element-by-element 122## with the FFT of second vector. The computational effort depends on 123## the length N of the vectors and is independent of the number of lags 124## requested. If you only need a few lags, the "direct sum" method may 125## be faster: 126## 127## Ref: Stearns, SD and David, RA (1988). Signal Processing Algorithms. 128## New Jersey: Prentice-Hall. 129 130## unbiased: 131## ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ [N:-1:N-k+1](:) 132## biased: 133## ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ N 134## 135## If length(x) == length(y) + k, then you can use the simpler 136## ( hankel(x(1:k),x(k:N-k)) * y ) ./ N 137 138function [R, lags] = xcorr (X, Y, maxlag, scale) 139 140 if (nargin < 1 || nargin > 4) 141 print_usage; 142 endif 143 144 ## assign arguments that are missing from the list 145 ## or reassign (right shift) them according to data type 146 if nargin==1 147 Y=[]; maxlag=[]; scale=[]; 148 elseif nargin==2 149 maxlag=[]; scale=[]; 150 if ischar(Y), scale=Y; Y=[]; 151 elseif isscalar(Y), maxlag=Y; Y=[]; 152 endif 153 elseif nargin==3 154 scale=[]; 155 if ischar(maxlag), scale=maxlag; maxlag=[]; endif 156 if isscalar(Y), maxlag=Y; Y=[]; endif 157 endif 158 159 ## assign defaults to missing arguments 160 if isvector(X) 161 ## if isempty(Y), Y=X; endif ## this line disables code for autocorr'n 162 N = max(length(X),length(Y)); 163 else 164 N = rows(X); 165 endif 166 if isempty(maxlag), maxlag=N-1; endif 167 if isempty(scale), scale='none'; endif 168 169 ## check argument values 170 if isempty(X) || isscalar(X) || ischar(Y) || ! ismatrix(X) 171 error("xcorr: X must be a vector or matrix"); 172 endif 173 if isscalar(Y) || ischar(Y) || (!isempty(Y) && !isvector(Y)) 174 error("xcorr: Y must be a vector"); 175 endif 176 if !isempty(Y) && !isvector(X) 177 error("xcorr: X must be a vector if Y is specified"); 178 endif 179 if !isscalar(maxlag) || !isreal(maxlag) || maxlag<0 || fix(maxlag)!=maxlag 180 error("xcorr: maxlag must be a single non-negative integer"); 181 endif 182 ## 183 ## sanity check on number of requested lags 184 ## Correlations for lags in excess of +/-(N-1) 185 ## (a) are not calculated by the FFT algorithm, 186 ## (b) are all zero; so provide them by padding 187 ## the results (with zeros) before returning. 188 if (maxlag > N-1) 189 pad_result = maxlag - (N - 1); 190 maxlag = N - 1; 191 %error("xcorr: maxlag must be less than length(X)"); 192 else 193 pad_result = 0; 194 endif 195 if isvector(X) && isvector(Y) && length(X) != length(Y) && ... 196 !strcmp(scale,'none') 197 error("xcorr: scale must be 'none' if length(X) != length(Y)") 198 endif 199 200 P = columns(X); 201 M = 2^nextpow2(N + maxlag); 202 if !isvector(X) 203 ## For matrix X, correlate each column "i" with all other "j" columns 204 R = zeros(2*maxlag+1,P^2); 205 206 ## do FFTs of padded column vectors 207 pre = fft (postpad (prepad (X, N+maxlag), M) ); 208 post = conj (fft (postpad (X, M))); 209 210 ## do autocorrelations (each column with itself) 211 ## -- if result R is reshaped to 3D matrix (i.e. R=reshape(R,M,P,P)) 212 ## the autocorrelations are on leading diagonal columns of R, 213 ## where i==j in R(:,i,j) 214 cor = ifft (post .* pre); 215 R(:, 1:P+1:P^2) = cor (1:2*maxlag+1,:); 216 217 ## do the cross correlations 218 ## -- these are the off-diagonal column of the reshaped 3D result 219 ## matrix -- i!=j in R(:,i,j) 220 for i=1:P-1 221 j = i+1:P; 222 cor = ifft( pre(:,i*ones(length(j),1)) .* post(:,j) ); 223 R(:,(i-1)*P+j) = cor(1:2*maxlag+1,:); 224 R(:,(j-1)*P+i) = conj( flipud( cor(1:2*maxlag+1,:) ) ); 225 endfor 226 elseif isempty(Y) 227 ## compute autocorrelation of a single vector 228 post = fft( postpad(X(:),M) ); 229 cor = ifft( post .* conj(post) ); 230 R = [ conj(cor(maxlag+1:-1:2)) ; cor(1:maxlag+1) ]; 231 else 232 ## compute cross-correlation of X and Y 233 ## If one of X & Y is a row vector, the other can be a column vector. 234 pre = fft( postpad( prepad( X(:), length(X)+maxlag ), M) ); 235 post = fft( postpad( Y(:), M ) ); 236 cor = ifft( pre .* conj(post) ); 237 R = cor(1:2*maxlag+1); 238 endif 239 240 ## if inputs are real, outputs should be real, so ignore the 241 ## insignificant complex portion left over from the FFT 242 if isreal(X) && (isempty(Y) || isreal(Y)) 243 R=real(R); 244 endif 245 246 ## correct for bias 247 if strcmp(scale, 'biased') 248 R = R ./ N; 249 elseif strcmp(scale, 'unbiased') 250 R = R ./ ( [ N-maxlag:N-1, N, N-1:-1:N-maxlag ]' * ones(1,columns(R)) ); 251 elseif strcmp(scale, 'coeff') 252 ## R = R ./ R(maxlag+1) works only for autocorrelation 253 ## For cross correlation coeff, divide by rms(X)*rms(Y). 254 if !isvector(X) 255 ## for matrix (more than 1 column) X 256 rms = sqrt( sumsq(X) ); 257 R = R ./ ( ones(rows(R),1) * reshape(rms.'*rms,1,[]) ); 258 elseif isempty(Y) 259 ## for autocorrelation, R(zero-lag) is the mean square. 260 R = R / R(maxlag+1); 261 else 262 ## for vectors X and Y 263 R = R / sqrt( sumsq(X)*sumsq(Y) ); 264 endif 265 elseif !strcmp(scale, 'none') 266 error("xcorr: scale must be 'biased', 'unbiased', 'coeff' or 'none'"); 267 endif 268 269 ## Pad result if necessary 270 ## (most likely is not required, use "if" to avoid unnecessary code) 271 ## At this point, lag varies with the first index in R; 272 ## so pad **before** the transpose. 273 if pad_result 274 R_pad = zeros(pad_result,columns(R)); 275 R = [R_pad; R; R_pad]; 276 endif 277 ## Correct the shape (transpose) so it is the same as the first input vector 278 if isvector(X) && P > 1 279 R = R.'; 280 endif 281 282 ## return the lag indices if desired 283 if nargout == 2 284 maxlag += pad_result; 285 lags = [-maxlag:maxlag]; 286 endif 287 288endfunction 289 290##------------ Use brute force to compute the correlation ------- 291##if !isvector(X) 292## ## For matrix X, compute cross-correlation of all columns 293## R = zeros(2*maxlag+1,P^2); 294## for i=1:P 295## for j=i:P 296## idx = (i-1)*P+j; 297## R(maxlag+1,idx) = X(:,i)' * X(:,j); 298## for k = 1:maxlag 299## R(maxlag+1-k,idx) = X(k+1:N,i)' * X(1:N-k,j); 300## R(maxlag+1+k,idx) = X(1:N-k,i)' * X(k+1:N,j); 301## endfor 302## if (i!=j), R(:,(j-1)*P+i) = conj(flipud(R(:,idx))); endif 303## endfor 304## endfor 305##elseif isempty(Y) 306## ## reshape X so that dot product comes out right 307## X = reshape(X, 1, N); 308## 309## ## compute autocorrelation for 0:maxlag 310## R = zeros (2*maxlag + 1, 1); 311## for k=0:maxlag 312## R(maxlag+1+k) = X(1:N-k) * X(k+1:N)'; 313## endfor 314## 315## ## use symmetry for -maxlag:-1 316## R(1:maxlag) = conj(R(2*maxlag+1:-1:maxlag+2)); 317##else 318## ## reshape and pad so X and Y are the same length 319## X = reshape(postpad(X,N), 1, N); 320## Y = reshape(postpad(Y,N), 1, N)'; 321## 322## ## compute cross-correlation 323## R = zeros (2*maxlag + 1, 1); 324## R(maxlag+1) = X*Y; 325## for k=1:maxlag 326## R(maxlag+1-k) = X(k+1:N) * Y(1:N-k); 327## R(maxlag+1+k) = X(1:N-k) * Y(k+1:N); 328## endfor 329##endif 330##-------------------------------------------------------------- 331