Calculate autocorrelation using FFT in Matlab

Question

I've read some explanations of how autocorrelation can be more efficiently calculated using the fft of a signal, multiplying the real part by the complex conjugate (Fourier domain), then using the inverse fft, but I'm having trouble realizing this in Matlab because at a detailed level.

Padlock

Is there some reason why you can't just use MATLAB's existing autocorrelation function ? (Homework perhaps ?) – Donofrio

@Paul R: xcorr is part of the signal-processing toolbox. – Tommyetommyrot

@Oli: OK - I guess the OP doesn't have the Signal-Processing Toolbox ? I use Octave rather than MATLAB and it seems to have xcorr. – Donofrio

I have the signal processing toolbox, but I'm just trying to understand the ACF better, particularly wrt any optimizations b/c I'll eventually port the algorithm I'm working on to C# (eek!). It's not HW, btw. :o) – Padlock

Answer 1

Just like you stated, take the fft and multiply pointwise by its complex conjugate, then use the inverse fft (or in the case of cross-correlation of two signals: Corr(x,y) <=> FFT(x)FFT(y)*)

x = rand(100,1);
len = length(x);

%# autocorrelation
nfft = 2^nextpow2(2*len-1);
r = ifft( fft(x,nfft) .* conj(fft(x,nfft)) );

%# rearrange and keep values corresponding to lags: -(len-1):+(len-1)
r = [r(end-len+2:end) ; r(1:len)];

%# compare with MATLAB's XCORR output
all( (xcorr(x)-r) < 1e-10 )

In fact, if you look at the code of xcorr.m, that's exactly what it's doing (only it has to deal with all the cases of padding, normalizing, vector/matrix input, etc...)

Brumley

Instead of taking the complex conjugate, you can also reverse one signal and then do the FFT of that. One might be easier than the other, depending on your program. – Wordy

why do you need to pad to 2^nextpow2(2*len-1) and not 2^nextpow2(len) ? – Ation

@Marius: it is to perform the FFT while padding the signal with zeros to double its size (preferably rounded up to next power of two as the implementation is faster that way). See this question for more explanations: dsp.stackexchange.com/q/741/679 – Brumley

If I compute the autocorrelation in this way the result will be symmetric with respect to the timelag and I don't need points 'r(end-len+2:end)'. Am I getting this right? – Ation

@Marius: yes the DFT of a real signal is conjugate-symmetric, so you can drop one-half if you want. The above code rearranges the correlation values to be symmetric with respect to the time lags (going from negative lags, to zero, to positive lags). Here is a nice explanation: blogs.uoregon.edu/seis/wiki/unpacking-the-matlab-fft – Brumley

@Brumley I don't understand why the result of the ifft is a real vector? I thought it should return a complex vector. When I try to pop in a random real vector in the ifft it returns a complex vector, so how come r is real? – Perdue

@nevos: the keyword is Hermitian symmetry. To quote help ifft: "ifft tests X to see whether vectors in X along the active dimension are conjugate symmetric. If so, the computation is faster and the output is real. An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x. " – Brumley

Thank you @Amro. Great answer - exactly what I needed for my application. BTW, I hope all is well. We haven't spoken in a while. – Contrary

Answer 2

By the Wiener–Khinchin theorem, the power-spectral density (PSD) of a function is the Fourier transform of the autocorrelation. For deterministic signals, the PSD is simply the magnitude-squared of the Fourier transform. See also the convolution theorem.

When it comes to discrete Fourier transforms (i.e. using FFTs), you actually get the cyclic autocorrelation. In order to get proper (linear) autocorrelation, you must zero-pad the original data to twice its original length before taking the Fourier transform. So something like:

x = [ ... ];
x_pad = [x zeros(size(x))];
X     = fft(x_pad);
X_psd = abs(X).^2;
r_xx = ifft(X_psd);

Tommyetommyrot