Algorytm Cooleya-Tukeya

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

Algorytm Cooleya-Tukeya nazwany imieniem J.W. Cooleya i Johna Tukeya, jest najbardziej rozpowszechnionym algorytmem 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.

Ponieważ algorytm Cooleya-Tukeya rozbija DFT na mniejsze DFT, może być połączony arbitralnie z dowolnym innym algorytmuem 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 | edytuj kod]

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 i w neo-łacinie)[1][2]. Jakkolwiek Gauss nie analizuje asymptotycznego czasu obliczeń. Różne ograniczone formy zostały również ponownie odkryte kilka razy w ciągu XIX i 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ł w którym ponownie wynaleźli algorytm i opisali jak wykonać go wygodnie na komputerze[3].

Tukey podobno wpadł na pomysł podczas spotkania amerykańskiej prezydenckiej komisji doradczej omawiania sposobów 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, który implementował go do innego (niesklasyfikowanego) problemu: analiza 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 co teraz nazywamy 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 różnym 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łozónych rozmiarów przez algorytm Cooleya–Tukeya)[6].

Przypadek Radix-2 DIT[edytuj | edytuj kod]

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, \ldots, x_{N-2}) i indeksach nieparzystych x_{2m+1} \ (x_1, x_3, \ldots, 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 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 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:


\begin{matrix} X_k= \underbrace{\sum \limits_{m=0}^{N/2-1} x_{2m}  e^{-\frac{2\pi i}{N/2} mk}}_{\mathrm{DFT\;of\;even-indexed\;part\;of\;} 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\;of\;odd-indexed\;part\;of\;} x_m} = E_k + e^{-\frac{2\pi i}{N}k} O_k.
\end{matrix}

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 \leq k z DFT długości N/2 są identyczne jak dane wyjściowe dla 0\leq k . To znaczy że E_{k + N/2} = E_k i O_{k + N/2} = O_k. Czynnik fazowy \exp[-2\pi i k/ N] jest posłuszny relacji: \exp[-2\pi i (k + N/2)/ N] = e^{-\pi i} \exp[-2\pi i k/ N] = -\exp[-2\pi i k/ N], przerzucania znaku O_{k + N/2} członu. Tak więc, całą DFT można obliczyć w następujący sposób:


\begin{matrix} X_k & =
& \left\{
\begin{matrix}
E_k + e^{-\frac{2\pi i}{N}k} O_k & \mbox{if } k < N/2 \\ \\
E_{k-N/2} - e^{-\frac{2\pi i}{N} (k-N/2)} O_{k-N/2} & \mbox{if }
k \geq N/2. \end{matrix} \right. \end{matrix}

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 i k/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 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] (wpływ pracy Runge's 1903[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 sekund na operację zmiennoprzecinkową, około 20% z nich to mnożenia).

Pseudokod[edytuj | edytuj kod]

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

Tutaj, ditfft2 (x ,N ,1), oblicza X =DFT(x ) nie na miejscu przez radix-2 DIT FFT, 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 i nie wymaga się dalszych permutacji odwracających bit; często wspominana konieczność oddzielnej fazy odwracania bitów powstaje jedynie dla pewnych algorytmów działających w miejscu jak opisane poniżej.)

Wysokowydajne implementacje FFT dokonują wielu zmian w realizacji takiego algorytmu w porównaniu z tym prostym pseudokodem. Na przykład, można użyć większego przypadku bazowego niż N = 1 do amortyzowania narzutu rekursji, Czynniki fazowe \exp[-2\pi i k/ N] można wcześniej obliczyć, a większe podstawy są często stosowane ze względu na pamięć cache; te i inne optymalizacje 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 lepszym umiejscowieniem pamięci[8][9]. ) Kilka z tych koncepcji jest opisane bardziej szczegółowo poniżej.

Ogólne faktoryzacje[edytuj | edytuj kod]

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.

Bardziej ogólnie, algorytmy Cooleya–Tukeya rekurencyjnie wyrażają 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, albo N 1 lub N 2 jest niewielkim czynnikiem (niekoniecznie pierwszym ), zwany podstawą (radix) (który może się różnić pomiędzy etapami rekursji). Jeśli N 1 jest podstawą, nazywa się algorytmem decymacji w dziedzinie czasu (DIT), podczas gdy N 2 jest podstawa, to decymacji w dziedzinie częstotliwości (DIF, zwany także algorytmem Sande-Tukeya). Wersja przedstawiona powyżej miała algorytm 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 znana jako motyl 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 ażeby osiągnąć najmniejszą znaną ilość operacji dla rozmiarów potęgi dwójki[10], chociaż ostatnie wariacje osiągnęły nawet niższy koszt[11][12]. (Na komputerach współczesnych, wydajność jest określana przez większą ilość pamięci podręcznej i kwestie potokowości procesora niż ściśle liczonej pracy, dobrze zoptymalizowane implementacje FFT często używają większe podstawy i/lub zakodowane na stałe znaczącej wielkości[13]. ) Innym sposobem patrzenia na algorytm 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. Wynikiem tych wszystkich transpozycji, dla algorytmu radix-2, odpowiada odwrócenie bitowe indeksów wejściowych (DIF) lub wyjściowych (DIT). Jeżeli zamiast użycia małej podstawy używa się podstawy około √ N i jawnej transpozycji matrycy wejścia / wyjścia , jest nazywany czteroetapowym algorytmem (lub sześcio krokowym, w zależności od liczby transpozycji), wstępnie zaproponowano poprawić położenie pamięci[14][15] np. do optymalizacji pamięci podręcznej lub algorytmu z pamięcią zewnętrzną działania, a później okazuje się optymalna nie baczącego na rozmiar cache algorytmu[16].

Ogólna faktoryzacja Cooleya-Tukeya przepisuje indeksy k oraz n k = N_2 k_1 + k_2 i n = N_1 n_2 + n_1 , odpowiednio, gdzie indeksy k i n zaczynają się od 0..N a-1 (dla a 1 lub 2). To jest, ona reindeksuje 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ą jak wspomniano powyżej. Kiedy to przeindeksowanie jest podstawione do formuły 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ą

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 }

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) mogą być stosowane, jak było pokazane przez Cooleya i Tukeya[3] jak również przez Gaussa (który dał 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 liczona złożoność dla podstawy r jest 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 optymalna podstawa znaleziona jest 3 (the najbliższa całkowita do e , która minimizuje r /log2r )[3][17]. Ta analiza była błędna, chociaż: podstawa obliczeń motylkowych jest również DFT i może być wykonany 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 ), i optymalne r jest określone przez bardziej skomplikowane rozważania. W praktyce dość duże r (32 lub 64) są ważne, aby skutecznie wykorzystać np. na dużą liczbę rejestrów procesora s na nowoczesnych procesorach[13], , a nawet nieograniczona podstawa r = √ N osiąga również złożoność O (N log N) i ma 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 | edytuj kod]

Choć streszczenie faktoryzacji DFT Cooley-Tukeya powyżej ma zastosowanie w jakiejś 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 z jego danych wyjściowych przy użyciu tylko O(1) dodatkowej pamięci.

Najbardziej znana techniką przeorganizowania obejmuje jawne odwrócenie bitu dla algorytmu radix-2 w miejscu. 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 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 iść do pierwszej i drugiej połowy tablicy wyjściowej, odpowiednio do najstarszego bitu 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 dostać dane wyjściowe w porządku wejściowym. (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.) Odpowiednio, jeśli wykonasz wszystkie czynności w odwrotnej kolejności, uzyskać algorytm radix-2 DIF z bitowym odwróceniem w postprocessingu (lub preprocessing, odpowiednio). Alternatywnie, niektóre zastosowania (takie jak splot) działa równie dobrze na danych z bitami odwróconymi, więc można wykonać transformacje w przód, przetwarzanie, a następnie przekształcić odwrotne wszystko bez odwracania bitów do wytwarzania końcowych wyników w naturalnym porządku.

Wielu użytkowników FFT, jednak preferuje dane wyjściowe w naturalnym porządku i oddzielna jawna faza odwrócenia bitów może mieć istotny wpływ na czas obliczeń[13], choć odwrócenie bitów można zrobić w czasie O(n) i było to przedmiotem wielu badań[18][19][20]. Ponadto, podczas gdy permutacja jest odwróceniem bitów w przypadku radix-2, jest to bardziej ogólnie w arbitralnym (mixed-base) odwrócenie cyfr dla przypadku mixed-radix i permutacji algorytmy stają się bardziej skomplikowane do zaimplementowania. Co więcej, jest pożądane na wielu architekturach sprzętowych uporządkować pośrednie fazy algorytmu FFT tak, by one operowały na kolejnych (lub co najmniej bardziej zlokalizowanych) elementach danych. W tym celu, wiele alternatywnych metod implementacji opracowano 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 jest poza miejscem: tablica wyjściowa różni się od tablicy wejściowej lub, równoważnie, równa wielkości pomocnicza tablica jest dostępna. Alorytm Stockhama automatycznego sortowania[21] wykonuje każdy etap FFT nie na miejscu, zwykle pisze się tam i z powrotem pomiędzy dwoma tablicami, transponując jedną "cyfrę" indeksów w każdym etapie, i było to szczególnie popularne w architekturze SIMD[22]. Nawet większy potencjał zalet SIMD (więcej kolejnych dostępów) zostały zaproponowane dla algorytmu Pease,[23] , który również zmienia kolejność nie w miejscu w każdej fazie, ale metoda ta wymaga 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, która produkuje naturalny porządek danych wyjściowych nie-w-miejscu bez osobnego kroku permutacji (jak w pseudokodzie powyżej) i można argumentować, że będą niezależne od wielkości cache korzyści z lokalizacji na systemach z hierarchiczną pamięcią[9][13][24].

Typowa strategia dla algorytmów działających w miejscu bez dodatkowej pamięci i bez oddzielnych przebiegó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, aby zmniejszyć liczbę przejść nad danymi[13][25][26][27][28].

Bibliografia[edytuj | edytuj kod]

  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. 2,0 2,1 2,2 2,3 2,4 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. 3,0 3,1 3,2 3,3 3,4 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. 8,0 8,1 8,2 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. 9,0 9,1 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. 10,0 10,1 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. 13,0 13,1 13,2 13,3 13,4 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. 14,0 14,1 Gentleman W. M., and G. Sande, "Fast Fourier transforms—for fun and profit," Proc. AFIPS 29 , 563–578 (1966).
  15. 15,0 15,1 Bailey, David H., "FFTs in external or hierarchical memory," J. Supercomputing 4 (1), 23–35 (1990)
  16. 16,0 16,1 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. Alan H. Karp. Bit reversal on uniprocessors. „SIAM Review”. 38 (1), s. 1–26, 1996. doi:10.1137/1038001. 
  19. Larry Carter, Kang Su Gatlin. Towards an optimal bit-reversal permutation program. , s. 544–553, 1998. doi:10.1109/SFCS.1998.743505. 
  20. A new superfast bit reversal algorithm. „Intl. J. Adaptive Control and Signal Processing”. 16 (10), s. 703–707, 2002. doi:10.1002/acs.718. 
  21. T. G. Stockham. High speed convolution and correlation. „Spring Joint Computer Conference, Proc. AFIPS”. 28, s. 229–233, 1966. 
  22. Vectorizing the FFTs. W: P. N. Swarztrauber: Parallel Computations. G. Rodrigue. New York: Academic Press, 1982, s. 51–83. ISBN 0-12-592101-2.
  23. 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. 
  24. 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
  25. H. W. Johnson, C.S. Burrus. An in-place in-order radix-2 FFT. „Proc. ICASSP”, s. 28A.2.1–28A.2.4, 1984. 
  26. C. Temperton. Self-sorting in-place fast Fourier transform. „SIAM J. Sci. Stat. Comput.”. 12 (4), s. 808–823, 1991. doi:10.1137/0912043. 
  27. 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. 
  28. 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 | edytuj kod]