Szybka transformacja Fouriera (FFT – Fast Fourier Transform)

Szybka transformacja Fouriera (FFT – Fast Fourier Transform)

Plan wykładu:

1. Transformacja Fouriera, iloczyn skalarny 2. DFT - dyskretna transformacja Fouriera 3. FFT – szybka transformacja Fouriera

a) algorytm PFA

b) algorytm Cooleya-Tukeya (radix-2) c) wielowymiarowe FFT

4. Przykłady zastosowań: mnożenie wielomianów, odszumianie sygnału, rozwiązanie równania Poissona, całkowanie.

2 Jeśli funkcja f(x) jest okresowa wówczas

zamiast wielomianów do jej interpolacji (aproksymacji) lepiej użyć wielomianów trygonometrycznych tj. rozwinąć funkcję w szereg Fouriera. Dla funkcji okresowej o okresie 2π:

Funkcję możemy też zapisać w postaci zespolonego szeregu Fouriera

f (x) = a

2 +

X

(a

cos(kx) + b

sin(kx))

a

= 1

¼

Z

¡¼

f (t)cos(kt)dt

b

= 1

¼

Z

¡¼

f (t)sin(kt)dt

Jeśli funkcja f(x) jest rzeczywista wówczas „zwykły”

szereg Fouriera jest częścią

rzeczywistą zespolonego szeregu Fouriera:

Dla ciągu współczynników

definiujemy

Co prowadzi do zależności pomiędzy szeregiem rzeczywistym i zespolonym

[a

]

[b

]

b

= 0 a

= a

b

= ¡b

f (x) »

X

f (k)e b

f (k) = b 1 2¼

Z

¡¼

f (t)e

dt

f (k) b = 1 2¼

Z

¡¼

f (t)[cos(kt) ¡ Isin(kt)]dt

= 1

2 (a

¡ Ib

)

k ¸ 0

c

= 1

2 (a

¡ Ib

)

a

2 +

X

(a

cos(kx) + b

sin(kx)) =

X

c

e

Funkcje

generują ciąg ortonormalnych funkcji w zespolonej przestrzeni Hilberta.

Iloczyn skalarny w tej przestrzeni:

Dysponując tablicą wartości funkcji f i g w węzłach siatki, iloczyn wewnętrzny

można zapisać w postaci dyskretnej:

hf; gi = 1 2¼

Z

¡¼

f (x)g

(x)dx

Własności iloczynu wewnętrznego/skalarnego:

oraz związek z normą euklidesową

Dla każdego

Dowód

hf; fi

¸ 0

hf; gi

= hg; fi

kf k

= p

hf; fi

hf; gi 1

X

N ¸ 1

h®f + ¯g; hi

= ® hf; hi

+ ¯ hg; hi

E

(x) = e

; k = 0; §1; §2; : : :

hE

; E

i = 1 2¼

Z

¡¼

E

(x)E

(x)dx

= 1

2¼

Z

¡¼

e

e

dx

= 1

2¼

Z

¡¼

e

dx

= 1

2¼

e

I(k ¡ n)

¯ ¯

¯ ¯

x=¼ x=¡¼

= 0

hE

; E

i

=

½ 1

2 Z Z Z 0

2 Z = Z Z

1 N

N

X

E

µ 2¼j N

¶ E

µ 2¼j N

¶

=

= 1 N

N

X

h e

i

4 Dla pozostałych przypadków można się posłużyć

wyrażeniem na sumę szeregu:

co daje

ze względu na postać licznika.

Funkcje E_k(x) tworzą ciąg ortogonalnych (ortonormalnych) jednomianów

eksponencjalnych, z których można utworzyć wielomian:

Wielomian eksponencjalny może posłużyć do interpolacji funkcji f(x).

N

X

¸

= ¸

¡ 1

¸ ¡ 1 ; ¸ 6= 1

Załóżmy, że jej wartości są określone na siatce zbudowanej z równoodległych węzłów:

Wielomian interpolujący ma wówczas postać

Współczynniki znajdziemy licząc iloczyny skalarne (lewa i prawa strona) z kolejnymi jednomianami E_m

Ciąg współczynników c_m wyznaczanych zgodnie z powyższym wzorem definiuje dyskretną

transformatę Fouriera (DFT – to wynik przekształcenia).

W metodzie najmniejszych kwadratów wielomian ten może posłużyć do

aproksymacji funkcji f(x) gdy stopień wielomianu aproksymującego jest mniejszy od N-1.

f (x) = P (x) =

N

X

c

E

(x) x

= 2¼j

N ; j = 0; 1; : : : ; N ¡ 1 e

= 1 , k ¡ m

N 2 Z Z Z

e

¡ 1

e

¡ 1 = 0

P (x) =

N

X

c

e

=

N

X

c

(e

)

=

X

c

E

(x)

hf; E

i =

N

X

c

hE

; E

i =

N

X

c

±

= c

f (x) = P (x) =

N

X

hf; E

iE

DFT można zapisać wykorzystując postać macierzową

Transformatę można znaleźć wykonując „tylko” mnożenie wektora przez macierz. Ale w ten sposób należy wykonać O(N²) operacji arytmetycznych.

Jednakże macierz E ma specyficzną postać – jej elementy są ze sobą ściśle powiązane co można wykorzystać w celu

zmniejszenia nakładu obliczeń.

Dzięki FFT liczba wykonywanych operacji może zmaleć do wartości O(Nlog₂N).

0 B B B @

E

(x

) E

(x

) : : : E

(x

) E

(x

) E

(x

) : : : E

(x

)

.. . .. . . .. .. .

E

(x

) E

(x

) : : : E

(x

) 1 C C C A

0 B B B @

f (x

) f (x

)

.. . f (x

)

1 C C C A =

0 B B B @

c

c

.. . c

1 C C C A Ef = c

N N

N ¢ log

(N )

1024 1048576 10240

4096 16777216 49152

16384 268435456 229375

6 FFT z rozkładem na czynniki pierwsze

(PFA - Prime Factor Algorithm)

Liczbę naturalną możemy zapisać jako

gdzie: r_i są liczbami pierwszymi przykład

Idea algorytmu PFA polega na zastąpieniu

obliczeń DFT w jednym wymiarze (skala N²), na obliczeniu iloczynu p transformat DFT

skalujących się jak

Algorytm PFA zmniejsza nakład obliczeń z N² do N(r₁+r₂+...+r_p)

Zakładamy, że funkcja f(x) jest stablicowana w n węzłach

Funcję rozwijamy w bazie wielomianów eksponencjalnych

N = r

r

: : : r

r

; r

; : : : ; r

; N = 56 = 2 ¢ 2 ¢ 2 ¢ 7

x

= 2¼j

n ; j = 0; 1; : : : ; n ¡ 1

Wykorzystując własności wielomianów, współczynniki rozwinięcia zapisujemy w postaci

Oznaczenia:

Rozkładamy liczbę węzłów na iloczyn liczb pierwszych:

i wprowadzamy kolejne oznaczenia:

n = r

r

: : : r

n

= r

r

: : : r

=

Y

r

º = 0; 1; : : : ; p ¡ 1

n

= 1

n = r

r

: : : r

n

F (x) =

n

X

c

e

c

= 1 n

n

X

f (x

)e

; k = 0; 1; : : : ; n ¡ 1

w = exp ¡

¡

¢ a

=

f ¡

n

¢

¾

c

=

n

X

a

w

Zmienną k zapisujemy w postaci:

i podobnie ułamek j/n:

k = ®

n

+ ®

n

+ : : : + ®

n

®

2 f0; 1; : : : ; r

¡ 1g

®

2 f0; 1; : : : ; r

¡ 1g .. .

®

2 f0; 1; : : : ; r

¡ 1g

l

2 f0; 1; : : : ; r

¡ 1g l

2 f0; 1; : : : ; r

¡ 1g

.. .

l

2 f0; 1; : : : ; r

¡ 1g

k

=

X

®

n

j

n

=

p¡1

X

i=º

l

n

k

< n

; º = 0; 1; : : : ; p

j

n = l

n

+ l

n

+ : : : l

n

M 2 Z Z Z k ¢ j

n =

Ã

X

i=1

®

n

! Ã

X

º=0

l

n

!

=

p¡1

X

º=0

l

n

Ã

X

i=º+1

®

n

!

+ M

=

p¡1

X

º=0

l

k

n

+ M

8 Wykorzystujemy uzyskany wynik do

obliczenia w^kj

Chcemy znaleźć wartość współczynników c_k

Wykorzystujemy teraz zależność pomiędzy wskaźnikiem j a l₁,l₂,...,l_p

w sumowaniu. Sumę po j możemy zapisać jako

ponieważ każdą wartość j realizuje odpowiednia kombinacja wskaźników l₁,l₂,...,l_p.

To przejście pozwala zapisać współczynnik c_k

jako iloczyn p transformat DFT jednowymiarowych.

j = l

+ l

r

+ : : : + l

r

r

: : : r

n¡1

X

j=0

!

r

X

r

X

: : :

r

X

w

= exp µ

¡2¼I kj n

¶

= exp Ã

¡2¼I

Ã

X

º=0

l

k

n

+ M

!!

=

p

Y

exp µ

¡2¼I l

k

n

¶

=

p

Y

w

c

=

n

X

f (x

) n exp

µ

¡Ik 2¼j n

¶

=

n

X

f (x

)

n w

Startujemy od obliczenia jednowymiarowego DFT dla wartości funkcji w węzłach, których indeksy zależą od aktualnych wartości l₁,l₂,...,l_p

Takich transformat będzie n/r_p , a wyznaczenie każdej z nich wiąże się z nakładem obliczeń rzędu (r_p)²

następnie obliczamy

Po wyznaczeniu w ten sposób p transformat dostajemy żądany współczynnik c_k (procedurę powtarzamy dla każdej wartości k).

c

(l

; l

; : : : ; ®

) =

r

X

c

(l

; l

; : : : ; l

)w

c

(l

; l

; : : : ; ®

; ®

) =

r_p¡1

X

c

(l

; l

; : : : ; l

; ®

)w

c

= c

(®

; ®

; : : : ; ®

; ®

) =

r

X

c

(l

; ®

; : : : ; ®

; ®

)w

c

= c(l

; l

; : : : ; l

)

=

r

X

r

X

: : :

r

X

c

(l

; l

; : : : ; l

)w

w

: : : w

c

= f (x

)

n ; j = l

+ l

r

+ : : : + l

r

r

: : : r

10 Algorytm radix-2

Najprostszy algorytm FFT to radix-2 (Cooley- Tukey) opracowany w latach 60 XX wieku w celu szybkiej analizy danych sejsmologicznych.

Naszym zadaniem jest obliczenie

współczynników transformaty Fouriera (DFT) c_k, ale wykonując jak najmniej obliczeń.

Zakładamy że całkowita liczba węzłów jest potęgą 2:

x

= 2¼ N j

j = 0; 1; 2; : : : ; N ¡ 1 N = 2

; r 2 N N N

Osobno grupujemy składniki a) parzyste

b) nieparzyste

j = 2m j = 2m + 1

c

= hE

; f i =

N

X

f (x

)E

(x

)

=

N

X

f (x

)exp ( ¡Ix

k)

=

N

X

f

exp µ

¡I 2¼ N jk

¶

c

=

N 2 ¡1

X

m=0

f

exp µ

¡I 2¼

N (2m)k

¶

+

N 2 ¡1

X

m=0

f

exp µ

¡I 2¼

N (2m + 1)k

¶

c

=

N 2 ¡1

X

m=0

f

exp µ

¡I 2¼

N=2 mk

¶

+exp µ

¡I 2¼ N k

¶

X

f

exp µ

¡I 2¼

N=2 mk

¶

c

= p

+ '

q

Korzystamy teraz z okresowości wyrazów p_k oraz q_k:

Natomiast czynnik fazowy ma następującą własność:

p

= p

q

= q

Uwagi:

a) współczynniki p_k oraz q_k można wyliczyć dzięki DFT nakładem O(N/2)²=O(N²/4)

b) dodatkowo oszczędzamy czas wyznaczając tylko współczynniki dla

ponieważ

Kolejnym krokiem w FFT jest podział sum w p_k oraz w q_k na sumy zawierające tylko elementy parzyste i nieparzyste. Po podziale liczba elementów w każdej z dwóch powstałych sum jest dwukrotnie mniejsza niż w elemencie macierzystym.

Proces rekurencyjnego podziału kończymy gdy liczba elementów jest równa 1.

k < N 2

c

= 8 <

:

p

+ '

q

; k <

p

2

¡ '

q

2

; k ¸

p

=

N 2 ¡1

X

m=0

f

exp µ

¡I 2¼

N=2 mk

¶

q

=

N 2 ¡1

X

m=0

f

exp µ

¡I 2¼

N=2 mk

¶

'

= exp µ

¡I 2¼ N k

¶

'

= exp µ

¡I 2¼ N

µ

k + N 2

¶¶

= exp µ

¡I 2¼ N k

¶ exp

µ

¡I 2¼ N

N 2

¶

= ¡exp µ

¡I 2¼ k

¶

= ¡'

12 Inne algorytmy FFT

1) Algorytm Winograda/Radera – DFT jest sformułowane w postaci cyklicznego splotu. Modyfikacja Winograda algorytmu FFT działa gdy N=p^m , p -liczba pierwsza.

2) Split-radix – modyfikacja algorytmu Cooleya-Tukeya. W każdym kroku DFT jest wyrażana jako suma DFT dla N/2 oraz dwóch DFT dla N/4. Jest to najszybszy algorytm FFT.

3) DST (discrete sine transform) oraz DCT (discrete cosine transform) – transformaty sinusowa i kosinusowa, opłaca się je stosować gdy transformację przeprowadzamy na funkcjach rzeczywistych. Unikamy w ten sposób operacji na liczbach zespolonych co jest kosztowne.

Wielowymiarowa FFT

Współczynnki wyznacza się stosując algorytm

jednowymiarowego FFT kolejno dla każdego z wymiarów.

N = N

N

: : : N

c

=

N₁

X

j1=0 N₂

X

j2=0

: : :

N_d

X

jd=0

f

exp µ

¡i2¼

µ j

k

N

+ j

k

N

+ : : : + j

k

N

¶¶

Zastosowania FFT:

1) interpolacja, aproksymacja 2) szybkie mnożenie

3) cyfrowe przetwarzanie sygnału (np. odszumianie – widmo częstotliwości) 4) kompresja danych

5) analiza sygnałów czasowych (korelacja, splot)

6) rozwiązywanie równań różniczkowych (rów. Poissona)

14 Szybkie mnożenie wielomianów przy

użyciu FFT

Chcemy obliczyć iloczyn dwóch wielomianów

Jeśli stopnie wielomianów są różne to je

wyrównujemy dodając do wielomianu niższego stopnia współczynniki równe 0.

Iloczyn wielomianów

Dokonujemy reindeksacji wskaźników

P (x) =

n

X

a

x

Q(x) =

n

X

b

x

i + j = k ! j = k ¡ i R(x) = P (x)Q(x) =

n

X

a

x

n

X

b

x

=

n

X

a

b

x

R(x) =

2n

X

n

X

a

b

x

=

2n

X

c

x

c

=

n

X

a

b

Jeśli współczynniki wielomianów a_i oraz b_i potraktujemy jako współrzędne wektorów

to wektor c jest ich splotem:

Korzystając z definicji tranformacji Fouriera dla splotu funkcji możemy zapisać

a a a = [a

; a

; : : : ; a

] bbb = [b

; b

; : : : ; b

]

ccc = a a a ¤ bbb

ccc = F F T

h

F F T (~ a ~ a a)F F T (~b~b~b) ~ i

~

a~ a ~ a = [a

; a

; : : : ; a

; a

; : : : ; a

]

~b~b~b = [b

; b

; : : : ; b

; b

; : : : ; b

] c

= 0

a

; b

= 0 , i > n ¡ 1

Filtracja sygnału

FFT

Dyskryminacja szumu

FFT^-1

f (x) = cos(x) + cos(2x) + cos(3x)

16 Rozwiązywanie równania Poissona (2D, 3D)

dokonujemy transformacji (FFT) całego równania (do przestrzeni odwrotnej)

skąd już łatwo wyznaczymy V_k

i gotowe rozwiązanie (w przestrzeni rzeczywistej)

Uwaga: musimy jeszcze uwzględnić warunki brzegowe (WB) .

1) jeśli WB są typu Dirichleta

to wykonujemy DST-FFT dla wnętrza obszaru, a WB uwzględniamy dokonując transformat

potencjału na brzegach – wyniki dodajemy 2) jeśli WB są typu Neumanna

r

V

= ½

k

V

= ½

V

= ½

k

V

= F F T

fV

g

V

j

= V

6= 0

@V

@~n j

= 0

to wówczas stosujemy DCT-FFT dla wnętrza obszaru. WB są automatycznie spełnione.

Całkowanie

Chcemy obliczyć całkę oddziaływania dwóch gęstości ładunku/materii

która zawiera osobliwość. Całkę możemy zapisać nieco inaczej

Potencjał V(r₁) jest splotem dwóch funkcji:

gęstości i funkcji f. Można więc wykorzystać tu twierdzenie o splocie i jego transformacie:

Transformata odwrotna ostatniego iloczynu daje poszukiwany potencjał.

F F T fV g = F F T f½

¤ fg

= F F T f½

g ¢ F F T f½

g j~r

¡ ~r

j

= f (~r

¡ ~r

) = f (~r

¡ ~r

)

C = Z

a

d~r

Z

a

d~r

½

(~r

)½

(~r

) j~r

¡ ~r

j

C = Z

a

d~r

½

(~r

)V (~r

) V (~r

) =

Z

d~r

½

(~r

)f (~r

¡ ~r

)