function [c,lags] = xcor(x,varargin) %XCORR Cross-correlation function estimates. % C = XCORR(A,B), where A and B are length M vectors (M>1), returns % the length 2*M-1 cross-correlation sequence C. If A and B are of % different length, the shortest one is zero-padded. C will be a % row vector if A is a row vector, and a column vector if A is a % column vector. % % XCORR produces an estimate of the correlation between two random % (jointly stationary) sequences: % C(m) = E[A(n+m)*conj(B(n))] = E[A(n)*conj(B(n-m))] % It is also the deterministic correlation between two deterministic % signals. % % XCORR(A), when A is a vector, is the auto-correlation sequence. % XCORR(A), when A is an M-by-N matrix, is a large matrix with % 2*M-1 rows whose N^2 columns contain the cross-correlation % sequences for all combinations of the columns of A. % The zeroth lag of the output correlation is in the middle of the % sequence, at element or row M. % % XCORR(...,MAXLAG) computes the (auto/cross) correlation over the % range of lags: -MAXLAG to MAXLAG, i.e., 2*MAXLAG+1 lags. % If missing, default is MAXLAG = M-1. % % [C,LAGS] = XCORR(...) returns a vector of lag indices (LAGS). % % XCORR(...,SCALEOPT), normalizes the correlation according to SCALEOPT: % 'biased' - scales the raw cross-correlation by 1/M. % 'unbiased' - scales the raw correlation by 1/(M-abs(lags)). % 'coeff' - normalizes the sequence so that the auto-correlations % at zero lag are identically 1.0. % 'none' - no scaling (this is the default). % % See also XCOV, CORRCOEF, CONV, COV and XCORR2. % Author(s): R. Losada % Copyright 1988-2002 The MathWorks, Inc. % \$Revision: 1.16 \$ \$Date: 2002/03/28 17:31:18 \$ % % References: % S.J. Orfanidis, "Optimum Signal Processing. An Introduction" % 2nd Ed. Macmillan, 1988. error(nargchk(1,4,nargin)); [x,nshift] = shiftdim(x); [xIsMatrix,autoFlag,maxlag,scaleType,msg] = parseinput(x,varargin{:}); error(msg); if xIsMatrix, [c,M,N] = matrixCorr(x,maxlag); else [c,M,N] = vectorXcorr(x,autoFlag,maxlag,varargin{:}); end % Force correlation to be real when inputs are real c = forceRealCorr(c,x,autoFlag,varargin{:}); lags = [-maxlag:maxlag]; % Keep only the lags we want and move negative lags before positive lags if maxlag >= M, c = [zeros(maxlag-M+1,N^2);c(end-M+2:end,:);c(1:M,:);zeros(maxlag-M+1,N^2)]; else c = [c(end-maxlag+1:end,:);c(1:maxlag+1,:)]; end % Scale as specified [c,msg] = scaleXcorr(c,xIsMatrix,scaleType,autoFlag,M,maxlag,lags,x,varargin{:}); error(msg); % If first vector is a row, return a row c = shiftdim(c,-nshift); %---------------------------------------------------------------- function [c,M,N] = matrixCorr(x,maxlag) % Compute all possible auto- and cross-correlations for a matrix input [M,N] = size(x); X = fft(x,2^nextpow2(2*M-1)); Xc = conj(X); C = []; [MX,NX] = size(X); C = zeros(MX,NX*NX); for n =1:N, C(:,(((n-1)*N)+1):(n*N)) = repmat(X(:,n),1,N).*Xc; end c = ifft(C); %---------------------------------------------------------------- function [c,M,N] = vectorXcorr(x,autoFlag,maxlag,varargin) % Compute auto- or cross-correlation for vector inputs x = x(:); [M,N] = size(x); if autoFlag, % Autocorrelation % Compute correlation via FFT X = fft(x,2^nextpow2(2*M-1)); c = ifft(abs(X).^2); else, % xcorrelation y = varargin{1}; y = y(:); L = length(y); % Cache the length(x) Mcached = M; % Recompute length(x) in case length(y) > length(x) M = max(Mcached,L); if (L ~= Mcached) & any([L./Mcached, Mcached./L] > 10), % Vector sizes differ by a factor greater than 10, % fftfilt is faster neg_c = conj(fftfilt(conj(x),flipud(y))); % negative lags pos_c = flipud(fftfilt(conj(y),flipud(x))); % positive lags % Make them of almost equal length (remove zero-th lag from neg) lneg = length(neg_c); lpos = length(pos_c); neg_c = [zeros(lpos-lneg,1);neg_c(1:end-1)]; pos_c = [pos_c;zeros(lneg-lpos,1)]; c = [pos_c;neg_c]; else if L ~= Mcached, % Force equal lengths if L > Mcached x = [x;zeros(L-Mcached,1)]; else, y = [y;zeros(Mcached-L,1)]; end end % Transform both vectors X = fft(x,2^nextpow2(2*M-1)); Y = fft(y,2^nextpow2(2*M-1)); % Compute cross-correlation c = ifft(X.*conj(Y)); end end %---------------------------------------------------------------- function [c,msg] = scaleXcorr(c,xIsMatrix,scaleType,autoFlag,... M,maxlag,lags,x,varargin) % Scale correlation as specified msg = ''; switch scaleType, case 'none', return case 'biased', % Scales the raw cross-correlation by 1/M. c = c./M; case 'unbiased', % Scales the raw correlation by 1/(M-abs(lags)). scale = (M-abs(lags)).'; scale(scale<=0)=1; % avoid divide by zero, when correlation is zero if xIsMatrix, scale = repmat(scale,1,size(c,2)); end c = c./scale; case 'coeff', % Normalizes the sequence so that the auto-correlations % at zero lag are identically 1.0. if ~autoFlag, % xcorr(x,y) % Compute autocorrelations at zero lag cxx0 = sum(abs(x).^2); cyy0 = sum(abs(varargin{1}).^2); scale = sqrt(cxx0*cyy0); c = c./scale; else, if ~xIsMatrix, % Autocorrelation case, simply normalize by c[0] c = c./c(maxlag+1); else, % Compute the indices corresponding to the columns for which % we have autocorrelations (e.g. if c = n by 9, the autocorrelations % are at columns [1,5,9] the other columns are cross-correlations). [mc,nc] = size(c); jkl = reshape(1:nc,sqrt(nc),sqrt(nc))'; tmp = sqrt(c(maxlag+1,diag(jkl))); tmp = tmp(:)*tmp; cdiv = repmat(tmp(:).',mc,1); c = c ./ cdiv; % The autocorrelations at zero-lag are normalized to % one end end end %---------------------------------------------------------------- function [xIsMatrix,autoFlag,maxlag,scaleType,msg] = parseinput(x,varargin) % Parse the input and determine optional parameters: % % Outputs: % xIsMatrix - flag indicating if x is a matrix % AUTOFLAG - 1 if autocorrelation, 0 if xcorrelation % maxlag - Number or lags to compute % scaleType - String with the type of scaling wanted % msg - possible error message % Set some defaults msg = ''; scaleType = ''; autoFlag = 1; % Assume autocorrelation until proven otherwise maxlag = []; errMsg = 'Input argument is not recognized.'; switch nargin, case 2, % Can be (x,y), (x,maxlag), or (x,scaleType) if ischar(varargin{1}), % Second arg is scaleType scaleType = varargin{1}; elseif isnumeric(varargin{1}), % Can be y or maxlag if length(varargin{1}) == 1, maxlag = varargin{1}; else autoFlag = 0; y = varargin{1}; end else % Not recognized msg = errMsg; return end case 3, % Can be (x,y,maxlag), (x,maxlag,scaleType) or (x,y,scaleType) maxlagflag = 0; % By default, assume 3rd arg is not maxlag if ischar(varargin{2}), % Must be scaletype scaleType = varargin{2}; elseif isnumeric(varargin{2}), % Must be maxlag maxlagflag = 1; maxlag = varargin{2}; else % Not recognized msg = errMsg; return end if isnumeric(varargin{1}), if maxlagflag, autoFlag = 0; y = varargin{1}; else % Can be y or maxlag if length(varargin{1}) == 1, maxlag = varargin{1}; else autoFlag = 0; y = varargin{1}; end end else % Not recognized msg = errMsg; return end case 4, % Must be (x,y,maxlag,scaleType) autoFlag = 0; y = varargin{1}; maxlag = varargin{2}; scaleType = varargin{3}; end % Determine if x is a matrix or a vector [xIsMatrix,m] = parse_x(x); if ~autoFlag, % Test y for correctness [maxlag,msg] = parse_y(y,m,xIsMatrix,maxlag); if ~isempty(msg), return end end [maxlag,msg] = parse_maxlag(maxlag,m); if ~isempty(msg), return end % Test the scaleType validity [scaleType,msg] = parse_scaleType(scaleType,errMsg,autoFlag,m,varargin{:}); if ~isempty(msg), return end %------------------------------------------------------------------- function [xIsMatrix,m] = parse_x(x) xIsMatrix = (size(x,2) > 1); m = size(x,1); %------------------------------------------------------------------- function [maxlag,msg] = parse_y(y,m,xIsMatrix,maxlag) msg = ''; [my,ny] = size(y); if ~any([my,ny] == 1), % Second arg is a matrix, error msg = 'B must be a vector (min(size(B))==1).'; return end if xIsMatrix, % Can't do xcorr(matrix,vector) msg = 'When B is a vector, A must be a vector.'; return end if (length(y) > m) & isempty(maxlag), % Compute the default maxlag based on the length of y maxlag = length(y) - 1; end %------------------------------------------------------------------- function [maxlag,msg] = parse_maxlag(maxlag,m) msg = ''; if isempty(maxlag), % Defaul hasn't been assigned yet, do so maxlag = m-1; else % Test maxlag for correctness if length(maxlag)>1 msg = 'Maximum lag must be a scalar.'; return end if maxlag < 0, maxlag = abs(maxlag); end if maxlag ~= round(maxlag), msg = 'Maximum lag must be an integer.'; end end %-------------------------------------------------------------------- function c = forceRealCorr(c,x,autoFlag,varargin) % Force correlation to be real when inputs are real forceReal = 0; % Flag to determine whether we should force real corr if (isreal(x) & autoFlag) | (isreal(x) & isreal(varargin{1})), forceReal = 1; end if forceReal, c = real(c); end %-------------------------------------------------------------------- function [scaleType,msg] = parse_scaleType(scaleType,errMsg,autoFlag,m,varargin) msg = ''; if isempty(scaleType), scaleType = 'none'; else scaleOpts = {'biased','unbiased','coeff','none'}; indx = strmatch(lower(scaleType),scaleOpts); if isempty(indx), msg = errMsg; return else scaleType = scaleOpts{indx}; end if ~autoFlag & ~strcmpi(scaleType,'none') &(m ~= length(varargin{1})), msg = 'Scale option must be ''none'' for different length vectors A and B.'; return end end % EOF