octave_packages/signal-1.1.3/xcorr.m

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