Algorytm FFT

e^2πilj/n= 1 n

n−1X

k=0

a_k

n−1X

j=0

e2πi(l−k)j/n. Jeśli k = l, to e2πi(l−k)j/n= e⁰= 1, zaś jeśli k 6= l, to liczby

e2πi(l−k)j/nsą pierwiastkami zespolonymi z 1; ich suma dla j∈ {0, . . . , n − 1} jest równa 0. Stąd wynika, że dl= a_l.

Interpolacja trygonometryczna i dyskretna transformata Fouriera występuje w wielu problemach związanych z analizą, transmisją i przetwarzaniem sygnałów (np. akustycznych lub obrazów), a także w rozwiązywaniu równań różniczkowych.

. . 361 . .

Algorytm FFT

Możemy zauważyć, że ciąg (b_j)_j∈Z, który jest transformatą ciągu (a_k)_k∈Z, składa się z wartości wielomianu stopnia n − 1

o współczynnikach a₀, . . . , a_n−1w punktach e^−2πij/n,

j = 0, . . . , n − 1. Dyskretną transformatę Fouriera można wyznaczyć za pomocą schematu Hornera; wyznaczenie pełnej transformaty kosztowałoby wtedy n²− nmnożeń i dodawań zespolonych. Okazuje się, że można to zadanie rozwiązać kosztem Θ(n(p₁+· · · + pr)) działań, gdzie p₁, . . . , prsą liczbami pierwszymi, takimi że

n = p₁· . . . · p^r. Odkrycia tego dokonali w 1952 r. Cooley i Tukey.

. . 362 .

Zauważmy, że ciąg okresowy (a_k)_k∈Z o okresie 1 jest ciągiem stałym i jest on identyczny ze swoją dyskretną transformatą Fouriera. Dalej, przypuśćmy, że liczba n jest podzielna przez p > 1. Oznaczmy w_j = e^−2πij/n. Wtedy wzór definiujący dyskretną transformatę Fouriera można przedstawić w postaci

b_j=

n/p−1

X

k=0

a_pkw^pk_j + w_j

n/p−1

X

k=0

a_pk+1w^pk_j +· · ·

+ w^p−1_j

n/p−1

X

k=0

a_pk+p−1w^pk_j .

Podzieliliśmy tu ciąg a₀, . . . , a_n−1na podciągi n/p-elementowe, wybierając do każdego z nich co p-ty element. Możemy dalej zauważyć, że sumy mnożone przez kolejne potęgi liczby w_j są wyrażeniami opisującymi transformaty tych podciągów, a dokładniej ich obustronnie nieskończonych rozszerzeń o okresie n/p.

Obliczenie dyskretnej transformaty Fouriera dla ciągu o okresie n może być zatem wykonane przez następujący algorytm rekurencyjny:

• Jeśli n = 1, to przyjmij b0 = a₀ (dla n = 1 przekształceniu poddajemy ciąg stały, którego obrazem jest ten sam ciąg).

• Jeśli n jest liczbą pierwszą, to zastosuj wzór podany jako definicja dyskretnej transformaty Fouriera i użyj schematu Hornera.

• Jeśli n > 1 jest podzielne przez liczbę pierwszą p < n, to podziel ciąg na p podciągów (zgodnie z opisem wyżej), oblicz

transformaty tych podciągów i „scal” je, stosując wzór podany wyżej i schemat Hornera.

Wzór opisujący transformatę odwrotną może być przekształcony podobnie; zamiast w_j= (cos^2πj_n , −sin^2πj_n )występuje w nim liczba w_j = (cos^2πj_n,sin^2πj_n). Możemy zatem użyć takiego samego algorytmu, zostawiając mnożenie wyniku działania procedury rekurencyjnej przez czynnik _n¹ na sam koniec. Koszt algorytmu w istotny sposób zależy od możliwości rozłożenia liczby n na czynniki.

Algorytm jest najbardziej efektywny, jeśli liczba n jest potęgą liczby 2 i często określenie FFT (od angielskiego Fast Fourier Transform) dotyczy takiego wariantu algorytmu. Zbadamy go dokładniej. Dla parzystej liczby n transformatę otrzymamy przez „scalenie”

transformat dwóch podciągów, złożonych odpowiednio z elementów parzystych i nieparzystych ciągu danego. Oznaczymy je symbolami (p_j)_j∈Z i (q_j)_j∈Z.

. . 365 . .

Transformaty te są ciągami obustronnie nieskończonymi, o okresie n/2, reprezentowanymi przez podciągi p₀, . . . , p_n/2−1

i q₀, . . . , q_n/2−1. Możemy napisać

b_j =

n/2−1

X

k=0

a_2kw^2k_j + w_j

n/2−1

X

k=0

a_2k+1w^2k_j = p_j+ w_jq_j. Podstawiając j + n/2 w miejsce j, i biorąc pod uwagę, że w_j+n/2= e−2πi(j+n/2)/n= e^−2πij/ne^{−2πin/(2n)} = −w_j oraz w^2k_j+n/2= w^2k_j , dostajemy

b_j+n/2=

n/2−1X

k=0

a_2kw^2k_j+n/2+ w_j+n/2

n/2−1X

k=0

a_2k+1w^2k_j+n/2

= p_j− w_jq_j.

. . 366 .

Implementacja algorytmu FFT w postaci procedury rekurencyjnej:

void rFFT ( int n, complex a[] ) {

complex *p, *q, u, w, t; int j;

if ( n > 1 ) {

p = malloc ( n*sizeof(complex) );

q = &p[n/2];

for ( j = 0; j < n/2; j++ ) { p[j] = a[2*j];

q[j] = a[2*j+1];

}

rFFT ( n/2, p );

rFFT ( n/2, q );

u = 1;

w = e^−2πi/n;

for ( j = 0; j < n/2; j++ ) { t = u*q[j];

a[j] = p[j] + t;

a[j+n/2] = p[j] - t;

u = u*w;

}

free ( p );

}

} /*rFFT*/

Widzimy, że choć procedura umieszcza transformatę w tej samej tablicy, w której są początkowo dane liczby a₀, . . . , a_n−1, potrzebuje ona sporo pamięci dodatkowej (w rzeczywistości potrzeba

dodatkowych tablic o sumarycznej długości 2n).

Można zaprojektować taką implementację, która wszystkie obliczenia wykonuje „w miejscu”, tj. która oprócz tablicy z danym ciągiem, który należy zastąpić przez jego transformatę, potrzebuje tylko niewielkiej ustalonej liczby zmiennych prostych. Aby otrzymać taką procedurę, nierekurencyjną i dodatkowo oszczędzającą pewne działania,

przyjrzymy się „przepływowi danych”, to znaczy zbadamy, od których współczynników zależą transformaty obliczane „po drodze”.

Dla n = 8 „przepływ danych” jest przedstawiony na rysunku (najlepiej go oglądać od prawej do lewej strony).

. . 369 . .

Krawędzie łączą dane z wynikami, tj. każda liczba (z wyjątkiem danych) jest obliczana na podstawie liczb znajdujących się

w kolumnie na lewo od niej, połączonych z nią kreskami. Widzimy, że w każdym przypadku obliczenie polega na zastąpieniu pary liczb przez inną parę, obliczoną tylko na jej podstawie (i dana para liczb nie jest do niczego innego potrzebna). Jeśli zatem ustawimy dane wejściowe w odpowiedniej kolejności, to można całe obliczenie wykonać bez potrzeby rezerwowania dodatkowej tablicy.

. . 370 .

Ostatnia transformata powstaje z transformat podciągów „parzystego”

i „nieparzystego”. Każda z tych dwóch transformat jest obliczana na podstawie transformat „parzystego” i „nieparzystego” podciągu odpowiedniego podciągu itd.; zatem ogólna reguła porządkowania danych wejściowych polega na ustawieniu ich w kolejności odwróconych bitów. Jeśli indeks j danego współczynnika a_j przedstawimy w układzie dwójkowym, przy użyciu l = log₂ncyfr dwójkowych (bitów), to indeks miejsca w tablicy, na którym ma się on znaleźć, otrzymamy wypisując te bity w odwrotnej kolejności.

Procedura FFT, która realizuje to obliczenie, ma postać void FFT ( int n, complex a[] )

{

complex t, u, w;

int i, j, k, l, m, p;

l = log₂n; m = n/2;

/* przestawianie danych w tablicy */

for ( i = 1, j = m; i < n-1; i++ ) {

/* obliczanie transformaty */

for ( k = 1; k <= l; k++ ) { m = 2^k; p = m/2;

u = 1; w = e^−πi/p;

for ( j = 0; j < p; j++ ) { i = j;

do {

t = a[i+p]*u;

a[i+p] = a[i]-t; a[i] = a[i]+t;

i += m;

} while ( i <= n );

u *= w;

} }

} /*FFT*/

. . 373 . .

Algorytm ten opublikowali w 1965 r. Cooley, Lewis i Welch. Pierwsza pętla, for ( i = ... ) ..., dokonuje przestawienia elementów w tablicy zgodnie z kolejnością odwróconych bitów. Kolejne przebiegi drugiej pętli, for ( k = ... ) ..., mają na celu obliczenie n/2 transformat podciągów o okresie 2, n/4 transformat podciągów o okresie 4, itd. Liczba e^−2πi/n(wartość zmiennej w) i jej potęgi, czyli liczby w_j (kolejne wartości zmiennej u) są obliczane tylko raz dla wszystkich transformat podciągów o tym samym okresie. W każdym przebiegu pętli for ( j = ... ) ... obliczane są pary

współczynników o numerach j oraz j + p we wszystkich transformatach podciągów o okresie m = 2p, ponieważ pętla najbardziej wewnętrzna (do...while) przebiega przez wszystkie te transformaty.

. . 374 .

Można udowodnić, że algorytm FFT, także w wersji ogólnej (dla dowolnego n), jest numerycznie stabilny, tj. istnieje stała K (zależna od n), taka że współczynniki ˜b_j obliczone przy użyciu arytmetyki zmiennopozycyjnej przybliżają dokładne współczynniki b_j dyskretnej transformaty Fouriera z błędem spełniającym nierówność

maxj |˜b_j− b_j| 6 Kνmax

j |b_j|,

gdzie ν = 2^−t. Dla liczby n będącej potęgą 2 można przyjąć K =^√

2log₂n + (log₂n − 1)(3 + 2ε)^√ n,

gdzie ε jest oszacowaniem błędu bezwzględnego obliczonych kosinusów i sinusów, tj. części rzeczywistych i urojonych liczb w_j.