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.
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.
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.
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 .
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...)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With