Tekst źródłowy strony Algorytm Cooleya-Tukeya

{{Algorytm infobox
 | nazwa            = Algorytm Cooleya-Tukeya
 | czas             = <math>O(n \log n)</math>
}}
'''Algorytm Cooleya-Tukeya''' – algorytm [[Szybka transformacja Fouriera|szybkiej transformacji Fouriera]] (FFT). Wyraża [[Dyskretna transformata Fouriera|dyskretną transformację Fouriera]] (DFT) o dowolnej [[Liczby złożone|złożonej]] wielkości <math>N = N_1 N_2</math> w członach mniejszych DFT wielkości <math>N_1</math> i <math>N_2,</math> [[Rekurencja|rekurencyjnie]], w celu ograniczenia czasu obliczeń do <math>\text{O}(N \log N),</math> szczególnie w przypadku <math>N</math> będącego liczbą wysoce złożoną ([[Liczby gładkie|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 [[James Cooley|J.W. Cooleya]] i [[John Tukey|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 [[Czynnik pierwszy|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 [[Carl Friedrich Gauss|Carla Friedricha Gaussa]], który używał go do interpolacji trajektorii [[Planetoida|planetoid]] [[(2) Pallas|Pallas]] i [[(3) Juno|Juno]], ale jego praca nie była powszechnie uznawana (opublikowana dopiero pośmiertnie w neo-łacinie)<ref>Gauss, Carl Friedrich, [http://lseet.univ-tln.fr/~iaroslav/Gauss_Theoria_interpolationis_methodo_nova_tractata.php ''Theoria interpolationis methodo nova tractata''], Werke, Band 3, s. 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866).</ref><ref name="Heideman84">M.T. Heideman, D.H. Johnson, C.S. Burrus, [http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1162257 ''Gauss and the history of the fast Fourier transform''], „IEEE ASSP Magazine”, 1(4), s. 14–21 (1984).</ref>. 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<ref name="Heideman84" />. Transformacje FFT stały się popularne po tym jak James Cooley (IBM) i John Tukey ([[Uniwersytet w Princeton|Princeton]]) opublikowali w 1965 roku artykuł o ponownym wynalezieniu algorytmu i opisaniu metody na jego użycie na komputerze<ref name="CooleyTukey65">{{cytuj pismo |nazwisko=Cooley |imię=James W. |nazwisko2=Tukey |imię2=John W. |tytuł=An algorithm for the machine calculation of complex Fourier series |czasopismo=[[Mathematics of Computation|Math. Comput.]] |wolumin=19 |strony=297–301 |rok=1965 |doi=10.2307/2003354}}</ref>.

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 [[Broń jądrowa|broni jądrowej]] w [[Związek Socjalistycznych Republik Radzieckich|ZSRR]]<ref>{{Cytuj pismo |nazwisko=Cooley |imię=James W. |nazwisko2=Lewis |imię2=Peter A.W. |nazwisko3=Welch |imię3=Peter D. |tytuł=Historical notes on the fast Fourier transform |czasopismo=IEEE Trans. on Audio and Electroacoustics |wolumin=15 |wydanie=2 |strony=76–79 |rok=1967}}</ref><ref>Rockmore, Daniel N., ''Comput. '' ''Sci. Eng.'' '''2''' (1), 60 (2000). [http://www.cs.dartmouth.edu/~rockmore/cse-fft.pdf The FFT – an algorithm the whole family can use] Special issue on „top ten algorithms of the century” [https://web.archive.org/web/20090407080904/http://amath.colorado.edu/resources/archive/topten.pdf].</ref>. 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<ref name="Heideman84" />. Ich artykuł cytuje jako inspirację jedynie pracę I.J. Good tym, obecnie nazywaną „algorytmem FFT czynników pierwszych” (PFA)<ref name="CooleyTukey65" />; 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)<ref>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).</ref>.

== Przypadek Radix-2 DIT ==
'''Radix-2''' [[Decymacja (cyfrowe przetwarzanie sygnałów)|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 <math>N</math> na dwie przeplatane DFT (stąd nazwa „radix-2”), o rozmiarze <math>N/2</math> w każdym etapie rekurencyjnym.

Dyskretna [[transformacja Fouriera]] (DFT) jest określona wzorem:
: <math>X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk},</math>

gdzie <math>k</math> jest liczbą całkowitą z zakresu od <math>0</math> do <math>N-1.</math>

Radix-2 DIT najpierw oblicza DFT dane wejściowe o indeksach parzystych <math>x_{2m} (x_0, x_2, \dots, x_{N-2})</math> i indeksach nieparzystych <math>x_{2m+1} (x_1, x_3, \dots, x_{N-1}),</math> a następnie łączy te dwa rezultaty do wytworzenia DFT całej sekwencji. Pomysł ten może następnie być przeprowadzony [[Rekurencja|rekurencyjnie]] do redukcji ogólnego czasu do <math>\text{O}(N \log N).</math> Ta uproszczona postać zakłada że <math>N</math> jest [[Potęga dwójki|potęgą dwójki]]; skoro liczba punktów próbki <math>N</math> może być zwykle wybrana przez aplikację, to jest często nieważne ograniczenie.

Algorytm Radix-2 DIT przestawia DFT funkcji <math>x_n</math> na dwie części: sumę dla indeksów parzystych <math>n={2m}</math> i sumę dla indeksów nieparzystych <math>n={2m+1}{:}</math>
: <math>\begin{matrix} X_k & =
& \sum \limits_{m=0}^{N/2-1} x_{2m} e^{-\frac{2\pi i}{N} (2m)k}
+ \sum \limits_{m=0}^{N/2-1} x_{2m+1} e^{-\frac{2\pi i}{N} (2m+1)k}.
\end{matrix}</math>

Można uwzględniać wspólny mnożnik <math>e^{-\frac{2\pi i}{N}k}</math> z drugiej sumy, jak przedstawiono w równaniu poniżej. Jest wtedy jasne że dwie sumy są DFT części o parzystych indeksach <math>x_{2m}</math> i DFT części o nieparzystych indeksach <math>x_{2m+1}</math> funkcji <math>x_n.</math> Oznaczmy DFT danych wejściowych o parzystych indeksach <math>x_{2m}</math> przez <math>E_k</math> a DFT danych o nieparzystych indeksach <math>x_{2m + 1}</math> przez <math>O_k</math> i otrzymujemy:
: <math>X_k
= \underbrace{\sum \limits_{m=0}^{N/2-1} x_{2m} e^{-\frac{2\pi i}{N/2} mk}}_{\mathrm{DFT\;element\dot{o}w\;parzystych\;} x_m}{}
+ e^{-\frac{2\pi i}{N}k}
\underbrace{\sum \limits_{m=0}^{N/2-1} x_{2m+1} e^{-\frac{2\pi i}{N/2} mk}}_{\mathrm{DFT\;element\dot{o}w\;nieparzystych\;} x_m}
= E_k + e^{-\frac{2\pi i}{N}k} O_k.</math>

Jednak te mniejsze DFT mają długość <math>N/2,</math> tak więc potrzebujemy obliczyć jedynie <math>N/2</math> danych wyjściowych: dzięki własności okresowości DFT, dane wyjściowe dla <math>N/2 \leqslant k</math> z DFT długości <math>N/2</math> są identyczne jak dane wyjściowe dla <math>0\leqslant k.</math> To znaczy, że <math>E_{k + N/2} = E_k</math> i <math>O_{k + N/2} = O_k.</math> Czynnik fazowy <math>\exp[-2\pi i k/N]</math> jest posłuszny relacji: <math>\exp[-2\pi i (k + N/2)/N] = e^{-\pi i} \exp[-2\pi i k/N] = -\exp[-2\pi i k/N],</math> przerzucania znaku <math>O_{k + N/2}</math> członu. Tak więc, całą DFT można obliczyć w następujący sposób:
: <math>X_k = \left\{\begin{matrix}
E_k + e^{-\frac{2\pi i}{N}k} O_k & \text{if } k < N/2 \\ \\
E_{k-N/2} - e^{-\frac{2\pi i}{N} (k-N/2)} O_{k-N/2} & \text{if }
k \geqslant N/2.
\end{matrix} \right.</math>

Ten wynik, wyrażając DFT o długości <math>N</math> rekurencyjnie w członach dwóch DFT wielkości <math>N/2,</math> 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ę <math>E_k</math> i <math>O_k \exp(-2\pi i k/N),</math> 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).
[[Plik:DIT-FFT-butterfly.png|thumb|300px|Diagram przepływu danych dla N = 8: [[Decymacja (cyfrowe przetwarzanie sygnałów)|decymacja]] w dziedzinie czasu radix-2 FFT rozdziela DFT długości N na dwie DFT długości <math>N/2</math> 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 [[Dziel i zwyciężaj|algorytmu dziel i zwyciężaj]], w wielu tradycyjnych realizacjach jednak jawna rekurencja jest unikana i zamiast niej przechodzi się obliczeniowe drzewo [[Przeszukiwanie wszerz|wszerz]].

Powyższe wyrażenie DFT rozmiaru <math>N</math> jako dwie DFT rozmiaru <math>N/2</math> jest czasami zwane [[lemat]]em '''Danielsona-[[Cornelius Lanczos|Lanczosa]]''' ponieważ ta tożsamość była zauważona przez tych dwóch autorów w 1942 roku<ref>G.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).</ref> (na co wpływ miała praca [[Carl David Tolmé Runge|C. Runge’a]] z 1903 roku<ref name="Heideman84" />). 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 [[Sumator (maszyna biurowa)|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ą [[Liczba zmiennoprzecinkowa|pojedynczą precyzją]], ~8 cyfr)<ref name="CooleyTukey65" />. 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 [[pseudokod]]zie powyższy proces można zapisać<ref name="Johnson08" />:

 ''X''<sub>0,...,''N'' −1</sub> ← '''ditfft2''' (''x'', ''N'', ''s''):                     ''DFT of (x''<sub>0</sub>, ''x''<sub>''s''</sub>, ''x''<sub>2''s''</sub>, ..., ''x''<sub>(''N'' -1)''s''</sub>):
     if ''N''  = 1 then
         ''X''<sub>0</sub> ← ''x''<sub>0</sub>                                      '''przypadek bazowy trywialnego rozmiaru DFT  = 1'''
     else
         ''X''<sub>0,...,''N''/2−1</sub> ← '''ditfft2''' (''x'', ''N''/2, 2''s'')              ''DFT (x''<sub>0</sub>, ''x''<sub>2''s''</sub>, ''x''<sub>4''s''</sub>, ...)
         ''X''<sub>''N''/2,...,''N'' −1</sub> ← '''ditfft2''' (''x'' +s, ''N''/2, 2''s'')            ''DFT (x''<sub>''s''</sub>, ''x''<sub>''s'' +2''s''</sub>, ''x''<sub>''s'' +4''s''</sub>, ...)
         for ''k''  = 0 to ''N''/2−1                           ''łączenie DFT dwóch połówek do pełnego DFT:''
             t ← ''X''<sub>''k''</sub>
             ''X''<sub>''k''</sub> ← t + exp(−2π''i''  ''k''/''N'') ''X''<sub>''k'' + ''N''/2</sub>
             ''X''<sub>''k'' +''N''/2</sub> ← t − exp(−2π ''i''  ''k''/''N'') ''X''<sub>''k'' +''N''/2</sub>
         endfor
     endif

W przedstawionym przykładzie, <code>'''ditfft2'''</code> (''x'', ''N'', 1), oblicza ''X'' =DFT(''x'') przez radix-2 DIT FFT (nie jest to algorytm typu [[Algorytm in situ|in situ]], gdzie <math>N</math> jest całkowitą potęgą dwójki 2, a <math>s = 1</math> jest rozstawem danych wejściowej [[Tablica (informatyka)|tablicy]] <math>x.</math> Oznaczenie <math>x + s</math> oznacza tablicę zaczynająca się od <math>x_s.</math>

(Wyniki są w poprawnej kolejności w <math>X,</math> 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ż <math>N=1</math> do amortyzowania narzutu rekursji, czynniki fazowe <math>\exp[-2\pi i k/N]</math> 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ęcej<ref name="Johnson08">S.G. Johnson, M. Frigo, [http://cnx.org/content/m16336/ ''Implementing FFTs in practice''], in ''Fast Fourier Transforms'' (C.S. Burrus, ed.), ch. 11, Rice University, Houston TX: Connexions, September 2008.</ref>. W wielu implementacjach podręcznikowych rekurencja depth-first jest całkowicie wyeliminowana na rzecz nierekurencyjnego podejścia wszerz, choć rekursja [[Przeszukiwanie w głąb|depth-first]] argumentuje za lepszą alokacją pamięci<ref name="Johnson08" /><ref name="Singleton67" />. Poniżej opisano szczegółowo kilka z wymienionych koncepcji.

== Ogólne faktoryzacje ==
[[Plik:Cooley-tukey-general.png|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 <math>N=N_1 N_2</math> jako<ref name="DuhamelVe90">P. Duhamel, M. Vetterli, ''Fast Fourier transforms: a tutorial review and a state of the art'', „Signal Processing” '''19''', s. 259–299 (1990).</ref>:
# Wykonaj <math>N_1</math> DFT wielkości <math>N_2.</math>
# Pomnóż przez zespolony pierwiastek z jedynki zwany czynnikami fazowymi.
# Wykonaj <math>N_2</math> DFT wielkości <math>N_1.</math>

Zazwyczaj <math>N_1</math> lub <math>N_2</math> jest niewielkim czynnikiem (niekoniecznie ''pierwszym''), zwanym '''podstawą''' (radix) (który może się różnić pomiędzy etapami rekursji). Jeśli podstawą będzie <math>N_1,</math> wówczas mówimy o algorytmie '''[[Decymacja (cyfrowe przetwarzanie sygnałów)|decymacji]] w dziedzinie czasu''' (DIT). Gdy z kolei <math>N_2</math> 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.

Istnieje wiele innych odmian algorytmu Cooley-Tukeya. Implementacje '''mixed-radix''' obsługują złożone rozmiary z różnymi (zwykle małymi) czynnikami oprócz dwójki, zwykle (ale nie zawsze) stosujące algorytm złożoności <math>\text{O}(N_1)</math> dla pierwszych przypadków bazowych rekurencji (można także zastosować algorytm o złożoności <math>N\log N</math> dla przypadków bazowych będących liczbami pierwszymi, taki jak algorytmy Radera lub Bluesteina). Algorytm split-radix łączy podstawy 2 i 4, z tym, że pierwsza transformacja podstawy 2 nie wymaga czynnika fazowego w celu osiągnięcia najmniejszej znanej liczby operacji dla rozmiarów potęgi dwójki<ref name="DuhamelVe90" />, chociaż ostatnie wariacje osiągnęły nawet niższą wartość<ref>{{Cytuj |autor=T. Lundy, J. Van Buskirk |tytuł=A new matrix approach to real FFTs and convolutions of length 2<sup><em>k</em></sup> |czasopismo=Computing |numer=80 |s=23–45 |data=2007}}</ref><ref>S.G. Johnson, M. Frigo, ''[http://www.fftw.org/newsplit.pdf A modified split-radix FFT with fewer arithmetic operations]'', „IEEE Trans. Signal Processing” '''55''' (1), s. 111–119 (2007).</ref>. Na komputerach współczesnych, wydajność jest określana przez większą ilość [[Pamięć podręczna procesora|pamięci podręcznej]] i kwestie potokowości procesora niż ścisłe obliczanie liczby wykonanych operacji; dobrze zoptymalizowane implementacje FFT często używają większych podstaw i/lub zakodowanych na stałe transformat o znaczących wielkościach<ref name="FrigoJohnson05" />). Inną cechą algorytmu Cooleya-Tukeya jest to, że wyraża rozmiar <math>N</math> jednowymiarowego DFT jako <math>N_1</math> przez <math>N_2</math> dwuwymiarową DFT (plus czynniki fazowe), gdzie macierz wyjściowa jest [[Macierz transponowana|transponowana]]. Wynik tych wszystkich transpozycji, dla algorytmu radix-2, odpowiada odwróceniu bitowemu indeksów wejściowych (DIF) lub wyjściowych (DIT). Jeżeli zamiast małej podstawy używa się podstawy o wielkości rzędu <math>\sqrt N</math> i jawnej transpozycji matrycy wejścia/wyjścia, mówimy o '''czteroetapowym''' algorytmie (lub ''sześciokrokowym'', w zależności od liczby transpozycji). Taki algorytm wstępnie został zaproponowany w celu poprawy alokacji pamięci<ref name="GenSande66">W.M. Gentleman, G. Sande, ''Fast Fourier transforms – for fun and profit'', „Proc. AFIPS” '''29''', 563–578 (1966).</ref><ref name="Bailey90">Bailey, David H., ''FFTs in external or hierarchical memory'', „J. Supercomputing” '''4''' (1), s. 23–35 (1990).</ref> (np. do optymalizacji pamięci podręcznej lub „operacji pozardzeniowych”, a w późniejszym czasie okazał się on algorytmem optymalnym, niezależnym od rozmiaru cache)<ref name="Frigo99">M. Frigo, C.E. Leiserson, H. Prokop, S. Ramachandran, ''Cache-oblivious algorithms'' in ''Proceedings of the 40th IEEE Symposium on Foundations of Computer Science'' (FOCS 99), s. 285–297. 1999. [http://ieeexplore.ieee.org/iel5/6604/17631/00814600.pdf?arnumber=814600 Extended abstract at IEEE], [http://citeseer.ist.psu.edu/307799.html at Citeseer].</ref>.

Ogólna faktoryzacja Cooleya-Tukeya przepisuje indeksy <math>k</math> oraz <math>n</math> w postaci <math>k = N_2 k_1 + k_2</math> i <math>n = N_1 n_2 + n_1,</math> gdzie indeksy <math>k</math> i <math>n</math> zaczynają się od <math>0\dots N_a - 1</math> (dla <math>a</math> 1 lub 2). Oznacza to, że reindeksuje ona wejście <math>(n)</math> i wyjście <math>(k)</math> jako <math>N_1</math> przez <math>N_2</math> dwuwymiarowe tablice odpowiednio w układzie kolumnowym i układzie rzędowym; różnica pomiędzy tym indeksowaniem jest transpozycją, wspomnianą powyżej. Kiedy to przeindeksowanie jest podstawione do równania DFT dla <math>nk,</math> <math>N_1 n_2 N_2 k_1</math> krzyż członów zanika (jego wykładnikiem jest jedność), a pozostałe człony dają

: <math>\begin{align}
X_{N_2 k_1 + k_2}
&= \sum_{n_1=0}^{N_1-1} \sum_{n_2=0}^{N_2-1} x_{N_1 n_2 + n_1}
   e^{-\frac{2\pi i}{N_1 N_2} \cdot (N_1 n_2 + n_1) \cdot (N_2 k_1 + k_2) } \\
&= \sum_{n_1=0}^{N_1-1}
   \left[ e^{-\frac{2\pi i}{N} n_1 k_2} \right]
   \left( \sum_{n_2=0}^{N_2-1} x_{N_1 n_2 + n_1}
   e^{-\frac{2\pi i}{N_2} n_2 k_2} \right)
   e^{-\frac{2\pi i}{N_1} n_1 k_1},
\end{align}</math>

gdzie każda wewnętrzna suma jest DFT rozmiaru <math>N_2,</math> każda zewnętrzna suma jest DFT rozmiaru <math>N_1,</math> a człony w nawiasie [...] są czynnikiem fazowym.

Arbitralna podstawa <math>r</math> (jak również mieszane podstawy) może być stosowana tak jak wykazali zarówno Cooley z Tukeyem<ref name="CooleyTukey65" />, jak i Gauss (który przedstawił przykłady kroków radix-3 i radix-6)<ref name="Heideman84" />. Cooley i Tukey pierwotnie przyjmowali że obliczenia motylkowe podstawy wymagają <math>\text{O}(r^2)</math> pracy i stąd obliczyli, iż złożoność dla podstawy <math>r</math> wynosi <math>\text{O}(r^2 N/r \log_r N) = \text{O}(N\log_2 (N) r/\log_2 r);</math> z liczenia wartości <math>r/\log_2 r</math> dla całkowitych wartości <math>r</math> od 2 do 12 odkryto, że optymalna podstawa wynosi 3 (najbliższa liczba całkowita do ''[[Podstawa logarytmu naturalnego|e]]'', która minimalizuje <math>r/\log_2 r</math>)<ref name="CooleyTukey65" /><ref>J.W. Cooley, P. Lewis, P. Welch, ''The Fast Fourier Transform and its Applications'', „IEEE Trans on Education” '''12''', 1, s. 28–34 (1969).</ref>. Ta analiza była błędna, chociaż podstawa obliczeń motylkowych jest również DFT i może być wykonana przez algorytm FFT w <math>\text{O}(r \log r)</math> operacjach, stąd podstawa <math>r</math> faktycznie anuluje złożoność <math>\text{O}(r\log(r) N/r \log_r N),</math> a optymalne <math>r</math> jest określone przez bardziej skomplikowane rozważania. W praktyce dość duże <math>r</math> (o wartości 32 lub 64) są potrzebne w celu skutecznego wykorzystania np. dużej liczby [[Rejestr procesora|rejestrów procesora]] na nowoczesnych procesorach<ref name="FrigoJohnson05" />, a nawet nieograniczona podstawa <math>r = \sqrt N</math> również osiąga złożoność <math>\text{O}(N \log N)</math> i wykazuje teoretyczne i praktyczne korzyści dla dużych <math>N,</math> jak wspomniano powyżej<ref name="GenSande66" /><ref name="Bailey90" /><ref name="Frigo99" />.

== Algorytmy zmiany kolejności danych, odwracania bitów i obliczeń w miejscu ==
Choć powyższe streszczenie [[Rozkład na czynniki|faktoryzacji]] DFT Cooley-Tukeya może być zastosowane w pewnej formie do wszystkich implementacji algorytmu, znacznie większa różnorodność istnieje w technikach porządkowania i dostępu do danych na każdym etapie FFT. Szczególnie interesujący jest problem opracowania algorytmu [[Algorytm in situ|w miejscu]], który nadpisuje dane wejściowe swoimi danymi wyjściowymi przy użyciu dodatkowej pamięci o wartości zaledwie O(1).

Najbardziej znana technika przeorganizowania obejmuje jawne '''odwrócenie bitu''' dla algorytmów radix-2 ''in situ''. Odwrócenie bitów jest [[Permutacja|permutacją]], gdzie dana o indeksie <math>n,</math> zapisana [[Dwójkowy system liczbowy|binarnie]] cyframi <math>b_4,</math> <math>b_3,</math> <math>b_2,</math> <math>b_1,</math> <math>b_0</math> (tj. 5 cyfr dla <math>N=32</math> danych wejściowych), jest przeniesiona do indeksu o odwróconych cyfrach <math>b_0,</math> <math>b_1,</math> <math>b_2,</math> <math>b_3,</math> <math>b_4.</math> Rozważmy ostatnią fazę algorytmu radix-2 DIT podobnego do tego prezentowanego powyżej, gdzie dane wyjściowe są nadpisane w miejscu na danych wejściowych: kiedy <math>E_k</math> i <math>O_k</math> są połączone z size-2 DFT, te dwie wartości są nadpisane przez dane wyjściowe. Jednak, dwie wartości wyjściowe powinny znaleźć się w pierwszej i drugiej ''połowie'' tablicy wyjściowej, co odpowiada ''najstarszemu'' bitowi <math>b_4</math> (dla <math>N=32</math>), podczas gdy dwie dane wejściowe <math>E_k</math> i <math>O_k</math> są przeplatane parzystymi i nieparzystymi elementami, co odpowiada ''najmłodszemu'' bitowi <math>b_0.</math> Tak więc, w celu uzyskania danych wyjściowych w odpowiednim miejscu, te dwa bity muszą być wymieniane. Jeśli na wszystkich rekurencyjnych fazach jest zawarty algorytm radix-2 DIT, ''wszystkie'' bity muszą być zamienione i w ten sposób trzeba przetwarzać wstępnie dane wejściowe („lub” przetwarzać na koniec dane wyjściowe) z odwróceniem bitu aby otrzymać dane wyjściowe w porządku wejściowym<ref>Jeżeli każda podtransformacja rozmiaru <math>N/2</math> operuje na przyległych danych, ''dane wejściowe'' DIT są wstępnie przetworzone przez odwrócenie bitów.</ref>. Analogicznie, jeżeli wszystkie czynności zostaną wykonane w odwrotnej kolejności, wynikiem będzie algorytm radix-2 DIF z bitowym odwróceniem w postprocessingu (lub odpowiednio w preprocessingu). Alternatywnie, niektóre zastosowania (takie jak splot) działają równie dobrze na danych z bitami odwróconymi, więc można wykonać transformacje w przód, przetwarzanie, a następnie odwrócić wszystkie transformaty bez odwracania bitów aby otrzymać końcowy wynik w naturalnym porządku.

Niemniej jednak wielu użytkowników FFT preferuje dane wyjściowe w naturalnym porządku, a oddzielna, jawna faza odwrócenia bitów może mieć istotny wpływ na czas obliczeń<ref name="FrigoJohnson05" /> (choć odwrócenie bitów można wykonać w czasie <math>\text{O}(n),</math> co było przedmiotem wielu badań)<ref>{{cytuj pismo |nazwisko=Karp |imię=Alan H. |tytuł=Bit reversal on uniprocessors |czasopismo=SIAM Review |wolumin=38 |wydanie=1 |strony=1–26 |rok=1996 |doi=10.1137/1038001}}</ref><ref>{{Cytuj pismo|nazwisko=Carter |imię=Larry |nazwisko2=Gatlin |imię2=Kang Su |tytuł=Towards an optimal bit-reversal permutation program |czasopismo=Proc. 39th Ann. Symp. on Found. of Comp. Sci. (FOCS) |strony=544–553 |rok=1998 |doi=10.1109/SFCS.1998.743505}}</ref><ref>{{Cytuj pismo |nazwisko=Rubio |imię=M. |nazwisko2=Gómez |imię2=P. |nazwisko3=Drouiche |imię3=K. |tytuł=A new superfast bit reversal algorithm |czasopismo=Intl. J. Adaptive Control and Signal Processing |wolumin=16 |wydanie=10 |strony=703–707 |rok=2002 |doi=10.1002/acs.718}}</ref>. Ponadto podczas gdy permutacja jest odwróceniem bitów w przypadku radix-2, jest to przeważnie arbitralnym (mixed-base) odwróceniem cyfr dla przypadku mixed-radix, a algorytmy permutacji stają się bardziej skomplikowane do zaimplementowania. Co więcej, na wielu architekturach sprzętowych zalecane jest uporządkowanie pośrednich faz algorytmu FFT tak, by one operowały na kolejnych (lub przynajmniej na lepiej zaalokowanych) elementach danych. W tym celu, opracowano wiele alternatywnych metod implementacji dla algorytmu Cooleya-Tukeya, które nie wymagają oddzielnego odwracania bitów i/lub angażowania dodatkowych permutacji w pośrednich fazach.

Problem jest znacznie uproszczony, jeśli roważymy algorytm '''nie będący in situ:''' tablica wyjściowa różni się od tablicy wejściowej lub, równoważnie, dostępna jest pomocnicza tablica o równym rozmiarze. '''Alorytm Stockhama automatycznego sortowania'''<ref>{{Cytuj pismo |nazwisko=Stockham |imię=T.G. |tytuł=High speed convolution and correlation |czasopismo=Spring Joint Computer Conference, Proc. AFIPS |wolumin=28 |wydanie= |strony=229–233 |rok=1966}}</ref> wykonuje każdy etap FFT nie na miejscu, zazwyczaj zapisując dane „w tył i w przód” pomiędzy dwiema tablicami, transponując jedną „cyfrę” indeksów w każdym etapie, co było szczególnie popularne w architekturze [[SIMD]]<ref>{{cytuj książkę |nazwisko=Swarztrauber |imię=P.N. |tytuł=Parallel Computations |url=https://archive.org/details/parallelcomputat0000unse |inni=G. Rodrigue |rozdział=Vectorizing the FFTs |wydawca=Academic Press |miejsce=New York |rok=1982 |strony=[https://archive.org/details/parallelcomputat0000unse/page/51 51]–83 |isbn=0-12-592101-2}}</ref>. Jeszcze większy potencjał zalet SIMD (więcej kolejnych dostępów) zostały zaproponowany dla algorytmu '''Pease'''<ref>{{Cytuj pismo |nazwisko=Pease |imię=M.C. |tytuł=An adaptation of the fast Fourier transform for parallel processing |czasopismo=J. ACM |wolumin=15 |wydanie=2 |strony=252–264 |rok=1968 |doi=10.1145/321450.321457}}</ref>, który również zmienia kolejność metodą „nie w miejscu” w każdej fazie, lecz wymaga ona odrębnego odwrócenia bitu/cyfry i <math>\text{O}(N \log N)</math> pamięci. Można również bezpośrednio zastosować definicję faktoryzacji Cooleya-Tukeya z jawną rekurencją ([[Przeszukiwanie w głąb|depth-first]]) i małymi podstawami, w wyniku której otrzymujemy dane wyjściowe w naturalnym porządku „nie na miejscu” bez osobnego kroku permutacji (jak w pseudokodzie powyżej). Jednak fakt iż będą korzyści związane z alokowaniem niezależnym od wielkości cache na systemach z hierarchiczną pamięcią, może być podważony<ref name="Singleton67">{{Cytuj pismo |nazwisko=Singleton |imię=Richard C. |tytuł=On computing the fast Fourier transform |czasopismo=Commun. of the ACM |wolumin=10 |wydanie=10 |rok=1967 |strony=647–654 |doi=10.1145/363717.363771}}</ref><ref name="FrigoJohnson05">{{Cytuj pismo |nazwisko=Frigo |imię=M. |nazwisko2=Johnson |imię2=S.G. |tytuł=The Design and Implementation of FFTW3 |url=http://fftw.org/fftw-paper-ieee.pdf |czasopismo=Proceedings of the IEEE |wolumin=93 |wydanie=2 |strony=216–231 |rok=2005 |doi=10.1109/JPROC.2004.840301}}</ref><ref>{{cytuj stronę |url=http://www.fftw.org/|tytuł=FFTW |nazwisko=Frigo |imię=Matteo |nazwisko2=Johnson |imię2=Steven G.}} A free ([[GNU General Public License|GPL]]) C library for computing discrete Fourier transforms in one or more dimensions, of arbitrary size, using the Cooley-Tukey algorithm.</ref>.

Typowa strategia dla algorytmów działających w miejscu bez dodatkowej pamięci i bez oddzielnych etapów odwracających cyfry pociąga za sobą transpozycję małych macierzy (które zamieniają poszczególne pary cyfr) na etapach pośrednich, które mogą być łączone z radix motylkowym w celu zredukowania liczby przejść nad danymi<ref name="FrigoJohnson05" /><ref>{{Cytuj pismo |nazwisko=Johnson |imię=H.W. |nazwisko2=Burrus |imię2=C.S. |tytuł=An in-place in-order radix-2 FFT |czasopismo=Proc. ICASSP |wolumin= |wydanie= |strony=28A.2.1–28A.2.4 |rok=1984}}</ref><ref>{{Cytuj pismo |nazwisko=Temperton |imię=C. |tytuł=Self-sorting in-place fast Fourier transform |czasopismo=SIAM J. Sci. Stat. Comput. |wolumin=12 |wydanie=4 |strony=808–823 |rok=1991 |doi=10.1137/0912043}}</ref><ref>{{Cytuj pismo |nazwisko=Qian |imię=Z. |nazwisko2=Lu |imię2=C. |nazwisko3=An |imię3=M. |nazwisko4=Tolimieri |imię4=R. |tytuł=Self-sorting in-place FFT algorithm with minimum working space |czasopismo=IEEE Trans. ASSP |wolumin=52 |wydanie=10 |strony=2835–2836 |rok=1994 |doi=10.1109/78.324749}}</ref><ref>{{Cytuj pismo |nazwisko=Hegland |imię=M. |tytuł=A self-sorting in-place fast Fourier transform algorithm suitable for vector and parallel processing |czasopismo=Numerische Mathematik |wolumin=68 |wydanie=4 |strony=507–547 |rok=1994 |doi=10.1007/s002110050074}}</ref>.

== Przypisy ==
{{Przypisy}}

== Bibliografia ==
* {{cytuj pismo |nazwisko=Cooley |imię=James W. |nazwisko2=Tukey |imię2=John W. |tytuł=An algorithm for the machine calculation of complex Fourier series |czasopismo=[[Mathematics of Computation|Math. Comput.]] |wolumin=19 |strony=297–301 |rok=1965 |doi=10.2307/2003354}}

== Linki zewnętrzne ==
* [http://www.librow.com/articles/article-10 prosty, pedagogiczny algorytm radix-2 Cooley-Tukeya FFT w C++.]
* [http://sourceforge.net/projects/kissfft/ KISSFFT] – prosta mixed-radix implementacja Cooleya-Tukeya w C (open source)

{{Kontrola autorytatywna}}

[[Kategoria:Algorytmy numeryczne]]