Idősorok zajszűrése gyors Fourier transzformációval
Készítette: Jakab Jenő Dátum: 2013.05.27.
Tartalom
Tartalom............................................................................................................................................ 2 1. Bevezető ......................................................................................................................................... 3 2. Matematika alapok ....................................................................................................................... 3 2.1. Alapfeltevések, definíciók ....................................................................................................... 3 2.2. A Fourier - együtthatók meghatározása ................................................................................ 4 2.2.1. a0 meghatározása ............................................................................................................. 4 2.2.2. ak és bk meghatározása .................................................................................................... 4 2.3. Általános, tetszőlegesen hosszú T szerint periodikus függvények vizsgálata ........................ 7 3. Fourier - transzformációban és Inverz Fourier - transzformációban rejlő lehetőségek ............. 9 4. Numerikus megoldás a Duna budapesti vízállás idősorán ........................................................... 9 4.1. Az diszkrét idősor optimalizálása az FFT számára ............................................................. 10 4.2. Az FFT algoritmus bemutatása............................................................................................ 15 4.3. Az FFT (és IFFT) alkalmazása a Duna vízállás idősorán ................................................... 18 5. Összefoglaló ................................................................................................................................. 23
2
1. Bevezető A vízmérnöki gyakorlatban a modern, digitális adatrögzítő mérőeszközök elterjedésével lehetőségünk van pontosabb, átfogóbb elemzések készítésére az időelőnyök növelése mellett. Ugyanakkor fel kell készülni - a kétségbevonhatatlan előnyök ellenére - az új megoldásokban rejlő hibák, hibaforrások szakszerű kezelésére is. Ennek az összefoglalónak éppen ez a célja: egy lehetséges megoldást szeretnénk szélesebb körben kipróbálni és ajánlani az adatregisztrálók hibájának mérséklésre, szűrésére. Ehhez jó eszköz a Fourier - analízis, pontosabban a Fourier - transzformáció, amely a villamosmérnöki – hírközlési, informatikai gyakorlatban elterjedt, rutinszerűen alkalmazott. Konkrét jó példája lehet a Fourier - transzformáció széles körű alkalmazásának a JPEG formátumra kódolt képek, MP3 formátumú zenék, tömörítők, vagy a Photoshop képszerkesztő több funkciója is.
2. Matematika alapok Ebben a vázlatpontban szeretném bemutatni a kollégáknak a későbbi alkalmazáshoz szükséges elméleti alapokat. A fejezet megírásában nagyrészt a „Matematika A2” egyetemi segédanyagra (BME, Sándor Csaba) támaszkodtam.
2.1. Alapfeltevések, definíciók Első lépésben vizsgáljuk meg a klasszikus, 2π szerint periodikus függvényeket. Ezekről a periodikus függvényekről feltesszük, hogy előállíthatók elemi 2π szerint periodikus függvények összegeként. Ilyen elemi függvények a szinusz és koszinusz szögfüggvények. Tehát a fenti állításunk (bizonyítás nélkül) a matematika nyelvén így néz ki: →
( )=
(
∙ cos( ∙ ) +
∙ sin( ∙ )),
( . )
ahol: - ak, bk amplitúdók (függőleges y tengely szerint nyújtja – zsugorítja, transzformálja az elemi függvényt). - k körfrekvencia, amely ebben az esetben egyenlő a frekvenciával. (vízszintes x tengely szerint transzformálja a függvényt). Esetünkben a „k”-nak egész számnak kell lenni, mert ha nem lenne egész, akkor már nem lenne 2π szerint periodikus. Tehát így belőle nem állítható elő pontosan akármelyik 2π szerint periodikus f(x) függvény sem. k=0 esetben a sin(0·x)=0, ezért az a0 tagot kiemelve általában k=1-től szokás felírni az összegképlet alakot. Jogosan felmerülhet (1.) alapján, hogy miért két fajta és miért ezt a két fajta elemi szögfüggvényt használjuk f(x) közelítésére. Bizonyítás nélkül fogadjuk el ebben az esetben, hogy a két elemi függvény független egymástól (geometriailag ez azt jelenti, hogy merőlegesek egymásra 1. ábra, vagy másként értelmezve a függetlenséget, skaláris szorzatuk 0) és a kettejük lineáris kombinációja a polinomokhoz hasonló univerzális összegző képességekkel (univerzális approximátor) rendelkezik. Érdekesség még, hogy a szinusz és koszinusz páratlan, illetve páros függvény (páratlan: sin(-x)=3
sin(x) => origóra középpontosan tükrös, páros: cos(-x)=cos(x) => y tengelyre tükrös). Magyarán azt a sejthető megállapítást fogalmaztuk meg, hogy tetszőleges periodikus függvények páros és páratlan elemi függvények összegeként előállíthatók. Az előző gondolatot másként megfogalmazva megvizsgálható, hogy adott periodikus függvény milyen mértékben páros és páratlan. Tovább fűzve a gondolatot logikailag előrevetíthető, hogy páros függvények Fourier - együtthatói csak a koszinuszos tagokból (ak) , páratlan függvények pedig csak szinuszosos (bk) tagokból állnak.
2.2. A Fourier - együtthatók meghatározása Az alcímben megfogalmazott cél közvetlenül következik a 2.1. pontban megfogalmazottakból: ha ismerjük az ak, bk együtthatókat elkészültünk a feladattal!
2.2.1. a0 meghatározása Tegyük fel, hogy a 2π szerint periodikus f(x) Riemann – integrálható a [0,2π] intervallumban, tehát: →
( )
=
(
+
∙ cos( ∙ ) +
(
∙ sin( ∙ )))
( . )
Az integrálást elvégezve az alábbi egyenletet kapjuk: →
( )
= [
∙ +
(
∙ sin( ∙ )
+
∙ cos( ∙ )
)] ( . )
Ha a (3.) -at megvizsgálva láthatjuk, hogy egyetlen nem 0 tagunk marad, mivel a szinuszos és a koszinuszos tagokat saját periódusukon, vagy azok többszörösén integráljuk. Tehát az egyenlet (4.) -re egyszerűsödik: ( )
= [
∙ ]
=
∙ 2 ( . )
Ebből már a0 kifejezhető: =
1 2
( )
( . )
2.2.2. ak és bk meghatározása Az ak és bk együtthatók számításához először nézzünk meg egy-két speciális integrált. Az integrálok argumentumaihoz mankóként alkalmazható a trigonometrikus függvények addíciós összefüggései, amelyeket az 1. ábrán grafikusan is nyomon követhetünk.
4
1. ábra trigonometrikus függvények addíciós összefüggései
Speciális – alkalmazott integrálok: k≠l esetben:
cos( ∙ ) ∙ cos( ∙ )
=
1 sin ( + ) ∙ 2 +
+
1 sin ( − ) ∙ 2 −
−
1 = 2
cos ( + ) 1 ∙ [− 2 +
= 0 ( . )
cos(( − ) ) − cos ( + )
sin ( + ) +
sin( ∙ ) ∙ cos( ∙ ) =
cos(( + ) ) + cos(( − ) )
sin ( − ) −
sin( ∙ ) ∙ sin( ∙ )
=
1 = 2
−
1 = 2
= 0 ( . )
sin(( + ) ) + sin(( − ) )
cos(( − ) ) ] −
= 0 ( . )
Látható, hogy a l => k-hoz speciális eset. A fenti képletekből nem derül ki, hogy mekkorák k=l esetén a fenti határozott integrálok!
5
k=l esetben: ( ∙ )
( ∙ )
1 + cos (2 ∙ 2
=
∙ )
1 − cos (2 ∙ 2
=
sin( ∙ ) ∙ cos( ∙ )
=
∙ )
=
sin (2 ∙ ∙ ) 2
=
1 sin(2 ∙ ∙ + 2 4∙
∙ )
1 sin(2 ∙ ∙ − 2 4∙
∙ )
cos(2 ∙ 4∙
∙ )
= −
= ( . )
= (
. )
= 0 (
. )
ak meghatározása: →
( ) ∙ cos (
∙ )
=
(
+
(
∙ cos( ∙ ) +
∙ sin( ∙ ))) ∙ cos (
∙ )
→
∙ cos(
= =
∙ ,
∙ )
+
(
∈ (
(
∙ cos( ∙ ) +
∙ sin( ∙ ))) ∙ cos (
)
. )
A (12.) képletnél (6.), (8.), (9.), (11.) felhasználva látható, hogy csak k=m esetnek van 0-tól különböző értéke az integrálnak. Tehát ak kifejezhető végül a következő alakban: =
1
( ) ∙ cos ( ∙ )
(
. )
bk meghatározása: →
( ) ∙ sin (
∙ )
=
(
+
(
∙ cos( ∙ ) +
∙ sin( ∙ ))) ∙ sin (
∙ )
→
∙ sin(
= =
∙
,
∙ )
+
(
∈ (
(
∙ cos( ∙ ) +
∙ sin( ∙ ))) ∙ sin (
)
. )
A (14.) képletnél (7.), (8.), (10.), (11.) felhasználva látható, hogy csak k=m esetnek van 0-tól különböző értéke az integrálnak. Tehát bk kifejezhető végül a következő alakban: =
1
( ) ∙ sin ( ∙ )
6
(
. )
2.3. Általános, tetszőlegesen hosszú T szerint periodikus függvények vizsgálata Ha a 2.1. és 2.2. pontokban felsorolt feltételek igazak, akkor az ott meghatározott állítások következményei kisebb formai módosításokkal tetszőlegesen hosszú T szerint periodikus függvényekre is alkalmazhatóak: Vegyük most a T szerint periodikus g(x) függvényt: →
( )=
(
∙ cos
2∙
∙
∙
+
∙ sin
2∙
∙
∙
),
(
.)
Vegyük észre, hogy a szinusz és koszinusz függvények argumentumában a fizikai tanulmányokból ismert körfrekvenciát figyelhetjük meg: ( )=
2∙
∙
=2∙
∙ ∙ (
.)
ahol, - T: periódusidő, de gyakran alkalmazzák formailag a periódusidő reciprokát: a frekvenciát is.
2. ábra Fourier Transzformált képzése g(x) függvényen ω1 körfrekvenciánál
7
3. ábra Fourier Transzformált képzése g(x) függvényen ω2 körfrekvenciánál
Ebben az esetben az elemi függvénykomponensek együtthatói (az amplitúdók) az alábbiak szerint alakulnak:
=
=
=
1
( )
(
. )
2
2∙ ( ) ∙ cos (
∙
∙ )
(
. )
2
2∙ ( ) ∙ sin (
∙
∙ )
(
. )
A 2. – 3. ábrákon grafikusan is megfigyelhetjük, hogyan kell értelmezni a Fourier transzformációt. Az ak együtthatókat a grafikonon jelölt zölddel satírozott előjeles területek összegeiként kapjuk, elosztva a periódusidő kétszeresével. A bk együtthatóknál pedig a pirossal satírozott területekkel kell hasonlóan eljárnunk.
8
3. Fourier - transzformációban és Inverz Fourier - transzformációban rejlő lehetőségek
A 2. pontban leírtak szerint látható, hogy bizonyos feltételek mellett függvények előállíthatók Fourier – sorba fejtéssel is tetszőleges pontossággal. Ilyen diszkrét időpontokban értelmezett függvény lehet akár egy vízállás idősor is. Tehát lehetőségünk van például a vízállás idősorainkat különböző frekvenciájú egységnyi harmonikus függvények különböző arányú összegeiként is értelmezni. Munkám elején élek azzal a feltételezéssel, hogy általánosságban a gyakorlatban előforduló vízmozgások valós változása lassabb, mint a mérőműszerek által belekevert zaj irányváltozása. Tehát a vízállás idősor frekvenciák szerinti un. spektrális felbontásánál (ez lényegében ak(k) és bk(k) függvények) van esély különböző, domináns frekvenciatartományok elkülönítésére, eredetüknek beazonosítására. Az ak(k) és bk(k) függvények meghatározását az irodalomban Fourier – transzformációnak hívják. Az integrálások során ak(k) és bk(k) függvények csak a k-tól (a frekvenciától) fognak függeni adott g(x) esetén, x-től (x az idősornál a független változó: az idő) nem. A frekvenciatérben tehát bizonyos frekvenciájú komponensek törlése lehetővé teszi a törölt frekvenciájú komponens törlését az idősorból is az inverz transzformáció után. Az inverz transzformáció az a számítási lépcső, amely során összegalakban ak(k) -t és bk(k) -t felhasználva előállítjuk g(x) -et, esetünkben az idősort. Elméletileg tehát a mérőműszer zajának spektrumát ismerve - kicsiny (cm-en belüli) hibával terhelten - jó esélyünk van idősorainkat elfogadható pontossággal zajmentesíteni.
4. Numerikus megoldás a Duna budapesti vízállás idősorán
A vízmérnöki gyakorlatban fontos állapotváltozók (pl.: vízszintek, vízhozamok, stb.) a mai kor méréstechnikájával gyakorlatilag tetszőlegesen sűrű időközönként megmérhetőek. Fontos - hogy céljaink szerint - méréseink információtartalma, számításigénye gazdaságos (optimális) legyen. Ez azt jelenti - hogy kedvezőtlen esetekben - a túl ritka mérési adatok között kénytelenek vagyunk az interpoláció bizonytalanságait felhasználni, valamint azt, hogy a túl sűrű adathalmazon kénytelenek vagyunk plusz erőforrás bevonással különböző elvű ritkításokat végrehajtani. Munkám további szakaszában élek azzal a feltételezéssel, hogy a mérőrendszerek általánosságban optimálisan üzemelnek. Tehát nem a rendszer felülbírálatára törekszem, hanem a rendszerben található adatsorok: idősorok optimális kihasználására. A vízállás idősorokat vizsgálva a Fourier – transzformációs zajszűrők gyakorlati megvalósításánál a gyors számítás és az adatszerkezet jellege miatt (diszkrét jellgű) a numerikus megoldók közül választottam. A lehetőségeket figyelembe véve a választott numerikus módszerem a gyors Fourier – transzformáció (angol rövidítés: FFT, visszalakítás inverz művelet: IFFT) lett. A módszer kiválasztásánál szempont volt, hogy az Országos Vízügyi Főigazgatóságban (OVF) homogén felbontású idősorok készítésére is van igény egyfajta feldolgozott késztermékként. 9
4. ábra Vízállás idősor a Dunán Budapestnél
Az elemzett idősorom - amelyen a módszertant ismertetem - a Duna budapesti vízmércéjén mért és a MAHAB-ban tárolt vízállás „hidrológiai idősor” a 2013.01.01. 7:00 – 2013.03.27 15:01 időszakra vonatkozóan. A későbbiekben bevezetésre kerülő relatív idő fogalmát az időintervallum kezdő időpontjához igazítom a rendszer alapegysége pedig az 1 perc lesz.
4.1. Az diszkrét idősor optimalizálása az FFT számára
Az FFT algoritmusnál fontos az idősort homogén felbontással megadni. Ez azt jelenti, hogy két szomszédos vízállás érték között dT idő kell, hogy elteljen egységesen. Az FFT számítási algoritmust részletesen a 4.2 fejezetben ismertetem. A kívánatos homogén időbeli felbontás viszont gyakorlati és technikai okokból kifolyólag is gyakran nem áll rendelkezésre. Tehát első lépésben a homogén felbontású idősor előállítását kell megoldani valamilyen (esetünkben lineáris) interpolációs eljárással. A vizsgált Duna vízállás idősor kb. órás időközű feldolgozatlanul. A homogén felbontást 1 perces finomsággal határoztam meg első lépésként. Az így létrejött idősorban viszont az interpolált vízszintek aránya túl nagy a mérésekhez képest és az adatsor mérete is jelentősen megnövekedett. Ez azt jelenti, hogy hígult az információ-sűrűség. A percnél durvább felbontás lehetőségét tehát az interpoláció fajlagos hibája és a számítási igény csökkentése mentén célszerű megvizsgálni. A kellően nagy információ-sűrűséggel rendelkező kívánatos adatstruktúra eléréséhez egy optimalizációs eljárást javasolok, amelynek szempontjai egyidejűleg: -
I . A lehető legdurvább időbeli felbontást kell alkalmazni
-
II. A lehető legtöbb mért adatot kell alkalmazni
-
III. A lehető legkevesebb interpolálat adatot kell felhasználni. (A II. és III. -ból levezethető, hogy legyen a lehető legnagyobb legyen a mért és interpolált adatok aránya!) 10
I-III. vázlatpontokat megfontolva tehát kaphatunk egy két szabadságfokú rendszert, amely változói az idő felbontása és a kezdő időpont (vagy az időrács) fázisa. A második szabadságfok azt jelenti, hogy ha van például homogén 2 perc felbontású mérésünk 0 hasznosul, ha a 2 perces időfelbontás a 2 perces mérések közé esik, de ha egybeesik a mérés a felbontás fázisával, akkor mindet hasznosítani tudjuk. A két szabadságfokú rendszer kimenete pedig legyen a felhasznált mérések aránya az összes méréshez (eredeti adat) képest! Ezt az arányszámot a továbbiakban hatásfoknak hívom. (
, )=
(
, )
(
.)
é
ahol - η: Hatásfok [%] - dT: Felbontás [min] - φ: Fázis [min] - Nmért: Megmért vízállások darabszáma az adott időintervallumban [-] - Nfelhaszn: Az interpolációval homogenizált idősorban a felbontás és a fázis függvényében a felhasznált eredeti mérések darabszáma [-]
5. ábra Hatásfok az időbeli felbontás és a fázis szerint
11
6. ábra Hatásfok fázis szerint maximált ábrája az időbeli felbontás szerint
Az 5. ábrán önkényesen az 1 – 300 perc közötti felbontásokkal vizsgálom meg a hatásfokokat 0 – 100 percig terjedő fázisok mellett. Az ábrán jól látható, hogy 1 perces felbontásnál minden megmért adat hasznosul, majd a felbontás csökkenésével a felbontás szerinti hatásfok trendje csökken. A fázis – hatásfok viszonya pedig jól mutatja, hogy kb. 60 perces periódusokban mért a regisztráló. A 6. ábrán az 5. ábrán bemutatott függvény egy speciális módosítását láthatjuk, ahol az egyes felbontások függvényében ábrázoljuk a hatásfokot úgy, hogy a fázis szerinti a maximális hatásfokot vesszük mindegyik felbontásnál. Ezek alapján a fent megfogalmazott (I. - III.) optimalizáció szempontjai szerint a 60 perces felbontás kedvezőnek tűnik. Annak érdekében, hogy az I. vázlatpontot megfelelően érvényre tudjuk juttatni egy súlyfüggvény és egy alapérték tényező bevezetését, illetve alkalmazását tartom kívánatosnak. Ez a kedvező időbeli felbontást felülreprezentáló súlyfüggvény első javaslatom szerint legyen a szinusz hullám első negyede: (
) = sin (
0.5 ∗ −
∙(
−
)) (
.)
ahol - S: Súly [-] - dT: időbeli felbontás [min] - dTmax: A választott legnagyobb alkalmazott felbontás, amely konstans [min] - dTmin: A választott legkisebb alkalmazott felbontás, amely konstans [min] 12
7. ábra Súlyfüggvény az időbeli felbontás függvényében
Tehát a kedvező megoldást a 6. ábrán már szemléltetett fázis szerint maximalizált hatásfok függvény súlyfüggvénnyel vett szorzatából és az alapérték tényező összegétől várom az alábbi alak szerint (a kapott függvényt fittségi (fitt, mint kedvező) függvénynek hívom): (
)
( → )
= ( ) ∙ ( ) ( → = ( ( )+ )∙ ( )
)
+ ∙ ( ( → ) (
) .)
( → )
ahol - ϑ: fittség. A függvény csak az időbeli felbontástól függ, olyan módon, hogy egy felbontáshoz a fázis szerint már maximalizált hatásfokot rendeljük. [-] - S: a súly [-] - η: a hatásfok [%] - λ: alapérték tényező, amely konstans. Választott értéke: λ=0.2 Szerepe: egyfajta kezdeti súly biztosítása. A 8. ábrán megszemlélhetjük a fittségi függvényt, amelynek maximuma van a kívánatos 60 perces felbontásnál. Ezek után a 60 perces felbontást kötött változóként kezelve a hatásfok függvényből megállapíthatjuk, hogy mekkora fázisnál érte el a fittségi függvény a maximát (a fittség és a hatásfok függvénynek a fázis szerint ugyanott van a maximuma).
13
A számítások menete: max
(
) á
max (
( → )
á
− → á
,
á
− →
)− →
,
( . ) ( . )
á
(
. )
8. ábra Fittség ábra az Időbeli felbontás szerint és az időbeli felbontás fázisa szerint maximalizálva
Az optimalizáció végeredményként a homogén időbeli felbontás 60 percre a kezdeti fázis 1 percre adódott. Ilyen paraméterek mellett a konkrét mérések 58.92% -át tudtam hasznosítani. Összességében a fent javasolt módszerrel egy jó becslés adható a két szabadságfokú rendszer paramétereire. Természetesen az olvasóban felmerülhet az a jogos kérdés, hogy miért pont ezt a formát, súlyfüggvényt és alapérték tényezőt alkalmaztam? Összességében érvként a mérnöki becslést tudom említeni, az eljárás alkalmazhatóságát, a programmodul finomítását a gyakorlat fogja meghatározni. A másik lényeges kérdés az, hogy a homogén időbeli felbontás mennyire jól alkalmazható? Egyfajta felbontás jól jellemzi-e a teljes vizsgált idősort? Ha csak arra gondolunk, hogy árvizes időszakokban a megnő az adatregisztrálás egy időegységben, célszerű lehet a teljes idősort különböző homogén felbontású rész idősorok sorozatásból előállítani és vizsgálni, esetenként a felbontás fázisát a szélső értékekhez igazítva. Ennek a problémának az orvoslására jó kezdő lépés lehet a mérések időközeinek változását megvizsgálni az idő függvényében.
14
4.2. Az FFT algoritmus bemutatása
Az FFT az alkalmazott homogén felbontású adatsoroknál az alábbi képlet alapján számítható:
( )=
∙
2
⎛ ∙⎜ ⎜
( ) ∙ cos
( ) ∙ sin
2∙
2∙
∙ ( − 1) ∙ ( − 1)
∙ ( − 1) ∙ ( − 1)
⎝
⎞ ⎟ ( ⎟
.)
⎠
ahol, - A: az amplitúdó vektor, amelynek első eleme ak, második bk [cm] - N: adatsor hossza, elemszáma [-] - x(j): j -dik mérés (ez lehet a diszkrét vízállás). j=1…N –ig változik [cm] -
∙
∙ ( − 1): Ez a körfrekvencia, amely k szerint változik. ={
0.5 , ℎ = 1 1 , ℎ > 1
Érdemes az amplitúdókat is k = 1…N -ig vizsgálni, így pont olyan hullámokra végezhetjük el a transzformációt, amelyek hullámhossza homogén módon sűrűsödve 0…N-1 –ig fedi le a vizsgált mintát (az első tagot az átlagból és a 0 -ból összetevődő vektor). A k > N fölötti hullámokat nem érdemes megvizsgálni, mert azokon a frekvenciatartományokon nincs megfelelő mintánk, ha pedig erősen durvítani akarjuk a transzformációt (k=1…M, M<
9. ábra FFT folyamatábra k=0 esetben cos –os összetevők és mérések 15
10. ábra FFT folyamatábra k=0 esetben sin –os összetevők és mérések
11. ábra FFT folyamatábra k=1 esetben cos –os összetevők és mérések
12. ábra FFT folyamatábra k=1 esetben sin –os összetevők és mérések
16
13. ábra FFT folyamatábra k=N esetben cos –os összetevők és mérések
14. ábra FFT folyamatábra k=N esetben sin –os összetevők és mérések
17
4.3. Az FFT (és IFFT) alkalmazása a Duna vízállás idősorán
Az FFT eljárás elvégezésével tehát előáll a jelösszetevők körfrekvencia szerinti függvénye a spektrum, amely a 15. ábrán megszemlélhető! Az ábrán kék szín jelöli a koszinuszos, piros a szinuszos függvényeket, az abszolút érték függvény pedig a zöld színű. Az abszolút érték függvényből jól látható, hogy a legnagyobb amplitúdójú jel 680 cm körüli, ellenőrzésként pedig megnézhető a 4. ábra maximális vízállása, amely kb. 580 cm. Ez első ránézésre reális lehet, az eltérés oka pedig a sok különböző jelkomponens összegzett hatása. | ( )| =
( ) +
( ) (
.)
ahol - |A(k)|: az amplitúdó nagysága [cm]
15. ábra A Duna budapesti 2013.01.01 7:00 - 2013.03.27 15:01 mért vízállás idősorának spektruma
A spektrum előállítása után a feladat a zavaró frekvenciák, frekvenciatartományok beazonosítása és törlése a spektrumból! A zaj periódusidejére, frekvenciájára az alábbiakat javasolom: =
∙
=
2048 ∙ 60 205
= 599.4
≈ 10ℎ →
=
1
1 = 0.1 ( ℎ
.)
ahol - T: a periódusidő [h]. k alapján a 10 óránál kisebb periódusidejű hullámokat törlöm a spektrumból! - N: homogén felbontású idősor elemszáma - kcut: az a körfrekvencia változó, amelynél nagyobb értékeket már törlök az adatsorból. k=205 –ös érték már a zaj nagyságrendjébe eső körfrekvenciát eredményez véleményem szerint, vagy ha nem is a 18
zaj frekvenciája a műszerzajból keletkező hamis frekvenciaként tekintek rá (Hamis frekvenciák akkor keletkezhetnek, ha nincs elég, legalább 2 mérésünk a frekvencia periódusideje alatt). - dTalk: az alkalmazott időfelbontás [min] A 16. ábrán a kis frekvenciájú összetevők láthatóak, az idősor jórészt a k=1…205 terjedő spektrumrészletből épül fel! Szemmel is jól látható, hogy a k>205 tartomány kicsiny arányba járul hozzá az idősor alakulásához. Ebben a frekvenciatartományban 1 cm–nél kisebb hullámok keletkeznek, amely esetünkben a 17. ábra alapján a zaj amplitúdójaként is azonosíthatók. Tehát ha töröljük a k > 205 spektrumot és visszaalakítjuk a módosított spektrumot az időtérbe, akkor egy jobb, a valóságot homogénebben jellemző idősort kaphatunk vissza. A spektrumból történő törlésnél viszont óvatosan kell eljárnunk, mert ha töröljük a k > 205 frekvenciasávot rossz eredményt kapunk. Ennek oka nem fizikai, hanem numerikus. A magyarázatához nézzük meg a 13. és 14. ábrákat, ahol látható az ábra elején, hogy ha a hullámhossz közel esik az időbeli felbontáshoz, akkor a diszkrét megoldásnál nem érvényesül a hullám teljes hatása, hanem csak az aktuális diszkrét érték. Ez felfogható úgy, mint egyfajta „numerikus lebegés (rezonancia)” a fizikából vett analógiával élve. Tehát ez a jelenség mindig az diszkrét alapon számolt spektrum végén jelentkezik az FFT-nél.
16. ábra A Duna budapesti 2013.01.01 7:00 - 2013.03.27 15:01 mért vízállás spektrumának részlete
A numerikus lebegéssel befolyásolt tartomány ugyanakkor nem tekinthető valós hibának, a jel előállításában fontos szerepet játszik, mert ugyan hamis frekvenciaként, de a kisebb frekvenciák hatását hordozza magában! Tehát a törlésre ajánlott tartomány: ( )
=
…
−
,
=
ezek alapján: ( )
= 0 (
19
.)
(
.)
ahol - I(k)clear: törlésre kijelölt frekvencia változók tartománya - klebeg: a numerikus lebegés által érintett frekvencia változó - kend: a legnagyobb értékű frekvencia változó. kend=N Ezek után alakítsuk át a módosított spektrumot az idő térbe az inverz transzformációval (IFFT) az alábbiak szerint:
( )=
⎛ ⎜ ⎜
( ) ∙ cos ( ) ∙ sin
2∙
2∙
⎝
∙ ( − 1) ∙ ( − 1)
⎞ ⎟ ( ⎟
.)
∙ ( − 1) ∙ ( − 1) ⎠
A 17. ábrán megszemlélhetjük a mért és a zajszűrt Duna budapesti mért és szűrt vízállás idősorokat. Ezek alapján szemre jó egyezés tapasztalható a két idősor között.
17. ábra A Duna budapesti 2013.01.01 7:00 - 2013.03.27 15:01 mért (piros) és zajszűrt (zöld) vízállásai
A 18. és 19. ábrákon a 2013-as év eleji időszakok budapesti mércén mért Duna árhullám csúcsai láthatóak, illetve a 20. ábrán a vizsgált időszak egyik jellemző kisvizes vízállása. Ezek alapján megfigyelhető, hogy a szélsőértékeknél jó – cm-belüli hibát okoz a szűrés. A szélsőértékeken kívül vizsgálandóak lehetnek az erőteljes áradó és apadó szakaszai a vízállás idősornak, de ezeken a szakaszokon sem tapasztaltam cm-es hibahatárt túllépő számított vízállás értéket.
20
18. ábra Árhullámcsúcs_1 a vízállás idősorban piros a mért, zöld a szűrt
19. ábra Árhullámcsúcs_2 a vízállás idősorban piros a mért, zöld a szűrt
20. ábra Kisvíz a vízállás idősorban piros a mért, zöld a szűrt
21
21. ábra Vízállás idősor eleje piros a mért, zöld a szűrt
22. ábra Vízállás idősor vége piros a mért, zöld a szűrt
Jelentős, cm-es hibatartományt meghaladó hibát az idősor elején és a végén találhatunk, amelyek a 21. és a 22. ábrákon is szerepelnek. A közelítő módszerből fakadó, az intervallum szélein megfigyelhető hullámosodást a matematikában Runge - jelenségnek hívják, amelynek hatását a zaj visszakeveréssel lehetne csökkenteni. Ha meggondoljuk viszont, hogy a teljes három hónapból a Runge - jelenség által jelentősebben befolyásolt időszak az intervallum elején és végén 2-2 nap, akkor látható, hogy ebben az esetben az idősor kb. 95.5% -án kapunk elfogadhatóan jó közelítést. Érdemes azt is megjegyezni, hogy hosszú idősoroknál a Runge - jelenség által meghatározott szakaszok arányaiban rövidülnek a teljes adatsor hosszához képest, tehát az FFT-t különösen hosszú idősoroknál lehet jó hatásfokkal alkalmazni.
22
5. Összefoglaló
Munkám eredményeként bizonyítható, hogy a Fourier - transzformáción alapuló FFT-s zajszűrés sikeresen alkalmazható vízállás idősorok műszerhibáinak csökkentésére. Eljárásom sikeres gyakorlati alkalmazásának és főként automatizálásának egyik sarokfeltétele lehet az általam definiált fittségi függvény paramétereinek általános érvényűsége. Az egyszerű alkalmazáshoz tehát elengedhetetlen egy teszt program készítése, amely szélesebb körű használatának tapasztalatai fontosak a paraméterhatárok tisztázásában, illetve az egyes modulok felhasználóbarát kialakításában. A dolgozatban látható grafikonok, eljárások nagyrészt Matlab környezetben valósultak meg, ahol egyrészt a Matlab belső függvényeit, másrészt saját fejlesztésű kódokat alkalmaztam. (Az alkalmazott kódokat szükség esetén az érdeklődők rendelkezésére tudom bocsájtani). Véleményem szerint a bemutatott eljárás teljes körűen megvalósítható önálló programként, illetve az OVF -ben található adatbázisokra épülő programok alkalmazásaként is. Másik fontos fejlesztési irányként az idősorok felbontásának inhomogenizálását tartom az FFT számára, úgy, hogy a teljes inhomogén idősor homogén részidősorok sorozatából épüljön fel. Természetesen minden numerikus eljárásnak a feladathoz, a kitűzött célokhoz, illetve a rendelkezésre álló adatokhoz kell igazodnia, tehát elképzelhetőek egyéb komolyabb numerikus megoldások fejlesztései is. Tudomásom szerint a hazai vízmérnöki gyakorlatban Fourier - transzformáción alapuló zajszűrést a Budapesti Műszaki és Gazdaságtudományi Egyetem Vízépítési és Vízgazdálkodási Tanszékén is alkalmaznak. A tanszék ilyen irányú tapasztalatai értékesek lehetnek az OVF -s fejlesztésekben is. A műszerhibák szűrésén kívül a bemutatott eljárás jó eszköze lehet a hidrológiai - hidraulikai analízisnek is. A Balaton vízállás idősorainál tesztelteltem a különböző lengésidejű hullámokhoz rendelt aluláteresztő (egy választott küszöbfrekvenciánál nagyobb komponenst töröl az idősorból) szűrőt, amely segítségével a nagyfrekvenciás hullámzás leválaszthatóvá vált az megmért idősorról, így jobban megvizsgálhatók lettek a tólengések a kis frekvenciájú vízállás komponensei. De a mérnöki feladatok sokrétűségét kiszolgálva az FFT átalakítható felüláteresztő (egy választott küszöbfrekvenciánál kisebb komponenst töröl az idősorból) szűrővé, vagy éppen sávszűrővé (két választott küszöbfrekvenciánál közé eső jeleket hagy meg az idősorból) is.
23
