Calcola autocorrelazione utilizzando FFT in MATLAB

Question

Ho letto alcune spiegazioni su come l'autocorrelazione può essere calcolata in modo più efficiente utilizzando il segnale di un segnale, moltiplicando la parte reale per il complesso coniugato (dominio di Fourier), quindi utilizzando l'inverso di fft, ma Ho difficoltà a realizzarlo in MATLAB perché ad un livello dettagliato, non so davvero cosa sto facendo. : o) Qualsiasi tipo di anime là fuori si preoccupa di condividere qualche codice e saggezza?

Grazie!

Answer 1

Proprio come lei ha dichiarato, prendere la FFT e moltiplicare puntuale con il suo complesso coniugato, quindi utilizzare la FFT inversa (o, nel caso di cross-correlazione di 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, che è esattamente quello che sta facendo (solo deve fare i conti con tutti i casi di imbottitura, normalizzazione, vettore/ingresso a matrice, ecc ...)

Answer 2

A Wiener–Khinchin theorem, la densità di potenza spettrale (PSD) di una funzione è la trasformata di Fourier dell'autocorrelazione. Per i segnali deterministici, il PSD è semplicemente la magnitudo al quadrato della trasformata di Fourier. Vedi anche lo convolution theorem.

Quando si tratta di trasformate di Fourier discreta (cioè utilizzando FFT), si ottiene effettivamente l'autocorrelazione ciclica. Per ottenere una corretta autocorrelazione (lineare), prima di prendere la trasformata di Fourier è necessario azzerare i dati originali al doppio della loro lunghezza originale. 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);