Calcolare autocorrelazione tramite FFT in Matlab

Question

Ho letto alcune spiegazioni di come autocorrelazione può essere calcolato in modo più efficiente utilizzando la FFT di un segnale, moltiplicando la parte reale da parte del complesso coniugato (dominio di Fourier), quindi utilizzando la FFT inversa, ma sto avendo difficoltà realizzando questo in Matlab, perché a livello dettagliato.

Answer 1

Come te affermato, prendere la FFT e puntuale moltiplicare per il suo complesso coniugato, quindi utilizzare la FFT inversa (o in caso di cross-correlazione dei due segnali: 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 realtà, se si guarda il codice di xcorr.m, questo è esattamente quello che sta facendo (solo deve fare i conti con tutti i casi di imbottitura, normalizzazione, vettore / ingresso a matrice, ecc ...)

Answer 2

Teorema di Wiener-Khinchin , la densità di potenza spettrale (PSD) di un funzione è la trasformata di Fourier dell'autocorrelazione. Per segnali deterministici, il PSD è semplicemente l'ampiezza-quadrato della trasformata di Fourier. Si veda anche la circonvoluzione teorema .

Quando si tratta di Trasformata discreta di Fourier (cioè utilizzando FFT), è effettivamente ottenere l'autocorrelazione ciclica. Al fine di ottenere una corretta autocorrelazione (lineare), è necessario azzerare-pad i dati originali a due volte la sua lunghezza originale prima di prendere la trasformata di Fourier. Quindi, qualcosa di simile a:

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