11 - Lecture notes 11 PDF

Preview

Summary

Lecture notes from lecture 11. Lecturer: Dr. Sütő Zoltán...

Description

A gyors Fourier transzformáció

11. AZ FFT ÉS JELGENERÁTOROK

11.1 A gyors Fourier transzformáció (Fast Fourier Transform, FFT) Az FFT a diszkrét idejű Fourier transzformáció (DFT) kiszámításának egy hatékony algoritmusa. Ha fn egy (általában komplex értékű) idősorozat, ahol n=0,1,2,…(N-1), akkor ennek DFT-je, (azaz a spektrum): N −1

Fk = ∑ f n e

−j

2π N

N −1

nk

= ∑ f nW Nnk

n= 0

k = 0,1,2,….,(N-1)

(11.1.)

n =0

−j

2π

Ahol: WN = e N (11.2.) Az algoritmus akkor lesz hatékony, ha N-et a 2-nek egész számú hatványaként választjuk meg. A (11.1.)-szerinti összegzést bontsuk fel két részre, mely részekben a páros (2n), ill. a páratlan (2n+1) sorszámú mintákat használjuk csak fel: N −1

Fk = ∑ f nW

N −1 2

∑f

=

nk N

n =0

n =0

Fk =

W

2n

2 nk N

(N / 2 )−1

∑ g W( n

nk N / 2)

N −1 2

+ ∑ f 2n +1W N( 2 n+1) k = n= 0

+ WNk

( N / 2) −1

∑ h W( n

n =0

nk N / 2)

= Gk + WNk H k

(11.3.)

n=0

A (11.3.)-ban Gk és Hk mennyiségeket felfoghatjuk úgy, mint az N/2 hosszúságú gn és hn sorozatok DFT-jét, ahol: n = 0,1,2,….,(N/2-1) (11.4.) hn = f 2n +1 g n = f 2n A Gk és Hk kiszámításánál a k index a k =(0,1,2,…,(N/2-1) ) értékeket veheti fel, míg a (11.1.) kiszámításához az indexnek a k =(0,1,2,…,(N-1) ) tartományon kell futnia. Ezért az (N/2) ill. az ennél nagyobb index értékekre a (11.1.)-et írjuk fel az alábbi alakban:

Fk + N / 2 =

N −1 2

n( k +

∑ g W( n

n= 0

N ) 2

N / 2)

k+

+ WN

N N 2 −1 2

n( k +

∑ h W( n

n =0

N ) 2

(11.5.)

N / 2)

ahol: k =(0,1,2,…,(N/2-1) ) A (11.5.) összefüggésben szereplő tényezők egyszerűbben is kifejezhetők: W Nk +N / 2 = e

W

k +N / 2 (N / 2 )

=e

− j

−j

N 2π ( k+ ) N 2

2π N

2(k +

N ) 2

=e

−j

=e

2π k − jπ N

−j

e

2π N

2k

e

= −e

− j 2π

1/14

−j

=e

2π k N

−j

2π N

= −W Nk 2k

= W ( kN/ 2 )

(11.6.) (11.7.)

A gyors Fourier transzformáció

Ezzel a spektrum felső fele: N −1 2

Fk + N / 2 = ∑ g n W

nk N/ 2

−W

k N

n= 0

N −1 2

∑h W n

nk N /2

= Gk − WNk H k

(11.8.)

n =0

A (11.3.) és a (11.8.) összefüggéseket egy jelfolyam gráfban ábrázolva kapjuk az u.n. "pillangót" (butterfly):

Gk

Hk

Fk W Nk Fk+N/2 11.1. ábra A pillangó

Ezzel az átalakítással az N-pontos DFT kiszámítását visszavezettük két (N/2) fokszámú DFT kiszámítására. g0 = f0 g1 = f2

gn = f2n

N/2

gN/2-1 = fN-2

F0

G1

F1

Gk

Fk

GN/2-1

FN/2-1

H0

h0 = f1 h1 = f3

hn = f2n+1

G00

N

WN0

F0+N/2

H1

N/2 Hk

hN/2-1 = fN-1

N/2

F1+N/2

WNk

N/2 Fk+N/2

HN/2-1

FN-1

11.2. ábra Az N-pontos DFT visszavezetése két (N/2)-pontos DFT kiszámítására.

2/14

A gyors Fourier transzformáció

Természetesen a fokszám redukálásos eljárás tovább folytatható, mivel az N a 2-nek egész számú hatványa. Az alábbi ábra egy N=8 pontos DFT teljes kifejtését mutatja: f0 f4 f2 f6 f1

a0 a1

W

b1

0 2

W

f3

d0 d1

A1

G1

B1

0 4

G2

1 4

G3

W W

C0

c0 c1

G0

B0

b0

f5

f7

0 2

A0

0 2

W

H0

C1

0 2

W

F0 F1 F2 F3 0 8

W

1 8

H1

W

D0

W40

H2

W82

D1

1 4

H3

W83

W

F4 F5 F6 F7

11.3. ábra Az N=8 pontos FFT kiszámításának folyamatábrája Mint az ábrából látható minden egyes fokozatban N/2 darab komplex szorzást, összeadást ill. kivonást kell elvégezni, míg a fokozatok száma log2N. A műveletek összes száma: N sz = log 2 N (11.9.) 2 Ha a DFT-t a (11.1.) szerint számítanánk ki, akkor N2 komplex szorzást kellene elvégeznünk. A megtakarítás nagy N-eknél jelentős, pld. N=1024 esetén az FFT művelet igénye 5120, míg a közvetlen számolással 10242 ~106 adódik. További egyszerűsítést jelent a 0-s indexű pillangók számításánál a: WN0 = e− j 0 = 1 azonosság használata. B0 f2

(11.10.)

B1

f6

11.4. ábra Az első fokozat pillangói Az algoritmus taglalását az idősorozat decimálásával kezdtük (Decimation In Time) ezért ezt DIT eljárásnak nevezik. A DIT módszerben a kimeneti sorozat rendezett (az indexek növekvő sorrendben helyezkednek el), míg a bemeneti sorozat első pillantásra rendezetlennek tűnik. [f0-7]: f0, f1, f2, f3, f4, f5, f6, f7 [g0-3]: f0, f2, f4, f6, [h0-3]: f1, f3, f5, f7, [a0-1]: f0, f4, [b0-1]: f2, f6, [c0-1]: f1, f5, [d0-1]: f3, f7 3/14

A gyors Fourier transzformáció

A bemeneti sorrend: f0, f4, f2, f6, f1, f5, f3, f7 Az alábbi táblázat első két oszlopában ezt a sorrendet tüntettük fel decimális és bináris formában.

ndec

nbin

rbin

rdec

0

000

000

0

4

100

001

1

2

010

010

2

6

110

011

3

1

001

100

4

5

101

101

5

3

011

110

6

7

111

111

7

A táblázat harmadik oszlopába a második oszlop bináris adatainak fordítottját (reverse) írtuk. Ezt a számot decimálissá alakítva a rend helyreáll. Itt nem bizonyítjuk, de ez a rend más pontszámokra is érvényes. A DIT eljárásban tehát az input vektor feltöltését bit reverse címzéssel kell elvégezni. A legtöbb DSP-ben lehetőség van az ilyen címzésre. DSP implementáció esetén a fontos szempont, hogy a számítások közben elkerüljük a túlcsordulást. Ha az fn sorozat formátuma pld 16 TC Q15, ( -1 ≤ real(fn )< 1, -1 ≤ imag(fn )< 1 ), akkor a részeredmények és kimenet adatai biztosan nem férnek el ebben a formátumban. Ezért a (11.1.)-el ellentétben az FFT esetében a pillangókban egy 2-vel való osztást is végrehajtunk, így az összeadások után a formátum mindig ugyanaz maradhat. Ez log2N számú fokozat után N-nel történő osztásnak felel meg. DFT esetén: N −1

Fk = DFT { f n } = ∑ f nW Nnk

(11.11.)

n= 0

FFT esetén: 1 N

Fk = FFT { f n } =

Az FFT-ben a pillangók: 1 Fk = (G k + W Nk H k ) 2

N −1

∑f W n

(11.12.)

nk N

n= 0

és

4/14

Fk+ N / 2 =

1 (G − WNk H k ) 2 k

(11.13.)

A gyors Fourier transzformáció

Ha egy sorozat inverz diszkrét Fourier transzformáltját (IDFT-jét) kell kiszámítanunk, akkor azt elvégezhetjük a rendelkezésünkre álló (11.12.) szerinti FFT eljárással. A feladat: 2

π j kn 1 N− 1 (11.14.) f n = IDFT {Fk } = ∑ Fk e N N k =0 Képezzük a (11.14.) szerinti kifejezés konjugáltját! ( Az összeg konjugáltja a tagok konjugáltjának összege, a szorzat konjugáltja a tényezők konjugáltjának szorzata.):

fn =

−j 1 N −1 Fk e ∑ N k= 0

2π kn N

=

1 N −1 ∑ F W kn = FFT Fk N k =0 k N

{ }

(11.15.)

Az input vektort a bemeneti sorozat konjugáltjával töltjük fel (bit reverse címzéssel), majd az FFT eredményét ismételten konjugáljuk, miáltal megkapjuk a sorozat IDFT-jét.

11.2 Jelgenerátorok A jelgenerátorok más néven az oszcillátorok feladata periodikus jelek előállítása. Jellemzőjük a hullámalak és a frekvencia. Bemenő jelük általában nincs, a kimenet maga az előállított jel. 11.2.1 Táblázatos módszer Tetszőleges hullámformájú periodikus jel legegyszerűbben úgy állítható elő, hogy a hullám egy periódusának mintáit táblázatban tároljuk. Egy állapotváltozó szükséges (az n index) amivel a táblázatot cirkulárisan címezzük. Ha a táblázat mérete N, akkor a generált jel frekvenciája a mintavételi frekvencia N-ed része:

ω0 =

ωs

N Leggyakrabban szinuszos jel előállítása a feladat. Ha a jel frekvenciáját változtatjuk, akkor is célszerű egyetlen táblázatot használni. Legyen a kívánt jel x(nT): x(nT ) = Asin ω 0 nT = Asin ϕ (n) A következő mintavételi ütemben érvényes fázis rekurzívan is kifejezhető:

ϕ ( n + 1) = ω 0 nT = ϕ ( n) + ω 0 T = ϕ (n) + 2πω 0 / ω s Praktikus okokból érdemes az állapotváltozónak felvett fázist nem radiánokban mérni, hanem az alábbiak szerint normálni: ϕ ( n + 1) ω f ( n + 1) = N = f ( n) + N 0 = f ( n) + Δ f (11.16) ωs 2π

5/14

A gyors Fourier transzformáció

A normált fázis növekménye: Δ f = N

ω0 ωs

Ha a jel frekvenciája a mintavételi frekvenciával az

ω0 = m

ωs

N viszonyban van (ahol m egész), akkor a fázis növekmény Δf = m is egész. Ebből következően az f(n) is egész lehet, amit a táblázat címzésére közvetlenül felhasználhatunk (f(n) az index). A táblázatot a: ⎛ 2π ⎞ T sin [k ] = A sin ⎜ k k = 0,1,2,…,(N-1) ⎟ ⎝ N ⎠ adatokkal feltöltve a jel generálása a memória olvasásával már nagyon egyszerű. A fázis frissítése a növekmény hozzáadásával, majd az eredmény modulo N végrehajtásával kell, hogy megtörténjen, biztosítva azt, hogy az index ne mutasson ki a táblázatból a következő ütemben sem. f ( n + 1) = [ f (n ) + m] mod N 0 ≤ f ≤ (N −1) A modulo (maradék) képzés akkor egyszerű, ha az N méretet a 2 egész számú hatványának választjuk és az összeget (N-1) –el a DSP programban “maszkoljuk” (AND). Pld.: ha N=256 f(n+1) = [f(n) + m] AND [255] (255dec=0000 0000 1111 1111bin). Gyakran igény a kétfázisú (sin-cos) jel előállítása. Ekkor nem kell még egy oszcillátort készíteni, hanem elég a táblázatot folytatni még egy negyed periódus hosszan. A táblázatot az fik idextről olvasva a szinuszos, az f+N /4 indexről a koszinuszos eredményt kapjuk.

11.2.2 Táblázatos módszer interpolációval (szinuszos jelekre) A táblázatos módszer kétségtelen hátránya, hogy használatával csak diszkrét frekvenciájú jelek állíthatók elő. Ha a (11.16) szerinti fázisnövekményben megengedünk tört részt is, ettől a hátránytól megszabadulhatunk. Bontsuk fel a normalizált (f) fázist (e) egész és (t) tört részére:

f(n) = e(n) + t(n) ahol: 0 ≤ t(n) < 1. Legyen most az előállítandó jel az:

(11.17.)

⎡ 2π ( )⎤ ⎡ 2π ( ( ) ( ))⎤ y (n ) = Acos [ϕ (n )] = Acos ⎢ f n ⎥ = Acos⎢ e n +t n ⎥= ⎣N ⎦ ⎣ N ⎦ ⎡ 2π ⎤ ⎡ 2π ⎤ ⎡ 2π ⎤ ⎡ 2π ⎤ = A cos⎢ e (n )⎥ cos ⎢ t (n )⎥ − A sin ⎢ e(n )⎥ sin ⎢ t (n )⎥ ⎣N ⎦ ⎣N ⎦ ⎣N ⎦ ⎣N ⎦

Felhasználva, a kis szögek szinuszára és koszinuszára vonatkozó közelítéseket:

6/14

(11.18)

A gyors Fourier transzformáció

⎡ 2π ⎤ 2π t (n )⎥ ≅ t (n )...