Algorytm Cooleya-Tukeya

Algorytm Cooleya-Tukeya – algorytm szybkiej transformacji Fouriera (FFT). Wyraża dyskretną transformację Fouriera (DFT) o dowolnej złożonej wielkości w członach mniejszych DFT wielkości i rekurencyjnie, w celu ograniczenia czasu obliczeń do szczególnie w przypadku będącego liczbą wysoce złożoną (liczbą gładką). Ze względu na znaczenie algorytmu, poszczególne warianty i style implementacji są znane pod własnymi nazwami, jak opisano poniżej. Algorytm został nazwany imieniem J.W. Cooleya i Johna Tukeya.

Ponieważ algorytm Cooleya-Tukeya rozbija DFT na mniejsze DFT, może być połączony arbitralnie z dowolnym innym algorytmem DFT. Na przykład algorytm Rädera lub Bluesteina może obsługiwać duże czynniki pierwsze, które nie mogą być rozkładane przez algorytm Cooley-Tukeya lub można wykorzystać algorytm czynników pierwszych dla większej wydajności wydzielania względnie pierwszych czynników.

Historia Algorytm ten, włączając jego rekurencyjne działanie, został wynaleziony około roku 1805 przez Carla Friedricha Gaussa, który używał go do interpolacji trajektorii planetoid Pallas i Juno, ale jego praca nie była powszechnie uznawana (opublikowana dopiero pośmiertnie w neo-łacinie)Gauss, Carl Friedrich, Theoria interpolationis methodo nova tractata, Werke, Band 3, s. 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866).M.T. Heideman, D.H. Johnson, C.S. Burrus, Gauss and the history of the fast Fourier transform, „IEEE ASSP Magazine”, 1(4), s. 14–21 (1984).. Jednak Gauss nie zajmował się analizą asymptotycznego czasu obliczeń. Różne ograniczone formy zostały ponownie odkryte kilka razy, zarówno w XIX w., jak i na początku XX wieku. Transformacje FFT stały się popularne po tym jak James Cooley (IBM) i John Tukey (Princeton) opublikowali w 1965 roku artykuł o ponownym wynalezieniu algorytmu i opisaniu metody na jego użycie na komputerze.

Potocznie uważa się, iż Tukey wpadł na ten pomysł podczas spotkania amerykańskiej prezydenckiej komisji doradczej, która obradowała nad sposobami wykrywania testów broni jądrowej w ZSRRRockmore, Daniel N., Comput. Sci. Eng. 2 (1), 60 (2000). The FFT – an algorithm the whole family can use Special issue on „top ten algorithms of the century” .. Inny uczestnik na tym spotkaniu, Richard Garwin z IBM, dostrzegł potencjał metody i skontaktował Tukeya z Cooleyem, jednocześnie mówiąc mu, iż ta metoda będzie mu potrzebna do analizy danych krystalograficznych 3D (patrz również: wielowymiarowe FFT). Następnie Cooley i Tukey opublikowali swój wspólny artykuł, który został szeroko przyjęty.

Mimo że Gauss opisał ten sam algorytm (chociaż nie analizując jego asymptotycznych kosztów), został on zrealizowany dopiero kilka lat po artykule Cooleya i Tukeya z 1965 roku. Ich artykuł cytuje jako inspirację jedynie pracę I.J. Good tym, obecnie nazywaną „algorytmem FFT czynników pierwszych” (PFA); chociaż początkowo błędnie sądzono, że algorytm Gooda jest odpowiednikiem algorytmu Cooleya-Tukeya, szybko zdano sobie sprawę, że PFA jest zupełnie innym algorytmem (działa jedynie dla rozmiarów które mają względnie pierwsze czynniki i opiera się na chińskim twierdzeniu o resztach, w przeciwieństwie do wspierania jakichkolwiek złożonych rozmiarów przez algorytm Cooleya-Tukeya)James W. Cooley, Peter A.W. Lewis, Peter W. Welch, Historical notes on the fast Fourier transform, „Proc. IEEE”, vol. 55 (no. 10), s. 1675–1677 (1967)..

Przypadek Radix-2 DIT Radix-2 decymacja w dziedzinie czasu (decimation-in-time) (DIT) FFT jest najprostszą i najbardziej powszechną formą algorytmu Cooleya-Tukeya, choć wysoce zoptymalizowane implementacje Cooleya-Tukeya zazwyczaj korzystają z innych form algorytmu jak opisano poniżej. Radix-2 DIT dzieli DFT o rozmiarze na dwie przeplatane DFT (stąd nazwa „radix-2”), o rozmiarze w każdym etapie rekurencyjnym.

Dyskretna transformacja Fouriera (DFT) jest określona wzorem:

gdzie jest liczbą całkowitą z zakresu od do

Radix-2 DIT najpierw oblicza DFT dane wejściowe o indeksach parzystych i indeksach nieparzystych a następnie łączy te dwa rezultaty do wytworzenia DFT całej sekwencji. Pomysł ten może następnie być przeprowadzony rekurencyjnie do redukcji ogólnego czasu do Ta uproszczona postać zakłada że jest potęgą dwójki; skoro liczba punktów próbki może być zwykle wybrana przez aplikację, to jest często nieważne ograniczenie.

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

Można uwzględniać wspólny mnożnik z drugiej sumy, jak przedstawiono w równaniu poniżej. Jest wtedy jasne że dwie sumy są DFT części o parzystych indeksach i DFT części o nieparzystych indeksach funkcji Oznaczmy DFT danych wejściowych o parzystych indeksach przez a DFT danych o nieparzystych indeksach przez i otrzymujemy:

Jednak te mniejsze DFT mają długość tak więc potrzebujemy obliczyć jedynie danych wyjściowych: dzięki własności okresowości DFT, dane wyjściowe dla z DFT długości są identyczne jak dane wyjściowe dla To znaczy, że i Czynnik fazowy jest posłuszny relacji: przerzucania znaku członu. Tak więc, całą DFT można obliczyć w następujący sposób:

Ten wynik, wyrażając DFT o długości rekurencyjnie w członach dwóch DFT wielkości jest rdzeniem szybkiej transformacji Fouriera radix-2 DIT. Algorytm zyskuje prędkość przez ponowne wykorzystanie wyników obliczeń pośrednich do obliczenia wielu danych wyjściowych DFT. Zauważmy, że ostateczne dane wyjściowe są otrzymane przez +/− kombinację i która jest po prostu rozmiarem-2 DFT (czasami zwana obliczeniami motylkowymi w tym kontekście); kiedy ona jest uogólniana do większych podstaw niż 2, poniżej, size-2 DFT jest zastąpione większym DFT (które samo może być obliczone przez FFT). thumb|300px|Diagram przepływu danych dla N = 8: decymacja w dziedzinie czasu radix-2 FFT rozdziela DFT długości N na dwie DFT długości następnie przez fazę łączącą składającej się z wielu DFT rozmiaru 2 zwanych operacjami „motylkowymi” (tzw. powodu kształtu diagramów przepływu danych).

Ten proces jest przykładem ogólnej techniki algorytmu dziel i zwyciężaj, w wielu tradycyjnych realizacjach jednak jawna rekurencja jest unikana i zamiast niej przechodzi się obliczeniowe drzewo wszerz.

Powyższe wyrażenie DFT rozmiaru jako dwie DFT rozmiaru jest czasami zwane lematem Danielsona-Lanczosa ponieważ ta tożsamość była zauważona przez tych dwóch autorów w 1942 rokuG.C. Danielson, C. Lanczos, Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids, „J. Franklin Inst.” 233, s. 365–380, 435–452 (1942). (na co wpływ miała praca C. Runge’a z 1903 roku). Zastosowali oni swój lemat we „wsteczny” rekurencyjny sposób, wielokrotnie „podwajając” rozmiar DFT aż spektrum transformaty skupiło się w jednym punkcie (choć najwyraźniej nie zdawali sobie sprawy z liniowo-logarytmicznej asymptotycznej złożoności którą osiągnęli). Praca Danielsona-Lanczosa wyprzedziła szeroką dostępność komputerów i wymagała obliczeń ręcznych (być może z pomocą mechaniczną, taką jak sumatory); relacjonowali czas obliczeń 140 minut dla rozmiaru DFT 64, działając na rzeczywistym wejściu do 3–5 znaczących cyfr. W artykule z 1965 roku Cooley i Tukey odnotowali czas pracy 0,02 minuty dla rozmiaru 2048 złożonej DFT na IBM 7094 (prawdopodobnie z 36-bitową pojedynczą precyzją, ~8 cyfr). Przeskalowując czas przez liczbę operacji, to odpowiada w przybliżeniu przyśpieszeniu o współczynnik około 800'000. (Aby umieścić w czasie obliczenia ręczne w perspektywie 140 minut dla rozmiaru 64 odpowiada to średnio co najwyżej 16 sekundom na każdą operację zmiennoprzecinkową, z których około 20% jest operacjami mnożenia).

Pseudokod W pseudokodzie powyższy proces można zapisać:

X0,...,N −1 ← ditfft2 (x, N, s): DFT of (x0, xs, x2s, ..., x(N -1)s): if N = 1 then X0 ← x0 przypadek bazowy trywialnego rozmiaru DFT = 1 else X0,...,N/2−1 ← ditfft2 (x, N/2, 2s) DFT (x0, x2s, x4s, ...) XN/2,...,N −1 ← ditfft2 (x +s, N/2, 2s) DFT (xs, xs +2s, xs +4s, ...) for k = 0 to N/2−1 łączenie DFT dwóch połówek do pełnego DFT: 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

W przedstawionym przykładzie, ditfft2 (x, N, 1), oblicza X =DFT(x) przez radix-2 DIT FFT (nie jest to algorytm typu in situ, gdzie jest całkowitą potęgą dwójki 2, a jest rozstawem danych wejściowej tablicy Oznaczenie oznacza tablicę zaczynająca się od

(Wyniki są w poprawnej kolejności w a dalsza permutacja odwracająca bit nie jest wymagana; często wspominana konieczność oddzielnej fazy odwracania bitów powstaje jedynie dla pewnych algorytmów działających in situ jak opisane poniżej.)

Wysokowydajne implementacje FFT dokonują wielu zmian w realizacji takiego algorytmu w porównaniu z przedstawionym pseudokodem. Na przykład można użyć większego przypadku bazowego niż do amortyzowania narzutu rekursji, czynniki fazowe można obliczyć wcześniej, a większe podstawy są często stosowane ze względu na pamięć cache; te i inne optymalizacje stosowane razem mogą poprawić wydajność o rząd wielkości lub więcejS.G. Johnson, M. Frigo, Implementing FFTs in practice, in Fast Fourier Transforms (C.S. Burrus, ed.), ch. 11, Rice University, Houston TX: Connexions, September 2008.. W wielu implementacjach podręcznikowych rekurencja depth-first jest całkowicie wyeliminowana na rzecz nierekurencyjnego podejścia wszerz, choć rekursja depth-first argumentuje za lepszą alokacją pamięci. Poniżej opisano szczegółowo kilka z wymienionych koncepcji.

Ogólne faktoryzacje thumb|500px|Podstawowy krok FFT Cooleya-Tukeya dla ogólnych faktoryzacji może być przedstawiony jako reinterpretacja jednowymiarowego DFT jako coś w rodzaju dwuwymiarowego DFT. Jednowymiarowa tablica wejściowa długości N = N1N2 jest reinterpretowana jako dwuwymiarowa macierz N1×N2 trzymana w pamięci w układzie kolumnowym. Wykonuje się mniejsze jednowymiarowe DFT wzdłuż kierunku N2 (nieprzylegająca kolejność), następnie mnoży przez czynniki fazowe, i na koniec wykonuje jednowymiarową DFT wzdłuż kierunku N1. Krok transpozycji może być wykonany w środku jak pokazano tutaj, lub na początku albo na końcu. To jest wykonywane rekurencyjnie dla mniejszych transformacji.

W ogólnym ujęciu, algorytmy Cooleya-Tukeya wyrażają w sposób rekurencyjny DFT złożonej wielkości jakoP. Duhamel, M. Vetterli, Fast Fourier transforms: a tutorial review and a state of the art, „Signal Processing” 19, s. 259–299 (1990).: Wykonaj DFT wielkości Pomnóż przez zespolony pierwiastek z jedynki zwany czynnikami fazowymi. Wykonaj DFT wielkości

Zazwyczaj lub jest niewielkim czynnikiem (niekoniecznie pierwszym), zwanym podstawą (radix) (który może się różnić pomiędzy etapami rekursji). Jeśli podstawą będzie wówczas mówimy o algorytmie decymacji w dziedzinie czasu (DIT). Gdy z kolei jest podstawą, wówczas mówimy o algorytmie decymacji w dziedzinie częstotliwości (DIF, zwanym także algorytmem Sande-Tukeya). Wersja przedstawiona powyżej była algorytmem radix-2 DIT; w końcowym wyrażeniu, faza mnożenia nieparzystej transformacji jest czynnikiem fazowym, a kombinacja „+/−” (obliczenia motylkowe) transformat parzystych i nieparzystych jest DFT rozmiaru 2. Podstawa małych DFT jest czasami określana mianem motyla z powodu kształtu schematu przepływu danych dla przypadku podstawy 2.