Algorytm Cooleya-Tukeya

Algorytm Cooleya-Tukeya – algorytm szybkiej transformacji Fouriera (FFT). Wyraża dyskretną transformację Fouriera (DFT) o dowolnej złożonej wielkości $N=N_{1}N_{2}$ w członach mniejszych DFT wielkości $N_{1}$ i $N_{2},$ rekurencyjnie, w celu ograniczenia czasu obliczeń do ${\text{O}}(N\log N),$ szczególnie w przypadku $N$ 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.

Algorytm Cooleya-Tukeya

Złożoność
Czasowa	$O(n\log n)$

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). 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 ZSRR. 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).

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 $N$ na dwie przeplatane DFT (stąd nazwa „radix-2”), o rozmiarze $N/2$ w każdym etapie rekurencyjnym.

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

$X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk},$

gdzie $k$ jest liczbą całkowitą z zakresu od $0$ do $N-1.$

Radix-2 DIT najpierw oblicza DFT dane wejściowe o indeksach parzystych $x_{2m}(x_{0},x_{2},\dots ,x_{N-2})$ i indeksach nieparzystych $x_{2m+1}(x_{1},x_{3},\dots ,x_{N-1}),$ 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 ${\text{O}}(N\log N).$ Ta uproszczona postać zakłada że $N$ jest potęgą dwójki; skoro liczba punktów próbki $N$ może być zwykle wybrana przez aplikację, to jest często nieważne ograniczenie.

Algorytm Radix-2 DIT przestawia DFT funkcji $x_{n}$ na dwie części: sumę dla indeksów parzystych $n={2m}$ i sumę dla indeksów nieparzystych $n={2m+1}{:}$

${\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}}$

Można uwzględniać wspólny mnożnik $e^{-{\frac {2\pi i}{N}}k}$ z drugiej sumy, jak przedstawiono w równaniu poniżej. Jest wtedy jasne że dwie sumy są DFT części o parzystych indeksach $x_{2m}$ i DFT części o nieparzystych indeksach $x_{2m+1}$ funkcji $x_{n}.$ Oznaczmy DFT danych wejściowych o parzystych indeksach $x_{2m}$ przez $E_{k}$ a DFT danych o nieparzystych indeksach $x_{2m+1}$ przez $O_{k}$ i otrzymujemy:

$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}.$

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

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

Ten wynik, wyrażając DFT o długości $N$ rekurencyjnie w członach dwóch DFT wielkości $N/2,$ 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ę $E_{k}$ i $O_{k}\exp(-2\pi ik/N),$ 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).

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 $N/2$ 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 $N$ jako dwie DFT rozmiaru $N/2$ jest czasami zwane lematem Danielsona-Lanczosa ponieważ ta tożsamość była zauważona przez tych dwóch autorów w 1942 roku (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 $N$ jest całkowitą potęgą dwójki 2, a $s=1$ jest rozstawem danych wejściowej tablicy $x.$ Oznaczenie $x+s$ oznacza tablicę zaczynająca się od $x_{s}.$

(Wyniki są w poprawnej kolejności w $X,$ 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ż $N=1$ do amortyzowania narzutu rekursji, czynniki fazowe $\exp[-2\pi ik/N]$ 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. 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

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 $N=N_{1}N_{2}$ jako:

1. Wykonaj $N_{1}$ DFT wielkości $N_{2}.$
2. Pomnóż przez zespolony pierwiastek z jedynki zwany czynnikami fazowymi.
3. Wykonaj $N_{2}$ DFT wielkości $N_{1}.$

Zazwyczaj $N_{1}$ lub $N_{2}$ jest niewielkim czynnikiem (niekoniecznie pierwszym), zwanym podstawą (radix) (który może się różnić pomiędzy etapami rekursji). Jeśli podstawą będzie $N_{1},$ wówczas mówimy o algorytmie decymacji w dziedzinie czasu (DIT). Gdy z kolei $N_{2}$ 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 ${\text{O}}(N_{1})$ dla pierwszych przypadków bazowych rekurencji (można także zastosować algorytm o złożoności $N\log N$ 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, chociaż ostatnie wariacje osiągnęły nawet niższą wartość. Na komputerach współczesnych, wydajność jest określana przez większą ilość 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). Inną cechą algorytmu Cooleya-Tukeya jest to, że wyraża rozmiar $N$ jednowymiarowego DFT jako $N_{1}$ przez $N_{2}$ dwuwymiarową DFT (plus czynniki fazowe), gdzie macierz wyjściowa jest 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 ${\sqrt {N}}$ 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 (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).

Ogólna faktoryzacja Cooleya-Tukeya przepisuje indeksy $k$ oraz $n$ w postaci $k=N_{2}k_{1}+k_{2}$ i $n=N_{1}n_{2}+n_{1},$ gdzie indeksy $k$ i $n$ zaczynają się od $0\dots N_{a}-1$ (dla $a$ 1 lub 2). Oznacza to, że reindeksuje ona wejście $(n)$ i wyjście $(k)$ jako $N_{1}$ przez $N_{2}$ 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 $nk,$ $N_{1}n_{2}N_{2}k_{1}$ krzyż członów zanika (jego wykładnikiem jest jedność), a pozostałe człony dają

${\begin{aligned}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{aligned}}$

gdzie każda wewnętrzna suma jest DFT rozmiaru $N_{2},$ każda zewnętrzna suma jest DFT rozmiaru $N_{1},$ a człony w nawiasie [...] są czynnikiem fazowym.

Arbitralna podstawa $r$ (jak również mieszane podstawy) może być stosowana tak jak wykazali zarówno Cooley z Tukeyem, jak i Gauss (który przedstawił przykłady kroków radix-3 i radix-6). Cooley i Tukey pierwotnie przyjmowali że obliczenia motylkowe podstawy wymagają ${\text{O}}(r^{2})$ pracy i stąd obliczyli, iż złożoność dla podstawy $r$ wynosi ${\text{O}}(r^{2}N/r\log _{r}N)={\text{O}}(N\log _{2}(N)r/\log _{2}r);$ z liczenia wartości $r/\log _{2}r$ dla całkowitych wartości $r$ od 2 do 12 odkryto, że optymalna podstawa wynosi 3 (najbliższa liczba całkowita do e, która minimalizuje $r/\log _{2}r$). Ta analiza była błędna, chociaż podstawa obliczeń motylkowych jest również DFT i może być wykonana przez algorytm FFT w ${\text{O}}(r\log r)$ operacjach, stąd podstawa $r$ faktycznie anuluje złożoność ${\text{O}}(r\log(r)N/r\log _{r}N),$ a optymalne $r$ jest określone przez bardziej skomplikowane rozważania. W praktyce dość duże $r$ (o wartości 32 lub 64) są potrzebne w celu skutecznego wykorzystania np. dużej liczby rejestrów procesora na nowoczesnych procesorach, a nawet nieograniczona podstawa $r={\sqrt {N}}$ również osiąga złożoność ${\text{O}}(N\log N)$ i wykazuje teoretyczne i praktyczne korzyści dla dużych $N,$ jak wspomniano powyżej.

Algorytmy zmiany kolejności danych, odwracania bitów i obliczeń w miejscu

Choć powyższe streszczenie 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 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 permutacją, gdzie dana o indeksie $n,$ zapisana binarnie cyframi $b_{4},$ $b_{3},$ $b_{2},$ $b_{1},$ $b_{0}$ (tj. 5 cyfr dla $N=32$ danych wejściowych), jest przeniesiona do indeksu o odwróconych cyfrach $b_{0},$ $b_{1},$ $b_{2},$ $b_{3},$ $b_{4}.$ 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 $E_{k}$ i $O_{k}$ 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 $b_{4}$ (dla $N=32$), podczas gdy dwie dane wejściowe $E_{k}$ i $O_{k}$ są przeplatane parzystymi i nieparzystymi elementami, co odpowiada najmłodszemu bitowi $b_{0}.$ 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. 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ń (choć odwrócenie bitów można wykonać w czasie ${\text{O}}(n),$ co było przedmiotem wielu badań). 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 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. Jeszcze większy potencjał zalet SIMD (więcej kolejnych dostępów) zostały zaproponowany dla algorytmu Pease, który również zmienia kolejność metodą „nie w miejscu” w każdej fazie, lecz wymaga ona odrębnego odwrócenia bitu/cyfry i ${\text{O}}(N\log N)$ pamięci. Można również bezpośrednio zastosować definicję faktoryzacji Cooleya-Tukeya z jawną rekurencją (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.

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.

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)