Algorytm Cooleya-Tukeya

Algorytm Cooleya-Tukeya

Algorytm Cooleya-Tukeya jest metodą szybkiej transformacji Fouriera (FFT), która przedstawia dyskretną transformację Fouriera (DFT) dla dowolnej złożonej wielkości N, wyrażonej jako:

N = N₁ × N₂,

gdzie człony N₁ i N₂ odpowiadają mniejszym DFT. Algorytm ten jest realizowany rekurencyjnie, co pozwala na zmniejszenie czasu obliczeń do O(N log N), szczególnie w przypadku, gdy N jest liczbą wysoce złożoną. Z uwagi na jego znaczenie, różne warianty i style implementacji tego algorytmu są znane pod różnymi nazwami, jak opisano poniżej. Nazwa algorytmu pochodzi od J.W. Cooleya i Johna Tukeya.

Algorytm Cooleya-Tukeya dzieli DFT na mniejsze DFT, co umożliwia łączenie go z innymi algorytmami DFT, takimi jak algorytm Rädera lub Bluesteina, które mogą obsługiwać duże czynniki pierwsze.

Historia

Algorytm, w tym jego rekurencyjne podejście, został wynaleziony około 1805 roku przez Carl Friedrich Gaussa, który stosował go do interpolacji trajektorii planetoid Pallas i Juno, jednak jego prace nie były szeroko znane aż do publikacji pośmiertnej. Gauss nie analizował jednak asymptotycznego czasu obliczeń. Formy algorytmu były na nowo odkrywane kilkukrotnie w XIX i XX wieku. Popularność transformacji FFT wzrosła po publikacji artykułu przez Jamesa Cooleya i Johna Tukeya w 1965 roku, w którym przedstawili oni metodę na użycie algorytmu na komputerach.

Uważa się, że Tukey wpadł na ten pomysł podczas posiedzenia amerykańskiej prezydenckiej komisji doradczej, badającej metody wykrywania testów broni jądrowej w ZSRR. Richard Garwin z IBM dostrzegł potencjał tej metody i skontaktował Tukeya z Cooleyem, sugerując, że metoda ta będzie przydatna w analizie danych krystalograficznych 3D. Ich wspólny artykuł zdobył dużą popularność.

Mimo że Gauss opisał ten sam algorytm, jego realizacja miała miejsce dopiero kilka lat po publikacji Cooleya i Tukeya. Ich artykuł odnosił się jedynie do pracy I.J. Gooda, znanej jako „algorytm FFT czynników pierwszych” (PFA), która jednak okazała się zupełnie innym algorytmem.

Przypadek Radix-2 DIT

Radix-2 decymacja w dziedzinie czasu (DIT) FFT jest najprostsza i najczęściej stosowana wersja algorytmu Cooleya-Tukeya. Radix-2 DIT dzieli DFT o rozmiarze N na dwie przeplatające się DFT, każda o rozmiarze N/2 w każdym etapie rekurencyjnym. DFT jest definiowana wzorem:

X_k = ∑_n=0^N−1 x_n e^{−(2πi/N)nk},

gdzie k jest liczbą całkowitą w zakresie od 0 do N−1.

Radix-2 DIT najpierw oblicza DFT dla danych wejściowych o indeksach parzystych x_2m oraz dla indeksów nieparzystych x_2m+1, a następnie łączy te dwa wyniki, aby uzyskać DFT całej sekwencji. Proces ten jest powtarzany rekurencyjnie, co redukuje czas do O(N log N).

Zakłada się, że N jest potęgą dwójki, co w praktyce jest często akceptowalne, ponieważ liczba punktów próbkowania N jest zazwyczaj ustalana przez aplikację.

Algorytm Radix-2 DIT dzieli DFT funkcji x_n na dwie części: sumę dla indeksów parzystych oraz sumę dla indeksów nieparzystych:

X_k = ∑_m=0^N/2−1 x_2m e^{−(2πi/N)(2m)k} + ∑_m=0^N/2−1 x_2m+1 e^{−(2πi/N)(2m+1)k}.

Można uwzględnić wspólny czynnik e^−(2πi/N)k z drugiej sumy. Wówczas obie sumy są DFT dla części o parzystych i nieparzystych indeksach funkcji x_n.

Pseudokod

W pseudokodzie proces ten można zapisać:

X0,...,N −1 ← ditfft2 (x, N, s):                     
    if N = 1 then
        X0 ← x0                                     
    else
        X0,...,N/2−1 ← ditfft2 (x, N/2, 2s)              
        XN/2,...,N −1 ← ditfft2 (x +s, N/2, 2s)            
        for k = 0 to N/2−1                          
            t ← Xk
            Xk ← t + exp(−2πi k/N) Xk + N/2
            Xk +N/2 ← t − exp(−2π i k/N) Xk +N/2
        endfor
    endif

Algorytm ditfft2 oblicza X = DFT(x) przez radix-2 DIT FFT, gdzie N jest całkowitą potęgą dwójki.

Algorytmy zmiany kolejności danych

Wiele technik reorganizacji danych istnieje w różnych implementacjach algorytmu Cooleya-Tukeya. Szczególnie interesujący jest problem opracowania algorytmu w miejscu, który nadpisuje dane wejściowe swoimi danymi wyjściowymi przy minimalnym użyciu pamięci.

Najbardziej znana technika to jawne odwrócenie bitów, które jest permutacją, w której dane o indeksie n przenoszone są do indeksu o odwróconych cyfrach. W przypadku algorytmu radix-2 DIT, dane wyjściowe są nadpisywane w miejscu na danych wejściowych, co wymaga zamiany wartości wyjściowych.

Opracowane zostały również różne alternatywne metody implementacji dla algorytmu Cooleya-Tukeya, które nie wymagają oddzielnego odwracania bitów i/lub angażowania dodatkowych permutacji w pośrednich fazach.

Bibliografia

James W. Cooley, John W. Tukey. An algorithm for the machine calculation of complex Fourier series. „Math. Comput.”. 19, s. 297–301, 1965. DOI: 10.2307/2003354.

Linki zewnętrzne

Prosty, pedagogiczny algorytm radix-2 Cooley-Tukeya FFT w C++.

KISSFFT – prosta mixed-radix implementacja Cooleya-Tukeya w C (open source).