Jakie są różnice między DFT i FFT, które sprawiają, że FFT jest tak szybki?

Próbuję zrozumieć FFT, oto co mam do tej pory:

Aby znaleźć wielkość częstotliwości w kształcie fali, należy je zbadać, mnożąc falę przez częstotliwość, której szukają, w dwóch różnych fazach (sin i cos) i uśredniając każdą z nich. Faza znajduje się w relacji do dwóch, a kod tego jest mniej więcej taki:

//simple pseudocode
var wave = [...];                //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...]  //all frequencies being tested for.  

function getMagnitudesOfSpectrum() {
   var magnitudesOut = [];
   var phasesOut = [];

   for(freq in spectrum) {
       var magnitudeSin = 0;
       var magnitudeCos = 0;

       for(sample in numSamples) {
          magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
          magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
       }

       magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
       phasesOut[freq] = //based off magnitudeSin and magnitudeCos
   }

   return magnitudesOut and phasesOut;
}

Aby to zrobić bardzo szybko dla bardzo wielu częstotliwości, FFT używają wielu sztuczek.

Jakie są sztuczki, dzięki którym FFT są znacznie szybsze niż DFT?

PS Próbowałem spojrzeć na kompletne algorytmy FFT w Internecie, ale wszystkie sztuczki są zwykle skondensowane w jeden piękny fragment kodu bez większego wyjaśnienia. Najpierw potrzebuję, zanim zrozumiem całą rzecz, jakieś wprowadzenie do każdej z tych skutecznych zmian jako koncepcji.

Dziękuję Ci.

Naiwna implementacja punktowego DFT to w zasadzie pomnożenie przez macierz N × N. Powoduje to złożoność O ( N 2 ) .$N$$N \times N$$\mathcal{O}(N^2)$

Jednym z najpopularniejszych algorytmów szybkiej transformaty Fouriera (FFT) jest algorytm FFT Radole-2 Cooley-Tukey Decimation-in-Time. To podstawowe podejście dziel i zwyciężaj.

Najpierw zdefiniuj „współczynnik twiddle” jako: gdziej≜

$$W_N \triangleq e^{-j\frac{2\pi}{N}}$$ jest jednostką urojoną, a DFTX[k]zx[n]jest dane przez X[k]= N - 1 Σ n = 0 x[n$j \triangleq \sqrt{-1}$$X[k]$$x[n]$ Jeśli N jest parzyste (i $$X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{kn} \, .$$ $N$ jest liczbą całkowitą), sumę można następnie podzielić na dwie sumy w następujący sposób X[k]= N / 2 - 1 ∑ n=0x[2n]W 2 k n N + N / 2 - 1 ∑ n=0x[2n+1]W k ( 2 n + 1 )$\tfrac{N}{2}$ $$X[k] = \sum_{n=0}^{N/2-1} x[2n]W_N^{2kn} + \sum_{n=0}^{N/2-1} x[2n+1]W_N^{k(2n+1)}$$ gdzie pierwsze sumowanie dotyczy parzystych próbek a drugie nieparzystych próbek x [ n ] . Zdefiniowanie x e [ n ] ≜ x [ 2 n ] i x o [ n ] ≜ x [ 2 n + 1 ] i wykorzystanie faktu, że$x[n]$$x[n]$$x_e[n] \triangleq x[2n]$$x_o[n] \triangleq x[2n+1]$

1. , i$W_N^{k(2n+1)} = W_N^{2kn}W_N^k$
2. ,$W_N^{2kn} = W_{N/2}^{kn}$

można to zapisać ponownie jako w którymXE[k]iXO[k],to

$$\begin{align} X[k] &= \sum_{n=0}^{N/2-1} x_e[n] W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1} x_o[n]W_{N/2}^{kn} \\ & = X_e[k] + W_N^k X_o[k] \end{align}$$ $X_e[k]$$X_o[k]$ punktowe transformaty DFT parzystych i nieparzystych próbek odpowiedniox[n]. Więc właśnie przekształciliśmy pojedynczypunktNDFT w dwa mniejsze$\tfrac{N}{2}$$x[n]$$N$ punktowe DFT. Zmniejsza to koszty obliczeniowe, ponieważ 2($\tfrac{N}{2}$ gdyN>2. $$2 \left( \frac{N}{2} \right)^2 + N < N^2$$ $N > 2$

Następnie możemy powtórzyć ten sam proces na tych dwóch mniejszych DFT. To podejście „dziel i rządź” pozwala osiągnąć złożoność , co jest znacznie lepsze niż O ( N 2 ), które mieliśmy z naiwną implementacją DFT (co dobrze ilustruje odpowiedź po lewej stronie ).$\mathcal{O}(N\log N)$$\mathcal{O}(N^2)$

http://nbviewer.jupyter.org/gist/leftaroundabout/83df89a7d3bdc24373ea470fb50be629

DFT, rozmiar 16

FFT, rozmiar 16

Różnica w złożoności jest dość oczywista, prawda?

---

Oto jak rozumiem FFT.

$\operatorname{FT} : \mathcal{L}^2(\mathbb{R}) \to \mathcal{L}^2(\mathbb{R})$

$\mathbb{R}\to\mathbb{C}$Najprostszym przypadkiem jest to, że twoja funkcja jest ciągła i dzielisz ją na tak małe regiony, że w zasadzie jest stała w każdym z nich. Następnie każda z STFT ma najsilniej zerowy wyraz. Jeśli zignorujesz (i tak rozkładające się) inne współczynniki, to każda domena będzie tylko jednym punktem danych. Ze wszystkich tych krótkoterminowych współczynników ograniczenia LF można przyjąć dyskretną transformatę Fouriera. W rzeczywistości to właśnie robisz, wykonując dowolną FT na zmierzonych rzeczywistych danych!

Zmierzone dane nie muszą jednak odpowiadać podstawowej wielkości fizycznej. Na przykład, gdy mierzysz natężenie światła , tak naprawdę mierzysz amplitudę fali elektromagnetycznej, której częstotliwość sama w sobie jest zbyt wysoka, aby można było próbkować za pomocą ADC. Ale wyraźnie można również obliczyć DFT próbkowanego sygnału natężenia światła i tanio, pomimo szalonej częstotliwości fali świetlnej.

Można to rozumieć jako najważniejszy powód, dla którego FFT jest tanie:

Nie zawracaj sobie głowy próbowaniem zobaczenia poszczególnych cykli oscylacyjnych z najwyższego poziomu. Zamiast tego przekształcaj tylko trochę informacji wysokiego poziomu, które zostały już wstępnie przetworzone lokalnie.

Ale to nie wszystko. Wspaniałą rzeczą w FFT jest to, że wciąż daje ci wszystkie informacje, które dałby pełny DFT . To znaczy wszystkie informacje, które uzyskasz również podczas próbkowania dokładnej fali elektromagnetycznej wiązki światła. Czy można to osiągnąć poprzez transformację sygnału fotodiody? - czy możesz zmierzyć z tego dokładną częstotliwość światła?

$\Delta \nu = 1/{\Delta t}$

Mając ogólnie dłuższy okres czasu, powinniśmy również być w stanie zawęzić niepewność częstotliwości. Jest to rzeczywiście możliwe, jeśli lokalnie mierzysz nie tylko częstotliwość szorstką, ale także fazę fali. Wiesz, że sygnał 1000 Hz będzie miał dokładnie tę samą fazę, jeśli spojrzysz na nią sekundę później. Podczas gdy sygnał 1000,5 Hz, chociaż jest nierozróżnialny w krótkiej skali, będzie miał odwróconą fazę sekundę później.

Na szczęście ta informacja o fazie może być bardzo dobrze przechowywana w pojedynczej liczbie zespolonej. I tak działa FFT! Zaczyna się od wielu małych, lokalnych przekształceń. Są tanie - z jednej strony oczywiście dlatego, że wykorzystują tylko niewielką ilość danych, ale po drugie dlatego, że wiedzą, że ze względu na krótki okres czasu nie są w stanie bardzo dokładnie ustalić częstotliwości - więc jest to nadal przystępne, nawet jeśli ty wykonaj wiele takich transformacji.

Rejestrują one jednak także fazę , dzięki czemu można dokładniej ustalić rozdzielczość częstotliwości na najwyższym poziomie. Wymagana transformacja jest znów tania, ponieważ sama nie przeszkadza w żadnych oscylacjach o wysokiej częstotliwości, ale tylko w przypadku wstępnie przetworzonych danych o niskiej częstotliwości.

---

^† _{Tak, moja argumentacja jest w tym momencie nieco okrągła. Nazwijmy to rekurencyjnym i nic nam nie jest ...}

^‡ _{Zależność ta nie jest mechaniką kwantową, lecz niepewnością Heisenberga ma ten sam podstawowy powód.}

$W_N^{nk}$$e^{j2\pi \frac{nk}{N}}$

Zwróć uwagę na pokazaną ścieżkę, a równanie poniżej pokazuje wynik dla przedziału częstotliwości X (1), zgodnie z równaniem Roberta.

Linie przerywane nie różnią się niczym od linii ciągłych, aby wyjaśnić, gdzie znajdują się połączenia sumowania.

zasadniczo, obliczając naiwny DFT bezpośrednio z podsumowania:

$$X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{j 2 \pi \frac{nk}{N}}$$

$N$$e^{j 2 \pi \frac{nk}{N}}$$N$$N-1$$X[k]$$k$$X[k+1]$ .

1. więc FFT przechowuje niektóre dane pośrednie.
2. FFT wykorzysta również faktoring współczynnika pośredniego nieco, aby ten sam współczynnik można zastosować do pośredniej kombinacji danych.

Jestem osobą wizualną. Wolę wyobrażać sobie FFT jako sztuczkę matrycową niż sztuczkę sumującą.

Aby wyjaśnić na wysokim poziomie:

Naiwny DFT oblicza każdą próbkę wyjściową niezależnie i wykorzystuje każdą próbkę wejściową w każdym obliczeniu (klasyczny algorytm N²).

Wspólna metoda FFT wykorzystuje symetrie i wzorce w definicji DFT, aby wykonać obliczenia w „warstwach” (warstwy N log), przy czym każda warstwa wymaga stałego czasu na próbkę, tworząc algorytm N log N.

Więcej szczegółów:

Jednym ze sposobów wizualizacji tych symetrii jest spojrzenie na DFT jako wejście macierzy 1 × N pomnożone przez macierz NxN wszystkich złożonych wykładników. Zacznijmy od przypadku „radix 2”. Podzielimy parzyste i nieparzyste wiersze macierzy (odpowiadające parzystym i nieparzystym próbkom wejściowym) i uznamy je za dwa oddzielne mnożenia macierzy, które sumują się, aby uzyskać ten sam końcowy wynik.

Spójrzmy teraz na te macierze: w pierwszej lewa połowa jest identyczna z prawą. Z drugiej strony prawa połowa to lewa połowa x -1. Oznacza to, że tak naprawdę musimy tylko użyć lewej połowy tych macierzy do pomnożenia i stworzyć prawą połowę tanio, mnożąc przez 1 lub -1. Następnie zauważ, że druga macierz różni się od pierwszej macierzy czynnikami, które są takie same w każdej kolumnie, więc możemy to wyliczyć i pomnożyć ją na wejściu, więc teraz zarówno parzyste, jak i nieparzyste próbki używają tej samej macierzy, ale wymagają mnożnika pierwszy. Ostatnim krokiem jest zaobserwowanie, że otrzymana macierz N / 2 × N / 2 jest identyczna z macierzą N / 2 DFT i możemy to robić wielokrotnie, aż dojdziemy do macierzy 1 × 1, gdzie DFT jest funkcją tożsamości.

Aby uogólnić poza podstawę 2, możesz spojrzeć na dzielenie co trzeci rząd i patrzeć na trzy fragmenty kolumn lub co 4 itd.

W przypadku danych wejściowych o podstawowej wielkości istnieje metoda poprawnego zerowania, FFT i obcięcia, ale jest to poza zakresem tej odpowiedzi.

Zobacz: http://whoiskylefinn.com/MatrixFFT.html

DFT dokonuje mnożenia macierzy brutalnej siły N ^ 2.

FFT wykonuje sprytne sztuczki, wykorzystując właściwości macierzy (zwyrodniając matrycę mnożąc) w celu zmniejszenia kosztów obliczeniowych.

Spójrzmy najpierw na mały DFT:

W = fft (oko (4));

x = rand (4,1) + 1j * rand (4,1);

X_ref = fft (x);

X = W * x;

aser (max (abs (X-X_ref)) <1e-7)

Świetnie, dlatego jesteśmy w stanie zastąpić wywołanie MATLAB do biblioteki FFTW małym mnożeniem macierzy 4x4 (złożonych) przez wypełnienie macierzy z funkcji FFT. Jak więc wygląda ta matryca?

N = 4

Wn = exp (-1j * 2 * pi / N),

f = ((0: N-1) '* (0: N-1))

f =

 0     0     0     0
 0     1     2     3
 0     2     4     6
 0     3     6     9

W = Wn. ^ F

W =

1 1 1 1

1 -i -1 i

1 -1 1 -1

1 i -1 -i

Każdy element to +1, -1, + 1j lub -1j. Oczywiście oznacza to, że możemy uniknąć pełnych złożonych multiplikacji. Co więcej, pierwsza kolumna jest identyczna, co oznacza, że ​​mnożymy pierwszy element x w kółko przez ten sam współczynnik.

Okazuje się, że produkty tensorowe Kroneckera, „współczynniki drgań” i macierz permutacji, w której indeks zmienia się zgodnie z odwróconą represantacją binarną, są zarówno zwarte, jak i dają alternatywne spojrzenie na sposób obliczania FFT jako zestawu rzadkich operacji macierzowych.

Poniższe wiersze to prosty FFT z dokładnością do dziesiętnej częstotliwości (DIF). Chociaż kroki mogą wydawać się niewygodne, wygodnie jest ponownie użyć FFT do przodu / do tyłu, radix4 / split-radix lub dziesiętnej w czasie, będąc jednocześnie rzetelną reprezentacją tego, w jaki sposób FFT na miejscu są zwykle wdrażane w prawdziwym świecie, Wierzę.

N = 4;

x = randn (N, 1) + 1j * randn (N, 1);

T1 = exp (-1j * 2 * pi * ([zera (1, N / 2), 0: (N / 2-1)]). '/ N),

M0 = kron (oko (2), fft (oko (2))),

M1 = kron (fft (oko (2)), oko (2)),

X = bitrevorder (x. '* M1 * diag (T1) * M0),

X_ref = fft (x)

aser (max (abs (X (:) - X_ref (:))) <1e-6)

CF Van Loan ma świetną książkę na ten temat.

Jeśli chcesz pić z Ognistego Węża Mądrości, proponuję:

„Szybkie transformacje - algorytmy, analizy, zastosowania” Douglas F. Elliott, K. Ramamohan Rao

Obejmuje FFT, Hartley, Winograd i aplikacje.

Jedną z mocnych stron jest to, że pokazuje, jak FFT jest zestawem rzadkich faktoryzacji macierzy z kolejnością odwracania bitów.