Efficiently calculating autocorrelation using FFTs

Question

I'm trying to calculate an autocorrelation on a platform where the only accelerated primitive I have available is the (I)FFT. I'm having a problem though.

I prototyped it in MATLAB. I am, however, slightly confused. I assumed that it works simply as follows (this is from memory so apologies if I've got it slightly wrong).

 autocorr = ifft( complex( abs( fft( inputData ) ), 0 ) )

However I get a different result than I get from using the xcorr function. Now I'm fully expecting not to get the left hand side of the auto correlation (as it's a reflection of the right hand side and thus not needed anyway). However, the problem is my right hand side appears to be, itself, reflected around the halfway point. Which effectively means that I get about half the amount of data that I'm expecting.

So I'm sure I must be doing something very simple wrong, but I just can't figure out what.

Answer 1

pichenettes is right, of course. The FFT implements a circular convolution while the xcorr() is based on a linear convolution. In addition you need to square the absolute value in the frequency domain as well. Here is a code snippet that handles all the zero padding, shifting & truncating.

%% Cross correlation through a FFT
n = 1024;
x = randn(n,1);
% cross correlation reference
xref = xcorr(x,x);

%FFT method based on zero padding
fx = fft([x; zeros(n,1)]); % zero pad and FFT
x2 = ifft(fx.*conj(fx)); % abs()^2 and IFFT
% circulate to get the peak in the middle and drop one
% excess zero to get to 2*n-1 samples
x2 = [x2(n+2:end); x2(1:n)];
% calculate the error
d = x2-xref; % difference, this is actually zero
fprintf('Max error = %6.2f\n',max(abs(d)));

Answer 2

Matlab's xcorr computes a $2N - 1$ FFT, where $N$ is the length of the input data (ie, the input is padded with $N - 1$ zeros). This avoids the circularity problem.

Answer 3

To elaborate a bit more on the previous answers, computing the auto-correlation of a signal of length $N$ results in a (sampled) auto-correlation of size $2N-1$. Well, actually, it should be infinite, but the auto-correlation outside $[-(N-1), N-1]$ equals $0$ anyway.

Now, you desire to use the discrete Fourier transform (DFT) to compute it, and the formula is indeed the inverse DFT of the squared magnitude of the DFT of your signal. But think about it: if we take it the other way around and compute the DFT of the auto-correlation, you end up with a spectrum of size $2N-1$, if you don't want to lose samples in the way! That spectrum therefore has to be of size $2N-1$, and that's the reason why you need to zero-pad your time-domain signal up to $2N-1$, compute the DFT (on $2N-1$ points), and proceed with it.

Another way of seeing this is to analyze what happens if you compute the DFT on $N$ points: this is equivalent to downsampling your (continuous frequency) discrete time Fourier transform (DTFT). Retrieving the auto-correlation, which should be of size $2N-1$, with an under-sampled spectrum of size $N$ therefore leads to time aliasing (the circularity pichenettes was talking about), which explains why you have this symmetrical pattern in the "right hand side" if your output.

Actually, the code provided by Hilmar also works, because as long as you zero-pad up to more than a size of $N-1$ (in his case, he computes an FT of size $N$), you "over-sample" your FT, and you still get your $2N-1$ "useful" samples (the other ones should be $0$s). So, for efficiency, just zero-pad to $2N-1$, that's all what you need (well, perhaps you'd better zero-pad up to the next power of 2 of $2N-1$, if you use FFTs).

In short: you should have done this (to be adapted to your programming language):

autocorr = ifft( complex( abs(fft(inputData, n=2*N-1))**2, 0 ) )

Or in MATLAB:

autocorr = ifft(abs(fft(inputData, 2*N-1)).^2)

Answer 4

The main reason for the desired output of xcorr function to be not similar to that of application of FFT and IFFT function is because while applying these function to signals the final result is circularly convoluted.

The main difference between Linear Convolution and Circular Convolution can be found in Linear and Circular Convolution.

The problem can be solved by Initially zero-padding the signal and truncating final output of IFFT.