Intuitive explanation of cross-correlation in frequency domain

Question

According to the cross-correlation theorem : the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal. After doing this, when we take the ifft of the product signal, we get a peak which indicates the shift between two signals.

I am not able to understand how this works? Why would i get a peak which indicates the shift between two signals. I got the math from : http://mathworld.wolfram.com/Cross-CorrelationTheorem.html but i am not able to understand as to what this means intuitatively. Can somebody please provide some explanation or point me to the right documents?

Thanks!

Answer 1

The concept is based on the convolution theorem, which states that for two signals $x(t)$ and $y(t)$, the product of their Fourier transforms $X(f)$ and $Y(f)$ is equal to the Fourier transform of the convolution of the two signals. That is:

$$ \mathcal{F}\{x(t) * y(t)\} = \mathcal{F}\{x(t)\}\mathcal{F}\{y(t)\} $$

You can read more on the derivation of this theorem at the above Wikipedia link. Now, convolution is a very important operation for linear systems in itself, so the theory on its properties is well-developed.

However, what you're looking for is the cross-correlation between $x(t)$ and $y(t)$. Here's the key: the cross-correlation integral is equivalent to the convolution integral if one of the input signals is conjugated and time-reversed. This allows you to utilize theory developed for evaluating convolutions (like frequency-domain techniques for calculating them quickly) and apply them to correlations.

In your example, you're calculating the following:

$$ \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^* $$

Recall that in the Fourier domain, complex conjugation is equivalent to time reversal in the time domain (this follows directly from the definition of the Fourier transform). Therefore, using the first equation given above, we can state that:

$$ \mathcal{F}\{x(t) * y^*(-t)\} = \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^* $$

If you then take the inverse Fourier transform of this equation, the signal you're left with is the cross-correlation between $x(t)$ and $y(t)$.

If you are working with real signals then we drop the complex conjugate in $y(t)$.

$$ \mathcal{F}\{x(t) * y(-t)\} = \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^* $$

And it is very easy to see that for real signals, cross correlation and convolution are equivalent if we flip one of the signals in time. In this case the convolution operation flip in time domain is compensated with another flip in $y(t)$ to yield the cross correlation on the left hand side of the last equation.

Answer 2

This can be confusing. I like Jason R's explanation, but found some code useful.

Given two random complex column vectors x and y of length M and N then in MATLAB you can verify that conv(x,flip(conj(y)) gives the same output as xcorr(x,y) except that the output of xcorr is controlled by the maxlag parameter, the default being max(M,N)-1 whereas the output of conv is length M+N-1. So, with the exception of some indexing differences, the two commands are the same.

Now, you can also verify that conv(x,y) = ifft(fft(x,L).*fft(y,L) where L = M+N-1.

So it follows that corr(x,y) = ifft(fft(x,L).*fft(flip(conj(y)),L)

(Again, taking care with the indexing of the lags.)

If L is big, you can also use a number bigger than L for the fft's (something that makes for an efficient FFT length) and trim the result to length L.

Answer 3

% Matlab function for frequency domain cross correlation
function [Lag,C]=xcorrf(X,Y,L)
% X, Y ---> Input vectors 
% L --->  maximum lag (must be less than minimum of (length of X, Y)
% C ---> correlation vector
% Lag ---> lag times  
X=X(:);
Y=Y(:);
s1=size(X);
s2=size(Y);
D=min(s1(1,1),s2(1,1));
for i=1:L
    X1=ifft(fft(X(1:D-i,:)).*conj(fft(Y(i+1:D,1))));
    C(i,1)=X1(1,1);
end

C=flipud(C);
X1=ifft(fft(X(1:D,:)).*conj(fft(Y(1:D,1))));
C(L+1,1)=X1(1,1);
for i=1:L
    X1=ifft(fft(Y(1:D-i,:)).*conj(fft(X(i+1:D,1))));
    C(i+L+1,1)=X1(1,1);
end
Lag=-L:1:L;
end