´ METODY PRO ME ˇ REN ˇ ´I VYBRANE ˇ EN ˇ ´I MEZI SIGNALY ´ ZPOZD EEG Ondˇrej Drbal K13131 - Katedra teorie obvod˚ u, Fakulta elektrotechnick´a ˇ e vysok´e uˇcen´ı technick´e v Praze Cesk´
Abstrakt V dneˇsn´ı dobˇe je vyuˇz´ıv´ ana v´ ypoˇcetn´ı technika prakticky ve vˇsech moˇzn´ ych oborech vˇedy a techniky, neurologie nevyj´ımaje. Poˇc´ıtaˇcem prov´ adˇen´a anal´ yza EEG z´aznamu (elektroencefalograf) velmi pˇrisp´ıv´a ke zpˇresnˇen´ı diagn´ ozy pacienta a zjednoduˇsen´ı pr´ace l´ekaˇr˚ u. Poˇc´ıtaˇc pˇri anal´ yz´ach m˚ uˇze napˇr. upozornit jen na rizikov´e u ´seky EEG a l´ekaˇr se tak nemus´ı prokous´ avat nˇekolikahodinov´ ymi z´aznamy. S dneˇsn´ım obrovsk´ ym v´ ypoˇcetn´ım potenci´alem poˇc´ıtaˇc˚ u maj´ı lid´e moˇznost archivovat velk´a mnoˇzstv´ı dat, navz´ajem je sd´ılet a d´ıky internetu je mohou analyzovat i l´ekaˇri doslova na druh´em konci svˇeta. Z toho pramen´ı velk´ a potˇreba spolupr´ ace odborn´ık˚ u r˚ uzn´ ych oblast´ı jako je medic´ına a informatika, kteˇr´ı se spoleˇcnˇe pod´ılej´ı na v´ yvoji program˚ u, anal´ yz a metod zpracov´an´ı biologick´ ych sign´al˚ u.
1
´ Uvod
V t´eto pr´aci se zab´ yv´ am porovn´ an´ım algoritm˚ u pro mˇeˇren´ı ˇcasov´eho zpoˇzdˇen´ı mezi sign´aly EEG. Sign´aly EEG jsou sn´ım´ any nejˇcastˇeji z povrchu skalpu a to 19 a 64 elektrodami, kter´e jsou uspoˇr´ad´any dle urˇcit´e konvence (nejˇcastˇeji uˇz´ıvan´ a je tzv. Montrealsk´ a konvence). Uspoˇr´ad´an´ı tvoˇr´ı p´arov´e a nep´arov´e sondy. Pomoc´ı elektroencefalografu z´ısk´ ame soubor o 19 aˇz 64 kan´alech sn´ıman´ ych souˇcasnˇe a pˇri znalosti p´arov´ ych kan´ al˚ u m˚ uˇzeme pomoc´ı anal´ yzy zpoˇzdˇen´ı mezi jednotliv´ ymi p´ary vytvoˇrit soubor ˇcasov´eho v´ yvoje zpoˇzdˇen´ı. Tato informace se jiˇz d´ ale d´a vyuˇz´ıt pro sledov´an´ı v´ yvoje epileptick´eho z´achvatu. Dalˇs´ı moˇznost´ı jak vyuˇz´ıt metody pro mˇeˇren´ı zpoˇzdˇen´ı mezi sign´aly EEG, je pro vˇcasnou detekci epileptick´eho z´ achvatu. Pokud se nach´ az´ı mozek v norm´aln´ım stavu, je zpoˇzdˇen´ı prakticky nemˇeˇriteln´e. Ovˇsem pˇri pˇr´ıchodu z´ achvatu se zasaˇzen´e ˇc´asti mozku sesynchronizuj´ı a zpoˇzdˇen´ı zaˇc´ın´a b´ yt dobˇre mˇeˇriteln´e. Probl´em vˇsak m˚ uˇze nastat, pokud i pˇri z´achvatu bude zpoˇzdˇen´ı nulov´e, nebo pˇri velmi kr´atk´ ych (interikti´aln´ıch) z´ achvatech.
1.1
P´ asma EEG sign´ alu
Elektrick´e charakteristiky mozkov´ ych vln rozdˇeluj´ı stavy vˇedom´ı do ˇctyˇr z´akladn´ıch hladin. Pro pojem hladiny se pouˇz´ıv´ a tak´e pojmu frekvenˇcn´ı p´asma nebo rytmy [3]. P´ asmo delta - pˇri tomto stavu jsou v´ yraznˇe utlumeny vˇsechny ˇzivotn´ı funkce. Nal´ez´ame se v nˇem bˇehem bezesn´eho sp´ anku, ale i bˇehem bezvˇedom´ı, zp˚ usoben´eho nemoc´ı ˇci u ´razem. P´ asmo theta - v tomto rytmu jde opˇet o v´ yrazn´ yu ´tlum funkc´ı. Nal´ez´ame se v nˇem pˇri ospalosti, us´ın´an´ı, uvolnˇen´ı, ale i pˇri lehk´e mozkov´e dysfunkci. P´ asmo alfa - v tomto rytmu tˇelo a mysl nereaguje na vnˇejˇs´ı podnˇety, ale ˇclovˇek je pˇri pln´em vˇedom´ı, jde tedy o klid a relaxaci. ˇ ım P´ asmo beta - v tomto stavu je ˇclovˇek zcela pˇri vˇedom´ı a je pˇripraven reagovat na vnˇejˇs´ı podnˇety. C´ vyˇsˇs´ı je dominantn´ı frekvence, t´ım v´ıce je ˇclovˇek podr´aˇzdˇen, ve stavu u ´zkosti, apod.
1.2
EEG a epilepsie
Sign´aly EEG zdrav´eho ˇclovˇeka maj´ı charakter ˇsumu (obr´ azek 1(a)), pˇriˇcemˇz nen´ı pozorovateln´a velk´a aktivita v p´asmu delta. Epileptick´ ych z´ achvat˚ u existuje nˇekolik druh˚ u, ale vˇetˇsinou je z´achvat charakteristick´ y tzv. komplexem hrot–vlna. Velk´ y n´ astup hrotu je zp˚ usoben souˇcasn´ ym v´ ybojem velk´eho mnoˇzstv´ı neuron˚ u. Tyto neurony jsou vlivem epilepsie sesynchronizov´any a pˇri takov´em z´achvatu nejsou schopny v mozku pˇren´ aˇset d˚ uleˇzit´e informace. Sign´al pˇri z´achvatu je na obr´azku 1(b). Je vidˇet, ˇze pˇri takov´em pr˚ ubˇehu je zv´ yˇsen´ a aktivita v p´ asmu delta, coˇz tak´e signalizuje poˇskozen´ı mozku. Z obr´azku je tak´e vidˇet, ˇze jsou zdeformovan´e i sign´ aly v ostatn´ıch p´ asmech. Jak je tato deformace z´avadn´a je jasn´e z popisu p´asem v odstavci 1.1. EEG a jeho základní rytmy
EEG a jeho základní rytmy
1000
10000
500
5000
0 0
−500 −1000
0
0.5
1 delta
1.5
2
2.5
1000 500 0 −500
0
1
2
3
4
5
3
3.5
4 theta
1
2 delta
6
1000 0
−100
−1000
−1000
−200
50 0
−100
−50 5
5
0
0
4
4
1000
0
1
2
3
4
5
−2000
0
2
4
beta
3
3
0
100
2
0
100
100
1
−5000
2000
alfa
0
5
2000
200
−200
4.5
200
−100
6
8
10
−2000
0
7
2
8 theta
4
alfa
9
10
6
8
10
6
8
10
beta
4000
2000 1000
2000
0 0
0
1
t [s]
2
3
4
t [s]
5
−2000
−1000 0
2
4
6
8
10
−2000
0
2
t [s]
(a) Norm´ aln´ı EEG
4 t [s]
(b) EEG pˇri z´ achvatu
Obr´ azek 1: EEG sign´al a jeho p´asma
2
Metody
Zd´anlivˇe jednoduch´ a u ´loha zjiˇstˇen´ı zpoˇzdˇen´ı mezi dvˇema totoˇzn´ ymi sign´aly pˇrech´az´ı v u ´lohu ne aˇz tak jednoduchou v pˇr´ıpadˇe mˇeˇren´ı zpoˇzdˇen´ı mezi sign´aly, kter´e jsou tvarovˇe odliˇsn´e. Sign´aly sn´ıman´e z p´arov´ ych sond mohou b´ yt ˇcasovˇe navz´ ajem zpoˇzdˇen´e a jsou si tvarovˇe pouze podobn´e, protoˇze kaˇzd´ y ze sign´al˚ u proch´ az´ı jinou ˇc´ ast´ı mozku a proto v kaˇzd´em kan´alu doch´az´ı k jin´e deformaci tvaru. Prost´ ym pˇrekr´ yv´an´ım jednotliv´ ych sign´ al˚ u pˇres sebe, aˇz je nalezena nejvˇetˇs´ı podobnost, vznikla metoda korelaˇcn´ı. Tato metoda m˚ uˇze obˇcas selhat (jeden z pˇr´ıpad˚ u je pˇr´ıtomnost nezanedbateln´eho ˇsumu), proto je tˇreba hledat dalˇs´ı metody. Tyto metody jsou pops´any v n´asleduj´ıc´ıch odstavc´ıch.
2.1
Korelaˇ cn´ı funkce
Jak jiˇz bylo naps´ ano, za touto metodou stoj´ı jednoduch´a myˇslenka vz´ajemn´eho porovn´av´an´ı sign´al˚ u EEG namˇeˇren´ ych r˚ uzn´ ymi sondami. Vz´ ajemn´ a korelaˇcn´ı funkce, pro kterou plat´ı vzorec R12 (τ ) =
∞ X
s1 (k + τ )s∗2 (k), τ ∈ Z,
k=−∞
vystihuje podobnost dvou sign´ al˚ u pro r˚ uzn´e vz´ajemn´e posuvy. Pokud vypoˇc´ıt´ame cel´ y vektor korelaˇcn´ıch koeficient˚ u R(τ ), tzn. pro vˇsechny moˇzn´e posuvy τ , bude jednomu urˇcit´emu τ odpov´ıdat nejvˇetˇs´ı korelaˇcn´ı koeficient R(τ )max . Potom toto τ je hledan´e zpoˇzdˇen´ı. Pro v´ ypoˇcet korelaˇcn´ı funkce existuje v Matlabu funkce xcorr, ale tato funkce vypoˇc´ıt´a cel´ y vektor korelaˇcn´ıch koeficient˚ u. Vzhledem ke skuteˇcnosti, ˇze zpoˇzdˇen´ı mezi sign´aly EEG se pohybuje mezi 0 aˇz 40 ms (tj. zpoˇzdˇen´ı 0 aˇz 10 vzork˚ u pˇri fs = 250 Hz), nen´ı tˇreba vypoˇc´ıt´avat cel´ y vektor, ale jen jeho ˇc´ ast, proto je vhodn´e (vzhledem k uˇsetˇren´emu v´ ypoˇcetn´ımu ˇcasu a pamˇet’ov´ ym n´arok˚ um) funkci pro v´ ypoˇcet
korelace nepatrnˇe upravit. Pˇr´ıklad vektoru korelaˇcn´ıch koeficient˚ u je na obr´azku 2. Na obr´azku 2(b) je zn´azornˇen detail vektoru okolo nulov´eho posuvu. Jen takov´a ˇc´ast vektoru staˇc´ı pro urˇcen´ı zpoˇzdˇen´ı. 120
120
100 100
80 80
60
40
60
20 40
0 20
−20
−40 −10
0
−8
−6
−4
−2
0 delay [s]
2
4
6
8
10
−0.15
−0.1
−0.05
(a) Cel´ y vektor
0 delay [s]
0.05
0.1
0.15
(b) Detail
Obr´ azek 2: Pr˚ ubˇeh korelaˇcn´ı funkce Tato metoda je efektivn´ı a velice jednoduch´a, z toho d˚ uvodu ji uˇz nebudu v´ıce popisovat, jej´ı v´ ysledky pouze pouˇziji pro porovn´ an´ı s dalˇs´ımi metodami a v pˇr´ıpadech sign´al˚ u bez ˇsumu, nebo jen s mal´ ym ˇsumem ji lze povaˇzovat za referenˇcn´ı (z d˚ uvodu nejvˇetˇs´ı pˇresnosti).
2.2
Metoda DFT spektra
2.2.1
Realizace metody
Tato metoda vyuˇz´ıv´ a vlastnost´ı diskr´etn´ı Fourierovy transformace a to, ˇze mimo jin´e nese informaci o f´azi. Diskr´etn´ı Fourierova transformace je definov´ana jako S(Ω) =
∞ X
s(k)e−jΩk .
k=−∞
Pokud chceme pomoc´ı DFT spektra nal´ezt zpoˇzdˇen´ı, je tˇreba vypoˇc´ıtat vz´ajemnou spektr´aln´ı hustotu [4], kter´ a je definov´ ana jako Sxy (ejΦ ) = Sx∗ (ejΦ )Sy (ejΦ ) = Sx (e−jΦ )Sy (ejΦ ), coˇz je opˇet komplexn´ı funkce re´ aln´e promˇenn´e a protoˇze ji lze zapsat tak´e n´asleduj´ıc´ım zp˚ usobem Sxy (f ) = Cxy (f ) + jQxy (f ) = |Sxy (f )|ejΦxy (f ) , tak je moˇzn´e pomoc´ı arkustangensu pod´ılu obou sloˇzek1 vypoˇc´ıtat funkˇcn´ı z´avislost vz´ajemn´e f´aze na frekvenci. V Matlabu je lze pro v´ ypoˇcet vz´ajemn´e spektr´aln´ı hustoty moˇzn´e pouˇz´ıt funkci csd (cross spectrum density), kter´ a pouˇz´ıv´ a Welchovu metodu v´ ypoˇctu a vyhlazen´ı spektra, coˇz je v´ yhodnˇejˇs´ı. Welchova metoda vych´ az´ı z principi´ aln´ı definice spektr´aln´ı v´ ykonov´e hustoty, ovˇsem spektra jsou pr˚ umˇerov´ana. Sign´ al je segmentov´ an na segmenty d´elky 2N , kde N je cel´e ˇc´ıslo a jednotliv´e segmenty se pˇrekr´ yvaj´ı, nejˇcastˇeji je voleno pˇrekryt´ı 50%. Kaˇzd´ y segment je v´ahov´an v´ahovac´ım oknem z d˚ uvodu potlaˇcen´ı jevu prosakov´ an´ı spektra. Z jednotliv´ ych segment˚ u se vypoˇc´ıt´a spektr´aln´ı hustota a hustoty pro jednotliv´e segmenty se zpr˚ umˇeruj´ı. L−1 1 X |Si [k]|2 Sbs [k] = L i=0 M V´ ysledn´ y pr˚ ubˇeh je d´ıky t´eto metodˇe vyhlazenˇejˇs´ı a l´epe se s n´ım pracuje. Samozˇrejmost´ı je zvolen´ı kompromisu mezi poˇctem segment˚ u a d´elkou segmentu, protoˇze velk´ y poˇcet oken a velk´e d´elky segmentu 1C
xy
je souf´ azov´ a a Qxy je kvadraturn´ı sloˇ zka vz´ ajemn´ eho spektra
znamen´a pˇr´ıliˇs velk´e ˇcasov´e zpr˚ umˇerov´ an´ı hodnot a kr´atk´e segmenty znamenaj´ı zhorˇsen´ı frekvenˇcn´ıho rozliˇsen´ı. Blokov´e sch´ema Welchovy metody je zn´azornˇeno na obr´azku 3.
x[n] xl [m]
w[m]
1 | . . . |2 N
DFT
pr˚ umˇer
Sˆx [k]
Obr´ azek 3: Welchova metoda v´ ypoˇctu spektra V pˇr´ıpadˇe, ˇze jiˇz zn´ ame pr˚ ubˇeh vz´ ajemn´e f´aze a nalezneme na kmitoˇctech, kter´e jsou urˇcuj´ıc´ı pro EEG sign´al (cca 0 − 40 Hz) t´emˇeˇr line´ arn´ı pr˚ ubˇeh f´aze, m˚ uˇzeme proloˇzit pr˚ ubˇeh pˇr´ımkou, pˇriˇcemˇz sklon t´eto pˇr´ımky ud´av´ a hledan´e zpoˇzdˇen´ı, kter´e lze vypoˇc´ıtat jako ∆t =
∆ϕ . ∆f · 2π
Pˇr´ıklad pr˚ ubˇeh spektra a f´ aze je zn´ azornˇen na obr´azku 4(a). 2.2.2
Koherence
Dosud byl pops´ an zp˚ usob hled´ an´ı zpoˇzdˇen´ı pomoc´ı proloˇzen´ı pr˚ ubˇehu f´aze pˇr´ımkou tak, ˇze ˇclovˇek musel oznaˇcit m´ısta, odkud kam se m´ a pr˚ ubˇeh prokl´adat. Pokud chceme odeˇcten´ı zpoˇzdˇen´ı zautomatizovat, je tˇreba pouˇz´ıt dalˇs´ı funkci, pro urˇcen´ı pˇribliˇznˇe line´arn´ı ˇc´asti pr˚ ubˇehu f´aze. Toto n´am poskytuje tzv. koherence, coˇz je m´ıra pˇresnosti vz´ajemn´e spektr´aln´ı hustoty a je urˇcena jako kvadr´at vz´ ajemn´e spektr´ aln´ı hustoty dˇelen´ y souˇcinem jednotliv´ ych spektr´aln´ıch hustot 2 γxy (f ) =
|Sxy (f )|2 2 , 0 ≤ γxy (f ) ≤ 1. Sx (f )Sy (f )
Ze vzorce je vidˇet, ˇze koherence m˚ uˇze nab´ yvat hodnot v rozmez´ı od 0 do 1, pˇriˇcemˇz 0 odpov´ıd´a nejmenˇs´ı m´ıˇre pˇresnosti a 1 maxim´ aln´ı m´ıˇre pˇresnosti vz´ajemn´e spektr´aln´ı hustoty. Pokud je koherence na vˇetˇs´ım frekvenˇcn´ım rozsahu dostateˇcnˇe velk´ a, m˚ uˇzeme ˇr´ıci, ˇze na stejn´em frekvenˇcn´ım rozsahu bude i line´arn´ı ˇc´ast pr˚ ubˇehu f´ aze. Pˇr´ıklad pr˚ ubˇehu koherence, jeˇz ukazuje na line´arn´ı ˇc´ast pr˚ ubˇehu f´aze je nakreslen na obr´azku 4(b). 4 20 0
2 φ [rad]
−40
0
S
xy
[dB]
−20
−60
−2
−80 −100
−4 0
20
40
60
80
100
120
140
0
20
40
60
80
100
120
140
80
100
120
140
f [Hz]
f [Hz]
1
4
0.8 koherence
Φ [rad]
2
0
−2
−4
0.6 0.4 0.2
0
20
40
60
80
100
120
140
0
0
20
40
60 f [Hz]
f [Hz]
(a) Spektr´ aln´ı hustota a f´ aze
(b) F´ aze a koherence
Obr´ azek 4: Pr˚ ubˇehy metody DFT spektra
2.3
Parametrick´ a metoda
Parametrick´ ych metod a zp˚ usob˚ u, jak s nimi pracovat, je mnoho. Jednoduˇse ˇreˇceno jde o hled´an´ı racion´aln´ı lomen´e funkce pˇrenosov´e funkce ˇc´ıslicov´eho filtru, kter´a popisuje parametrick´y model. Parametrick´ y model m´a sv˚ uj ˇr´ad M a ve spektru potom modeluje M ˇspiˇcek. Tento ˇr´ad modelu je tˇreba citlivˇe nastavit (odhadnout [1]), protoˇze podhodnocen´ y resp. nadhodnocen´ y model vˇetˇsinou funguje nespr´avnˇe. Z uveden´eho vypl´ yv´ a, ˇze nalezen´ı parametrick´eho modelu ˇr´adu M je aproximace spektr´aln´ı hustoty racion´aln´ı lomenou funkc´ı stupnˇe M metodou nejmenˇs´ıch ˇctverc˚ u. Je zˇrejm´e, ˇze doch´az´ı ke kompresi (redukci) dat, protoˇze sign´ al o d´elce N je nahrazen M < N parametry. Tak´e je vidˇet, ˇze vytv´aˇret parametrick´ y model o ˇr´ adu M → N nem´ a smysl (Pro M = N m´a napˇr. LPC spektrum stejn´ y tvar jako standardnˇe vypoˇc´ıtan´e spektrum). 2.3.1
Typy model˚ u
ARMA je autoregresn´ı model klouzav´ ych souˇct˚ u (autoregressive moving average model). Jde o nejpˇresnˇejˇs´ı model, ale nalezen´ı jeho parametr˚ u je nejsloˇzitˇejˇs´ı. Pˇrenosov´a funkce syntetizuj´ıc´ıho filtru m´a nuly i p´oly. MA je model klouzav´ ych pr˚ umˇer˚ u (moving average model). Pˇrenosov´a funkce syntetizuj´ıc´ıho filtru m´a pouze nulov´e body. Nalezen´ı parametr˚ u vede na soustavu neline´arn´ıch rovnic, stejnˇe jako v pˇr´ıpadˇe modelu ARMA. AR je autoregresn´ı model (autoregressive model). Pˇrenosov´a funkce syntetizuj´ıc´ıho filtru m´a pouze p´ oly a nalezen´ı jeho parametr˚ u vede na soustavu line´arn´ıch rovnic, proto je tento typ modelu pouˇz´ıv´an nejˇcastˇeji. V m´em pˇr´ıpadˇe jsem pouˇzil tak´e AR model. 2.3.2
Estimace sign´ alu
Podstatou metody je estimace sign´ alu (Wienerova filtrace [2]), kterou lze definovat jako filtraci sign´alu line´arn´ım filtrem tak, abychom z´ıskali sign´ al poˇzadovan´ y. Obecn´ y estim´ator je zn´azornˇen na obr´azku 5. Na tomto obr´ azku je zn´ azornˇen filtr s impulzovou odezvou g[n], kter´ y generuje na sv´em v´ ystupu sign´al d[n] u[n]
g[n]
ˆ d[n]
e[n]
Obr´ azek 5: Obecn´ y estim´ator ˆ filtrac´ı sign´ d[n] alu u[n]. Sign´ al d[n] je analyzovan´ y, od kter´eho odeˇc´ıt´ame sign´al filtrovan´ y. V pˇr´ıpadˇe, ˇze ˆ je tento sign´al aproximov´ an nejl´epe (nalezen ide´aln´ı model), bl´ıˇz´ı se sign´al e[n] = d[n] − d[n] b´ıl´emu ˇsumu. Po nalezen´ı parametr˚ u modelu lze pomoc´ı tohoto modelu a nekorelovan´eho ˇsumu modelovat analyzovan´ y sign´al. 2.3.3
ˇ sen´ı Wienerovy filtrace Reˇ
ˇ sen´ı odvodil pan Norbert Wiener v podobˇe Wienerovy-Hopfovy rovnice, kter´a pouˇz´ıv´a jako krit´eria Reˇ minimalizaci stˇredn´ı kvadratick´e chyby. Tento pˇr´ıstup je statistick´ y, proto je tˇreba, aby sign´aly byly stacion´arn´ı. Chyba estimace je d´ ana vztahem ˆ e[n] = d[n] − d[n] a n´asledn´ ymi u ´pravami, kter´e zde nepopisuji, z´ısk´ame hledanou Wienerovu-Hopfovu rovnici Ru G = −Rud → G = −R−1 u Rud , kde matice Ru je symetrick´ a ekvidiagon´ aln´ı matice autokorelaˇcn´ıch koeficient˚ u, vektor G je vektor koeˇ sen´ı t´eto ficient˚ u hledan´e impulzn´ı odezvy a vektor Rud je vektor vz´ajemn´ ych korelaˇcn´ıch koeficient˚ u. Reˇ rovnice lze prov´est r˚ uzn´ ymi metodami (Gaussova eliminace apod.), v Matlabu pro inverzn´ı matici existuje pˇr´ıkaz inv, pot´e je jiˇz ˇreˇsen´ı velice jednoduch´e.
2.3.4
Estimace pouˇ zit´ a pro mˇ eˇ ren´ı zpoˇ zdˇ en´ı mezi sign´ aly EEG
V pˇr´ıpadˇe mˇeˇren´ı zpoˇzdˇen´ı mezi sign´ aly nepotˇrebujeme synt´ezu, ale vystaˇc´ıme si pouze s anal´ yzou. Jako vstupn´ı sign´aly u[n] a d[n] pouˇzijeme dva EEG sign´aly, mezi kter´ ymi chceme mˇeˇrit zpoˇzdˇen´ı. Pˇred zapoˇcet´ım u ´konu hled´ an´ı parametrick´eho modelu je dobr´e nejprve filtrovat vstupn´ı sign´aly doln´ı propust´ı. Sign´ aly EEG obsahuj´ı frekvence do cca 40 Hz a pokud bychom hledali model pro cel´ y rozsah frekvenc´ı, kter´ y n´ am poskytuje z´ aznam z pˇr´ıstroj˚ u, snaˇzil by se algoritmus modelovat sign´al v cel´em spektru. Protoˇze ˇr´ ad AR modelu pro EEG sign´aly se prakticky pohybuje v rozmez´ı M = 4 ∼ 16, mohlo by se st´at, ˇze by na ˇc´ ast spektra, kde se vyskytuje EEG sign´al, zbylo jen m´alo spektr´aln´ıch ˇspiˇcek. To by vedlo na chybn´e urˇcen´ı f´ aze, potaˇzmo zpoˇzdˇen´ı. Dalˇs´ım krokem je jiˇz nalezen´ı impulzn´ı odezvy analyzuj´ıc´ıho filtru. Pomoc´ı t´eto odezvy lze nal´ezt frekvenˇcn´ı charakteristiku FIR2 filtru, coˇz je komplexn´ı funkce re´aln´e promˇenn´e. D´ıky tomu m´ame opˇet moˇznost vypoˇc´ıtat f´ azi v z´ avislosti na frekvenci. V Matlabu je to zajiˇstˇeno opˇet intern´ı funkc´ı, a to funkc´ı freqz. Tato f´aze je vz´ ajemn´ a a lze ji pouˇz´ıt stejn´ ym zp˚ usobem jako v pˇr´ıpadˇe metody DFT spektra, tzn. ˇze na relevantn´ıch kmitoˇctech hled´ ame pˇribliˇznˇe line´arn´ı ˇc´ast pr˚ ubˇehu a proloˇz´ıme ho pˇr´ımkou, jej´ıˇz sklon ud´av´a velikost zpoˇzdˇen´ı. V tomto pˇr´ıpadˇe je ale u ´loha hled´ an´ı line´arn´ı ˇc´asti jednoduˇsˇs´ı, nebot’ v pˇr´ıpadˇe parametrizace jde vlastnˇe o odhad spektra, tud´ıˇz je pr˚ ubˇeh spektra i f´aze vyhlazen. Jak takov´ y pr˚ ubˇeh vypad´a je zn´azornˇeno na obr´azku 6(a). Pˇr´ımku lze proloˇzit velmi jednoduˇse a to v poˇc´ateˇcn´ı oblasti frekvenc´ı aˇz do takov´ ych frekvenc´ı, dokud je pr˚ ubˇeh f´aze dostateˇcnˇe line´ arn´ı. To jsem zaˇr´ıdil postupn´ ym prokl´ad´an´ım a zkracov´an´ım prokl´adac´ı pˇr´ımky (pouˇzil jsem funkci polyfit) tak, aby odchylka regrese byla pod urˇcitou empiricky stanovenou mez. Pˇr´ıklad nalezen´e ˇc´ asti je vyobrazen na obr´azku 6(b). 4
2
H [−]
1.5
3
1
2
0.5
0
20
40
60
80
100
120
140
f [Hz]
φ [rad]
1
0
0
4 −1
Φ [rad]
2 −2
0 −3
−2
−4
−4
0
20
40
60
80
100
120
140
0
20
40
60
(a) Pˇrenosov´ a funkce a f´ aze
80
100
120
140
f [Hz]
f [Hz]
(b) Proloˇ zen´ a f´ aze
Obr´ azek 6: Pr˚ ubˇehy parametrick´e metody
2.3.5
Pr˚ ubˇ eˇ zn´ e metody Wienerovy filtrace
Pr˚ ubˇeˇzn´e metody jsou pouˇzity totoˇznˇe, jako byl dosavadn´ı popis, pouze k z´ısk´an´ı frekvenˇcn´ı charakteristiky resp. f´azov´e charakteristiky filtru nedojde blokovˇe, ale adaptivnˇe. Pr˚ ubˇeˇznˇe jsem Wienerovu filtraci navrhl jednak gradientn´ım stochastick´ ym algoritmem LMS (Least Mean Squares) a RLS algoritmem (Recursive Least Squares) s pˇredv´ ahov´ an´ım. U blokov´eho ˇreˇsen´ı doch´ az´ı k naˇcten´ı kompletn´ıch sign´al˚ u a jejich zpracov´an´ı. U pr˚ ubˇeˇzn´ ych metod je naˇcteno pouze nˇekolik vzor˚ u sign´ al˚ u a po jejich zpracov´an´ı je naˇcten dalˇs´ı vzorek a podle nˇej upraveny parametry modelu. Takto doch´ az´ı k u ´prav´ am parametr˚ u, dokud jsou k dispozici data. Pr˚ ubˇeˇzn´e metody vykazuj´ı tak´e dobr´e v´ ysledky, ale v mnoha pˇr´ıpadech urˇcuj´ı zpoˇzdˇen´ı ˇspatnˇe. To je zp˚ usobeno pravdˇepodobnˇe komplexem hrot–vlna, kter´ y parametry modelu rozhod´ı. Pr˚ ubˇeˇzn´e metody proto nepouˇz´ıv´ am a aby bylo moˇzn´e je pouˇz´ıt, bylo by tˇreba upravit uˇc´ıc´ı algoritmus tak, aby fungoval l´epe, neˇz se mi zat´ım podaˇrilo realizovat. 2 finite
impulse response
3
Testov´ an´ı a modely sign´ al˚ u
Pomoc´ı programu Matlab jsem mˇel moˇznost vˇsechny zde popsan´e metody naprogramovat a pouˇz´ıt re´aln´a a modelovan´a data. Je nutno podotknout, ˇze jsem nemohl hned testovat algoritmy na re´aln´ ych sign´alech, protoˇze u takov´ ych sign´ al˚ u nezn´ am skuteˇcnou hodnotu zpoˇzdˇen´ı. Proto jsem musel vymodelovat takov´ y sign´al, kter´ y by se spektr´ alnˇe i ˇcasovˇe co nejv´ıce podobal sign´alu re´aln´emu. Hlavn´ı vˇec, na kterou jsem se zamˇeˇril, bylo vystihnout tvar komplexu hrot–vlna. To se mi celkem povedlo a potom i spektrum bylo rozprostˇreno hlavnˇe na frekvenc´ıch stejn´ ych jako v pˇr´ıpadˇe re´aln´ ych sign´al˚ u. Teprve potom, kdy jsem otestoval algoritmy na modelovan´ ych sign´alech a vˇsechny tˇri vykazovaly pˇribliˇznˇe stejn´e (spr´ avn´e) hodnoty, viz. tabulka 1 a 2, jsem pˇreˇsel na testov´an´ı sign´al˚ u re´aln´ ych. Po mnoha pokusech jsem doˇsel k z´ avˇeru, ˇze vˇsechny tˇri metody funguj´ı na re´aln´ ych sign´alech spr´avnˇe, totiˇz, ˇze vˇsechny metody vykazovaly pˇribliˇznˇe stejn´ ych hodnot zpoˇzdˇen´ı.
4
Porovn´ an´ı metod
V tomto odstavci pop´ıˇsi, jak´ ym zp˚ usobem jsem porovn´aval jednotliv´e algoritmy. Jako testovac´ı sign´aly jsem pouˇzil moje modely, jejichˇz vz´ ajemn´e zpoˇzdˇen´ı je 20ms. Nejprve jsem testoval, jak se algoritmy chovaj´ı pro r˚ uzn´e d´elky EEG sign´al˚ u. Nejpˇresnˇejˇs´ı v tomto smˇeru je metoda korelaˇcn´ı, kter´ a vesmˇes pro vˇsechny moˇzn´e d´elky vyk´azala spr´avn´e hodnoty. Metoda DFT spektra je nepouˇziteln´ a pro sign´ aly kratˇs´ı neˇz 256 vzork˚ u, coˇz je zp˚ usobeno podstatou Fourierovy transformace. Parametrick´ a metoda pˇri extr´emnˇe kr´atk´ ych sign´alech jiˇz tak´e nefunguje, ale funguje s dostateˇcnou pˇresnost´ı pro podstatnˇe kratˇs´ı sign´aly, neˇz metoda DFT spektra. D´elka sign´ al˚ u [vzorky] 2420 1210 605 302 151 75 37 18 13
Korelaˇcn´ı funkce ∆t [ms] ε [%] 20 0 20 0 20 0 16 −20 16 −20 20 0 20 0 20 0 16 −20
DFT spektrum ∆t [ms] ε [%] 18,71 −6,44 18,67 −6,62 18,92 −5,39 2,02 −89,90 — — — — — — — — — —
Parametrick´a metoda ∆t [ms] ε [%] 18,08 −9,59 18,44 −7,82 17,72 −11,40 17,05 −14,74 15,21 −23,96 21,47 7,37 21,38 6,92 3,19 −84,07 3,20 −84,01
Tabulka 1: V´ ysledky mˇeˇren´ı pro r˚ uzn´e d´elky sign´al˚ u V druh´em pˇr´ıpadˇe jsem pˇriˇc´ıtal k sign´al˚ um b´ıl´ y ˇsum a to se SNR=0dB (ˇsum se stejn´ ym v´ ykonem jako uˇziteˇcn´ y sign´ al) aˇz SNR=10dB. Jak je z tabulky vidˇet, metody DFT spektra a parametrick´a funguj´ı t´emˇeˇr bezchybnˇe i pˇri nejvˇetˇs´ım ˇsumu. Naopak metoda korelaˇcn´ı pro SNR < 5dB jiˇz zaˇc´ın´a chybovat. Nicm´enˇe uˇz to je tak velk´ y ˇsum, kter´ y by se snad v laboratorn´ıch podm´ınk´ach nemˇel vyskytovat. SNR [dB] 0 1 2 3 4 5 6 7 8 9 10
Korelaˇcn´ı funkce ∆t [ms] ε [%] 79 295 40 100 36 80 25 25 21 5 21 5 20 0 20 0 20 0 20 0 20 0
DFT spektrum ∆t [ms] ε [%] 22,0 10,0 18,8 −6,0 21,0 5,0 20,6 3,0 19,5 −2,5 17,5 −12,5 16,0 −20,0 16,7 −16,5 17,0 −15,0 18,7 −6,5 20,0 0,0
Parametrick´a metoda ∆t [ms] ε [%] 20,56 2,81 20,89 4,43 21,00 5,01 20,58 2,90 20,55 2,73 20,68 3,42 20,84 4,19 19,70 −1,51 20,00 0,01 18,81 −5,94 18,34 −8,31
Tabulka 2: V´ ysledky mˇeˇren´ı pro r˚ uzn´e SNR Na obr´azku 7 jsou vyobrazeny pr˚ ubˇehy pro porovn´an´ı chybovosti jednotliv´ ych metod.
80
25
Korelace DFT spektrum AR modelovani
20
70
15
60 10
∆t [ms]
∆t [ms]
50 5
Korelace DFT spektrum AR modelování
0
40
−5
30 −10
20 −15
−20
0
500
1000
1500
2000
2500
10
0
1
2
3
4
délka [vzorky]
(a) V z´ avislosti na d´ elce sign´ al˚ u
5 SNR [dB]
6
7
8
9
10
(b) V z´ avislosti na SNR
Obr´ azek 7: Porovn´ an´ı namˇeˇren´ ych zpoˇzdˇen´ı r˚ uzn´ ymi metodami
5
Z´ avˇ er
Pˇredloˇzen´e algoritmy funguj´ı vˇsechny s dostateˇcnou pˇresnost´ı, ale kaˇzd´ y z algoritm˚ u dosahuje r˚ uzn´ ych v´ ysledk˚ u dle konkr´etn´ıch podm´ınek (d´elka sign´alu a zaˇsumˇen´ı). Pokud uvaˇzujeme sign´aly bez aditivn´ıho ˇsumu, je nejpˇresnˇejˇs´ı metoda korelaˇcn´ı, kter´a ale zaˇc´ın´a chybovat pˇri zvˇetˇsuj´ıc´ım se ˇsumu. Pokud bychom uvaˇzovali promˇennou d´elku sign´ al˚ u i ˇsum, tak pr˚ umˇernˇe vych´az´ı nejspolehlivˇejˇs´ı metoda parametrick´ a. Uvaˇzujeme-li rychlost algoritm˚ u, nejrychlejˇs´ı je metoda korelaˇcn´ı - potˇrebuje N (2N − 1) souˇct˚ u a N souˇcin˚ u, kde N je d´elka sign´ al˚ u. Pokud pouˇzijeme upravenou korelaˇcn´ı metodu a vypoˇc´ıt´ame pouze ξ korelaˇcn´ıch koeficient˚ u, potˇrebuje metoda pouze N (2ξ − 1) souˇct˚ u a N souˇcin˚ u. V´ ypoˇcet f´ azov´eho DFT spektra potˇrebuje N 2 komplexn´ıch souˇcin˚ u a N 2 komplexn´ıch souˇct˚ u, proto je vyuˇz´ıv´ana FFT (Fast Fourier Transformation), kter´a je prakticky dvakr´at rychlejˇs´ı. U AR modelov´ an´ı je tˇreba sestavit matici o rozmˇeru M × M , kde M je ˇr´ad modelu, vyrobit jej´ı inverzi a vyn´asobit s vektorem o velikosti M . Z v´ ysledku je vypoˇc´ıt´ana frekvenˇcn´ı a f´azov´a charakteristika. Tento postup je tˇreba opakovat tak dlouho, dokud nen´ı nalezen optim´aln´ı ˇr´ad, coˇz algoritmus ponˇekud zpomaluje, ale i tak je rychlejˇs´ı neˇz metoda DFT spektra.
5.1
Dalˇ s´ı pr´ ace
Touto prac´ı se jiˇz d´ ale zab´ yvat nebudu, ale v´ yzkum v oblasti epilepsie neopouˇst´ım. Nyn´ı se zab´ yv´ am zkoum´an´ım dostupn´ ych algoritm˚ u pro predikci epileptick´ ych z´achvat˚ u a do budoucnosti se budu snaˇzit nˇekter´ y z nich realizovat v Matlabu.
Reference ˇ [1] Cmejla R. Krit´eria pro urˇcen´ı ˇr´ adu ar modelu pˇri zpracov´an´ı ˇreˇcov´ ych sign´al˚ u. Akustick´e listy, 22:4–7, 2000. ˇ [2] Sovka P., Poll´ ak P. Vybran´e metody ˇc´ıslicov´eho zpracov´ an´ı sign´ al˚ u. CVUT, 1. vyd´an´ı, 2001. ˇ [3] Svatoˇs J. Biologick´e sign´ aly I - Geneze, zpracov´ an´ı a anal´yza. CVUT, 2. vyd´an´ı, 1998. ˇ ıslicov´e zpracov´ ˇ [4] Uhl´ıˇr J., Sovka P. C´ an´ı sign´ al˚ u. CVUT, 2. vyd´an´ı, 2002. Kontaktn´ı informace: Ondˇrej Drbal ˇ ˇ a republika CVUT FEL, K13131, Technick´ a 2, 166 27 Praha 6, Cesk´ tel: (+420) 22435 2820 e-mail:
[email protected]
