Algorytm Cooleya-Tukeya

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj
Algorytm Cooleya-Tukeya
Złożoność
Czasowa О(n log2(n))

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

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[edytuj]

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)[1][2]. 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[2]. 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[3].

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.[4][5] 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[2]. Ich artykuł cytuje jako inspirację jedynie pracę I. J. Good tym, obecnie nazywaną "algorytmem FFT czynników pierwszych (PFA)"[3]; 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)[6].

Przypadek Radix-2 DIT[edytuj]

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:

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 O (N log N). Ta uproszczona postać zakłada że N jest potęgą dwójki; skoro ilość 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 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ść 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 z DFT długości N/2 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 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ę 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).

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 w 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[7] (na co wpływ miała praca C. Runge'a z 1903 roku.[2]). 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)[3]. 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[edytuj]

W pseudokodzie powyższy proces można zapisać:[8]

X 0,...,N −1ditfft2 (x , N , s ):             DFT of (x 0, x s , x 2s , ..., x (N -1)s ):
    if N  = 1 then
        X 0x 0                                      ' przypadek bazowy trywialnego rozmiaru DFT  = 1
    else
        X 0,...,N /2−1ditfft2 (x , N /2, 2s )             DFT (x 0, x 2s , x 4s , ...)
        X N /2,...,N −1ditfft2 (x +s, N /2, 2s )           DFT (x s , x s +2s , x s +4s , ...)
        for k  = 0 to N /2−1                           łączenie DFT dwóch połówek do pełnego DFT:
            t ← X k 
            X k  ← t + exp(−2πi  k /N ) X k +N /2
            X k +N /2 ← t − exp(−2πi  k /N ) X k +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 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[8]. 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[8][9]. Poniżej opisano szczegółowo kilka z wymienionych koncepcji.

Ogólne faktoryzacje[edytuj]

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 1N 2 jako:[10]

  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óczas 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óznymi (zwykle małymi) czynnikami oprócz dwójki, zwykle (ale nie zawsze) stosujące algorytm złożoności O(N 2) 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 ilości operacji dla rozmiarów potęgi dwójki[10], chociaż ostatnie wariacje osiągnęły nawet niższą wartość[11][12]. 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 ilości 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[13]. ) 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 √ N i jawnej transpozycji matrycy wejścia / wyjścia , mówimy o czteroetapowym algorytmie (lub sześcio krokowym, w zależności od liczby transpozycji). Taki algorytm wstępnie został zaproponowany w celu poprawy alokacji pamięci[14][15](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.)[16].

Ogólna faktoryzacja Cooleya-Tukeya przepisuje indeksy k oraz n w postaci i , gdzie indeksy k i n zaczynają się od 0..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 , krzyż członów zanika (jego wykładnikiem jest jedność), a pozostałe człony dają

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 Tukey'em[3], jak i Gauss (który przedstawił przykłady kroków radix-3 i radix-6)[2]. Cooley i Tukey pierwotnie przyjmowali że obliczenia motylkowe podstawy wymagają O(r 2) pracy i stąd obliczyli, iż złożoność dla podstawy r wynosi O(r 2 N /r  logr N ) = O(N  log2(Nr /log2r ); z liczenia wartości r /log2r 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 /log2r )[3][17]. Ta analiza była błędna, chociaż podstawa obliczeń motylkowych jest również DFT i może być wykonana przez algorytm FFT w O(r log r ) operacjach, stąd podstawa r faktycznie anuluje złożoność O(r  log(rN /r  logr 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[13], a nawet nieograniczona podstawa r = √ N również osiąga złożoność O (N log N) i wykazuje teoretyczne i praktyczne korzyści dla dużych N jak wspomniano powyżej[14][15][16].

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

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ęcio 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 4b 3b 2b 1b 0 (t.j. 5 cyfr dla N =32 danych wejściowych), jest przeniesiona do indeksu o odwróconych cyfrach b 0b 1b 2b 3b 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 i 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 i 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.[18] 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ń[13] (choć odwrócenie bitów można wykonać w czasie O(n) co było przedmiotem wielu badań)[19][20][21]. 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[22] wykonuje każdy etap FFT nie na miejscu, zazwyczaj zapisując dane "w tył i w przód" pomiędzy dwoma tablicami, transponując jedną "cyfrę" indeksów w każdym etapie, co było szczególnie popularne w architekturze SIMD[23]. Jeszcze większy potencjał zalet SIMD (więcej kolejnych dostępów) zostały zaproponowany dla algorytmu Pease,[24], który również zmienia kolejność metodą "nie w miejscu" w każdej fazie, lecz wymaga ona odrębnego odwrócenia bitu/cyfry i 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[9][13][25].

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[13][26][27][28][29].

Bibliografia[edytuj]

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

Przypisy

  1. Gauss, Carl Friedrich, "Theoria interpolationis methodo nova tractata", Werke, Band 3, 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866)
  2. a b c d e Heideman, M. T., D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine, 1, (4), 14–21 (1984)
  3. a b c d e 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. 
  4. Cooley, Peter A. W. Lewis, Peter D. Welch. Historical notes on the fast Fourier transform. „IEEE Trans. on Audio and Electroacoustics”. 15 (2), s. 76–79, 1967. 
  5. Rockmore, 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 "[1]
  6. James W. Cooley, Peter A. W. Lewis, and Peter W. Welch, "Historical notes on the fast Fourier transform," Proc. IEEE , vol. 55 (no. 10), p. 1675–1677 (1967).
  7. Danielson, G. C., and C. Lanczos, "Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids," J. Franklin Inst. 233 , 365–380 and 435–452 (1942).
  8. a b c S. G. Johnson and M. Frigo, “Implementing FFTs in practice,” in Fast Fourier Transforms (C. S. Burrus, ed.), ch. 11, Rice University, Houston TX: Connexions, September 2008.
  9. a b Richard C. Singleton. On computing the fast Fourier transform. „Commun. of the ACM”. 10 (10), s. 647–654, 1967. DOI: 10.1145/363717.363771. 
  10. a b Duhamel, P., and M. Vetterli, "Fast Fourier transforms: a tutorial review and a state of the art," Signal Processing 19 , 259–299 (1990)
  11. Lundy, T., and J. Van Buskirk, "A new matrix approach to real FFTs and convolutions of length 2k ," Computing 80 , 23-45 (2007).
  12. Johnson, S. G., and M. Frigo, "A modified split-radix FFT with fewer arithmetic operations," IEEE Trans. Signal Processing 55 (1), 111–119 (2007).
  13. a b c d e M. Frigo, S. G. Johnson. The Design and Implementation of FFTW3. „Proceedings of the IEEE”. 93 (2), s. 216–231, 2005. DOI: 10.1109/JPROC.2004.840301. 
  14. a b Gentleman W. M., and G. Sande, "Fast Fourier transforms—for fun and profit," Proc. AFIPS 29 , 563–578 (1966).
  15. a b Bailey, David H., "FFTs in external or hierarchical memory," J. Supercomputing 4 (1), 23–35 (1990)
  16. a b M. Frigo, C.E. Leiserson, H. Prokop, and S. Ramachandran. Cache-oblivious algorithms. In Proceedings of the 40th IEEE Symposium on Foundations of Computer Science (FOCS 99), p.285-297. 1999. Extended abstract at IEEE, at Citeseer.
  17. Cooley, J. W., P. Lewis and P. Welch, "The Fast Fourier Transform and its Applications", IEEE Trans on Education 12 , 1, 28-34 (1969)
  18. Jeżeli każda podtransformacja rozmiaru N /2 operuje na przyległych danych, dane wejściowe DIT są wstępnie przetworzone przez odwrócenie bitów.
  19. Alan H. Karp. Bit reversal on uniprocessors. „SIAM Review”. 38 (1), s. 1–26, 1996. DOI: 10.1137/1038001. 
  20. Larry Carter, Kang Su Gatlin. Towards an optimal bit-reversal permutation program. , s. 544–553, 1998. DOI: 10.1109/SFCS.1998.743505. 
  21. A new superfast bit reversal algorithm. „Intl. J. Adaptive Control and Signal Processing”. 16 (10), s. 703–707, 2002. DOI: 10.1002/acs.718. 
  22. T. G. Stockham. High speed convolution and correlation. „Spring Joint Computer Conference, Proc. AFIPS”. 28, s. 229–233, 1966. 
  23. Vectorizing the FFTs. W: P. N. Swarztrauber: Parallel Computations. G. Rodrigue. New York: Academic Press, 1982, s. 51–83. ISBN 0-12-592101-2.
  24. M. C. Pease. An adaptation of the fast Fourier transform for parallel processing. „J. ACM”. 15 (2), s. 252–264, 1968. DOI: 10.1145/321450.321457. 
  25. Matteo Frigo, Steven G. Johnson: FFTW. A free (GPL) C library for computing discrete Fourier transforms in one or more dimensions, of arbitrary size, using the Cooley–Tukey algorithm
  26. H. W. Johnson, C.S. Burrus. An in-place in-order radix-2 FFT. „Proc. ICASSP”, s. 28A.2.1–28A.2.4, 1984. 
  27. C. Temperton. Self-sorting in-place fast Fourier transform. „SIAM J. Sci. Stat. Comput.”. 12 (4), s. 808–823, 1991. DOI: 10.1137/0912043. 
  28. Qian, C. Lu, M. An, R. Tolimieri. Self-sorting in-place FFT algorithm with minimum working space. „IEEE Trans. ASSP”. 52 (10), s. 2835–2836, 1994. DOI: 10.1109/78.324749. 
  29. M. Hegland. A self-sorting in-place fast Fourier transform algorithm suitable for vector and parallel processing. „Numerische Mathematik”. 68 (4), s. 507–547, 1994. DOI: 10.1007/s002110050074. 

Linki zewnętrzne[edytuj]