Calcular autocorrelación usando FFT en matlab

Question

He leído algunas explicaciones de cómo la autocorrelación puede calcularse más eficientemente utilizando la fft de una señal, multiplicando la parte real por el conjugado complejo (dominio de Fourier), luego usando el fft inverso, pero Tengo problemas para darme cuenta de esto en Matlab porque, en un nivel detallado, realmente no sé lo que estoy haciendo. : o) ¿Alguna alma amable por ahí se preocupa por compartir algún código y sabiduría?

Gracias!

Answer 1

Al igual que usted dijo, tomar la FFT y se multiplica por puntos por su complejo conjugado, a continuación, utilizar la FFT inversa (o en el caso de correlación cruzada de dos señales: 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) 

de hecho, si nos fijamos en el código de xcorr.m, eso es exactamente lo que está haciendo (sólo tiene que ocuparse de todos los casos de relleno, la normalización, el vector/matriz de entrada, etc ...)

Answer 2

Por el Wiener–Khinchin theorem, la densidad espectral de potencia (PSD) de una función es la transformada de Fourier de la autocorrelación. Para las señales deterministas, el PSD es simplemente la magnitud cuadrada de la transformada de Fourier. Vea también el convolution theorem.

Cuando se trata de transformaciones discretas de Fourier (es decir, mediante el uso de FFT), se obtiene la autocorrelación cíclica. Para obtener una autocorrelación (lineal) adecuada, debe rellenar con cero los datos originales al doble de su longitud original antes de tomar la transformada de Fourier. Así que algo como:

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