Skutecznie oblicza autokorelację za pomocą FFT

Próbuję obliczyć autokorelację na platformie, na której jedyną dostępną przyspieszoną operacją podstawową jest (I) FFT. Mam jednak problem.

Prototypowałem go w MATLAB . Jestem jednak nieco zdezorientowany. Założyłem, że działa po prostu w następujący sposób (pochodzi z pamięci, więc przepraszam, jeśli się trochę mylę).

 autocorr = ifft( complex( abs( fft( inputData ) ), 0 ) )

Jednak uzyskuję inny wynik niż przy użyciu tej xcorrfunkcji. Teraz w pełni oczekuję, że nie otrzymam lewej strony automatycznej korelacji (ponieważ jest to odbicie prawej strony, a zatem i tak nie jest potrzebne). Problemem jest jednak to, że moja prawa ręka wydaje się być odbita w połowie odległości. Co faktycznie oznacza, że ​​otrzymuję około połowy danych, których oczekuję.

Jestem więc pewien, że robię coś bardzo prostego, ale nie wiem, co.

oczywiście feniksy mają rację. FFT implementuje splot kołowy, podczas gdy xcorr () opiera się na splotie liniowym. Ponadto należy również obliczyć wartość bezwzględną w dziedzinie częstotliwości. Oto fragment kodu, który obsługuje wypełnianie zerami, przesuwanie i obcinanie.

%% Cross correlation through a FFT
n = 1024;
x = randn(n,1);
% cross correlation reference
xref = xcorr(x,x);

%FFT method based on zero padding
fx = fft([x; zeros(n,1)]); % zero pad and FFT
x2 = ifft(fx.*conj(fx)); % abs()^2 and IFFT
% circulate to get the peak in the middle and drop one
% excess zero to get to 2*n-1 samples
x2 = [x2(n+2:end); x2(1:n)];
% calculate the error
d = x2-xref; % difference, this is actually zero
fprintf('Max error = %6.2f\n',max(abs(d)));

Xcorr Matlaba oblicza FFT , gdzie jest długością danych wejściowych (tzn. Dane wejściowe są wypełnione zerami ). Pozwala to uniknąć problemu okrągłości.$2N - 1$$N$$N - 1$

Aby rozwinąć nieco więcej na temat poprzednich odpowiedzi, obliczenie autokorelacji sygnału o długości skutkuje (próbkowaną) autokorelacją wielkości . Właściwie powinno być nieskończone, ale autokorelacja poza tak wynosi .$N$$2N-1$$[-(N-1), N-1]$$0$

Teraz chcesz użyć dyskretnej transformaty Fouriera (DFT) do jej obliczenia, a formuła jest w rzeczywistości odwrotnym DFT kwadratowej wielkości DFT twojego sygnału. Ale pomyśl o tym: jeśli weźmiemy to na odwrót i obliczymy DFT autokorelacji, uzyskasz spektrum wielkości , jeśli nie chcesz stracić próbek! Widmo to musi więc mieć rozmiar i dlatego musisz zerować swój sygnał w dziedzinie czasu do , obliczyć DFT (w punktach ) i kontynuować.2 N - 1 2 N - 1 2 N - $2N-1$$2N-1$$2N-1$$2N-1$

Innym sposobem na to jest analiza tego, co się stanie, jeśli obliczymy DFT w punktach : jest to równoważne z próbkowaniem w dół twojego dyskretnego czasu (ciągłej częstotliwości) transformaty Fouriera (DTFT). Odzyskanie autokorelacji, która powinna być wielkości , z niedopróbkowanym spektrum wielkości prowadzi zatem do aliasingu w czasie (o których mówiono o fenenetach o okrągłości), co wyjaśnia, dlaczego masz ten symetryczny wzór po prawej stronie hand side ”, jeśli Twój wynik.2 N - 1 $N$$2N-1$$N$

W rzeczywistości kod podany przez Hilmara również działa, ponieważ tak długo, jak zero-pad do rozmiaru większego niż (w jego przypadku oblicza FT o rozmiarze ), „przesada” próbkę FT , i nadal otrzymujesz „przydatne” próbki (pozostałe powinny wynosić s). Tak więc, dla wydajności, od zera do , to wszystko, czego potrzebujesz (cóż, być może lepiej jest zerować do następnej potęgi 2 z , jeśli używasz FFT).N 2 N - 1 0 2 N - 1 2 N - $N-1$$N$$2N-1$$0$$2N-1$$2N-1$

W skrócie: powinieneś to zrobić (aby dostosować się do swojego języka programowania):

autocorr = ifft( complex( abs(fft(inputData, n=2*N-1))**2, 0 ) )

Lub w MATLAB:

autocorr = ifft(abs(fft(inputData, 2*N-1)).^2)

Głównym powodem, dla którego pożądane wyjście funkcji xcorr nie jest podobne do zastosowania funkcji FFT i IFFT, jest to, że podczas stosowania tych funkcji do sygnałów wynik końcowy jest zawijany kołowo .

Główną różnicę między zwijaniem liniowym i zwojem kołowym można znaleźć w zwojach liniowych i kołowych .

Problem można rozwiązać, wypełniając początkowo zerowanie sygnału i obcinając końcowe wyjście IFFT .