Calculate autocorrelation using FFT in Matlab

23
votes

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.

2
matlabsignal-processingfftcorrelation
Is there some reason why you can't just use MATLAB's existing autocorrelation function ? (Homework perhaps ?) – Paul R
@Paul R: xcorr is part of the signal-processing toolbox. – Oliver Charlesworth
@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. – Paul R
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) – skj

2 Answers

32
votes

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...)

31
votes

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);