Elég régen írtam az előző cikket, ezért most a bevezetőt nem fogom hosszú sorokra ereszteni. Az előző fejezetben megígértem, hogy egy kicsit körüljárom a Fourier-transzformáció algoritmusát és velejáróit. Vessük hát bele magunkat, egyébként sem egy sakk bajnokság jegyzőkönyvét olvashatja a Kedves Olvasó. (Helyreigazítás: A szerző egyáltalán nem gondolja a sakkot unalmasnak.)

Figyelem!
Az alábbi cikk tömegesen tartalmaz matematikai okfejtéseket,
melyek az olvasó lelkivilágára ártalmasak lehetnek!

Akkor matekozzunk egy kicsit!

Először is egy kicsit beszéljünk a komplex számokról. Majd egy idő után el is engedhetjük az információt (bár egy műszaki ember azért nem árt, ha birtokában van) de annak okán hogy ne a levegőbe beszéljek a következőkben, tisztázok egy-két dolgot. Természetesen a matematikai teljesség igénye nélkül fogok írni a dolgokról, így ha valakinek egy kicsit bántja a szemét, akkor azért elnézést kérek.

A komplex számokat a legegyszerűbben úgy tudjuk elképzelni, mint két összetartozó valós számot. Ebben az aspektusban felfoghatjuk úgy is, mint egy kétdimenziós vektort. Bevezetünk egy úgynevezett képzetes egységet teljesen abszurd és látszólag indokolatlan módon, mint aminek a négyzete mínusz egy \(\left(j^{2}=-1\right)\). A filozófiáról elég is ennyi, ami most fontos számunkra, hogy hogyan jelöljük és milyen alapvető műveleteket tudunk végezni a komplex számokkal.

Mint mondottam két része van a komplex számnak a valós (reális) és képzetes (imaginárius). Három alakja létezik a komplex számoknak ahhoz, hogy utóbbi kettőt megértsük, el kell fogadnunk, hogy a komplex számoknak létezik hossza (norma) és szöge (argumentuma). A listánkban az első alak az úgynevezett algebrai alak.

\begin{align}
z&=a+jb\\
\Re z&=a\\
\Im z&=b
\end{align}

Ahol a komplex szám valós és képzetes része \(a\) és \(b\). Ebből a kettőből képezhető a norma (\(r\)) és a szög (\(\varphi\)).

\begin{align}
r&=\sqrt{a^{2}+b^{2}}\\
\varphi&=\textrm{atan2}\left(a,b\right)
\end{align}

Ha ehhez még Leonhard Euler egy zseniális észrevételét hozzátesszük, miszerint \(e^{j\varphi}=\cos\varphi+j\sin\varphi\), akkor további két alakhoz jutunk. A trigonometrikus és a kényelmesebb Euler-alakhoz.

\begin{align}
z&=r\left(\cos\varphi+j\sin\varphi\right)\\
z&=re^{j\varphi}
\end{align}

A komplex számokon végezhető műveletek, amik számunkra fontosak, azok az összeadás, számmal való szorzás és esetleg a konjugálás.

\begin{align}
z&=a+jb\\
w&=c+jd\\
z+w&=\left(a+c\right)+j\left(b+d\right)\\
\mu z&=\mu a+j\mu b\\
z^{*}&=a-jb
\end{align}

Egyelőre hagyjuk abba a matematikát, mielőtt megüli a gyomrunkat, majd később még lesz belőle részünk bőven (mint említettem).

Az algoritmusok előkészítése

Mivel később komplex számokon kell majd matematikai műveleteket végeznünk, ezeket jobb, ha előkészítjük. A programkódokat C-ben fogjuk írni. Amelyhez (csak hogy ráérezzünk) magunk fogjuk megírni az apparátust. Ehhez definiáljuk a típust:

A manipulációs függvények prototípusa legyen az alábbi (csak mert):

A függvények definícióját a cikk végén található linken találja az Olvasó. Gyakorlatilag összegyűjtöttem őket egy “RynComplex.h” header fájlba. Nem csinálnak mást, mint a fentiekben bemutatott műveleteket. Ezt fogom használni a későbbiekben.

A Fourier-transzformáció az időmilliomosoknak való

A definíció szerinti algoritmus az adatsor hosszának négyzetével arányos, ugyanis az \(N\) darab együtthatót kell kiszámolni és minden egyes együtthatóhoz \(N\) darab szorzást kell elvégezni.

Ha az \(x[k]\) adatsor \(N\) hosszú, akkor az \(X_{p}^{F}\) spektrum szintén \(N\) hosszú és a \(p\)-edik Fourier együttható (igen, oka van annak, hogy ott az \(F\) betű):

\begin{align}
X_{p}^{F}=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x_{k}e^{-j\frac{2\pi}{N}pk}
\end{align}

A Fourier-transzformáció drága részét a trigonometrikus függvények kiszámítása jelenti. A képletben szereplő exponenciális az Euler-formula szerint számolható két szögfüggvény kiszámításával. Először is a szokásjoghoz tartva magunkat vezessük be az alábbi jelölést Joseph L. Walsh-ra emlékezve:

\begin{align}
e^{-j\frac{2\pi}{N}pk}=W^{pk}_{N}=\cos\left(\frac{2\pi}{N}pk\right)-j\sin\left(\frac{2\pi}{N}pk\right)
\end{align}

Egyelőre ne legyen drága semmi és írjuk meg a fenti exponenciálist kiszámító függvényt. Használva a már megírt header fájlt.

A spektrum kiszámításához használjuk tehát a definíciót:

Most már bizonyára mindenki tűkön ül, hogy működik-e a dolog! Ez akkor fog kiderülni, ha kipróbáljuk és ehhez írjunk is egy kis példaprogramot. Ahhoz, hogy ellenőrizni tudjuk az eredményeinket, hívjuk elő a Matlab-ot (gyakran fogom ezt tenni). Adjunk meg egy adatsort és számítsuk ki a spektrumát!

A példaprogramunk pedig ennek megfelelően az alábbi eredményt adja:

 

Ha az Olvasó szemfülesebb, akkor észre vehette, hogy az általunk kapott eredmény nem pontosan egyezik meg a Matlab által produkáltakkal. Ennek az oka nagyon egyszerű. A Fourier-transzformációnál az exponenciális szorzó előjele az általunk definiáltak szerint negatív. A visszatranszformálásnál pedig pozitív. A matematikusok ezt általában fordítva használják. Ha a bemeneti adatsor valós, akkor az előjel felcserélése ezt a kis különbséget okozza.

Gondok

Ha egy kicsit lemászunk az elefántcsonttoronyból, akkor rájöhetünk, hogy ezzel több baj is van.

1. A lépésszám probléma elég nagy gond, főleg akkor, ha egy ügyes csellel a feldolgozandó minták számát többszörösére növeljük. Már körülbelül ezer minta esetén (ami valljuk be nem olyan sok), az elvégezendő műveletek száma milliós nagyságrendbe esik. Ezzel jó lenne kezdeni valamit, hiszen a mi eszközünk az adatokat “real time” dolgozza fel, tehát két minta között el kell végezni a teljes transzformációt és a kommunikációt is a számítógéppel. Szerencsére a hangfrekvenciás sávban ez annyira nem drasztikus és mint mondtam a minőségen tudunk mesterségesen javítani.
2. A Fourier-transzformáció a valós adatsorra egy komplex adatsort ad válaszul. Mint említettem ez azt jelenti, hogy az eredmény eltárolásához kétszer annyi hely kell, mint a bemenethez. Ez persze nem teljesen igaz, mert matematikai fondorlatosságunknak köszönhetően egyébként is elég csak a felét átvinni. Ennek ellenére kicsit megbütyköljük a transzformációt, elnevezzük valaki másról és kis extra tulajdonságot is szerzünk ezzel.
3. A mi kis eszközünk nincs arra felkészítve, hogy trigonometrikus függvényeket számolgassunk vele, ezért valahogyan meg kell oldanunk azt, hogy a trigonometrikus függvények értékei a rendelkezésünkre álljanak, méghozzá minél hamarabb.

Cooley-Tukey algoritmus

Vegyünk észre egy-két érdekességet az exponenciális egységvektorokkal kapcsolatban:

\begin{align}
W_{N}^{pk}&=e^{-j\frac{2\pi}{N}pk}\\
W_{N}^{p\left(2k+1\right)}&=W_{N}^{2kp}W_{N}^{p}\\
W_{N}^{2kp}&=W_{\frac{N}{2}}^{kp}\\
W_{N}^{pk}&=-W_{N}^{N-pk}
\end{align}

Nézzük meg még egyszer a Fourier-transzformáció definícióját (kicsit egyszerűsítettem a jelölést, hogy kényelmesebb legyen a szemnek):

\begin{align}
X_{p}^{F}&=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x_{k}W_{N}^{pk}
\end{align}

Szedjük szét a summát \(k\) páros, illetve páratlan volta szerint.

\begin{align}
X_{p}^{F}&=\frac{1}{\sqrt{N}}\left(\sum_{k=0}^{\frac{N}{2}-1}x_{2k}W_{N}^{2pk}+\sum_{k=0}^{\frac{N}{2}-1}x_{2k+1}W_{N}^{p(2k+1)}\right)\\
X_{p}^{F}&=\frac{1}{\sqrt{N}}\left(\sum_{k=0}^{\frac{N}{2}-1}x_{2k}W_{N}^{2pk}+W_{N}^{p}\sum_{k=0}^{\frac{N}{2}-1}x_{2k+1}W_{N}^{2pk}\right)\\
X_{p}^{F}&=\frac{1}{\sqrt{N}}\left(\sum_{k=0}^{\frac{N}{2}-1}x_{2k}W_{\frac{N}{2}}^{pk}+W_{N}^{p}\sum_{k=0}^{\frac{N}{2}-1}x_{2k+1}W_{\frac{N}{2}}^{pk}\right)\\
\end{align}

Ha jobban megnézzük az utolsó képletet, akkor kiderül, hogy az eredeti adatsor Fourier-transzformáltját ki lehet számolni az adatsor páros és páratlan elemeiből alkotott adatsor transzformáltjaiból. Jelöljük a páros elemekből alkotott sor transzformáltját \(E_{p}^{F}\)-vel, a páratlanok transzformáltját pedig \(O_{p}^{F}\)-vel. Ekkor a kifejezésben lévő a summákra bevezetett jelölésekkel a transzformáció:

\begin{align}
X_{p}^{F}&=\frac{1}{\sqrt{2}}\left(E_{p}^{F}+W_{N}^{p}O_{p}^{F}\right)
\end{align}

Ez önmagában még nem alkalmas az érvágáson kívül másra, azonban ne csüggedjünk, nézzük tovább! Az adatsor transzformáltjának első \(\frac{N}{2}\) elemét el fogom most nevezni \(X_{p}^{F\textrm{(bal)}}\), a másodikat \(X_{p}^{F\textrm{(jobb)}}\) jelöléseknek. Ennek köszönhetően (és kihasználva még egy fent említett azonosságot a Walsh vektorokra) egy érdekességet vehetünk észre:

\begin{align}
X_{p}^{F\textrm{(bal)}}&=\frac{1}{\sqrt{2}}\left(E_{p}^{F}+W_{N}^{p}O_{p}^{F}\right)\\
X_{p}^{F\textrm{(jobb)}}&=\frac{1}{\sqrt{2}}\left(E_{p}^{F}-W_{N}^{p}O_{p}^{F}\right)
\end{align}

Most már érdekesebb a dolog, hiszen azt látjuk, hogy a transzformáció “drága” részét kétszer is fel tudjuk használni, ha egyszer kiszámoltuk. Mindemellett egy rekurzív képlet is kirajzolódik a szemünk előtt.

A fenti műveletre egy kis ábrát szokás rajzolni. Amikor \(A\) és \(B\) mennyiségekből a fenti módon képzünk két mennyiséget:

Ezt az ábrát pillangó diagramnak nevezzük. Ezzel egy egyszerű ábra a nyolc elemű adatsor Fourier-transzformációjára:

 

Az rekurzív algoritmus tesztelése némiképp felesleges lenne, mert nem azt fogjuk használni. A rekurziót meg lehet valósítani iteráció (ciklus) segítségével is, ami a memóriafelhasználást csökkenti. Ehhez azonban a bemenetet (mivel nincs rekurzió) rendezni kell egy úgynevezett bit-reversal algoritmus segítségével.

Egy kis számolás után kiderül, hogy a definíció szerinti DFT számítási idejével \(\mathcal{O}(N^{2})\) szemben az Cooley-Tukey algoritmus lépésszáma \(\mathcal{O}(N\log N)\). Ez például egy \(1024\) elemből álló adatsor esetében \(99{,}5\%\)-os megtakarítást jelent. Ez már annyira gyors, hogy ezt az algoritmust is elnevezték Gyors Fourier-transzformációnak (FFT). Itt megemlíteném, hogy az FFT egy gyűjtő fogalom. Olyan algoritmusoknak a csoportja, amelyek a Fourier-transzformációt gyorsabban végrehajtják az eredeti elképzelésnél. Szerepelhet még erre az algoritmusra a DIT (Decimation In Time). Ez arra utal, hogy az időbeli koordináták páros-páratlan csoportosításával történt a lépésszám redukálása. Lehetne olyan megoldást is találni, hogy az időtartománybeli elemeket középen bal-jobb csoportra osztjuk. Ekkor úgynevezett DIF (Decimation In Freqency) algoritmust kapunk. Fontos, hogy ekkor az eredmény páros-páratlan csoportokra lesz osztva, de ezzel most mi nem fogunk foglalkozni. Egy ilyen algoritmus például a Sande-Tukey algoritmus.

Elég bőségesre sikeredett eddig az írás, ezért inkább itt most elvágom és hagyom az Olvasót egy picit lenyugodni. A következő részekben még egy kicsit folytatom a matematikai eszmefuttatást és ígérem megoldást találunk a fent felsorolt problémákra.

A cikkben megjelenő kódrészletek az alábbi linken elérhetők:
https://github.com/patrikradvanyi/SEM_DSP