1 . Dyskretna transformata Fouriera
Problem: f jest funkcją o okresie a. Znamy jej wartości w N punktach
równomiernie rozłożonych na odcinku [0, a):
f
k
a
N
!
= y
k
dla k = 0, 1, 2, . . . , N − 1.
Mówiąc inaczej, sygnał f jest próbkowany w regularnych odstępach czasu
o długości
a
N
. Na podstawie tych informacji chcemy aproksymować współ-
czynniki c
n
szeregu Fouriera funkcji f .
Będziemy również zakładać, że szereg Fouriera funkcji f jest punktowo
zbieżny do f i w punktach nieciągłości t zachodzi równość
f (t) =
f (t+) + f (t−)
2
.
Mając dane N punktów będziemy wyznaczać N współczynników Fo-
uriera. Wiedząc, że c
n
→ 0 przy n → ∞, wyznaczymy współczynniki dla
n = −
N
2
, . . . ,
N
2
− 1 (lub n = −
N −1
2
, . . . ,
N −1
2
, jeśli N jest nieparzyste). Te
założenia pozwalają wyznaczyć przybliżoną wartość całki
c
n
=
1
a
Z
a
0
f (t)e
−2iπn
t
a
dt.
Sposób 1. Stosujemy metodę trapezów:
c
0
n
=
1
a
N −1
X
k=0
f
k
a
N
!
e
−2iπn
k
N
(k + 1)
a
N
− k
a
N
!
=
1
N
N −1
X
k=0
y
k
e
−2iπn
k
N
lub inaczej
c
0
n
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
,
gdzie ω
N
= e
2iπ
1
N
.
Sposób 2. Możemy wyznaczyć współczynniki Fouriera c
N
n
wielomianu try-
gonometrycznego
P (t) =
N
2
−1
X
n=−
N
2
c
N
n
e
2iπn
t
a
,
1
który interpoluje funkcję f w punktach k
a
N
, k = 0, 1, 2, . . . , N − 1. Musimy
rozwiązać układ równań liniowych rzędu N:
N
2
−1
X
n=−
N
2
c
N
n
ω
nk
N
= y
k
,
k = 0, 1, 2, . . . , N − 1.
Przesuwamy dla wygody ujemne indeksy o N w prawo. Możemy to zrobić,
gdyż wyrażenie ω
k
N
jest N -okresowe (ze względu na zmienną k). Stąd
−1
X
n=−
N
2
c
N
n
ω
nk
N
=
N −1
X
p=
N
2
c
N
p−N
ω
k(p−N )
N
=
N −1
X
p=
N
2
c
N
p−N
ω
kp
N
=
N −1
X
n=
N
2
c
N
n−N
ω
kn
N
.
Definiujemy:
Y
n
=
c
N
n
dla 0 ¬ n ¬
N
2
− 1
c
N
n−N
dla
N
2
¬ n ¬ N − 1
i wówczas układ równań ma postać
N −1
X
n=0
Y
n
ω
nk
N
= y
k
,
k = 0, 1, 2, . . . , N − 1.
Niech 0 ¬ p ¬ N − 1. Mnożymy obie strony równań przez ω
−kp
N
i doda-
jemy równania stronami:
N −1
X
k=0
y
k
ω
−kp
N
=
N −1
X
k=0
N −1
X
n=0
Y
n
ω
k(n−p)
N
=
N −1
X
n=0
Y
n
N −1
X
k=0
ω
k(n−p)
N
.
Zauważmy, że
N −1
X
k=0
ω
k(n−p)
N
=
0,
jeśli n 6= p
N,
jeśli n = p,
a więc
N −1
X
k=0
y
k
ω
−kp
N
= N Y
p
i stąd nieznane Y
n
są dane wzorami
Y
n
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
,
n = 0, 1, 2, . . . , N − 1.
2
Wniosek: Całkując metodą trapezów, otrzymujemy N przybliżonych współ-
czynników Fouriera c
N
n
, równych współczynnikom wielomianu trygonome-
trycznego interpolującego funkcję f w punktach t
k
= k
a
N
. Mamy dwie rów-
noważne formuły:
y
k
=
N −1
X
n=0
Y
n
ω
nk
N
,
k = 0, 1, 2, . . . , N − 1,
Y
n
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
,
n = 0, 1, 2, . . . , N − 1.
Przybliżone współczynniki szeregu Fouriera mają postać
c
n
≈ c
N
n
=
Y
n
,
gdy 0 ¬ n <
N
2
Y
n+N
,
gdy −
N
2
¬ n < 0.
Definicja: Dyskretną transformatą Fouriera (DFT) ciągu y
0
, y
1
, . . . , y
N −1
nazywamy ciąg liczbowy Y
0
, Y
1
, . . . , Y
N −1
dany wzorem
Y
n
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
,
n = 0, 1, 2, . . . , N − 1.
Piszemy
Y = (Y
0
, Y
1
, . . . , Y
N −1
) = F
N
((y
0
, y
1
, . . . , y
N −1
)).
Wniosek: F
N
: C
N
→ C
N
jest odwracalnym odwzorowaniem liniowym.
Odwzorowanie odwrotne dane jest wzorem
y
k
=
N −1
X
n=0
Y
n
ω
nk
N
,
k = 0, 1, 2, . . . , N − 1.
F
−1
N
zadaje się macierzą
Ω
N
= (ω
nk
N
) =
1
1
1
. . .
1
1
ω
N
ω
2
N
. . .
ω
N −1
N
1
ω
2
N
ω
4
N
. . . ω
2(N −1)
N
...
...
...
...
1 ω
N −1
N
ω
2(N −1)
N
. . . ω
(N −1)
2
N
.
3
Macierzą odwzorowania F
N
jest
Ω
−1
N
=
1
N
Ω
N
.
Uwaga: Pamiętając, że y
k
= f (k
a
N
) i f jest funkcją o okresie a, wygodnie
jest rozważać ciąg (y
k
) jako ciąg N -okresowy, określony dla k ∈ Z. Wzór
Y
n
=
1
N
M +N −1
X
k=M
y
k
ω
−nk
N
zadaje okresowe rozszerzenie ciągu Y
n
niezależnie od wyboru M :
Y
N l+n
=
1
N
M +N −1
X
k=M
y
k
ω
−(N l+n)k
N
=
1
N
M +N −1
X
k=M
y
k
ω
−nk
N
=
1
N
N j+s+N −1
X
k=N j+s
y
k
ω
−nk
N
=
=
1
N
N −1
X
k=0
y
N j+s+k
ω
−n(N j+s+k)
N
=
1
N
N −1
X
k=0
y
s+k
ω
−n(s+k)
N
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
= Y
n
.
Wynika stąd, że także współczynniki c
N
n
tworzą ciąg okresowy. Należy jed-
nak pamiętać, że c
N
n
aproksymują c
n
tylko dla −
N
2
¬ n <
N
2
, bowiem c
n
→ 0
przy n → ∞.
DFT ciągów rzeczywistych
W przypadku ciągów rzeczywistych przy wyznaczaniu DFT można zre-
dukować liczbę wykonywanych operacji o połowę, stosując do dwóch takich
ciągów pojedynczą transformatę zespoloną.
Chcemy wyznaczyć DFT dwóch ciągów rzeczywistych (x
k
) i (y
k
):
F
N
:(x
k
) → (X
n
),
F
N
:(y
k
) → (Y
n
).
Wiemy, że
X
N −n
=
X
n
i
Y
N −n
= Y
n
.
Niech z
k
= x
k
+ iy
k
i niech (Z
n
) oznacza transformatę ciągu (z
k
):
F
N
: (z
k
) → (Z
n
).
Na mocy liniowości
Z
n
= X
n
+ iY
n
,
4
przy czym X
n
i Y
n
niekoniecznie muszą być rzeczywiste, i wtedy
X
n
=
1
2
(Z
n
+ Z
N −n
),
Y
n
=
1
2
(Z
n
− Z
N −n
).
Wystarczy wyznaczyć te wartości dla n = 0, 1, . . . ,
N
2
, ponieważ wartości
dla n pomiędzy
N
2
i N − 1 pojawią się jako sprzężenia.
Zależność pomiędzy rzeczywistymi a przybliżonymi współczynni-
kami Fouriera
Załóżmy dla ułatwienia, że funkcja okresowa f daje się przedstawić w
postaci
f (t) =
∞
X
n=−∞
c
n
e
2iπn
t
a
i szereg współczynników jest absolutnie zbieżny:
∞
X
n=−∞
|c
n
| < +∞.
Możemy wówczas zmienić kolejność sumowania i sumować najpierw wyrazy
z indeksami mającymi ustaloną resztę modulo N, a potem sumować po tych
resztach. Biorąc t = k
a
N
, otrzymujemy:
f (k
a
N
) = y
k
=
∞
X
m=−∞
c
m
ω
mk
N
=
N −1
X
n=0
∞
X
q=−∞
c
n+qN
ω
nk
N
,
gdzie m = n + qN , n ∈ {0, 1, . . . , N − 1}, q ∈ Z. Otrzymujemy:
c
N
n
=
∞
X
q=−∞
c
n+qN
.
Powyższa równość wyraża przybliżone współczynniki poprzez współczynni-
ki właściwe i daje wzór na błąd aproksymacji:
c
N
n
− c
n
=
X
q6=0
c
n+qN
.
Widzimy więc, że dla ustalonego N im szybciej współczynniki c
n
zbiegają
do zera przy n → ∞, tym dokładniejsze będzie przybliżenie.
5
Własności DFT
Fakt: Jeśli F
N
: (y
k
) → (Y
n
), to
a) F
N
: (y
−k
) → (Y
−n
),
b) F
N
: (y
k
) → (Y
n
),
c) F
N
: (y
−k
) → (Y
n
).
Dowód:
a) Niech F
N
((y
−k
)) = (Y
0
n
). Wówczas
Y
0
n
=
1
N
N −1
X
k=0
y
−k
ω
−nk
N
=
1
N
0
X
k=−N +1
y
k
ω
nk
N
= Y
−n
.
b) Niech F
N
((y
k
)) = (Y
0
n
). Wówczas
Y
0
n
=
1
N
N −1
X
k=0
y
k
ω
−nk
N
=
1
N
N −1
X
k=0
y
k
ω
nk
N
= Y
−n
.
c) Niech F
N
((y
−k
)) = (Y
0
n
). Wówczas
Y
0
n
=
1
N
N −1
X
k=0
y
−k
ω
−nk
N
=
1
N
N −1
X
k=0
y
−k
ω
nk
N
=
1
N
0
X
k=−N +1
y
k
ω
nk
N
= Y
n
.
Wniosek: Jeśli F
N
: (y
k
) → (Y
n
), to zachodzą własności:
a) (y
k
) jest parzysty (nieparzysty) ⇐⇒ (Y
n
) jest parzysty (nieparzysty),
b) (y
k
) jest rzeczywisty ⇐⇒ Y
−n
= Y
n
dla każdego n ∈ Z,
c) (y
k
) jest parzystym ciągiem liczb rzeczywistych ⇐⇒ (Y
n
) jest parzystym
ciągiem liczb rzeczywistych,
d) (y
k
) jest nieparzystym ciągiem liczb rzeczywistych ⇐⇒ (Y
n
) jest nie-
parzystym ciągiem liczb czysto urojonych.
6
Fakt: Niech (x
k
) i (y
k
) będą dwoma zespolonymi ciągami o okresie N i niech
(X
k
) i (Y
k
) oznaczają ich DFT.
a) Transformata splotu cyklicznego, tzn. ciągu zdefiniowanego wzorem
z
k
=
N −1
X
q=0
x
q
y
k−q
,
k ∈ Z,
ma postać
F
N
: (z
k
) → (Z
n
= N X
n
Y
n
).
b) Transformata iloczynu ciągów (x
k
) i (y
k
) ma postać
F
N
: (p
k
= x
k
y
k
) → (P
n
=
N −1
X
q=0
X
q
Y
n−q
).
Dowód:
a) Z definicji
Z
n
=
1
N
N −1
X
k=0
x
q
y
k−q
ω
−nk
N
=
1
N
N −1
X
q=0
x
q
ω
−nq
N
N −1
X
k−0
y
k−q
ω
−n(k−q)
N
= N X
n
Y
n
.
b) Wyznaczamy odwrotną transformatę ciągu (P
n
):
p
k
=
N −1
X
n=0
N −1
X
q=0
X
q
Y
n−q
ω
nk
N
= x
k
y
k
.
Fakt: Jeśli F
N
: (y
k
) → (Y
n
), to
N −1
X
k=0
|y
k
|
2
= N
N −1
X
n=0
|Y
n
|
2
.
Dowód:
N
N −1
X
n=0
|Y
n
|
2
= N
Y
T
Y = Y
T
Ω
N
T
Ω
T
N
Y = (Ω
N
Y )
T
(Ω
N
Y ) = y
T
y =
N −1
X
k=0
|y
k
|
2
.
7
2 . Szybka transformata Fouriera
Wyznaczenie ciągu (Y
0
, Y
1
, . . . , Y
N −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 ω
j
N
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 F
N
, której ele-
mentami są pierwiastki z jedynki. Koszt FFT jest rzędu N log N .
Niech N = 2m, m ∈ N.
Y
n
=
1
2m
2m−1
X
k=0
y
k
ω
−nk
N
=
1
2m
m−1
X
k=0
y
2k
ω
−2nk
N
+ y
2k+1
ω
−n(2k+1)
N
=
=
1
2
1
m
m−1
X
k=0
y
2k
ω
−2nk
N
+ ω
−n
N
1
m
m−1
X
k=0
y
2k+1
ω
−2nk
N
,
ale
ω
−2nk
N
= ω
−2nk
2m
= e
2iπ
1
2m
(−2nk)
=
e
2iπ
1
m
−nk
= ω
−nk
m
,
czyli
Y
n
=
1
2
1
m
m−1
X
k=0
y
2k
ω
−nk
m
+ ω
−n
2m
1
m
m−1
X
k=0
y
2k+1
ω
−nk
m
=
1
2
P
n
+ ω
−n
2m
I
n
.
Zauważmy, że
P
m+n
=
1
m
m−1
X
k=0
y
2k
ω
−(m+n)k
m
=
1
m
m−1
X
k=0
y
2k
ω
−nk
m
= P
n
i podobnie
I
n+m
= I
n
oraz
ω
−(m+n)
2m
= e
−2iπ
1
2m
(m+n)
= ω
−n
2m
· e
−iπ
= −ω
−n
2m
.
Równości te prowadzą do algorytmu:
1
Krok 1. Obliczyć P
k
i ω
−k
N
I
k
.
Krok 2. Utworzyć Y
k
=
1
2
(P
k
+ ω
−k
N
I
k
).
Krok 3. Wyznaczyć Y
k+m
=
1
2
(P
k
− ω
−k
N
I
k
).
Kroki te należy wykonać dla k = 0, 1, . . . , m − 1.
P
k
i I
k
są dodatkowo dwoma niezależnymi DFT rzędu m =
N
2
:
F
N
2
:(y
0
, y
2
, . . . , y
2m−2
) → (P
0
, P
1
, . . . , P
m−1
),
F
N
2
:(y
1
, y
3
, . . . , y
2m−1
) → (I
0
, I
1
, . . . , I
m−1
).
Przy założeniu, że m jest parzyste, można więc algorytm powtórzyć. Jeśli
N = 2
p
, to iterujemy ten proces tak długo, aż dojdziemy do DFT rzędu 2.
Ma ona postać:
Y
0
=
y
0
+ y
1
2
,
Y
1
=
y
0
− y
1
2
.
Przykład: N = 8.
1. (y
0
, y
1
, . . . , y
7
) dzielimy na dwa ciągi długości 4:
(y
0
, y
2
, y
4
, y
6
),
(y
1
, y
3
, y
5
, y
7
).
2. Ciągi te dzielimy na ciągi długości 2:
(y
0
, y
4
),
(y
2
, y
6
),
(y
1
, y
5
),
(y
3
, y
7
).
3. Obliczamy dla każdego z tych ciągów P
0
i I
0
oraz wyznaczamy (Y
0
, Y
1
).
4. Przy użyciu formuł
Y
k
=
1
2
(P
k
+ ω
−k
N
I
k
)
Y
k+m
=
1
2
(P
k
− ω
−k
N
I
k
).
przechodzimy od wektora długości m do wektora długości 2m.
2
Dla N = 2
p
ogólny koszt tego algorytmu to
1
2
N (log
2
N − 2) + 1 dodawań
i N log
2
N 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];
3
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.
(x
n
)
n∈Z
, (h
q
)
q∈Z
- ciągi o wartościach zespolonych, okresowe o okresie N .
Splot cykliczny tych ciągów zdefiniowany jest wzorem
y
n
=
N −1
X
q=0
h
q
x
n−q
=
N −1
X
q=0
h
n−q
x
q
.
Jest to ciąg okresowy o okresie N . Przy ustalonym (h
q
) przekształcenie
zadane powyższym wzorem jest liniowym przekształceniem C
N
w siebie:
X → Y = HX,
X = (x
0
, x
1
, . . . , x
N −1
)
T
,
Y = (y
0
, y
1
, . . . , y
N −1
)
T
,
H =
h
0
h
N −1
h
N −2
. . . h
1
h
1
h
0
h
N −1
. . . h
2
h
2
h
1
h
0
. . . h
0
...
h
N −1
h
N −2
h
N −3
. . . h
0
.
Macierz tą nazywamy macierzą cykliczną.
Wyznaczenie splotu cyklicznego bezpośrednio z definicji wymaga wyko-
nania:
mnożenia zespolonego -
N
2
razy,
dodawania zespolonego - N (N − 1) razy.
Możemy jednak rozważać DFT (X
k
), (H
k
) i (Y
k
) odpowiednio ciągów (x
n
),
(h
n
) i (y
n
). Równanie zadające splot będzie miało wówczas postać
Y
k
= N H
k
X
k
,
k = 0, 1, 2, . . . , N − 1.
4
Jeśli założymy, że N = 2
p
, to wykonujemy teraz następujące kroki:
Kolejne kroki
Koszt
1. Wyznaczamy transformaty
[N (p − 2); 2N p]
F
N
: (x
n
) → (X
k
)
F
N
: (h
n
) → (H
k
)
2. Obliczamy iloczyn
[N ; 0]
Y
k
= N H
k
X
k
, k = 0, 1, . . . , N − 1
3. Wyznaczamy transformatę odwrotną
h
N
2
(p − 2); N p
i
.
F
−1
N
: (Y
k
) → (y
n
)
Koszt całkowity to
mnożenie zespolone –
N
2
(3 log
2
N − 4) razy,
dodawanie zespolone – 3N log
2
N 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
=
Y
k
, możemy zastąpić to obli-
czaniem transformaty rzędu
N
2
.
Całkowity koszt wynosił będzie
h
N
4
(3 log
2
N − 5);
N
2
(3 log
2
N − 1)
i
opera-
cji zespolonych, czyli
h
N (3 log
2
N − 5);
N
2
(9 log
2
N − 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 (x
n
)
n∈Z
, (h
n
)
n∈Z
będą dwoma sy-
gnałami nieokresowymi o zwartych nośnikach. W szczególności niech
x
n
= 0 jeśli n < 0 lub n M ,
h
n
= 0 jeśli n < 0 lub n Q (Q M ).
Chcemy wyznaczyć
y
n
=
Q−1
X
q=0
h
q
x
n−q
.
5
y
n
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 (x
n
) jest praktycznie „nieskończony”, a
nośnik filtru (h
n
) relatywnie mały. Są jednak metody pozwalające obniżyć
ten koszt, rozważa się w nich ciąg (x
n
) podzielony na mniejsze części.
6