2 . Szybka transformata Fouriera
Wyznaczenie ciągu (Y0, Y1, . . . , YN −1) przy użyciu DFT wymaga wykonania:
– mnożenia zespolonego – (N − 1)2 razy, – dodawania zespolonego – N (N − 1) razy,
przy założeniu, że wartości ωNj są już dane. Najczęściej przyjmowane war- tości N są rzędu 1000, co daje około miliona operacji każdego rodzaju.
Naturalnym stało się więc poszukiwanie sposobów wyznaczania DFT, które pozwoliłyby na obniżenie kosztu tej metody. Algorytm taki został opisany w 1965 roku przez dwóch matematyków amerykańskich J. W. Cooleya i J.
W. Tukeya i nosi nazwę szybkiej transformaty Fouriera (FFT). Algorytm ten wykorzystuje specjalną postać macierzy przekształcenia FN, której ele- mentami są pierwiastki z jedynki. Koszt FFT jest rzędu N log N .
Niech N = 2m, m ∈ N.
Yn = 1 2m
2m−1 X k=0
ykωN−nk = 1 2m
m−1 X k=0
y2kωN−2nk+ y2k+1ωN−n(2k+1)
=
=1 2
1 m
m−1 X k=0
y2kωN−2nk + ωN−n 1 m
m−1 X k=0
y2k+1ωN−2nk
, ale
ωN−2nk = ω2m−2nk = e2iπ2m1 (−2nk) =
e2iπm1−nk = ωm−nk, czyli
Yn = 1 2
1 m
m−1 X k=0
y2kωm−nk + ω2m−n 1 m
m−1 X k=0
y2k+1ωm−nk
= 1 2
Pn + ω2m−nIn. Zauważmy, że
Pm+n = 1 m
m−1 X k=0
y2kωm−(m+n)k = 1 m
m−1 X k=0
y2kωm−nk = Pn
i podobnie
In+m = In oraz
ω2m−(m+n) = e−2iπ2m1 (m+n) = ω−n2m· e−iπ = −ω−n2m. Równości te prowadzą do algorytmu:
Krok 1. Obliczyć Pk i ω−kN Ik.
Krok 2. Utworzyć Yk = 12(Pk + ωN−kIk).
Krok 3. Wyznaczyć Yk+m = 12(Pk − ωN−kIk).
Kroki te należy wykonać dla k = 0, 1, . . . , m − 1.
Pk i Ik są dodatkowo dwoma niezależnymi DFT rzędu m = N2: FN
2 :(y0, y2, . . . , y2m−2) → (P0, P1, . . . , Pm−1), FN
2 :(y1, y3, . . . , y2m−1) → (I0, I1, . . . , Im−1).
Przy założeniu, że m jest parzyste, można więc algorytm powtórzyć. Jeśli N = 2p, to iterujemy ten proces tak długo, aż dojdziemy do DFT rzędu 2.
Ma ona postać:
Y0 =y0 + y1 2 , Y1 =y0 − y1
2 . Przykład: N = 8.
1. (y0, y1, . . . , y7) dzielimy na dwa ciągi długości 4:
(y0, y2, y4, y6), (y1, y3, y5, y7).
2. Ciągi te dzielimy na ciągi długości 2:
(y0, y4), (y2, y6), (y1, y5), (y3, y7).
3. Obliczamy dla każdego z tych ciągów P0 i I0 oraz wyznaczamy (Y0, Y1).
4. Przy użyciu formuł
Yk =1
2(Pk + ωN−kIk) Yk+m =1
2(Pk − ωN−kIk).
przechodzimy od wektora długości m do wektora długości 2m.
Dla N = 2p ogólny koszt tego algorytmu to 12N (log2N − 2) + 1 dodawań i N log2N mnożeń. Dla N = 1024 daje to nam koszt 250 razy mniejszy od kosztu DFT.
Program:
Procedure FFT(n,w,y,Y);
begin
if n=1 then Y[0]:=y[0] else begin
m:=n div 2;
for k:=0 to m-1 do begin
b[k]:=y[2*k];
c[k]:=y[2*k+1]
end;
w2:=w*w;
FFT(m,w2,b,B);
FFT(m,w2,c,C);
wk:=1;
for k:=0 to m-1 do begin
X:=B[k]; T:=wk*C[k];
Y[k]:=(X+T)/2;
Y[k+m]:=(X-T)/2;
wk:=wk*w end
end end.
Zastosowania FFT do obliczeń numerycznych Przykład 1.: Splot cykliczny.
1. Ciągi o wartościach zespolonych.
(xn)n∈Z, (hq)q∈Z - ciągi o wartościach zespolonych, okresowe o okresie N . Splot cykliczny tych ciągów zdefiniowany jest wzorem
yn =
N −1 X q=0
hqxn−q =
N −1 X q=0
hn−qxq.
Jest to ciąg okresowy o okresie N . Przy ustalonym (hq) przekształcenie zadane powyższym wzorem jest liniowym przekształceniem CN w siebie:
X → Y = HX,
X = (x0, x1, . . . , xN −1)T, Y = (y0, y1, . . . , yN −1)T,
H =
h0 hN −1 hN −2 . . . h1 h1 h0 hN −1 . . . h2
h2 h1 h0 . . . h0 ...
hN −1 hN −2 hN −3 . . . h0
.
Macierz tą nazywamy macierzą cykliczną.
Wyznaczenie splotu cyklicznego bezpośrednio z definicji wymaga wyko- nania:
mnożenia zespolonego - N2 razy,
dodawania zespolonego - N (N − 1) razy.
Możemy jednak rozważać DFT (Xk), (Hk) i (Yk) odpowiednio ciągów (xn), (hn) i (yn). Równanie zadające splot będzie miało wówczas postać
Yk = N HkXk, k = 0, 1, 2, . . . , N − 1.
Jeśli założymy, że N = 2p, to wykonujemy teraz następujące kroki:
Kolejne kroki Koszt
1. Wyznaczamy transformaty [N (p − 2); 2N p]
FN : (xn) → (Xk) FN : (hn) → (Hk)
2. Obliczamy iloczyn [N ; 0]
Yk = N HkXk, k = 0, 1, . . . , N − 1
3. Wyznaczamy transformatę odwrotną hN2(p − 2); N pi. FN−1 : (Yk) → (yn)
Koszt całkowity to
mnożenie zespolone – N2(3 log2N − 4) razy, dodawanie zespolone – 3N log2N razy.
Dla N = 64 daje to koszt [448; 1152], zamiast [4096, 4032].
2. Ciągi o wartościach rzeczywistych.
W tym przypadku krok 1. może być zastąpiony wyznaczeniem poje- dynczej transformaty zmiennych zespolonych zamiast dwóch transformat zmiennych rzeczywistych. W kroku 3. oblicza się transformatę odwrotną rzędu N ciągu, który spełnia warunek Y−k = Yk, możemy zastąpić to obli- czaniem transformaty rzędu N2.
Całkowity koszt wynosił będzie hN4(3 log2N − 5);N2 (3 log2N − 1)i opera- cji zespolonych, czyli hN (3 log2N − 5);N2 (9 log2N − 7)i operacji rzeczywi- stych. Przy N = 64 zmniejsza to liczbę wykonywanych operacji z [16 384;
16 256] do [832, 1 504].
Przykład 2.: Splot niecykliczny. Niech (xn)n∈Z, (hn)n∈Z będą dwoma sy- gnałami nieokresowymi o zwartych nośnikach. W szczególności niech
xn = 0 jeśli n < 0 lub n M ,
hn = 0 jeśli n < 0 lub n Q (Q M ).
Chcemy wyznaczyć
yn =
Q−1 X q=0
hqxn−q.
yn jest równe 0, jeśli n < 0 lub n M + Q − 1. Niech N będzie naj- mniejszą potęga liczby 2 taką, że N M + Q − 1. Tworząc z oryginalnego ciągu ciąg okresowy o okresie N , wracamy do problemu wyznaczenia splotu cyklicznego z użyciem FFT.
Koszt takiego postępowania może się jednak okazać wyższy niż koszt bezpośredniego rachunku w przypadkach, gdy długości tych dwóch sygnałów różnią się znacznie, np. gdy ciąg (xn) jest praktycznie „nieskończony”, a nośnik filtru (hn) relatywnie mały. Są jednak metody pozwalające obniżyć ten koszt, rozważa się w nich ciąg (xn) podzielony na mniejsze części.
