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