Az előző cikkben megoldottuk, hogy gyorsabban ki tudjuk számítani a spektrumot egyelőre csak rekurzív módon. Mielőtt átírnánk a rekurziót iterációra oldjuk meg, hogy ne kelljen a komplex számokkal bajlódni.

Ralph Hartley transzformációja

A Fourier-transzformáció a valós bemeneti adatsor esetén is komplex eredményt produkál. Mint láttuk, egy komplex szám eltárolásához két valós számot kell összetartozóan eltárolni. A kimeneti adatsor így kétszer annyi adatmennyiség, mint a bemeneti. Persze ki lehet használni, hogy valós bemenet esetén elég csak az adatsor első felét átvinni, de még mindig nehézség az összetartozó adatok tárolása, mozgatása és a közöttük történő műveletvégzés. A fordító ugyanis nem támogatja a komplex számok kezelését, így azt nekünk kell megoldani szoftveresen egy nem objektum orientált nyelven.

Ralph Hartley egy olyan transzformációt javasolt, ami valós bemenet esetén valós kimenetet eredményez, továbbá öninvertáló. A Hartley-transzformáció olyan, hogy ha az adatsoron egymás után kétszer elvégezzük minden változtatás nélkül, akkor visszakapjuk a kiindulási adatokat.

A definíciója nagyon hasonlít a Fourier-transzformációhoz:

\begin{align}
X_{p}&=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x_{k}\left(\cos\left(\frac{2\pi kp}{N}\right)+\sin\left(\frac{2\pi kp}{N}\right)\right)
\end{align}

A transzformáció örökli ihletője tulajdonságait, többek között a DIT és DIF algoritmusoknak (melyeket az előző cikkben láthatott az Olvasó) is beadja a derekát. Ezeket az FFT algoritmusok mintájára FHT algoritmusoknak nevezik. Sajnos a trigonometrikus számításokat nem spórolhatjuk meg vele, de könnyebb bánni a valós számokkal a már említett okok miatt.

Természetesen a transzformáció eredménye nem a spektrum, de van lehetőség a spektrum meghatározására. Könnyen belátható az alábbi két formula (talán az előző cikkben elég volt a matematikai túlbuzgásból):

\begin{align}
\Re{X_{p}^{F}}&=\frac{1}{2}\left(X_{p}^{H}+X_{N-p}^{H}\right)\\
\Im{X_{p}^{F}}&=-\frac{1}{2}\left(X_{p}^{H}-X_{N-p}^{H}\right)
\end{align}

A továbbiakban elhagyom a Hartley-transzformáltról az azt jelző felső indexet. Felmerül a kérdés, hogy feltétlenül ki kell-e számolni a Fourier-transzformáltat? Természetesen nem tettem volna fel a kérdést így, ha a válasz igen lenne. Nekünk igazából az amplitúdóspektrumra van szükségünk, amit a Fourier-transzformáltak normája. Most már viszket a bőröm egy kis levezetéshez, de most sem fogom túlzásba vinni:

\begin{align}
||X_{p}^{F}||&=\sqrt{\left(\Re{X_{p}^{F}}\right)^{2}+\left(\Im{X_{p}^{F}}\right)^{2}}=\dots=\frac{1}{\sqrt{2}}\sqrt{X_{p}^{2}+X_{N-p}^{2}}
\end{align}

A mi kis eszközünknek a feladata, hogy ezt a normát megszorozza egy pozitív valós számmal (negatívval nincs értelme, mert a norma mindig nagyobb vagy egyenlő, mint nulla).

\begin{align}
\mu||X_{p}^{F}||&=\mu\sqrt{\left(\Re{X_{p}^{F}}\right)^{2}+\left(\Im{X_{p}^{F}}\right)^{2}}=\dots=\frac{1}{\sqrt{2}}\sqrt{\left(\mu X_{p}\right)^{2}+\left(\mu X_{N-p}\right)^{2}}
\end{align}

Ez a sor azért nagyon jó, mert igazából a spektrumot tudjuk úgy is kezelni, hogy megkerüljük a Fourier-transzformációt. Ezzel teljesen kiküszöbölve a komplex számokat, továbbá a transzformáció irányát (ez Fourier-transzformáció esetén az exponenciális vektor kitevőjének előjele) sem kell figyelembe venni, mivel (mint említettem) a Hartley-transzformáció öninvertáló. Természetesen a trigonometrikus számításokkal még foglalkozni kell, de nem kell mindent egyszerre megoldani.

Csakhogy legyen már végre valami sikerélményünk, csináljuk meg a lassú DHT-t (Diszkrét Hartley-transzformáció).

Az előző cikkben szereplő adatsorra teszteljük le a transzformációt. Lássuk milyen eredményt szolgáltat.

Mint látható a Hartley-transzformációból számított spektrum megegyezik a Fourier-transzformációból számolttal.

Áttérés a rekurzív algoritmusról az iteratív algoritmusra

A Hartley-transzformáció örökölte a Cooley-Tukey algoritmussal való megvalósíthatóságot. Ebben az esetben a pillangók egy kicsit másként alakulnak. A hosszas bizonyítgatásokat ezúttal nem részletezem (most nem olyan egyszerű, mint a Fourier-transzformációban), azonban a forrásfájlok között található a bizonyítás, ha érdekli az Olvasót. A rekurzív formulákra az alábbiak adódnak:

\begin{align}
X_{p}^{\textrm{(bal)}}&=\frac{1}{\sqrt{2}}\left[E_{p}+O_{p}\cos\left(\frac{2\pi p}{N}\right)+O_{\frac{N}{2}-p}\sin\left(\frac{2\pi p}{N}\right)\right]\\
X_{p}^{\textrm{(jobb)}}&=\frac{1}{\sqrt{2}}\left[E_{p}-O_{p}\cos\left(\frac{2\pi p}{N}\right)-O_{\frac{N}{2}-p}\sin\left(\frac{2\pi p}{N}\right)\right]
\end{align}

A pillangók ezúttal első ránézésre nem öntenek olyan formát, mint az eddigiekben. Bontsuk két részre! Iktassunk be az alábbi műveletet a pillangók előtt:

\begin{align}
O_{p}’&=O_{p}\cos\left(\frac{2\pi p}{N}\right)+O_{\frac{N}{2}-p}\sin\left(\frac{2\pi p}{N}\right)
\end{align}

Most ezt én önkényesen elnevezem “Node-shift” műveletnek és valóban egy csomóval fogom ábrázolni. Ebből a pillangó egyenletei:

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

Az ábra kísértetiesen hasonló a Fourier-transzformációnál kapott pillangókhoz.

A fenti rekurzió bázisfeltétele, az elemszám. Ha a bemeneti vektor hossza egy, akkor a kimenet megegyezik a bemenettel.

Szétválogatjuk a bemenetet páros és páratlan indexek szerint.

Ezek után transzformáljuk a két (ezúttal fele akkora hosszú vektort).

Megcsináljuk a “Node-shift” műveletet.

Végül a Fourier-transzformációból is ismert pillangó műveletet.

Így a teljes függvény:

A program eredménye természetesen azonos:

Lássuk ez iteratívan hogyan néz ki! Rajzoljuk fel pillangókkal a dolgot (például \(N=8\)-ra).

Az első ami feltűnhet az Olvasónak, hogy a bemeneti adatsor picit furán van rendezve. Ez annak köszönhető, hogy páros és páratlan indexekre szedtük szét kétszer egymás után. A sorrend elég furán néz ki a tízes számrendszerben. Próbáljuk meg a kettes számrendszerben!

\begin{align*}
0 & \rightarrow 0 & 000 \rightarrow 000\\
1 & \rightarrow 4 & 001 \rightarrow 100\\
2 & \rightarrow 2 & 010 \rightarrow 010\\
3 & \rightarrow 6 & 011 \rightarrow 110\\
4 & \rightarrow 1 & 100 \rightarrow 001\\
5 & \rightarrow 5 & 101 \rightarrow 101\\
6 & \rightarrow 3 & 110 \rightarrow 011\\
7 & \rightarrow 7 & 111 \rightarrow 111\\
\end{align*}

Ha jobban megnézzük, akkor az adott számhoz hozzá rendeljük a bitek sorrendjének megfordításából kapott számot. Ez az úgynevezett “bit-reversal” rendezés. Ehhez elő kell állítani a megfordítottat. Vegyünk például egy nyolc hosszú szót! Csoportosítsuk kettesével a számjegyeket és csoporton belül cseréljük fel őket. Ezek után csoportosítsuk négyesével a számjegyeket és az előző kettes csoportokat cseréljük fel a négyeseken belül. Tovább menve egy nyolcas csoportot (tehát az egész szót) nézzük és a benne lévő négyes csoportokat cseréljük fel.

Eredményképpen a nyolc hosszú szám fordított sorrendben van felírva. A gondolat az, hogy veszünk egy olyan maszkot, ami minden második bitet hagyja csak meg és annak az inverzét (ami minden első bitet hagyja csak meg). Ezekkel “ÉS” kapcsolatot végzünk, az egyiket balra, a másikat jobbra toljuk és a végén “össze VAGY-oljuk”. Az első maszkolás folyamatát szemlélteti a következő ábra:

Javaslom az Olvasónak, hogy próbálja meg papíron az alábbi kódot a fenti nyolc hosszú szóra.

A “length” paraméterben kell megadni a mi általunk használt szó bitekben mért hosszát. Ugyanis ha nem harminckét bites a szó (például számok nullától nyolcig csak három bitet vesznek igénybe), akkor korrigálni kell, vissza kell tolni. Az általunk rendezni kívánt tömbök hossza szinte kivétel nélkül kettőnek pozitív egész kitevős hatványai, így a tologatás könnyen kivitelezhető. A rendező függvény:

A tesztprogram a rendezésre az alábbi eredményt adja:

Az első feladat tehát, hogy jól rendezzük a bemeneti tömböt. Ezek után a fenti képen látszik (mondjuk úgy, hogy érezzük meg, de belátható matematikailag) hogy annyiszor kell megcsinálni a pillangót, ahányszor az adatsor kettéosztható. Ehhez vegyük a kettes alapú logaritmusát az adatsor hosszának. Hatékonyabban járunk el, ha a függvényünkkel eleve a logaritmust tudatjuk és majd ha kell, akkor kiszámítjuk a hosszt.

Ezek után szükség lesz egy ciklusra, ami a fenti ábrán elkülönülő részeken végig megy. Mint kiderítettük ez \(\log_{2}N\) lépést vesz majd igénybe. Minden ilyen iterációban meg kell csinálni a “Node-shift”-et a szükséges helyeken és a pillangókat.

Ahol \(M\), az éppen aktuális ciklusban a bemeneti adatok száma. Ez rendre \(2,4,8,16\dots\). Most tovább példálózva a fenti képpel, az első iterációban megcsináljuk a négy darab kettes pillangót. A következő ciklusban két darab négyes pillangót kéne elvégezni, de előtte a két darab “Node-shift”-et csináljuk meg. Végül az egy darab nyolcas pillangót csináljuk meg, előtte a “Node-shift”-tel.

Ki kell számolni, hogy az aktuális ciklusban mekkora az M. Ezek után végig kell menni az M-es csoportokon egymás után. Gyakorlatilag minden körben elvégzünk \(\frac{N}{M}\) darab transzformációt.

A \(j\)-edik \(M\)-es csoporton belül a feladat a “Node-shift”, ami \(\frac{M}{4}-1\) iteráció.

A “Node-shift” után meg kell csinálni az \(M\)-es csoporton belül a pillangókat. Ez \(\frac{M}{2}\) ismétlést vesz igénybe.

Ezek egymás után következnek, így a két ciklus együtt \(\left(\frac{M}{4}-1+\frac{M}{2}\right)\sim M\) időt vesz igénybe. A teljes függvény időigénye ezáltal:

\begin{align*}\require{cancel}
\cancel{M}\cdot \frac{N}{\cancel{M}}\cdot\log_{2}N&=N\log_{2}N
\end{align*}

Az iteratív FHT függvény teljes egészében:

Az eredmény természetesen megegyezik a rekurzív függvény eredményével.

Ez az iromány is elég hosszúra sikeredett, de remélem tanulságos volt. A következő cikkben átevezünk egy kicsit a fizikai síkra, elkezdünk mintákat venni a világból.

Github: https://github.com/patrikradvanyi/SEM_DSP