25

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.

NoDataDumpNoContribution
  • 10,591
  • 9
  • 64
  • 104
skj
  • 275
  • 1
  • 5
  • 7
  • 1
    Is there some reason why you can't just use MATLAB's existing autocorrelation function ? (Homework perhaps ?) – Paul R Oct 16 '10 at 16:41
  • @Paul R: `xcorr` is part of the signal-processing toolbox. – Oliver Charlesworth Oct 16 '10 at 19:33
  • 1
    @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 Oct 16 '10 at 20:28
  • 4
    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 Oct 17 '10 at 03:05

2 Answers2

36

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

John Clements
  • 16,895
  • 3
  • 37
  • 52
Amro
  • 123,847
  • 25
  • 243
  • 454
  • 2
    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. – endolith May 15 '12 at 18:16
  • 1
    why do you need to pad to `2^nextpow2(2*len-1)` and not `2^nextpow2(len)` ? – Marius Apr 11 '14 at 15:22
  • 1
    @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: http://dsp.stackexchange.com/q/741/679 – Amro Apr 11 '14 at 20:41
  • 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? – Marius Apr 14 '14 at 07:26
  • @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: http://blogs.uoregon.edu/seis/wiki/unpacking-the-matlab-fft/ – Amro Apr 14 '14 at 12:52
  • @Amro 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? – nevos Feb 11 '16 at 11:27
  • 2
    @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. " – Amro Feb 11 '16 at 22:07
  • 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. – rayryeng Nov 13 '17 at 01:43
32

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);
Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680