Skutecznie oblicza autokorelację za pomocą FFT

12

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.

— Goz
źródło

1

Bądź ostrożny. O ile dane nie są deterministyczne, zazwyczaj możemy jedynie oszacować sekwencję autokorelacji. Istnieją dwie popularne wersje szacunków autokorelacji: stronnicze i obiektywne. Bezstronne wyniki w szacunkach autokorelacji, które są statystycznie bezstronne. Jednak wariancja może być bardzo duża dla opóźnień wysokiego rzędu, co prowadzi do problemów, jeśli oszacowanie autokorelacji jest stosowane na przykład w odwróceniu macierzy. Tendencyjne próbki wykazują błąd statystyczny, ale z mniejszą wariancją (i średnim błędem kwadratowym). Oba są statystycznie spójne. Powyżej masz nietypowe, tendencyjne oszacowanie.

— Bryan

16

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)));

— Hilmar
źródło

Wow, które działało pięknie. Prosta wersja C (jednowątkowa, bez SIMD) mojego trackera pitch działała w 0,8 sekundy przy użyciu powyższej metody, w przeciwieństwie do wersji opartej na prymitywie inteligencji, która działała w 0,4 sekundy. To jest wspaniałe! Dzięki

— Goz,

7

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$

— fenenety
źródło

3

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ć. $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. $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-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)

— Jean-Louis Durrieu
źródło

0

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 .

— Venu
źródło