X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fxcorr.m;fp=octave_packages%2Fsignal-1.1.3%2Fxcorr.m;h=15207d02f160f365d89a7c0373c63e5e46a786fc;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/signal-1.1.3/xcorr.m b/octave_packages/signal-1.1.3/xcorr.m new file mode 100644 index 0000000..15207d0 --- /dev/null +++ b/octave_packages/signal-1.1.3/xcorr.m @@ -0,0 +1,283 @@ +## Copyright (C) 1999-2001 Paul Kienzle +## Copyright (C) 2004 +## Copyright (C) 2008,2010 Peter Lanspeary +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + +## usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale]) +## +## Estimate the cross correlation R_xy(k) of vector arguments X and Y +## or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X, +## for a range of lags k specified by argument "maxlag". If X is a +## matrix, each column of X is correlated with itself and every other +## column. +## +## The cross-correlation estimate between vectors "x" and "y" (of +## length N) for lag "k" is given by +## R_xy(k) = sum_{i=1}^{N}{x_{i+k} conj(y_i), +## where data not provided (for example x(-1), y(N+1)) is zero. +## +## ARGUMENTS +## X [non-empty; real or complex; vector or matrix] data +## +## Y [real or complex vector] data +## If X is a matrix (not a vector), Y must be omitted. +## Y may be omitted if X is a vector; in this case xcorr +## estimates the autocorrelation of X. +## +## maxlag [integer scalar] maximum correlation lag +## If omitted, the default value is N-1, where N is the +## greater of the lengths of X and Y or, if X is a matrix, +## the number of rows in X. +## +## scale [character string] specifies the type of scaling applied +## to the correlation vector (or matrix). is one of: +## 'none' return the unscaled correlation, R, +## 'biased' return the biased average, R/N, +## 'unbiased' return the unbiassed average, R(k)/(N-|k|), +## 'coeff' return the correlation coefficient, R/(rms(x).rms(y)), +## where "k" is the lag, and "N" is the length of X. +## If omitted, the default value is "none". +## If Y is supplied but does not have the ame length as X, +## scale must be "none". +## +## RETURNED VARIABLES +## R array of correlation estimates +## lag row vector of correlation lags [-maxlag:maxlag] +## +## The array of correlation estimates has one of the following forms. +## (1) Cross-correlation estimate if X and Y are vectors. +## (2) Autocorrelation estimate if is a vector and Y is omitted, +## (3) If X is a matrix, R is an matrix containing the cross- +## correlation estimate of each column with every other column. +## Lag varies with the first index so that R has 2*maxlag+1 +## rows and P^2 columns where P is the number of columns in X. +## If Rij(k) is the correlation between columns i and j of X +## R(k+maxlag+1,P*(i-1)+j) == Rij(k) +## for lag k in [-maxlag:maxlag], or +## R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)). +## "reshape(R(k,:),P,P)" is the cross-correlation matrix for X(k,:). +## + +## The cross-correlation estimate is calculated by a "spectral" method +## in which the FFT of the first vector is multiplied element-by-element +## with the FFT of second vector. The computational effort depends on +## the length N of the vectors and is independent of the number of lags +## requested. If you only need a few lags, the "direct sum" method may +## be faster: +## +## Ref: Stearns, SD and David, RA (1988). Signal Processing Algorithms. +## New Jersey: Prentice-Hall. + +## unbiased: +## ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ [N:-1:N-k+1](:) +## biased: +## ( hankel(x(1:k),[x(k:N); zeros(k-1,1)]) * y ) ./ N +## +## If length(x) == length(y) + k, then you can use the simpler +## ( hankel(x(1:k),x(k:N-k)) * y ) ./ N + +function [R, lags] = xcorr (X, Y, maxlag, scale) + + if (nargin < 1 || nargin > 4) + print_usage; + endif + + ## assign arguments that are missing from the list + ## or reassign (right shift) them according to data type + if nargin==1 + Y=[]; maxlag=[]; scale=[]; + elseif nargin==2 + maxlag=[]; scale=[]; + if ischar(Y), scale=Y; Y=[]; + elseif isscalar(Y), maxlag=Y; Y=[]; + endif + elseif nargin==3 + scale=[]; + if ischar(maxlag), scale=maxlag; maxlag=[]; endif + if isscalar(Y), maxlag=Y; Y=[]; endif + endif + + ## assign defaults to missing arguments + if isvector(X) + ## if isempty(Y), Y=X; endif ## this line disables code for autocorr'n + N = max(length(X),length(Y)); + else + N = rows(X); + endif + if isempty(maxlag), maxlag=N-1; endif + if isempty(scale), scale='none'; endif + + ## check argument values + if isempty(X) || isscalar(X) || ischar(Y) || ! ismatrix(X) + error("xcorr: X must be a vector or matrix"); + endif + if isscalar(Y) || ischar(Y) || (!isempty(Y) && !isvector(Y)) + error("xcorr: Y must be a vector"); + endif + if !isempty(Y) && !isvector(X) + error("xcorr: X must be a vector if Y is specified"); + endif + if !isscalar(maxlag) || !isreal(maxlag) || maxlag<0 || fix(maxlag)!=maxlag + error("xcorr: maxlag must be a single non-negative integer"); + endif + ## + ## sanity check on number of requested lags + ## Correlations for lags in excess of +/-(N-1) + ## (a) are not calculated by the FFT algorithm, + ## (b) are all zero; so provide them by padding + ## the results (with zeros) before returning. + if (maxlag > N-1) + pad_result = maxlag - (N - 1); + maxlag = N - 1; + %error("xcorr: maxlag must be less than length(X)"); + else + pad_result = 0; + endif + if isvector(X) && isvector(Y) && length(X) != length(Y) && \ + !strcmp(scale,'none') + error("xcorr: scale must be 'none' if length(X) != length(Y)") + endif + + P = columns(X); + M = 2^nextpow2(N + maxlag); + if !isvector(X) + ## For matrix X, correlate each column "i" with all other "j" columns + R = zeros(2*maxlag+1,P^2); + + ## do FFTs of padded column vectors + pre = fft (postpad (prepad (X, N+maxlag), M) ); + post = conj (fft (postpad (X, M))); + + ## do autocorrelations (each column with itself) + ## -- if result R is reshaped to 3D matrix (i.e. R=reshape(R,M,P,P)) + ## the autocorrelations are on leading diagonal columns of R, + ## where i==j in R(:,i,j) + cor = ifft (post .* pre); + R(:, 1:P+1:P^2) = cor (1:2*maxlag+1,:); + + ## do the cross correlations + ## -- these are the off-diagonal colummn of the reshaped 3D result + ## matrix -- i!=j in R(:,i,j) + for i=1:P-1 + j = i+1:P; + cor = ifft( pre(:,i*ones(length(j),1)) .* post(:,j) ); + R(:,(i-1)*P+j) = cor(1:2*maxlag+1,:); + R(:,(j-1)*P+i) = conj( flipud( cor(1:2*maxlag+1,:) ) ); + endfor + elseif isempty(Y) + ## compute autocorrelation of a single vector + post = fft( postpad(X(:),M) ); + cor = ifft( post .* conj(post) ); + R = [ conj(cor(maxlag+1:-1:2)) ; cor(1:maxlag+1) ]; + else + ## compute cross-correlation of X and Y + ## If one of X & Y is a row vector, the other can be a column vector. + pre = fft( postpad( prepad( X(:), length(X)+maxlag ), M) ); + post = fft( postpad( Y(:), M ) ); + cor = ifft( pre .* conj(post) ); + R = cor(1:2*maxlag+1); + endif + + ## if inputs are real, outputs should be real, so ignore the + ## insignificant complex portion left over from the FFT + if isreal(X) && (isempty(Y) || isreal(Y)) + R=real(R); + endif + + ## correct for bias + if strcmp(scale, 'biased') + R = R ./ N; + elseif strcmp(scale, 'unbiased') + R = R ./ ( [ N-maxlag:N-1, N, N-1:-1:N-maxlag ]' * ones(1,columns(R)) ); + elseif strcmp(scale, 'coeff') + ## R = R ./ R(maxlag+1) works only for autocorrelation + ## For cross correlation coeff, divide by rms(X)*rms(Y). + if !isvector(X) + ## for matrix (more than 1 column) X + rms = sqrt( sumsq(X) ); + R = R ./ ( ones(rows(R),1) * reshape(rms.'*rms,1,[]) ); + elseif isempty(Y) + ## for autocorrelation, R(zero-lag) is the mean square. + R = R / R(maxlag+1); + else + ## for vectors X and Y + R = R / sqrt( sumsq(X)*sumsq(Y) ); + endif + elseif !strcmp(scale, 'none') + error("xcorr: scale must be 'biased', 'unbiased', 'coeff' or 'none'"); + endif + + ## Pad result if necessary + ## (most likely is not required, use "if" to avoid uncessary code) + ## At this point, lag varies with the first index in R; + ## so pad **before** the transpose. + if pad_result + R_pad = zeros(pad_result,columns(R)); + R = [R_pad; R; R_pad]; + endif + ## Correct the shape (transpose) so it is the same as the first input vector + if isvector(X) && P > 1 + R = R.'; + endif + + ## return the lag indices if desired + if nargout == 2 + maxlag += pad_result; + lags = [-maxlag:maxlag]; + endif + +endfunction + +##------------ Use brute force to compute the correlation ------- +##if !isvector(X) +## ## For matrix X, compute cross-correlation of all columns +## R = zeros(2*maxlag+1,P^2); +## for i=1:P +## for j=i:P +## idx = (i-1)*P+j; +## R(maxlag+1,idx) = X(i)*X(j)'; +## for k = 1:maxlag +## R(maxlag+1-k,idx) = X(k+1:N,i) * X(1:N-k,j)'; +## R(maxlag+1+k,idx) = X(k:N-k,i) * X(k+1:N,j)'; +## endfor +## if (i!=j), R(:,(j-1)*P+i) = conj(flipud(R(:,idx))); endif +## endfor +## endfor +##elseif isempty(Y) +## ## reshape X so that dot product comes out right +## X = reshape(X, 1, N); +## +## ## compute autocorrelation for 0:maxlag +## R = zeros (2*maxlag + 1, 1); +## for k=0:maxlag +## R(maxlag+1+k) = X(1:N-k) * X(k+1:N)'; +## endfor +## +## ## use symmetry for -maxlag:-1 +## R(1:maxlag) = conj(R(2*maxlag+1:-1:maxlag+2)); +##else +## ## reshape and pad so X and Y are the same length +## X = reshape(postpad(X,N), 1, N); +## Y = reshape(postpad(Y,N), 1, N)'; +## +## ## compute cross-correlation +## R = zeros (2*maxlag + 1, 1); +## R(maxlag+1) = X*Y; +## for k=1:maxlag +## R(maxlag+1-k) = X(k+1:N) * Y(1:N-k); +## R(maxlag+1+k) = X(k:N-k) * Y(k+1:N); +## endfor +##endif +##--------------------------------------------------------------