Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Calculate autocorrelation using FFT in Matlab

Tags:

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.

like image 503
skj Avatar asked Oct 16 '10 14:10

skj


People also ask

How does MATLAB calculate autocorrelation?

The autocorrelation function measures the correlation between the univariate time series yt and yt+k, where k = 0,...,K and yt is a stochastic process. c k = 1 T ∑ t = 1 T − k ( y t − y ¯ ) ( y t + k − y ¯ ) . c0 is the sample variance of the time series.

How does fft in MATLAB work?

Y = fft( X ) computes the discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm. If X is a vector, then fft(X) returns the Fourier transform of the vector. If X is a matrix, then fft(X) treats the columns of X as vectors and returns the Fourier transform of each column.

What does autocorrelation mean in MATLAB?

Autocorrelation is the linear dependence of a variable with itself at two points in time. For stationary processes, autocorrelation between any two observations depends only on the time lag h between them. Define Cov(yt, yt–h) = γh. Lag-h autocorrelation is given by. ρ h = C o r r ( y t , y t − h ) = γ h γ 0 .


1 Answers

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

like image 156
Amro Avatar answered Sep 20 '22 18:09

Amro