Dekonvolúció és spektrumbecslés

femtokémiai mérések elemzésében

Diplomamunka

Írta: Bánhegyesi Gábor

Alkalmazott matematikus szak

Témavezető:

Dr. Keszei Ernő, egyetemi tanár Fizikai Kémiai Tanszék

Eötvös Loránd Tudományegyetem, Természettudományi Kar

Konzulens:

Pröhle Tamás Valószínűségelméleti és Statisztika Tanszék

Eötvös Loránd Tudományegyetem, Természettudományi Kar

Eötvös Loránd Tudományegyetem Természettudományi Kar

Tartalomjegyzék

1. Bevezető 1 1.1. Konvolúció és dekonvolúció a femtokémiában . . . . . . . . . . . . . . . . . . . . . 1 1.2. A dekonvolúció nehézségei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3. Elterjedt dekonvolúciós módszerek áttekintése . . . . . . . . . . . . . . . . . . . . . 3

1.3.1. Lineáris módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3.2. Nemlineáris módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.4. Impulzus- és lépcső-függvények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2. Gans és Nahman módszere lépcsőfüggvények kezelésére 6

2.1. A Gans-Nahman módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1. Folytonos lépcsőfüggvény esete . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2. Diszkrét lépcsőfüggvények esete . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.3. Egy általánosabb megközelítés . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2. Samulon és Nicolson módszerei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1. Általános alapfeltevések . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2. Samulon eljárása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.3. Nicolson módszere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.4. Ekvivalencia Gans eljárásával . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.5. Következmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3. Spektrumbecslés 14 3.1. Bayes modell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2. A spektrum becslése Kálmán szűrő segítségével . . . . . . . . . . . . . . . . . . . . 15

3.2.1. Algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2.2. Paraméterek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2.3. A torzítás kiküszöbölése . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.3. Függvény előállítása a becsült spektrum alapján . . . . . . . . . . . . . . . . . . . 18 3.3.1. Általános tapasztalatok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.2. A függvény visszaállítása egy konkrét esetben . . . . . . . . . . . . . . . . . 19

3.4. Következmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4. Rejtett periodicitások 23 4.1. Fisher-teszt: ismeretlen frekvenciájú rejtett periodicitások keresése . . . . . . . . . 23 4.2. Kolmogorov–Szmirnov teszt alkalmazása . . . . . . . . . . . . . . . . . . . . . . . . 24 4.3. Nem ekvidisztáns mintavételezés esete . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.4. Periodicitások nem ekvidisztáns minta esetén . . . . . . . . . . . . . . . . . . . . . 26 5. Periodicitások a mért adatsorokban 28

5.1. Ekvidisztáns mintavételezésű adatsorok . . . . . . . . . . . . . . . . . . . . . . . . 28 5.1.1. Periodogramok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.1.2. Fisher-teszt és a Kolmogorov–Szmirnov próba . . . . . . . . . . . . . . . . . 29

5.2. Nem ekvidisztáns mintavételezésű adatsorok . . . . . . . . . . . . . . . . . . . . . . 31 5.2.1. Adatsorok és Lomb-féle periodogramjaik . . . . . . . . . . . . . . . . . . . . 31 5.2.2. Lomb-periodogram és frekvenciák . . . . . . . . . . . . . . . . . . . . . . . . 32

6. Dekonvolúció és spektrumbecslés 34 6.1. Harmonikus összetevők . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 6.2. A valódi és az érzékelt függvény spektruma . . . . . . . . . . . . . . . . . . . . . . 35 6.3. A dekonvolúció végrehajthatósága . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 6.4. Következmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

7. Összefoglaló és további tennivalók 39 7.1. Összefoglaló . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 7.2. További vizsgálati lehetőségek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

A. Forráskódok 42 A.1. Spektrumbecslő algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 A.2. Lomb periodogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.3. Periodogram, Fisher-teszt és Kolmogorov–Szmirnov-próba . . . . . . . . . . . . . . 46

Ábrák jegyzéke

1.1. Az s függvény (vékony vonal), a valódi o négyszögfüggvény és kettejük s ∗ o konvolúciója (vastag vonal). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1. Az f_{3} függvény grafikonja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2. Egy f_{0} diszkrét lépcsőfüggvény grafikonja . . . . . . . . . . . . . . . . . . . . . . . 8 2.3. Egy diszkrét lépcsőfüggvény a Gans által javasolt átalakítás után . . . . . . . . . . 9

3.1. Egy véletlenszerűen generált 128 adatpontból álló adatsor FFT-je (vékony vonal) és a tárgyalt módszerrel becsült spektruma. Az adatsor 0,05 szórásnégyzetű fehér zajjal terhelt, a periodikus komponensek frekvenciái rendre 10,5, 31,3, illetve 47,7 Hz. Az adatsor ekvidisztáns mintavételezésű. . . . . . . . . . . . . . . . . . . . . . 16

3.2. Egy nem ekvidisztánsan mintavételezett adatsor FFT-je (vékony vonal) és a becsült spektrum grafikonja. A folyamatot terhelő fehér zaj szórásnégyzete 0,01, az adatpontok húsz százalékát figyelmen kívül hagytuk, a periodikus komponensek frekvenciái 10,5, 23,3 és 31,7 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.3. Mesterségesen generált adatsor grafikonja. . . . . . . . . . . . . . . . . . . . . . . . 20 3.4. A spetrum becslése az algoritmus egyszeri végrehajtása után. . . . . . . . . . . . . 20 3.5. A valódi (folytonos), illetve a becsült spektrum alapján számolt adatsor (szaggatott). 22

4.1. Egy mesterségesen generált, több harmonikus összetevőt is tartalmazó, fehér zajjal

terhelt adatsor periodogramja és a különbőző k_{α}-khoz tartozó határok. . . . . . . . 25 4.2. Ekvidisztánsan mintavételezett adatsor Lomb-féle periodogramja (folytonos) és FFT-

je (szaggatott). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5.1. Ekvidisztánsan mintavételezett szimulált adatsorok. . . . . . . . . . . . . . . . . . 29 5.2. Az 5.1 ábrán látható adatsorok periodogramjainak logaritmusai. . . . . . . . . . . 30 5.3. Az 5.1 ábra adatsorainak alapján számolt sztenderdizált kumulatív periodogramok.

A szaggatott vonalak a k_{0,01} = 1, 63-hoz tartozó Kolmogorov–Szmirnov-határokat jelölik. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5.4. Valódi mért adatsorok, 2080nm és 558nm nem ekvidisztáns, 800nm ekvidisztáns

mintavételezésű. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.5. Valódi adatsorok Lomb-féle periodogramjainak logaritmusai. A vízszintes vonal a

13,92-nak megfelelő kritikus értéket ábrázolja. . . . . . . . . . . . . . . . . . . . . . 32

6.1. Az inF2.dat adatsorhoz tartozó Lomb-periodogram és az első menetben becsült spektrum logaritmusai. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 6.2. Az eredeti o függvény becsült spektrumának logaritmusa és a detektált mintából és

az impulzusból származtatott becslése. . . . . . . . . . . . . . . . . . . . . . . . . . 36 6.3. A 6.2 ábrán látható spektrum egyik kritikus szakasza, a súlyok logaritmusaival… 36 6.4. Súlyok a 6.2 ábrán látható spektrum egy kritikus szakaszán. . . . . . . . . . . . . . 37 6.5. Az eredeti adatsor (oF2), és a becslés inverziójával előállított függvény. . . . . . . . 38 6.6. Az inF2 adatsor alapján készített becslés első lépésben. . . . . . . . . . . . . . . . 38

1. fejezet

Bevezető

A dolgozatban áttekintést nyújtunk a femtokémiában napjainkban használt dekonvolúciós eljárásokról valamint ismertetünk néhány módszert, melyek a jelfeldolgozás más területein használatosak. Célunk, hogy ezen módszerek segítségével egy olyan dekonvolúciós eljárást dolgozzunk ki, mely bizonyos esetekben pontosabb eredményhez vezet a jelenleg használatos módszereknél.

A fő probléma lépcsőszerű függvények spektrumának becslése, illetve ennek segítségével zajszűrés és a dekonvolúció elvégzése. Ezzel kapcsolatban ismertetjük Gans és Nahman eljárását, melyet már sikerrel alkalmaznak a femtokémiában, és amely lépcsőfüggvények periodikus függvényekké történő átalakítására használható. Ismertetünk ezen kívül egy olyan módszert, mely alkalmazható nem stacionárius és nem ekvidisztánsan mintavételezett adatsorok esetén is. E Bayes becslés különösen akkor alkalmazható jó hatásfokkal, ha valamiféle előzetes információval rendelkezünk a spektrumról. Legtöbbször ez a helyzet, hiszen sok megoldás kémiai-fizikai szempontból eleve nem megengedett. Megmutatjuk, hogy a Gans-Nahman módszerrel átalakított adatasorokra a spektrum ilyetén becslése jobb megoldást ad, mint a nyers adatsorok esetén. Többek között csökkenti a szélhatást, mely tulajdonsága igen értékes.

A Bayes spektrumbecslést a későbbiekben a dekonvolúciós formulákban mintegy a gyors Fourier transzformáció hellyettesítéseként használjuk, így kerülve el a Fourier-transzformáció okozta pontatlanságokat, illetve ezáltal kiterjeszthetők az eljárások olyan esetekre is, melyekre sem a klasszikus FFT, sem általánosításai nem használhatók.

A módszerek alkalmazhatóságához azonban először arról is meg kell győződnünk, hogy az adatsorok valóban tartalmaznak-e periodikus komponenseket. Három módszert ismertetünk, melyek közül kettő ekvidisztáns mintavételezésű adatsorok esetén alkalmazható tesztet nyújt, míg a harmadik módszer nem ekvidisztáns mintavételezésű adatsor esetén vizsgálja a rejtett periodikus komponensek jelenlétét. Látni fogjuk, hogy a periodikus komponensek jelenléte a vizsgált adatsorok mindegyikében magas szignifikancia-szint mellett igazolható.

1.1. Konvolúció és dekonvolúció a femtokémiában

Femtokémiai folyamatokban, ahol a vizsgált jelenségek femtoszekundumokban mérhető időskálákon zajlanak le (így a folyamatok megfigyelése – egyenlőre – nem lehetséges elekronikus eszközökkel), a reakciók észlelése és rögzítése csak fényimpulzusok segítségével lehetséges.

A fénnyel történő leképezés úgy valósul meg, hogy a kémiai reakciót egy fényimpulzus indítja be és egy néhány femtoszekundummal későbbi hasonló impulzus képezi le a reakció aktuális állapotát.

1.1. ábra. Az s függvény (vékony vonal), a valódi o négyszögfüggvény és kettejük s ∗ o konvolúciója (vastag vonal).

Az ilyen típusú mintavételezés azonban több szempontból is problematikus. Egyfelől a leképezett jel zajjal terhes, amihez hozzáadódik a fényjelek felfogására használt érzékelők tökéletlenségéből adódó zavaró hatás, valamint az a tény, hogy az egyes reakciók beindításához hosszabb időtartamú megvilágítás szükséges, mint amennyi a leképezés szempontjából optimális lenne. A reakciót gerjesztő és leképező fényimpulzusok csak távolról hasonlíthatók az optimálisnak tekinthető Dirac δ függvényhez.

Midezeket a torzító hatásokat a kémikusok úgy veszik figyelembe, hogy a megfigyelésként kapott adatsorokat a valódi jel és az azt beindító illetve leképező fény-impulzusok konvolváltjaként jellemzik, azaz ha a valódi folyamat egy o(t) függvény, az észlelést megvalósító fénnyalábot pedig egy s(t) függvény jellemzi, akkor az észleléskor rögzített függvény

i(t) = ∫_{−∞}^{∞} s(τ)o(t − τ)dτ.

A konvolúció műveletét a továbbiakban a ∗ jellel fogjuk jelölni. Ezen jelöléssel a fenti egyenletet az i(t) = s(t) ∗ o(t) alakba írhatjuk.

Így az i(t) és az s(t) függvények ismeretében az o(t) meghatározása dekonvolúcióval lehetséges. A dekonvolúció azonban rendkívül érzékeny a zajokra, és egyéb torzító hatásokra, ami megnehezíti alkalmazását.

1.2. A dekonvolúció nehézségei

Mint ismertes, a konvolúció fenti alakja a frekvenciatérben a következő formát veszi föl

I(ω) = S(ω)O(ω),

ahol az I, S, O függvények az i, s, o függvények Fourier-transzformáltjai. A frekvenciatérben ennek alapján a dekonvolúció matematikailag jól definiálható és – elvileg – könnyen el is végezhető, legalábbis ahol az S(ω) nem tűnik el, hiszen

O(ω) = \frac{I(ω)}{S(ω)}

előállítja a kérdéses o függvény Fourier-transzformáltját, amiből az ineverz transzformáció segítségével o meghatározható. Ehhez azonban mind az O, mind pedig az S függvény pontos ismerete szükséges, amely feltétel általában nem teljesül. A dekonvolúciónak ebben a formájában megjelenő osztás okozza a legnagyobb problémát.

Jelölje N(ω) az I függvény zaj-összetevőjét, és legyen I_{v} = I − N. Ekkor az előző egyenlet felhasználásával O becslésére a

Ô(ω) = \frac{I_{v}(ω) + N(ω)}{S(ω)}

kifejezés adódik. A vizsgált femtokémiai mérésekben a zaj tipikusan fehér zaj, tehát magas frekvenciájú komponeneseket is tartalmaz, az S függvény azonban gyakorlatilag eltűnik magasabb frekvenciák esetén, ami olykor azzal jár, hogy az N/S hányados nagyságrendekkel nagyobb az I_{v}/S-nél. Így az O-ra kapott becslés használhatatlan, hiszen a valódi jelet teljesen elnyomja a magas frekvenciájú zaj.

1.3. Elterjedt dekonvolúciós módszerek áttekintése

A dekonvolúció a fentebb ismertetett probléma miatt ebben a formában gyakorlatilag soha nem kivitelezhető, a viszonylag jó eredményre vezető módszerek többnyire különféle approximációs eljárások.

Ebben a szakaszban rövid áttekintést adunk az irodalomban gyakran említett és a gyakorlatban is gyakran használt dekonvolúciós módszerekről.

1.3.1. Lineáris módszerek

Diszkrét esetben a konvolúciós formulában az integrálás összegzéssé alakul. Ha i_{n}, o_{n},-nel, illetve s_{n}-nel jelöljük a vonatkozó függvények n időpillanathoz tartozó értékeit, akkor a konvolúciós formula:

i_{n} = ∑_{m} s_{n−m}o_{n},

ahol az összegzés az összes olyan m-re végzendő, amelyekhez a megfelelő s_{n−m} nem zérus. Ha i-t és o-t oszlopvektorokkal reprezentáljuk, úgy az egyenlet mátrix-alakja

i = So

lesz, ahol S egy N × N-es mátrix, és ahol N az i és o vektorok mérete, általában az i_{n} adatsor hossza. Az S mátrix elemei rendre az s_{n,m} = s_{n−m} számok. Az S ilyen definíciójából adódik, hogy a mátrix sorai egymás cirkuláris eltoltjai, tehát a mátrix Toeplitz-féle.

A fenti alakból az ismeretlen o vektor mátrix-inverzió segítségével számolható. Azonban a Toeplitz mátrix invertálása numerikus szempontból instabil művelet, mivel a sorok alig különböznek egymástól.

Az egyenlet azonban közelítő iterációval is megoldható. A Van-Cittert-féle eljárás lényege a szukcesszív approximáció. A keresett függvényt a

o^{(k+1)} = o^{(k)} + λ (i − s ∗ o^{(k)})

iterációval közelítjük, ahol λ a relaxációs konstans. Az iteráció az o^{(0)} = i feltétellel indítható. Az iteráció konvergenciája néhány speciális eset kivételével, amikor némi módosításra szorul, biztosítható [15].

A Van-Cittert módszer hibája, hogy az iteráció során a zaj minden lépésben erősödik, ami jelentősen torzíthatja az eredményt.

Lineáris módszernek tekinthető még az előző szakaszban tárgyalt frekvencia-térben adott megoldás is, így a többi lineáris módszerhez hasonlóan ez is komoly problémákkal küzd.

Másrészt azonban az

Ô(ω) = \frac{I_{v}(ω)}{S(ω)} + \frac{N(ω)}{S(ω)}

kifejezésből leolvasható egy igen egyszerű módja a zaj csökkentésének. Ha I_{v} és N frekvenciái jelentősen különböznek egymástól, azaz egyes frekvencia-tartományokban csak egyikük domináns, akkor a zaj egy alkalmasan választott aluláteresztő szűrővel eltávolítható Ô-ból. Ha viszont a két tartománynak van közös része, akkor szűréskor az eredmény is torzulni fog.

1.3.2. Nemlineáris módszerek

A lineáris módszerekkel összehasonlítva sokkal jobb eredmények érhetők el a dekonvolúcióban, ha a jelenséget leíró o függvényről előzetes ismereteink vannak, például tudjuk, hogy exponenciális, vagy trigonometrikus függvény. Ebben az esetben ugyanis a becslés során esetleg csak néhány ismeretlen paramétert kell megbecsülni, és általában a becslés elég jól használható lesz.

Jansson módszerének lényege, hogy a Van-Cittert eljárás λ relaxációs faktorát nem konstansnak választja, hanem egy olyan r súlyfüggvénynek, melynek tulajdonságait fizikai törvényszerűségek alapján határozza meg. Konkrétan az

o^{(k+1)} = o^{(k)} + r (o^{(k)}) (i − s ∗ o^{(k)})

iterációs lépésben szereplő r függvényt a

r (o^{(k)}) = r_{0} (1 − \frac{|o_{max} − o_{min} − 2o^{(k)}|}{o_{max} − o_{min}})

egyenlőséggel definiálja.

A nemlineáris módszerek egy másik családjában az additív korrekció helyett multiplikatívat szokás használni. A következő eljárást elsőként Gold javasolta. Az iterációs lépés Gold szerint az

o^{(k+1)} = o^{(k)} \frac{i}{s ∗ o^{(k)}}

legyen.

Megjegyzendő, hogy a Gold eljárás a Jansson módszer egy speciális esetének tekinthető, méghozzá abban az esetben, ha relaxációs függvénynek az

r (o^{(k)}) = \frac{o^{(k)}}{s ∗ o^{(k)}}

függvényt választjuk. Érdemes még megemlíteni egy Bayes eljárást melyben az iterációs lépésben az s függvény y-tengelyre vonatkozó tükörképe is szerepel. Az iterációt a

o^{(k+1)}(t) = o^{(k)}(t)s(−t) ∗ \frac{i(t)}{s(t) ∗ o^{(k)}(t)}

összefüggés definiálja. Ezen módszer erőssége azon kívül, hogy hatékonyan nyomja el a zajt az, hogy az iteráció folyamán minden lépésben pozitív értékeket kapunk szemben például a Jansson módszerrel, ami a megengedett tartományon kívül negatív korrekciós együtthatókat használ.

1.4. Impulzus- és lépcső-függvények

A címben említett két függvénytípus gyakran előfordul a kémiában és a fizikában. A dolgozat által vizsgált jelenségek vizsgálatában is mindkét függvénytípussal találkozunk.

Lépcsőfüggvények jellemeznek tulajdonképpen minden kémiai és fizikai folyamatot. Egy anyagi reakciót többnyire adott mennyiségű energia közlésével lehet beindítani, majd kellő idő elmúltával (a reakció végére) az anyag egy másik (stabil) energia állapotba jut. A folyamat során az anyag válaszfüggvényét egy lépcsőfüggvény írja le. A dolgozatban általában ezt jelöljük o-val.

Az energiát a reakcióban résztvevő anyagoknak a mi esetünkben egy fényimpulzus adja. Mivel azonban technikailag nem megoldható olyan fényimpulzust előállítani, mely a Dirac-δ-hoz közelítene, azaz az energia nem adható át a reakciónak egy adott végtelen kis kiterjedésű időpillanatban, az energiát hordozó fényimpulzus torzítja magát az o függvényt is. A reakció leképezését egy hasonló függvény oldja meg.

Gans és Nahman módszere

lépcsőfüggvények kezelésére

Gans és Nahman [3] egy olyan módszert javasol, mely lépcsőfüggvényeket alakít át véges tartójú függvényekké, alkalmassá téve őket FFT számítására. A módszer előnye, hogy könnyen implementálható és jól kezelhetők a bejövő új frekvenciák. Az irodalomban két másik módszer található ugyanerre a problémára. Samulon egy lineáris függvénnyel vágja le a lépcsőfüggvény „végét,” Nicolson pedig deriválással éri el, hogy a függvény impulzus-jellegűvé váljon.

A fejezetben ismertetjük azt az eredményt, mely szerint a három – látszólag különböző – módszer ugyan ahhoz a numerikus eljáráshoz vezet. Gans módszerét a könnyű implementáció miatt választottuk. A leírásban a jelöléseket Gans és Nahman eredeti cikkének megfelelően használjuk.

2.1. A Gans-Nahman módszer

2.1.1. Folytonos lépcsőfüggvény esete

Tekintsünk egy folytonos lépcsőfüggvényt, amely rendelkezik azzal a tulajdonsággal, hogy egy megfelelően választott t_{1} indextől kezdve a függvény értékeinek különbsége egy előre meghatározott ε állandónál kisebb. Ha ε kellően kicsiny, a függvény tekinthető – fizikai értelemben – olyannak, ami az index után már nem változtatja értékét, azaz megadható az

f_{0}(t) = \{ \begin{matrix}0, & t < 0, \\ f_{1}(t), & 0 < t < t_{1}, \\ K, & t > t_{1} \end{matrix}

hozzárendeléssel, ahol f_{1} valamilyen tetszőleges függvény. Ezt a függvényt alakítjuk olyan alakúvá, hogy egy periodikus függvény egy periódusának legyen tekinthető és így alkalmassá váljon FFTszámításokra.

Legyen tehát

f_{3}(t) = f_{0}(t) − f_{0}(t − t_{1}).^{1}

1Megjegyezzük, hogy az eredeti cikkben az összegzés második tényezője még be van szorozva a (t_{1}, 2t_{1}) intervallum karakterisztikus függvényével, ez azonban fölösleges, hiszen a definíció szerint az f_{0} függvény nullánál kisebb értékekre eltűnik.

FEJEZET 2. GANS ÉS NAHMAN MÓDSZERE LÉPCSŐFÜGGVÉNYEK KEZELÉSÉRE 7

2.1. ábra. Az f_{3} függvény grafikonja

Az ilyen módon nyert f_{3} függvény Fourier-transzformáltja az

F_{3}(ω) = ∫_{−∞}^{∞} f_{3}(t)e^{−jωt}dt = ∫_{0}^{2t_{1}}f_{3}(t)e^{−jωd}dt = F_{0}(ω) (1 − e^{−jωt_{1}})

alakban áll elő, ahol j a komplex egységet jelöli. Ezt egyszerűbb formába írva

F_{3}(ω) = F_{0}(ω)F_{4}(ω).

A továbbiakban F_{4}-et vizsgáljuk. Amikor cos ωt_{1} = −1, a körfrekvenciák ω_{n} = nπ/t_{1}, ahol n páratlan egész szám, így sin ω_{n}t_{1} eltűnik. Következésképpen

F_{3}(ω_{n}) = 2F_{0}(ω_{n}), n = ±1, ±3, ±5, …

Másrészről, amikor cos ωt_{1} = 1, a körfrekvenciák ω_{m} = mπ/t_{1} alakban állnak elő, ahol m páros szám, így a képzetes rész ismét eltűnik. Páros egészekre tehát a

F_{3}(ω_{n}) = 0, n = ±0, ±2, ±4, …

formula áll.

Látható tehát, hogy a diszkrét ω_{n} frekvenciákra az F_{3} függvény abszolút értéke éppen kétszerese a keresett F_{0}(ω) függvény abszolút értékének:

|F_{3}(ω)| = |F_{0}(ω)||F_{4}(ω)| = 2|F_{0}(ω)|| sin ωt_{1}/2|.

2.1.2. Diszkrét lépcsőfüggvények esete

Diszkrét lépcsőfüggvény átalakítása hasonlóan megy, mint a folytonos lépcsőfüggvények átalakítása. Ebben az esetben a periodikussá alakított függvény

f_{3}(nT) = f_{0}(nT) − f_{0}(nT − nt),

2.2. ábra. Egy f_{0} diszkrét lépcsőfüggvény grafikonja

ahol T a lépésköz, t pedig az az index, amitől kezdve az eredeti f_{0} lépcsőfüggvény konstansnak tekinthető. Ennek a függvénynek a diszkrét Fourier transzformáltja

F_{3}(kω) = F_{0}(kω) − (e^{−jωkNT} F_{0}(kω)) ,

ahol NT = t, és k = 0, 1, …, 2N − 1. Az egyenletet némileg átalakítva kapjuk, hogy

F_{3}(ωk) = ∑_{n=0}^{2N−1} f_{0}(nT)e^{−jωT nk} (1 − e^{−jωT kN}) .

A szummában szereplő (1 − e^{−jωT kN}) kifejezés megfelel a folytonos esetben bevezetett F_{4}(ω) függvénynek. Jelöljük most ezt a függvényt is F_{4}-gyel. A függvény abszolút értékét

|F_{4}(ωk)| = |1 − e^{−jωT kN} | = 2| sin πk/2|

alapján határozhatjuk meg. Ebből látszik, hogy az F_{4} függvény éppen úgy súlyozza a spekrumot a diszkrét esetben, mint ahogy tette azt a folytonosban.

2.1.3. Egy általánosabb megközelítés

Gans módszerével a folytonos esetben arra jutottunk, hogy a kiszámítható F_{3} függvény a keresett F_{0} függvény és egy F_{4} szorzataként áll elő. Megpróbáljuk ezt az összefüggést felhasználni arra, hogy a spektrumot ne csak a Gans és Nahman által javasolt esetekben tudjuk visszaállítani.

Tekintsük tehát a

F_{3}(ω) = F_{0}(ω) (1 − e^{iωt_{1}}) = F_{0}(ω)F_{4}(ω)

2.3. ábra. Egy diszkrét lépcsőfüggvény a Gans által javasolt átalakítás után

egyenletet, ahol t_{1} a szokásos indexet jelöli. Az F_{4} függvény olyan, hogy csak két valós értéket vehet föl, nevezetesen

F_{4}(ω) = \{ \begin{matrix}2, & ω = \frac{(2k+1)π}{t_{1}} , \\ 0, & ω = \frac{2kπ}{t_{1}} . \end{matrix}

Gans csak azzal az első esettel foglalkozik, amikor a keresett F_{0} értéke éppen fele az F_{3}-nak.

Az F_{3} = F_{0}F_{4} komplex egyenletet e^{−ix} = cos x − i sin x összefüggés felhasználásával átalakítva kapjuk, hogy

F_{3}(ω) = F_{0}(ω) [1 − (cos ωt_{1} − i sin ωt_{1})] .

Ebből az egyenletből a minden x, y komplex számra érvényes

Rxy = RxRy − IxIy,
Ixy = RxIy + RyIx

azonosság felhasználásával az

RF_{3}(ω) = (1 − cos ωt_{1})RF_{0}(ω) − sin ωt_{1}IF_{0}(ω),
IF_{3}(ω) = (1 − cos ωt_{1})IF_{0}(ω) + sin ωt_{1}RF_{0}(ω)

egyenletrendszer adódik. Ezt a könnyebb átláthatóság kedvéért mátrixos alakban írjuk föl:

( \begin{matrix}1 − cos(ωt_{1}) & − sin(ωt_{1}) \\ sin(ωt_{1}) & 1 − cos(ωt_{1}) \end{matrix} ) ( \begin{matrix}RF_{0}(ω) \\ IF_{0}(ω) \end{matrix} ) = ( \begin{matrix}RF_{3}(ω) \\ IF_{3}(ω) \end{matrix} ) .

Az egyenletet balról a négyzetes mátrix inverzével szorozva éppen F_{0} egy előállítását kapjuk, ahol az inverz létezik. Az inverz létezésének feltétele, hogy a szóban forgó mátrix determinánsa ne tűnjön el. Némi számolással adódik, hogy az inverz az ω = \frac{2kπ}{t_{1}} helyeken nem létezik, azaz ott nincs inverz, ahol az eredeti módszerben a szorzó zérus. Ha tehát ω nem áll elő ebben az alakban, akkor

( \begin{matrix}RF_{0}(ω) \\ IF_{0}(ω) \end{matrix} ) = \frac{1}{2} ( \begin{matrix}1 & \frac{sin(ωt_{1})}{1−cos(ωt_{1})} \\ \frac{− sin(ωt_{1})}{1−cos(ωt_{1})} & 1 \end{matrix} ) ( \begin{matrix}RF_{3}(ω) \\ IF_{3}(ω) \end{matrix} ) .

Ha azokban az esetekben, amikor a mellékátló elemei indefinitté válnak, a mellékátlót nullának tekintjük, éppen Gans eredményét kapjuk. Azonban ez a módszer nem csak ezekben az esetekben teszi lehetővé a spektrum előállítását.

A megközelítés szépséghibája, hogy nem alkalmazható automatikusan, mivel a mellékátló elemei nem konvergálnak a kritikus esetekben, konkrétan:

lim_{ω→ \frac{2lπ}{t_{1}}} \frac{± sin ωt_{1}}{1 − cos ωt_{1}} = lim_{ω→ \frac{2lπ}{t_{1}}} \frac{± cos ωt_{1}}{sin ωt_{1}} = ±∞,

ahol l tetszőleges egész szám.

2.2. Samulon és Nicolson módszerei

Az alábbiakban két további módszert ismertetünk, melyek lépcsőfüggvények Fourier transzformáltjának meghatározására használhatók, és megmutatjuk, hogy a két módszer numerikus szempontból mind egymással, mind pedig Gans módszerével egyenértékű.

2.2.1. Általános alapfeltevések

Legyen a továbbiakban f(t) egy lépcsőfüggvény, ami csak egy T időpontig mutat változásokat, azaz f(t) = f(T), ha t ≥ T. Feltesszük, hogy az f függvényt Δt időközönként mintavételezzük, összesen N mintát véve. Δt-t úgy választjuk, hogy NΔt = T teljesüljön. f-et periodikus függvénynek tekintjük T periódussal, így a spektrum diszkrét lesz. A harmonikusok frekvenciái a következők lehetnek

ω_{k} = kω_{0} = k \frac{1}{T} , k = 0, …, N − 1.

Legyen a továbbiakban f_{i} = f(iΔt), illetve F_{k} = F (kω_{0}). Ezek között a DFT

F_{k} = ∑_{i=0}^{N−1} f_{i}e^{\frac{−jik2π}{N}} , k = 0, …, N − 1

illetve az inverz DFT

f_{i} = \frac{1}{N} ∑_{k=0}^{N−1} F_{k}e^{\frac{jik2π}{N}} , i = 0, …, N − 1

teremt kapcsolatot, ahol j a komplex egységet jelöli.

Mivel az alábbi módszerekben lépcsőfüggvények impulzus-szerű függvényekké való átalakításáról van szó, így a zéró frekvencia esete nem számolható, tehát a számítások csak a k = 1, …, N − 1 esetekre vonatkoznak.

2.2.2. Samulon eljárása

Samulon [5] módszerének lényege abban áll, hogy az eredeti adatsor deriváltját tekintjük, az ugyanis már impulzus-jellegű függvény lesz. Ebben a megközelítésben szükség van az egymás utáni adatpontok amplítúdókülönbségeinek kiszámítására. A transzformációt az

F (ω) = \frac{Δt}{1 − e^{−j2πωΔt}} ∑_{i=0}^{∞} (f_{i} − f_{i−1}) e^{−ji2πωΔt}

összefüggés adja meg. Az ω helyére az

ω_{k} = kω_{0} \frac{1}{T} = k \frac{1}{NΔt}, k = 1, …, N − 1

értékeket helyettesítve az egyenlet az

F_{k}^{(S)} = \frac{Δt}{1 − e^{− \frac{jk2π}{N}}} ∑_{i=0}^{N−1} (f_{i} − f_{i−1}) e^{− \frac{jik2π}{N}}

alakba írható át, ahol (S) Samulon nevére utal. Az FFT formulákkal való konzisztencia kedvéért a Δt faktort eltávolítjuk. Ebben a formában a jobb oldalon szerplő összegzést két komponensre bontva az

F_{k}^{(S)} = \frac{1}{1 − e^{− \frac{jk2π}{N}}} \{ ∑_{i=0}^{N−1} f_{i}e^{− \frac{jik2π}{N}} − ∑_{i=0}^{N−1} f_{i−1}e^{− \frac{jik2π}{N}} \}

kifejezést kapjuk. Legyen S_{k} = ∑_{i=0}^{N−1} f_{i}e^{− \frac{jik2π}{N}} . Ekkor

F_{k}^{(S)} = \frac{1}{1 − e^{− \frac{jk2π}{N}}} \{S_{k} − ∑_{i=0}^{N−1} f_{i−1}e^{− \frac{jik2π}{N}} \} .

A szummából e^{− \frac{jk2π}{N}} -et kiemelve kapjuk, hogy

F_{k}^{(S)} = \frac{1}{1 − e^{− \frac{jk2π}{N}}} \{S_{k} − e^{− \frac{jk2π}{N}} ∑_{l=0}^{N−2} f_{l}e^{− \frac{jlk2π}{N}} \} ,

amiből viszont

F_{k}^{(S)} = \frac{1}{1 − e^{− \frac{jk2π}{N}}} \{S_{k} − e^{− \frac{jk2π}{N}} (S_{k} − f_{N−1}e^{− \frac{j(N−1)k2π}{N}} )\} ,

következik. Mivel e^{−jk2π} = 1,

F_{k}^{(S)} = S_{k} + \frac{f_{N−1}}{1 − e^{− \frac{jk2π}{N}}} , k = 1, …, N − 1.  (2.1)

2.2.3. Nicolson módszere

Tekintsük az alábbi függvényt

r(iΔt) = i\frac{f_{N−1}}{N} , i = 0, …, N − 1.

Nicolson [6] megmutatta, hogy ezt a függvényt kivonva az eredeti lépcsőfüggvényből egy impulzusjellegű függvény adódik, melynek harmonikusai között az r függvény harmonikusai nem jelennek meg, ha a keletkezett függvényt T = NΔt-től zérónak tekintjük. Így

F_{k}^{(N)} = Δt ∑_{i=0}^{N−1} (f_{i} − r_{i}) e^{\frac{−jik2π}{N}} .

Samulon esetéhez hasonlóan itt is eltekintünk a Δt faktortól, és a két szummát ismét külön számoljuk:

F_{k}^{(N)} = Δt ∑_{i=0}^{N−1} f_{i}e^{\frac{−jik2π}{N}} − \frac{f_{N−1}}{N} ∑_{i=0}^{N−1} ie^{\frac{−jik2π}{N}} = S_{k} − \frac{f_{N−1}}{N} R_{k},

ahol

R_{k} = ∑_{i=0}^{N−1} ie^{\frac{−jik2π}{N}} .

Felhasználva, hogy R_{k} − e^{− \frac{jk2π}{N}} R_{k} = ∑_{l=1}^{N} e^{− \frac{jlk2π}{N}} − Ne^{− \frac{jNk2π}{N}} , valamint azt, hogy ∑_{i=1}^{N} e^{\frac{jik2π}{N}} = 0,

illetve e^{−jk2π} = 1, következik, hogy

R_{k} = \frac{−N}{1 − e^{\frac{−jk2π}{N}}} ,

és így

F_{k}^{(N)} = S_{k} + \frac{f_{N−1}}{1 − e^{− \frac{jk2π}{N}}} , k = 1, …, N − 1.  (2.2)

Fontos megjegyezni, hogy az állítás csak akkor teljesül, ha az eredeti lépcsőfüggvény nulláról indul, azaz f_{0} = 0. Ha ez nem teljesül, akkor F_{k}^{(N)} előállításánál egy

\frac{−f_{0}}{1 − e^{− \frac{jk2π}{N}}}

hibatag jön be, mint azt Laguna [7] kimutatta.

2.2.4. Ekvivalencia Gans eljárásával

Az egységesség kedvéért jelölje most a Gans cikkében f_{3}-mal jelölt függvényt g, azaz g(t) = f(t) − f(t − T). Ekkor tehát

g_{i} = g(iΔt) = \{ \begin{matrix}f_{i}, & i = 0, …, N − 1 \\ f_{N−1} − f_{i−N}, & i = N, …, 2N − 1 \\ 0, & egyébként \end{matrix} .

Ekkor, Gans alapján, legyen

G(ω) = F (ω) (1 − e^{−j2πωT} ) .

Ha ω_{1} = \frac{1}{2NΔt}, akkor

G(mω_{1}) = G_{m} = ∑_{i=0}^{2N−1} g_{i}e^{− \frac{jim2π}{2N}} ,

amit tovább alakítva

G_{m} = ∑_{i=0}^{N−1} f_{i}e^{− \frac{jim2π}{2N}} + ∑_{i=N}^{2N−1}(f_{N−1} − f_{i−N})e^{− \frac{jim2π}{2N}} , m = 1, …, 2N − 1

adódik. Mivel a ω_{k} = kω_{0} = k \frac{1}{2NΔt} harmonikus frekvenciákat keressük, a páros m-ekkel számolunk tovább:

G(2kω_{1}) = G(kω_{0}) = G_{k} = ∑_{i=0}^{N−1} f_{i}e^{− \frac{jim2π}{N}} + ∑_{i=N}^{2N−1}(f_{N−1} − f_{i−N})e^{− \frac{jim2π}{N}} , k = 1, …, N − 1.

A második összegzést szétválasztva

G_{k} = S_{k} + ∑_{i=N}^{2N−1} f_{N−1}e^{− \frac{jim2π}{N}} − ∑_{i=N}^{2N−1} f_{i−N}e^{− \frac{jim2π}{N}} ,

ahonnan

G_{k} = (S_{k} + \frac{f_{N−1}}{1 − e^{− \frac{jk2π}{N}}} ) (1 − e^{−jk2π}) .

Ekkor G_{k} helyére annak eredeti definícióját helyettesítve:

F_{k} (1 − e^{−jk2π}) = (S_{k} + \frac{f_{N−1}}{1 − e^{− \frac{jk2π}{N}}} ) (1 − e^{−jk2π}) .

Mivel k egész értékeire a (1 − e^{−jk2π}) kifejezés eltűnik, az egyenlet valódi információt hordozó része vizsgálhatatlanná válik. Ha azonban feltesszük, hogy k egészhez közeli érték, akkor már joggal írható, hogy

F (ω) = F_{ω/ω_{0}} = S_{ω/ω_{0}} + \frac{f_{N−1}}{1 − e^{− \frac{jω/ω_{0}2π}{N}}} .

Ezt az eredményt k egész értékeire kiterjesztve

F_{k}^{(G)} = S_{k} + \frac{f_{N−1}}{1 − e^{− \frac{jk2π}{N}}} , k = 1, …, N − 1.  (2.3)

2.2.5. Következmények

Megmutattuk, hogy a lépcsőfüggvények átalakítására használt három módszer ugyanahhoz a numerikus módszerhez vezet, hiszen (2.1), (2.2), és (2.3) egyenletek ekvivalensek egymással, következésképpen mindhárom módszer ugyanolyan érzékeny a zajra és a Nyquist frekvenciánál magasabb frekvenciák esetében mutatott torzításuk is megegyezik.

Samulon eljárásából az is látható, hogy a kiindulási függvény alapvonala egyik esetben sem kell, hogy a nulla legyen.

Spektrumbecslés

Spektrum becslésénél általános feltételnek számít, hogy a megfigyelt folyamat legyen stacionárius, másrészt, hogy a mintavételi időpontok ekvidisztánsan helyezkedjenek el. A stacionaritás legtöbbször valamilyen ablakolási eljárással érhető el, míg az ekvidisztáns időosztást további időpontok beszúrásával és ott a folyamat becslésével lehet megoldani. Yuan Qi és mások [8] megoldást javasolnak mindkét problémára.

A javasolt megoldás egy nem-stacionárius Kálmán-szűrő segítségével egyszerre becsli az összes koefficienst, míg Bretthorst általánosított DFT eljárása és a Lomb-Scargle periodogram egyszerre csak egy frekvenciát képes becsülni, így ezeket a módszereket sokszor kell egy adatsorra alkalmazni, hogy jó eredményre jussunk. Az fejezetben bemutatott eljárás további előnye, hogy teljesen független a mintavételezési időpontok között eltelt időtől. A szerzők azt is kimutatják, hogy a módszer pontosabb eredményeket ad, mint a Lomb-Scargle eljárás és más klasszikus módszerek.

3.1. Bayes modell

Ebben a szakaszban ismertetjük a cikkben [8] javasolt modellt, melyre a becslést végző algoritmus épül. A modell nem tételezi föl a jel stacionaritását, szemben a klasszikus modellekkel. Az eljárás ekvidisztánsan és nem ekvidisztánsan mintavételezett adatsor esetén egyaránt működik.

Tegyük fel, hogy a folyamatot a t_{i} időpontokban figyeljük meg (i = 1, …, N), és az előzetes feltételezéseknek megfelelően t_{i} − t_{i−1} nem feltétlenül állandó. A t_{i} időpontokban megfigyelt értékeket jelöljük x_{i}-vel. Ilyen módon kapunk egy x és egy t vektort, mely a megfigyelés minden információját tartalmazza.

Az x_{i}-ket a következő összefüggéssel modellezzük:

x_{i} = a_{i0} + ∑_{j=1}^{M} a_{ij} sin(2πf_{j}t_{i}) + b_{ij} cos(2πf_{j}t_{i}) + v_{i}, i = 1, …, N,

ahol v_{i} a zajt jelzi. Az M és f_{j} változók tetszőlegesen választhatók. A szerzők a cikkben az f_{j} frekvenciákat ekvidisztánsnak választják. A trigonometrikus függvények a_{ij}, illetve b_{ij} együtthatóit valósaknak tételezzük föl. Megjegyzendő, hogy nem-stacionárius folyamat esetén mind ezek az együtthatók, mind pedig a zajtényező függ a mintavételezés időpontjaitól.

Mivel mindkét trigonometrikus alapfüggvényt felhasználjuk, ezáltal elérhető, hogy a folyamat amplitúdója is időfüggő legyen és változzon a fázis is, nevezetesen:

A = \sqrt{a_{ij}^{2} + b_{ij}^{2}}, φ = arctan (\frac{b_{ij}}{a_{ij}} )

a j-edik frekvencia-sáv esetén. Az egyszerűség kedvéért bevezetjük az alábbi jelöléseket:

s_{i} = [a_{i0}, …, a_{iM}, b_{i1}, …, b_{iM}]^{T}  (3.1)
c_{i} = [1, sin(2πf_{1}t_{i}), …, sin(2πf_{M}t_{i}), cos(2πf_{1}t_{i}), …, cos(2πf_{M}t_{i})] .  (3.2)

E jelölések segítségével, és feltéve, hogy az x_{1}, …, x_{N} megfigyelések egy Markov-lánc által generált idősorból származnak, a modell az

s_{i} = s_{i−1} + w_{i}  (3.3)
x_{i} = c_{i}s_{i} + v_{i}  (3.4)

alakba írható, ahol w_{i} a folyamat zaja a t_{i} időpontban, v_{i} pedig ugyanebben az időpontban a megfigyelés zaja. A számítások egyszerűsítése végett a szerzők felteszik, hogy a zaj normális eloszlású.

E modell esetén a folyamat rejtett állapotainak együttes eloszlására

p(s_{1:N}, x_{1:N}) = p(s_{1})p(x_{1}|s_{1}) ∏_{i=2}^{N} p(s_{i}|s_{i−1})p(x_{i}|s_{i}),

teljesül, ahol s_{1:N} = [s_{1}, …, s_{N}]^{T} , és s_{1:N} = [s_{1}, …, s_{N}]^{T} .

A p(s_{i}|x_{1:i}) valószínűségek a

p(s_{i}|x_{1:i−1}) = ∫_{s_{i−1}} p(s_{i}|s_{i−1})p(s_{i−1}|x_{1:i−1})  (3.5)
p(s_{i}|x_{1:i}) = \frac{p(x_{i}|s_{i})p(s_{i}|x_{1:i−1})}{p(x_{i}|x_{1:i−1})}  (3.6)

formulákkal becsülhetők. Így a spektrum becslése a p(s_{i}|x_{1:i}) átlagainak összegeként áll elő.

3.2. A spektrum becslése Kálmán szűrő segítségével

3.2.1. Algoritmus

Ha a lineáris Gauss modellt (3.3 és 3.4 egyenletek) használjuk azzal a feltevéssel, hogy a modellben szereplő zajok normális fehér zajok (GWN), azaz

w_{i} ∼ N(0, Γ_{i})
v_{i} ∼ N(0, σ),

akkor a keresett p(s_{i}|x_{1:i−1}) is normális eloszlású, és Kálmán szűrővel hatékonyan számíthatók ezek a valószínűségek. A nem ekvidisztáns mintavétel okozta probléma kezelésére bevezetjük a

Γ_{i} = Z(t_{i} − t_{i−1})

mátrixokat, ahol Z egy előre definiált konstans mátrix, és Γ_{0} = 0.

Jelölje a p(s_{i}|x_{1:i}) kovariancia-mátrixát V_{i}, a várható értékét pedig m_{i}. Ezekre ekkor

m_{i} = m_{i−1} + K_{i}(x_{i} − c_{i}m_{i−1})  (3.7)
V_{i} = (I − K_{i}c_{i})P_{i−1},  (3.8)

ahol

P_{i−1} = V_{i−1} + Γ_{i−1}
K_{i} = P_{i−1}c_{i}^{T} (c_{i}P_{i−1}c_{i}^{T} + σ)^{−1} .

3.1. ábra. Egy véletlenszerűen generált 128 adatpontból álló adatsor FFT-je (vékony vonal) és a tárgyalt módszerrel becsült spektruma. Az adatsor 0,05 szórásnégyzetű fehér zajjal terhelt, a periodikus komponensek frekvenciái rendre 10,5, 31,3, illetve 47,7 Hz. Az adatsor ekvidisztáns mintavételezésű.

Mivel egy nem stacionárius Kálmán szűrőről van szó, a c_{i} és Γ_{i−1} vektorok időfüggőek. A rekurzióhoz a kezdőértékeket az alábbiaknak választhatjuk:

m_{1} = m_{0} + K_{1}(x_{1} − c_{1}m_{0})
V_{1} = (I − K_{1}c_{1})V_{0}
K_{1} = V_{0}c_{1}^{T} (c_{1}V_{0}c_{1}^{T} + σ)^{−1} ,

ahol m_{0} és V_{0} előre megadott paraméterek, melyek a p(s_{0}) eloszlását jellemzik.

3.2.2. Paraméterek

Mivel az algoritmus bayes-i alapokon nyugszik lehetőség van rá, hogy a priori ismereteinket felhasználhassuk. Ha semmi információval nem rendelkezünk, akkor a Z mátrixot Z = zI alakúnak választhatjuk, ahol a z valós értékű amplitúdó az egyes frekvenciákhoz tartozó amplitúdók szórását befolyásolja. Használható az információt nem hordozó a priori p(s_{0}) eloszlás is. Ehhez V_{0}-t egy megfelelően skálázott egységmátrixnak, m_{0}-t pedig nullvektornak kell választani.

Ha az adatsor csak néhány ismert frekvenciából áll, a becslés végezhető kizárólag ezekre a frekvenciákra is. Ekkor azonban nincs rá mód, hogy az esetleges nem ismert frekvenciákhoz tartozó amplitúdókat is becsüljük. Célravezetőbb lehet ezért a V_{0} mátrixot diagonálisnak választani úgy, hogy a nem ismert frekvenciákhoz tartozó értékeket sem nullának, hanem valamilyen kis értéknek definiáljuk. Ilyen módon nem zárjuk ki a többi frekvencia megjelenését sem.

3.2. ábra. Egy nem ekvidisztánsan mintavételezett adatsor FFT-je (vékony vonal) és a becsült spektrum grafikonja. A folyamatot terhelő fehér zaj szórásnégyzete 0,01, az adatpontok húsz százalékát figyelmen kívül hagytuk, a periodikus komponensek frekvenciái 10,5, 23,3 és 31,7 Hz.

3.2.3. A torzítás kiküszöbölése

Bretthorst [11] kimutatta, hogy a DFT egy olyan általánosításával, amely lehetővé teszi alkalmazását nem ekvidisztáns mintavételezésből származó adatsorokra is, nagyobb effektív sávszélesség érhető el. Az előbbiekben ismertetett módszerre hasonló állítás igaz: az a kritikus frekvencia, ami fölött jelntős torzítással kell számolni nem ekvidisztánsan mintavételezett adatok esetén akár nagyon nagy is lehet.

Ekvidisztánsan mintavételezett adatsor esetén

t = nΔt, n = 1, …, N
f_{s} = \frac{1}{Δt}
f_{Ny} = \frac{f_{s}}{2} = \frac{1}{2Δt},

ahol Δt az egymást követő mintavételek közötti időtartam, f_{s} a mintavételezési frekvencia, f_{Ny} pedig a Nyquist frekvencia. Jól ismert, hogy ebben az esetben a torzítások a mintavételezési frekvencia többszöröseinél jelentkeznek.

Ezzel szemben jelölje nem ekvidisztáns mintavételezés esetén Δt az összes mintavételi intervallum legnagyobb közös osztóját, azaz

Δt = max {t : t|(t_{i} − t_{i−1}), ∀ i} .

Ekkor

t_{n} = k_{n}Δt, n = 1, …, N,

ahol k_{n} megfelelően választott természetes szám.

Ebben az esetben a Nyquist frekvencia f′_{Ny} az egyenletesen mintavételezett esettel összehasonlítva lényegesen nagyobb is lehet, hiszen az

f′_{Ny} = \frac{1}{2Δt}

formulában a Δt nem nagyobb a legrövidebb mintavételezési intervallumnál.

3.3. Függvény előállítása a becsült spektrum alapján

Ha a spektrum becslése helyes akkor az eredeti függvény könnyen visszaállítható. Tapasztalataink szerint azonban az eljárás nem becsli kellő potossággal a spektrumot, legalábbis abban az esetben, amikor nem rendelkezünk semmiféle előzetes információval.

Megmutatjuk, hogy ebben az esetben az algoritmus kétszeri alkalmazásával érhető el jó eredmény úgy, hogy az első alkalmazás eredményeit a priori ismeretként használjuk fel a második alkalmazáskor.

A módszer használhatóságát úgy demonstráljuk, hogy egy adatsorra elvégezzük a spektrumbecslést, majd a kapott spektrum segítségével próbáljuk visszaállítani az eredeti függvényt.

3.3.1. Általános tapasztalatok

Az algoritmust számos szintetikus és valódi adatsorra alkalmazva az a tapasztalati tény állapítható meg, hogy a tárgyalt eljárás adatsorra történő egyszeri alkalmazása alapvetően helytelenül becsli meg a spektrumban az egyes frekvenciákhoz tartozó súlyokat, amennyiben a becsült értékek olykor nagyságrendi eltérést mutatnak a valódi értékekhez képest. Ezek a súlyok azonban szinte kivétel nélkül a helyes frekvenciákhoz tartoznak és élesen elkülönülnek a spektrum többi részétől, ugyanis az adatsor periodikus összetevőihez tartozó frekvenciák súlya legalább két, de gyakran három, vagy négy nagyságrenddel nagyobb, mint a többi frekvenciáé.

Technikai szempontpól az okoz problémát, hogy az algorimust végrehajtó program – nyilvánvalóan – csak diszkrét frekvenciákra becsli meg a súlyokat. Ebben az esetben pedig könnyen előfordulhat (sőt ez az általános eset), hogy egy valódi frekvencia két vizsgált frekvencia közé esik. Ekkor a becslés pontossága attól függ, hogy a valódi frekvencia a vizsgált frekvenciáktól milyen messze helyezkedik el. Legrosszabb a helyzet akkor, amikor a felezőpont közelében van. Ilyen esetekben a szomszédos pontokban a becsült súly közel azonos. A valódi frekvencia helyére a szomszédos frekvenciákhoz tartozó csúcsok arányából lehet következtetni. Elképzelhető természetesen az is, hogy szomszédos frekvenciákhoz valódi csúcsok tartoznak. Ebben az esetben a frekvencia ilyetén kiválasztása az eljárás hibáját növelheti, ezért célszerű ezeket az eseteket réyzletesebben megvizsgálni.

Természetes gondolat, hogy a becslés precizitását a vizsgált frekvenciák számának növelésével, azaz (lévén, hogy adott időintervallumokról van szó) a frekvencia-felosztás finomításával próbáljuk elérni. Mivel azonban itt egy, a vizsgált frekvenciák számában másodfokú algoritmusról van szó (ami a tárigényt illeti), a számítási idő gyorsan növekszik.

Az első menetben így nagy hatékonysággal meghatározhatók a lényeges frekvenciák, feltéve, hogy azok nincsenek egymáshoz „túl közel.”

Célszerű ekkor az algoritmust újra alkalmazni, méghozzá úgy, hogy csak az előbb kiválasztott frekvenciákhoz tartozó súlyok becslését végezze el.^{1} Az ekkor kapott eredmények kellően precízek ahhoz, hogy a szóban forgó függvény visszaállítása az ekkor kapott spektrum alapján már sikeres legyen.

1Természetesen a folyamat tetszőleges sok lépésre bontható szét, ekkor minden lépésben egy intervallumot vá-

Érdemes azt is megjegyezni, hogy az algoritmus újbóli alkalmazásának időigénye többnyire nem jelentős, hiszen míg az első futtatáskor akár frekvenciák százaihoz tartozó súlyokat kell becsülni, addig a második alkalommal ezt már csak a kiválasztott frekvenciák esetén kell megtenni, melyek száma legtöbbször legfeljebb tizes nagyságrendű.

3.3.2. A függvény visszaállítása egy konkrét esetben

A visszaállíthatóság demonstrálásához először generáltunk egy adasort, mely három periodikus komponens és zaj összegeként áll elő. Feltesszük, hogy n számú adattal rendelkezünk, melyek a t_{i} időpontokhoz tartoznak. Ezeket az időpontokat, valamint a frekvenciákat ebben az esetben mi határoztuk meg. E fejezet jelöléseivel

x = cs + v,

ahol v fehér zaj, melynek szórásnégyzete most σ^{2} = 0, 01 és

s_{i} = [0 2 1 1 1 1 1]^{T} ,

c = ( \begin{matrix}1 & sin(2πf_{1}t_{1}) & sin(2πf_{2}t_{1}) & sin(2πf_{3}t_{1}) & cos(2πf_{1}t_{1}) & cos(2πf_{2}t_{1}) & cos(2πf_{3}t_{1}) \\ 1 & sin(2πf_{1}t_{2}) & sin(2πf_{2}t_{2}) & sin(2πf_{3}t_{2}) & cos(2πf_{1}t_{2}) & cos(2πf_{2}t_{2}) & cos(2πf_{3}t_{2}) \\ ⋮ & & & & & & ⋮ \\ 1 & sin(2πf_{1}t_{i}) & sin(2πf_{2}t_{i}) & sin(2πf_{3}t_{i}) & cos(2πf_{1}t_{i}) & cos(2πf_{2}t_{i}) & cos(2πf_{3}t_{i}) \\ ⋮ & & & & & & ⋮ \\ 1 & sin(2πf_{1}t_{n}) & sin(2πf_{2}t_{n}) & sin(2πf_{3}t_{n}) & cos(2πf_{1}t_{n}) & cos(2πf_{2}t_{n}) & cos(2πf_{3}t_{n}) \end{matrix} ) ,

a folyamat zajától (w) most eltekintünk. Az adatsor tizenöt százalékát ezután kitöröltük a hozzájuk tartozó időpontokkal együtt, hogy a nem-ekvidisztáns időosztást modellezzük. A frekvenciák f_{1} = 10, 5, f_{2} = 23, 3, illetve f_{3} = 31, 7.

A nyers adatsorra (3.3. ábra) ekkor végrehajtottuk az algoritmust 1/4 Hz-es lépésenként vizsgálva a frekvenciákat. A becsült spektrum (3.4 ábra) alapján nyilvánvaló, hogy a folyamatban három periodikus komponens található, a kapott numerikus eredmények alapján pedig célszerűnek tűnhet a 9,25 Hz-től 11,25 Hz-ig, a 22,75 Hz-től 24 Hz-ig, valamint a 31 Hz-től 32,5 Hz-ig terjedő frekvencia-intervallumok részletesebb vizsgálata. Ezekre az intervallumokra újra végrehajtottuk az algoritmust, ám ezúttal 0,05 Hz-es lépésközökkel. Az alábbiakban az érdekesebb eredmények közül közlünk néhányat:

                              freq.  wsinˆ 2  wcosˆ 2
                              10.40  0.0148  0.0142
                              10.45  0.0219  0.0093
                              10.50  0.0280  0.0045
                              10.55  0.0314  0.0012
                              10.60  0.0315  0.0000

lasztunk ki és ezt a frekvencia-intervallumot finomítjuk és vizsgáljuk a következő lépésekben. Ezzel a többlépcsős módszerrel gyorsan érhető el bizonyos intervallumok finomítása. Alkalmazni többnyire olyankor ajánlatos, amikor szomszédos pontokhoz összemérhető súlyok tartoznak.

FEJEZET 3. SPEKTRUMBECSLÉS 20

3.3. ábra. Mesterségesen generált adatsor grafikonja. 3.4. ábra. A spetrum becslése az algoritmus egyszeri végrehajtása után.

Ezen adatok alapján nem lehet egyértelműen eldönteni, hogy melyik ferkvencia áll legközelebb a valódihoz, egyáltalán az sem látható, hogy csak egy frekvencia van jelen. Ha azonban élünk azzal a feltételezéssel, hogy egy ekkora frekvencia-intervallumba legfeljebb egy periodikus komponens frekvenciája esik, akkor célszerű a 10,45 Hz mellett dönteni, hiszen ehhez a frekvenciához tartozik a legnagyobb összesített súly. Technikai okokból ugyanis a táblázatban nem a súlyok, hanem a súlyok négyzetei jelennek meg. Még egy további finomítását is végrehajtottuk a frekvenciaintervallumoknak és erre még egyszer alkalmaztuk az algoritmust. Az ekkor kapott eredmények alapján, ha a legnagyobb összsúllyal rendelkező frekvencia-jelölteket választjuk, akkor a 10,48 Hz, 23,31 Hz, illetve a 31,69 Hz-es frekvenciák mellett döntünk, azaz a frekvenciák becsült vektora

f̂ = [10, 48 23, 31 31, 69].

Negyedszer alkalmazva az algoritmust, már csak ezekre az értékekre, a következő súlyok adódnak:

                              freq.  wsinˆ 2  wcosˆ 2
                              10.48  3.8094  1.3788
                              23.31  1.0971  0.6563
                              31.69  0.8854  1.1905

                               const.: 3.2291e-005,

ahol ezúttal a konstans vektor (a c mátrix első oszlopa) együtthatóját is feltüntettük. Ennek alapján a becsült vektor tehát

ŝ = [3, 2291 · 10^{−5} 1, 9518 1, 0474 0, 9410 1, 1742 0, 8101 1, 0911]^{T} .

Az együtthatókra vonatkozó becslés hibája (e = s − ŝ) ekkor

e = [−0, 0057 0, 0482 − 0, 0474 0, 0590 − 0, 1742 0, 1899 − 0, 0911]^{T} ,

így már sikeresen becsülhetjük a kiindulási adatsort (x̂ = cŝ) annak ellenére, hogy a koszinuszos tagok együtthatói esetében a hiba némileg nagyobb, mint a szinuszos tagok esetében (ld. 3.5. ábra).

A jobb összehasonlíthatóság kedvéért megjegyezzük, hogy abban az esetben, ha az eljárással az igazi frekvenciákhoz tartozó súlyokat becsüljük meg, azaz ha a bemeneti paraméterként megadott frekvenciák vektora

f = [10, 5 23, 3 31, 7],

akkor a becsült súlyokra az

                              freq.  wsinˆ 2  wcosˆ 2
                              10.5   4.2852  0.9122
                              23.31  1.0026  0.7534
                              31.69  1.0012  1.0845

                                  const.: 0.0000

eredmények adódnak. Így tehát ebben az esetben

ŝ = [0, 0052 2, 0701 1, 0013 1, 0006 0, 9551 0, 8680 1, 0414]^{T} ,

vagyis a becslés hibája

e = [−0, 0052 0, 0701 − 0, 0013 0, 0006 − 0, 0449 − 0, 1320 − 0, 0414]^{T} .

3.5. ábra. A valódi (folytonos), illetve a becsült spektrum alapján számolt adatsor (szaggatott).

3.4. Következmények

Megmutattuk, hogy a tárgyalt módszerrel valóban lehetséges egyfelől a spektrum hatékony becslése, másrfelől példát mutattunk arra, hogy a becsült spektrum alapján valóban lehetséges az eredeti adatsor hatékony visszaállítása. Láttuk, hogy a módszer sok esetben lényegesen hatékonyabb, mint az FFT, hiszen a becsült spektrumban a csúcsok jobban elkülönülnek, emellett sok olyan esetben is alkalmazható, amikor a gyors Fourier transzformáció csak korlátozottan, vagy egyáltalán nem.

Komoly problémát jelent azonban az adatsor spektruma alapján a valódi frekvenciák megfelelő kiválasztása, ráadásul a kiválasztás mechanizmusa nehezen automatizálható. További probléma a számítás művelet- és tárigénye, ami lényegesen nagyobb erőforrások alkalmazását teszi szükségessé, mint az FFT, vagy akár a Lomb-Scargle periodogram számítása.

Rejtett periodicitások

Ebben a fejezetben az idősorokban fellelhető periodikus komponensek keresésére használt statisztikai eljárásokat ismertetünk. Elsőként Fisher tesztjét tárgyaljuk. Fisher az adatsor periodogramjából származtatott valószínűségi változók eloszlásának pontos meghatározásával éri el, hogy a periodikus komponensek szignifikancia-szintje meghatározható.

A másodikként ismertetett Kolmogorov-Szmirnov próba egy változatának használhatósága is az adatsor periodogramján alapszik.

A klasszikus periodogram azonban csak ekvidisztánsan mintavételezett adatsorra számítható, így további próbák szükségesek annak eldöntéséhez, hogy a nem ekvidisztánsan mintavételezett adatsorokban mennyire szignifikánsak a periodikus komponensek. Ez utóbbi problémára egy – a Lomb–Scargle-periodogramon alapuló – tesztet ismertetünk.

4.1. Fisher-teszt: ismeretlen frekvenciájú rejtett periodicitá-

sok keresése

Fisher tesztje a folyamat periodogramjának vizsgálatán alapszik. Az x_{t} folyamat periodogramján az

I(ω_{j}) = n^{−1} |∑_{t=1}^{n} x_{t}e^{−itω_{j}} | ,

mennyiséget értjük, ahol ω_{j} = \frac{2πj}{n} .

Ha x_{t} normális eloszlású fehér zaj σ^{2} szórásnégyzettel, akkor jól ismert, hogy

\frac{I(ω_{k})}{σ^{2}} ∼ \frac{χ^{2}(2)}{2}, k = 1, …, q,

ahol

q = ⌊\frac{n − 1}{2}⌋ .

Mivel V_{k} sűrűségfüggvénye e^{−x}I_{[0,∞)}(x), így a V_{1}, …, V_{q} változók együttes sűrűségfüggvénye

f_{V_{1},...,V_{q}}(v_{1}, …, v_{q}) = ∏_{i=1}^{q} I_{[0,∞)}(v_{i}).

Ebből következik az ismert tétel, miszerint, ha x_{t} normális eloszlású fehér zaj, akkor az

Y_{i} = \frac{∑_{k=1}^{i} I(ω_{k})}{∑_{k=1}^{q} I(ω_{k})}, i = 1, …, q − 1

változók függetlenek és eloszlásaik egyenletesek a [0, 1] intervallumon, azaz

Y_{i} ∼ U(0, 1), ∀ i.

Ugyanezen feltételek mellett, ha Y_{0}-t zérusnak, Y_{q}-t pedig egynek definiáljuk valamint bevezetjük az

M_{q} = max_{1≤i≤q}(Y_{i} − Y_{i−1}) = max_{1≤i≤q} \frac{I(ω_{i})}{∑_{i=1}^{q} I(ω_{i})}

állandót, akkor

P (M_{q} ≤ a) = ∑_{j=0}^{q}(−1)^{j}\binom{q}{j} (1 − \frac{ja}{q})_{+}^{q−1} ,  (4.1)

ahol z_{+} = max(z, 0).

Fisher ezen eredmények alapján alkotta meg tesztjét, mellyel azt a hipotézist vizsgálta, hogy egy adott normális eloszlású fehér zajban vajon vannak-e ismeretlen frekvenciájú periodikus komponenesek. A teszt lényege abban áll, hogy akkor utasítja el a nullhipotézist, azaz azt, hogy a folyamat mentes periodikus komponensektől, ha a periodogram az átlagnál szignifikánsan magasabb értékeket is tartalmaz.

Bevezetve a ξ_{q} változót, melyet a

ξ_{q} = \frac{max_{1≤i≤q} I(ω_{i})}{q^{−1} ∑_{k=1}^{q} I(ω_{i})} = qM_{q}

összefüggés definiál, és felhasználva (4.1) egyenletet kapjuk, hogy

P (ξ_{q} ≥ a) = 1 − ∑_{j=0}^{q}(−1)^{j}\binom{q}{j} (1 − \frac{ja}{q})_{+}^{q−1} .  (4.2)

Ha ez a valószínűség α-nál kisebb, akkor a nullhipotézist α szinten elvetjük.

4.2. Kolmogorov–Szmirnov teszt alkalmazása

A Kolmogorov–Szmirnov próba is alkalmazható a korábban definiált nullhipotézis tesztelésére. Tekintsük ugyanis az Y_{i} változók tapasztalati eloszlásfüggvényét, és ellenőrizzük, hogy az mennyire tér el az egyenletes eloszlás F (x) = x eloszlásfüggvényétől. q > 30 esetén egy α szintű Kolmogorov– Szmirnov tesztben elvetjük a nullhipotézist, ha a tapasztalati eloszlásfüggvény kilép az

y = x ± k_{α}(q − 1)^{− \frac{1}{2}} , 0 < x < 1

tartományból, ahol k_{0,05} = 1, 36 és k_{0,01} = 1, 63.

Ezzel az eljárással ekvivalens, ha a normalizált kumulatív peridogramot,

C(x) = \{ \begin{matrix}0, & x < 1, \\ Y_{i}, & i ≤ x < i + 1, i = 1, …, q − 1, \\ 1, & x ≥ q, \end{matrix}

használva akkor utasítjuk el a nullhipotézist α szinten, ha bármely x ∈ [1, q]-re a C függvény kívül esik az

y = \frac{x − 1}{q − 1} ± k_{α}(q − 1)^{− \frac{1}{2}}

4.1. ábra. Egy mesterségesen generált, több harmonikus összetevőt is tartalmazó, fehér zajjal terhelt adatsor periodogramja és a különbőző k_{α}-khoz tartozó határok.

által rögzített határokon.

Fontos megjegyezni, hogy sem a Kolmogorov–Szmirnov, sem pedig a Fisher-féle teszt nem szolgáltat információt arra nézve, hogy milyen frekvenciájúak az esetleges periodikus összetevők. A frekvenciák megkeresésére ezért más módszereket kell keresni. A periodogram általában nem alkalmas ezek megállapítására, ugyanis zaj jelenlétében a valódi frekvenciákhoz tartozó csúcsokat gyakran hamis csúcsok nyomják el.

4.3. Nem ekvidisztáns mintavételezés esete

A fentiekben ismertetett elméletek alkalmazhatóságának egyik feltétele, hogy a mintavételezés időpontjai között eltelt idő állandó legyen. Ez a feltétel azonban ritkán teljesül. Sok esetben az észlelést végző berendezés hibája miatt esnek ki megfigyelési pontok, más esetekben pedig (pédául bizonyos csillagászati mérések esetén) a megfigyelés tárgya nem teszi lehetővé az ekvidisztáns mintavételezést.^{1}

A problémára egyfajta megoldást jelent a hiányzó adatpontok interpolációval történő becslése, azonban ez a módszer nem, vagy csak nagyon korlátozottan alkalmazható olyan esetekben, amikor sok egymást követő adatpont hiányzik. Különösen az olyan periodikus komponensek meghatározásánál jelent ez problémát, melyek hullámhossza összemérhető a hiányzó adatpontok intervallumának hosszával. Ezért a femtokémiában gyakori a logartimikus időosztású mintavételezés is.

Másfajta megoldást kínál a Lomb [9] féle sztenderdizált periodogram, melynek lényege, hogy nem frekvencia-intervallumok súlyát számítja ki, hanem pontonként rendel súlyokat az egyes frekvenciákhoz. A módszer ezen tulajdonsága miatt gyakran pontosabb eredménnyel szolgál, mint az FFT. A továbbiakban [14] jelöléseit használjuk.

1A dolgozat alapproblémáját jelentő femtokémiai mérések esetén a leggyakrabban nem ekvidisztáns mintavételezés vezet optimális információtartalomhoz.

A Lomb periodogram kiszámításához először a vizsgált adatsor számtani közepét és korrigált tapasztalati szórásnégyzetét kell kiszámítani. A szokásos jelölésekkel

x̄ = \frac{1}{N} ∑_{i=1}^{N} x_{i}, σ^{2} = \frac{1}{N − 1} ∑_{i=1}^{N}(x_{i} − x̄)^{2},

ahol x_{i}, i = 1, …, N a detektált értékek. A normalizált Lomb periodogramot az ω = 2πf körfrekvenciák esetén a

P_{N}(ω) = \frac{1}{2σ^{2}} (\frac{(∑_{j}(x_{j} − x̄) cos ω(t_{j} − τ))^{2}}{∑_{j} cos^{2} ω(t_{j} − τ)} + \frac{(∑_{j}(x_{j} − x̄) sin ω(t_{j} − τ))^{2}}{∑_{j} sin^{2} ω(t_{j} − τ)} )  (4.3)

összefüggés határozza meg. A képletben szereplő τ mennyiségnek fontos szerepe van, amennyiben a formulát invariánssá teszi a mintavételezési időpontokat tartalmazó t vektor eltolásaival szemben. További, mélyebb összefüggések megvilágításában is fontos szerephez jut. Lomb kimutatta, hogy ha τ-t a

tan(2πτ) = \frac{∑_{j} sin 2πωt_{j}}{∑_{j} cos 2πωt_{j}}

formulából származtatjuk, akkor a (4.3) egyenlet megegyezik azzal a legkisebb négyzetes becsléssel, mely egy adott ω-hoz tartozó harmonikus komponens

x(t) = A cos ωt + B sin ωt

modellhez történő illesztésével adódik.

4.4. Periodicitások nem ekvidisztáns minta esetén

A periodogram – akárcsak a korábbi esetben – lehetőséget nyújt a mintában jelenlévő periodikus komponensek kiszűrésére. Igen gyakori az az eset, amikor a mért adatsor harmonikus összetevőkből és fehér zajból áll; a korábbi modellek is mind ezzel a feltételezéssel éltek. Ahhoz azonban, hogy egy spektrumban szereplő csúcsról megállapítható legyen valamiféle szignifikancia, ismerni kell a periodogram értékeinek az átlagtól vett eltérését, illetve az eltérések eloszlását.

Scargle [10] megmutatta, hogy a szokásos nullhipotézis teljesülése mellett a P_{N} valószínűségi változó eloszlása egy egység várható értékű exponenciális eloszlás lesz. Ebből következik, hogy M független frekvenciát megvizsgálva annak a valószínűsége, hogy az egyik z-nél nagyobb értéket ad (1 − e^{−z})^{M}. Így tehát a másodfajú hiba valószínűsége

P (> z) = 1 − (1 − e^{−z})^{M}.  (4.4)

A másodfajú hiba kis valószínűsége szignifikáns csúcsot jelez.

A szignifikancia kvantitatív meghatározásához M ismerete szükséges. A fontosnak számító (kis) értékekre a fenti valószínűség jól közelíthető a jobb oldal Taylor-sorával, azaz

P (> z) ≈ Me^{−z},  (4.5)

4.2. ábra. Ekvidisztánsan mintavételezett adatsor Lomb-féle periodogramja (folytonos) és FFT-je (szaggatott).

tehát kis értékek esetén a valószínűség jó közelítéssel az M lineáris függvénye. Szokásos szignifikanciaszintek esetén akár a becslés ötven százalékos hibája is elfogadható, hiszen a szintek akár nagyságrendben is eltérhetnek. Mindebből az a következtetés vonható le, hogy M-et nem kell „túl” pontosan meghatározni.

Horne és Baliunas [13] kutatásai alapján elmondható, hogy M elsősorban az adatpontok számától és a mintavételezési intervallumtól függ. Számos kísérlet alapján arra a megállapításra jutottak, hogy ekvidisztáns mintavételezés esetén az optimális M az adatpontok számával közel megegyezik, sőt ez a helyzet abban az esetben is, ha az adatpontok véletlenszerűen oszlanak el a vizsgált időintervallumon. Jelentős eltérés a vizsgálat szerint csak akkor mutatkozik, ha a megfigyelési időpontok szoros csoportokat alkotnak, ez az eset azonban jelen dolgozat esetében nem meghatározó fontosságú.^{2}

2M értékét a számítások egyszerűségének megtartásáért mindig a minta nagyságával megegyezőnek választottuk.

Periodicitások a mért adatsorokban

A fejezetben megvizsgáljuk, hogy a mérések adatsorai alapján statisztikailag igazolható-e periodikus komponensek jelenléte. Az ekvidisztánsan mintavételezett adatsorok esetén elsőként a periodogramot számítjuk ki és annak alapján végezzük el a Fisher, illetve a Kolmogorov–Szmirnovpróbát. Egyenetlen mintavételezésű adatsorokra pedig a Lomb–Scargle periodogram alapján teszteljük a nullhipotézist, mely szerint nincsenek periodikus komponensek a vizsgált mintában.

Kimutatjuk, hogy mindkét esetben a periodikus komponensek jelenléte statisztikai értelemben magas szignifikancia-szint mellett igazolható a vizsgált adatsorok esetén.

5.1. Ekvidisztáns mintavételezésű adatsorok

Elsőként ekvidisztáns mintavételezésű adatsorok (5.1. ábra.) periodogramjait vizsgáljuk. A kiszámított periodogramok alapján elvégezzük a Fisher tesztet, majd – az eredményt ellenőrizendő – összehasonlítjuk a Kolmogorov–Szmirnov-próba eredményeivel. Látni fogjuk, hogy a tesztek azonos jelentésű eredményeket adnak.^{1}

Ebben az alfejezetben olyan adatsorokat vizsgálunk, melyek nem valódi mérésekből, hanem valódi mérések szimulációiból^{2} származnak. Így lehet biztosítani például a mintavételezési intervallumok állandóságát.

5.1.1. Periodogramok

Az előző fejezetben leírtak alapján kiszámoltuk az adatsorok periodogramjait. A periodogramokat (5.2. ábra.) a nullához közeli frekvenciákhoz tartozó magas értékek dominálják. Ebből azt a következtetést lehet levonni, hogy az adatsorokban is az alacsony frekvenciájú összetevők a meghatározó fontosságúak.

A periodogramok alapján várható, hogy a tesztek statisztikailag igazolják a periodikus komponensek jelenlétét. Fisher valószínűségeire zérushoz közeli értékek várhatók, a Kolmogorov–Szmirnovtesztben pedig az, hogy a tapasztalati eloszlásfüggvény, illetve a sztenderdizált kumulatív periodogram elhagyja a kis α szinthez tartozó sávokat is.

1A számításokat a period.m, fish.m, illetve a ks.m-ben definiált MATLAB programokkal végeztük.

2A három vizsgált adatsort a mellékelt inF1.DAT, inF2.DAT és az inF3.DAT fájlok tartalmazzák.

5.1. ábra. Ekvidisztánsan mintavételezett szimulált adatsorok.

5.1.2. Fisher-teszt és a Kolmogorov–Szmirnov próba

A periodogramok alapján Fisher módszeréhez híven meghatározzuk a szóban forgó (4.2) valószínűségeket. A három vizsgált adatsor esetén ezek

                               adatsor   P értéke
                                  1.         0
                                  2.     3.253e-014
                                  3.         0.

A táblázatban szereplő zérus értékek azt jelzik, hogy a vizsgált valószínűség nullától való távolsága kisebb, mint a MATLAB-pontosság.^{3} Ez egyben azt is jelenti, hogy a periodikus komponensek jelenléte mindhárom esetben szinte bármilyen értelmes α szint mellett igazolható.

Ezen megállapítást megerősíti a Kolmogorov–Szmirnov-próba is, hiszen az adatsorok tapasztalati eloszlásfüggvénye már a vizsgált intervallum legelején kilép még az α = 0, 01 szinthez tartozó k_{α} = 1, 63 által meghatározott

y = x ± k_{α}(q − 1)^{− \frac{1}{2}} , 0 < x < 1

határok közül is (5.3. ábra.).

Mindkét teszt pozítív eredményt hozott, így a nullhipotézist, miszerint az adatsorok periodikus összetevőktől mentesek, el kell vetnünk.

3A tesztkonfiguráción a dupla MATLAB pontosság körülbelül 10^{−325}-ös nagyságrendű volt.

FEJEZET 5. PERIODICITÁSOK A MÉRT ADATSOROKBAN 30

5.2. ábra. Az 5.1 ábrán látható adatsorok periodogramjainak logaritmusai.

5.3. ábra. Az 5.1 ábra adatsorainak alapján számolt sztenderdizált kumulatív periodogramok. A szaggatott vonalak a k_{0,01} = 1, 63-hoz tartozó Kolmogorov–Szmirnov-határokat jelölik.

5.4. ábra. Valódi mért adatsorok, 2080nm és 558nm nem ekvidisztáns, 800nm ekvidisztáns mintavételezésű.

5.2. Nem ekvidisztáns mintavételezésű adatsorok

Ebben az alfejezetben nem ekvidisztáns mintavételezésű adatsorokat^{4} vizsgálunk. A rejtett periodicitások fellelésére először kiszámítjuk az adatsorok Lomb-periodogramját, majd ebből Scargle módszerével határozzuk meg az egyes szignifikancia-szintekhez tartozó kritikus értékeket.^{5}

Akárcsak az előző alfejezetben, itt is kimutatjuk, hogy a periodicitások jelenléte statisztikailag igazolható, magas szignifikancia-szint mellett.

5.2.1. Adatsorok és Lomb-féle periodogramjaik

A vizsgált adatsorok (5.4. ábra.) tehát valódi mérésekből származnak. A mintavételezés két esetben logaritmikus, pontosabban mindkét esetben a mintavételezési időpontok logaritmikusan ritkulnak, a harmadik minta pedig ekvidisztáns mintavételezésű. A Lomb-periodogramra és különösen a Scargle-módszerre vonatkozó feltételek ebben az esetben teljesülnek. Emlékeztetőül: a csúcsok szignifikanciájára vonatkozó becslés abban az esetben válik bonyolulttá, ha a mintavételezési időpontok szoros csoportokat alkotnak nagy „foghíjakkal.”

Az adatsorok Lomb-periodogramjából (ld. 4.3) nem csak a periodikus komponensek jelenlétére, hanem azok frekvenciájára is következtetni lehet. A periodogramokat (5.5. ábra.) vizsgálva nyilvánvalónak tűnik, hogy vannak az adatsorban periodikus komponensek. Az adatsorokhoz és az egyes szintekhez tartozó kritikus értékeket táblázatban foglaltuk össze:

4A szóban forgó adatsorokat a 2080nm.dat, 800nm.dat, illetve az 558nm.dat fájlok tartalmazzák. Ezek valódi mérésekből származó adatsorok.

5Az ebben az alfejezetben felhasznált MATLAB-programok listája: lomb.m, lombtest.m.

5.5. ábra. Valódi adatsorok Lomb-féle periodogramjainak logaritmusai. A vízszintes vonal a 13,92nak megfelelő kritikus értéket ábrázolja.

                           dataset   level      value
                           2080nm    0.05 0.01 0.001     7.1901 8.8197 11.1268

                            558nm    0.05 0.01 0.001     7.2191 8.8487 11.1558

                            800nm    0.05 0.01 0.001     7.6800 9.3097 13.9198

Látható, hogy egyfelől a táblázat értékei nem különböznek számottevően, ami az adatsorok alaki hasonlóságának^{6} következménye, másfelől pedig, hogy a Lomb-féle periodogramban a csúcsok túlnyomó többsége messze meghaladja még a 0,001 szinthez tartozó kritikus értékeket is. Ilyen szinten vethető el tehát a nullhipotézis, miszerint az adatsor mentes periodikus komponensektől.

5.2.2. Lomb-periodogram és frekvenciák

A későbbiekben hasznát vesszük annak a ténynek, hogy a Lomb-féle periodogram olyan esetekben is használható a periodikus komponensek frekvenciáinak meghatározására, amikor a gyors Fourier transzformáció nem számolható.

6Az alaki hasonlóság annak a következménye, hogy mindhárom adatsorban az időbeli változásokat ugyanaz a fizikai folyamat szabja meg.

A Bayes spektrumbecslés hibája, hogy nehéz a megfelelő frekvenciákat kiválasztani. Hasznos ezért a konkrét esetekben a Lomb-periodogramot is kiszámolni, hogy a frekvenciák meghatározása egyszerűbb legyen.

Ha a frekvenciák viszonylag jó pontossággal meghatározásra kerültek, már hatékonyan alkalmazható a spektrumbecslés és felhasználható arra, hogy egy rögzített függvény spektrumának és az impulzus függvény spektrumának becsléséből egy használható közelítést kapjunk a folyamatot ténylegesen leíró függvény spektrumára, lehetővé téve ezzel a dekonvolúciót.

Dekonvolúció és spektrumbecslés

Ebben a fejezetben megkíséreljük végrehajtani a dekonvolúciót a korábban ismertetett, jelfeldolgozásban használt módszerek segítségével pontosabban, mint ahogy az a hagyományos módszerekkel lehetséges volt.

Elsőként meghatározzuk a folyamatokat összetevő harmonikus függvények frekvenciáit a Bayes spektrumbecslés, illetve a Lomb-periodogram segítségével, majd megbecsüljük ezen komponensek súlyait. Hasonló analízist hajtunk végre az impulzusfüggvényen. A becslés pontosságának növeléséhez némely esetben az adatsort az analízis előtt Gans javaslatának megfelelően alakítjuk át.

Végül eddigi eredményeink segítségével megbecsüljük az eredeti folyamat spektrumát. Példát mutatunk arra, hogy ez a becslés igen pontos lehet még egyszerű módszerek esetén is. Utoljára a becsült spektrum alapján megpróbáljuk előállítani az eredeti függvényeket.

6.1. Harmonikus összetevők

Mint azt már korábban is említettük, a dekonvolúció elvégzésekor a harmonikus összetevők frekvenciáinak meghatározása okozza a legnagyobb problémát. A vizsgált esetekben ezért a Bayes spektrumbecslés és a Lomb-peiodogram segítségével próbáljuk majd a frekvenciákat meghatározni.

Hogy az eredményeket ellenőrizni tudjuk, a módszert olyan esetekben alkalmazzuk, amikor az eredeti függvény is rendelkezésre áll, azaz szintetikus adatsorokat használunk.^{1}

Kiindulásképpen megvizsgáljuk az adatsor (példánkban az inF2.dat fájl adatsora) Lomb-peri-odogramját és a becsült spektrumot (ld. 6.1. ábra.). Az ábrából látható, hogy mind a két módszer nagyjából egyforma arányban súlyozza az egyes frekvenciákat, a nullához közeli frekvenciák nyilvánvalóan dominálják a folyamatot. Ebben az esetben igen nehéz egyes konkrét frekvenciákat kiválasztani a spektrumból, mint a folyamatot alkotó frekvenciát, ezért nem is tesszük ezt. Megmutatjuk viszont, hogy az impulzus- és az észlelt függvény spektrumainak segítségével jó becslést kapunk a valódi függvény spektrumára. Nevezetesen azt látjuk be, hogy az

O = \frac{I}{S} \frac{|S|}{|S| + α},

egyenlőség jó közelítéssel fennál, ahol I, S, O most nem az i, s, o függvények Fourier transzformáltjait, hanem ezen függvények Bayes spektrumbecsléseit jelölik, α pedig megfelelően választott konstans. Így ez egy gyakran használt, egyszerű dekonvolúciós módszer olyan változata, ahol spektrumbecslést alkalmazunk.

1Ilyen adatsorokat tartalmaznak az inF1-3.DAT, oF1-3.DAT, illetve a Pulse30_100.dat fájlok.

6.1. ábra. Az inF2.dat adatsorhoz tartozó Lomb-periodogram és az első menetben becsült spektrum logaritmusai.

Az O függvény ilyen meghatározásában a fenti szorzat második tényezője mint egy aluláteresztő szűrő működik. Alacsony frekvenciák esetén, ahol az impulzusfüggvényből spektrumbecsléssel származtatott S függvény értékei viszonylag nagyok, a szorzótényező értéke α megfelelő választásával egyhez közeli, míg magasabb frekvenciákon, ahol S értéke kicsi, nullához van közel. Így α megválasztásával az áteresztés frekvenciája szabályozható.

Célunk példával illusztrálni, hogy az általunk javasolt módszer annyiban jobb a hagyományos Fourier-transzformáción alapuló eljárásnál, hogy még α kis értékeinél sem jelentkeznek a spektrum magas frekvenciájú részein a Fourier-transzformációt alkalmazó módszerek esetén jellemző igen magas csúcsok, melyek lehetetlenné teszik a keresett függvény előállítását.

Hogy célunkat elérjük, meghatározzuk a fényimpulzust jellemző függvény, valamint az eredeti függvény spektrumát is. Megjegyezzük, hogy a felhasznált szintetikus adatsorok ekvidisztáns mintavételezésűek, így jó alapot szolgáltatnak az összehasonlításhoz, hiszen végezhetők rajtuk klasszikus FFT-számítások is, ha előbb Gans módszerével átalakítjuk őket.

6.2. A valódi és az érzékelt függvény spektruma

Tekintsük a 6.2 ábrát! Jól látható, hogy az O függvény és az \frac{I}{S} \frac{|S|}{|S|+α} függvény a teljes becsült frekvencia-tartományon azonos nagyságrendű (ezúttal α = 1 helyettesítéssel). Látható az is, hogy a spektrum elején a két függvény gyakorlatilag megegyezik. Mivel ezek az értékek a dominánsak a spektrumban, ezért joggal feltételezhető, hogy az itteni egyezés elegendő a dekonvolúció kellőképpen pontos elvégzéséhez.

FEJEZET 6. DEKONVOLÚCIÓ ÉS SPEKTRUMBECSLÉS 36 6.2. ábra. Az eredeti o függvény becsült spektrumának logaritmusa és a detektált mintából és az impulzusból származtatott becslése.

6.3. ábra. A 6.2 ábrán látható spektrum egyik kritikus szakasza, a súlyok logaritmusaival.

6.4. ábra. Súlyok a 6.2 ábrán látható spektrum egy kritikus szakaszán.

6.3. A dekonvolúció végrehajthatósága

Az ábrák alapján kijelenthetjük, hogy a becsült spektrumokra mondott

O = \frac{I}{S} \frac{|S|}{|S| + α},

egyenlőség tartható. Az az állítás is igaz, miszerint az eredeti torzítatlan függvény becsült spektrumából az eredeti függvény jó közelítéssel előállítható (ld. 6.5. ábra).

Egyedüli kérdésnek marad, hogy vajon elő lehet-e állítani az impulzusfüggvény és a megfigyelt függvény spektruma alapján az eredeti adatsort. A tapasztalat azt mutatja, hogy akár az első becslés után is (tehát a vizsgált frekvencia-tartomány szűkítése nélkül) jó eredmények adódnak (ld. 6.6. ábra). Az előállított függvény valóban értékes információt hordoz, nem válik olyan zajossá, mint amilyen a Fourier-transzformáció alkalmazása tenné.

Fontos továbbá azt is megjegyezni, hogy bár a spektrumból a frekvenciák kiválasztása nehéz, nagyban javítható vele a becslés, amennyiben a zaj nagy része kiszűrhető, hiszen a tárgyalt esetekben az alacsony frekvenciák dominánsak.

6.4. Következmények

Kimutattuk tehát, hogy a javasolt dekonvolúciós módszer végrehajtható lépcsőfüggvények esetében is, anélkül, hogy a más módszerek által keltett hibák megjelennének. A becslés ráadásul rugalmasan és hatékonyan tovább javítható, ha sikerül meghatározni a folyamatot domináló frekvenciákat.

FEJEZET 6. DEKONVOLÚCIÓ ÉS SPEKTRUMBECSLÉS 38

6.5. ábra. Az eredeti adatsor (oF2), és a becslés inverziójával előállított függvény. 6.6. ábra. Az inF2 adatsor alapján készített becslés első lépésben.

Összefoglaló és további tennivalók

7.1. Összefoglaló

A dolgozatban egy olyan módszert ismertettünk, melynek segítségével a Fourier-transzformáción alapuló dekonvolúciós eljárások az eddiginél pontosabban végezhetők.

Betekintést nyújtottunk a probléma femtokémiai hátterébe, ismertettük, hogy milyen jelenségek kapcsán merül fel a konvolúció és a dekonvolúció problémája, valamint hogy ezek során milyen függvények vizsgálata szükséges. Különösen fontosak ezek közül a lépcsőszerű függvények, melyek a legtöbb kémiai reakciót jellemzik.

A lépcsőfüggvényekre azonban nem alkalmazhatók közvetlenül a Fourier-transzformációt használó algoritmusok, hiszen a gyors Fourier vagy a diszkrét Fourier-transzformáció alkalmazhatóságának előfeltétele, hogy a szóban forgó függvény periodikus legyen. Fourier-transzformáció alkalmazása nem periodikus függvények esetében a spektrum torzulásához vezet, ugyanis a spektrumban megjelennek a tökéletes lépcsőfüggvényhez tartozó frekvenciák. Három olyan módszert ismertettünk, mely lépcsőfüggvények periodikus függvényekké történő átalakítását célozza. Kiemelten foglalkoztunk és némileg általánosítottuk a Gans-Nahman módszert, amelyet a későbbiekben a spektrum becslésekor használtunk fel. A másik két, Samulon és Nicolson által jegyzett, módszer szintén széles körben használatos. Ismertettük Waldmeyer eredményét, mely szerint a három módszer ugyanarra a numerikus megoldásra vezet.

A függvények spektrumát egy olyan Bayes becslést használó algoritmussal becsültük, mely alkalmazható olyan esetkben is, amikor a Fourier-transzformált nem számolható. A Kálmán-szürőt használó eljárás az idő szempontjából nem ekvidisztansan elhelyezkedő adatsorok esetén is alkalmazható, sőt a folyamatnak nem szükséges stacionáriusnak sem lennie. Az algoritmus szinte valamennyi paramétere a priori módosítható, tág teret biztosítva ezáltal a legkülönfélébb javítások számára. Mi kizárólag a bemeneti frekvencia-intervallumok tekintetében aknáztuk ki ezen lehetőségeket.

Példát mutattunk arra, hogy a javasolt módszer által becsült spektrumból az eredeti függvény jól visszaállítható. Megemlítettük a módszer néhány gyenge pontját is. Nem elég az algoritmust egyszer alkalmazni egy adott adatsorra, mivel a tapasztalat szerint az elsőre becsült spektrum (tehát, amikor nem rendelkezünk semmiféle a priori ismerettel) csak arányaiban jó becslés, az egyes frekvenciák valódi súlyaira azonban pontatlan eredményt ad. Az első alkalmazás így inkább csak a kritikus frekvenciák meghatározására használható. Ráadásul – lévén az algoritmus négyzetes tárigényű a vizsgált frekvenciák számában (ami viszont az adtpontok számának függvénye) – az algoritmus futási ideje jóval nagyobb annál, mint amit a Fourier-transzformált számítása igényel. További problémát jelent a valódi frekvenciák kiválasztása a spektrum alapján, ha ezek egymáshoz közel helyezkednek el. Egy-két esetben szélhatás is jelentkezik, a spektrum magas frekvenciájú részein jelenhetnek meg hamis csúcsok. Tapasztalataink szerint az adatsorok Gans-féle módosítása jelentősen csökkenti a jelentkező szélhatásokat.

Megvizsgáltuk azt is, hogy vajon statisztikai értelemben igazolható-e egyáltalán periodikus komponensek jelenléte a vizsgált adatsorokban. Több módszert ismertettünk, melyek ekvidisztáns és nem ekvidisztáns mintavételezésű adatsorok esetén szolgáltatnak próbákat e probléma eldöntésére. Ekvidisztáns mintavételezésű adatsorok esetén használható Fisher tesztje, mely a vizsgált adatsor periodogramjából származtatott valószínűségi változók eloszlását használja fel, valamint a Kolmogorov–Szmirnov-próba egy változata. Nem ekvidisztáns mintavételezésű adatsorok esetén a Lomb periodogram segítségével konstruálható próba. A Lomb-periodogram segítségével Scargle módszerével adtunk meg kritikus értékeket.

Ezekkel a módszerekkel vizsgáltuk meg az összes olyan adatsort, melyekkel dolgoztunk. Mind az ekvidisztáns, mind pedig a nem ekvidisztáns mintavételezésű adatsorok esetén magas szignifikanciaszint mellett volt igazolható a harmonikus összetevők jelenléte.

Legvégül egy Fourier-transzformáción alapuló dekonvolúciós eljárást hajtottunk végre azzal a módosítással, hogy a Fourier-transzformáltak helyett a képfüggvény és az impulzus-függvény becsült spektrumait használtuk fel. Előbb összehasonlítottuk ezt a becsült spektrumot az eredeti adatsor alapján becsült spektrummal. Mivel az egyezés jónak volt mondható különösen a kritikus (nullához közeli) frekvencia-tartományok esetén, joggal volt várható, hogy a dekonvolúció sikerülhet. Előállítottuk a becsült spektrum alapján az eredeti függvény közelítését és elég jó egyezést kaptunk. Fontos megjegyezni, hogy az egyezés annak ellenére kielégítőnek mondható, hogy a spektrumot csak egy lépésben, a priori ismeretek nélkül végeztük. Várható tehát, hogy ennél jobb közelítések is viszonylag egyszerűen előállíthatók, például úgy, hogy a Lomb-periodogram és a becsült spektrum segítségével első körben a vizsgát frekvenciák körét szűkítve javítjuk a becslést.

7.2. További vizsgálati lehetőségek

A javasolt módszer igen érzékeny a kiválasztott frekvenciák pontosságára. Sokban segítené az adatsorok elemzését, ha automatikus módszert lehetne találni a tényleges frekvenciák megtalálására, ugyanis ez jelenleg szubjektív döntés alapján történik.

Fontos további cél a spektrumbecslés gyorsítása. Az algoritmus négyzetes a vizsgált frekvenciák számában, azonban a spektrumot ezek közül csak néhány, legfeljebb néhány tucat dominálja, azaz rengeteg fölösleges számítást végzünk. Javítható lenne az algoritmus a frekvencia-felosztás lépésenkénti sűrítésével. Így ugyan többször kell az algoritmust végrehgajtani, azonban néhány előzetes kísérlet alapján elmondható, hogy valószínűleg jelentős időmegtakarítás érhető el.

Elképzelhető, hogy a Gans-átalakítás alapötletét felhasználva, vagyis azt, hogy a mért adatsort önmagával összefűzve állítsunk elő periodikus mintát, az eredeti módszer által előállított adatsorokhoz képest jobb tulajdonságokkal rendelkező adatsorok állíthatók elő. Ígéretes átalakításnak látszik az eredeti adatsor függőleges tengely menti tükörképének az eredeti adatsorhoz történő hozzáfűzése.

Irodalomjegyzék

[1] Á. Bányász, E. Mátyus, E. Keszei, „Deconvolution of Ultrafast Kinetic Data with Inverse Filtering,” Radiat. Phys. Chem., accepted for publication.

[2] Á. Bányász, G. Dancs, and E. Keszei, „Direct Deconvolution of Femtosecond Kinetic Data”

[3] W. L. Gans and N. S. Nahman, „Continous and Discrete Fourier Transforms of Steplike Waveforms,” IEEE Trans. Instrum. Meas., vol. IM-31, pp. 97-101, Jun. 1982.

[4] Jürg Waldmeyer, „Fast Fourier Trasnsform for Step-Like Functions: The Synthesis of Three Apparently Different Methods,” IEEE Trans. Instrum. Meas., vol. IM-29, pp. 36-39, Mar. 1980.

[5] H. A. Samulon, „Spectrum Analysis of Transient Response Curves,” Proc. IRE, vol. 39, pp. 175-186, 1951.

[6] A. M. Nicolson, „Forming the Fast Fourier Transform of a Step Response in Time-Domain Metrology,” Electron. Lett., vol. 9, pp. 317-318, 1973.

[7] J. R. Laguna, communication presented at L’école d’été du Trégor (Lannion, France, July 1978)

[8] Yuan Qi, T. P. Minka and R. W. Picard, „Bayesian Spectrum Estimation of Unevenly Sampled Nonstationary Data,” IEEE, ICASSP, May 2002.

[9] N. R. Lomb, „Least-squares Frequency Analysis of Unevenly Spaced Data,” Astrophysical and Space Science, pp. 447-462, 1976.

[10] J. D. Scargle, „Studies in Astronomical Time Series Analysis ii. Statistical Aspects of Spectral Analysis of Unevenly Sampled Data,” Astrophisical Journal, pp. 835-853, 1982.

[11] G. L. Bretthorst, „Bayesian Spectrum Analysis and Parameter Estimation,” Lecture Notes in Statistics, 48. Springer-Verlag, 1988.

[12] P. J. Brockwell and R. A. Davis, „Time Series: Theory and Methods,” pp. 320-332, Springer-Verlag, 1987.

[13] J. H. Horne and S. L. Baliunas Astrophysical Journal, vol. 302, pp. 757-763, 1986.

[14] „Numerical Recipes in C: The Art of Scientific Computing,” pp. 575-579, 1992.

[15] A. Bennia and S. M. Riad, „Filtering Capabilities and Convergence of the Van-Cittert Deconvolution Technique,” IEEE Trans. Instrum. Meas. vol. 41, pp. 246-250, April 1992.

A. Függelék

Forráskódok

A.1. Spektrumbecslő algoritmus

Az alábbi algoritmus számolja ki egy x adatsorhoz a becsült spektrumot. Bemeneti paraméterként az adatsort, a hozzá tartozó időpontokat tartalmazó t vektort, illetve a vizsgálandó frekvenciák f vektorát kell megadni. Az idővektor és a frekvenciavektor adható nullának, ekkor a program a mintát ekvidisztáns mintavételezésűként kezeli, méghozzá úgy, mintha a megfigyelés egységnyi ideig tartott volna, a frekvenciavektor pedig ebben az esetben úgy áll elő, hogy a benne szereplő frekvenciák egyenletesen helyezkednek el a zéró frekvencia és a között a frekvencia között, melyet az adatpontok számának felében határozunk meg.^{1}

%%% bse.m %%%
% the function returns with the estimated spectrum (MM), the corresponding
% frequency-vector (freq), and the matrix that consists of trigonometric
% functions (c)
function [MM,freq,c]=bse(x,t,f)
if size(x,1)==1 x=x'; n=size(x,1); else n=size(x,1); end;
if t == 0 t=1/n:1/n:1; elseif size(t,2)==1 t=t'; end; % generates t
% calculates the max. freq. that is sane to check
mint = Inf;
for i=2:n;
    if ((t(i)-t(i-1))==0)
    elseif((t(i)-t(i-1))<mint) mint=t(i)-t(i-1);
    end;
end;
maxf = 1/2/mint;
k = 4; % determines the fineness of the partition of the freq. interval
N = size(x,1); % same as n previously (remained here for historical reasons)
if f==0 M = k*round(N/2); else M = length(f); end; % def. length of f

1Ez az eljárás – ha nem is precízen – a Nyquist frekvenciáig ellenőrzi a frekvenciákat. A frekvencia-intervallum felosztásának finomsága a k állandóval szabályozható.

%uncomment the following line if the sample is not equidistant
%M=round(k*maxf);
if f==0 for j=1:M; freq(1,j) = j/k; end; % generates f if it is not
else freq = f; end;                      % given in the input
sigma   =   1; % actually this is sigma^2
z       =   1; % scale parameter
% initiate variables
V       =   zeros(2*M+1,2*M+1);
Gamma   =   zeros(2*M+1,2*M+1);
P       =   zeros(2*M+1,2*M+1);
K       =   zeros(2*M+1,1);
m       =   zeros(2*M+1,1);
c       =   zeros(N,2*M+1);
% calculation of the c matrix the columns of which are trigonometric
% functions, the first column is constant. note that it is an overkill to keep
% c permanently in the memory, each vector could be generated on-the-fly
% when needed, however the code is simpler and more convenient this way.
c(:,M+2:2*M+1) =    cos(2*pi*(t(1:N))'*freq);
c(:,2:M+1) =        sin(2*pi*(t(1:N))'*freq);
c(:,1) =            ones(size(c(:,1)));
% for the sake of integrity of notations we introduce the following variables,
% which is a must, for MATLAB does not cope with the zero index
m0      =   zeros(2*M+1,1);
V0      =   eye(2*M+1);
Z       =   z*eye(2*M+1);
Gamma0  =   zeros(2*M+1);
P0      =   V0 + Gamma0;
% the following four lines can not be included in the below for-loop because
% of the MATLAB indexing-problem
K(:,1)    =   V0*c(1,:)'*inv(c(1,:)*V0*c(1,:)'+sigma);
V(:,:)    =   (eye(2*M+1)-K(1)*c(1))*V0;
m = m0 + K(1)*(x(1)-c(1)*m0);
Gamma(:,:)=   Z;
% the actual updates follow here:
for i=2:N;
    P      = V  +   Gamma;
    Gamma  = Z*(t(i)-t(i-1));

    K      = P*c(i,:)'*inv(c(i,:)*P*c(i,:)'+sigma);
    V      = (eye(2*M+1)-K*c(i,:))*P;
    m      = m + K*(x(i)-c(i,:)*m); % m holds the values we are interested in
end;
% to obtain the estimated power of the spectrum the square of the following
% has to be calculated
MM = m;

A program igen tárigényes, sok nagyméretű mátrix van permanensen a memóriában. Javítható a helyzet azzal, ha a trigonometrikus függvényeket tartalmazó c mátrixot nem számítjuk előre ki, csak az aktuális oszlopot használjuk, amikor szükség van rá. Más mátrixok, mint a V, a Gamma, vagy a P, gyakran ritka mátrixok, így sokkal kisebb tárigényűek, mint amekkorát tényleges méretük indokolttá tesz.

A.2. Lomb periodogram

Az alábbi kód egy adatsor Lomb-féle periodogramját számítja ki, illetve – Scargle módszere alapján – meghatározza az előre definiált α-szintekhez tartozó kritikus értékeket. Bemeneti paraméterként az adatsor értékeit tartalmazó x, illetve a mintavételezési időpontokat tartalmazó t vektort kell megadni. Mint a spektrumbecslés esetén, itt is adható a t = 0 érték a bemenetben. Ekkor a program az időintervallumot egység hosszúnak tekinti, a mintavételezést pedig ekvidisztánsnak.

%%% ls.m %%%
% returns the Lomb-periodogram (PN), and a vector that contains the critical
% values (level) for pre-defined alpha-levels.
function [PN,level]=ls(x,t)
if size(x,1)==1 x=x'; end; % a column vector is exprcted
n = size(x,1);
if t==0 t=1/n:1/n:1; end;
if size(t,2)==1 t=t'; end; % a row vector is expected
mn = mean(x);
s2 = cov(x)/(n-1);
% determine the max freq. that is sane to calculate
mint = Inf;
for i=2:n;
    if ((t(i)-t(i-1))==0)
    elseif((t(i)-t(i-1))<mint) mint=t(i)-t(i-1);
    end;
end;
% different methods of calculating the max freq. we use the static value since
% this method gave the most sane results

%maxf = round(1/2/mint);
%maxf = 4*maxf;
maxf = 300;
% calculation of the periodogram itself
for w=1:maxf;
    tau = 1/(2*w)*atan(sum(sin(2*w*t))/sum(cos(2*w*t)));
    PN(w) =
    1/(2*s2)*((sum((x-mn)'.*cos(w*(t-tau)))).^2/sum((cos(w*(t-tau))).^2) +
    (sum((x-mn)'.*sin(w*(t-tau)))).^2/sum((sin(w*(t-tau))).^2));
end;
% this program invokes lombtest.m at this point to calculate the critical
% values. note that lombtest.m is a script, not a function.
lombtest;

Az alábbiakban az ls.m által meghívott script kódja található, mely a kritikus értékeket határozza meg.

%%% lombtest.m %%%
% this script calculates the critical values for the Lomb-Scargle-test.
lombM = size(x,1);
% pre-defined alpha-levels
c1=.05; c2=.01; c3=.001; c4=.0001;
% the critical values
level(1) = -log(1-(1-c1)^(1/lombM));
level(2) = -log(1-(1-c2)^(1/lombM));
level(3) = -log(1-(1-c3)^(1/lombM));
level(4) = -log(1-(1-c4)^(1/lombM));

A kritikus értékeket Scrgle módszere alapján határozzuk meg a (4.4) egyenlet alapján. Emlékeztetőül:

α = P (> z) = 1 − (1 − e^{−z})^{M},

amiből a programban szereplő

z = − ln (1 − (1 − α)^{\frac{1}{M}} )

összefüggés elemi átalakításokkal következik.

Némileg egyszerűbb formulát kapunk, ha a valószínűségre álló egyenlet helyett annak Taylorsorát (4.5) használjuk fel. Ekkor a kritikus z értékre az alábbi összefüggés adódik

z = − ln \frac{α}{M} .

A.3. Periodogram, Fisher-teszt és Kolmogorov–Szmirnov-próba

Ebben a szakaszban az ekvidisztáns mintavételezésű adatsorok esetén számolható periodogram kiszámítását végző kódot közöljük. Bemeneti adatként itt is az adatsor x vektorát és az időpontokat tartalmazó t vektort kell megadni. A program kiszámítja a Fisher-féle tesztet és a Kolmogorov– Szmirnov-próbát is.

%%% period.m %%%
% returns the periodogram (II), Fisher's p (P), and the Kolmogorov-Smirnov
% standardized cumulative periodogram (SCP->Y)
function [II,P,Y]=period(x,t)
if size(x,1)==1 x=x'; end;
n = size(x,1);
I = zeros(n,1);
if t==0 t=1/n:1/n:1; end;
for j=1:n;
    for ttt=1:n;
        e(ttt) = exp(-complex(0,1)*ttt*2*pi*j/n);
    end;
    I(j) = 1/n*(abs(e*x))^2;
end;
II = I(1:round(n/2)); % get the important part
% invoking fish.m, and ks.m to calculate Fisher's p, and the SCP
P=fish(I);
Y=ks(x,I);

A következő kódrészlet Fisher tesztje rejtett periodicitások keresésére. Számításához bemeneti paraméterként mindössze az adatsor periodogramjára van szükség.

%%% fish.m %%%
% returns Fisher's p for a given periodogram
function [P]=fish(I)
n=length(I);
q = floor((n-1)/2);
II=I(1:q);
zeta = max(II)/(1/q*sum(II));
test = zeros(1,q+1);
for j=1:(q+1); test(j) = (-1)^(j-1)*choose(q,j-1)*max((1-(j-1)*zeta/q),0)^(q-1); end;
% note that 'choose' is not a MATLAB function. it is a piece of code written
% by Szlávik, Árpád at Ericsson. this is necessarry since MATLAB cannot
% compute this efficiently especially for large values. (the MATLAB function

% 'factorial' is accurate only for n<=21 according to the documentation.)
P = sum(test);
P = 1-P;

Az alábbi program a Kolmogorov–Szmirnov-próbához szükséges sztandardizált kumulatív periodogramot számítja ki az adatsor és periodogramja alapján.

%%% ks.m %%%
% calculates the SCP (Y) from a given dataset and it's periodogram
function Y=ks(x,I)
n=length(x);
q=floor((n-1)/2);
Y=zeros(q-1,1);
s=0;
s_k=0;
sf=1; %starting freq. Actually this is startfreq+1
% note that setting sf to greater values it is possible to see wether low
% or high frequency components dominate the spectrum (we did not do that)
for i=sf:q; s=s+I(i); end; %\sum_{k=1}^q I(\omega_k)
% SCP
for i=sf:(q-1);
    s_k=0;
    for j=sf:i; s_k=s_k+I(j); end;
    Y(i)=s_k/s;
end;

A programban lehetőség van arra, hogy a kezdőfrekvenciát (sf) megahtározzuk. Ekkor a sztenderdizált kumulatív periodogramot csak ennél nagyobb frekvenciákra számítja ki. Ha az adatsorban nem dominánsak a magasabb frekvenciájú komponensek, akkor a kezdőfrekvencia emelésével a grafikon egyre lineárisabb lesz, bizonyos sf fölött már nem lép ki a kritikus tartományból. Ezzel a módszerrel az adatsor legmagasabb frekvenciájú periodikus komponensének frekvenciájára lehet következtetni. Ezt az elképzelést nem ültettük át a gyakorlatba.