ANALÝZA A ZPRACOV`N˝ ØE¨OVÝCH A BIOLOGICKÝCH SIGN`LÙhynek/citing_papers/... · 2011. 1....

ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE

Fakulta elektrotechnickáKatedra teorie obvodů

ANALÝZA A ZPRACOVÁNÍ

ŘEČOVÝCH A BIOLOGICKÝCH SIGNÁLŮ

SBORNÍK PRACÍ 2010

Editoři sborníkuDoc. Ing. Petr Pollák, CSc.Doc. Ing. Roman Čmejla, CSc.

Prosinec 2010

ANALÝZA A ZPRACOVÁNÍŘEČOVÝCH A BIOLOGICKÝCH SIGNÁLŮSBORNÍK PRACÍ 2010

Editoři:Doc. Ing. Petr Pollák, CSc.Doc. Ing. Roman Čmejla, CSc.

[email protected]@fel.cvut.cz

Katedra teorie obvodůhttp://amber.feld.cvut.czvedoucí: Prof. Ing. Pavel Sovka, CSc.http://noel.feld.cvut.cz/speechlab - Laboratoř zpracování řečihttp://amber.feld.cvut.cz/bio - LaBiS - Laboratoř biologických signálů

Foniatrická klinika 1.LF UK a VFNhttp://fonja.lf1.cuni.czvedoucí: Doc. MUDr. Olga Dlouhá, CSc.

Poděkování:Tato publikace vznikla za podpory grantu GAČR 102/08/0707 „Rozpoznávání mluvenéřeči v reálných podmínkáchÿ, GAČR 102/08/H008 „Analýza a modelováníbiomedicínských a řečových signálůÿ a výzkumných záměrů MSM 210000012„Transdisciplinární výzkum v oblasti biomedicínského inženýrstvíÿ a MSM 212300014„Výzkum v oblasti informačních technologií a komunikacíÿ.

Vydalo nakladatelství ČVUT, Zikova 4, 166 36 Praha 6, v roce 2010.

ISBN: 978-80-01-04680-7

Ediční poznámka

Předložený sborník je souhrnem prací realizovaných doktorandy katedry teorie obvodův oblasti číslicového zpracování signálů a aplikačním zaměřením na zpracování biomedi-cínských a řečových signálů a navazuje na sborníky vydávané od roku 2005.

Sborník dává přehled o jednotlivých výzkumných aktivitách řešených ve skupině zpraco-vání signálů na katedře teorie obvodů. Prezentované příspěvky jsou shrnující a podrobnějšíinformace o řešených problémech lze nalézt v odkazovaných pramenech.

V Praze 26. listopadu 2010

Doc. Ing. Petr Pollák, CSc.Doc. Ing. Roman Čmejla, CSc.editoři sborníku

Obsah

Jan Bartošek: Porovnávání algoritmů pro detekci základní frekvence se zaměřenímna řečové signály 1

Tomáš Bořil: Vícerozměrné autoregresní modelování kauzálních vztahů v EEG 9

Vladyslava Demchenko: Analýza srdeční variability v závislosti na dýchání 14

Jaromír Doležal: Systém pro zpracování EEG v reálném čase 22

Radek Janča: Možnosti detekce SOZ v intrakraniálním EEG signálu 29

Jan Janda: Odhad logopedického věku z řeči dítěte 34

Ján Janík: Úvod do selektívnych Zolotarevových kosínov 40

Robert Krejčí: Optimalizace algoritmů rozpoznávačů řeči s využitím architekturvícejádrových signálových procesorů 44

Ondřej Kučera: Zpětnovazební regulace mikromechanických experimentů 52

Tomáš Lustyk: Hodnocení koktavosti a experimenty s adaptivním prahem u baye-sovského spektrálního detektoru 57

Martina Nejepsová: Analýza subjektivního hodnocení dětského věku dle promluv 63

Václav Procházka: Příprava a analýza Českého Web 1T 5-gram korpusu pro použitív jazykovém modelu 67

Barbora Richtrová: Princip aplikace speciální arteterapeutické techniky 74

Jan Rusz: Hodnocení dysfonie u neléčené Parkinsonovy nemoci 78

Jan Sova: Detekce náhlých změn v intrakraniálním EEG pomocí vlastních čísel ko-relační matice 85

Adam Stráník: Klasifikace mezi /s/ a /š/ na základě parametrizace vstupního sig-nálu pomocí LSF 92

Daniel Špulák: Vliv průměrování na možnosti odhadu krevního tlaku z EKG a PPG 99

Václav Turoň: Zolotarevova transformace a spektrální analýza 104

Porovnavanı algoritmu pro detekci zakladnı

frekvence se zamerenım na recove signaly

Jan Bartosek

Ceske vysoke ucenı technicke v Praze, Fakulta [email protected]

Abstrakt: Prıspevek volne navazuje na [1] a zabyva se predevsım objek-tivnım porovnanım implementovanych algoritmu pro detekci zakladnı frek-vence signalu (F0). Pro tento ucel byl navrzen a realizovan hodnotıcı fra-mework vyuzıvajıcı pro porovnanı referencnı databazi a byla ustavena sadahodnotıcıch kriteriı. Z vysledku vyplyva hlavnı nedostatek vetsiny algoritmuv nızke uspesnosti detekce znelych a neznelych useku promluvy. Zaroven jediskutovana otazka optimalnıho casoveho rozlisenı algoritmu.

1. Uvod

Intonace (zmena zakladnı frekvence v case) lidskeho hlasu je jednım z hlavnıch prozo-dickych rysu nası reci. Extrakce prubehu intonace muze mıt ve zpracovanı reci nemalyvyznam [1]. Uloha nenı bohuzel zcela snadna, protoze na rozdıl od zpevu ci prodlouzenefonace ne vsechny useky promluvy majı vlastnost znelosti (tj. jsou tvoreny s vyuzitımhlasivkovych pulzu). Tudız je treba, aby algoritmus detekujıcı zakladnı frekvenci (PitchDetection Algorithm, PDA) delal toto s maximalnı moznou presnostı, ale zaroven spravnerozlisil znele a neznele useky. Objektivnıho vzajemneho porovnanı PDA mezi sebou lzedosahnout pouzitım F0 referencnı databaze s vhodnou sadou kriteriı. Prace predstavujenavrh a realizaci takoveho hodnotıcıho prostredı. Dale popisuje vıce do hloubky nekterenove PDA metody (zameruje se zejmena na MNBFC) a na zaklade vysledku se zabyvatake realizacı jednoducheho detektoru znelych a neznelych useku reci. Krome shrnutıvysledku je poslednı cast venovana zhodnocenı pozadavku na planovany interpunkcnıhodetektoru.

2. Testovane PDA algoritmy

2.1. Standardnı implementovane algoritmy

Byla implementovana vetsina PDA metod teoreticky popsanych v [1], jmenovite kromejiz drıve implementovane autokorelacnı funkce ve frekvencnı oblasti (ACF freq) take au-tokorelace v casove oblasti (ACF time), Average Magnitude Difference Function (AMDF)a kepstralnı metoda (Ceps). Pro pripomenutı jsou rovnicemi (1), (2), (3) a (4) vyjadrenyvztahy pro jejich vypocet.

Jan Bartošek 1

ACFtime(τ) =1

N

N−n−1∑

n=0

x(n)x(n+ τ) (1)

AMDF (τ) =1

N

N−1∑

n=0

|x(n)− x(n+ τ)| (2)

ACFfreq(n) = IFFT{[abs(FFT (x(k)))]2} (3)

Ceps(n) = IFFT{log(abs(FFT (x(k))))} (4)

2.2. Real-time time domain pitch tracking using wavelets

Metoda je detailne popsana v [5] a byla teoreticky opet predstavena jiz v [1]. Jejı uvadenevysledky vypadaly velmi slibne. Behem implementace vsak bylo zjisteno, ze avizovananekolikaurovnova vlnkova transformace je v tomto prıpade ve vysledku v kazde urovnizplostena na pouhou dolnı propust signalu s naslednym podvzorkovanım (5) a testemkandidata na F0 (peak-picking a hledanı centralnı modu casovych diferencı). Pokud nenınalezen, takto transformovany signal postupuje do dalsı urovne. Soucastı algoritmu bymel byt detektor znelych/neznelych useku na zaklade pomeru energiı pocıtanych z tretinaktualnıho okna. Takovy detektor se vsak svymi vysledky ani nepriblizoval realite a protobyl tento PDA vedome testovan bez samostatne faze detekce Voiced/Unvoiced (avsakalgoritmus samotny ma take schopnost ohodnotit usek v krajnıch prıpadech jako neznely).

a(n) = [x(2 ∗ n− 1) + x(2 ∗ n)]/2 (5)

2.3. Merged Normalized Forward-Backward Correlation (MNFBC)

Tato metoda zpracovanı signalu pracujıcı v casove oblasti byla popsana v [3] jako zakladvelmi komplexnıho PDA. Jejım jadrem je pocıtanı dvou proti sobe jdoucıch (auto)korelacı.Rovnice (7)/(8) ukazuje vypocet dopredne/zpetne normalizovane korelace, kde konstantaMAX PER odpovıda casove periode nejnizsı detekovane frekvence, vzdy se pocıtajı zokna delky 4*MAX PER. Prubehy obou funkcı pro referencnı znelou cast promluvy jsouvykresleny na obrazku 1a. Obe funkce jsou nasledne jednocestne usmerneny (NFC’ aNFC’) a pouzity pro vypocet normalizovane

”spojene“ korelace MNFBC (9), jejı prubeh

je na obrazku 1b. Rovnice (6) ukazuje formalnı predpis pro pouzity korelacnı clen.

< xwk[n], xwl

[n] > =2∗MAX PER−1∑

n=0

xw[n + k]xw[n + l] (6)

NFC[t] =< xw0 [n], xwt [n] >√

< xw0 [n], xw0 [n] >< xwt [n], xwt [n] >(7)

NBC[t] =< xw2MAX PER [n], xw2MAX PER−t [n] >√

< xw2MAX PER [n], xw2MAX PER [n] >< xw2MAX PER−t [n], xw2MAX PER−t [n] >(8)

MNFBC[t] =< xw0 [n]], xw0 [n]] > (NFC ′[t])2+ < xw2MAX P ER

[n], xw2MAX P ER[n] > (NBC ′[t])2

< xw0 [n], xw0 [n] > + < xw2MAX P ER[n], xw2MAX P ER

[n] >(9)

2 Jan Bartošek

0 50 100 150 200 250−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1Funkce NFC(t) a NBC(t)

NFC(t)NBC(t)

Casovy interval [vzorku]

Hodnot

anor

mal

izov

ane

kore

lace

(a) Funkce NFC a NBC

0 50 100 150 200 2500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1MNFBC(t)

MNFBC(t)

Casovy interval [vzorku]

Hodnot

anor

mal

izov

ane

kore

lace

(b) Funkce MNFBC

Obrazek 1: Prubehy korelacnıch funkcı na znelem referencnım useku promluvy

2.4. Direct Frequency Estimation (DFE)

Metoda pracuje ciste v casove oblasti. Je detailne popsana v [2] a algoritmus byl prevzatv binarnı podobe. Obsahuje fazi pro rozhodnutı o znelosti daneho useku. Je doposudpouzıvana v mnoha vedeckych pracıch zabyvajıcıch se analyzou hlasu na nasem pracovisti.

3. Otazka optimalnıho casoveho rozlisenı PDA

Pro zodpovezenı teto otazky je treba uvedomit si nektera biologicka fakta o hlasovemustrojı. Casove rozlisenı algoritmu udava, jak casto je pocıtan novy odhad F0 a je danozejmena delkou prekryvu analyzovanych oken signalu. Na jedne strane je nasım cılemzıskat co nejdetailnejsı informaci o prubehu F0, na strane druhe jsou zde prave biofy-zikalnı limitace traktu, dıky kterym vyssı casove rozlisenı od urcite hranice postrada smysl,protoze neprinası zadnou novou informaci a vede jen k narustu vypocetnı narocnosti, kterahraje klıcovou roli zejmena pri nasazenı aplikacı pracujıcıch v realnem case. Dale v textuzmınena F0 referencnı databaze pracuje s casovym rozlisenım 1ms, avsak je k rozhodnutı,zda-li jiz nenı takto vysoke rozlisenı nad limitem hlasoveho ustrojı.Odpoved’ lze nalezt naprıklad ve studii [6] zabyvajıcı se rychlostı zmeny vysky hlasiv-koveho tonu a podle ktere nejvetsı detekovana mluvena intonacnı rychlost v danstine byla50 pultonu za sekundu (50 centu pultonu za 10ms). Lze predpokladat, ze tato rychlostnenı vyrazneji zavisla na jazyku a pro cestinu bude dosahovat podobne hodnoty. Dodejmejen, ze se jedna o limitnı hodnoty, kterych se zrıdkakdy v realnych promluvach dosahuje.Vysledky studie potvrzujı i nasledujıcı ilustrace se dvema ukazkovymi prubehy F0 naceloznelem useku promluvy detekovanymi testovanym algoritmem ACF freq. Obrazek 2aznazornuje jeden z moznych velmi rychlych intonacnıch prubehu (melodemu) tazacı vetys casovym rozlisenım 23ms (FS=11025Hz, posun oken 256 vzorku). V useku nejvetsıhovzestupu je videt, ze takove casove rozlisenı ne zcela dostacuje. Proto bylo v dalsı pracipouzıvano rozlisenı 16ms, jehoz aplikaci na rozechveny zpevovy signal (vibrato) lze videtna obrazku 2b. Tato hodnota se zda byt dostacujıcı a byla proto zvolena jako vychozı ipro testovanı vsech F0 algoritmu.

Jan Bartošek 3

(a) Rozlisenı F0 23ms, velmi rychla otazka (b) Rozlisenı F0 16ms, vibrato zpev

Obrazek 2: Vliv casoveho rozlisenı F0 na detekovany prubeh intonace

4. F0 referencnı databaze

K tomu, abychom dokazali objektivne ohodnotit PDA, potrebujeme znat spravne vysledkypro vstupnı data. Proto byla v praci pouzita rucne znackovana cast Speecon SpanishDatabase, jejız referencnı cast byla vytvorena jako vysledek pitch-marking1 algoritmu [4]a dale rucne prekontrolovana. Referencnı F0 jsou potom zıskany jako prevracene hodnotycasovych vzdalenostı sousednıch pitch-marku.Databaze byla nahrana s nasledujıcımi parametry: RAW audio format se vzorkovacıfrekvencı 16kHz, bitova hloubka 2B/vzorek, linearnı kodovanı, mono. V nahravkach jeprıtomno 60 mluvcıch (30 zen a 30 muzu). Celkova doba zaznamu presahuje jednu ho-dinu, coz predstavuje v prumeru vıce nez 1 minutu na mluvcıho. Databaze byla nahravanasoucasne ctyrmi mikrofony umıstenymi v ruznych vzdalenostech od mluvcıho a lisıcıchse tak odstupem uzitecneho signalu od sumu (SNR). Promluvy jsou take odlisitelneprostredım, ve kterych byly nahrany (automobil, kancelar a verejna mısta). Krome re-ferencnıch F0 hodnot s intervalem 1ms obsahuje take okamziky zmınenych pitch-markua informaci o znelosti daneho useku.

5. PDA hodnotıcı framework

5.1. Motivace

Motivacı pro vznik hodnotıcıho prostredı pitch-detection algoritmu je nejenom moznostjejich vzajemne objektivnıho porovnanı oproti zname referenci, ale naprıklad i hledanıoptimalnıch parametru v nastavenı daneho PDA. Rozlicna hodnotıcı kriteria a kategoriepote umoznı vybrat vhodny PDA pro potreby aplikace (napr. chybovost urcenı znelostiversus presnost urcenı F0).

5.2. Typy souboru na vystupu PDA

Algoritmy detekujıcı zakladnı frekvenci mıvajı obvykle svuj vystup v jednom ze trı pouzıvanychformatu:Typ 1 oznacuje nativnı format .pda souboru referencnı databaze. V souborech nenı uve-dena specialnı informace o casovem kroku (predpoklada se, ze je znama a-priori), a proto

1Pitch-mark je presne definovany okamzik hlasivkoveho cyklu

4 Jan Bartošek

soubor zacına prımo detekovanou frekvencı, kazda dalsı je pak na novem radku. V hod-notach je zakodovana informace o znelosti useku (hodnota 0 znacı ticho, hodnota 1 neznelyusek a hodnota vetsı nez 1 validnı F0).Typ 2 je velmi podobny typu 1, avsak na prvnım radku obsahuje informaci o casovemkroku v milisekundach.Typ 3 neodpovıda zadnemu z predchozıch, protoze nema konstantnı casovy krok. Najednom radku obsahuje dve cısla udavajıcı jednak detekovanou zakladnı frekvenci a takedelku jejıho trvanı (napr. v sekundach). Tento typ ulozenı dat komprimuje informaci,protoze eliminuje delsı sekvence stejnych hodnot a je tak pamet’ove nejuspornejsı.

5.3. Navrh a implementace hodnotıcıho frameworku

Obrazek 3 predstavuje zakladnı blokove schema architektury frameworku. Jadrem je F0referencnı databaze, ktera obsahuje jednak testovacı promluvy, tak jejich referencnı .pdasoubory s prubehy F0. Seznam testovacıch souboru je predan spoustecımu bloku, ktery jeodpovedny za zavolanı daneho PDA na kazdem ze zadanych testovacıch souboru zvlast’.Je schopen volat ruzne typy implementacı PDA (binarnı forma operacnıho systemu nebovyvolanı prostredı Matlabu z prıkazoveho radku pro PDA ve forme m-filu). Blok takezarucuje prıpadny prevod vstupnıho souboru do formatu vhodneho pro dany PDA. Vprıpade, ze algoritmus nema fazi detekce znelych/neznelych useku, lze mu predradit blokV/UV, ktery mu dale ke zpracovanı preda jen useky vyhodnocene jako znele. Nasledne jesjednocen typ vystupnıch .pda souboru a ty jsou jednotlive porovnany s referencemi.Vysledkem je jednotlivy hodnotıcı soubory pro kazdou promluvu zvlast’. Nad temitovysledky jsou dale pocıtany globalnı hodnocenı a dale hodnocenı naprıc ruznymi kate-goriemi vstupnıch promluv (kanaly, prostredı, pohlavı) a jejich kombinacemi (napr. jenkanal 1 v prostredı automobilu).Framework byl implementovan pod operacnım systemem UNIXoveho typu jako kombi-nace skriptu okolo referencnı databaze napsanych v multiplatformnım interpretovanemjazyce PERL (vlastnı logika) a BASH (jednoduche souborove operace). Cele prostredı jetedy velmi snadno prenositelne i na ostatnı platformy. Mısto pouzite referencnı databazemuze byt s minimalnım usilım pouzita jina F0 referencnı databaze a cely framework bypak mohl najıt upotrebenı i mimo recove systemy (naprıklad v hudebnım odvetvı).

5.4. Hodnotıcı kriteria

Existuje sada ustalenych kriteriı pouzıvanych na poli hodnocenı PDA, cılem prace vsakbylo take smysluplne navrhnout nektera kriteria nova. Chyby znelosti (Voiced Errors -VE) se vyskytnou, pokud algoritmus znely usek prohlası za neznely. Unvoiced Errors (UE)znacı chybu opacneho typu. Gross Error High (Low) je prıpad chyby, kdy detekovana frek-vence v Hz prekrocı 20% toleranci shora (zdola) a nese oznacenı GEH (GEL). Jelikoz tatotolerance je vzhledem k vyssım pozadavkum dnesnı doby k algoritmum pomerne benevo-lentnı, byla nove zavedena take tolerance 10% (oznacenı GEH10 a GEL10). Nekdy jsoutake pouzıvany soucty UE+VE a GEH+GEL. Dale byly nove zavedeny tzv.

”halving HE

(doubling DE)“ chyby, ke kterym dochazı pri odhadu dvojnasobne (polovicnı) frekvence.Gross Error chyby s obema tolerancnımi pasmy jsou take evidovany v ramci jednotlivychneprekryvajıcıch se frekvencnıch pasem (napr. pro recove signaly pet 2/3 oktavovychpasem v rozmezı cca 60-560Hz). V literature lze take nalezt statisticka kriteria (abso-lutnı rozdıl mezi strednımi hodnotami reference a odhadu a take absolutnı rozdıl jejichsmerodatnych odchylek) pocıtana v Hz pres cely soubor promluv, avsak tyto hodnotynemajı dıky logaritmickemu prubehu vnımanı frekvence prılisnou vypovıdajıcı hodnotu.Proto je pouzito modifikovanych statistickych kriteriı [2] - strednı odchylka ∆% (10) a

Jan Bartošek 5

Obrazek 3: Blokove schema architekrury PDA hodnotıcıho frameworku

smerodatna odchylka δ% (11) pocıtanych v centech pultonu:

∆% =1200

N

N∑

n=1

log2

Fest(n)

Fref (n)(10)

δ% =

√√√√ 1

N

N∑

n=1

[1200 log2

Fest(n)

Fref (n)−∆%]2 (11)

Pro vysvetlenı jeste dodejme, ze kriterium VE+UE nenı pouhym souctem obou chybo-vostı, ale je to soucet vsech chybne klasifikovanych useku ku celkovemu poctu useku.Aby nedochazelo k akumulaci chyb, postupujı do faze hodnocenı presnosti jen takoveznele useky, ktere byly spravne klasifikovany jako znele. To muze v urcitych prıpadechzvyhodnit v kriteriu tolerance F0 algoritmy s vysokou hodnotou VE, protoze takove me-tode se presnost F0 pocıta jen z useku, ktere oklasifikovala jako znele a problematickeuseky, ktere mohou snizovat presnost napr. algoritmum bez V/UV faze, jiz v hodnocenıpresnosti prıtomny nejsou.

5.5. Vysledky

Testovany byly vsechny algoritmy uvedene v kapitole 2.. Vysledky testovane sady algo-ritmu na kanalech 0 (nejvyssı hodnota SNR, close-talk mikrofon) a 1 (klopovy mikrofon)uvadı tabulky 1 a 2. Nektere z algoritmu (ACF time, AMDF, CEPS a MNBFC) byly naim-plementovany bez rozhodovacıch prahu pro znele a neznele useky a nebyl jim predrazenani zadny samostatny detektor. Proto detekujı vsechny neznele useky promluv jako znelea kriterium Unvoiced Error u nich dosahuje hodnoty 100 %. Zajımave je, ze algoritmyAMDF a CEPS dosahli, ackoliv pracujı na zcela odlisnem principu, temer totoznychvysledku ve vsech sledovanych kriteriı. Z pohledu vysledku se tedy zdajı ekvivalentnı,ackoliv lze nalezt informace o tom, ze nektere komplexnejsı PDA jsou zalozeny pravena jejich komplementarite, ktera se v nasem testu vubec neprojevila. Metoda MNBFCpodala bohuzel v presnosti vysledky horsı nez ACF freq. Potvrdilo se, ze rostoucı zastou-penı sumu v signalu je velkym problemem (napr. obrovske zhorsenı GEH pro Acf freq).

6 Jan Bartošek

PDA VE UE VE+UE GEH GEL GEH10 GEL10 DE HE[%] [%] [%] [%] [%] [%] [%] [%] [%]

ACF freq 44,4 23,5 31,6 1,2 0,1 1,5 0,18 0,4 0,06ACF time 0 100 61,9 4,7 2,3 6,2 3,5 0,8 1,3

AMDF 0 100 61,9 0,6 27,2 1,4 28,3 0,1 16,2CEPS 0 100 61,9 0,6 27,1 1,4 28,1 0,1 16,0DFE 26,6 15,5 20,4 8,4 4,2 16,5 8,9 0,2 1,3

Wavelets 67,7 11,3 32,7 2,5 4,9 3,7 6,0 1,1 3,9MNBFC 0 100 61,9 4,8 4,4 6,6 6,6 0,4 2,8

Tabulka 1: Souhrnne vysledky na kanalu 0

PDA VE UE VE+UE GEH GEL GEH10 GEL10 DE HE[%] [%] [%] [%] [%] [%] [%] [%] [%]

ACF freq 52,7 34,1 41,3 23,3 0,1 23,5 0,2 3,2 0,03ACF time 0 100 61,9 28,8 2,5 29,8 3,4 3,6 1,5

AMDF 0 100 61,9 10,8 44 11,3 45,2 1,3 21,4CEPS 0 100 61,9 10,1 43,4 10,5 44,7 1,3 21,4DFE 45,4 11,1 25,9 8,5 8,1 17,9 13,1 0,05 4,3

Wavelets 70,4 9,5 32,6 14,3 9,9 17,4 11,6 4,3 6,7MNBFC 0 100 61,9 29,1 4,9 30,4 6,5 2,1 3,1

Tabulka 2: Souhrnne vysledky na kanalu 1

Nejrobustnejsı chovanı naopak vykazuje DFE.

5.6. Pozadavky pro interpunkcnı detektor

Pro ucely pouzitı s dalsımi cıli vyzkumu (interpunkcnı detektor) budeme hledat algo-ritmus s minimalnı chybou kriteria UV (chybne detekovane F0 na neznelych usecıch bymohly vyrazne narusit prubeh melodemu), naopak urcitou chybovost VE lze tolerovat(jen chybejıcı bod melodemu). Celkove vysledky presnosti (GEH) algoritmu na kanalu 0by mely ucelu postacovat.

5.7. Dodatecny experimentalnı V/UV blok

Na zaklade predchozıch vysledku byl implementovan dle [4] jednoduchy detektor znelych aneznelych useku vychazejıcı z pomeru energie signalu (E) a poctu pruchodu nulou (ZCR).Energie je pocıtana z predzpracovaneho okna, kdy jsou zvyrazneny periodicke struk-tury znelych segmentu pomocı kratkodobe energeticke obalky. Predpis pro jeho vystupnıfunkci, ktera se nakonec prahuje empiricky zjistenou hodnotou, je vyjadren vztahem (12).Lze se totiz domnıvat, ze pro znele useky je hodnota funkce EZR vysoka, protoze energiesignalu je pomerne velka a pocet pruchodu nulou je ve srovnanı se sumem nızky. Naopak,pro neznele useky s vysokym ZCR a nızkou energiı signalu bude vysledna hodnota mala.Hodnotıcı framework byl obohacen o modul umoznujıcı evaluaci i samotnych V/UV algo-ritmu a takto implementovany detektor dosahl na kanalu 0 vysledku o neco horsıch nezmetoda DFE (UE+VE: 24,3 % pro EZR versus 20,4 % DFE) a mohl by tak byt uspesnepredrazen algoritmum bez detekce znelych useku a vylepsit tak jejich celkove vysledky.

Jan Bartošek 7

EZR[m] =E[m]

ZCR[m](12)

6. Zavery

Krome implementace zakladnıch PDA byly zkoumany i pokrocilejsı metody. Experi-mentalne bylo overeno, ze casove rozlisenı 16ms je pro sledovanı prubehu intonace recivhodnou hodnotou. Pro objektivnı hodnocenı F0 detekcnıch algoritmu byl navrzen a im-plementovan hodnotıcı framework obalujıcı F0 referencnı databazi. Byla navrzena sadaruznorodych kriteriı umoznujıcı detailnı analyzu chovanı algoritmu. Podle ocekavanı klesauspesnost vsech algoritmu s klesajıcım SNR. Metoda detekce zalozena na sjednocenychnormalizovanych korelacıch (MNBFC) bohuzel nepodala ocekavane vysledky co do presnostiurcenı F0. Vysledky take ukazujı, ze nejvetsı slabinou je obecne faze detekce znelych aneznelych useku.

Podekovanı

Tento vyzkum byl podporovan z grantu GACR 102/08/H008 “Analyza a modelovanıbiomedicınskych a recovych signalu” a GACR 102/08/0707 “Rozpoznavanı mluvene reciv realnych podmınkach”.

Reference

[1] Bartosek, J. Prozodie, zjistenı a vyuzitı zakladnıho tonu v rozpoznavanı reci. Seminarekatedry teorie obvodu, analyza a zpracovanı recovych a biologickych signalu - sbornıkpracı 2009 (2009), 1–8.

[2] Boril, H.; Pollak, P. Direct time domain fundamentalfrequency estimation of speechin noisy conditions. in Proceedings of EUSIPCO 2004 (European Signal ProcessingConference, Vol. 1) (2004), 1003–1006.

[3] Kotnik, B.; et al. Noise robust f0 determination and epoch-marking algorithms. SignalProcessing 89. (2009), 2555–2569.

[4] Kotnik, B.; Hoge, H.; Kacic, Z. Evaluation of pitch detection algorithms in adverseconditions. Proc. 3rd International Conference on Speech Prosody, Dresden, Germany(2006), 149–152.

[5] Larson, E. Real-time time domain pitch tracking using wavelets. Journal of theAcoustical Society of America (2005), 111(4).

[6] Xu, Y.; Sun, X. Maximum speed of pitch change and how it may relate to speech.Journal of Acoustical Society of America, Vol. 111, No. 3 (2002), 1399–1413.

8 Jan Bartošek

Multivariate Autoregressive Modelling of Causal

Connections in EEG

Tomas Boril

Czech Technical University in Prague, Faculty of Electrical [email protected]

Abstract: Conditional Granger causality and related methods have beenused in neurophysiology in recent years for revealing causal relations in brain.Multivariate autoregressive models can model dynamic relations in electroen-cephalography and derived estimators can measure not only the strength ofconnections but also their directions which are important for understandingbrain function during specific tasks. This paper propose an extension of theconditional Granger causality in order to improve a differentiation betweenstrong and weak causal connections.

1. Introduction

In the late ‘60s of the 20th century, econometrics started examination of causal relationsamong time series in order to create economic models and predict a behavior of variableslike agricultural prices. The basic idea of causality was formulated by Wiener (1956)and formalized later in the scope of autoregressive models by Granger in 1969 [4]. Ifprediction of one time series could be improved by knowledge of past values of anotherone, then we say the second series has a causal influence on the first one. Although thistechnique received a great acceptance, some studies like [7] warn against inconsiderate useof Granger causality which can lead to incorrect answers. For instance, the thoughtlessuse can find statistically significant causal relations which do not have its match inthe real word sense. It is necessary to analyse data with autoregressive nature, theexamined segment must be stationary and sufficiently long and all time series with possibleinfluence must be included into the model. While meeting of such requirements can bedifficult in econometrics, a whole new area opens for Granger causality in the form ofneurophysiology at the end of 20th century. Functional magnetic resonance imaging(fMRI), electroencephalography (EEG) and magnetoencephalography (MEG) bring datawhich meets conditions of causality analysis very well and together with progress ofcomputational technology it is possible to find causal relations in large data of brainactivity records, e.g. [5, 2].EEG processing is the oldest method of human brain function analysis. Among advan-tages, good availability should be mentioned (short waiting periods for investigation),low price and simplicity (no problems with patients suffering claustrophobia like withfMRI). However, the main advantage is a high temporal resolution of measured data(important for causal connectivity analysis) and direct relationship with electric activity

Tomáš Bořil 9

of brain neurons. Surface data recorded with a small number of electrodes gives onlyvery inaccurate spatial location of sources of such an activity. Electric voltage potentialsmeasured on one electrode are results of summation of electric activity of large amount ofneurons. High resolution surface EEG recordings bring a higher spatial resolution thanconventional cerebral electromagnetic measures.

2. Renormalized Granger Causality

Let’s assume k stationary stochastic processes v1 . . .vk. Process v1 can be fitted intomultivariate autoregressive (MVAR) model of order p:

v1 (t) =p∑

j=1

a1,jv1 (t− j) +p∑

j=1

a2,jv2 (t− j) + · · ·+p∑

j=1

ak,jvk (t− j) + ε1 (t) (1)

where t stands for an index in a discrete time instance and the prediction error ε1 is awhite noise. To estimate the MVAR model coefficients one can minimize the variance ofthe prediction error via MVAR estimators discussed in [6]. Model order can be estimatedby Akaike Information Criterion (AIC) or using a method published in [1].Let’s define notation x(t−1) for vector of all values of x excluding the last sample x (t).Then we can define conditional variance [3]:

var(v1|v1(t−1), . . . ,vk(t−1)

)= var (ε1) , (2)

signifying the variance of prediction error when actual sample of v1 is expressed viaprevious values of v1 . . .vk.Conditional Granger causality (CGC) is often defined for three variables like causalityfrom v1 to v2 conditional on v3 [3]:

Fv1→v2|v3 = lnvar

(v2|v2(t−1),v3(t−1)

)

var(v2|v1(t−1),v2(t−1),v3(t−1)

) (3)

which measures the causality from v1 to v2 but by the reason of including also v3 intothe model a direct connection is distinguished from an indirect one when the informationflow from the first variable to the second variable is completely mediated by the thirdvariable.We can easily generalize the equation (3) to a form with k variables:

Fv1→v2|v3,...,vk= ln

var(v2|v2(t−1), . . . ,vk(t−1)

)

var(v2|v1(t−1), . . . ,vk(t−1)

) . (4)

The numerator denotes variance of prediction error of target variable in MVAR model notincluding the source variable, whereas the denominator includes the source variable. Thisdefinition corresponds with the original causal idea definition perfectly – if the inclusion ofthe previous values of the source variable decreases the variance of the prediction error, theconditional Granger causality is a non-zero positive value and we say the source variableGranger-causes the target variable.Let’s define maximal prediction ratio:

MPRvi=

var (vi)

var(vi|v1(t−1), . . . ,vk(t−1)

) (5)

10 Tomáš Bořil

determining the ratio of how many times the variance of the variable is reduced in the formof residual prediction error when it is modelled via MVAR model including all variables.In the case of conditional Granger causality, strong relations can be falsely detected ifpredictions of variables are overall low and small changes in the value can lead to largeresult of the ratio. Let’s define renormalized Granger causality (RGC):

RGCv1→v2 = Fv1→v2|v3,...,vk

k∏

i=1

MPRvi(6)

which takes into account the rate of prediction of each variable included in the MVARmodel.

3. Experiments and Results

We have generated 2000 samples of three MVAR signals of order 5 with causal influencesv1 → v2 → v3 corresponding to [1] in order to compare the results:

v1(t) = 0.9v1(t− 1)− 0.3v1(t− 4) + ε(t),v2(t) = 0.8v2(t− 1)− 0.5v2(t− 2)+

+ 0.16v1(t− 1)− 0.2v1(t− 2)++ 0.2v1(t− 5) + η(t),

v3(t) = −0.2v3(t− 2)− 0.4v3(t− 5)−− 0.27v2(t− 1) + 0.1v2(t− 3) + γ(t)

(7)

where ε, η and γ are Gaussian white noises with zero means and variances var (ε) = 1,var (η) = 0.7 and var (γ) = 0.4.In the next step, Gaussian white noises with zero means were mixed with test signalswith different signal-to-noise ratio levels (SNR) in order to analyse a behavior of CGCand RGC in noise conditions. In figures 1 to 3, we can see the recognition of causalitiesfor different SNR levels and signal lengths, the recognition decreases rapidly with lowerSNR. It is obvious with shorter signal segments the noise immunity of CGC and RGCdecreases. However, RGC behaves much better in low SNR conditions, it suppresses falsecausality detections when they cannot be estimated from the data segment well.

4. Conclusions

A new approach of causality analysis is presented, a normalization extension of well-accepted conditional Granger causality has been proposed. Presented results show therenormalized Granger causality can distinguish better a case when the prediction is overalllow from the case when data have highly predictive character. This property can behelpful in real data processing when the signal-to-noise ration value is unknown andthe conditional Granger causality can estimate high causality because both values innumerator and denominator are small and little statistically insignificant changes canlead to great change in the ratio result. Renormalized Granger causality improves thisproperty and a demonstration on artificial data is shown in this paper.

Tomáš Bořil 11

−200

0 20 40 60 80 100

0.05

0.1

0.15

0.2

0.25

SNR (dB)

F

v1 → v2v2 → v3v1 → v3v2 → v1v3 → v1v3 → v2

−200

0 20 40 60 80 100

0.5

1

1.5

2

2.5

3

3.5

SNR (dB)

RG

C

v1 → v2

v2 → v3

v1 → v3

v2 → v1

v3 → v1

v3 → v2

Figure 1: Additive white noise, signal length 2000 samples, CGC and RGC indexes

−200

0 20 40 60 80 100

0.05

0.1

0.15

0.2

0.25

SNR (dB)

F

v1 → v2v2 → v3v1 → v3v2 → v1v3 → v1v3 → v2

−200

0 20 40 60 80 100

0.5

1

1.5

2

2.5

3

3.5

SNR (dB)

RG

C

v1 → v2

v2 → v3

v1 → v3

v2 → v1

v3 → v1

v3 → v2


−200

0 20 40 60 80 100

0.05

0.1

0.15

0.2

0.25

SNR (dB)

F

v1 → v2v2 → v3v1 → v3v2 → v1v3 → v1v3 → v2

−200

0 20 40 60 80 100

0.5

1

1.5

2

2.5

SNR (dB)

RG

C

v1 → v2

v2 → v3

v1 → v3

v2 → v1

v3 → v1

v3 → v2


12 Tomáš Bořil

Acknowledgments

Research described in this paper has been supported by SGS10/176/OHK3/2T/13 “Brainactivity mapping and analysis”, Czech Grant Agency under grant No. GD102/08/H008“Analysis and modelling of biomedical and speech signals” and by research program“Transdisciplinary Research in Biomedical Engineering” No. MSM6840770012 of CzechTechnical University in Prague.

References

[1] Boril, T. Revealing of Relations in EEG via Granger Causality. 3th InternationalStudent Conference on Electrical Engineering (2009), 1–4.

[2] Brovelli, A.; Ding, M.; Ledberg, A.; Chen, Y.; Nakamura, R.; Bressler, S. L.Beta oscillations in a large-scale sensorimotor cortical network: directional influencesrevealed by Granger causality. Proceedings of the National Academy of Sciences of theUnited States of America 101, 26 (June 2004), 9849–9854.

[3] Chen, Y.; Bressler, S. L.; Ding, M. Frequency decomposition of conditional Grangercausality and application to multivariate neural field potential data. Journal ofNeuroscience Methods 150, 2 (2006), 228–237.

[4] Granger, C. W. J. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 3 (July 1969), 424–38.

[5] Liang, H.; Ding, M.; Nakamura, R. R.; Bressler, S. L. Causal influences in primatecerebral cortex during visual pattern discrimination. Neuroreport 11 (2000), 2875–2880.

[6] Schlogl, A. A comparison of multivariate autoregressive estimators. Signal Process.86, 9 (2006), 2426–2429.

[7] Ziemer, R. F.; Collins, G. S. Granger causality and U.S. crop and livestock prices.Southern Journal of Agricultural Economics 16, 01 (1984).

Tomáš Bořil 13

Analýza variability srdečního rytmuv závislosti na dýchání

Vladyslava Demchenko

České vysoké učení v Praze, Fakulta elektrotechnická[email protected]

Abstrakt: Tématem příspěvku je záznam signálu srdeční variability (HRV) a jejích zpracování s důrazem na analýzu v časové oblasti a zjištění podobnosti průběhu HRV s dýcháním. V rámci této práce se prováděli testy srdeční variability a také závislost změny srdečního rytmu během spontánního dýchání.

1. Úvod

Práce popsána v tomto článku se zabývá signálem, vzniklým z faktu, že srdeční tep není ani během dýchání ani po dobu zadržení dechu konstantní. Tato srdeční variabilita (HRV) je řízená zejména aktivitou autonomní nervové soustavy (ANS) a frekvencí dýchání. V rámci vlivu ANS na srdeční variabilitu můžeme také vyjmenovat takové důležité faktory jako jsou stres, fyzická aktivita a silné emoce. Dýchání má různě silný vliv na HRV podle toho, zdá je řízené (až 50% vliv) nebo spontánní (až 10% vliv). Přičemž pří řízeném dýchání zároveň klesá vliv ANS [1]. Parametry srdeční variability se zkoumají v časové a frekvenční oblasti a používají se v klinické praxi pro kontrolu stavu kardiovaskulárního systému. Mohou posloužit k včasné predikci náhlé smrti, diagnostice srdečních nemocí, komplikujících například diabetes melitus apod [2]. Pro posouzení srdeční variability se používají různé testy: ortostatický, valsalvův, test zadržení dechu a test hlubokého dýchání (poslední dva jmenované byly použity v této práci). Křivka srdeční variability, která byla obdržená po provedení zmíněných testů, byla zpracována v časové a frekvenční oblasti. Protože zpracování ve frekvenční oblasti se zabývá poměrně hodně studii v této práci jsme se zabývali zejména zpracování signálu srdečné variability časové oblasti a také zpracování signálu dýchání.

14 Vladyslava Demchenko

2. Databáze signálů

Pro účely této práce byla naměřená databáze signálů od 14 dobrovolníků. Všichni dobrovolnici byli zdraví a mladí, pouze jeden pravidelně vykonává různá cvičení hlubokého dýchání pro zlepšení srdeční variability. Mezi dobrovolníky byla také jedna žena.

2.1 Druhy naměřených signálů

Protože měření bylo prováděno pomocí přístroje Biopac MP35, který je schopen měřit současně na čtyřech kanálech, měřili jsme 4 signály pro každého dobrovolníka. Tyto signály jsou postupně:

1. elektrokardiogram (EKG) - měřeno na druhém svodu2. fotopletysmogram (PPG) - tento signál nebyl v této části zpracování využit3. signál dýchání měřený z pásu, obepínajícího hrudník4. signál dechu měřený pomocí termistoru, umístěného u nosu

2.2 Rozdělení dýchání na řízené a spontánní

Měření se skládalo z dvou částí. Během první části se měřili signály pří klidném neřízeném dýchání celkovou délkou trvání 5min. V druhé fáze bylo měřeno řízené dýchání pří různých dýchacích frekvencích. Délka trvání jednotlivých fázi byla stanovena na cca 2 min (s tím, že vždy fáze končila výdechem). Frekvence dýchání byly zvolené následovně (notace: [doba nádechu - doba zadržení dechu v nádechu - doba výdechu - doba zadržení dechu ve výdechu], doby jsou udávané v sekundách).

1. 4 - 1 - 4 - 12. 5 - 1 - 5 - 1 (test hlubokého dýchání před apnoe)3. fáze apnoe (zadržení dechu) - doba trvání 62s4. 5 - 1 - 5 - 1 (test hlubokého dýchání po apnoe)5. 4 - 0 - 4 - 0 6. 4 - 1 - 8 - 1

Obrázek 1: Signály naměřené během řízeného dýchání. Frekvence dýchání se měnili po cca 2 minutách a byly řízené programem. Apnoetická fáze trvala 62 sekundy.

Vladyslava Demchenko 15

2.2 Program pro řízení dechu

Pro přesné řízení dýchání měřených jedinců byl navržen program figB.m, který měřil čas jednotlivých fázi dýchání a pomocí nabíhajícího (resp. sebíhajícího) sloupcového grafu reprezentujícího plíce naznačoval jak má dotyčný dýchat. Program lze spustit pomocí příkazu v prostředí MatLab. Důvodem napsaní takového programu byla hlavně nutná možnost zadávání více frekvenci dýchání do jednoho měřicího cyklu. Hlavní obrazovka je uvedená na obrázku 2.

Obrázek 2: Hlavní obrazovka programu pro řízení dechu

3. Metody zpracování

Analýza signálu byla rozdělená do dvou hlavních části a to do analýzy v časové oblasti a analýzy v oblasti frekvenční. Pro analýzu v časové oblasti byly použité jak lineární tak nelineární metody zpracování.

3.1 Analýza v časové oblasti

Z mnoha hodnot, které byly vypočítané v rámci použití lineárních metod zpracování signálu v časové oblasti stoji za zmínění zejména následující.Hodnota RRmean nebo střední hodnota délky intervalu mezi sousedními R-spičkami po dobu měření jednotlivých dechových sekvencí.Hodnota NN50, udávající počet RR intervalů kratších než 50 ms (často se používá pro diagnostiku kvality HRV).Směrodatná odchylka od střední hodnoty všech RR intervalů naměřených během jedné dechové frekvence (SDNN).K diagnostice v klinické praxi se také používá hodnota rMSSD (root mean square of successive diferences), která se počítá podle vztahu:

rMSSD= 1N−1 ∑i=1

N−1

RRi1−RRi2 [ms; -, ms, ms] (1)

kde N je počet všech RR intervalů za daný period měření.


Nelineární metodou, která byla v této práci použita je Poincarého graf, který udává závislost intervalu na následujícím. Na obrázku 3 je vidět příklad použití tohoto grafu na reálných datech z naměřené databáze. Poincarého graf je dobrým nástrojem pro posouzení kvality srdeční variability. Pokud je HRV pravidelné na grafu by se mela vytvořit elipsoidní tvar a body grafu by měli ležet blízko u os.

(a) (b)Obrázek 3: Zobrazení Poincarého grafu pro 5 minut trvající záznam (a) neřízeného dýchání a (b) řízeného dýchání během testu hlubokého dýchání

3.1 Analýza ve frekvenční oblasti

Druhou částí analýzy byla analýza signálů ve frekvenční oblasti. Během řízeného dýchání srdeční variabilita se přizpůsobí dechu a křivka poměrně dobře připomíná sinusový tvar, proto se také nejčastěji zkoumá ve frekvenční oblasti. Zkoumají se zejména parametry výkonu ve tří hlavních frekvenčních pásmech: velmi nízkém (pod 40 mHz), nízkém (mezi 40 a 150 mHz) a vysokém (mezi 150 a 400 mHz). Sledovali jsme změnu těchto parametrů během v závislosti na frekvenci dýchání.

(a)


(b)Obrázek 4. Výkonové spektrum HRV při (a) neřízeném dýchání a (b) testu hlubokého dýchání

4. Výsledky

Mladí a zdraví jedinci, kteří netrénují řízenou respiraci, měli při spontánním dýchání prakticky periodický průběh HRV, nekoreloval však přesně s průběhem dýchání. Stále však lze pozorovat vliv dýchání na srdeční variabilitu. Pří řízeném dýchání se změnili některé parametry HRV - NN50 se zmenšil a to jak absolutně tak i relativně, body Poincarého grafu se nashromáždili blíže k osám, tak že začalo být patrné eliptické uspořádání. Parametr RRmean se však významně nezměnil. Při změně frekvence dýchání na nepřirozenou se parametry srdeční variability zhoršili (směrem k hodnotám, které jsme pozorovali pří neřízeném dýchání). Nepřirozenou frekvenci rozumíme frekvenci bez zadržení dechu v nádechu a ve výdechu a také frekvenci 4 - 1 - 8 - 1.nic

(a)


(b)nic

Obrázek 5: Vývoj lineárních parametrů (a) RRmean, (b) NN50, popisujících srdeční variabilitu během změny frekvence dýchaní

Na spektru HRV, který byl pořízen během analýzy signálu ve frekvenční oblasti je dobře vidět rapidní zvýraznění spektrální složky odpovídající dýchání. Což odpovídá nárůstu vlivu dechu na regulaci srdeční variability. Na spektrech některých signálů jsme schopni rozeznat vyšší harmonické složky frekvence dýchání. Byla také provedená segmentace signálu dechu z termistoru, tak že jsme vždy rozdělili dech na části nádechu, výdechu a zadržení dechu v nádechu a ve výdechu. Algoritmus segmentace však nefungoval zcela spolehlivě a je potřeba ho do budoucna vylepšit. I přes špatnou segmentaci však můžeme pozorovat určité zpoždění signálu HRV za signálem dechu.

(a)


(b)Obrázek 6: (a) Ukázka segmentace signálu dechu a zobrazení jedné periody dýchání a srdeční variability, odpovídající této periodě pro sekvenci dýchání 4 - 1 - 4 - 1.

Obrázek 7: Ukázka srdeční variability pří neřízeném dýchání

Obrázek 8: Ukázka srdeční variability pří řízeném dýchání s frekvenci 5 - 1 - 5 - 1 po apnoe


Ne zcela zpracovnou části však zůstává fáze apnoe. Bylo potvrzeno, že apnoetická fáze se vyznačuje některými charakteristickými úseky. Každý z těchto úseků trval cca 20s a postupně je lze popsat následovně: úsek stejné srdeční variability (kdy HRV se některou dobu prakticky neměnilo své vlastnosti), úsek nulové variability (amplituda HRV byla nepatrná) a úsek poloviční amplitudy (kdy najednou se srdeční variabilita zase byla dobře pozorovatelná).

Obrázek 9: HRV během zadržení dechu4. Závěr

V tomto článku byly uvedené parametry závislosti srdeční variability na frekvenci dýchání. Pro detekci byly použité různé metody v časové a frekvenční oblasti. Pro jedince, kteří necvičí řízené dýchání není závislost HRV na dýchání pří přirozeném dechu úplně jednoznačná, HRV měřené osoby, která dýchání cvičila i při neřízené respiraci vykazuje hodnoty, které jsou normální pro respiraci řízenou.

5. Poděkování

Tento výzkum byl podporován z výzkumného záměru MSM6840770012 "Transdisciplinární výzkum v oblasti biomedicínského inženýrství" a z grantů GAČR102/08/Х008 !Analýza a modelování biomedicínských a řečových signálů" a SGS10/273/OHK3/3T/13 "Analýza signálů mechanické aktivity srdce"

Reference

[1] Javorka, K.; a kolektív: Variabilita frekvencie srdca, Osveta, 2008

[2] Kumari, S; Kyriacou, P.A.; Shafqat, K; Time-frequency analysis of HRV data from locally anesthetized patients. IEEE EMBS, pages 1824-1827, 2009.

[3] Hotoleanu, C; Agoston L.C.; Zdrenghea, D.; Dumitrascu, DL.; Rusu, LD; Poant. L; Heart rate variability assessment physiological and pathological aspects. IEEE, 2008.


Systém pro zpracování EEG v reálném čase

Jaromír Doležal

České vysoké učení v Praze, Fakulta elektrotechnická [email protected]

Abstrakt: Naše skupina se zabývá výzkumem v oblasti rozpoznávání pohybové aktivity v EEG za účelem konstrukce rozhraní mozek-stroj. V současnosti pracujeme na off-line klasifikaci předem separovaných realizací pohybového EEG na základě jejich časového vývoje. Nyní se soustředíme na převedení našich algoritmů pro zpracování kontinuálního EEG. Pro jejich ověření jsme navrhli modulární systém pro zpracování EEG v reálném čase. Systém je poskládán ze samostatných modulů komunikujících přes síťové rozhraní. Provádění experimentů v reálném čase navíc umožní využít zpětnou vazbu, která zesiluje požadované projevy EEG aktivity. Systém pomáhají realizovat také studenti naší katedry, kteří jsou plně integrováni do vývojového týmu. Funkčnost systému byla úspěšně ověřena v experimentech s alfa aktivitou.

1. Úvod

Koncept rozhraní Brain Computer Interface (BCI) umožňuje přímou komunikaci mozek-stroj. V našich pracích se zabýváme rozpoznáváním pohybové aktivity na základě časového vývoje EEG signálu [1]. Mezi výhody pohybové aktivity v obecné rovině patří především přirozené ovládání systému, jelikož touto cestou obvykle ovládáme naše okolí. BCI založený na pohybové aktivitě lze využít pro náhradu ztracených motorických funkcí pomocí protéz nebo pro jejich rychlejší obnovu rehabilitací, například po mrtvici [2]. Experimenty klasifikace pohybové aktivity jsme dosud prováděli pouze off-line, tedy na předem nahraných datech a ručně separovaných realizacích [1]. Nyní se soustředíme na převedení algoritmů pro zpracování kontinuálního EEG. Pro ověření jejich funkce jsme navrhli systém pro zpracování EEG v reálném čase. Provádění experimentů v reálném čase navíc umožní využít zpětnou vazbu, která zesiluje požadované projevy EEG aktivity [2] a tím zlepšuje celkovou funkci systému. 2. Návrh systému

Protože naše algoritmy jsou stále ve vývoji, rozhodli jsme se celý systém navrhnout modulárně a jednotlivé moduly spojovat přes síťové rozhraní. Systém jsme od začátku koncipovali jako distribuovaný a otevřený. Pro implementaci celého systému byl zvolen jazyk JAVA díky jeho nezávislosti na operačním systému a hardware.

2.1 Týmová práce

Projekt je implementován ve spolupráci se studenty, proto byla nejdříve navržena specifikace systému [3] která slouží jako zadání pro implementaci. Pro kontrolu implementace byla vyvinuta metodika automatického testování. Moduly jsou průběžně testovány a nové moduly jsou do systému integrovány co nejdříve. Týmová práce je také distribuovaná, proto jsme pro její organizaci použili následující podpůrné technologie:

- Subversion server [4] pro sdílení zdrojových kódů. - Mantis bugtracker [5] pro úkolování a kontrolování dílčích kroků implementace.

22 Jaromír Doležal

2.2 Vlastnosti systému

Modulární architektura nám umožňuje jednoduše měnit a přidávat další algoritmy a funkce. Systém je díky tomu možno jednoduše konfigurovat pro různé třídy experimentů. Pro konkrétní experiment se systém poskládá na míru z odpovídajících modulů. Díky propojení modulů pomocí síťového rozhraní je náš systém také distribuovaný. Distribuovaná architektura zjednodušuje provádění experimentů na různých pracovištích, kdy stačí na snímací stanici spustit jen odpovídající modul zachytávání dat z EEG přístroje. Vlastní výpočetně náročné zpracování pak může probíhat na vyhrazeném zařízení či přenosném počítači. Distribuování systému mimo jiné také umožňuje efektivní rozložení výpočtů na nejen na pracovní stanice ale i jejich jednotlivé procesory. Pro prezentační vrstvu lze využít i tenkého klienta, jako je např. chytrý mobilní telefon či jiné zařízení pro integraci s assistivními technologiemi. Otevřenost našeho systému spočívám v tom, že veškeré nastavení je ukládáno v konfiguračních textových souborech. Ve zdrojové kódech jsou implementovány algoritmy obecně a jejich konkrétní použití je definováno až při konfiguraci systému pro zamýšlený experiment.

2.3 Síťový přenos EEG

Telemedicínský přenos EEG dat po síti je podobný přenosu multimédií jako je např. přenos hlasu a video konference. Při práci v reálném čase je třeba zajistit co nejmenší odezvu (latenci) i za cenu ztráty či nepoužití pozdě doručených dat. Proto jsme pro komunikaci mezi moduly vybrali protokol RTP (Real Time tranfer Protocol), který dosahuje nejnižší odezvy díky prioritě na síťových prvcích kde datagramy RTP předbíhají ostatní provoz. Pro podporu protokolu RTP byla použita knihovna LibRTP [6]. Pro naše potřeby jsme navrhli vlastní strukturu přenášených dat. Rozšíření bylo provedeno tak aby pakety s daty bylo možné normálně zpracovávat libovolným softwarem či hardwarem podporující přenos pomocí protokolu RTP. Obsah paketu je rozdělen do tří hlavních částí:

- Povinná hlavička obsahuje základní parametry společné pro všechny pakety, je to například vzorkovací frekvence, počet kanálů, počet vzorků a číslo bloku dat.

- Volitelná hlavička obsahuje další parametry které se v paketu můžou ale nemusí vyskytovat. Tyto parametry jsou konfigurovatelné a jejich výčet je uložen v textovém souboru konfigurace. Přenášené parametry jsou jednak příkazy pro moduly, stavy systému a další zprávy.

- Datová část obsahuje vlastní data – vzorky EEG signálu a časové průběhy vypočtených parametrizací.

2.4 Zachytávání dat z přístroje BIOPAC

První verze systému byla navržena pro práci s přístrojem AlienEEG v laboratoři evokovaných potenciálů na lékařské fakultě Karlovy University v Hradci Králové. Pro potřeby testování a ladění systému v prostorách naší laboratoře bylo třeba získávat v reálném čase data z dostupného přístroje BIOPAC MP 35. Jelikož originální software dodávaný k zařízení předávání dat v reálném čase neumožňuje, přistoupili jsme k vývoji systému pro zachytávání komunikace na USB. Navržené řešení je schematicky zobrazeno na obrázku 1.

Jaromír Doležal 23

Obrázek 1: Schéma zachytávání dat z přístroje BIOPAC

Námi implementované řešení má tři části:

- Pro nastavení parametrů přístroje a jeho ovládání jsme využili stávající software dodávaný k přístroji. Zařízení se nakonfiguruje pro konkrétní experiment a jeho nastavení se uloží do standardního souboru záznamu. Zařízení se ovládá originálním SW takže je k dispozici jeho plná funkčnost.

- Transparentní filtr na úrovni jádra Windows se zavede do cesty datové komunikace na sběrnici USB mezi stávajícím software a přístrojem a na základě dekódovaného protokolu vybírá z komunikace vzorky signálu a dále je předává modulu pro zachytávání dat BiopacRTPBridge.

- Modul BiopacRTPBridge načítá soubory s konfigurací přístroje a na jejich základě pak ze zachytávaných vzorků vytváří a odesílá datové pakety k dalšímu zpracování.

Specifika zachytávání komunikace na sběrnici USB a dekódovaný protokol lze najít v samostatné výzkumné zprávě [7]. 2.5 Konfigurace systému

Konfigurace systému je uložena v samostatných textových souborech. Konfigurace má čtyři hlavní části:

- Soubor definic parametrů obsahuje výčet parametrů přenášených po síti a rozšiřuje definici protokolu pro přenos dat.

- Soubor konfigurace parametrů obsahuje výčet hodnot používaných algoritmy systému jako jsou například konkrétní hodnoty použitých filtrů.

- Soubory konfigurace experimentů obsahují informace které algoritmy se mají použít pro daný experiment.

- Spouštěcí skripty experimentů pak obsahují informace o propojení modulů, jejich rozmístění v síti a zajišťují spuštění modulů s odpovídající konfigurací.


Tento přístup umožňuje jednoduché rozšiřování systému i protokolu pro přenos dat bez úprav zdrojových kódů. 2.6 Architektura modulů

Moduly jsou implementovány s využitím vláken a paralelizmu, což umožňuje využít podporu paralelismu poskytovanou moderním vícejádrovými procesory a technologie hyper-threading. Společné funkce modulů jsou implementovány v abstraktních třídách balíku. Jedná se o jednotlivá výpočetní vlákna, jejich základní funkce a metody pro mezi-vláknovou asynchronní synchronizaci. Použity jsou čtyři hlavní vlákna:

- Vlákno přijímaní dat spravuje příchozí komunikaci, provádí třídění paketů, detekci ztracených paketů a ukládání přijatých paketů do fronty pro zpracování. Samostatné vlákno pro přijímaní dat je nutné aby nedocházelo ke ztrátám paketů, protože protokol RTP sám doručení dat nezajišťuje.

- Vlákno zpracování provádí dekódování paketů a tvorbu všech odpovídajících datových objektů obsažených v paketech. V tomto vláknu je pak implementován vlastní algoritmu který má modul provádět.

- Vlákno logování zajišťuje ukládání veškeré komunikace do souborů za účelem automatického testování a ladění systému.

- Vlákno odesílání zajišťuje tvorbu datových paketů a jejich rozesílání na zadané adresy.

Tento přístup umožňuje jednoduše implementovat další moduly, kdy se v novém modulu aktivují požadovaná vlákna a implementuje se jen vlastní funkčnost (algoritmus). 2.7 Řetězec BCI

Návrh řetězce pro zpracování v reálném čase úzce reflektuje architekturu BCI systému, kde jednotlivým blokům BCI systému odpovídají samostatné moduly v řetězci pro zpracování EEG. Typické zapojení modulů je schematicky naznačeno na obrázku 2. Ve schématu nejsou uvedeny moduly určené pro testování systému – modul pro ukládání surových dat do souborů a modul generování dat ze souborů. Detailní informace o všech modulech lze získat ve specifikaci systému [3].

Obrázek 2: Typické zapojení modulů při experimentu se zpětnou vazbou.


2.7.1 Moduly pro zachytávání dat V současné době jsou podporovány dva EEG přístroje. Moduly pro zachytávání dat (zeleně na obrázku 2) slouží zachycení vzorků signálu, tvorbě datových paketů RTP a jejich odeslání k dalšímu zpracování. Tyto moduly je třeba instalovat na pracovní stanici ke které je připojeno zařízení pro snímání dat. 2.7.2 Ovládací modul Ovládací modul slouží k ovládání systému za běhu experimentátorem. Jedná se především o přepínání systému mezi stavy učení/klasifikace a vysílání příkazů pro další modulu jako je například resetování či načítání a ukládání jejich stavů (např. naučeného klasifikátoru). Modul dále slouží k úpravě parametrů dalších modulů za běhu jako je například rychlosti učení nebo hodnot použitých pro prahování. Modul je zařazen v řetězci pro zpracování jako první, tak aby se parametry a příkazy vložené do paketů dostali ke všem zbývajícím modulům. 2.7.3 Detekce aktivity Modul detekce aktivity slouží k trénování systému v asynchronním režimu experimentu, tj. kdy si experimentální osoba sama volí čas pro provádění detekovaných akcí. Pro pohybovou aktivity to znamená zpracování výstupů EMG kanálů. V případě synchronního přístupu by funkce detekce aktivity nahradil modul zajištující zvukové či obrazové znamení. 2.7.4 Povrchové filtrace Protože EEG signál prochází přes různé vrstvy, mozkomíšní mok, lebku, je aktivita na povrchu hlavy rozmazána a úzce lokalizovaná EEG aktivita je díku tomu méně zřejmá. Modul povrchové filtrace proto implementuje diskrétní laplacovské filtry, které umožní požadovanou aktivitu zvýraznit. 2.7.5 Výpočet příznaků Modul výpočtu příznaků provádí výpočet příznaků ze signálu pro klasifikaci. Pro naše první experimenty byly implementovány algoritmy FIR filtrů a navrženy filtry pro detekci ERS u pohybové aktivity a alfa aktivity u experimentů s otevřenýma a zavřenýma očima. 2.7.6 Klasifikace Modul klasifikace provádí klasifikaci vypočtených příznaků. Modul má implementované funkce učení a vlastní klasifikaci. Natrénovaný klasifikátor lze uložit do souboru a využít v dalších experimentech. Pro první experimenty byl implementován perceptronový algoritmus. 2.7.7 Visualizační modul Visualizační modul zajišťuje na základě výsledků klasifikace zpětnou vazbu pro experimentální osobu. V prvních experimentech byl použit vznášející se míček a rozšiřující se obdélník. 2.7.8 Monitorovací modul Monitorovací modul slouží ke sledování funkce systému za běhu a ladění systému. Modul umožňuje zobrazit všechny druhy dat v grafickém prostředí (obrázek 3). K dispozic jsou údaje o experimentální osobě, zapojení elektrod a stavu systému. Příchozí data jsou časově synchronizována a tak je možné hned vidět nad sebou průběhy surového EEG, filtrovaného EEG a vypočtených příznaků pro klasifikaci. Události jsou zobrazeny barevnými svislými


čarami, jedná se například o nedoručená a nahrazená data (červeně) a o detekovanou či klasifikovanou aktivitu (zeleně).

Obrázek 3: Grafické rozhraní monitorovacího modulu.

V textovém okně jsou pak zobrazeny všechny zprávy mezimodulové komunikace. 3. Experimenty a další kroky

Funkčnost systému byla ověřena při experimentech na našem pracovišti. Snímek pracoviště je na obrázku 4. Pro experiment s alfa aktivitou bylo použito zapojení modulů schematicky uvedené na obrázku 2.

Obrázek 4: Snímek pracoviště při experimentu s alfa aktivitou.


Systém se úspěšně naučil klasifikovat alfa aktivitu při otevřených a zavřených očích. Největším problémem se kterým se potýkáme na našem pracovišti je malý odstup signál-šum, který komplikuje klasifikace pohybové aktivity. Pohybová aktivita u extezních a flexních pohybů ukazováčku kterou se zabýváme [1] je výrazně slabší a více lokalizovaná než alfa aktivita a proto je v našich podmínkách hůře detekovatelná. V dalším experimentech se proto zaměříme na aktivitu při déle představovaných pohybech větších části těla tak, aby bylo možné využít zpětné vazby. Experimenty na našem pracovišti poslouží k odladění systému, metodiky provádění experimentů a nastavení pro jednotlivé experimenty. Takto připravené experimentální procedury poté zopakujeme v laboratoři evokovaných potenciálů na lékařské fakultě Karlovy University v Hradci Králové kde je k dispozici kvalitnější vícekanálový EEG přístroj. 4. Závěry

Byl navržen otevřený modulární systém pro zpracování EEG. Jednotlivé moduly odpovídají částem BCI systému, systém lze poskládat a nakonfigurovat na míru konkrétnímu experimentu. Moduly komunikují pomocí síťového rozhraní a systém není závislý na konkrétním EEG přístroji či operačním systému. Navržená architektura umožňuje volně distribuovat části systému po síti a na plno využít paralelismu a vícejádrových procesorů. Byly implementovány základní algoritmy pro jednotlivé moduly a funkčnost systému byla úspěšně ověřena při skutečných experimentech s alfa aktivitou v reálném čase. 5. Poděkování

Tento výzkum byl podporován z grantu GAČR č. 102/08/H008 "Analýza a modelování biologických a řečových signálů" a výzkumného záměru MŠMT MSM6840770012 "Transdisciplinární výzkum v oblasti biomedicínského inženýrství 2" a grantem studentské grantové soutěže ČVUT č. SGS 10/178/OHK3/2T/13. Poděkování patří také doc. Janu Kremláčkovi za umožnění nahrávání a spolupráci při experimentech v laboratoři evokovaných potenciálů na lékařské fakultě Karlovy University v Hradci Králové.

Reference

[1] J. Doležal, BCI založený na manifestaci pohybové aktivity v EEG II, semináře katedry teorie obvodů, analýza a zpracování řečových a biologických signálů - sborník prací 2009, str. 38-43, 2009.

[2] Chista Neuper, Feedback-regulated mental imagery in BCI applications and NIRS signals, BBCI Workshop 2009, Advances in Neurotechnology, Berlin, 2009.

[3] Doležal J, Šťastný J. Real Time EEG Processing software, system specification. ČVUT, FEL, 2010.

[4] Software pro sdílení zdrojových kódů Subversion server. http://subversion.tigris.org/

[5] Software pro úkolování a řízení týmu Mantis bugtracker. http://www.mantisbt.org/

[6] Knihovna LibRTP library. http://sourceforge.net/projects/jlibrtp/

[7] Jan Kubák, Jaromír Doležal, zachytávání dat z přístroje BIOPAC MP35, výzkumná zpráva #Z10-01, ČVUT, FEL, 2010


Možnosti detekce SOZ v intrakraniálním EEG signálu

Radek Glajcar, Radek Janča

České vysoké učení v Praze, Fakulta elektrotechnická [email protected], [email protected]

Abstrakt: Příspěvek představuje možnosti detekce počátku epileptických záchvatů za účelem lokalizace epileptoformních ložisek z intrakraniálního EEG. Shrnuje několik metod postavených na extrakci parametrů z časově-frekvenční analýzy, transformaci signálů neuronovou sítí, korelační analýze a hlouběji rozvádí metodu na bázi multidimenzionálních autoregresních modelů. Dále jsou představeny prvotní výsledky zmíněné metody implementované v prostředí Matlab.

1. Úvod

U neurologických onemocnění mozku je jedním ze standardních vyšetření analýza skalpového EEG. Epilepsie patří mezi záchvatovité neurologické onemocnění a jedním z příznaků je patologický nález v těchto signálech mozkové aktivity. Hodnocení patologií je ze strany lékařů více či méně subjektivní a v běžné praxi zahrnuje nalezení artefaktů v signálu, jako jsou například výkyvy amplitudy, zvrat fáze, hroty či aktivita s vyšší frekvencí. Velikost amplitudy, čas výskytu a četnost těchto artefaktů vede v kombinaci s dalšími vyšetřeními k lokalizaci pravděpodobných epileptoformních center. U pacientů nereagujících na medikamentózní léčbu se mnohdy přistupuje k chirurgické resekci postižených oblastí. Nicméně ne vždy lze lokalizovat ložiska neinvazivními metodami před samotným zákrokem. Skalpové EEG je málo odolné vůči rušení, lebka funguje jako dolní frekvenční propust a způsobuje prostorový rozptyl signálů. Proto se přistupuje k dvouplánovým operacím. Při první operaci se přímo na povrch mozku umisťují elektrody, čímž se eliminují nedostatky skalpového EEG. Po týdenním monitoringu a analýze intrakraniálního EEG a zpřesnění lokalizace probíhá druhá resekující operace. Přesné určení epileptoformních center a jejich odstranění má zásadní vliv na úspěšnost léčby a na následné postižení pacienta. Spolupráce s Klinikou dětské neurologie FN Motol má za cíl vyvinout metodu, která objektivně napomůže lékařům identifikovat patologie v intrakraniálních EEG signálech a tím i přispět k přesnější lokalizaci ložisek. 2. Metody detekce SOZ

Určení oblasti začátku záchvatu (SOZ – Seizure Onset Zone) je důležité k detekci samotného záchvatu výstražnými systémy, primárně však slouží k určení ložiska a směru šíření. Detekce SOZ není omezena pouze na počátky záchvatů, ale k lokalizaci mohou posloužit i např. vysokofrekvenční oscilace (ripples) v pásmu 80-250 Hz, vyskytující se výhradně v místech SOZ. Podstatou většiny metod je časově-frekvenční analýza signálů se společnými prvky extrakce časových a spektrálních příznaků, častěji však z frekvenčního spektra, a jejich následné klasifikace. 2.1 Časově-frekvenční analýza

T. Tzallase a M. G. Tsipourase uvádějí přehled bázových funkcí časově-frekvenčních distribucí k výpočtu spektrální výkonové hustoty [1]. Z výpočtů jsou extrahovány příznaky klasifikované dopřednou neuronovou sítí. Pro porovnání autoři uvádějí ještě další

Radek Janča 29

klasifikátory a sice bayessovský, k-NN (k-nearest neighbour) a rozhodovací stromy. Klasifikátor dělí vstupní data do dvou tříd, a to do třídy příznaků s nebo bez záchvatu. Autoři se primárně věnují problematice přesnosti klasifikace, nicméně z výsledků lze určit čas, kdy dochází k počátku záchvatu.

Metoda autorů v čele s Alexandrem M. Chanem rovněž zjišťuje čas počátku záchvatů [2]. Signál je klouzavě segmentován, pro každý segment je periodogramem odhadnuta výkonová spektrální hustota. Z té je vypočten příznakový vektor integrací ve frekvenčních pásmech korespondujících s frekvenčními rozsahy vln běžného EEG (pásmo θ, α, β, γ-1, γ-2). Příznakové vektory jsou klasifikovány klasifikátorem SVM (Support Vector Machine) natrénovaným daty obsahující počátek záchvatu nebo naopak daty bez počátku záchvatu. K přesné detekci času SOZ se z výsledných hodnot klasifikace používá regresní analýzy.

Detektor výstražného systému pro nemocniční personál je sestaven z dvojstupňové filtrace [3]. První filtrace je v pásmu 5 - 45 Hz zajištěna pásmovou propustí FIR filtru sestaveného na základě koeficientů vlnky třetí úrovně z rodiny DAUB4 vlnkové transformace. Kvadrát výsledného signálu je následně filtrován tentokrát klouzavým mediánovým filtrem, který je odolný proti statisticky odchýleným vzorkům. Mediánový filtr zachovává rozdíly ve skokovém přechodu v signálech, čehož se využívá při detekci přechodu signálu bez záchvatu a se záchvatem. Různou úpravou výstupu mediánového filtru se vytvoří „pozadí“ a „popředí“ signálu z prvního filtru. Prahovaný poměr „popředí“ a „pozadí“ spouští varovný signál.

Vývoj detekčního algoritmu je v případě autorů Shoeb a spol. motivován zkrácením doby podání radiofarmaka od začátku záchvatu při vyšetření jednofotonovou emisní tomografií (SPECT). Příznakový vektor je tvořen energií čtyř různých časových měřítek signálu po průchodu bankou filtrů vlnkové transformace v rozsahu 0,5 – 25 Hz [4]. Každý kanál je segmentován dvousekundovým oknem, pro které jsou spočteny 4 hodnoty. Pro dané časové okno je sestaven příznakový vektor skládající se z příspěvků 4 hodnot každého kanálu. Příznakový vektor je klasifikován pomocí SVM do tříd „záchvat / bez záchvatu“. Pokud jsou tři po sobě jdoucí příznakové vektory klasifikovány do třídy „záchvat“, je aplikováno radiofarmakum.

Gotman a spol. využívají modifikovaný klasifikátor NN (nearest neighbour) [5]. Jako příznakový vektor je využíván parametrizovaný úsek 2,5 sekundového signálu. Šest prvků příznakového vektoru reprezentují parametry: průměrná amplituda vlny, průměrná doba trvání vlny, variační koeficient doby trvání vlny, dominantní frekvence, průměrný výkon a lokalizační příznak. 2.2 Využití neuronové sítě

Neuronová síť je využita v práci R. Batese a kol. [6]. Jedná se o třívrstvou rekurentní síť s architekturou 5-10-5. Na vstupy sítě jsou přiváděny signály kanálů intrakraniálního EEG a síť je trénována na cílovou hodnotu nula a jedna dle signálu se záchvatem a bez záchvatu. Při testování jsou výstupy převedeny regresní analýzou na koeficienty R. Pokud se R blíží nule, je signál označen s největší pravděpodobností jako signál bez záchvatu a naopak. 2.3 Korelační analýza

V současné době se ukazuje nadějný směr detekce záchvatu pomocí korelační analýzy. Např. K. Schindler a kol. filtroval multikanálové signály v pásmu 80 – 200 Hz [7]. Ze signálů vypočetli kovarianční matici s celkovou silou korelace (TSC - Total Correlation Strength).

30 Radek Janča

Z matic lze vysledovat nárůst korelace při počátku záchvatu, nicméně metoda nebyla ověřena na dostatečném počtu pacientů. 2.4 Metody na bázi autoregresních (AR) modelů

B. Swiderski a spol. využili k analýze EEG signálů směrovou přenosovou funkci (DTF-Directed Transfer Function) [8]. Popisovaná parametrická metoda je založena na multidimenzionálních AR modelech popisujících směr a sílu toku elektrické aktivity mozku. Signály kanálů jsou segmentovány po 2 sekundách s 50% překryvem. Mezi časově odpovídajícími segmenty z různých kanálů je vypočtena DTF, která tvoří příznakový vektor každého segmentu. Příznakový vektor je klasifikován pomocí OC-SVM (One Class-SVM) do třídy „záchvat / bez záchvatu“. Ne všechny klasifikované segmenty spadají do správné třídy, proto se z několika po sobě jdoucích klasifikovaných segmentech počítá odhad apriorní pravděpodobnosti výskytu záchvatu (p-estimate). 3. Implementace metody s DTF

Implementace v prostředí Matlab vychází z SOZ metody postavené na bázi AR modelů popisované v 2.4 [8]. Z m-kanálového signálu x, je vypočten multidimenzionální parametrický AR model,

i

N

iiti exA =∑

=−

0

(1)

A i je matice koeficientů modelu, xt je vektor vytvořený z m EEG signálů v čase t, N odpovídá řádu modelu a et je vektor bílého nekorelovaného šumu. Řád AR modelu je počítán dle informačního Akaike kritéria individuálně pro každého pacienta [9]. Vzhledem ke krátkodobé stacionaritě signálů je provedena segmentace oknem odpovídající jedné sekundě záznamu s překryvem 75 %. Matice A i je časově závislá a může být odhadnuta pomocí multikanálové kovarianční analýzy. Použitím Fourierovy transformace se definuje matice přenosové funkce H(f),

1

0

)2exp()(−

=

−= ∑

N

i si f

fijf πAH (2)

kde f je frekvence, fs je vzorkovací kmitočet a 1−=j . Na základě DTF )(2 fijγ popisuje

„tok“ z kanálu j do kanálu i,

∑=

= m

kik

ij

ij

fH

fHf

1

2

2

2

)(

)()(γ (3)

kde Hij(f) je element matice H(f) v řádku i a sloupci j. )(2 fijγ je tedy normalizovaná DTF.

Původní signál je filtrován do 8 subpásem FIR filtrem řádu sf3 , pásma jsou popsána

v tabulce 1.

Označení δδδδ θθθθ αααα1 αααα2 ββββ1 ββββ2 γγγγ1 γγγγ2 Pásmo [Hz] 2-4 4-8 8-10 10-12 12-18 18-25 25-48 52-85

Tabulka 1.: Subpásma EEG signálu

Radek Janča 31

Pro každé subpásmo je definován střední „tok“ ( ))()( kijij fmeanpásmo γγ = z kanálu j do

kanálu i, kde fk značí diskrétní frekvence v daném subpásmu. „Tok“ aktivity z kanálu j do všech kanálů je definován rovnicí 4.

∑≠=−

=m

iji

ijj pásmom

DF1

)(1

1 γ (4)

[ ])( ), DF( ), DF( ), DF( ), DF( ), DF( ), DF( ), DF(DF jjjjjjjjj 212121 γγββααθδ=y (5)

Vzniklý vektor yj reprezentuje parametrizovaný segment kanálu j. Vektory yj jsou klasifikovány pomocí OC-SVM do třídy „baseline / záchvat“ (+1 / −1). Klasifikátor je trénován množinou parametrizovaných segmentů signálů bez záchvatové aktivity; pro potřeby použitého klasifikátoru bylo nutno snížit dimenzi yj pomocí PCA (Principal Component Analysis). OC-SVM hledá při trénování rozhodovací vektor nadroviny tak, aby 10 % dat leželo mimo třídu „baseline“. To má za následek 10% chybu klasifikace na trénovací množině, ale zato bez potřeby definovat trénovací data třídy „záchvat“. Pro N po sobě jdoucích segmentů je vypočtena apriorní pravděpodobnost P četnosti výskytu klasifikací „záchvat“ N−, tzv. p-odhad. Nárůst pravděpodobnosti P na jednotlivých kanálech indikuje počátek záchvatu a tím i ložisko.

N

NP

−

= (6)

Výsledná data jsou vizualizována do kortikální mapy přibližně odpovídající skutečnému rozložení elektrod na povrchu mozku. Hodnota P na jednotlivých elektrodách je reprezentována barevnou škálou. Časový vývoj mozkové aktivity lze dobře pozorovat ve videosekvenci.

5.51 5.512 5.514 5.516

x 104

0

0.2

0.4

0.6

0.8

1

P

t [s]

low activity

high activity

Obrázek 1.: Příklad vizualizace p-odhadu

v kortikální mapě barevnou škálou. Obrázek 2.: Příklad rozdílu průběhu

p-odhadu mezi dvěma kanály v čase začátku záchvatu.

32 Radek Janča

4 Závěr

Implementovaná metoda využívající DTF sleduje na rozdíl od jiných algoritmů změny napříč všemi kanály. Dává tak komplexnější informaci o začátku a průběhu epileptického záchvatu. V současné době probíhá zkoumání vlivu velikosti a překryvu segmentace na lokalizaci ložiska a počátku záchvatu. Výsledky analýzy signálů zkoumaných pacientů ukazují na velkou shodu s nálezy lékařů. Nicméně k velmi malému počtu pacientů s různorodou anamnézou je předčasné hodnotit kvalitu metody. V budoucnu se plánuje porovnání výsledků pomocí analýzy jinými metodami. Modifikací metody a za použití jiných parametrizací bude snaha detekovat i rychlejší změny v EEG signálech reprezentované např. mezizáchvatovými (interiktálními) výboji. 5 Poděkování

Tato práce je podporována granty IGA NT11460-4/2010 Komplexní analýza intrakraniálního EEG záznamu a identifikace epileptogenní zóny u pacientů s nelezionální farmakorezistentní epilepsií, SGS 10/272/OHK4/3T/13 Analýza intrakraniálních EEG záznamů výzkumný záměr MSM6840770012 Transdisciplinární výzkum v biomedicínském inženýrství. Reference

[1] A. T. Tzallas, M. G. Tsipouras, and D. I. Fotiadis. Epileptic seizure detection in EEGs using time-frequency analysis. IEEE transactions on information technology in biomedicine, 13(5): 703–710, 2009.

[2] A. M. Chan, F. T. Sun, E. H. Boto, and B. M. Wingeier. Automated seizure onset detection for accurate onset time determination in intracranial EEG. Clinical Neuro-physiology, 119: 2687–2696, 2008.

[3] I. Osorio, M. Frei, and et al. Real-time automated detection and quantitative analysis of seizures and short-term prediction of clinical onset. Epilepsia, 39(6): 615–627, 1998.

[4] A. Shoeb, H. Edwards, J. Connolly, B. Bourgeois, and et al. Patient-specific seizure onset detection. Epilepsy & behaviour, 5: 483–498, 2004.

[5] H. Qu and J. Gotman. A patient-specific algorithm for the detection of seizure onset in long-term EEG monitoring: Possible use a warning device. IEEE transactions on biomedical engineering, 44(2): 115–122, 1997.

[6] R. R. Bates, S. Mingui, M. L. Scheue, and R. J. Sclabassi. Detection of seizure foci by recurrent neural networks. Proceeings of the 22nd Annual EMBS international conference, pages 1377–1379, 2000.

[7] K. Schindler, F. Amor, H. Gast, and et al. Peri-ictal correlation dynamics of high frequency (80-200 Hz) intracranial EEG. Epilepsy research, 2009. doi:10.1016/j.eplepsyres.2009.11.006, stav z 18. 10. 2010.

[8] B. Swiderski, S. Osowski, A. Cichocki, and A. Rysz. Single-class SVM nad directed transfer function approach to the localization of the region containing epileptic focus. Neurocomputing, 72: 1575–1583, 2009.

[9] R. Čmejla. Kritéria pro určení řádu AR modelu při zpracování řečových signálů.

Akustické listy, 22: 4–7, 2000.

Radek Janča 33

Odhad logopedického věku z řeči dítěte

Jan Janda

České vysoké učení technické v Praze, Fakulta elektrotechnická[email protected]

Abstrakt: Tato práce ukazuje a hodnotí vybrané akusticko-fonetické a logo-pedické parametry dětské řeči. Tyto parametry byly získány z analýz dětskýchpromluv v nově pořízené databázi a jsou porovnány podle míry jejich věkovézávislosti.

1. Úvod

Vývoj řeči v dětství je po akustické stránce zapříčiněn částečně růstem a anatomickoupřestavbou řečového traktu. Od narození do dospělosti se délka vokálního traktu prodloužípřibližně na dvojnásobek délky. Významně se však také mění geometrické proporce jed-notlivých měkkých a tvrdých tkání relativně k délce vokálního traktu. Jedná se napříklado postupné zakřivení vokálního traktu do pravého úhlu v oblasti nosohltanu, sestoupeníhrtanu, jazylky a příklopky hrtanové a pokles zadní části jazyka tak, aby tvořil přednístěnu hltanu [6].V důsledku těchto anatomických změn rostou jednotlivé kostní a měkko-tkáňové strukturyv oblasti ústní dutiny a hltanu různou rychlostí a svých dospělých rozměrů dosahují vširokém rozmezí od 7 do 18 let věku.V práci [6] byla pomocí moderních zobrazovacích metod (zejména magnetické rezonance)provedena kvantitativní měření anatomických charakteristik jednotlivých částí vokálníhotraktu během prvních 20 let života. Délkou vokálního traktu se zde rozumí délka křivkyv mediální rovině vedené od středu hlasivek po její průsečík s tečnou rtů. Během vývojese délka této křivky zvětší z přibližně 8 cm u novorozenců po asi 18 cm u dospělých mužů(obr. 2).Práce dále předkládá věkové závislosti dalších anatomických charakteristik. Jedná se otloušťku maxilárního a mandibulárního rtu, délku tvrdého a měkkého patra, délku ja-zyka, délku mandibuly, a vzdálenost hrtanu, příklopky hrtanové a jazylky od spina nasa-lis posterior. Tyto charakteristiky vykazují podobný rostoucí trend, liší se však ve věku,ve kterém dochází k největšímu růstu. Každé charakteristice je pak přiřazen po částechlineární model růstu. Ze závěrů anatomických měření můžeme předpokládat některé akus-tické důsledky. Jednak s rostoucí délkou celého vokálního traktu můžeme předpokládatočekávaný pokles formantů [4, 2], dále díky absenci pohlavního dimorfismu v uvedenýchcharakteristikách do 6 let věku můžeme do jisté míry očekávat i pohlavní nezávislost sou-visejících akustických parametrů. Také vedle zvyšující se artikulační zručnosti bude mítvliv na koartikulaci i zvětšující se prostor pro pohyb jazyka.

34 Jan Janda

Obrázek 1: Výstup magnetické rezonance v mediální rovině. Vlevo 7 měsíční dívka, vpravodospělá žena. Převzato z [6]

Obrázek 2: Věková závislost délky vokálního traktu dětí a dospělých (prázdné trojúhelníkydolů pro muže a stínované nahoru pro ženy). Převzato z [6]

Jan Janda 35

Vzhledem k tomu, že se během dospívání mění nejen rozměry, ale i vzájemné uspořádáníjednotlivých částí vokálního traktu, můžeme předpokládat, že i další foneticko-akusticképarametry dětské řeči budou vykazovat věkovou závislost.

2. Databáze

V rámci této studie byla na základě dosavadních poznatků a zkušeností nově pořízenadatabáze dětských promluv. Obsah je již cíleně navržen tak, aby se v databázi vysky-tovaly především typy promluv, na nichž se dařilo věkově závislé parametry naměřit.Databáze je oproti té předchozí jednotná pro děti předškolního a školního věku. Slovaděti vyslovují pouze podle toho, co vidí na obrázku, a nikoliv už pouhým opakováním sly-šeného. Z izolovaných slov byla po konzultaci s Foniatrickou klinikou 1.LF UK vybránaslova: babička, máma, čokoláda, sluníčko, popelnice, košile, silnice, Rákosníček, hambur-ger, velryba, ucho, ježek, ředkvičky, fotbalista. Pro posouzení základních fonetických para-metrů byly do databáze zařazeny prodloužené fonace hlásek /a/, /e/, /i/, /o/, /u/,/s/, /ss/. Dále byla zařazena říkanka „Ententýky dva špalíky. . .ÿ. Pro následné měřenídiadochokinetických parametrů se mluvčí měli pokusit o co nejrychlejší opakování sekvencí/pa/-/ta/-/ka/ a /ba/-/da/-/ga/. Každé nahrávací sezení pak dítě ukončilo krátkýmspontánním projevem – vyprávěním příběhu podle předložených obrázků.Databáze v tomto okamžiku obsahuje promluvy od 250 dětí ve věku od 4 do 15 let.

3. Věkově závislé akustické charakteristiky

3.1. Analýza samohlásek

3.1.1. Základní frekvence hlasivkového tónu F0

závisí na velikosti hrtanu a délce hlasivek. Jedná se o nejčastěji uváděnou charakteristikuv souvislosti s věkem člověka. Nabývá hodnot od 500Hz u nejvyšších dětských hlasů a svěkem může u mužů klesnout až na hodnoty kolem 80Hz. Analýza základního hlasivkovéhotónu byla provedena na prodloužené fonaci (cca 5 s) jednotlivých samohlásek.Analýza byla provedena autokorelační metodou v programu Praat v. 5.0.15 s parametrytime step=0.0, pitch floor=75Hz a pitch ceiling=600Hz. Výsledné hodnoty byly ověřenyv programu Wavesurfer v. 1.8.5 a případně ručně modifikovány.Pro statistické potvrzení věkové závislosti F0 uvažujme nulovou hypotézu H0, která tutozávislost popírá.H0 můžeme zamítnout na základě výsledku t-testu pro korelovaná měření.V našem případě lze H0 zamítnout na hladině p < 0, 001, n = 250. Sílu korelace můžemevyjádřit Pearsonovým korelačním koeficientem r. Pro věkovou závislost F0 samohlásky/a/ dostaneme r = −0, 57, tedy středně silnou, uspokojivou korelaci.Na obrázku 3 můžeme vidět průměrné F0 jednotlivých věkových skupin. Pro znázorněnívlivu mutace u chlapců jsou vyneseny zvlášť věkové závislosti F0 pro chlapce a pro dívky.

3.1.2. Formanty F1, F2

Věková závislost je u formantů méně zřetelná než u F0. Pro F2 pro samohlásku /a/naměřen korelační koeficient r = −0, 40. Obrázek 4 ukazuje posun a mírné natočenívokalického trojúhelníka /a/-/i/-/u/.

36 Jan Janda

2 4 6 8 10 12 14 1650

100

150

200

250

300F0

vek

Hz

dívky

chlapci

Obrázek 3: Věková závislost F0 pro hlásku /a/

300 400 500 600 700 800 900 1000800

1000

1200

1400

1600

1800

2000

2200

2400

2600

2800Formantove pole

F1(Hz)

F2(

Hz)

I

U A

Obrázek 4: Věková závislost vokalického trojúhelníka. Věk roste ve směru šipek od 4 do15 let.

Jan Janda 37

3.2. Analýza slov

3.2.1. Akumulovaná vzdálenost řečových parametrizací

Při posuzování srozumitelnosti (resp. patlavosti) dětské řeči použijeme vedle analyzovanépromluvy i promluvu referenční stejného obsahu, precizně vyřčenou. V matici vzdálenostíjednotlivých segmentů v prostoru dané řečové parametrizace nalezneme křivku DTW.Kumulativni vzdálenost podél křivky DTW bude značně korelovat s nesrozumitelnostizkoumané promluvy.Z provedených experimentů bylo ověřeno, že tato kumulativní vzdálenost s věkem klesáa promluvy starších dětí jsou tedy srozumitelnější. Největší věkovou závislost vykazovalakumulovaná vzdálenost v prostoru RASTA-PLP a kepstrálních LPC koeficientů (obr. 5).

4 6 8 10 12 14 160

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

vek

Cumsum DTW − slovo "Hamburger"

Cum

sum

CLPCRASTA−PLP

Obrázek 5: Věková závislost kumulované vzdálenosti DTW funkce slova hamburger vprostoru RASTA-PLP a CLPC koeficientů.

4. Přehled věkově závislých charakteristik

Tabulka 1 shrnuje doposud zkoumané akustické a logopedické charakteristiky. Jednot-livé příznaky jsou seřazeny podle míry korelace s věkem. Sloupec Ho obsahuje hodnotyhladin významnosti, na který je teoreticky možné zamítnout nulovou hypotézu o věkovénezávislosti parametru.

5. Závěr

Doposud řečové charakteristiky vykazují různě velikou závislost na věku. Nejčastěji uvá-děné charakteristiky založené na základní hlasivkové frekvenci vykazují korelaci s věkemokolo 0,66. Navržená metoda měření srozumitelnosti pomocí kumulované vzdálenosti valgoritmu DTW dokonce vykazuje korelaci s věkem až 0,72. Oproti výsledkům pořízenýmna předchozí datábazi zde nevykazují nijak významnou věkovou závislost spektrální cha-rakteristiky sykavek. Začíná se však ukazovat, že tyto charakteristiky jsou značně odlišnéu chlapců a u dívek a budou dále zkoumány. V dalším výzkumu budou analyzovány dalšípotencionálně věkově závislé parametry a z nich pak natrénován klasifikátor pro strojovéurčení věku.

38 Jan Janda

Charakteristika r H0

DTW fotbalista -0,72 0,001DTW popelnice -0,71 0,001DTW čokoláda -0,70 0,001DTW velryba -0,69 0,001DTW hamburger -0,66 0,001F0 /I/ -0,66 0,001F0 /U/ -0,65 0,001F0 /A/ -0,57 0,001DTW Rákosníček -0,53 0,001F2 /A/ -0,40 0,001F1 /U/ -0,38 0,001F1 /A/ -0,34 0,001F1 /I/ -0,30 0,001F2 /I/ 0,20 0,005F2 /U/ 0,18 0,005

Tabulka 1: Přehled věkově závislých charakteristik.

Poděkování

Tento výzkum je podporován z grantu GD102/08/H008 - Analysis and modeling biome-dical and speech signals.

Reference

[1] GEROSA M.,LEE S. et al.: Analyzing Children’s Speech: an Acoustic Study of Con-sonants and Consonant-Vowel Transition. In Proc. of the International Conferenceon Acoustic, Speech and Signal Processing,(Tolouse, France), May 2006

[2] HUBERT J. E., STATHOPOULOS E. T et al.: Formants of children, women, andmen: the effects of vocal intensity variation. Journal of the Acoustical Society ofAmerica 1999 Sep;106(3 Pt 1):1532-42.

[3] LEE S., POTAMIANOS A., NARAYANAN S.: Acoustic of children’s speech: Deve-lopmental changes of temporal and spectral parameters. Journal of the AcousticalSociety of America, pp. 1455-1468, Mar. 1999

[4] MENARD L., SCHWARTZ J.-L., Articulatory–acoustic relationships during vocaltract growth for French vowels: Analysis of real data and simulations with an articu-latory model. Journal of Phonetics Volume 35, Issue 1, January 2007, Pages 1-19

[5] POTAMIANOS, A.; NARAYANAN, S.: A review of the acoustic and linguistic pro-perties of children’s speech Multimedia Signal Processing, 2007. MMSP 2007. IEEE9th Workshop on Volume , Issue , 1-3 Oct. 2007 pp. 22 - 25

[6] VORPERIAN K. H., KENT R. D. et al.: Development of vocal tract length duringearly childhood: A magnetic resonance imagining study. Journal of the AcousticalSociety of America 117, January 2005, pp. 338-350

Jan Janda 39

Introduction to Selective Zolotarev-Cosines

Ján Janík, Miroslav Vlček+), Pavel Sovka

Czech Technical University in Prague, Faculty of Electrical Engineering +) Czech Technical University in Prague, Faculty of Transportation Sciences

[email protected], [email protected], [email protected]

Abstract: Several time-frequency transforms have been introduced that claim better performance in the field of signal classification and compression. One of them is the discrete cosine transform (DCT), which has become a basic tool for speech and image signal processing. The abilities of DCT are limited by the cosine basis functions, which do not achieve adequate results in detecting a non-stationary signal. We present the selective Zolotarev-cosines which are able to extend the operation of DCT by substituting the cosine basis functions.

1. Introduction

Yegor Ivanovich Zolotarev was a Russian mathematician who lived in the second half of the 19th century. He stated and solved four problems in approximation theory. The first and second problems concern polynomials which deviate least from zero in two disjoint intervals

w−− ;1 and 1;−w . Closed form solutions of these approximation problems led to the set of

polynomials which forms fundamentals for the discrete Zolotarev transform (DZT) [1], [2]. Our objective is to substitute basis cosine function of the discrete cosine transform for the selective Zolotarev-cosine function with intention to enhance its properties. Recently we have developed two recursive algorithms for both even, and odd selective cosines. 2. Discrete Cosines Transform

A discrete cosine transform [3] is an orthogonal transform related to DFT. The purpose of DCT is to decorrelate a signal to the cosine series; it therefore gives an N-point real spectrum for an N-point real signal. This property has been employed in many applications in the area of signal processing. 2.1 Definition of DCT There are several variants of DCT with a slightly modified definition. The most common definition of DCT is (1) with coefficient (2), which is derived from DFT of a 2N-point even extension of series x(n).

( ) ( ) ( )∑−

=

+=1

0 2

12cos

N

n

kN

nnxkkC πε (1)

( )

≠

==

02

01

kforN

kforNkε (2)

40 Ján Janík

2.2 DCT basis function As the first we have to take a look at DCT basis function (3). There are generally defined two parameters. The parameter n defines the n-th sample within the interval 1;0 −N while the

parameter k defines degree and symmetry of the function which is composed of k-multiple of one half of standard cosine period. We can define the even and odd cosine function satisfying particular symmetry relations. If k is even, the cosine is even and vice versa, see Fig. 1. This specific definition is necessary for section 3.

+ πkN

n

2

12cos (3)

It is worth noting that the basis cosine functions are actually a solution of the Chebyshev approximation problem - a polynomial approximation of a constant in the interval 1;1− .

3. Selective Zolotarev cosine Functions

The selective cosines zcos(n,k,w) are understood as solutions of the first and second Zolotarev's approximation problem. The zcos function is determined by parameters n, k and wp, ws . The parameters n and k stand for the same as in (3). The parameters wp, ws determine a shape of the central lobe or lobes according to the function symmetry. Therefore it is essential to define the even and odd cosine in section 2. A selective cosine of the N-th order can be expressed as a weighted sum of cosines of the same parity with degrees not exceeding N.

Figure 1: Standard and selective odd and even cosines.

Ján Janík 41

3.1 Even Selective Cosine The even selective cosine is generally a cosine wave with ability to enlarge central lobe. It is a solution of the approximation differential equation (4)

( )( ) ( )2222

222 141 Zwmdw

dZwww p −=

−− , (4)

where Z in the interval 1;1−=w stands for an even zcos (5). The parameter wp, defines

width of the central lobe, see Fig. 1, and is the most important parameter regulating selectivity of the function. In case when pp ww −= , zcos function is equal to a standard cos. We have

developed the recursive evaluation of the coefficients a2µ(wp)

( ) ( )∑=

+=m

pp N

nwawmnz

02 2

2

12cos,2,cos

µµ µπ . (5)

3.2 Odd Selective Cosine The odd selective cosine is more complicated case. Since the centre of symmetry has zero value, we have to work with two neighboring central lobes, see Fig. 1. Therefore odd zcos introduces three interconnected parameters wm, wp, and ws defining a shape of the central lobes. For the DZCT the most attractive parameter regarding regulation of selectivity is ws which is the most similar to parameter wp of the even selective cosine. Relation of the parameters wm, wp and ws to the parameter ws can be seen in Figure 2. These parameters appear in an approximation differential equation (6)

( )( )( ) ( ) ( )( )222222

22222 1121 msp wwZmdw

dZwwwww −−+=

−−− . (6)

In eq. (6) the function Z stands for an odd zcos (7) within the interval 1;1−=w . We have

developed the recursive evaluation of the coefficients a2µ+1(ws)

( ) ( ) ( )∑=

+

++=+m

ss N

nwawmnz

012 12

2

12cos,12,cos

µµ πµ . (7)

42 Ján Janík

Figure 2: Relation between parameters ws, wm and wp.

4. Conclusion Recently we have performed several experiments with our new DZCT and the results are promising. In short future we intend to develop a complete mathematical theory dealing with the selective cosine basis and use it for DZCT applications. We will perform this research within the project "Novel selective transforms for non-stationary signal processing" of the Grant Agency of the Czech Republic in 2011-2013. 4. Acknowledgements

This work was supported by the Grant Agency of the Czech Technical University in Prague, grant No. SGS10/181/OHK3/2T/13 “Spectral properties of Zolotarev transform”, by the Grant Agency of the Czech Republic GD102/08/H008 “Analysis and modeling biomedical and speech signals” and the research program Transdisciplinary Research in Biomedical Engineering II No. MSM6840770012 of the Czech University in Prague. References [1] R. Špetík, The Discrete Zolotarev Transform, Doctoral Thesis. CTU in Prague, Czech

Republic, 2009. [2] M. Vlček, R. Unbehauen. Zolotarev Polynomials and Optimal FIR Filters, IEEE

Transactions on Signal Processing, vol. 47, no. 3, March 1999, pp. 717-730. [3] N. Ahmed, T. Natarajan and K. R. Rao. Discrete Cosine Transform. IEEE

Transactions on Computers, January 1974, pp. 90-93.

Ján Janík 43

Optimalizace algoritmů rozpoznávačů řečis využitím architektur vícejádrových signálových procesorů

Robert Krejčí


Abstrakt: Tento článek popisuje hardwarové a softwarové možnosti vytvoření rozpoznávače řeči založeného na platformě signálových procesorů TMS320C64x+. Je popsáno jedno z možných řešení výpočetně náročných úloh zpracování signálů formou návrhu koncepce DSP serveru, který je schopen obsloužit požadavky na rychlý výpočet úloh číslicového zpracování signálů. Je popsána myšlenka „blízkého“ a „vzdáleného“ DSP serveru s využitím vícejádrových signálových procesorů. Dále jsou popsány některé optimalizační metody pro zrychlení vybraných částí rozpoznávače řeči na zvolené hardwarové platformě: parametrizace a výstupní pravděpodobnostní funkce.

1. ÚvodExistuje mnoho typů úloh zpracování signálů, které jsou výpočetně velmi náročné, např. rozpoznávání nebo syntéza řeči, zpracování obrazů a videa, analýza biologických signálů. Existují aplikace, ve kterých není možné použít „běžnou“ platformu PC např. z důvodu požadavku na přenositelnost, miniaturní rozměry nebo malou spotřebu. V některých aplikacích také může být problematická značná produkce tepla, která může vést ke snížení spolehlivosti systému. V takových případech je vhodné použít jinou platformu, která by lépe vyhovovala požadavkům DSP systému.

Mnoho výzkumných týmů používá pro vývoj algoritmů rozpoznávání řeči běžnou platformu PC. Při vývoji rozpoznávačů řeči pro jiné platformy však často není efektivní pouze převzít algoritmy optimalizované pro platformu PC. Je vhodnější vytvářet spíše nové optimalizační metody a algoritmy, s využitím znalostí o zvolené hardwarové architektuře, které vedou k efektivnějšímu provádění výpočtů a zpracování dat.

Tento článek popisuje jednak hardwarové uspořádání a také softwarové optimalizace použité při vývoji rozpoznávače řeči s využitím vícejádrových signálových procesorů.

2. Koncept DSP serveru

2.1. Trend k vícejádrovým systémůmSoučasná nabídka významných výrobců signálových procesorů (Texas Instruments, Freescale) jednoznačně směřuje k rozvoji vícejádrových architektur, ať už homogenních (více stejných jader na jednom čipu), nebo heterogenních (jádro procesoru pro všeobecné použití se signálovým jádrem na jednom čipu). [TICOM] Tato skutečnost nás vede k vytvoření myšlenky DSP serveru a jeho využití pro rozpoznávání řeči. [KRE10P]1

1 Příspěvek oceněný 2. místem na konferenci POSTER 2010 v kategorii Electronics and Instrumentation.

44 Robert Krejčí

2.2. Architektura signálových procesorů Texas Instruments TMS320C64x+Jako jedna z velmi vhodných architektur, které jsou na trhu běžně dostupné, se ukazuje např. architektura signálových procesorů od firmy Texas Instruments s označením TMS320C64x+ pracující s pevnou řádovou čárkou. [TMS06] Je to jedna z nejvýkonnějších architektur signálových procesorů a lze ji v současné době považovat za průmyslový standard2. Jedná se o 32bitovou architekturu typu VLIW3 s osmi paralelními jednotkami a s instrukční sadou umožňující rychlé provádění operací typu SIMD4.

2.3. Koncept DSP serveruServer je obecně elektronické zařízení, které je schopno obsloužit požadavek na zpracování dat. Běžně používáme HTTP server, který generuje webové stránky, SQL server, který vyhledává v databázích nebo do nich ukládá data, e-mailový server, který slouží k elektronické komunikaci.

Na podobném principu lze vybudovat DSP server, který bude schopen obsloužit požadavky na rychlé zpracování signálových dat. Jak si ukážeme v následujícím textu, může se jednat o softwarový server v rámci jednoho čipu, který jsme nazvali „blízký (near) DSP server“, nebo to také může být zařízení připojené k systému externě pomocí síťového (nebo jiného – sériového, paralelního) rozhraní, to jsme pojmenovali jako „vzdálený (far) DSP server“. Klíčová slova near a far jsme zvolili na základě analogie s klíčovými slovy v programovacím jazyce C, která řídí umístění proměnné v paměti s rychlejším přístupem (near), nebo v obecné paměti (far).

2.4. Blízký DSP serverV současné době je běžné řešit výpočetně náročné úlohy pomocí vícejádrových procesorových systémů. V našem případě používáme pro výzkum a testování algoritmů rozpoznávačů řeči dvoujádrový heterogenní procesor řady OMAP s velmi nízkým příkonem. Procesor lze nakonfigurovat tak, že na jádru ARM běží operační systém Linux a na jádru TMS320C674x funguje realtimový systém označovaný jako DSP/BIOS. Obě jádra mezi sebou komunikují pomocí softwarové vrstvy označované jako DSP/BIOS Link. Komunikace je založena na sdílené paměti.

2.4.1. Rozdělení procesůCelý proces zpracování signálů je řízen procesorem pro všeobecné použití ARM. Ten se stará o komunikaci s uživatelem, obsluhu rozhraní, pravidelný odběr vzorků z A/D převodníku atd. Jádro signálového procesoru je naprogramováno jako DSP server, který je schopen relativně rychle obsloužit požadavek na výpočet časově náročné operace zpracování signálů. Vzhledem k umístění na stejném čipu nazvěme tuto konfiguraci jako blízký DSP server. V případě, že je potřeba provést nějakou výpočetně náročnou operaci, procesor pro všeobecné použití se obrátí s požadavkem a příslušnými daty na signálový procesor. Ten je, díky architektuře přizpůsobené pro provádění operací typických pro zpracování signálů, schopen provést výpočet mnohem rychleji než procesor pro všeobecné použití. Po dokončení výpočtů vrátí výsledky zpět jádru procesoru pro všeobecné použití.

Při zadávání požadavku signálovému jádru však je potřeba počítat s určitou časovou prodlevou pro sestavení komunikace a předání dat (řádově desítky až stovky mikrosekund).

2 Několik týdnů před zveřejněním této publikace byla vydána nová ještě výkonnější architektura TMS320C66x pracující s plovoucí řádovou čárkou.

3 VLIW = Very Long Instruction Word. Velmi dlouhé instrukční slovo.4 SIMD = Single Instruction Multiple Data. Současné zpracování více dat během jedné instrukce.

Robert Krejčí 45

Např. v případě výpočtu komplexní FFT do 64 vzorků se nevyplatí spouštět proces na signálovém jádru – výpočet na řídicím jádru ARM je rychlejší. Při 128 vzorcích trvají obě varianty stejně dlouho a od 256 vzorků je použití signálového jádra jednoznačně rychlejší. [TMS10]

2.4.2. Vývojové nástrojeDonedávna bylo nutné vyvíjet programy pro každé jádro zvlášť a poměrně složitě řešit komunikaci mezi oběma jádry. Ve druhé polovině roku 2010 Texas Instruments uvedl nový vývojový nástroj C6Run, který zajišťuje komunikaci mezi oběma jádry procesoru automaticky. Existují dvě varianty tohoto nástroje. [TMS10]

Aplikace C6RunLib umístí výpočetně náročné funkce do oddělené knihovny. Hlavní program běží na jádru ARM, oddělená knihovna je spouštěna na signálovém jádru. Tento nástroj je vhodný zejména pro linuxové vývojáře, kteří nyní mohou jednoduše využívat výhod jádra signálového procesoru.

Aplikace C6RunApp zajistí komunikaci „z opačného pohledu“. Jádro signálového procesoru pomocí této aplikace získá přístup k souborovému systému, k obrazovce a klávesnici atd., což bylo dosud v podstatě nemyslitelné. Tato aplikace je naopak určena primárně vývojářům na platformě signálových procesorů, kteří nyní mají snadnější přístup k periferiím procesoru.

2.5. Vzdálený DSP server V případě, že je potřeba ještě větší výpočetní výkon, nabízí se možnost začlenit do systému další jednotky založené na vícejádrových homogenních procesorech. Jedná se o jednotky připojené k řídicímu systému externě pomocí vhodného komunikačního rozhraní, proto nazvěme tuto konfiguraci jako vzdálený DSP server. Počet připojených jednotek (modulů) odpovídá požadavkům na výpočetní výkon.

V našem případě pracujeme s 6jádrovými procesory rovněž od výrobce Texas Instruments s označením TMS320C6472. V případě použití součástky s taktovací frekvencí 500 MHz výrobce uvádí spotřebu 3,4 W při plném výpočetním výkonu 24 GMACS5. Nepatrný příkon spolu s efektivní hardwarovou architekturou jsou jedny z hlavních důvodů, proč jsme zvolili toto řešení.

2.6. Rozdělení výpočetní zátěžeRozhodnutí o umístění výpočtu na blízký, nebo vzdálený DSP server je zcela v rukou programátora. Zatím neexistuje žádný nástroj, který by to dělal automaticky. Přitom se lze rozhodovat podle jednoduchého kritéria: 1) Rychlé výpočty, jejichž doba je srovnatelná s prodlevou pro komunikaci mezi oběma jádry, se počítají na řídicím jádru. 2) Středně náročné výpočty se umístí na blízký DSP server, pokud je výpočet včetně komunikace rychlejší než na řídicím jádru a pokud stihne provést v dostatečném časovém intervalu. 3) Velmi náročné výpočty se umístí na vzdálený DSP server; navíc je nutné optimalizovat jejich rovnoměrné rozložení na všechna jádra.

V případě programování rozpoznávače řeči může být rozdělení procesů následující: odběr vzorků z A/D převodníku zajišťuje řídicí jádro procesoru pro všeobecné použití, parametrizaci počítá blízký DSP server a další výpočetně velmi náročné algoritmy, jako je např. výstupní pravděpodobnostní funkce, nebo algoritmus token passing, počítá vzdálený DSP server.

5 GMACS = Giga Multiply and Accumulate Operations per Second. Miliard operací násobení s akumulací za sekundu.

46 Robert Krejčí

Dále je potřeba řešit rovnoměrné rozdělení úkolů mezi jednotlivá jádra signálových procesorů. Např. výstupní pravděpodobnostní funkce, v literatuře běžně označovaná jako funkce b(o), která vyhodnocuje akustickou podobnost mezi vstupním signálem a natrénovanými daty, se skládá z několika set (v našem případě 435) nezávislých smyček. Přitom výpočet každé smyčky lze umístit do jednoho jádra, takže v jednom okamžiku se počítá současně šest smyček, resp. jejich násobky v případě použití více šestijádrových procesorů.

2.7. Porovnání příkonuVýsledky našich testů algoritmů rozpoznávání řeči ukazují, že pro dosažení výpočetního výkonu moderního PC jsou potřeba tři šestijádrové procesory TMS320C6472 s taktovací frekvencí 500 MHz. Rozdíl v příkonu však je značný. Zatímco PC při plném výpočetním zatížení odebírá příkon řádově 100 W, v případě použití procesorů s architekturou TMS320C64x+ odhadujeme příkon na 10 až 15 W.

Ukazuje se tedy, že tato navržená konfigurace je vhodná pro výkonné nízkopříkonové přenosné systémy pro zpracování signálů a rozpoznávání řeči v reálném čase.

3. Optimalizace algoritmů rozpoznávačů řeči pro platformu signálových procesorů TMS320C64x+

Představili jsme možné uspořádání hardwarového systému pro rozpoznávání řeči založeného na vícejádrových signálových procesorech řady TMS320C64x+. Nyní popíšeme některé možnosti optimalizace algoritmů rozpoznávání řeči pro tuto platformu.

Nejprve je potřeba zmínit některé specifické vlastnosti této hardwarové architektury, které je potřeba při optimalizaci algoritmů mít na paměti. [TMS06]

➢ Maximální propustnost datových sběrnic je 128 bitů pro čtení nebo zápis do interní (on-chip) paměti v každém instrukčním cyklu.

➢ Instrukční sada umožňuje ve velké míře využívat operace SIMD pro datové typy 16bitový signed a 8bitový unsigned. Podpora operací SIMD je poněkud více omezená pro datové typy 16bitový unsigned a 8bitový signed.

➢ Lze zpracovávat paralelně až 8 operací, avšak např. pouze 2 operace typu násobení s akumulací současně.

➢ Na rozdíl od platformy PC zde není tak markantní rozdíl mezi taktovací frekvencí jádra procesoru a externí paměti. Tato zdánlivě nenápadná vlastnost však vede k vytváření zcela odlišných algoritmů (viz dále).

➢ Jedna podskupina této architektury (označovaná jako TMS320C674x) také poskytuje možnost výpočtů ve formátu operandů s plovoucí řádovou čárkou. Možnosti optimalizací jsou však poněkud omezené oproti výpočtům v pevné řádové čárce. Nevýhoda instrukční sady pro výpočty v plovoucí řádové čárce na této architektuře spočívá v tom, že neobsahuje instrukce s využitím principů SIMD. Výsledky výpočtů jsou sice velmi přesné, ale pomalejší oproti výpočtům v pevné řádové čárce s plným využitím možností hardwarové architektury.

Robert Krejčí 47

3.1. Optimalizace výpočtu parametrizace MFCCJedním z prvních bloků rozpoznávače řeči je parametrizace. Cílem parametrizace je extrahovat ze vstupního řečového signálu takové parametry, které dostatečným způsobem vypovídají o jeho řečových vlastnostech a zároveň potlačují redundantní informace jako např. základní tón řeči, emoce mluvčího a další. Existuje více algoritmů parametrizace, v našem případě používáme parametrizaci pomocí mel-kepstrálních koeficientů. [UHL07]

Součástí parametrizace pomocí mel-kepstrálních koeficientů je diskrétní kosinová transformace:

c i= 2P ∑j=1

P

g jcosiP

j−0,5 ,1iM−1,

kde ci je i-tý koeficient parametrizace, gj je logaritmus energie v j-tém spektrálním pásmu, P je počet pásem (počet trojúhelníkových filtrů) a M je požadovaný počet koeficientů MFCC, přičemž obvykle bývá M < P.

Základem této funkce je výpočet kosinu. Pokud se na tento vztah podíváme pozorněji, vidíme, že argument kosinu vůbec nezávisí na neznámé proměnné gp. V argumentu se vyskytují pouze indexy, které souvisejí se základní konfigurací, která je známá během kompilace: počet pásem melovských filtrů P a počet koeficientů diskrétní kosinové transformace M. Tyto kosiny v algoritmu figurují jako konstanty, a můžeme je tedy spočítat předem. Nový algoritmus se redukuje pouze na operaci typu skalární součin:

c i=∑j=1

P

g j⋅costab[ i , j ]=DOTPROD g ,costab[ i ] ,

kde costab je dvourozměrná matice s předem spočítanými koeficienty. V našem případě má tato datová struktura 32×16 prvků.

Násobení s akumulací se, na rozdíl od čistě matematického přístupu, dá pokládat za elementární operaci. Signálové procesory obvykle mohou takovou operaci provést během jednoho instrukčního cyklu. Některé hardwarové architektury navíc umožňují počítat v každém instrukčním cyklu operace typu DOT PRODUCT, což znamená paralelní násobení vektoru operandů a součet výsledků do jedné sumy.

Výsledky vlivu vybraných optimalizačních metod na dobu výpočtu parametrizace MFCC jsou zobrazeny v grafu na Obrázku 1. Zatímco bez optimalizace trval výpočet každého segmentu signálu přes 55 ms, po aplikaci vhodných optimalizačních metod došlo přibližně k desetinásobnému zrychlení výpočtu na 5,5 ms, přičemž největší přínos měla právě popsaná metoda vyhledávací tabulky.

Obrázek 1: Vliv vybraných optimalizačních metod na dobu výpočtu parametrizace MFCC

48 Robert Krejčí

Tato optimalizační metoda ukazuje typický rozdíl oproti platformě PC. Ukazuje se, že na platformě PC je rychlejší přímo spočítat kosinus než čekat na relativně pomalé přečtení tabulky dat z paměti. V případě naší hardwarové architektury tomu je přesně naopak.

Další dvě metody – optimalizace výpočtů energie v pásmech pomocí vhodně připravených datových struktur a začlenění kruhového FIFO bufferu – měly už menší přínos, ale ukazuje se, že rovněž vedou ke zrychlení výpočtů.

3.2. Modifikace výstupní pravděpodobnostní funkceSoučástí akustického modelu rozpoznávače řeči je výpočet hustoty výstupní pravděpodobnosti. Tato funkce znamená ohodnocení zvukové shody mezi vstupním segmentem a modelovanou referencí [UHL07] a počítá se v každém segmentu pro každý emitující stav skrytých Markovových modelů:

b jo t=∏s=1

S [∑m=1

M s

c jsm N ost ; jsm , jsm]s

;

N o ; ,=1

2n∣∣e−1

2 o− ' −1o− ,

kde S je počet proudů (skupin parametrů), γs je váha proudu, Ms je počet směsí v proudu, cjsm

je váha m-té složky, N(·; μ, Σ) je vícerozměrné Gaussovo rozdělení s vektorem středních hodnot μ a kovarianční maticí Σ. [YOU09]

Jedná se o výpočetně velmi náročnou funkci. Důvodem výpočetní náročnosti není ani tak extrémní složitost algoritmu, ale spíše obrovské množství dat, se kterými se opakovaně pracuje.

Pokud na tuto funkci aplikujeme některé optimalizační metody, získáme tento tvar, který pro naše experimenty budeme považovat za referenční:

ln b jo t = ln ∑m=1

M s≈40

expX jm∑k=1

K≈51

[otk− jmk ⋅y jmk ]2

3.2.1. Optimaizovaná funkce – TVAR 1 – MACV publikaci [KRE09S] jsme představili první modifikaci této funkce vhodnou pro architektury signálových procesorů. Pokud se s původní funkcí b(o) provede následující substituce, v jádru funkce bude eliminováno odečítání a druhá mocnina:

pk = ok2 , z jmk = yk

2 , v jmk =−2 jmk y jmk2 , X jm1 = X jm∑

k=1

K≈51

y jmk2 ⋅ jmk

2

ln b jo t = ln∑i=1

M

exp X jm1∑k=1

K≈51

z jmk⋅p jmkv jmk⋅o jmkZůstaly nám pouze operace typu „násobení s akumulací“ (Multiply and Accumulate – MAC), což pokládáme za elementární operaci, jak již bylo řečeno. V jádře funkce jsou dvě elementární operace typu MAC, které jsou navíc proveditelné paralelně.

3.2.2. Optimaizovaná funkce – TVAR 2 – ADD-MACNašli jsme ještě jiný tvar funkce b(o), který vrací numericky stejné výsledky, obsahuje stejný počet elementárních operací jako původní tvar, ale je rychleji proveditelný na některých

Robert Krejčí 49

hardwarových architekturách. [KRE10CG] Pokud se provede následující substituce, dostaneme další tvar modifikované funkce b(o):

z k , m, i = yk , m ,i2 , t k , m,i =−2k ,m ,i , U k , m = X k , m∑

i=1

51

yk ,m ,i2 ⋅k , m, i

2

ln bk o = ln ∑m=1

M ≈40

expU k , m∑i=1

K≈51

z k , m ,i⋅o i⋅oit k , m, iTato funkce obsahuje násobení a součet. Tyto dvě operace mohou být provedeny paralelně, pokud to příslušná hardwarová architektura umožňuje. Následuje násobení a akumulace. Další průběh je stejný jako v předchozích případech.

Bitový rozsah nově přepočítaných koeficientů je menší než u prvního tvaru modifikované funkce. Proto také můžeme dosáhnout přesnějších výsledků v případě výpočtů se sníženou bitovou přesností.

3.2.3. Zpřesnění výpočtů pásmovým škálovánímV obou modifikacích funkce b(o) lze dosáhnout větší přesnosti, pokud veškerá data rozdělíme do více pásem, která nemusí nutně odpovídat uspořádání statických, diferenciálních a akceleračních parametrů [KRE09S], [KRE10CG]. V tomto případě jsme podle dynamického rozsahu zvolili tři pásma s 16, 16 a19 koeficienty. Každé pásmo se nyní lépe vejde do rozsahu 16bitových čísel. Počet trénovacích parametrů, které v důsledku podtečení budou vynulovány, je podstatně menší. Výpočet probíhá v každém pásmu s optimálním bitovým rozsahem. Výsledky jednotlivých smyček v jádře funkce jsou srovnány na stejnou bitovou úroveň a další výpočet, který už není tak časově kritický, probíhá stejně jako v referenční funkci.

3.2.4. Porovnání obou modifikovaných funkcíVýsledky měření přesnosti a rychlosti na signálových procesorech s architekturou TMS320C64x+ ukazují, že v obou případech lze dosáhnout podstatně menší chyby, když výpočet probíhá v jednotlivých pásmech s optimální bitovou přesností, přičemž se ukázalo, že výpočetní čas spotřebovaný na uspořádání jednotlivých smyček je minimální. Průměrováním výsledků ze 100 realizací jsme zjistili, že algoritmy s globální redukcí přesnosti vedou k odchylce kolem 3%. Naproti tomu v případě, že výpočet probíhá ve třech pásmech, přičemž každé pásmo má svou optimální bitovou přesnost, lze dosáhnout poměrně přesného výsledku vzhledem k původní funkci b(o), v případě nové funkce i pod 0.5%. Dále se ukazuje, že modifikovaná funkce s operací MAC (TVAR 1), je stále nejrychlejší. Přitom funkce s operacemi ADD-MAC (TVAR 2) je rovněž rychlejší než referenční funkce.

Obrázek 2: Blokové schéma modifikované funkce b(o) - TVAR 2 - ADD-MAC

50 Robert Krejčí

4. ZávěryTento článek prezentuje některé výsledky našeho výzkumu zaměřeného na vývoj miniaturních, nízkopříkonových či přenosných systémů rozpoznávání řeči pracujících v reálném čase. K tomu využíváme kombinaci moderních vícejádrových signálových procesorů navrženou pro řešení výpočetně náročných úloh číslicového zpracování signálů se zaměřením na algoritmy rozpoznávání řeči, kterou lze zobecnit pod pojmem DSP server. Pro urychlení algoritmů rozpoznávače řeči využíváme modifikované algoritmy, které efektivně využívají možností zvolené hardwarové architektury. Přestože předmětem našeho dosavadního výzkumu je především optimalizace algoritmů rozpoznávačů řeči pro skupinu hardwarových architektur rodiny TMS320C64x+, ukazuje se, že popsané optimalizace jsou v mnoha případech obecnější, a tedy aplikovatelné i na jiné podobné architektury.

5. PoděkováníTento výzkum byl podporován z grantu GAČR 102/08/H008 “Analýza a modelování biomedicínských a řečových signálu”, GAČR 102/08/0707 “Rozpoznávání mluvené řeči v reálných podmínkách“ resp. výzkumného záměru MŠMT MSM6840770014 “Výzkum perspektivních informačních a komunikačních technologií”.

Reference[KRE09S] KREJČÍ, Robert: Optimalizace výpočetně náročné části rozpoznávače řeči se zaměřením na hardwarovou platformu OMAP. In Analýza a zpracování řečových a biologických signálů – sborník prací 2009. Praha: České vysoké učení technické v Praze, 2009, s. 50-57. ISBN 978-80-01-04474-2.

[KRE10P] KREJČÍ, Robert: Building DSP Server on TMS320C64x+ Platform. In POSTER 2010 - Proceedings of the 14th International Conference on Electrical Engineering. Praha: ČVUT v Praze, FEL, 2010. ISBN 978-80-01-04544-2.

[KRE10CG] KREJČÍ, Robert: Methods for Faster Calculations of Output Probability Density. In 20th Czech-German Workshop on Speech Processing [CD-ROM]. Praha: Institute of Photonics and Electronics AS CR, 2010.

[TICOM] Texas Instruments [online]. 2010 [cit. 2010-09-07]. Digital Signal Processors & ARM Microprocessors. URL: <http://focus.ti.com/dsp/docs/dsphome.tsp?sectionId=46&DCMP=TIHomeTracking&HQS=Other+OT+home_p_dsp>.

[TMS06] Texas Instruments: TMS320C6000 Programmer's Guide [online]. 2006 [cit. 2010-09-07]. URL: <http://www.ti.com/litv/pdf/spru198i>.

[TMS10] Allred Daniel: C6Run DSP Software Development Tool [online]. 2010 [cit. 2010-11-24]. URL: <http://www.ti.com/litv/pdf/spry143>.

[UHL07] UHLÍŘ, Jan, et al. Technologie hlasových komunikací. Praha: Nakladatelství ČVUT, 2007. 276 s. ISBN 978-80-01-03888-8.

[YOU09] Young Steve. et al., The HTK Book. cambridge University Engineering Department, 2006. [online] URL: <http://htk.eng.cam.ac.uk/ftp/software/htkbook.pdf.zip>.

Robert Krejčí 51

Zpětnovazební regulace mikromechanickýchexperimentů

Ondřej KučeraČeské vysoké učení technické v Praze, Fakulta elektrotechnická

Ústav fotoniky a elektroniky, Akademie věd Č[email protected]

Abstrakt: Scanning Probe Microscopy, a large family of microscopies, is asensitive technique with unprecedented resolution. However, the accuracy ofthe measurement is strongly influenced by the setting and calibration of theapparatus. In this article, we present system analysis approach to the problemof the feedback adjustment and probe’s stiffness selection.

1. Úvod

Princip mikroskopie skenující sondou [5, 1, 3, 4, 6], spočívá v silové interakci mezi velmimalou sondou (typicky s poloměrem křivosti v řádu desítek, nebo jednotek nanometrů) avzorkem. Typy interakcí jsou zejména:

• atomové síly, jako síly van der Waalsovy, Pauliho odpuzování, atd.,

• tunelující proud,

• elektrostatické síly,

• optické síla,

• elektromagnetický odraz,

• chemické síly,

• a mnoho dalších.

Vzhledem k tomu, že mikroskopie skenující sondou je schopna měřit pouze na místníúrovni, topografie se rekonstruuje skenováním povrchu vzorku. Mikroskopie skenující son-dou je také široce používaná jako nezobrazovací technika umožňující měření lokální vlast-ností vzorku, jako je Youngův modul, permitivita, hustota náboje, přilnavost, atd. Proudržení požadované přesnosti měření se používá zpětná vazba. V tomto článku uvádímeanalýzu problému nastavení zpětné vazby a přizpůsobení tuhosti sondy.

52 Ondřej Kučera

2. Lineární model

Vzhledem k tomu, že se lineární model založený na tlumeném oscilátoru [2] (viz Obr.1) mnohdy ukazuje jako dostatečně přesný, můžeme reprezentovat systém sonda-vzoreknásledující rovnicí pohybu

md2z

dt2+ β

dz

dt+ (k + κ)z = κzo, (1)

kde m je efektivní hmotnost, k je tuhost, a β je tlumení sondy, zo je topografie nezatíženéhovzorku, κ je tuhost povrchu vzorku a z je vzdálenost sondy od povrchu. Chceme-li býtpřesní, musíme konstatovat, že κ obsahuje i tuhost interakce, ale to není v lineárnímrežimu důležité.

F

κ

z

zo

m

k β

Sample

Cantilever

Obrázek 1: Model systému (vpravo) a odpovídající prvky Atomic Force Microscopy(vlevo).

Zpětná vazba (zpravidla proporcionálně-integrální) se používá k udržení žádaného pra-covního bodu a zvýšení citlivosti v případě vysokého poměru k/κ. Můžeme ji do rovnice1 zavést následujícím způsobem

md2z

dt2+ β

dz

dt+ (k + κ)z = κ

zo − Pz − I

t∫

0

z(τ)dτ

︸︷︷︸zp.vazba

, (2)

kde P a I jsou proporcionální a integrální zisk zpětné vazby. Přenos zpětnovazebně řízenéhosystému sonda-vzorek můžeme psát jako

H(s) =sκ

s3m + s2β + s(k + κ + Pκ) + Iκ, (3)

kde s je komplexní Laplaceův parametr. V ideálním případě L−1H → 0 a informace otopografii je obsažená v signálu zpětné vazby s přenosem

Hh(s) =κ(sP + I)

s3m + s2β + s(k + κ + Pκ) + Iκ, (4)

a odpovídající amplitudovou frekvenční charakteristikou

|Hh(jω)| =√√√√ κ2 (P 2ω2 + I2)

(Iκ− ω2β)2 + (ω (k + κ + Pκ)− ω3m)2 , (5)

kde ω je úhlová frekvence.

Ondřej Kučera 53

3. Nelineární model

0

z

F(z

)

Obrázek 2: Lennard-Jonesova síla v závislosti na vzdálenosti z.

Síla působící mezi sondou a vzorkem, nemůže být v některých případech linearizovanáa nelineární povaha interakce κ = κ(z) musí být brána v úvahu. Vzhledem k tomu,že interakce sleduje obvykle derivaci Lennard-Jonesova potenciálu (viz Obr. 2), tuhostinterakce je pak

κ(z) = κo∥4ε

z∇(

ζ

z

)12

−(

ζ

z

)6 , (6)

kde ε je konstanta (hloubka potenciálové jámy), ζ představuje konečnou (!) vzdálenostnulového potenciálu a κo je tuhost samotného vzorku. Existují samozřejmě i další ne-linearity v systému, ale lze je mnohdy kompenzovat, či zanedbat. Výsledné rovnice jetřeba řešit numericky. Pro nejpřesnější analýzu lze použít model ve stavovém prostoru,kteý navrhl Stark et al. [7]. Tento model používá pohybovou rovnici vetknutého nosníkunamísto tlumeného oscilátoru.

4. Výsledky a závěr

Správně nastavená zpětná vazba zvyšuje citlivost mikroskopie skenující sondou. Přestomůže být citlivost naopak významně zhoršená špatným nastavením zpětné vazby, a tozejména pro vysoký poměr κ/k. Hodnota |Hh(jω)| (rov. 5) vyjádřená v procentech před-stavuje vertikální přesnost měření. Cílem nastavení zpětné vazby je upravit P a I ziskyzpůsobem umožňujícím získat nejvyšší citlivost, a zároveň neuvádějící systém do oscilacízpůsobených časovým zpožděním z-piezo pohonu. Nesprávné nastavení vede k potlačenívertikálních informací obrazu a nakonec ke ztrátě detailů. Sonda se pak více či méněboří do vzorku a může mechanicky narušit nebo vyčerpat zkoumaný proces. To má velkývýznam především pro biologické experimenty.

Poděkování

Práce byla částečně podpořena z grantů č. 102/08/H081 GA ČR a SGS10/179/OHK3/2T/13GA ČVUT.

54 Ondřej Kučera

0.01 0.1 1 10 100

0

0.25

0.5

0.75

1

Relative sample stiffness κ/k [−]

Tra

nsfe

r H [−

]

Obrázek 3: Efekt nastavení zpětné vazby na přenos v propustném pásmu v závislosti napoměru κ/k. Modrá křivka reprezentuje případ s parametry zpětné vazby P = 6, I = 4.Šedé křivky znázorňují hodnoty parametrů zpětné vazby logaritmicky distribuované mezi0,1 a 100. Červená křivka značí případ bez použití zpětné vazby.

100

102

104

−100

−50

0

50

Frequency f [kHz]

Tra

nsfe

r H [d

B]

Obrázek 4: Vliv nastavení zpětné vazby na amplitudovou frekvenční charakteristiku.Parametry sondy jsou m = 3.909−12 kg, k = 0.05 N/m, and β = 1.418−8 Ns/m. Plnáa čárkovaná křivka ukazuje případ se zpětnou vazbou (P = 6,I = 4), respektive bez ní.Solid and dashed curves represent cases with and without feedback (P = 6,I = 4), respec-tively. Jsou ukázány tři odlišné hodnoty poměru tuhostí sondy a vzorku: k = κ = 0.05N/m zelená, k > κ = 0.005 N/m modrá a k < κ = 5000 N/m červená barva.

Reference

[1] Giessibl, F. Advances in atomic force microscopy. Reviews of modern physics 75, 3(2003), 949–983.

[2] Gittes, F.; Schmidt, C. Thermal noise limitations on micromechanical experiments.European Biophysics Journal 27, 1 (1998), 75–81.

[3] Grier, D. A revolution in optical manipulation. Nature 424, 6950 (2003), 810–816.

[4] Parot, P.; Dufrêne, Y.; Hinterdorfer, P.; Le Grimellec, C.; Navajas, D.; Pellequer, J.;Scheuring, S. Past, present and future of atomic force microscopy in life sciences andmedicine. Journal of Molecular Recognition 20, 6 (2007), 418–431.

Ondřej Kučera 55

[5] Poggi, M.; Gadsby, E.; Bottomley, L.; King, W.; Oroudjev, E.; Hansma, H. Scanningprobe microscopy. Anal. Chem 76, 12 (2004), 3429–3444.

[6] Salapaka, S.; Salapaka, M. Scanning probe microscopy. IEEE Control Systems Mag-azine 28, 2 (2008), 65–83.

[7] Stark, R.; Drobek, T.; Heckl, W. Thermomechanical noise of a free v-shaped cantileverfor atomic-force microscopy. Ultramicroscopy 86, 1-2 (2001), 207–216.

56 Ondřej Kučera

Hodnocení koktavosti a experimenty s adaptivním prahem u bayesovského spektrálního detektoru

Tomáš Lustyk


22. listopadu 2010

Abstrakt: V příspěvku je uvedeno stručné shrnutí výsledků výzkumu automa-tického hodnocení neplynulosti řeči. K hodnocení jsou užity parametry jak v ča-sové tak frekvenční oblasti. Výsledky hodnocení jsou porovnávány s hodnocením lékařů. Parametr v časové oblasti průměrná délka ticha dosahuje korelačního ko-eficientu s hodnocením lékařů 0,793, parametr ve frekvenční oblasti počet maxim bayesovského detektoru -0,782. Také je představen experiment s adaptivním prahem pro získání parametru počet maxim bayesovského detektoru, nejlepší vý-sledek tohoto parametru je korelační koeficient -0,726.

1. Úvod

Koktavost je jednou z poruch plynulosti řeči, má mnoho podob a mnoho různých příčin. Koktavost (balbuties) se projevuje opakováním hlásek či slabik (repetice), prodlužováním hlásek (prolongace), častými pauzami a přerušeními v řeči. Oproti jiným poruchám řeči (např. brebtavosti) si koktavý svůj řečový handicap uvědomují. Stres spojený s tím, že si je mluvčí vědom své poruchy, má velký vliv na osobnost člověka a může vést až ke strachu před řečí (logofobii).Správné posouzení tíže poruchy a na ni navazující léčebný postup jsou obtížné úkoly. Určení tíže poruchy je založeno na subjektivním posouzení dle promluvy a chování pacienta, které dnes provádí specialisté na logopedii. Metoda, která by umožnila automaticky a objektivně posoudit tíži poruchy by byla přínosem. Její možné využití by bylo např. při: 1) Určení tíže poruchy. 2) Hodnocení výsledků léčby. 3) Porovnání různých léčebných přístupů.Tento příspěvek krátce popisuje současný stav problematiky hodnocení neplynulosti řeči po-mocí analýzy audio nahrávek pacientů. Je zde také představen experiment modifikující postup pro získání parametru počet změn bayesovského detektoru pomocí adaptivního prahu.

2. Databáze promluv

Databáze signálů vznikala v průběhu posledních několika let na Foniatrické klinice 1. Lékař-ské fakulty Univerzity Karlovy a VFN v Praze. Obsahuje nahrávky od mluvčích různého věku, s různou vážností poruchy plynulosti řeči. Promluvy jsou čtené i volně formulované, nahrávané se zpožděnou akustickou vazbou i bez zpětné vazby. Přibližně 30 signálů jsou kontrolní signály zdravých mluvčích. Zpožděná akustická vazba (DAF - Delayed Auditory Feedback) se pohybuje v rozmezí od 10 do 110 ms. Vzorkovací frekvence při nahrávání byla 44 kHz, pro další experimenty byly signály převzorkovány na 16 kHz.

Tomáš Lustyk 57

Čtené promluvy vznikaly čtením úryvku Podzim na Starém Bělidle z knihy Boženy Němcové Babička, text obsahuje přibližně 70 slov. Volně formulované (spontánní) promluvy vznikaly popisem situace, která byla nakreslena na obrázku.Velmi důležité je, že v průběhu roku 2008 byly všechny promluvy ohodnoceny dvěma foniat-ry, kteří tíži poruchy popsali pomocí 5-ti stupňové klasifikace na základě relativní četnosti ne-plynulých slov, od 0 - žádné projevy koktavosti až 4 - velmi těžká koktavost. Pro každého mluvčího jsou tak k dispozici dvě známky hodnotící neplynulost promluvy. Tato data slouží jako kontrolní pro všechny navrhované algoritmy.

3. Současný stav řešené problematiky

V práci [1] je popsáno velké množství parametrů v časové a frekvenční oblasti umožňujících odlišení plynulé a neplynulé řeči. Parametry jsou navrženy tak, aby zohledňovali různé projevy koktavosti. V následujících dvou podkapitolách budou uvedeny pouze dva parametry reprezentující obě oblasti zpracování. Všechny experimenty jsou prováděny na čtených promluvách bez DAF.

3.1 Parametr v časové oblasti

V úvodu bylo uvedeno, že promluvy balbutiků obsahují velké množství ticha. Lze předpoklá-dat, že parametr zkoumající právě dobu trvání ticha v promluvě by mohl být užitečný k odli-šení neplynulé a plynulé řeči.Parametr průměrná délka ticha vychází právě z tohoto předpokladu. V promluvě jsou pomocí VAD (Voice Activity Detector) nalezeny úseky řeči a ticha. Poté je měřena průměrná délka úseků ticha v celé promluvě. Výsledky parametru měřícího pouze délku ticha nebyly přesvěd-čivé.Neplynulá řeč může obsahovat mnoho repetic a ty jsou také obklopeny tichem. Pokud se tyto krátké úseky řeči odstraní, rozdíl mezi plynulou a neplynulou řečí se zvýrazní. Postup odstra-ňování úseků řeči kratších než určitý práh je demonstrován na obrázku 1.Výsledky tohoto parametru jsou spolu s výsledky dalších parametrů uvedeny v kapitole 3.4. Dalšími mírami v časové oblasti jsou poměr ticha a řeči, počet úseků řeči a ticha a parametry zkoumající energii signálu.

Obrázek 1: Postup zpracování signálu při určení průměrné délky ticha. Nahoře: Průběh signá-lu řeč/ticho. Dole: Průběh signálu řeč/ticho po odstranění úseků kratších než 150 ms.

58 Tomáš Lustyk

3.2 Parametr ve frekvenční oblasti

Parametr zkoumající neplynulou řeč ve frekvenční oblasti, se snaží zohlednit přítomnost pro-longací. Prolongace je nepřirozené prodloužení hlásky, během jejího trvání nedochází ke změně nastavení mluvidel a tedy ani ke změnám ve spektru signálu. Dlouhé úseky ticha se chovají velmi podobně. Počet spektrálních změn by tedy měl ukazovat na tíži poruchy.Příklad výstupu bayesovského detektoru spektrálních změn ukazuje obrázek 2. Výstup ob-sahuje významné i méně významné změny. Pro jejich odlišení je třeba zavést práh. Práh je od-vozen jako zlomek vybraného maxima ve výstupu detektoru (např. pátého nejvyššího). Para-metr počet změn bayesovského detektoru je určen počtem maxim nad prahem.Zkoumán nemusí být pouze počet změn ale i rozestupy maxim ve výstupu detektoru. Tuto ob-last zkoumají další parametry ve frekvenční oblasti jako např. směrodatná odchylka z interva-lů mezi změnami. Dalším přístupem může být také použití jiného detektoru spektrálních změn, např. detektor založený na GLR vzdálenosti, nebo použití rozpoznávače řeči HTK. Vý-sledky všech parametrů ve frekvenční oblasti jsou uvedeny v následující kapitole.

Obrázek 2: Ukázka výstupu bayesovského detektoru a pevného prahu. Nahoře: Výstup detek-toru pro celý signál. Dole: Část výstupu se zobrazeným prahem, křížky jsou označeny lokální maxima, kroužkem maxima nad prahem.

3.4 Kombinace více měr a dosažené výsledky pro čtené promluvy

Úspěšnost všech parametrů byla testována na 120 signálech čtených promluv s dobrou tech-nickou kvalitou nahrávky. Výsledky byly srovnávány s hodnocením lékařů. K porovnání sloužily tři hlediska: korelační koeficient, Wilcoxonův test a odchylka od hodnocení lékařů. Dosažené hodnoty korelačních koeficientů jednotlivých parametrů s hodnocením lékařů jsou ve velké většině případů vyšší jak 0,75. Konkrétně výše uvedené parametry mají korelační koeficient 0,793 pro průměrnou délku ticha a -0,782 pro počet změn bayesovského detektoru.Hodnoty parametrů jsou reálná čísla. Chceme-li provést přímé srovnání hodnocení lékařů s výsledky parametrů je nutné tyto reálné hodnoty diskretizovat do stejné škály jako používají lékaři. To umožňuje diskriminační analýza, jež byla použita a umožnila stanovit odchylku vý-sledků od hodnocení lékařů. V tabulce 1 jsou zobrazeny průměrné odchylky pro tři nejú-spěšnější parametry.Parametry lze různým způsobem kombinovat. Pokud je použito kombinace více parametrů, které zohledňují různé projevy koktavosti, je takto vytvořená míra úspěšnější než jednotlivé parametry samostatně. Průměrná odchylka tohoto systému od hodnocení lékařů je uvedena v tabulce 1. Systém kombinující více různých měr správně odhaduje stupeň neplynulosti pro 63% jedinců, přičemž chyba odhadu o více než jednu třídu se vyskytuje pouze u 1,7% mluvčích.

0 1 0 2 0 3 0 4 0 5 00

1

2

3

BD

1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 90

1

2

3

c a s [ s ]

BD

Tomáš Lustyk 59

průměrná odchylka od hodnocení lékařůkombinace 11 parametrů v čas. i frek. oblasti 0,297směrodatná odchylka vzdálenosti maxim - BACD 0,446pravidelnost energie 0,454průměrný počet maxim - GLR 0,454

Tabulka 1: Průměrná odchylka od hodnocení lékařů.

4. Modifikace určení počtu spektrálních změn pomocí adaptivního prahu

4.1 Algoritmus nastavení adaptivního prahu

Motivací pro užití adaptivního prahu pro určení počtu spektrálních změn může být: 1) Maxi-ma ve výstupu detektoru nemají stálou velikost. 2) Z experimentů uvedených v předchozí kapitole bylo vyřazeno 33 signálů, které měly špatnou technickou kvalitu, způsobenou převážně změnami hlasitosti v průběhu nahrávání.Možností určování prahu adaptivně je několik, v práci [3] byly popsány dvě možnosti, jejich nejlepší výsledky však výrazně zaostávali za pevným prahem použitým v práci [1]. Nejvyšší korelační koeficient 0,709, oproti -0,782.Další způsob určení adaptivního prahu vychází z postupu určení prahu pevného. Výstup de-tektoru je zpracováván po oknech. V okně je vybráno maximum (např. 5-té nejvyšší) z nějž je práh odvozen. Každé okno tak má jiný práh pro oddělení významných maxim od méně vý-znamných, více obrázek 3.

Obrázek 3: Výstup Bayesova detektoru, kroužkem označené maximum je maximum, sloužící k odvození hodnoty prahu pro dané okno.

Takto neošetřený algoritmus aplikovaný na výstup bayesovského detektoru způsobí, že práh kolísá a výrazně ovlivní hodnotu parametru počet změn bayesovského detektoru (počet maxim výrazně vzroste). V místech, kde jsou maxima, která nepředstavují výrazné spektrální změny, práh poklesne a určí jako výrazné spektrální změny i ty, které ve skutečnosti výrazné nejsou Vše ilustruje obrázek 4, zakroužkované části.Obrázek také ukazuje hodnotu směrodatné odchylky výstupu Bayesova detektoru a maxim v jednotlivých oknech. Můžeme si všimnout, že v místech, kde dochází ke zmíněnému zkres-lení výsledků, směrodatná odchylka kolísá. Právě tohoto kolísání je využito pro regulaci hodnoty adaptivního prahu. Po použití směrodatné odchylky k úpravě hodnoty prahu, vidíme, že hodnota prahu „vystřelí“ nahoru pokud je směrodatná odchylka velmi malá, naopak není-li, zůstává v normálních mezích.

0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 20

0 . 2

0 . 4

0 . 6

0 . 8

1

60 Tomáš Lustyk

Obrázek 4: Shora dolů: Výstup Bayesova detektoru a práh s regulací pomocí směrodatné od-chylky (červeně jsou označena všechna maxima, zeleně maxima nad prahem, černě je vyzna-čen práh), výstup Bayesova detektoru a práh bez regulace, směrodatná odchylka, řečový signál.

4.2 Výsledky experimentů s adaptivním prahem

Oba způsoby byly vyzkoušeny na 120 čtených promluvách s dobrou technickou kvalitou na-hrávky, které byly součástí výše uvedených experimentů. První způsob (bez regulace) neměl dobré výsledky. Proto bylo od tohoto způsobu upuštěno.Práh s regulací pomocí směrodatné odchylky má výrazně lepší výsledky. Na obrázku 5 je po-mocí funkce Matlabu boxplot zobrazen výsledek nejlepšího nastavení. V záhlaví jsou uvedeny korelační koeficienty a výsledky Wilcoxonova testu pro dvě hladiny významnosti.Wilcoxonův test testuje hypotézu, že mediány mezi jednotlivými skupinami jsou stejné. Např. zápis „0 0 1 0“ značí, že hypotézu o shodných mediánech zamítáme mezi daty pro známky 2 a 3. Jednoduše řečeno, zkoumaný parametr by mohl být užitečný pro rozeznávání, zda mluvčí patří do skupiny 2 nebo 3. Pro ideální parametr třídící do pěti skupin bychom měli dostat osm jedniček. Vyzkoušeno bylo zatím jen velmi málo nastavení. Nejlepších výsledků dosáhlo nastavení s délkou okna 64000 vzorků (4 s). Hodnota prahu je určena jako 0,25 násobek čtvrtého nej-vyššího maxima. Dosažený korelační koeficient -0,726 a sedm jedniček u Wilcoxonova testu u prvního lékaře. Tento způsob určení prahu nedosahuje kvalit pevného prahu (korelační ko-eficient -0,782), ale jeho výsledky jsou lepší než u způsobů uvedených v [3].

5. Závěr

Práce shrnuje dosavadní výsledky při řešení problematiky koktavosti a neplynulosti na ka-tedře teorie obvodů ČVUT FEL ve spolupráci S Foniatrickou klinikou 1. LF UK a VFN v Praze a modifikaci spektrálního detektoru při použití adaptivního prahu.Všechny experimenty analyzující neplynulou řeč jsou prováděny na čtených promluvách bez zpožděné akustické vazby. Pro analýzu jsou použity parametry v časové i frekvenční oblasti. Výsledky všech parametrů jsou srovnávány s hodnocením dvou lékařů.Z každé oblasti zpracování je přestavena jedna míra reprezentující danou oblast. Navržené pa-rametry zohledňují různé projevy koktavosti a neplynulosti. V časové oblasti je to průměrná délka ticha. Tento parametr dosahuje korelačního koeficientu s hodnocením lékařů 0,793.

Tomáš Lustyk 61

Obrázek 5: Výsledky nejlepšího nastavení adaptivního prahu.

Ve frekvenční oblasti je představen parametr průměrný počet změn bayesovského detektoru. Korelační koeficient s hodnocením lékařů je -0,782.V práci [1] uvedeno více než dvacet měr pro analýzu neplynulých promluv. Systémy kombi-nující více parametrů dosahují lepších výsledků než jednotlivé míry samostatně. Nejlepšího výsledku dosahuje systém kombinující 11 měr v časové i frekvenční oblasti. Systém správně odhaduje stupeň neplynulosti pro 63% jedinců, přičemž chyba odhadu o více než jednu třídu se vyskytuje pouze u 1,7% mluvčích.Také je představena modifikace určování počtu změn bayesovského detektoru pomocí adap-tivního prahu. K regulaci hodnoty prahu je použita směrodatná odchylka maxim ve výstupu bayesovského detektoru. Nejvyšší dosažený korelační koeficient s hodnocením lékařů je -0,726.

Poděkování

Děkujeme MUDr. M. Hrbkové a Dr. Ing. J. Vokřálovi z Foniatrické kliniky 1. LF UK a VFN za poskytnutí signálů. Tento výzkum byl podporován z výzkumného záměru MŠMTMSM6840770012 „Transdisciplinární výzkum v oblasti biomedicínského inženýrství“, a grantů GAČR 102/08/H008 “Analýza a modelování biomedicínských a řečových signálu”, GAČR 102/08/0707 “Rozpoznávání mluvené řeči v reálných podmínkách“ a grantem Stu-dentské grantové soutěže ČVUT č. SGS10/180/OHK3/2T/13 „Hodnocení poruch hlasu a řeči“.

Reference

[1] Bergl P., Objektivizace poruch plynulosti řeči. Disertační práce, Fakulta elektrotech-nická, ČVUT v Praze, 2010.

[2] Bergl P., Čmejla R., Hrbková M., Černý L.: Systém pro automatické hodnocení koktavosti. Automatizace 2010, roč. 53, č. 1-2, s. 49-52. ISSN 0005-125X.

[3] Lustyk T., Analýza neplynulé řeči. Diplomová práce, Fakulta elektrotechnická, ČVUT v Praze, 2010.

62 Tomáš Lustyk

Analýza subjektivního hodnocení dětského věku dle promluv

Martina Nejepsová

České vysoké učení v Praze, Fakulta elektrotechnická [email protected]

Abstrakt: Subjektivní hodnocení dětského věku bylo prováděno na základě pořízených zvukových záznamů od dětí ve věku 3-12 let. Hodnocení bylo prováděno logopedy i laiky z běžné populace a porovnáváno s objektivním hodnocením pomocí regresní analýzy.

1. Úvod

Subjektivní poslechové testy se používají především k hodnocení kvality zvukových záznamů, akustické kvality hudebních nástrojů, testování přenosové kvality elektroakustických zařízení, atd. [1]. Tyto testy jsou prováděny individuálně nebo skupinově při zajištění stejných poslechových podmínek. Dále je třeba zajištění opakovatelnosti testu. Jsou založeny na psychoakustice, subjektivním vnímání a posuzování, respektive na zkušenostech hodnotitelů s posuzovanými subjekty. Již ze zmíněných faktů vyplývá, že příprava a následně i vyhodnocování subjektivních poslechových testů je velmi náročné. Z toho důvodu, pokud je to možné, se přechází na objektivní hodnocení dle daných kritérií, které je nezávislé na opakování. V tomto projektu jsou vzájemně porovnávány subjektivní poslechové testy a objektivní hodnocení pří klasifikaci zvukových nahrávek a jejich kategorizaci. 2. Databáze

Pro analýzu byl použita databáze nahrávek dětí ve věku 3-12 let pořízených v mateřských a základních školách v Praze, která obsahuje promluvy od 103 chlapců a 90 dívek bez poruch hlasu či řeči.Každý z nich namluvil pro databázi slova uvedená v Tabulce 1. Histogram věkových kategorií je na Obrázku 1.

Věk [rok]

Obrázek 1: Počet promluv v jednotlivých věkových kategoriích

2 3 4 5 6 7 8 9 10 11 12 130

10

20

30

Martina Nejepsová 63

předškolní děti školní děti babička silnice čokoláda ptáček košile škola hodiny sluníčko letadlo zmrzlina koník vlak maluje prší zpívá

Tabulka 1: Slova v databázi

3. Objektivní hodnocení

V tomto projektu bylo objektivní hodnocení prováděno algoritmem M5 v prostředí WEKA [4]. Atributy, které byly použity v regresním stromu jsou fundamentální frekvence, první a druhý formant u samohlásek A,E a O, spektrální těžiště, spektrální směrodatná odchylka, spektrální zešikmení a špičatost u sykavek S a Š [3]. Fundamentální frekvence a první dva formanty byly získány v programu PRAAT [5] a kontrolovány v programu WAVESURFER [6] při analýze samohlásek z difónů LA, LE a LO vyjmutých ze slov v databázi. Spektrální koeficienty jsou počítány vzorci pro výpočet spektrálních momentů v prostředí Matlab [7] z vyjmutých sykavek ze slov „silnice“a „košile“ v databázi. 4. Subjektivní poslechové testy Pro subjektivní hodnocení byly připraveny dva poslechové testy TEST 1 a TEST 2, které každý hodnotitel provedl individuálně [2]. V testech byly užity nahrávky mluvčích ve věku 3-7 let z databáze. Pro oba testy byly nahrávky shodné, pouze se zaměněným pořadím. Hodnotitelé po jednorázovém vyslechnutí nahrávky zapsaly do předpřipraveného formuláře odhadovaný věk mluvčího z nahrávky. Klasifikace byla prováděna na škále v rozmezí věkových kategorií nahrávek z databáze – tedy 3-7 let. Testy TEST 1 a TEST 2 byly prováděny v časovém odstupu jeden měsíc. Subjektivní poslechový test provedlo 6 hodnotitelů – 4 dobrovolníci jako výběr z běžné populace a 2 zástupci z Foniatrické kliniky v Praze viz Tabulka 2.

4 dobrovolníci z různých oborů a zaměření Hodnotitel 2-I žena 33 let Hodnotitel 2-II muž 40 let Hodnotitel 2-III žena 59 let Hodnotitel 2-IV žena 41 let

2 zástupci z Foniatrické kliniky v Praze Hodnotitel 3-I žena 32 let Hodnotitel 3-II žena 34 let

Tabulka 2: Hodnotitelé 5. Statistické hodnocení Věk mluvčích byl znám jako počet roků a měsíců od narození do doby pořízení nahrávky. Maximální chyba od biologického věku je tedy ±15 dní. Pro hodnocení byl skutečný věk mluvčích dán počtem roků zaokrouhleným na 2 desetinná místa.

64 Martina Nejepsová

Hodnotitelé odhadovali věk mluvčích v počtu celých roků, výjimečně s upřesněním na ±0,25 či ±0,5 roku. Objektivní hodnocení je označeno číslem hodnotitele 1-I, hodnotitelé subjektivních poslechových testů jsou označeny dle Tabulky 2. 5.1 Odchylky hodnocení od skutečného věku mluvčího První porovnání úspěšnosti v hodnocení věku dětí dle jimi namluvených promluv bylo sledování odchylek od skutečné hodnoty. Z průměrné i jednostranných odchylek odhadovaného věku od skutečného (viz Tabulka 3 a Tabulka 4) je patrné, že hodnotitelé klasifikují promluvy do nižších věkových kategorií než jsou mluvčí ve skutečnosti. Pouze při objektivním hodnocení získáme kladnou průměrnou odchylku. Průměrná absolutní odchylka je v Tabulce 5.

číslo hodnotitele 1-I 2-I 2-II 2-III 2-IV 3-I 3-II

ø odchylka [rok] 1,49 -0,90 -0,16 -0,56 -0,22 -1,24 -0,88 pořadí 7 5 1 3 2 6 4

Tabulka 3: Hodnocení průměrné odchylky v TEST 1


ø kladná odchylka [rok] 1,83 1,27 1,18 1,17 0,91 0,72 0,58 ø záporná odchylka [rok] -0,85 -1,38 -1,31 -1,10 -1,00 -0,85 -1,52

rozdíl ø odchylek [rok] 0,98 -1,18 -0,13 -0,53 -0,57 -2,67 -2,18 pořadí 4 5 1 2 3 7 6

Tabulka 3: Hodnocení kladných a záporných odchylek samostatně v TEST 1


ø absolutní odchylka [rok] 1,71 1,26 1,06 0,90 0,74 1,35 1,02 pořadí 7 5 4 2 1 6 3

Tabulka 5: Hodnocení absolutní průměrné odchylky v TEST 1 Pro jednotlivé věkové kategorie jsme sledovali úspěšnost jednotlivých skupin hodnotitelů (viz Tabulka 6).

skupiny hodnotitelů 1 2 3

3roky - ø odchylka [rok] 2,72 0,50 0,50 4roky - ø odchylka [rok] 2,20 -0,21 0,00 5let - ø odchylka [rok] 1,43 -0,67 -0,84 6let - ø odchylka [rok] 1,19 -0,91 -1,48 7let - ø odchylka [rok] 1,15 -0,21 -1,89

Tabulka 6: Hodnocení průměrné odchylky jednotlivými skupinami v TEST 1 5.2 Míra korelace odhadů Dále byla sledována míra korelace jednotlivých skupin hodnotitelů se skutečným věkem mluvčích (viz Tabulka 7) a hodnotitelů subjektivních poslechových testů s objektivním hodnocením (viz Tabulka 8).

Martina Nejepsová 65


míra korelace 0,45 0,80 0,68 Tabulka 7: Míra korelace se skutečným věkem v TEST 1


míra korelace XXX 0,42 0,43 Tabulka 8: Míra korelace hodnotitelů s objektivním hodnocením v TEST 1

5.3 Odchylka opakovaného hodnocení Při vzájemném porovnávání výsledků z opakování subjektivního poslechového testu je patrná velmi nízká vzájemná odchylka odhadů zástupců z Foniatrické kliniky v Praze (viz Tabulka 9). Od dalších hodnotitelů zatím nebylo uskutečněno druhé hodnocení nahrávek, neboť neuplynula dostatečná doba od hodnocení prvního testu. číslo hodnotitele 1-I 2-I 2-II 2-III 2-IV 3-I 3-II

ø odchylka [rok] 0,00 -0,37 XXX XXX XXX -0,08 0,10 ø absolutní odchylka [rok] 0,00 0,66 XXX XXX XXX 0,23 0,21

Tabulka 9: Odchylka hodnocením v TEST 1 a v TEST 2 6. Závěry Cílem projektu bylo zhodnocení úspěšnosti objektivního hodnocení vůči subjektivnímu hodnocení logopedy i laiky. Po porovnání výsledků je možné konstatovat, že pokud bude pro objektivní hodnocení regresní analýzou dostatek vstupních dat a vhodných klasifikačních atributů, bude algoritmus klasifikovat s velkou přesností. V opačném případě dosahují účastníci subjektivních poslechových testů větší přesnosti. Objektivní hodnocení je nezávislé na opakování. Při opakování subjektivního poslechového testu byla sledována přesnost hodnocení jednotlivými hodnotiteli. S průměrnou odchylkou ±0,1 roku se shodovaly výsledky obou poslechových testů u zástupců z Foniatrické kliniky v Praze. Z řad dobrovolníků provedl TEST 2 pouze jeden uchazeč a to s odchylkou ±0,4 roku. 7. Poděkování

Tento výzkum byl podporován z grantu GAČR 102/08/H008 “Analýza a modelování biomedicínských a řečových signálu”, SGS10/180/OKH3/20/13 „Hodnocení poruch hlasu a řeči“ a výzkumného záměru MSM6840770012 „Transdisciplinární výzkum v oblasti biomedicínského inženýrství II“. Reference [1] Otčenášek, Z.: O subjektivním hodnocení zvuku, AMU, Praha 2008 [2] Melka, A: Základy experimentální akustiky, AMU, Praha 2005 [3] Janda, J.: Studie věkově závislých akustických parametrů v dětské řeči, Studie k odborné rozpravě, ČVUT, Praha 2010 [4] http://www.cs.waikato.ac.nz/ml/weka/ [5] http://www.fon.hum.uva.nl/praat/ [6] http://www.speech.kth.se/wavesurfer/ [7] http://www.mathworks.com/

66 Martina Nejepsová

Příprava a analýza Českého Web 1T 5-gram korpusupro použití v jazykovém modelu

Václav Procházka


Abstrakt: V této práci je popsán postup analýzy českého Web 1T 5-gramkorpusu. Korpus byl analyzován a byly vyhodnoceny jeho základní charak-teristiky před a v průběhu zpracování. Při zpracování byl slovník korpusufiltrován různými metodami, tak aby pokud možno obsahoval pouze smys-luplná slova. Z pročištěného korpusu byly vygenerovány jazykové modely proLarge Vocabulary Continuous Speech Recognition (LVCSR) a spočítána jejichperplexita. Pro srovnání stejnými filtrovacími postupy byl také zpracovaný 5-gramový korpusu založený na SYN2006 korpusu který sestavil Český národníkorpus (ČNK).

1. Úvod

Tvorba jazykových modelů je běžnou součástí tvorby rozpoznávače spojité řeči (LargeVocabulary Continuous Speech Recognition (LVCSR) LVCSR). Tento článek dále popisujeanalýzu nedávno zveřejněného zdroje českých webových n-gramů [6] a výsledky porovnávás daty získanými z Českého národního korpusu (ČNK) [2].V současnosti asi nejčastějším zdrojem textů pro jazykové modelování jsou knihy, noviny,časopisy či jiná publicistická tvorba. Získávání textů z těchto zdrojů je obtížné a zdlou-havé, a to hlavně ze dvou důvodů: texty nejsou volně k dispozici nebo jsou chráněnévlastnickými právy a není možně je uchovávat. Alternativní možností, jak získat velkémnožství textových dat, je potom internet. Množství dostupných textů se spolu s růstempočtu webových stránek neustále zvětšuje, a tím také roste atraktivita internetu jako jejichzdroje, přes mnohé komplikace při zpracování, které internet jako zdroj textů přináší.

2. Český Web 1T 5-gram korpus

Český Web 1T 5-gram korpus [6] je souborem unigramů až 5-gramů sestavených firmouGoogle z webových stránek vydaný v Linguistic Data Consortium (LDC) [4]. Použitímtohoto korpusu je možné ušetřit mnoho náročné práce se sběrem korpusových dat z WEBo-vých stránek. Velké množství dat (∼136 miliard slov) je již předzpracováno do podobyn-gramů s následujícími vlastnostmi [6]:

• kódování n-gramů je sjednocené (UTF8),

• obsahuje n-gramy 1 až 5 řádu,

Václav Procházka 67

• (X)HTML značky jsou odstraněné,

• slova obsahují malá i velká písmena,

• speciální token <UNK> pro neznámá nebo neplatná slova,

• speciální token <S> resp. </S> pro začátek resp. konec věty,

• neobsahuje n-gramy jdoucí přez hranici věty, tzn. <S> se vyskytuje pouze na za-čátku a </S> pouze na konci n-gramu,

• čísla o více než 16 číslicích a slova delší než 32 znaků jsou nahrazena tokenem<UNK>,

• slova se znaky z neevropských jazyků, neplatnými UTF8 znaky a kontrolními ASCIIznaky jsou také nahrazena tokenem <UNK>,

• slova vyskytující se méně než 40 krát jsou nahrazena tokenem <UNK>

• a všechny n-gramy které po předchozích operacích měly výskyt menší než 40 byliz korpusu odstraněny úplně.

Prvotní náhled do slovníku ukázal, že přes výše uvedené rozsáhlé předzpracování sev tomto korpusu vyskytuje velké množství nevhodných nebo dokonce nesmyslných slovjako jsou slova obsahující písmena i číslice, slova z jiných evropských jazyků, slova sklá-dající se z jiných znaků než písmen a číslic (interpunkce), či náhodné shluky písmen. Proúčely tvorby jazykových modelů pro rozpoznávání řeči je nutné tato slova odstranit. Ab-sence n-gramů jdoucích přez hranice vět a ořezání nejméně častých n-gramů, v některýchpřípadech omezuje nebo komplikuje možnosti dalšího zpracování.

3. Referenční SYN2006PUB 5-gram korpus

Jako referenční korpus byl použit SYN2006PUB 5-gram korpus, což je soubor n-gramůvygenerovaný ze SYN2006PUB [2] od ČNK [3]. Zdrojový korpus SYN2006PUB je syn-chronní korpus psané publicistiky o rozsahu 300 milionů textových slov. V porovnánís Web 1T 5-gram korpusem je tento korpus výrazně čistější. Tento n-gramový korpus mánásledující vlastnosti:

• obsahuje n-gramy 1 až 5 řádu,

• kódování ISO-8859-2,

• slova obsahují malá i velká písmena,

• žádné ořezání (cutoff), i slova s jedním výskytem jsou zachována,

• speciální token </s> označující hranici věty,

• obsahuje n-gramy jdoucí přes hranice vět

• a neobsahuje interpunkční znaménka.

Z nevhodných slov vyskytujících se v korpusu Web 1T 5-gram obsahuje tento korpusve větší míře pouze číselné výrazy, jiná nevhodná slova jako například náhodné shlukypísmen jsou méně časté.

68 Václav Procházka

původní korpus pročištěný korpuspo prvním kole po druhém kole

unigramy 9 786 424 1 804 682 954 771trigramy 117 264 988 47 376 975 8 567 409

Tabulka 1: Počet unikátních n-gramů českého Web 1T 5-gram korpusu.

původní korpus pročištěný korpuspo prvním kole po druhém kole

unigramy 2 554 028 1 799 0057 859 403trigramy 189 152 100 170 243 851 169 671 342

Tabulka 2: Počet unikátních n-gramů českého SYN2006PUB 5-gram korpusu.

4. Zpracování korpusů

Zpracování korpusů se skládá z několika na sobě více či méně závislých krocích. Některékroky byly aplikovány jen na jeden z korpusů nebo pouze v druhém kole zpracování.Jednotlivé kroky zpracování byly aplikovány v tomto pořadí:

• slova byla převedena na malá písmena,

• vybraná nevhodná slova nebo tokeny nahrazena tokenem <UNK>:

– všechna slova, která obsahovala jiné znaky než písmena než ty z české abecedya interpunkce,

– čísla, číselné výrazy (např.: +4 nebo 1,5) a data,– slova označená při kontrole pravopisu nástrojem Aspell [1] jako nevhodná (jen

v druhém kole),– ručně vybraná nevhodná slova o třech a méně znacích, například zkratky nebo

náhodné kombinace vyskytující se mezi nejčastějšími 60 tisíci slovy (jen v dru-hém kole),

• interpunkce byla vynechána bez náhrady (jen Web 1T 5-gram, SYN2006PUB 5-gram interpunkci neobsahuje),

• nahrazení tokenu </s> tokeny </s> a <s> (jen SYN2006PUB 5-gram , Web 1T5-gram již oba tokeny obsahuje),

• totožné n-gramy vzniklé v předchozích krocích byly sjednoceny a jejich četnostisečteny.

Operace vynechání interpunkce ve svém důsledku způsobuje, že řád zpracovávaného n-gramu je snížen. Aby byla zachována možnost tvorby trigramových jazykových modelů,na vstupu byly 5-gramy. N-gram který byl zredukován více než na trigram byl zahozenzcela. Při zpracování SYN2006PUB 5-gram byl tento krok zbytečný a byl zcela vynechán,čímž bylo možné od začátku pracovat pouze s trigramy, a tím zrychlit celý proces.


bigramový model trigramový modelslovník OOV rate cut off perplexity bigramy perplexity trigramy60 000 17.76% 40 49 507 11 437 940 162 509 35 166 960120 000 12.14% 40 126 001 14 904 654 255 144 42 125 506180 000 9.44% 40 153 682 16 522 260 664 866 43 882 665240 000 7.63% 40 206 831 17 666 327 935 077 45 443 605

Tabulka 3: Parametry jazykových modelů natrénovaných z pročištěného českého Web 1T5-gram korpusu po prvním kole.

cut off 1 bigramový model trigramový modelslovník OOV rate cut off perplexity bigramy perplexity trigramy60 000 10.93% 1 932 17 587 483 928 65 783 406120 000 6.36% 1 1 226 29 928 179 1 197 96 243 374180 000 4.18% 1 1 432 35 987 677 1 400 108 421 863240 000 3.06% 1 1 571 39 963 558 1 521 115 381 912

Tabulka 4: Parametry jazykových modelů natrénovaných z pročištěného SYN2006PUB5-gram korpusu po prvním kole.

SYN2006PUB 5-gram korpus obsahuje pouze jeden token pro hranici věci. Jelikož dalšízpracování pomocí Hidden Markov Model Toolkit (HTK) vyžaduje token pro začátek i ko-nec věty, byl tento token nahrazen dvojicí tokenů </s> a <s>. Tím došlo ke zvýšení řádun-gramu resp. k jeho následnému rozdělení na dva n-gramy se stejným řádem a počtemvýskytů dle původního n-gramu.Po prvním kole zpracování, tedy bez filtrování slovníku nástrojem Aspell a ručního odstra-nění krátkých slov, ve slovníku Web 1T 5-gram korpusu stále zbylo mnoho špatných nebonevhodných slov, například cizí slova, slova s překlepy, slova s odstraněnou diakritikounebo i jen pouze nesmyslné shluky písmen (kód výrobku, jméno firmy). Tyto problémyjsou řešeny v rámci druhého kola filtrování.Trigramy a bigramy připravené v prvním kole pak byly použity jako základ pro dalšíanalýzy a pro tvorbu jazykových modelů pomocí HTK [11]. Počty vybraných unikátníchn-gramů jsou shrnuté v tabulce 1. Perplexity na jazykových modelech postavených z tri-gramů každého korpusu jsou shrnuty v tabulkách 3 a 4 převzatých z [9]. Perplexita bylapočítána na podmnožině přepisů z České SPEECON databáze [8], [7]. Tato podmnožinaobsahuje 148 557 slov v 14 914 větách. Perplexity napočítané na jazykových modelechz SYN2006 5-gramů mají typický průběh a pohybují se v rozumných rozsazích. Jazykovémodely postavené z českého Web 1T 5-gram korpusu mají velmi velké perplexity a také vy-kazují výrazně větší četnost mimoslovníkových slov (OOV - Out Of Vocabulary (OOV)).Podrobnější výsledky tohoto kola zpracování jsou shrnuty v článku [9].Pro takto vysoké perplexity Web 1T 5-gram korpusu existuje několik možných vysvětlení.Po provedeném filtrování v korpusu zůstalo mnoho cizích slov či slov s chybějící diakri-tikou. Typické vlastnosti internetových stránek mohou být dalšími důvody. Například jemožné se setkat se stránkami s částečné identickým kódem (hlavičky, patičky nebo menu)nebo dokonce stránkami téměř identickými viz. stránky internetových obchodů. Odstra-nění prvních dvou zmíněných problémů, cizích slov a slov s chybějící diakritikou byloprovedeno v druhém kole zpracování.


bigramový model trigramový modelslovník OOV rate cut off perplexity bigramy perplexity trigramy60 000 13.02% 40 337 226 1 512 233 120 476 1 765 425120 000 8.14% 40 570 781 1 797 817 249 281 1 991 280180 000 5.75% 40 525 619 1 920 150 380 179 2 079 416240 000 4.48% 40 1 184 050 1 985 640 507 689 2 123 533

Tabulka 5: Parametry jazykových modelů natrénovaných z pročištěného českého Web 1T5-gram korpusu po druhém kole.

cut off 1 bigramový modelslovník OOV rate cut off perplexity bigramy60 000 10.77% 1 1 727 30 205 968120 000 6,71% 1 2 215 38 181 993180 000 4,78% 1 2 578 41 981 434240 000 3,77% 1 2 848 44 146 507

Tabulka 6: Parametry jazykových modelů natrénovaných z pročištěného SYN2006PUB5-gram korpusu po druhém kole.

V druhém kole zpracování byly kromě operací nových zachovány i všechny operace před-chozího kola. Na rozdíl od prvního kola, všechny pročišťovací operace, které ze své pod-staty nemusí pracovat na n-gramech požadovaného řádu, byly provedeny na unigramech.Jedná se především o všechny operace odstranění nevhodných slov a tokenů, které jsouv n-gramech nahrazeny speciálním tokenem <UNK>. Zpracování a čištění slovníku v prv-ním kroku s následující kontrolou existence slova ve vyčištěném slovníku je při filtrovánín-gramů nejen univerzálnější, ale také řádově rychlejší než provádění obou operací sou-časně na n-gramech. Výsledná velikost slovníku po druhém kole filtrování je v tabulce 1.Počty n-gramů a napočítané perplexity počítané na stejných datech jako v kole prvnímjsou v tabulkách 5 a 6. U modelů se slovníkem o velikosti 60000 slov došlo podle oče-kávání k poklesu OOV. U Web 1T 5-gram modelů byl pokles výrazný, u SYN2006PUB5-gram modelů jen nepatrnému. Pro slovníky ostatních velikostí OOV modelů Web 1T5-gram korpusu rovnoměrně klesalo, oproti tomu u modelů ze SYN2006PUB 5-gram OOVmírně vzrostlo, důvodem je pravděpodobně příliš agresivní filtrování, které odstranilo isprávná slova. Perplexity trigramových modelů postavených z Web 1T n-gramů se trochuzlepšily, nicméně úrovně perplexit modelů ze SYN2006PUB 5-gram zdaleka nedosáhly.Perplexity bigramových modelů z Web 1T n-gramů výrazně vzrostly, to je pravděpo-dobně způsobeno vlastnostmi korpusu: např. ořezáním, obsahem interpunkce a krokemvypuštění interpunkce, kdy pak 5-gramy zredukované na bigramy již dostatečné dobřenereprezentují původní data. Drobný vzrůst perplexit bigramů ze SYN2006PUB 5-gramkorpusu má stejnou příčinu jako vzrůst OOV, tj. příliš agresivní filtrování.

5. Popis použitých nástrojů

Pro stavbu modelů a počítání perplexit a OOV byly použity veřejně dostupné balíčkynástrojů Hidden Markov Model Toolkit (HTK) a Stanford Research Institute LanguageModeling Toolkit (SRILM).


HTK [11] je jeden z balíčků nástrojů pro práci se skrytými Markovovými modely (HMM).Především je určen pro aplikaci v oblasti rozpoznávání řeči, ale používá se i pro různédalší aplikace jako je syntéza řeči, rozpoznávání textu, nebo další neřečové aplikace např.sekvencování DNA [11]. V oblasti rozpoznávání řeči je tento balíček nástrojů předevšímzaměřen na vytváření akustických modelů založených na Hidden Markov Model (HMM)avšak také obsahuje nástroje pro vyváření jazykových modelů. Nástroji HTK je možnénatrénovat akustické a jazykové modely a vyhodnotit jejich kvalitu a funkčnost v LVCSRnebo samostatně. Část nástrojů pro práci s jazykovými modely byla použita předevšímv prvním kole zpracování n-gramů.SRILM [10] [5] je balíček nástrojů pro tvorbu a použití jazykových modelů, předevšímv oblastech rozpoznávání řeči, statistického značkování a segmentace či strojového pře-kladu. Pro práci s jazykovými modely je tento nástroj silnější než HTK, mimo jiné umož-ňuje pracovat přímo s n-gramy v textové podobě a podporuje více vyhlazovacích technik.Jazykové modely ve formátu ARPA vytvořené tímto nástrojem je možné přímo použíti v HTK. Tyto nástroje byly požity pro tvorbu a analýzu jazykových modelů v druhémkroku.

6. Závěr

V této práci byla popsána analýza a další zpracování Web 1T 5-gram korpusu od Google.Výsledky těchto analýz byly průběžně porovnávány s výsledky na SYN2006PUB 5-gramkorpusu, který byl použit jako referenční.

• Český Web 1T n-gram corpus se zdá být dobrým počáteční zdrojem pro prácí s datyz webu. Základní problémy spojené se získáváním textů z webu jako sjednocenékódování, odstranění (X)HTML značek a jednoduché filtrování jsou v tomto korpusujiž vyřešené.

• Perplexita modelů vytvořených z n-gramů českého Web 1T 5-grams po první fázebyla velmi vysoká, více než 105 v porovnání s výsledky perplexit SYN2006PUB 5-gram korpusu které byly mezi 900 a 1600. Takto vysoké perplexity českého Web 1Tkorpusu znamenají, že filtrování provedené v prvním kole nebylo dostatečné.

• V druhém kole bylo filtrování rozšířeno o filtrování pomocí slovníků pro automatic-kou kontrolu pravopisu. Při tomto filtrování došlo k významnému poklesu slov veslovníku. OOV a perplexity trigramových modelů z českého Web 1T 5-grams kor-pusu se zlepšily, nicméně stále zůstávají příliš vysoké. Stejné charakteristiky modelůze SYN2006PUB korpusu se převážně mírně zhoršily, což odpovídá příliš agresiv-nímu automatickému filtrování.

• Další pokračování těchto pokusů s čištěním n-gramových korpusů a tvorby jazyko-vých modelů se budou soustředit na získání společného referenčního slovníku a ná-strojů pro jeho snadné rozšiřování.

Poděkování

Tento výzkum byl podporován z grantů GAČR 102/08/0707 “Rozpoznávání mluvené řečiv reálných podmínkách” resp. výzkumného záměru MŠMT MSM6840770014 “Výzkumperspektivních informačních a komunikačních technologií”.


Reference

[1] GNU Aspell. accesed: 17.11.2010, http://aspell.net/.

[2] Český národní korpus (Czech National Corpus) - SYN2006PUB. Institute of theCzech National Corpus FF UK, Prague, 2006. http://www.korpus.cz.

[3] Czech National Corpus, 2010. http://www.korpus.cz.

[4] Linguistic Data Consortium, 2010. http://www.ldc.upenn.edu.

[5] Stanford Research Institute Language Modeling Toolkit (SRILM), 2010. http://www.speech.sri.com/projects/srilm/.

[6] Brants, T.; Franz, A. Web 1T 5-gram, 10 European languages, version 1. LinguisticData Consortium, Web page, Philadelphia, 2009. http://www.ldc.upenn.edu.

[7] Iskra, D.; Grosskopf, B.; Marasek, K.; van den Heuvel, H.; Diehl, F.; Kiessling, A.SPEECON - speech databases for consumer devices: Database specification and va-lidation. In Proc. of LREC’02 May 2002.

[8] Pollák, P.; Černocký, J. Czech SPEECON adult database, Nov 2003. accesed:17.11.2010, http://www.speechdat.org/speecon/index.html.

[9] Procházka, V.; Pollák, P. Analysis of Czech Web 1T 5-gram Corpus and Its Compa-rison with Czech National Corpus Data. In LNAI 6231, (TSD 2010) 2010, P. Sojkaet al., Ed., Springer-Verlag Berlin Heidelberg, pp. 181–188.

[10] Stolcke, A. Srilm—an extensible language modeling toolkit. In In Proceedings of the7th International Conference on Spoken Language Processing (ICSLP 2002) 2002,pp. 901–904.

[11] Young et al., S. The Hidden Markov Model Toolkit (HTK), Version 3.4.1, 2009.accesed: 9.6.2010, http://htk.eng.cam.ac.uk.


Princip aplikace speciální arteterapeutické techniky

u dětí s vývojovou dysfázií

Barbora Richtrová

Univerzita Karlova v Praze, 1. lékařská fakulta

[email protected]

Abstrakt: Cílem výzkumného projektu je potvrzení hypotézy, že arteterapeutické techniky a

strategie zmírňují nebo eliminují symptomy vývojové dysfázie, především obtíže v oblasti

zpracování řečových signálů. Tato vývojová porucha řeči je „charakterizována specifickým

řečovým vývojem, který je především aberantní“. Široká symptomatologie zasahuje osobnost

celostně a zároveň jedinečně, kdy u každého klienta se symtomatologie projevuje v různých

oblastech a v různé míře.

1. Úvod

1.1 Arteterapie

„…..arteterapie je teoreticky zakotvené působení na člověka jako celek v jeho tělesné,

psychické a duševní realitě, v jeho vědomém i nevědomém snažení, sociálních a ekologických

vazbách, plánované ovlivňování postojů a chování pomocí umění a technik z umění

odvozených, s cílem léčby či zmírnění nemoci a integrace či obohacení osobnosti“ (Petzold,

2001, s. 86) Jedna z mnoha artetrapeutických definic, plně vystihuje předpoklad, že vhodně

zvolené arteterapeutické techniky a strategie působí nejen na konkrétní symptom, ale mají

velký vliv na člověka celkově.

1.2 Vývojová dysfázie

Vývojová dysfázie je centrální porucha komunikační schopnosti. Jedná se o specificky

narušený řečový vývoj s dvojí patofyziologií [3]: A) specifická centrální sluchová porucha,

B) všeobecná nespecifická korová porucha. Může se projevovat neschopností nebo sníženou

schopností verbálně komunikovat, a to i přesto, že podmínky pro vytvoření schopnosti

verbálně komunikovat jsou dobré. To znamená, že se: a) nevyskytují se závažné

psychiatrické, b) neurologické nálezy, c) inteligence je přiměřená, d) nevyskytuje se závažná

sluchová, zraková vada a e) sociální prostředí je podnětné.

Symptomatologie u vývojové dysfázie se projevuje ve všech jazykových úrovních, kdy

základní příčinou je porucha centrálního sluchového zpracování řeči. Rozvoj aktivní slovní

zásoby je pomalý, v řeči se objevují dysgramatismy. Obtíže se objevují v percepci,

diskriminaci, syntéze a analýze, v paměti (především krátkodobé fonologické), diferenciaci,

produkci atd. Nesnáze se objevují i ve vývoji orofaciální motoriky.

74 Barbora Richtrová

2. Výzkum

2.1 Cíl výzkumu

Cílem mého výzkumného projektu je potvrzení hypotézy, že arteterapeutické techniky a

strategie zmírňují nebo eliminují symptomy vývojové dysfázie, především obtíže v oblasti

zpracování řečových signálů. Tato vývojová porucha řeči je „charakterizována specifickým

řečovým vývojem, který je především aberantní“. Široká symptomatologie zasahuje osobnost

celostně a zároveň jedinečně, kdy u každého klienta se symtomatologie projevuje v různých

oblastech a v různé míře.

K uvedenému výzkumu mě vede jednak skutečnost, že neustále roste počet klientů s

diagnózou vývojová dysfázie. Terapie je velmi náročná a je otázkou, jestli terapeutické

působení by nebylo vhodné rozšířit o prostředky uměleckých terapií, které působí na klienta

celostně. Arteterapie je dalším možným nástrojem k rehabilitaci a celostnímu rozvoji

osobnosti, kdy osobní prožitek a vlastní zkušenost výrazně determinují rozvoj kognitivních

procesů. Právě tato oblast je pro klienty s vývojovou dysfázií zásadní vzhledem k jejich

řečovému vývoji.

2.2 Koherentní arteterapeutické systémy

Koherentní arteterapeutické systémy odvozují svou odlišnost od plně využívaného procesu

výtvarné tvorby: „Arteterapeutická teorie vychází z teorie umění a psychoterapeutických škol.

Arteterapie je součástí jasně definovaného léčebného (psychoterapeutického) procesu. Je

třeba odlišovat psychoterapii užívající artetechnik od arteterapie. Zatímco v psychoterapii

jsou artetechniky zařazovány cíleně a izolovaně, zpravidla proto, aby byl získán materiál pro

zpracování určitého tématu, v arteterapii jde o využití plnohodnotného kanálu pro komunikaci

a introspekci. Neverbální tvořivá činnost zde slouží nejen pro otevírání, ale i pro zpracování

témat.“ (Stiburek in Slavík, 2001, s.35).

Tvořivý proces přináší rozvoj osobnosti celkově. Na základě vlastního prožitku a zážitku

vznikají nové zkušenosti, tudíž bychom mohli říct, že tvůrčí proces stimuluje kognitivní

stránku jedince. Zároveň je člověku posilována i emocionální, sociální a komunikakční oblast.

Z výše uvedeného vyplývá předpoklad, že tvůrčí proces je vzhledem k symptomatologii

vývojové dysfázie velmi důležitý a účinný. Na základě symptomatologie vývojové dysfázie se

jeví účinější vést terapii multisenzorialně, kdy stěžejním terapeutickým kamenem pro můj

výzkumný projekt je propojení audio – hapticko – vizuálního cítění.

Ve svém projektu pracuji s technikami právě na tomto principu. Klíčovým se mi zdá být

haptické vnímání. Přičemž haptické vnímání neproudí jen skrze ruce, ale skrze jejich pohyb,

přes motoriku. Pohyb rukou a doteky rukou jsou důležité, zvláště dotýkání se sebe sama a

dotýkání se svého okolí. Když se nás něco dotýká nebo my se něčeho dotýkáme, tak nás to

vede k určitým reakcím, jsou to impulsy, které nás mohou „rozhýbat “.

Arteterapeutické techniky, které stojí na základě audio-hapticko-vizuálního spojení, nás

vedou k pohybu. Jsou vhodné pro stimulaci vývoje a rozvoje nejen intaktní populace, ale

především pro celostní vývoj a sebenahlížení klientů s vývojovou dysfázií, u kterých je toto

„rozhýbání“ nesmírně důležité vzhledem k jejich vývoji.

Barbora Richtrová 75

2.3 Technika práce na „Hliněném poli“

S technikou práce na „Hliněném poli“ jsem se seznámila v ateliéru Extraart při víkendovém

semináři, který vedl arteterapeut M. Weigele. Protože mě technika velmi zaujala a chtěla jsem

jít více do hloubky, absolvovala jsem čtrnáctidenní stáž u M. Weigeleho v Německu.

Hlavním představitelem této techniky je v Německu arteterapeut Heinz Deuser. Technika je

určena pro děti i dospělé, kteří mají různé psychologické obtíže. Stáž, kterou jsem

absolvovala v Německu vede žák Heinze Deusera, M. Weigele. Pracuje především na klienty

s ADHD syndromem. Zásluhou profesora H. Deusera je, že práce na Hliněném poli má již

ustálenou a srozumitelnou metodiku. V roce 1996 byl vytvořen „Spolek pro tvarování“, který

mimo mnoho jednotlivých projektů uspořádal už tři sympozia věnovaná této technice.

Doposud největší sympozium se konalo v říjnu 2007 v Freiburg/Breisgau s názvem „Im

Greifen sich begreifen“.

Základním principem práce na Hliněném poli je haptické dění, kdy se něčeho dotýkáme a to

se nás vnitřně dotýká. Zažíváme sami sebe skrze něco jiného a zažíváme ostatní skze nás. V.v

Weizsäcker nazval tuto propojenost jako „Gestaltkreis“ a J.v Uexküll jako „Funktionskreis“.

V haptickém dění se omezuje představa o světě i představa o sobě samých na naše bazální

myšlení. Díky tomuto bazálnímu myšlení tvarujeme a organizujeme sami sebe.

Arteterapeutické techniky a strategie v plné míře využívají princip zkušenosti, prožívání, hry,

tvořivosti. Činnostní pojetí a vlastní zážitek klienta posilují osobnost v mnoha stránkách.

Vzhledem k projektu vyhledáváme, upravujeme a hledáme zobecnění právě pro ty strategie a

techniky, které především stimulují zpracování řečového signálu a kognitivní stránky

osobnosti důležité pro komunikační schopnosti a dovednosti (percepce, diferenciace,

produkce atd.).

2.4 Plán výzkumu

Vytvořila jsem následující skupiny:

1) Děti ve věku 3-5 let

2) Děti ve věku 5-7 let.

Ke každému věkovému intervalu jsou vytvořeny tři podskupiny:

1) klienti s vývojovou dysfázií – terapie „klasická“

2) klienti s vývojovou dysfázií – terapie „arteterapeutická“

3) kontrolní skupina.

V současné fázi výzkumu pokračuje sběr dat. Metodika i konkrétní techniky se v průběhu

projektu upravují. Vzhledem k danému projektu vyhledáváme, upravovujeme a zobecňujeme

právě ty strategie a techniky, které především stimulují zpracování řečového signálu a

kognitivní stránky osobnosti důležité pro komunikační schopnosti a dovednosti (percepce,

diferenciace, produkce atd.).

76 Barbora Richtrová

3. Závěry

Neustále roste počet klientů s diagnózou vývojová dysfázie. Terapie je velmi náročná a je

otázkou, jestli terapeutické působení by nebylo vhodné rozšířit o prostředky uměleckých

terapií, které působí na klienta celostně. Na základě symptomatologie vývojové dysfázie se

jeví účinější vést terapii multisenzorialně, kdy stěžejním terapeutickým kamenem pro můj

výzkumný projekt je propojení audio – hapticko – vizuálního cítění. Klíčovým se mi zdá být

haptické vnímání. Přičemž haptické vnímání neproudí jen skrze ruce, ale skrze jejich pohyb,

přes motoriku. Výsledkem projektu by v ideálním případě mělo být potvrzení hypotézy, že

vybrané arteterapeutické techniky a strategie zmírňují nebo eliminují symptomy vývojové

dysfázie.

4. Poděkování

Tento výzkum byl podporován z grantu GAČR 102/08/H008 “Analýza a modelování

biomedicínských a řečových signálů”.

Reference

[1] DEUSER, H. a kol. Der haptische Sinn. Hamburg: Verein für Gestaltbildung e. V., 2009.

[2] DEUSER, H. Im Greifen sich begreifen. Keutschach: Gerhild Tschachler-Nagy, 2007.

[3] DLOUHÁ, O. Vývojové poruchy řeči. Praha: Publisher, 2003, ISBN 80-239-1832-X.

[4] MIKULAJOVÁ, M; RAFAJDUSOVÁ, I. Vývinová dysfázia. Bratislava: 1993.

[5] SLAVÍK, J. Umění zážitku, zážitek umění. I. a II. díl. Praha: UK Pedagogická fakulta,

2001 a 2004.

[6] ŠKODOVÁ, E.; JEDLIČKA, I. Základy klinické logopedie. Praha: Portál, 2004.

[7] VEREIN FÜR GESTALTBILDUNG e. V. (Hrsg.) Der haptische Sinn. Hamburg, 2009.

www.artefiletika.cz

www.arteterapie.cz

http://is.muni.cz/th/10957/pedf_d/007_kapitola_4.doc

www.tonfeld.de

www.tonfeld-ammerbuch.de

Barbora Richtrová 77

Evaluation of dysphonia in untreated Parkinson’s disease

Jan Rusz Czech Technical University, Faculty of Electrical Engineering

[email protected]

Abstract: Parkinson’s disease (PD) is a chronic neurodegenerative disorder characterized by progressive lost of dopaminergic neurons in the substantia nigra. In addition to the most ostensible motor symptoms such as tremor, many patients with PD develop non-motor deficits in speech known as hypokinetic dysarthria. The disorders of voice and speech in Parkinson’s disease result from the involvement of one or more speech subsystems including respiration, phonation, articulation, and prosody. The deficits in phonation characterized as dysphonia is one of the most salient features of PD speech impairment. This study thus investigates the signs of dysphonia in untreated patients with PD in comparison with healthy control participants. Speech data was collected using sustained phonations from 46 Czech native speakers, 24 with PD. Then, 7 representative measurements of dysphonia including jitter, shimmer, harmonic-to-noise ratios, recurrence period density entropy, detrended fluctuation analysis, and pitch period entropy were selected. The acoustic analyses for its assessment were performed using Praat and Matlab. Significant differences between both groups were found in all measurements except detrended fluctuation analysis. The results of this study confirm that dysphonia is sign of PD-related vocal impairment from the early stage of disease.

1. Introduction The progressive lost of dopaminergic neurons in Parkinson’s disease (PD) is associated with a variety of motor deficits, and non-motor deficits such as disorders of speech, mood, behaviour, thinking, and sensation [1]. As the second most common neurodegenerative disorder after Alzheimer’s disease, PD affects a large part of worldwide population. PD is estimated to affect the population over the age of 50; only approximately 10% of patients report symptoms before the age of 40 [2]. Moreover, PD affects 1.6% of persons after the age of 65 [3]. Moreover, statistics of the number of persons with PD are expected to increase because the worldwide population is growing older [4]. Several studies have shown that deficiencies in speech affect approximately 75-90% people with PD [5, 6]. The most salient features of PD speech impairment include deficits in the production of vocal sounds (dysphonia), and problems with motor speech disorder (dysarthria) [6, 7]. There are several vocal tests using various speech samples to assess the extent of these PD-related symptoms. These among others include sustained phonation where the participant produces a single vowel and holds it at a comfortable pitch and loudness as constantly and long as possible [8]. Previous research have found statistical differences between PD patients in comparison to healthy control (HC) participants using traditional dysphonia measurements of jitter, shimmer, harmonic-to-noise ratio (HNR), and noise-to-harmonics ratio (NHR) [9, 10]. Recently, novel non-standard measurement methods have been proposed to assess dysphonic symptoms. As randomness and noise are inherent to vocal production, methods based on nonlinear dynamical system theory such as recurrence period density entropy (RPDE) and

78 Jan Rusz

detrended fluctuation analysis (DFA) are good candidates for their ability to detect general voice disorders [11]. Another measurement of dysphonia called pitch period entropy (PPE), is a robust measurement sensitive to observed changes in speech specific to PD [12]. The aim of this study was to analyze the extent of dysphonia occurrence, using the various traditional and novel measurements, in early stages of PD, before symptomatic pharmacotherapy treatment is administered. The remainder of this paper is organized as follows: The section “Methods” presents the subjects, recording and data, measurement method, and statistics used. In the section “Results” the result obtained are shown. The sections “Discussion and Conclusion” discuss the findings and contain a summary of results for future work. 2. Methods 2.1. Subjects From 2007 to 2009, a total of 46 Czech native participants were recruited for this study, 24 of these subjects (20 men and 4 women) were diagnosed with PD. They were examined in the drug-naive state, before the symptomatic treatment was started. The age of PD subjects ranged from 34 to 83 years ([mean age in years 60.92 (±SD 11.24)]), and the duration of PD prior to recording ranged from 6 months to 7 years ([mean duration of PD in months 31.29 (±SD 22.25)]). All these patients fulfilled diagnostic criteria for PD, and were evaluated with the Unified PD Rating Scale III (UPDRS III), and the modified Hoehn and Yahr (HY) staging scale ([mean UPDRS III score 17.42 (±SD 7.14), mean modified HY scale 2.19 (±SD 0.48)]).

Parameter Description

Jitter:local (%) Average absolute difference between consecutive periods, divided by the average period.

Jitter:RAP (%) Relative Average Perturbation, the average absolute difference between a period and the

average of it and its two neighbours, divided by the average period.

Jitter:PPQ5 (%) Five-point Period Perturbation Quotient, the average absolute difference between a period

and the average of it and its four closest neighbours, divided by the average period.

Jitter:DDP (%) Average absolute difference between consecutive differences between consecutive periods,

divided by the average period.

Shimmer:local (%) Average absolute difference between the amplitudes of consecutive periods, divided by the

average amplitude.

Shimmer:local (dB) Average absolute base-10 logarithm of the difference between the amplitudes of consecutive

periods, multiplied by 20.

Shimmer:APQ3 (%) Three-point Amplitude Perturbation Quotient, the average absolute difference between

the amplitude of a period and the average of the amplitudes of its neighbours, divided by

the average amplitude.

Shimmer:APQ5 (%) Five-point Amplitude Perturbation Quotient, the average absolute difference between the

amplitude of a period and the average of the amplitudes of it and its four closest neighbours,

divided by the average amplitude.

Shimmer:APQ11 (%) Eleven-point Amplitude Perturbation Quotient, the average absolute difference between the

amplitude of a period and the average of the amplitudes of it and its ten closest neighbours,

divided by the average amplitude.

Shimmer:DDA (%) Average absolute difference between consecutive differences between the amplitudes of

consecutive period.

NHR (-) Noise-to-Harmonics-Ratio, the amplitude of noise relative to tonal components.

HNR (dB) Harmonics-to-Noise-Ratio, the amplitude of tonal relative to noise components.

RPDE (-) Recurrence Period Density Entropy, the extension of the concept of periodicity.

DFA (-) Detrended Fluctuation Analysis, the stochastic self-similarity of the turbulent noise.

PPE (-) Pitch Period Entropy, the inefficiency of voice frequency control.

Table 1. Overview of measurements used as a features applied to acoustic signals recorded from each subject.

Jan Rusz 79

Only three patients were diagnosed in HY stage of 3 (mild to moderate stage of disease) with some postural instability; the remaining 21 patients were diagnosed in early to mild stage of disease (modified HY 1-2.5). No patient had been under voice therapy. Detailed descriptions of PD patients that participated in this study are summarized in Table 1. As a control group, 22 persons with no history of neurological disorders (15 men and 7 women) were matched for the respective age which ranged from 40 to 91 years ([mean age in years 58.73 (±SD 14.61)]). 2.2. Recording and speech data The speech data was recorded in a quiet room using an external condenser microphone placed at approximately 15 cm from the mouth and coupled to a Panasonic NV-GS 180 video camera. The voice signals were sampled at 48 kHz, with 16-bit resolution. The video material was not used. All subjects were recorded at the time during a single session with a speech pathologist. Each participant was familiarized with vocal tasks and recording procedure. Each participant was instructed to perform on one breath sustained phonation of /i/ at a comfortable pitch and loudness as constant and long as possible, at least 5 seconds. 2.3. Measurement methods Calculations of traditional measures for detecting PD-related dysphonia are performed using the algorithms supplied in the software package Praat [13]. These measures are usually based on application of a short-time autocorrelation for determining the frequency and location of each vibration of the vocal folds (pitch marks) [14]. The jitter and measures of period perturbation are the cycle-to-cycle variability of the period duration of the acoustic signal coming from voice production, derived by taking successive absolute differences between frequencies of each vocal cycle and averaging over a varying number of cycles, optionally normalizing by the overall average. The shimmer and measures of amplitude perturbation are the cycle-to-cycle variability of the maximum extent of the period amplitude of vocal fold vibrations. The average difference of this sequence is taken as a measure of the deviation between cycle amplitudes. The NHR and HNR are derived from the signal-to-noise estimates from the autocorrelation of each cycle. For practical reasons, in this study, only traditional measures that are not affected by individual differences such as age are used include Jitter:local, Jitter:RAP, Jitter:PPQ5, Jitter:DDP, Shimmer:local, Shimmer:APQ3, Shimmer:APQ5, Shimmer:APQ11, Shimmer:DDA, NHR, and HNR. From non-standard measures, recurrence period density entropy addresses the ability of the vocal folds to sustain simple vibration, quantifying the deviations from exact periodicity. It is determined from the normalized entropy Hnorm of the distribution of the signal recurrence periods P(T), representing the uncertainty in the measurement of the exact period in the signal [11]. An increase in RPDE value represents greater hoarseness in the voice. Detrended fluctuation analysis characterizes the extent of turbulent noise in the speech signal, quantifying the stochastic self-similarity of the noise caused by turbulent air-flow in the vocal tract. The DFA algorithm calculates the extent of amplitude variations F(L) of the speech signal over a range of time scales L, and the self-similarity of the speech signal is quantified by the slope α of a straight line on a log-log plot of L versus F(L) [11]. Breathiness or similar dysphonias caused by incomplete vocal fold closure can increase the DFA value. Pitch period entropy measures the impaired control of the stable pitch during sustained phonation [12]. This measure uses a logarithmic pitch scale and subsequently estimates relative logarithmic pitch sequence variation rp using a linear whitening filter. Finally, PPE computes entropy from probability distribution P(rp) of this rp sequence. An increase in the PPE measure better

80 Jan Rusz

reflects the variations over and above natural healthy variations in pitch observed in healthy speech production [12]. See Table 1 for the list of all measurements used in this study. 2.4. Statistics For obtaining statistically significant differences between the groups, we compare the measure of rhythm by using the non-parametric two-sided Wilcoxon rank sum test against the null hypothesis of equal medians, at a significance probability of 0.05.

6.7 6.72 6.74

−0.5

0

0.5

pitch period

signal amplitude

x

t (s)

(a)

PD

1.42 1.425 1.43 1.435

−0.5

0

0.5

t (s)x

HC

pitch period

signal amplitude (b)

0 100 2000

0.05

0.1

0.15

P(T

)

T (samples)

Hnorm = 0.42 (c)

0 100 2000

0.05

0.1

0.15

P(T

)

T (samples)

Hnorm = 0.32 (d)

2 3 40

0.5

1α norm = 0.64

log

F(L

)

log L

(e)

2 3 40

0.5

1α norm = 0.56

log

F(L

)

log L

(f)

−2 0 20

0.5

1

P(r

p)

rp (semitones)

PPE = 0.88 (g)

−2 0 20

0.5

1

P(r

p)

rp (semitones)

PPE = 0.23 (h)

Figure 1. Details of measurement results performed on sustained phonation. The left panel is for a person with PD, the right panel is for the HC subject; (a-b) the sustained vowel phonation signal over five periods zoomed in. The amplitude stability influences the measures of shimmer and the pitch period stability influences the measures of jitter. As can be seen, the healthy signal is more harmonic, which is captured by noise-to-harmonics measures, (c-d) recurrence period density P(T) for recurrence times T, (e-f) log-log plot of scaling window sizes L against fluctuation amplitudes F(L) in detrended fluctuation analysis measure, (g-h) probability densities P(rp) of residual pitch period rp. Pitch period entropy measures the entropy of these probability densities. See main text for detailed description of the algorithm used to calculate these results.

Jan Rusz 81

3. Results Statistically significant differences between the both groups were found in all measures except DFA. Figure 1 demonstrates the example of results of calculating the RPDE, DFA, and PPE values. It also shows the impaired pitch and amplitude stability and more noise addition in signal. The RPDE measure shows higher peaks in the healthy voice, while for PD they are spread over a wide range of values, which can indicate that the vocal folds are not oscillating at regular intervals. The PPE measures the irregular semitone variations in pitch production. The higher peak represents a stable healthy pitch sequence, and it could be picked up by the large entropy value. Table 2 shows the means, standard deviations, differences between PD and HC groups for all measurements in this study.

Subjects Differences

Measurement PD HC between

Mean (SD) Mean (SD) groups

Jitter:local (%) 1.77 1.78 0.62 0.58 P < .001

Jitter:RAP (%) 1.04 1.08 0.34 0.39 P < .001 Jitter:PPQ5 (%) 0.89 0.83 0.31 0.25 P < .001

Jitter:DDP (%) 3.11 3.24 1.02 1.16 P < .001 Shimmer:local (%) 8.05 5.03 3.43 2.58 P < .001

Shimmer:local (dB) 0.76 0.44 0.32 0.23 P < .001 Shimmer:APQ3 (%) 4.07 2.56 1.73 1.43 P < .001

Shimmer:APQ5 (%) 4.56 2.98 1.91 1.44 P < .001 Shimmer:APQ11 (%) 6.41 3.86 2.78 1.78 P < .001

Shimmer:DDA (%) 12.22 7.69 5.19 4.29 P < .001 NHR (-) 0.21 0.28 0.03 0.04 P < .001

HNR (dB) 14.73 6.99 22.37 5.21 P < .001 RPDE (-) 0.34 0.09 0.27 0.07 P < .01

DFA (-) 0.61 0.02 0.61 0.02 P = .086

PPE (-) 0.48 0.28 0.30 0.12 P < .01

Table 2. List of results of all measurements with mean values, standard deviation values, and statistical significances between the measurements.

0 1 2 30

0.02

0.04

0.06

P(J

itter

)

Jitter (%)0 1 2

0

0.05

0.1

P(S

him

mer

)

Shimmer (dB)0 0.5 1

0

0.05

0.1

P(N

HR

)

NHR (−)0 20 40

0

0.01

0.02

0.03

0.04

P(H

NR

)

HNR (dB)

0 0.2 0.4 0.60

0.01

0.02

0.03

P(R

PD

E)

RPDE (−)0.5 0.6 0.70

0.05

0.1

P(D

FA

)

DFA (−)0 0.5 1 1.5

0

0.01

0.02

0.03

0.04

P(P

PE

)

PPE (−)

HC

PD

Figure 2. The probability densities for all measures. The vertical axes are the probability densities P(‘measure’) of the normalized features values of each measure. The dashdot lines are for HC speakers, the solid lines for Parkinsonian’s.

82 Jan Rusz

Figure 2 shows distributions estimated by using the Gaussian kernel density method for all measures; it is retained only one representative measure from all kind of jitter and shimmer feature. From the first view, the measures of shimmer and HNR better separate the PD and HC groups than the other measures.

4. Conclusions

Summarized, significant differences in phonation were found using almost each of the dysphonia measures. This founding are in agreement with other studies that report phonatory impairment as the most salient features of PD speech [6, 7]. The acoustic findings of increased values in jitter, shimmer, and NHR/HNR that may be clinically interpreted as hypophonia, voice hoarseness, and tremolo are in agreement with a previous report on untreated patients with PD [15]. On the other hand, in PD patients treated by dopaminergic drugs, only the jitter values were increased in PD patients compared to controls while shimmer values were similar to those of controls, and NHR/HNR findings were controversial [16, 17]. The results of this paper can ease the clinical monitoring of speech progression and can be also used as feedback in speech treatment. Acknowledgement The study was supported by Grant Agency of the Czech Technical University in Prague – project SGS 10/180/OHK3/2T/13, and Czech Science Foundation - project GACR 102/08/H008. We are obliged to doctors Hana Ruzickova, Evzen Ruzicka, Jan Roth, Jiri Klempir, Veronika Majerova, and Jana Picmausova for provision of clinical data. References [1] O. Hornykiewicz, “Biochemical aspects of Parkinson’s disease,” Neurology Suppl. 2, 51, S2-S9, (1998). [2] M. Hoehn, “The natural history of Parkinson’s disease in the pre-levodopa and post-levodopa eras,” Neurologic Clinics, 10, 331-339, (1992). [3] M. C. de Rijk, C. Tzourio, M. M. Breteler, J. F. Dartiques, L. Amaducci, S. Lopez-Pousa, J. M. Manubens-Bertran, A. Alpérovitch and W. A. Rocca, “Prevalence of parkinsonism and Parkinson's disease in Europe: the EUROPARKINSON Collaborative Study. European Community Concerted Action on the Epidemiology of Parkinson's disease,” J. Neurol. Neurosurg. Psychiatry, 62, 10-15, (1997). [4] S. K. Van Den Eeden, C. M. Tanner, A. L. Bernstein, R. D. Fross, A. Leimpeter, D. A. Bloch, and L. M. Nelson, “Incidence of Parkinson’s disease: Variation by Age, Gender, and Race/Ethnicity,” Am. J. Epidem., 157, 1015-1022, (2003). [5] A. K. Ho, R. Iansek, C. Marigliani, J. Bradshaw, and S. Gates, “Speech impairment in large sample of patients with Parkinson’s disease,” Behav. Neurol., 11, 131-137, (1998).

Jan Rusz 83

[6] J. A. Logemann, H. B. Fisher, B. Boshes, and E. R. Blonsky, “Frequency and coocurrence of vocal tract dysfunction in the speech of a large sample of Parkinson patients,” J. Speech Hear. Disord., 43, 47-57, (1978). [7] C. Ludlow, N. Connor, and C. Bassich, “Speech timing in Parkinson’s and Huntington’s Disease,” Brain. Lang., 32, 195-214, (1987). [8] P. H. Dejonckere, P. Bradley, P. Clemente, G. Cornut, L. Crevier-Buchman, G. Friedrich, P. Van De Heyning, M. Remacle, and V. Woisard, "A basic protocol for functional assessment of voice pathology, especially for investigating the efficacy of (phonosurgical) treatments and evaluating new assessment techniques. Guideline elaborated by the Committee on Phoniatrics of the European Laryngological Society (ELS)," Eur Arch Otorhinolaryngol, vol. 258, pp. 77-82, 2001. [9] A. M. Goberman, C. Coelho, “Acoustic analysis of Parkinsonian speech I: Speech characteristics and L-Dopa therapy,” Neurorehab., 17, 237-246, (2002). [10] J. Rusz, R. Cmejla, H. Ruzickova, E. Ruzicka, “Quantitative acoustic measurements for characterization of speech and voice disorders in early untreated Parkinson’s disease”, J. Acoust. Soc. Am., (2011), in press. [11] M. A. Little, P. E. McSharry, S. J. Roberts, D. A. Costello, and I. M. Moroz, "Exploiting Nonlinear recurrence and Fractal scaling properties for voice disorder detection," Biomedical Engineering Online, vol. 6, pp. -, 2007. [12] M. A. Little, P. E. McSharry, E. J. Hunter, J. Spielman, and L. O. Ramig, “Suitability of dysphonia measurement for telemonitoring of Parkinson’s disease,” IEEE Trans. Biomed. Eng., 56, 1015-1022, (2009). [13] P. Boersma, and D. Weenink, “Praat, a system for doing phonetics by computer,” Glot International, 5, 341-345, (2001). [14] P. Boersma, “Accurate short-term analysis of the fundamental frequency and harmonics-to-noise ratio of a sampled sound,” In Proceedings of the Institute of Phonetics Sciences, (Amsterdam, 17, 97-110, 1993). [15] F. J. Jimenez-Jimenez, J. Gamboa, A. Nieto, J. Guerrero, M. Orti-Pareja, J. A. Molina, E. Garcia-Albea, and I. Cobeta, “Acoustic voice analysis in untreated patients with Parkinson’s disease.” Parkinsonism. Relat. D., 3, 111-116, (1997). [16] J. Gamboa, F. J. Jimenez-Jimenez, A. Nieto, J. Montojo, M. Orti-Pareja, J. A. Molina, E. Garcia-Albea, and I. Cobeta, “Acoustic voice analysis in patients with Parkinson’s disease treated with dopaminergic drugs.” J. Voice, 11, 314-320, (1997) . [17] I. Hertrich, and H. Ackermann, “Gender-specific vocal dysfunction in Parkinson’s disease: electroglottographic and acoustic analysis.” Ann. Otol. Rhinol. Laryngol., 104, 197-202, (1995).

84 Jan Rusz

Detekce náhlých změn v intrakraniálním EEGpomocí vlastních čísel korelační matice

Jan Sova


Abstrakt: Práce ukazuje možnosti využití korelační matice a vlastních číselkorelační matice při detekci počátků epileptických záchvatů a klasifikaci jejichšíření. Epilepsie je zde, v souladu s jedním z aktuálních paradigmat epileptolo-gie, považována za globální děj, do kterého jsou zapojeny i zdánlivě upozaděné,do šíření záchvatů nezapojené, struktury mozku.

1. Úvod

Epilepsie1 je záchvatové onemocnění mozku. Jedná se tedy o onemocnění neurologické a vpopulaci se udává jeho výskyt 1%, ale bude pravděpodobně vyšší. Epilepsie a epileptickésyndromy lze dělit z několika hledisek. Nejčastěji se uvádí dělení dle klinického obrazuzáchvatů. V lékařské veřejnosti se uvádí klasifikace dle oblastí mozku, kde záchvaty vzni-kají (frontální epilepsie, meziotemporopolární epilepsie apod.), jak se šíří a jak velkoučást mozku zabírají (jednoduché parciální záchvaty, komplexní parciální záchvaty, gene-ralizované záchvaty bez křečí, generalizované záchvaty s křečemi) apod. Nebezpečné jsoutzv. farmakorezistentní epilepsie, tj. epilepsie, u nichž se pomocí medikamentů nedaří do-sáhnout uspokojivého stavu pacienta (dlouhá bez záchvatová období či úplné vymizenízáchvatů) a pacient se tak stává kandidátem neurochirurgického výkonu. Cílem mé práceje navrhnout takové analýzy signálů intrakraniálního2 EEG3, které ve výsledku pomohoulékaři při rozhodování o typu a závažnosti onemocnění a volbě následné léčby. Motivacepro tento výzkum je diskutována v kapitole 2.

1.1. Analyzovaná data

Data poskytli Doc. MUDr. Pavel Kršek, Ph.D. a MUDr. Alena Jahodová z Kliniky dětskéneurologie 2.LF UK FN Motol. Jedná se o záznamy intrakraniálního EEG, které byly poří-zeny během týdenního předoperačního vyšetření kandidátů neurochirurgického výkonu sfarmakorezistentní epilepsií. Záznamy obsahují 88 - 122 svodů se vzorkovací frekvencífs = 200 Hz nebo fs = 1000 Hz s iktální, interiktální i klidovou aktivitou diagnostikova-ných pacientů.

1padoucnice, padoucí nemoc, latinsky eufemicky také morbus sacer nebo morbus divinus2vnitrolebečního3Nutno dodat, že intrakraniální EEG je jen jednou z řady diagnostických nástrojů epiletologa.

Jan Sova 85

2. Motivace

Jako motivaci pro svou práci mohu uvést dva problémy, s kterými se setkávají epileptolo-gové v klinické praxi [4] a které tak mohou být i konkretizovanými požadavky na výsledkyanalýzy EEG signálů. První problém můžeme označit jako tzv. pseudotemporální epilep-sie. V případě pseudotemporální epilesie se jedná o možnou záměnu extratemporální epi-lepsie za temporální na základě vyšetření pomocí intrakraniálního (interiktálního a někdyi iktálního) EEG4. Předpokládá se, že k chybné diagnóze dochází proto, že interiktálníEEG výboje v temporálním laloku jsou odrazem šíření z oblasti extratemporální a totošíření nebylo detekováno. Druhým problémem je nedostatečnost paradigmatu epilepto-genního ložiska, pro některé typy a projevy epilepsie. Proto se i kliničtí epileptologovépřiklánějí k paradigmatu „epileptickýchÿ sítí.

2.1. Pseudotemporální epilepsie

„V klinické praxi se opakovaně setkáváme s pacienty, u nichž nás charakter epileptickýchzáchvatů vede k chybné diagnóze epilepsie temporálního5 laloku (TLE). MRI nález u po-stižených sice zpravidla neprokáže strukturální patologii v temporálních lalocích, nicméněnavržené diagnóze většinou odpovídá nález při interiktálním, a někdy i iktálnám EEGvyšetření. Méně charakteristické jsou už výsledky dalších vyšetření, jako neuropsycholo-gie, interiktálního PET a/nebo iktálního SPECT vyšetření. Ve skutečnosti mohou býtinteriktální EEG výboje v temporálním laloku odrazem šíření z oblasti extratemporální.ÿ

2.2. „Epileptickáÿ síť

„Jeden z novějších konceptů zabývajících se patogenezí vzniku epilepsií (včetně TLE)navrhuje představu dynamicky se chovajících „epileptickýchÿ sítí či okruhů, zahrnujícíchvětší množství potenciálních zdrojů iktálních výbojů. Jednotlivé kortikální a subkortikálnístruktury, jež jsou součástí zmíněných okruhů, se navzájem ve své aktivitě významně ovli-vňují, čímž je modifikována citlivost jednotlivých participujících struktur k záchvatovéaktivitě. Podle tohoto konceptu by bylo možné, že záchvatová aktivita může být gene-rována v různých částech rozsáhlého neuronálního okruhu v závislosti na chování jehodalších částí.ÿ

3. Zpracování intrakraniálního EEG signálu

Zvolená metoda analýzy vlastních čísel korelační matice [12] umožňuje detekci synchro-nizačních a desynchronizažních změn, ke kterým dochází před i během epileptického zá-chvatu. Od zvolené metody (v kombinaci s dalšími metodami) si slibuji především: roz-poznávání typu záchvatů (parciální, generalizované, sekundárně generalizované), segmen-taci záchvatu na jednotlivé významné okamžiky (počátek, šíření, generalizace, sekundárnígeneralizace apod.), detekci ložiska a analýzu vlivu jednotlivých částí mozku na šířenízáchvatu.

1. Segmentace: Signál je segmentován na okna velikosti 2 - 5 s ve všech n kanálech,viz obrázek 1. Tím vzniká matice s n segmenty mezi nimiž jsou počítány vzájemnékorelace (jejichž celkový počet je n(n − 1)). Výsledkem je korelační matice R (vizdále).

4částečně charakteristické jsou i výsledky neuropsychologie, PET, SPECT apod.5spánkového

86 Jan Sova

Obrázek 1: Segmentace n kanálového EEG signálu.

2. Výpočet korelační matice: Korelační matice R má strukturu

R =

ρ11 ρ12 · · · ρ1n

ρ21 ρ22 · · · ρ2n...

.... . .

...ρn1 ρn2 · · · ρnn

, (1)

kde n je počet svodů (kanálů) intrakraniálního EEG a ρij je korelace mezi i-tým aj-tým svodem, kterou lze určit vztahem:

ρij =cov(i, j)

σiσj

=E((X − µx)(Y − µy))

σiσj

. (2)

Matice R je symetrická (RT = R) a pro její koeficienty platí

ρij = ρji, (3)

neboť korelace ρij mezi i-tým a j-tým svodem je stejná jako korelace ρji mezi j-týma i-tým svodem. Navíc jsou prvky na diagonále matice jednotkové:

ρij = 1 pro ∀i, j : i = j, (4)

neboť ρii je korelace signálu se sebou samým.

Pozn. U jednotlivých korelačních koeficientů při výpočtu vlastních čísel není uvažo-váno znaménko. Pro zjednodušení bude také dále uvažovano, že maticeR neobsahujenulové prvky.

3. Výpočet vlastních čísel: Nechť A je čtvercová matice řádu n. Pokud existuje čísloλ ∈ C a vektor ~u ∈ Rn, pro které platí

A~u = λ~u, (5)

pak λ se nazývá vlastní číslo (též charakteristické číslo) matice A a ~u vlastní vektornáležící vlastnímu číslu λ. Rovnice (8) se obvykle zapisuje jako homogenní soustavarovnic

Jan Sova 87

(A− λE)~u = ~0, (6)

kde E je jednotková matice řádu n. Tato homogenní soustava rovnic má netriviálnířešení právě tehdy, když je determinant matice soustavy roven nule, tj. pokud platí

|A− λE| = 0. (7)

Determinant A(λ) = |A− λE| nazýváme charakteristický polynom matice A -jedná se o polynom stupně n v proměnné λ, který má v oboru komplexních číseln kořenů. Rovnici A(λ) = 0 nazýváme charakteristická rovnice matice A a jejímikořeny jsou vlastní čísla matice A.

Dále bez důkazu zmíním pro další výklad důležité vlastnosti vlastních čísel a vlast-ních vektorů pro případ kladné reálné symetrické matice (což je právě případ kore-lační matice R).

Věta 1. Je-li A reálná symetrická matice řadu n, potom jsou všechny kořeny cha-rakteristické rovnice reálné, tj. všechna vlastní čísla reálné symetrické matice jsoureálná.

Věta 2. Dva vlastní vektory odpovídající různým vlastním číslům matice A jsounavzájem ortogonální.

Věta 3. Ke každé symetrické matici A existuje ortonormální matice P taková, žePBPT = B je diagonální matice - prkvy bii na diagonále matice B jsou všechnavlastní čísla λi matice A (počítána i s jejich násobností) a sloupcové vektory ma-tice P jsou jednotkové vzájemně ortogonální6 vlastní vektory matice A příslušné kvlastním číslům λi.

Definujeme-li dále stopu matice

tr(A) =n∑

i=1

aii, (8)

pak můžeme s využitím věty 37 ukázat, že pro symetrickou maitici A platí

tr(A) = tr(PTBP) = tr(APPT ) = tr(B). (9)

Protože

tr(R) =n∑

i=1

ρii =n∑

i=1

1 = tr(B) =n∑

i=1

λi, (10)

dostáváme

n∑

i=1

λi = n. (11)

6Pro každou ortogonální matici P platí PTP = E, což využijeme dále.7A skutečnosti, že tr(AB) = tr(BA).

88 Jan Sova

Tedy součet všech vlastních čísel bude roven počtu svodů EEG.

Shrnutí: Po těchto úvahách jsme došli k závěru, že jako výsledek výpočtu vlastníchčísel a vlastních vektorů korelační matice R řádu n tedy vždy dostaneme n reálnýchkladných vlastních čísel, přičemž bude platit

∑ni=1 λi = n, a jim odpovídajících

navzájem ortogonálních n vlastních vektorů.

Výkon je pro jednotlivá vlastní čísla rozdělen v poměru velikosti daného čísla k suměvšech vlastních čísel

P (λi) =λi∑n

j=1 λj

=λi

n. (12)

4. Výsledky

Dále budu diskutovat výsledky aplikace výše popsané metody na signál s epileptickou ak-tivitou, který byl poskytnut z Kliniky dětské neurologie 2. LF UK FN Motol. Při výpočtukorelační matice se ukazuje její zásadní odlišnost pro stav bez záchvatu, pro počátek zá-chvatu, pro šíření záchvatu, a pro stav generalizace. Zatímco pro normální nezáchvatovouEEG aktivitu je korelační matice bez významných shluků lokální synchronizace nebo de-synchronizace (viz obrázek 2a)), objevuji se prakticky bezprostředně po vzniku záchvatumísta se silnou lokální synchronizací i desynchronizaci (viz obrázek 2b)). Během trváníepileptického záchvatu se synchronost šíří do dalších částí mozku (viz obrázek 2c)), ažv konečném stádiu dochází ke globální synchronizaci a ke konvergenci velké části mozkuprakticky k totožnému signálu (výsledný stav generalizace, viz obrázek 2d)). Významnéjsou také změny v hodnotách vlastních čísel. Zatímco v případě nezáchvatového EEG jevýkon rozdělen mezi přibližně 10% vlastních čísel, dochází během vzniku a šíření k přesunuvýznamné části energie mezi zbylých 90% vlastních čísel. Tento přesun energie se projevízmenšením velikosti prvních 10% vlastních čísel a zvětšení velikosti zbylých vlastních čí-sel. Při generalizaci dochází naopak k přesunu 95% veškeré energie na hlavní vlastní číslo.Čím větší část energie je během generalizace záchvatu soustředěna do jednoho vlastníhočísla, tím větší část mozku konverguje k jedné synchronní aktivitě.Statistická významnost změny vlastních čísel korelační matice byla diskutována v [9] navelké skupině pacientů s epileptogenním ložiskem v různých částech mozku.

5. Závěr

Byla prezentována metoda globální analýzy intrakraniálního EEG signálu a bylo diskuto-váno její předpokládané využití v epileptologii (včetně klasifikace záchvatu). Do budoucnabych rád jednak do této globální metody zahrnul lokální vhled (bez toho nebude možnéhledat ložiska záchvatů ani sledovat šíření záchvatů), jednak tuto metodu zkombinoval sdalšími metodami detekce náhlých změn v číslicových signálech.Podstatnou otázkou také je, jaké výsledky bychom dostali při aplikaci metody na čistěparciální záchvaty. Překážkou těmto analýzám je skutečnost, že vyšetření pomocí intra-kraniálního EEG se většinou u pacientů s dobře lokalizovaným epileptogenním ložiskema s negeneralizovanými (parciálními) záchvaty nedělá.

Jan Sova 89

Obrázek 2: a) korelační matice normálního EEG, b) korelační matice těsně na začátkuzáchvatu, c) korelační matice během záchvatu, d) korelační matice již generalizovanéhozáchvatu

Poděkování

Tato práce je podporována granty IGA NT11460-4/2010 Komplexní analýza intrakra-niálního EEG záznamu a identifikace epileptogenní zóny u pacientů s nelezionální far-makorezistentní epilepsií, SGS 10/272/OHK4/3T/13 Analýza intrakraniálních EEGzáznamů výzkumný záměr MSM6840770012 Transdisciplinární výzkum v biomedicín-ském inženýrství.

Reference

[1] Arhuis, M.; Valton, L.; Régis, J. Impaired consciousness during temporal lobe sei-zure is related to increased long-distance cortical-subcortical synchronization. Brain(2009), 2091–2101.

[2] Blume, W. T.; Ociepa, D.; Kander, V. Frontal lobe seizure propagation: Scalp andsubdural eeg studies. Epilepsia (2001), 491–503.

90 Jan Sova

[3] Brabec, J. Vybrané kapitoly z teorie matic. Vydavatelství ČVUT, Praha, 1975.

[4] Brázdil, M.; Marusič, P. Epilepsie temporálního laloku. TRITON, Praha, 2006.

[5] Conlon, T.; Ruskin, H. J.; Crane, M. Seizure characterisation using frequency-dependmultuvariate dynamics. Computers in Biology and Medicine (2009), 760–767.

[6] Krajník, E. Základy maticového počtu. Vydavatelství ČVUT, Praha, 2006.

[7] Netoff, T. I.; Schiff, S. J. Decreased neuronal synchronization during experimentalseizure. The Journal of Neuroscience (2002), 7297–7307.

[8] P., S.; P., P. Vybrané metody číslicového zpracování signálů. Vydavatelství ČVUT,Praha, 2003.

[9] Schnindler, K.; Leung, H.; Elger, C. E.; Lehnertz, K. Assessing seizure dynamics byanalysing the correlation structure of multichannel intracranial eeg. Brain (2007),65–77.

[10] Stam, C. J. Nonlinear dynamical analysis of eeg and meg: Review of an emergingfield. Clinical Neurophysiology (2005), 2266–2301.

[11] Uhlíř, J.; Sovka, P. Číslicové zpracování signálů. Vydavatelství ČVUT, Praha, 2002.

[12] Wendling, F.; Bartolomei, F.; Bellanger, J. J.; Bourien, J.; Chauvel, P. Epileptic fastintracerebral eeg activity: evidence for spatial decorrelation at seizure onset. Brain(2003), 1449–1459.

[13] Zdeněk, V. EEG v epileptologii dospělých. GRADA, Praha, 2004.

Jan Sova 91

Klasifikace mezi /s/ a /š/ na základě parametrizacevstupního signálu pomocí LSF

Adam Stráník


Abstrakt: V běžné řeči jsou sykavky posluchačem velmi často podvědoměhodnoceny. Při tvorbě sykavek může vznikat několik jevů, které negativněovlivňují výsledný vjem. Jedním z nich je hvízdání, tzv. stridence. Jedná se oostrý, jasně ohraničený vysokofrekvenční zvuk, který vzniká nevhodným při-blížením mluvidel při tvorbě /s/. Oproti tomu se často vyskytuje jev, kdy jsousykavky tzv. tupé a mohou přecházet až do hlásky /š/. V příspěvku je před-stavena možnost detekce hvízdání nebo naopak tupých sykavek na základěpopisu signálu pomocí LSF.

1. Úvod

Frikativní hlásky, do kterých sibilanty /s/ a /š/ patří, vznikají těsným přiblížením mlu-videl. Na tomto zúžení dochází k tření proudícího vzduchu o okraje mluvidel a tím vznikácharakteristický šum. Při tvorbě hlásky /s/ zúžení vzniká těsným přitisknutím podélnýchokrajů jazyka k horní dásni. Štěrbina zůstává mezi hřbetem špičky jazyka a přední částídásní. Retní štěrbina bývá poměrně úzká. Podobně u hlásky /š/ je zúžení vytvořeno při-tisknutím okrajů jazyka po stranách k horní dásni, hlavní zúžení je ovšem posunuté vícedozadu, tzn. dále od zubů. Štěrbina zůstává mezi hřbetem přední části jazyka a zadníčásti dásňového výstupku. Hrot jazyka bývá také často skloněn dolů.Hvízdání při syčení, tzv. stridence, je jev, kdy vlivem nevhodného přiblížení mluvidel přitvorbě hlásky /s/ dojde k nárůstu energie na vyšších frekvencích. [1].Line Spectral Pairs (LSP) nebo Line Spectral Frequencies (LSF) je metoda parametrizacesignálu používaná k jednoduchému popisu spektrálních vlastností signálu na základělineární prediktivní analýzy (LPC - Linear Predictive Coding). Koeficienty LSP a lineárníprediktivní analýzy je možné mezi sebou poměrně jednoduše přepočítávat. V této prácibudou ukázány vlastnosti LSP, způsob výpočtu a dále možnost popisu sibilantů pomocíLSF.

2. Princip LSP

Popis hlasového signálu pomocí LSP je založen na válcovém modelu hlasového traktu– trakt je možné modelovat sadou válců stejné délky, ale různého průměru, které jsoudo sebe zasunuté. Pokud bude takovými válci proudit vzduch, bude docházet k různýmrezonancím v závislosti na tom, zda bude tento model na konci otevřený nebo uzavřený

92 Adam Stráník

(otevření a uzavření hlasivek, úst atp.), přičemž počet rezonančních frekvencí závisí napočtu válců, kterými je hlasový trakt modelován, tedy na řádu modelu. Šířku a pozicitěchto rezonančních frekvencí je možné popsat právě pomocí LSF, které přímo souvisí sLSP.LSP je polynom, jehož kořeny jsou LSF. Pokud tyto kořeny seřadíme, dostaneme páryčísel1 (v tomto případě frekvencí), které popisují rezonanční frekvence hlasového traktu.

3. Výpočet a vlastnosti LSP

LSP je odvozeno z koeficientů LPC.

3.1. Výpočet LPC

Vlastnosti LPC vychází z toho, že je možné aproximovat n tý vzorek signálu jako lineárníkombinaci M předchozích vzorků

x[n] =M∑

m=1

amx[n−m], (1)

kde x[n] je odhad n tého vzorku a am je váha daného předchozího vzorku.Výpočet vah vychází z toho, že chceme minimalizovat výkon chyby predikce, která jedána vztahem

e[n] = x[n]− x[n]. (2)

Po derivaci vztahu (2) podle všech vah am a položení derivace rovné nule se dostanemeaž ke vztahu

R(0) R(1) . . . R(M − 1)R(1) R(2) . . . R(M − 2)

. . . . . .R(M − 1) R(M − 2) . . . R(0)

a1

a2

.aM

=

R(1)R(2)

.R(M)

, (3)

kde R je autokorelační funkce a am jsou váhy vzorků. Matici koeficientů am pak jednodušedopočítáme.Vzhledem k tomu, že matice autokorelačních koeficientů je symetrická a pozitivně definitní,je možné k výpočtu použít efektivní rekurzivní Levinson-Durbinův algoritmus — díkyněmu není nutné počítat inverzní matici, získáme odhad chyby pro všechny předchozířády a díky rekurzi získáme koeficienty všech předchozích řádů označované jako PAR-COR (PARtial CORrelation coefficients), viz např. [2].

3.2. Výpočet LSP

Pokud známe LPC koeficienty am, můžeme zavést polynom

A(z) = 1 + a1z + a2z2 + . . . + aMzM , (4)

což je LPC model M tého řádu hlasového traktu pro daný signál.Pokud zavedeme dva polynomy M + 1 řádu P (z) a Q(z), které jsou antisymetrické2 amají k polynomu A(z) následující vztah

A(z) =P (z) + Q(z)

2, (5)

1Sudý a lichý kořen.2Polynom A#(z) je antisymetrický k polynomu A(z), jestliže A#(z) = −A(z).

Adam Stráník 93

získáme LSP polynomy. Pokud nalezneme kořeny těchto polynomů, získáme LSF.Polynomy P (z) a Q(z) lze získat z následujících vztahů:

P (z) = A(z)− z−(M+1)A(z−1), (6)

Q(z) = A(z) + z−(M+1)A(z−1). (7)

3.3. Vlastnosti LSP

Lze dokázat ([3] nebo [4]), že komplexní kořeny Θk LSP polynomů P (z) a Q(z), tedy LSF,leží v z rovině na jednotkové kružnici a jsou navzájem proložené.Samotné frekvence ωk lze z komplexních kořenů Θk polynomů P (z) a Q(z) vypočítatnásledovně

ωk = tan−1

(<Θk

=Θk

). (8)

Dále je dokázáno, že vzhledem k tomu, že polynomy jsou navzájem antisymetrické, odpovídájeden polynom modelu zavřeného hlasového traktu a druhý otevřenému. Pro matematickáodvození a zpětné přepočty mezi LPC a LSP viz [3] nebo [4].Pozice LSF ukazují, kde jsou v signálu rezonanční frekvence, viz obr. 1 – pokud jsou siodpovídající páry frekvenčně blízké, je mezi nimi významnější nárůst energie (např. 1. a 2.LSF nebo 5. a 6. LSF) a naopak, čím jsou od sebe vzdálenější, tím větší lokální minimumenergie mezi nimi je (3. a 4. LSF).

500 1000 1500 2000 2500 3000 3500 400010

−1

100

101

102

→ frequency [Hz]

H(f

) [d

B]

odhad LPCspektra

6. LSF

4. LSF3. LSF

5. LSF

1. LSF 2. LSF

Obrázek 1: Ukázka odhadu LPC spektra 6. řádu pro hlásku /a/ s LSF.

94 Adam Stráník

4. Experimenty

Pro experimenty byly použity nahrávky hlásek /s/ a /š/ jejichž typické spektrum jezobrazené na obr. 2. Nahrávky jsou pořízené s vzorkovací frekvencí 44,1 kHz, 16 bitů navzorek, lineární kvantování. Signál byl segmentován okem délky 20 ms, okna se z 50%překrývala.

Time

Fre

quen

cy

s with stridence

0.5 1 1.5 20

0.5

1

1.5

2

x 104

(a)

Time

Fre

quen

cy

s

1 2 30

0.5

1

1.5

2

x 104

(b)

Time

Fre

quen

cy

sh

0.2 0.4 0.6 0.8 1 1.20

0.5

1

1.5

2

x 104

(c)

Obrázek 2: Typické spektrogramy analyzovaných sibilantů: (a) hláska /s/ s pískáním, (b)hláska /s/, (c) hláska /š/.

0.5 1 1.5 2

x 104

10−1

100

101

→ frequency [Hz]

H(f

) [d

B]

s with stridence

(a)

0.5 1 1.5 2

x 104

100

101

→ frequency [Hz]

H(f

) [d

B]

s

(b)

0.5 1 1.5 2

x 104

10−1

100

101

→ frequency [Hz]

H(f

) [d

B]

sh

(c)

Obrázek 3: Typická LPC spektra a zobrazené LSF v analyzovaných signálech: (a) hláska/s/ s pískáním, (b) hláska /s/, (c) hláska /š/.

Pro parametrizaci byly použity následující míry vypočtené z LSF koeficientů:

4.1. Diference LSF – dLSFa,b

Parametr dLSFa,b je dán vztahem

dLSFa,b = ωi − ωi−1, (9)

kde ω je frekvence vypočtená ze vztahu (8) a i je přirozené kladné sudé číslo, a = i− 1 ab = i.Tento parametr jednoduchým způsobem popisuje, jak je významný nárůst energie vefrekvenční oblasti mezi těmito frekvenčními páry – čím je diference menší, tím je nárůstvýznamnější a naopak. Pokud je diference příliš velká, může se jednat o lokální minimumenergie.

Adam Stráník 95

4.2. Průměr LSF – mLSFa,b

Parametr mLSFa,b je dán vztahem

mLSFa,b =1

2(ωi−1 + ωi) , (10)

kde ω je frekvence vypočtená ze vztahu (8) a i je přirozené kladné sudé číslo, a = i− 1 ab = i.Tento parametr zhruba odpovídá pozici lokálního maxima případně minima energie meziLSP – o minimum se jedná, pokud je diference mezi LSP příliš velká.

5. Výsledky

Na obr. 4 je zobrazen výsledek parametrizace sibilantů /s/ a /š/ pomocí LSF. Z obr. 4(a)a obr. 4(b) je patrné, že existuje významný rozdíl mezi těmito hláskami. Z obr. 4(b) jedále patrné, že je možné oddělit hlásku /s/ s hvízdáním a hlásku /s/ bez hvízdání – hláska/s/ s hvízdáním má jednak menší diferenci 3. a 4. LSF a jednak větší průměr těchto LSF.To odpovídá tomu, že lalok energie mezi těmito LSF je užší a je situován na vyššíchfrekvencích.

2000 3000 4000 5000 6000 7000 80000

1000

2000

3000

4000

5000

LSF means [Hz]

LSF

diff

s [H

z]

1st and 2nd LSF

s with stridence s sh silence

(a)

6000 7000 8000 9000 10000 11000 12000 130000

1000

2000

3000

4000

5000

LSF means [Hz]

LSF

diff

s [H

z]

3rd and 4th LSF


(b)

1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7

x 104

0

1000

2000

3000

4000

5000

6000

LSF means [Hz]

LSF

diff

s [H

z]

5th and 6th LSF


(c)

Obrázek 4: Zobrazení parametrizace hlásek /s/ a /š/: (a) parametrizace z 1. a 2. LSF, (b)parametrizace z 3. a 4. LSF, (c) parametrizace z 5. a 6. LSF.

Na základě parametrizace byl nalezen klasifikátor, který je schopný rozdělovat vstupnísignál do následujících skupin:

96 Adam Stráník

• /s/ s hvízdáním (Sstridence),• /s/ (S),• /š/ (Sh),• ticho (silence).

Při hledání vhodného klasifikátoru byl použit nástroj pro dolování dat WEKA [5]. Vs-tupní data byla rozdělena na 10 skupin se stejným poměrovým zastoupením každé třídya hledané klasifikátory byly testovány pomocí crossvalidace. Výsledný klasifikátor je za-ložený na rozhodovacím stromu, který je zobrazen na obr. 5. V tab. 1 je konfúzní maticetohoto klasifikátoru.

Obrázek 5: Vizualizace rozhodovacího stromu. mean12 – průměr z 1. a 2. LSF; mean56– průměr z 5. a 6. LSF; diff34 – rozdíl 3. a 4. LSF.

Sstridence S Sh silence ← klasifikováno jako267 27 0 0 Sstridence2 384 0 0 S0 1 1525 1 Sh0 0 0 391 silence

Tabulka 1: Konfúzní matice nalezeného klasifikátoru.

6. Závěry

Byl proveden rozbor možností využití LSP (resp. LSF) pro popis sibilantů /s/ a /š/, resp.možnost detekce hvízdání při syčení, tzv. stridence. Hlavní rozdíly mezi hláskou /s/ a /š/při popisu pomocí LSF jsou patrné z diferencí odpovídajících si LFS, kdy hláska /s/ mámenší diferenci než hláska /š/. To je způsobeno rozdílným rozložením energií ve spektruu těchto hlásek. Detekce hvízdání je možná při popisu vstupního signálu LPC modelem6. řádu, ze kterého jsou následně vypočteny diference LSF. Pokud je hvízdání v signálupřítomné, sníží se diference mezi 3. a 4. LSF, viz obr. 4(b).Na základě parametrizace signálu byl nalezen pomocí nástroje WEKA [5] klasifikátor,který je schopný s přesností 98% rozhodnout, zda je vstupní signál /s/ s hvízdáním, /s/,/š/ nebo ticho. Klasifikátor je reprezentován rozhodovacím stromem zobrazeným na obr. 5a jeho konfúzní matice je v tab. 1.

Adam Stráník 97

Poděkování

Tato práce je podporována z: GACR102/08/H008 Analýza a modelování biomedicín-ských a řečových signálů, SGS10/180/OHK3/2T/13 Hodnocení poruch hlasu a řeči,MSM6840770012 Transdisciplinární výzkum v oblasti biomedicínského inženýrství II.

Reference

[1] JOVICIC, S, PUNISIC, S, SARIC, Z. Time-frequency detection of stridence infricatives and africates. In Acoustics. 2008. s. 5137-5141.

[2] PSUTKA, Josef, et al. Mluvíme s počítačem česky. Praha : Academica, 2006. 752 s.

[3] MCLOUGHLIN, Ian Vince. LineSpectral pairs. Signal processing. 2008, 88, s. 448-467. Dostupný také z WWW: <www.sciencedirect.com>.

[4] BÄCKSTRÖM, Tom; MAGI, Carlo. Properties of line spectrum pair polynomi-als: A review. Signal processing. 2006, 88, s. 3286-3298. Dostupný také z WWW:<www.sciencedirect.com>.

[5] HALL M. et. al., The WEKA Data Mining Software: An Update; SIGKDD Explo-rations, 2009, Volume 11, Issue 1.

98 Adam Stráník

Vliv průměrování na možnosti odhadukrevního tlaku z EKG a PPG

Daniel Špulák


Abstrakt: V posledních letech jsou předmětem intenzivního výzkumu nové metody pro neinvazivní kontinuální monitorování krevního tlaku. Oblastí našeho zkoumání jsou zejména postupy založené na využití elektrokardiogramu (EKG) a fotopletysmogramu (PPG). V článku presentujeme dílčí výsledky týkající se vlivu průměrování veličin na možnosti odhadu krevního tlaku.

1. Úvod

Měření krevního tlaku patří mezi běžnou součást lékařské diagnostiky. Ve většině případů se uplatňuje auskultační nebo oscilometrická měřicí metoda. V obou případech je nedílnou součástí aparatury vzduchová manžeta, jíž se zaškrcuje tok krve v přilehlé části těla. Měření je poměrně zdlouhavé a výsledkem je pouze jeden údaj o systolickém a diastolickém tlak.Ve speciálních případech se vyžaduje průběžné monitorování krevního tlaku. Nejpřesnější metodou je invazivní měření, při němž se část aparatury zavádí přímo do krevního oběhu pacienta. Vzhledem k náročnosti této metody a určitým rizikům musí být pro její nasazení zvlášť závažné důvody.Rozšířenou neinvazivní alternativou je metoda odtížení arterie vyvinutá J. Peňázem ([PEN73]), nevýhodu však představuje jednak omezená přesnost, jednak pro pacienta nepříjemné neustálé změny tlaku v zaškrcovací manžetě. Další, zatím převážně experimentální metody, využívají korelaci mezi krevním tlakem a rychlostí šíření pulsní vlny v krevním řečišti.

2. Netradiční neinvazivní měření krevního tlaku

Mezi obvyklé postupy patří snímání EKG a PPG z prstu nebo z ušního lalůčku. Měří se prodleva mezi R-špičkou EKG a stanoveným charakteristickým bodem průběhu PPG (pulse arrival time, PAT) a na jejím základě se odhaduje krevní tlak.V některých článcích se objevují různé modifikace tohoto přístupu. Článek [PAR09] uvádí odkazy na literaturu popisující bezkontaktní snímání EKG i PPG skrz oblečení. Snímače přitom mohou být na židli, v posteli, na klozetu, na počítačové myši apod. Místo PPG může být použit také mechanický pletysmogram. Pramen [KAN06] popisuje využití mechanického snímače tvořeného bimetalovým páskem, jehož indukčnost se mění při ohýbání.Určitý podklad k odhadu průběhu krevního tlaku lze získat i z balistokardiogramu, tedy záznamu mechanických pohybů srdce. Ke snímání se mohou využít modifikované fonokardiografické snímače, různé plošné tlakové snímače přikládané přímo na tělo měřené osoby či na židli nebo postel apod., viz odkazy v [PAR09]. Článek [LIM07] popisuje v souvislosti s monitorováním krevního tlaku měření prodlevy mezi R-špičkou EKG a stahem srdečních komor (pre-ejection period, PEP) za použití EKG a signálu z tlakového snímače na opěradle židle.Mezi nekonvenční přístupy patří i odvozování průběhu krevního tlaku na základě záznamu EKG a fonokardiogramu. [WON06b] rozvíjí úvahy o použití prodlevy mezi R-špičkou EKG a druhou srdeční ozvou a uvádí střední korelaci se systolickým tlakem o velikosti –0,85. V

Daniel Špulák 99

pramenu [ZHA06] je uveden teoretický rozbor závislosti spektrálního složení druhé srdeční ozvy na systolickém tlaku.

3. Experimentální měření krevního tlaku s použitím EKG a PPG

3. 1 Databáze pacientů a naměřených signálů

Od prvních měření prováděných v rámci diplomové práce ([SPU09]) jsme získali několik dalších náměrů. Zatímco první měření probíhala ve Fakultní nemocnici Motol (FNM), novější soubory pochází z nemocnice Na Homolce (NH), jak shrnuje tabulka 1.

V tabulce 1 jsou uvedeny pouze kompletní, pro další zpracování použitelné záznamy. Poslední sloupec přibližuje rozsah změřeného krevního tlaku (malý do 15 mmHg, střední do 30 mmHg, velký nad 30 mmHg): obecně jsou pro naše účely vhodné záznamy s co nejvyšším rozsahem naměřených hodnot.

3.2 Vliv průměrování

Některé prameny se zmiňují o hysterezním charakteru závislosti mezi PAT a krevním tlakem ([KAN06]). Uvádí se, že zlepšení přesnosti odhadu krevního tlaku je možné dosáhnout průměrováním veličin přes několik srdečních cyklů (resp. R-R intervalů).Při našem zkoumání jsme volili průměrování přes 3, 6, 11 a 16 srdečních cyklů a srovnávali korelace mezi PAT a systolickým nebo diastolickým tlakem. PAT byl u všech subjektů definován jako čas mezi R-špičkou EKG a následujícím minimem PPG. Výsledky shrnují tabulky 2 a 3 a grafy 1 a 2.

Tabulka 1: Soupis náměrů; FNM – Fakultní nemocnice Motol, NH – nemocnice Na Homolce

Datum Místo Označení Rozsah tlaku

18.3.2009 FNM 2 180 velký23.4.2009 FNM 4 224 velký23.4.2009 FNM 5 211 střední23.4.2009 FNM 6 210 malý16.3.2010 NH 16 306 malý16.3.2010 NH 17 319 střední13.4.2010 NH 25 248 střední13.4.2010 NH 26 304 malý13.4.2010 NH 27 257 střední13.4.2010 NH 28 345 malý

Délka záznamu/s

Tabulka 2: Absolutní hodnoty korelačních koeficientů mezi systolickým tlakem a PAT

Subjekt

2 16 17 25 26 27 28

1 0,25 0,47 0,11 0,14 0,56 0,00 0,033 0,15 0,54 0,25 0,56 0,72 0,14 0,066 0,16 0,51 0,27 0,53 0,76 0,27 0,1611 0,13 0,30 0,29 0,52 0,80 0,37 0,2616 0,16 0,10 0,33 0,59 0,86 0,43 0,32

Korelační koeficienty (abs. h.): systolický tlak vs. čas do minima PPG

Počet R-R intervalů pro průměrování

100 Daniel Špulák

Tabulka 3: Absolutní hodnoty korelačních koeficientů mezi diastolickým tlakem a PAT

Subjekt

2 16 17 25 26 27 28

1 0,20 0,45 0,08 0,29 0,53 0,28 0,003 0,01 0,56 0,16 0,53 0,69 0,07 0,026 0,13 0,57 0,15 0,53 0,72 0,12 0,1311 0,34 0,45 0,14 0,51 0,76 0,17 0,2316 0,46 0,35 0,17 0,57 0,81 0,22 0,28

Korelační koeficienty (abs. h.): diastolický tlak vs. čas do minima PPG


Graf 1: Závislost korelačních koeficientů na průměrování pro jednotlivé subjekty; systolický tlak

1 3 6 11 160,00

0,20

0,40

0,60

0,80

1,00

Absolutní hodnota korelačního koeficientuSystolický tlak vs. čas do minima PPG

2161725262728


Kor

el. k

oef.

(abs

. hod

n.)

Graf 2: Závislost korelačních koeficientů na průměrování pro jednotlivé subjekty; diastolický tlak

1 3 6 11 160

0,2

0,4

0,6

0,8

1

Absolutní hodnota korelačního koeficientuDiastolický tlak vs. čas do minima PPG

2161725262728


Kor

el. k

oef.

(abs

. hod

n.)

Daniel Špulák 101

Korelace se obecně zvyšuje s rostoucím počtem srdečních cyklů uplatněných pro průměrování. Výjimku tvoří pouze subjekt 16 (u systolického i diastolického tlaku) a částečně i subjekt 2 (u systolického tlaku). Absolutní hodnoty korelačních koeficientů mezi PAT a systolickým tlakem se pohybují od 0,00 do 0,56 bez průměrování až po 0,10 až 0,86 při průměrování přes 16 srdečních cyklů. V případě diastolického tlaku jsou hodnoty obecně nižší, od 0,00 do 0,53 bez průměrování až po 0,17 až 0,81 při průměrování přes 16 R-R intervalů.

4. Diskuse a závěr

U většiny subjektů se absolutní hodnota korelačního koeficientu zvyšovala s počtem R-R intervalů použitých pro průměrování. Je však zřejmé, že s rostoucí délkou průměrovacího okna klesá časové rozlišení odhadu tlaku a tedy i schopnost zachytit rychlé změny krevního tlaku. Z tohoto hlediska se jeví jako rozumné maximum přibližně 16 srdečních cyklů, což odpovídá obvykle zhruba 10 až 15 sekundám.Celkově jsou námi zjištěné korelace mezi tlakem a dobou šíření pulsní vlny nižší, než uvádějí některé prameny: například článek [WON10] udává u systolického tlaku absolutní hodnotu kolem 0,81, ve [WON06b] je publikována dokonce průměrná absolutní hodnota 0,95. Domníváme se, že příčinou odlišných výsledků je zejména odlišný zdravotní stav měřených osob, neboť ve zmíněných studiích se měřilo převážně na zdravých osobách nízkého až středního věku, zatímco námi prováděná měření probíhala na starších pacientech s různými chorobami oběhové soustavy.

5. Poděkování

Tento výzkum je podporován výzkumným záměrem MSM6840770012 „Transdisciplinární výzkum v oblasti biomedicínského inženýrství II“ a granty GAČR 102/08/H008 „Analýza a modelování biomedicínských a řečových signálu” a SGS10/273/OHK3/3T/13 „Analýza signálů mechanické aktivity srdce“.

Reference

[DEB07] S. Deb, C. Nanda, et al., Cuff-less Estimation of Blood Pressure using Pulse Transit Time and Pre-ejection Period, 2007 International Conference on Convergence Information Technology, 2007

[KAN06] E. Kaniusas, H. Pfützner, et al., Method for Continuous Nondisturbing Monitoring of Blood Pressure by Magnetoelastic Skin Curvature Sensor and ECG, IEEE Sensors Journal, Vol. 6., No. 3, 2006

[LIM07] Y. G. Lim, K. H. Hong, et al., Mechanocardiogram Measured at the Back of Subject Sitting in a Chair as a Non-intrusive Pre-ejection Period Measurement, 2007

[PEN73] J. Peňáz, Photoelectric measurement of blood pressure, volume and flow in the finger, Digest of the 10th International Conference on Medical, and Biological Engineering, International Federation for Medical and Biological Engineering, s. 904, 1973.

102 Daniel Špulák

[SPU09] D. Špulák, R. Čmejla, Monitorování systolického tlaku z údajů EKG a PPG, diplomová práce, ČVUT FEL, 2009

[PAR09] K. S. Park, Nonintrusive Measurement of Biological Signals for Ubiquitous Healthcare, 31st Annual International Conference of the IEEE EMBS, 9/2009

[WON06b] M. Y. M. Wong, C. C. Y. Poon, et al., Can the Timing-Characteristics of Phonocardiographic Signal be Used for Cuffless Systolic Blood Pressure Estimation?, Proceedings of the 28th IEEE EMBS Annual International Conference, 9/2006

[WON10] M. Y. M. Wong, E. Pickwell-MacPherson, et al., The effects of pre-ejection period on post-exercise systolic blood pressure estimation using the pulse arrival time technique, Eur J Appl Physiol, 8/2010

[ZHA06] X. Y. Zhang, Y. T. Zhang, Model-based Analysis of Effects of Systolic Blood Pressure on Frequency Characteristics of the Second Heart Sound, Proceedings of the 28th IEEE EMBS Annual International Conference, 9/2006

Daniel Špulák 103

Zolotarevova transformace a spektrální analýza

Václav Turoň, Radim Špetík, Pavel Sovka, Miroslav Vlček

České vysoké učení v Praze, Fakulta elektrotechnická [email protected], [email protected], [email protected], vlcek [email protected]

Abstrakt: Tento článek se zabývá využitím nové Zolotarevovy transformace (ZT) při krátkodobé spektrální analýze nestacionárních signálů. ZT je frekvenčně-časová transformace, která využívá adaptivní bázi složenou ze Zolotarevových polynomů.

1. Úvod

Při zpracování signálů se velmi často analyzují nestacionární signály, jako jsou například biologické signály, řeč nebo různé diagnostické signály z mechanických strojů. K tomuto účelu se využívá mnoho typů transformací - Krátkodobá Fourirova transformace (Short Time Fourier Transform - STFT), vlnková transformace (Wavelet Transform - WT), Hilbert-Huangova transformace (HT) nebo analýza hlavních komponent (Principal Spectral Component - PCA). Zolotarevova transformace (ZT) je nová frekvenčně-časová transformace využívající adaptivní bázi měnící s mírou nestacionarity analyzovaného signálu. 2. Teorie

Na první pohled se může zdát, že ZT je velmi podobná Fourierově transformaci (FT), protože stejně jako FT hledá maximální korelaci mezi vstupním signálem a danou bází transformace. Mezi FT a ZT transformací je však jeden velký rozdíl. ZT využívá adaptivní bázi, která se mění s mírou nestacionarity vstupního signálů. Tato báze je složena ze dvou selektivních Zolotarevových polynomů prvního a druhého druhu

( )

( ) ( ) ,,2sin2cos

)2zsin()2zcos(2zexp

02

02 Zlmtbimta

ltiltltil

mm

l

mm ∈′+′=

+=

∑∑==

ππ

πππ (1)

kde ma2′ , mb2′ jsou koeficienty Zolotarevových polynomů a l je jejich řád. Báze FT může být

zapsána pomocí komplexní exponenciály jako ( ) Zlltiltlti ∈+= ),2sin()2cos(2zexp πππ , (2)

kde l představuje rovněž řád dané báze. Jak je vidět z porovnání obou bází FT a ZT (obr. 1), je možné vyjádřit bázi FT pomocí báze ZT. Výška centrálního laloku zcos a dvou centrálních laloků zsin se nastavuje eliptickým modulem k , který je úměrný míře nestacionarity vstupního signálu. Jestliže je vstupní signál stacionární, eliptický modul k je roven nule a báze ZT odpovídá bázi FT.

104 Václav Turoň

Obrázek 1: a) Trigonometrická funkce cos b) Zolotarevův polynom zcos c) Trigonometrická

funkce sin d) Zolotarevův polynom zsin 3. Aplikace 3.1 Krátkodobá spektrální analýza nestacionárního signálu Harmonický signál se skokovou změnou frekvence je vygenerován podle následujícího vztahu

( ) 2000,10,2

cos =−=

= NNn

f

nfns

S

K

π (3)

Frekvence signálu f je změněna z 600Hz na 3100Hz podle podmínky

−+∈

∈=

1,12/,3100

2/,0,600

NNn

Nnf . (4)

Na obrázku 2 je zobrazen vstupní signál společně s krátkodobou Fourierovou transformací (tzv. spektrogram) a krátkodobou Zolotarevovou transformací (tzv. zologram), kterou lze získat zaměněním báze FT za bázi ZT při výpočtu krátkodobé Fourierovy transformace. Spektrogram je vytvořen pomocí Hammingova okna a Zologram pomocí obdélníkového okna. Krok segmentace je roven jedné a je pro obě transformace stejný. Nestacionarita vstupního signálu je vytvořena skokovou změnou frekvence v polovině délky signálu (obr. 2a). Spektrogram s délkou okna segmentace 512 vzorků tuto nestacionaritu zachytí, ale její časová lokalizace je dosti nepřesná (obr. 2b). Použitím kratšího okna (128 vzorků) u Fourierovy transformace získáme lepší časové rozlišení, ale zároveň se tím zhorší frekvenční

Václav Turoň 105

rozlišení (obr. 2c). Na druhou stranu zologram tuto skokovou změnu frekvence dokáže velmi přesně časově lokalizovat a to s přesností na 3 vzorky na rozdíl od uvedených spektrogramů (obr. 2d). Délka okna pro vytvoření zologramu je rovna 512 a jak je z obrázku (obr. 2d) vidět, zologram má zachované jak dobré frekvenční, tak i časové rozlišení.

Obrázek 2: a) Harmonický signál se skokovou změnou frekvence generovaný podle

vztahu (3), b) Logaritmické zobrazení spektrogramu vytvořeného Hammingovým oknem o délce 512 vzorků, c) Logaritmické zobrazení spektrogramu vytvořeného Hammingovým

oknem o délce 128 vzorků, d) Logaritmické zobrazení zologramu vytvořeného obdélníkovým oknem o délce 512 vzorků

3.2 Krátkodobá spektrální analýza řeči Obrázek 3 zobrazuje reálný signál (slovo „šest“) společně s jeho spektrogramem a zologramem. Spektrogram je vytvořen pomocí Hammingova okna o délce 256 vzorků. Tato délka je zvolena jako kompromis mezi dobrým frekvenčním a časovým rozlišením (obr. 3b). Zologram je vytvořen pomocí obdélníkového okna o délce 512 vzorků (obr. 3c). Krok segmentace je v obou případech roven jednomu vzorku. Jak je na zobrazeném zologramu vidět, krátkodobá Zolotarevova transformace dokáže zachytit pulsování mezi jednotlivými frekvenčními složkami řeči na rozdíl od spektrogramu, který tuto informaci neposkytuje.

106 Václav Turoň

Obrázek 3: a) Reálný signál – slovo „šest“, b) Spektrogram vytvořený Hammingovým oknem

o délce 256 vzorků, c) Zologram vytvořený obdélníkovým oknem o délce 512 vzorků

V logaritmickém zobrazení (obr. 4) výše uvedeného zologramu lze vidět, že zologram dokáže zaznamenat glotální rázy, které vznikají při tvorbě hlasu. Glotální rázy jsou zobrazeny jako spektrum několika dirakových impulsů (obr. 4c – oblast mezi segmenty 1600÷2100 a 2600÷3200). Spektrogram tyto glotální rázy zaznamenat nedokáže, protože časové rozlišení FT se zvyšuje zkracováním délky okna segmentace. U ZT se časové rozlišení zvyšuje zvyšováním samotných bázových polynomů.

Obrázek 4: a) Reálný signál – slovo „šest“, b) Logaritmické zobrazení spektrogramu vytvořeného Hammingovým oknem o délce 256 vzorků, c) Logaritmické zobrazení

hologramu vytvořeného obdélníkovým oknem o délce 512 vzorků

Václav Turoň 107

4. Závěr Tento článek představuje novou časově-frekvenční Zolotarevovu transformaci společně s jejím možným využitím při krátkodobé spektrální analýze. Jak je z výše uvedených výsledků patrné, Zolotarevova transformace se zdá být dobrým nástrojem pro analýzu nestacionárních signálu. 5. Poděkování

Tento výzkum byl podporován z grantu GAČR 102/08/H008 “Analýza a modelování biomedicínských a řečových signálu”, SGS10/181/OHK3/2T13 „Spektrální vlastnosti Zolotarevovy transformace“ a výzkumného záměru MŠMT MSM6840770012 “Transdisciplinární výzkum v biomedicínckém inženýrství 2”. Reference [1] Špetík, R. The Discrete Zolotarev Transform, Czech Technical University in Prague,

Faculty of Electrical Engineering, Prague, 2009 [2] M. Vlček and R. Unbehauen. Zolotarev Polynomials and Optimal FIR Filters, IEEE

Transaction on Signal processing, vol. 47, no. 3, March 1999, pp. 717 – 730, ISSN 1053-587X

[3] P. Sovka, P. Pollák, Vybrané metody číslicového zpracování signálů, České vysoké

učení technické v Praze, Fakulta elektrotechnická, Praha, 2003

108 Václav Turoň

Date post:	12-Oct-2020
Category:	Documents
Upload:	others
View:	0 times
Download:	0 times

ANALÝZA A ZPRACOV`N˝ ØE¨OVÝCH A BIOLOGICKÝCH SIGN`LÙhynek/citing_papers/... · 2011. 1....

Documents