Tính tự tương quan bằng cách sử dụng FFT trong MATLAB

Question

Tôi đã đọc một số giải thích về cách tự tương quan có thể được tính toán hiệu quả hơn bằng cách sử dụng fft của tín hiệu, nhân phần thực bởi liên hợp phức tạp (miền nhiễu), sau đó sử dụng nghịch đảo fft, Tôi gặp khó khăn khi nhận ra điều này trong MATLAB bởi vì ở mức độ chi tiết, tôi thực sự không biết mình đang làm gì. : o) Có linh hồn nào ngoài kia quan tâm để chia sẻ một số mã và sự khôn ngoan?

Cảm ơn!

Answer 1

Cũng giống như bạn đã nêu, đi FFT và nhân pointwise bởi liên hợp phức tạp của nó, sau đó sử dụng FFT ngược (hoặc trong trường hợp tương quan chéo của hai tín hiệu: 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) 

trong thực tế, nếu bạn nhìn vào mã của xcorr.m, đó là chính xác những gì nó làm (chỉ nó có để đối phó với tất cả các trường hợp đệm, bình thường hóa, vector/ma trận đầu vào, vv ...)

Answer 2

Bằng số Wiener–Khinchin theorem, mật độ phổ công suất (PSD) của hàm là biến đổi Fourier của tự tương quan. Đối với các tín hiệu xác định, PSD đơn giản là độ lớn bình phương của biến đổi Fourier. Xem thêm convolution theorem.

Khi nói đến các biến đổi Fourier rời rạc (tức là sử dụng FFT), bạn thực sự nhận được sự tương quan tự tuần hoàn. Để có được sự tương quan đúng đắn (tuyến tính), bạn phải tạo 0-pad dữ liệu ban đầu gấp hai lần chiều dài ban đầu của nó trước khi chuyển đổi Fourier. Vì vậy, một cái gì đó như:

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