ANALÝZA A ZPRACOVÁNÍ SBORNÍK PRACÍ 2010

Č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íku Doc. 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] [email protected]

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

Foniatrická klinika 1.LF UK a VFN http://fonja.lf1.cuni.cz vedoucí: 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í biomedicí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ě zpracová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ím na ř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 architektur ví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 bayesovské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 korelační matice 85 Adam Stráník: Klasifikace mezi /s/ a /š/ na základě parametrizace vstupního signá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

Jan Bartošek

1

Porovn´ av´ an´ı algoritm˚ u pro detekci z´ akladn´ı frekvence se zamˇ eˇ ren´ım na ˇ reˇ cov´ e sign´ aly Jan Bartoˇsek ˇ e vysok´e uˇcen´ı technick´e v Praze, Fakulta elektrotechnick´a Cesk´ [email protected] Abstrakt: Pˇr´ıspˇevek volnˇe navazuje na [1] a zab´ yv´a se pˇredevˇs´ım objektivn´ım porovn´an´ım implementovan´ ych algoritm˚ u pro detekci z´akladn´ı frekvence sign´alu (F0). Pro tento u ´ˇcel byl navrˇzen a realizov´an hodnot´ıc´ı framework vyuˇz´ıvaj´ıc´ı pro porovn´an´ı referenˇcn´ı datab´azi a byla ustavena sada hodnot´ıc´ıch krit´eri´ı. Z v´ ysledk˚ u vypl´ yv´a hlavn´ı nedostatek vˇetˇsiny algoritm˚ u v n´ızk´e u ´spˇeˇsnosti detekce znˇel´ ych a neznˇel´ ych u ´sek˚ u promluvy. Z´aroveˇ n je diskutov´ana ot´azka optim´aln´ıho ˇcasov´eho rozliˇsen´ı algoritm˚ u.

1.

´ Uvod

Intonace (zmˇena z´akladn´ı frekvence v ˇcase) lidsk´eho hlasu je jedn´ım z hlavn´ıch prozodick´ ych rys˚ u naˇs´ı ˇreˇci. Extrakce pr˚ ubˇehu intonace m˚ uˇze m´ıt ve zpracov´an´ı ˇreˇci nemal´ y ´ v´ yznam [1]. Uloha nen´ı bohuˇzel zcela snadn´a, protoˇze na rozd´ıl od zpˇevu ˇci prodlouˇzen´e fonace ne vˇsechny u ´seky promluvy maj´ı vlastnost znˇelosti (tj. jsou tvoˇreny s vyuˇzit´ım hlasivkov´ ych pulz˚ u). Tud´ıˇz je tˇreba, aby algoritmus detekuj´ıc´ı z´akladn´ı frekvenci (Pitch Detection Algorithm, PDA) dˇelal toto s maxim´aln´ı moˇznou pˇresnost´ı, ale z´aroveˇ n spr´avnˇe rozliˇsil znˇel´e a neznˇel´e u ´seky. Objektivn´ıho vz´ajemn´eho porovn´an´ı PDA mezi sebou lze dos´ahnout pouˇzit´ım F0 referenˇcn´ı datab´aze s vhodnou sadou krit´eri´ı. Pr´ace pˇredstavuje n´avrh a realizaci takov´eho hodnot´ıc´ıho prostˇred´ı. D´ale popisuje v´ıce do hloubky nˇekter´e nov´e PDA metody (zamˇeˇruje se zejm´ena na MNBFC) a na z´akladˇe v´ ysledk˚ u se zab´ yv´a tak´e realizac´ı jednoduch´eho detektoru znˇel´ ych a neznˇel´ ych u ´sek˚ u ˇreˇci. Kromˇe shrnut´ı v´ ysledk˚ u je posledn´ı ˇca´st vˇenov´ana zhodnocen´ı poˇzadavk˚ u na pl´anovan´ y interpunkˇcn´ıho detektoru.

2. 2.1.

Testovan´ e PDA algoritmy Standardn´ı implementovan´ e algoritmy

Byla implementov´ana vˇetˇsina PDA metod teoreticky popsan´ ych v [1], jmenovitˇe kromˇe jiˇz dˇr´ıve implementovan´e autokorelaˇcn´ı funkce ve frekvenˇcn´ı oblasti (ACF freq) tak´e autokorelace v ˇcasov´e oblasti (ACF time), Average Magnitude Difference Function (AMDF) a kepstr´aln´ı metoda (Ceps). Pro pˇripomenut´ı jsou rovnicemi (1), (2), (3) a (4) vyj´adˇreny vztahy pro jejich v´ ypoˇcet.

2

Jan Bartošek

1 ACFtime (τ ) = N AM DF (τ ) =

2.2.

N −n−1 X

x(n)x(n + τ )

(1)

n=0

−1 1 NX |x(n) − x(n + τ )| N n=0

(2)

ACFf req (n) = IF F T {[abs(F F T (x(k)))]2 }

(3)

Ceps(n) = IF F T {log(abs(F F T (x(k))))}

(4)

Real-time time domain pitch tracking using wavelets

Metoda je detailnˇe pops´ana v [5] a byla teoreticky opˇet pˇredstavena jiˇz v [1]. Jej´ı uv´adˇen´e v´ ysledky vypadaly velmi slibnˇe. Bˇehem implementace vˇsak bylo zjiˇstˇeno, ˇze avizovan´a nˇekolika´ urovˇ nov´a vlnkov´a transformace je v tomto pˇr´ıpadˇe ve v´ ysledku v kaˇzd´e u ´rovni zploˇstˇena na pouhou doln´ı propust sign´alu s n´asledn´ ym podvzorkov´an´ım (5) a testem kandid´ata na F0 (peak-picking a hled´an´ı centr´aln´ı m´odu ˇcasov´ ych diferenc´ı). Pokud nen´ı nalezen, takto transformovan´ y sign´al postupuje do dalˇs´ı u ´rovnˇe. Souˇc´ast´ı algoritmu by mˇel b´ yt detektor znˇel´ ych/neznˇel´ ych u ´sek˚ u na z´akladˇe pomˇeru energi´ı poˇc´ıtan´ ych z tˇretin aktu´aln´ıho okna. Takov´ y detektor se vˇsak sv´ ymi v´ ysledky ani nepˇribliˇzoval realitˇe a proto byl tento PDA vˇedomˇe testov´an bez samostatn´e f´aze detekce Voiced/Unvoiced (avˇsak algoritmus samotn´ y m´a tak´e schopnost ohodnotit u ´sek v krajn´ıch pˇr´ıpadech jako neznˇel´ y). a(n) = [x(2 ∗ n − 1) + x(2 ∗ n)]/2 2.3.

(5)

Merged Normalized Forward-Backward Correlation (MNFBC)

Tato metoda zpracov´an´ı sign´al˚ u pracuj´ıc´ı v ˇcasov´e oblasti byla pops´ana v [3] jako z´aklad velmi komplexn´ıho PDA. Jej´ım j´adrem je poˇc´ıt´an´ı dvou proti sobˇe jdouc´ıch (auto)korelac´ı. Rovnice (7)/(8) ukazuje v´ ypoˇcet dopˇredn´e/zpˇetn´e normalizovan´e korelace, kde konstanta MAX PER odpov´ıd´a ˇcasov´e periodˇe nejniˇzˇs´ı detekovan´e frekvence, vˇzdy se poˇc´ıtaj´ı z ubˇehy obou funkc´ı pro referenˇcn´ı znˇelou ˇca´st promluvy jsou okna d´elky 4*MAX PER. Pr˚ vykresleny na obr´azku 1a. Obˇe funkce jsou n´aslednˇe jednocestnˇe usmˇernˇeny (NFC’ a NFC’) a pouˇzity pro v´ ypoˇcet normalizovan´e spojen´e“ korelace MNFBC (9), jej´ı pr˚ ubˇeh ” je na obr´azku 1b. Rovnice (6) ukazuje form´aln´ı pˇredpis pro pouˇzit´ y korelaˇcn´ı ˇclen. < xwk [n], xwl [n] > =

2∗M AX XP ER−1

xw [n + k]xw [n + l]

(6)

n=0

N F C[t] = p

< xw0 [n], xwt [n] > < xw0 [n], xw0 [n] >< xwt [n], xwt [n] >

< xw2M AX

[n], xw2M AX

(7)

[n] >

N BC[t] = q

< xw2M AX

M N F BC[t] =

< xw0 [n]], xw0 [n]] > (N F C 0 [t])2 + < xw2M AX P ER [n], xw2M AX P ER [n] > (N BC 0 [t])2 < xw0 [n], xw0 [n] > + < xw2M AX P ER [n], xw2M AX P ER [n] >

P ER

[n], xw2M AX

P ER

P ER

P ER−t

[n] >< xw2M AX

P ER−t

[n], xw2M AX

P ER−t

[n] >

(8)

(9)

Jan Bartošek

3

Funkce NFC(t) a NBC(t)

MNFBC(t)

1

1 NFC(t) NBC(t)

0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8

MNFBC(t)

0.9

Hodnota normalizovan´e korelace

Hodnota normalizovan´e korelace

0.8

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0

50

100

150

ˇ Casov´ y interval [vzork˚ u]

200

250

(a) Funkce NFC a NBC

0

0

50

100

150

ˇ Casov´ y interval [vzork˚ u]

200

(b) Funkce MNFBC

Obr´azek 1: Pr˚ ubˇehy korelaˇcn´ıch funkc´ı na znˇel´em referenˇcn´ım u ´seku promluvy 2.4.

Direct Frequency Estimation (DFE)

Metoda pracuje ˇcistˇe v ˇcasov´e oblasti. Je detailnˇe pops´ana v [2] a algoritmus byl pˇrevzat v bin´arn´ı podobˇe. Obsahuje f´azi pro rozhodnut´ı o znˇelosti dan´eho u ´seku. Je doposud pouˇz´ıv´ana v mnoha vˇedeck´ ych prac´ıch zab´ yvaj´ıc´ıch se anal´ yzou hlasu na naˇsem pracoviˇsti.

3.

Ot´ azka optim´ aln´ıho ˇ casov´ eho rozliˇ sen´ı PDA

Pro zodpovˇezen´ı t´eto ot´azky je tˇreba uvˇedomit si nˇekter´a biologick´a fakta o hlasov´em ˇ u ´stroj´ı. Casov´ e rozliˇsen´ı algoritmu ud´av´a, jak ˇcasto je poˇc´ıt´an nov´ y odhad F0 a je d´ano zejm´ena d´elkou pˇrekryvu analyzovan´ ych oken sign´alu. Na jedn´e stranˇe je naˇs´ım c´ılem z´ıskat co nejdetailnˇejˇs´ı informaci o pr˚ ubˇehu F0, na stranˇe druh´e jsou zde pr´avˇe biofyzik´aln´ı limitace traktu, d´ıky kter´ ym vyˇsˇs´ı ˇcasov´e rozliˇsen´ı od urˇcit´e hranice postr´ad´a smysl, protoˇze nepˇrin´aˇs´ı ˇz´adnou novou informaci a vede jen k n´ar˚ ustu v´ ypoˇcetn´ı n´aroˇcnosti, kter´a hraje kl´ıˇcovou roli zejm´ena pˇri nasazen´ı aplikac´ı pracuj´ıc´ıch v re´aln´em ˇcase. D´ale v textu zm´ınˇen´a F0 referenˇcn´ı datab´aze pracuje s ˇcasov´ ym rozliˇsen´ım 1ms, avˇsak je k rozhodnut´ı, zda-li jiˇz nen´ı takto vysok´e rozliˇsen´ı nad limitem hlasov´eho u ´stroj´ı. Odpovˇed’ lze nal´ezt napˇr´ıklad ve studii [6] zab´ yvaj´ıc´ı se rychlost´ı zmˇeny v´ yˇsky hlasivkov´eho t´onu a podle kter´e nejvˇetˇs´ı detekovan´a mluven´a intonaˇcn´ı rychlost v d´anˇstinˇe byla 50 p˚ ult´on˚ u za sekundu (50 cent˚ u p˚ ult´onu za 10ms). Lze pˇredpokl´adat, ˇze tato rychlost nen´ı v´ yraznˇeji z´avisl´a na jazyku a pro ˇceˇstinu bude dosahovat podobn´e hodnoty. Dodejme jen, ˇze se jedn´a o limitn´ı hodnoty, kter´ ych se zˇr´ıdkakdy v re´aln´ ych promluv´ach dosahuje. V´ ysledky studie potvrzuj´ı i n´asleduj´ıc´ı ilustrace se dvˇema uk´azkov´ ymi pr˚ ubˇehy F0 na celoznˇel´em u ´seku promluvy detekovan´ ymi testovan´ ym algoritmem ACF freq. Obr´azek 2a zn´azorˇ nuje jeden z moˇzn´ ych velmi rychl´ ych intonaˇcn´ıch pr˚ ubˇeh˚ u (melod´em˚ u) t´azac´ı vˇety s ˇcasov´ ym rozliˇsen´ım 23ms (FS=11025Hz, posun oken 256 vzork˚ u). V u ´seku nejvˇetˇs´ıho vzestupu je vidˇet, ˇze takov´e ˇcasov´e rozliˇsen´ı ne zcela dostaˇcuje. Proto bylo v dalˇs´ı pr´aci pouˇz´ıv´ano rozliˇsen´ı 16ms, jehoˇz aplikaci na rozechvˇen´ y zpˇevov´ y sign´al (vibrato) lze vidˇet na obr´azku 2b. Tato hodnota se zd´a b´ yt dostaˇcuj´ıc´ı a byla proto zvolena jako v´ ychoz´ı i pro testov´an´ı vˇsech F0 algoritm˚ u.

250

4

Jan Bartošek

(a) Rozliˇsen´ı F0 23ms, velmi rychl´a ot´ azka

(b) Rozliˇsen´ı F0 16ms, vibrato zpˇev

Obr´azek 2: Vliv ˇcasov´eho rozliˇsen´ı F0 na detekovan´ y pr˚ ubˇeh intonace

4.

F0 referenˇ cn´ı datab´ aze

K tomu, abychom dok´azali objektivnˇe ohodnotit PDA, potˇrebujeme zn´at spr´avn´e v´ ysledky pro vstupn´ı data. Proto byla v pr´aci pouˇzita ruˇcnˇe znaˇckovan´a ˇc´ast Speecon Spanish Database, jej´ıˇz referenˇcn´ı ˇc´ast byla vytvoˇrena jako v´ ysledek pitch-marking1 algoritmu [4] a d´ale ruˇcnˇe pˇrekontrolov´ana. Referenˇcn´ı F0 jsou potom z´ısk´any jako pˇrevr´acen´e hodnoty ˇcasov´ ych vzd´alenost´ı sousedn´ıch pitch-mark˚ u. Datab´aze byla nahr´ana s n´asleduj´ıc´ımi parametry: RAW audio form´at se vzorkovac´ı frekvenc´ı 16kHz, bitov´a hloubka 2B/vzorek, line´arn´ı k´odov´an´ı, mono. V nahr´avk´ach je pˇr´ıtomno 60 mluvˇc´ıch (30 ˇzen a 30 muˇz˚ u). Celkov´a doba z´aznamu pˇresahuje jednu hodinu, coˇz pˇredstavuje v pr˚ umˇeru v´ıce neˇz 1 minutu na mluvˇc´ıho. Datab´aze byla nahr´av´ana souˇcasnˇe ˇctyˇrmi mikrofony um´ıstˇen´ ymi v r˚ uzn´ ych vzd´alenostech od mluvˇc´ıho a liˇs´ıc´ıch se tak odstupem uˇziteˇcn´eho sign´alu od ˇsumu (SNR). Promluvy jsou tak´e odliˇsiteln´e prostˇred´ım, ve kter´ ych byly nahr´any (automobil, kancel´aˇr a veˇrejn´a m´ısta). Kromˇe referenˇcn´ıch F0 hodnot s intervalem 1ms obsahuje tak´e okamˇziky zm´ınˇen´ ych pitch-mark˚ u a informaci o znˇelosti dan´eho u ´seku.

5.

PDA hodnot´ıc´ı framework

5.1.

Motivace

Motivac´ı pro vznik hodnot´ıc´ıho prostˇred´ı pitch-detection algoritm˚ u je nejenom moˇznost jejich vz´ajemn´e objektivn´ıho porovn´an´ı oproti zn´am´e referenci, ale napˇr´ıklad i hled´an´ı optim´aln´ıch parametr˚ u v nastaven´ı dan´eho PDA. Rozliˇcn´a hodnot´ıc´ı krit´eria a kategorie pot´e umoˇzn´ı vybrat vhodn´ y PDA pro potˇreby aplikace (napˇr. chybovost urˇcen´ı znˇelosti versus pˇresnost urˇcen´ı F0). 5.2.

Typy soubor˚ u na v´ ystupu PDA

Algoritmy detekuj´ıc´ı z´akladn´ı frekvenci m´ıvaj´ı obvykle sv˚ uj v´ ystup v jednom ze tˇr´ı pouˇz´ıvan´ ych form´at˚ u: Typ 1 oznaˇcuje nativn´ı form´at .pda soubor˚ u referenˇcn´ı datab´aze. V souborech nen´ı uvedena speci´aln´ı informace o ˇcasov´em kroku (pˇredpokl´ad´a se, ˇze je zn´ama a-priori), a proto 1

Pitch-mark je pˇresnˇe definovan´ y okamˇzik hlasivkov´eho cyklu

Jan Bartošek

5

soubor zaˇc´ın´a pˇr´ımo detekovanou frekvenc´ı, kaˇzd´a dalˇs´ı je pak na nov´em ˇr´adku. V hodnot´ach je zak´odov´ana informace o znˇelosti u ´seku (hodnota 0 znaˇc´ı ticho, hodnota 1 neznˇel´ y u ´sek a hodnota vˇetˇs´ı neˇz 1 validn´ı F0). Typ 2 je velmi podobn´ y typu 1, avˇsak na prvn´ım ˇra´dku obsahuje informaci o ˇcasov´em kroku v milisekund´ach. Typ 3 neodpov´ıd´a ˇz´adn´emu z pˇredchoz´ıch, protoˇze nem´a konstantn´ı ˇcasov´ y krok. Na jednom ˇr´adku obsahuje dvˇe ˇc´ısla ud´avaj´ıc´ı jednak detekovanou z´akladn´ı frekvenci a tak´e d´elku jej´ıho trv´an´ı (napˇr. v sekund´ach). Tento typ uloˇzen´ı dat komprimuje informaci, protoˇze eliminuje delˇs´ı sekvence stejn´ ych hodnot a je tak pamˇet’ovˇe nej´ uspornˇejˇs´ı. 5.3.

N´ avrh a implementace hodnot´ıc´ıho frameworku

Obr´azek 3 pˇredstavuje z´akladn´ı blokov´e sch´ema architektury frameworku. J´adrem je F0 referenˇcn´ı datab´aze, kter´a obsahuje jednak testovac´ı promluvy, tak jejich referenˇcn´ı .pda soubory s pr˚ ubˇehy F0. Seznam testovac´ıch soubor˚ u je pˇred´an spouˇstˇec´ımu bloku, kter´ y je odpovˇedn´ y za zavol´an´ı dan´eho PDA na kaˇzd´em ze zadan´ ych testovac´ıch soubor˚ u zvl´aˇst’. Je schopen volat r˚ uzn´e typy implementac´ı PDA (bin´arn´ı forma operaˇcn´ıho syst´emu nebo vyvol´an´ı prostˇred´ı Matlabu z pˇr´ıkazov´eho ˇr´adku pro PDA ve formˇe m-filu). Blok tak´e zaruˇcuje pˇr´ıpadn´ y pˇrevod vstupn´ıho souboru do form´atu vhodn´eho pro dan´ y PDA. V pˇr´ıpadˇe, ˇze algoritmus nem´a f´azi detekce znˇel´ ych/neznˇel´ ych u ´sek˚ u, lze mu pˇredˇradit blok V/UV, kter´ y mu d´ale ke zpracov´an´ı pˇred´a jen u ´seky vyhodnocen´e jako znˇel´e. N´aslednˇe je sjednocen typ v´ ystupn´ıch .pda soubor˚ u a ty jsou jednotlivˇe porovn´any s referencemi. V´ ysledkem je jednotliv´ y hodnot´ıc´ı soubory pro kaˇzdou promluvu zvl´aˇst’. Nad tˇemito v´ ysledky jsou d´ale poˇc´ıt´any glob´aln´ı hodnocen´ı a d´ale hodnocen´ı napˇr´ıˇc r˚ uzn´ ymi kategoriemi vstupn´ıch promluv (kan´aly, prostˇred´ı, pohlav´ı) a jejich kombinacemi (napˇr. jen kan´al 1 v prostˇred´ı automobilu). Framework byl implementov´an pod operaˇcn´ım syst´emem UNIXov´eho typu jako kombinace skript˚ u okolo referenˇcn´ı datab´aze napsan´ ych v multiplatformn´ım interpretovan´em jazyce PERL (vlastn´ı logika) a BASH (jednoduch´e souborov´e operace). Cel´e prostˇred´ı je tedy velmi snadno pˇrenositeln´e i na ostatn´ı platformy. M´ısto pouˇzit´e referenˇcn´ı datab´aze m˚ uˇze b´ yt s minim´aln´ım u ´sil´ım pouˇzita jin´a F0 referenˇcn´ı datab´aze a cel´ y framework by pak mohl naj´ıt upotˇreben´ı i mimo ˇreˇcov´e syst´emy (napˇr´ıklad v hudebn´ım odvˇetv´ı). 5.4.

Hodnot´ıc´ı krit´ eria

Existuje sada ust´alen´ ych krit´eri´ı pouˇz´ıvan´ ych na poli hodnocen´ı PDA, c´ılem pr´ace vˇsak bylo tak´e smysluplnˇe navrhnout nˇekter´a krit´eria nov´a. Chyby znˇelosti (Voiced Errors VE) se vyskytnou, pokud algoritmus znˇel´ yu ´sek prohl´as´ı za neznˇel´ y. Unvoiced Errors (UE) znaˇc´ı chybu opaˇcn´eho typu. Gross Error High (Low) je pˇr´ıpad chyby, kdy detekovan´a frekvence v Hz pˇrekroˇc´ı 20% toleranci shora (zdola) a nese oznaˇcen´ı GEH (GEL). Jelikoˇz tato tolerance je vzhledem k vyˇsˇs´ım poˇzadavk˚ um dneˇsn´ı doby k algoritm˚ um pomˇernˇe benevolentn´ı, byla novˇe zavedena tak´e tolerance 10% (oznaˇcen´ı GEH10 a GEL10). Nˇekdy jsou tak´e pouˇz´ıv´any souˇcty UE+VE a GEH+GEL. D´ale byly novˇe zavedeny tzv. halving HE ” (doubling DE)“ chyby, ke kter´ ym doch´az´ı pˇri odhadu dvojn´asobn´e (poloviˇcn´ı) frekvence. Gross Error chyby s obˇema toleranˇcn´ımi p´asmy jsou tak´e evidov´any v r´amci jednotliv´ ych nepˇrekr´ yvaj´ıc´ıch se frekvenˇcn´ıch p´asem (napˇr. pro ˇreˇcov´e sign´aly pˇet 2/3 okt´avov´ ych p´asem v rozmez´ı cca 60-560Hz). V literatuˇre lze tak´e nal´ezt statistick´a krit´eria (absolutn´ı rozd´ıl mezi stˇredn´ımi hodnotami reference a odhadu a tak´e absolutn´ı rozd´ıl jejich smˇerodatn´ ych odchylek) poˇc´ıtan´a v Hz pˇres cel´ y soubor promluv, avˇsak tyto hodnoty nemaj´ı d´ıky logaritmick´emu pr˚ ubˇehu vn´ım´an´ı frekvence pˇr´ıliˇsnou vypov´ıdaj´ıc´ı hodnotu. Proto je pouˇzito modifikovan´ ych statistick´ ych krit´eri´ı [2] - stˇredn´ı odchylka ∆% (10) a

6

Jan Bartošek

Obr´azek 3: Blokov´e sch´ema architekrury PDA hodnot´ıc´ıho frameworku smˇerodatn´a odchylka δ% (11) poˇc´ıtan´ ych v centech p˚ ult´on˚ u: ∆% =

δ% =

v u u t

N Fest (n) 1200 X log2 N n=1 Fref (n)

N Fest (n) 1 X [1200 log2 − ∆% ]2 N n=1 Fref (n)

(10)

(11)

Pro vysvˇetlen´ı jeˇstˇe dodejme, ˇze krit´erium VE+UE nen´ı pouh´ ym souˇctem obou chybovost´ı, ale je to souˇcet vˇsech chybnˇe klasifikovan´ ych u ´sek˚ u ku celkov´emu poˇctu u ´sek˚ u. Aby nedoch´azelo k akumulaci chyb, postupuj´ı do f´aze hodnocen´ı pˇresnosti jen takov´e znˇel´e u ´seky, kter´e byly spr´avnˇe klasifikov´any jako znˇel´e. To m˚ uˇze v urˇcit´ ych pˇr´ıpadech zv´ yhodnit v krit´eriu tolerance F0 algoritmy s vysokou hodnotou VE, protoˇze takov´e metodˇe se pˇresnost F0 poˇc´ıt´a jen z u ´sek˚ u, kter´e oklasifikovala jako znˇel´e a problematick´e u ´seky, kter´e mohou sniˇzovat pˇresnost napˇr. algoritm˚ um bez V/UV f´aze, jiˇz v hodnocen´ı pˇresnosti pˇr´ıtomny nejsou. 5.5.

V´ ysledky

Testov´any byly vˇsechny algoritmy uveden´e v kapitole 2.. V´ ysledky testovan´e sady algoritm˚ u na kan´alech 0 (nejvyˇsˇs´ı hodnota SNR, close-talk mikrofon) a 1 (klopov´ y mikrofon) uv´ad´ı tabulky 1 a 2. Nˇekter´e z algoritm˚ u (ACF time, AMDF, CEPS a MNBFC) byly naimplementov´any bez rozhodovac´ıch prah˚ u pro znˇel´e a neznˇel´e u ´seky a nebyl jim pˇredˇrazen ani ˇz´adn´ y samostatn´ y detektor. Proto detekuj´ı vˇsechny neznˇel´e u ´seky promluv jako znˇel´e a krit´erium Unvoiced Error u nich dosahuje hodnoty 100 %. Zaj´ımav´e je, ˇze algoritmy AMDF a CEPS dos´ahli, aˇckoliv pracuj´ı na zcela odliˇsn´em principu, t´emˇeˇr totoˇzn´ ych v´ ysledk˚ u ve vˇsech sledovan´ ych krit´eri´ı. Z pohledu v´ ysledk˚ u se tedy zdaj´ı ekvivalentn´ı, aˇckoliv lze nal´ezt informace o tom, ˇze nˇekter´e komplexnˇejˇs´ı PDA jsou zaloˇzeny pr´avˇe na jejich komplementaritˇe, kter´a se v naˇsem testu v˚ ubec neprojevila. Metoda MNBFC podala bohuˇzel v pˇresnosti v´ ysledky horˇs´ı neˇz ACF freq. Potvrdilo se, ˇze rostouc´ı zastoupen´ı ˇsumu v sign´alu je velk´ ym probl´emem (napˇr. obrovsk´e zhorˇsen´ı GEH pro Acf freq).

Jan Bartošek

PDA ACF freq ACF time AMDF CEPS DFE Wavelets MNBFC

7

VE UE VE+UE [%] [%] [%] 44,4 23,5 31,6 0 100 61,9 0 100 61,9 0 100 61,9 26,6 15,5 20,4 67,7 11,3 32,7 0 100 61,9

GEH [%] 1,2 4,7 0,6 0,6 8,4 2,5 4,8

GEL GEH10 [%] [%] 0,1 1,5 2,3 6,2 27,2 1,4 27,1 1,4 4,2 16,5 4,9 3,7 4,4 6,6

GEL10 [%] 0,18 3,5 28,3 28,1 8,9 6,0 6,6

DE [%] 0,4 0,8 0,1 0,1 0,2 1,1 0,4

HE [%] 0,06 1,3 16,2 16,0 1,3 3,9 2,8

Tabulka 1: Souhrnn´e v´ ysledky na kan´alu 0 PDA ACF freq ACF time AMDF CEPS DFE Wavelets MNBFC

VE UE VE+UE [%] [%] [%] 52,7 34,1 41,3 0 100 61,9 0 100 61,9 0 100 61,9 45,4 11,1 25,9 70,4 9,5 32,6 0 100 61,9

GEH [%] 23,3 28,8 10,8 10,1 8,5 14,3 29,1

GEL GEH10 [%] [%] 0,1 23,5 2,5 29,8 44 11,3 43,4 10,5 8,1 17,9 9,9 17,4 4,9 30,4

GEL10 [%] 0,2 3,4 45,2 44,7 13,1 11,6 6,5

DE HE [%] [%] 3,2 0,03 3,6 1,5 1,3 21,4 1,3 21,4 0,05 4,3 4,3 6,7 2,1 3,1

Tabulka 2: Souhrnn´e v´ ysledky na kan´alu 1 Nejrobustnˇejˇs´ı chov´an´ı naopak vykazuje DFE. 5.6.

Poˇ zadavky pro interpunkˇ cn´ı detektor

Pro u ´ˇcely pouˇzit´ı s dalˇs´ımi c´ıli v´ yzkumu (interpunkˇcn´ı detektor) budeme hledat algoritmus s minim´aln´ı chybou krit´eria UV (chybnˇe detekovan´e F0 na neznˇel´ ych u ´sec´ıch by mohly v´ yraznˇe naruˇsit pr˚ ubˇeh melod´emu), naopak urˇcitou chybovost VE lze tolerovat (jen chybˇej´ıc´ı bod melod´emu). Celkov´e v´ ysledky pˇresnosti (GEH) algoritm˚ u na kan´alu 0 by mˇely u ´ˇcelu postaˇcovat. 5.7.

Dodateˇ cn´ y experiment´ aln´ı V/UV blok

Na z´akladˇe pˇredchoz´ıch v´ ysledk˚ u byl implementov´an dle [4] jednoduch´ y detektor znˇel´ ych a neznˇel´ ych u ´sek˚ u vych´azej´ıc´ı z pomˇeru energie sign´alu (E) a poˇctu pr˚ uchodu nulou (ZCR). Energie je poˇc´ıt´ana z pˇredzpracovan´eho okna, kdy jsou zv´ yraznˇeny periodick´e struktury znˇel´ ych segment˚ u pomoc´ı kr´atkodob´e energetick´e ob´alky. Pˇredpis pro jeho v´ ystupn´ı funkci, kter´a se nakonec prahuje empiricky zjiˇstˇenou hodnotou, je vyj´adˇren vztahem (12). Lze se totiˇz domn´ıvat, ˇze pro znˇel´e u ´seky je hodnota funkce EZR vysok´a, protoˇze energie sign´alu je pomˇernˇe velk´a a poˇcet pr˚ uchod˚ u nulou je ve srovn´an´ı se ˇsumem n´ızk´ y. Naopak, pro neznˇel´e u ´seky s vysok´ ym ZCR a n´ızkou energi´ı sign´alu bude v´ ysledn´a hodnota mal´a. Hodnot´ıc´ı framework byl obohacen o modul umoˇzn ˇuj´ıc´ı evaluaci i samotn´ ych V/UV algoritm˚ u a takto implementovan´ y detektor dos´ahl na kan´alu 0 v´ ysledk˚ u o nˇeco horˇs´ıch neˇz metoda DFE (UE+VE: 24,3 % pro EZR versus 20,4 % DFE) a mohl by tak b´ yt u ´spˇeˇsnˇe pˇredˇrazen algoritm˚ um bez detekce znˇel´ ych u ´sek˚ u a vylepˇsit tak jejich celkov´e v´ ysledky.

8

Jan Bartošek

EZR[m] =

6.

E[m] ZCR[m]

(12)

Z´ avˇ ery

Kromˇe implementace z´akladn´ıch PDA byly zkoum´any i pokroˇcilejˇs´ı metody. Experiment´alnˇe bylo ovˇeˇreno, ˇze ˇcasov´e rozliˇsen´ı 16ms je pro sledov´an´ı pr˚ ubˇehu intonace ˇreˇci vhodnou hodnotou. Pro objektivn´ı hodnocen´ı F0 detekˇcn´ıch algoritm˚ u byl navrˇzen a implementov´an hodnot´ıc´ı framework obaluj´ıc´ı F0 referenˇcn´ı datab´azi. Byla navrˇzena sada r˚ uznorod´ ych krit´eri´ı umoˇzn ˇuj´ıc´ı detailn´ı anal´ yzu chov´an´ı algoritm˚ u. Podle oˇcek´av´an´ı kles´a u ´spˇeˇsnost vˇsech algoritm˚ u s klesaj´ıc´ım SNR. Metoda detekce zaloˇzen´a na sjednocen´ ych normalizovan´ ych korelac´ıch (MNBFC) bohuˇzel nepodala oˇcek´avan´e v´ ysledky co do pˇresnosti urˇcen´ı F0. V´ ysledky tak´e ukazuj´ı, ˇze nejvˇetˇs´ı slabinou je obecnˇe f´aze detekce znˇel´ ych a neznˇel´ ych u ´sek˚ u.

Podˇ ekov´ an´ı ˇ 102/08/H008 “Anal´ Tento v´ yzkum byl podporov´an z grant˚ u GACR yza a modelov´an´ı ˇ biomedic´ınsk´ ych a ˇreˇcov´ ych sign´al˚ u” a GACR 102/08/0707 “Rozpozn´av´an´ı mluven´e ˇreˇci v re´aln´ ych podm´ınk´ach”.

Reference [1] Bartoˇsek, J. Prozodie, zjiˇstˇen´ı a vyuˇzit´ı z´akladn´ıho t´onu v rozpozn´av´an´ı ˇreˇci. Semin´aˇre katedry teorie obvod˚ u, anal´yza a zpracov´an´ı ˇreˇcov´ych a biologick´ych sign´al˚ u - sborn´ık prac´ı 2009 (2009), 1–8. [2] Boˇril, H.; Poll´ak, P. Direct time domain fundamentalfrequency estimation of speech in noisy conditions. in Proceedings of EUSIPCO 2004 (European Signal Processing Conference, Vol. 1) (2004), 1003–1006. [3] Kotnik, B.; et al. Noise robust f0 determination and epoch-marking algorithms. Signal Processing 89. (2009), 2555–2569. [4] Kotnik, B.; H¨oge, H.; Kacic, Z. Evaluation of pitch detection algorithms in adverse conditions. 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 the Acoustical 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.

Tomáš Bořil

9

Multivariate Autoregressive Modelling of Causal Connections in EEG Tom´aˇs Boˇril Czech Technical University in Prague, Faculty of Electrical Engineering [email protected] Abstract: Conditional Granger causality and related methods have been used in neurophysiology in recent years for revealing causal relations in brain. Multivariate autoregressive models can model dynamic relations in electroencephalography and derived estimators can measure not only the strength of connections but also their directions which are important for understanding brain function during specific tasks. This paper propose an extension of the conditional Granger causality in order to improve a differentiation between strong and weak causal connections.

1.

Introduction

In the late ‘60s of the 20th century, econometrics started examination of causal relations among time series in order to create economic models and predict a behavior of variables like 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]. If prediction of one time series could be improved by knowledge of past values of another one, then we say the second series has a causal influence on the first one. Although this technique received a great acceptance, some studies like [7] warn against inconsiderate use of Granger causality which can lead to incorrect answers. For instance, the thoughtless use can find statistically significant causal relations which do not have its match in the real word sense. It is necessary to analyse data with autoregressive nature, the examined segment must be stationary and sufficiently long and all time series with possible influence must be included into the model. While meeting of such requirements can be difficult in econometrics, a whole new area opens for Granger causality in the form of neurophysiology at the end of 20th century. Functional magnetic resonance imaging (fMRI), electroencephalography (EEG) and magnetoencephalography (MEG) bring data which meets conditions of causality analysis very well and together with progress of computational technology it is possible to find causal relations in large data of brain activity records, e.g. [5, 2]. EEG processing is the oldest method of human brain function analysis. Among advantages, good availability should be mentioned (short waiting periods for investigation), low price and simplicity (no problems with patients suffering claustrophobia like with fMRI). However, the main advantage is a high temporal resolution of measured data (important for causal connectivity analysis) and direct relationship with electric activity

10

Tomáš Bořil

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

2.

Renormalized Granger Causality

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

p 

j=1

a1,j v1 (t − j) +

p 

j=1

a2,j v2 (t − j) + · · · +

p 

j=1

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

(1)

where t stands for an index in a discrete time instance and the prediction error ε1 is a white noise. To estimate the MVAR model coefficients one can minimize the variance of the prediction error via MVAR estimators discussed in [6]. Model order can be estimated by 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 via previous values of v1 . . . vk . Conditional Granger causality (CGC) is often defined for three variables like causality from v1 to v2 conditional on v3 [3]: Fv1 →v2 |v3 = ln



var 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 into the model a direct connection is distinguished from an indirect one when the information flow from the first variable to the second variable is completely mediated by the third variable. 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 not including the source variable, whereas the denominator includes the source variable. This definition corresponds with the original causal idea definition perfectly – if the inclusion of the previous values of the source variable decreases the variance of the prediction error, the conditional Granger causality is a non-zero positive value and we say the source variable Granger-causes the target variable. Let’s define maximal prediction ratio: MPRvi =



var (vi )

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



(5)

Tomáš Bořil

11

determining the ratio of how many times the variance of the variable is reduced in the form of 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 if predictions of variables are overall low and small changes in the value can lead to large result of the ratio. Let’s define renormalized Granger causality (RGC): RGCv1 →v2 = Fv1 →v2 |v3 ,...,vk

k 

MPRvi

(6)

i=1

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

3.

Experiments and Results

We have generated 2000 samples of three MVAR signals of order 5 with causal influences v1 → 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 signals with different signal-to-noise ratio levels (SNR) in order to analyse a behavior of CGC and RGC in noise conditions. In figures 1 to 3, we can see the recognition of causalities for different SNR levels and signal lengths, the recognition decreases rapidly with lower SNR. It is obvious with shorter signal segments the noise immunity of CGC and RGC decreases. However, RGC behaves much better in low SNR conditions, it suppresses false causality 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 wellaccepted conditional Granger causality has been proposed. Presented results show the renormalized Granger causality can distinguish better a case when the prediction is overall low from the case when data have highly predictive character. This property can be helpful in real data processing when the signal-to-noise ration value is unknown and the conditional Granger causality can estimate high causality because both values in numerator and denominator are small and little statistically insignificant changes can lead to great change in the ratio result. Renormalized Granger causality improves this property and a demonstration on artificial data is shown in this paper.

12

Tomáš Bořil

3.5

0.25

3 0.2 2.5

F

RGC

0.15 v1 v2 v1 v2 v3 v3

0.1

0.05

0 −20

0

20

60 40 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

2 v1 v2 v1 v2 v3 v3

1.5 1 0.5 0 −20

100

0

20

60 40 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

100

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

3.5

0.25

3 0.2 2.5

F

RGC

0.15 v1 v2 v1 v2 v3 v3

0.1

0.05

0 −20

0

20

40 60 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

2 v1 v2 v1 v2 v3 v3

1.5 1 0.5 0 −20

100

0

20

40 60 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

100

Figure 2: Additive white noise, signal length 1000 samples, CGC and RGC indexes

2.5

0.2

2

0.15

1.5

F

RGC

0.25

v1 v2 v1 v2 v3 v3

0.1

0.05

0 −20

0

20

40 60 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

v1 v2 v1 v2 v3 v3

1

0.5

100

0 −20

0

20

40 60 SNR (dB)

→ v2 → v3 → v3 → v1 → v1 → v2 80

100

Figure 3: Additive white noise, signal length 300 samples, CGC and RGC indexes

Tomáš Bořil

13

Acknowledgments Research described in this paper has been supported by SGS10/176/OHK3/2T/13 “Brain activity 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 Czech Technical University in Prague.

References [1] Boˇril, T. Revealing of Relations in EEG via Granger Causality. 3th International Student 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 influences revealed by Granger causality. Proceedings of the National Academy of Sciences of the United States of America 101, 26 (June 2004), 9849–9854. [3] Chen, Y.; Bressler, S. L.; Ding, M. Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. Journal of Neuroscience Methods 150, 2 (2006), 228–237. [4] Granger, C. W. J. Investigating causal relations by econometric models and crossspectral methods. Econometrica 37, 3 (July 1969), 424–38. [5] Liang, H.; Ding, M.; Nakamura, R. R.; Bressler, S. L. Causal influences in primate cerebral cortex during visual pattern discrimination. Neuroreport 11 (2000), 2875– 2880. [6] Schl¨ogl, 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).

14

Vladyslava Demchenko

Analýza variability srdečního rytmu v 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í.

Vladyslava Demchenko

2.

15

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 svodu 2. fotopletysmogram (PPG) - tento signál nebyl v této části zpracování využit 3. signál dýchání měřený z pásu, obepínajícího hrudník 4. 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 - 1 2. 5 - 1 - 5 - 1 (test hlubokého dýchání před apnoe) 3. fáze apnoe (zadržení dechu) - doba trvání 62s 4. 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.

16

2.2

Vladyslava Demchenko

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:



N −1

1 ∑  RRi1−RRi 2 [ms; -, ms, ms] N −1 i =1 kde N je počet všech RR intervalů za daný period měření. rMSSD=

(1)

Vladyslava Demchenko

17

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)

18

Vladyslava Demchenko

(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)

Vladyslava Demchenko

19

(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)

20

Vladyslava Demchenko

(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

Vladyslava Demchenko

21

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í dechu

4.

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.

22

Jaromír Doležal

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.

Jaromír Doležal

23

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.

24

Jaromír Doležal

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í.

Jaromír Doležal

25

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.

26

Jaromír Doležal

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

Jaromír Doležal

27

č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.

28

Jaromír Doležal

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

Radek Janča

29

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ší

30

Radek Janča

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).

Radek Janča

31

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 (DTFDirected 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, N

∑A x i =0

i

t −i

= ei

(1)

Ai 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 Ai 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),

 N f  H ( f ) =  ∑ A i exp(− j 2πi )  fs   i =0

−1

(2)

kde f je frekvence, fs je vzorkovací kmitočet a j = − 1 . Na základě DTF γ ij2 ( f ) popisuje „tok“ z kanálu j do kanálu i,

γ (f)= 2 ij

H ij ( f )

2

m

∑H k =1

ik

(f)

(3) 2

kde Hij(f) je element matice H(f) v řádku i a sloupci j. γ ij2 ( f ) je tedy normalizovaná DTF. Původní signál je filtrován do 8 subpásem FIR filtrem řádu 3 f s , pásma jsou popsána v tabulce 1.

Označení Pásmo [Hz]

δ 2-4

θ α1 α2 β1 4-8 8-10 10-12 12-18 Tabulka 1.: Subpásma EEG signálu

β2 18-25

γ1 25-48

γ2 52-85

32

Radek Janča

Pro každé subpásmo je definován střední „tok“ γ ij ( pásmo) = mean(γ ij ( f k ) ) 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. DF j =

1 m ∑ γ ij ( pásmo) m − 1 i=1

(4)

j ≠i

[

]

y j = DF j (δ ), DF j (θ ), DF j ( α 1 ), DF j ( α 2 ), DF j ( β 1 ), DF j ( β 2 ), DF j ( γ 1 ), DF j ( γ 2 )

(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. P=

N− N

(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.

1 low activity high activity

0.8

P

0.6 0.4 0.2 0 5.51

Obrázek 1.: Příklad vizualizace p-odhadu v kortikální mapě barevnou škálou.

5.512 t [s]

5.514

5.516 4

x 10

Obrázek 2.: Příklad rozdílu průběhu p-odhadu mezi dvěma kanály v čase začátku záchvatu.

Radek Janča

4

33

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 Neurophysiology, 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.

34

Jan Janda

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 logopedické parametry dětské řeči. Tyto parametry byly získány z analýz dětských promluv 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 anatomickou př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 jednotlivých měkkých a tvrdých tkání relativně k délce vokálního traktu. Jedná se například o 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é struktury v 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ího traktu během prvních 20 let života. Délkou vokálního traktu se zde rozumí délka křivky v mediální rovině vedené od středu hlasivek po její průsečík s tečnou rtů. Během vývoje se 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 o tloušťku maxilárního a mandibulárního rtu, délku tvrdého a měkkého patra, délku jazyka, délku mandibuly, a vzdálenost hrtanu, příklopky hrtanové a jazylky od spina nasalis 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 částech lineární model růstu. Ze závěrů anatomických měření můžeme předpokládat některé akustické důsledky. Jednak s rostoucí délkou celého vokálního traktu můžeme předpokládat očekávaný pokles formantů [4, 2], dále díky absenci pohlavního dimorfismu v uvedených charakteristikách do 6 let věku můžeme do jisté míry očekávat i pohlavní nezávislost souvisejících akustických parametrů. Také vedle zvyšující se artikulační zručnosti bude mít vliv na koartikulaci i zvětšující se prostor pro pohyb jazyka.

Jan Janda

35

Obrázek 1: Výstup magnetické rezonance v mediální rovině. Vlevo 7 měsíční dívka, vpravo dospě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íky dolů pro muže a stínované nahoru pro ženy). Převzato z [6]

36

Jan Janda

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řízena databáze dětských promluv. Obsah je již cíleně navržen tak, aby se v databázi vyskytovaly 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. Slova dě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ána slova: babička, máma, čokoláda, sluníčko, popelnice, košile, silnice, Rákosníček, hamburger, velryba, ucho, ježek, ředkvičky, fotbalista. Pro posouzení základních fonetických parametrů 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ým spontá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. 3.1.

Věkově závislé akustické charakteristiky 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 charakteristiku v souvislosti s věkem člověka. Nabývá hodnot od 500 Hz u nejvyšších dětských hlasů a s věkem může u mužů klesnout až na hodnoty kolem 80 Hz. Analýza základního hlasivkového tó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 parametry time step=0.0, pitch floor =75 Hz a pitch ceiling=600 Hz. Výsledné hodnoty byly ověřeny v 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á tuto zá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ůžeme vyjá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/.

Jan Janda

37

F0 300

250

Hz

200

150 dívky chlapci

100

50

2

4

6

8

10

12

14

16

vek

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

Formantove pole 2800 2600 2400

I 2200

F2(Hz)

2000 1800 1600 1400

U

1200

A

1000 800 300

400

500

600

700

800

900

1000

F1(Hz)

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

38

3.2.

Jan Janda

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 nesrozumitelnosti zkoumané 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 vykazovala kumulovaná vzdálenost v prostoru RASTA-PLP a kepstrálních LPC koeficientů (obr. 5). Cumsum DTW − slovo "Hamburger" 0.4 CLPC RASTA−PLP

0.35 0.3

Cumsum

0.25 0.2 0.15 0.1 0.05 0

4

6

8

10 vek

12

14

16

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

4.

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

Tabulka 1 shrnuje doposud zkoumané akustické a logopedické charakteristiky. Jednotlivé příznaky jsou seřazeny podle míry korelace s věkem. Sloupec Ho obsahuje hodnoty hladin 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ěkem okolo 0,66. Navržená metoda měření srozumitelnosti pomocí kumulované vzdálenosti v algoritmu DTW dokonce vykazuje korelaci s věkem až 0,72. Oproti výsledkům pořízeným na předchozí datábazi zde nevykazují nijak významnou věkovou závislost spektrální charakteristiky 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.

Jan Janda

39

Charakteristika DTW fotbalista DTW popelnice DTW čokoláda DTW velryba DTW hamburger F0 /I/ F0 /U/ F0 /A/ DTW Rákosníček F2 /A/ F1 /U/ F1 /A/ F1 /I/ F2 /I/ F2 /U/

r -0,72 -0,71 -0,70 -0,69 -0,66 -0,66 -0,65 -0,57 -0,53 -0,40 -0,38 -0,34 -0,30 0,20 0,18

H0 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,001 0,005 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 biomedical and speech signals.

Reference [1] GEROSA M.,LEE S. et al.: Analyzing Children’s Speech: an Acoustic Study of Consonants and Consonant-Vowel Transition. In Proc. of the International Conference on Acoustic, Speech and Signal Processing,(Tolouse, France), May 2006 [2] HUBERT J. E., STATHOPOULOS E. T et al.: Formants of children, women, and men: the effects of vocal intensity variation. Journal of the Acoustical Society of America 1999 Sep;106(3 Pt 1):1532-42. [3] LEE S., POTAMIANOS A., NARAYANAN S.: Acoustic of children’s speech: Developmental changes of temporal and spectral parameters. Journal of the Acoustical Society of America, pp. 1455-1468, Mar. 1999 [4] MENARD L., SCHWARTZ J.-L., Articulatory–acoustic relationships during vocal tract growth for French vowels: Analysis of real data and simulations with an articulatory 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 properties of children’s speech Multimedia Signal Processing, 2007. MMSP 2007. IEEE 9th 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 during early childhood: A magnetic resonance imagining study. Journal of the Acoustical Society of America 117, January 2005, pp. 338-350

40

Ján Janík

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 − 1;− w and w;−1 . 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). N −1  2n + 1  C (k ) = ε (k )∑ x(n ) cos kπ  (1)  2N  n =0  1 for k = 0  ε (k ) =  N (2) 2  for k ≠ 0  N

Ján Janík

41

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 0; N − 1 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.  2n + 1  cos kπ   2N 

(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.

42

Ján Janík

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)

(1 − w )(

2

)

(

)

 dZ  2 2 2 w −w   = 4m w 1 − Z ,  dw 

2

2

2 p

(4)

where Z in the interval w = − 1;1 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 w p = − w p , zcos function is equal to a standard cos. We have developed the recursive evaluation of the coefficients a (wp) 2µ

m  2n + 1  2µπ  . z cos(n,2m, w p ) = ∑ a2 µ (w p )cos  2N  µ =0

(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)

(1 − w )(w 2

2

−w

2 p

)(

)

2

(

)(

)

 dZ  2 2 2 2 2 w −w   = (2m + 1) 1 − Z w − wm .  dw  2

2 s

(6)

In eq. (6) the function Z stands for an odd zcos (7) within the interval w = − 1;1 . We have developed the recursive evaluation of the coefficients a2µ+1(ws) m  2n + 1 (2µ + 1)π  . z cos(n,2m + 1, ws ) = ∑ a2 µ +1 (ws )cos  2N  µ =0

(7)

Ján Janík

43

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.

44

Robert Krejčí

Optimalizace algoritmů rozpoznávačů řeči s využitím architektur vícejádrových signálových procesorů Robert Krejčí České vysoké učení v Praze, Fakulta elektrotechnická [email protected] 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. Úvod Existuje 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ům Souč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.

Robert Krejčí

45

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ý standard 2. 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 serveru Server 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 server V 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.

46

Robert Krejčí

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ástroje Donedá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ěže Rozhodnutí 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.

Robert Krejčí

47

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říkonu Vý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.

48

Robert Krejčí

3.1. Optimalizace výpočtu parametrizace MFCC Jední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=



P





2 i g j cos  j−0,5 ,1iM −1, ∑ P j=1 P

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: P

c i=∑ g j⋅costab[i , j ]=DOTPROD  g , costab[i]  , j =1

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.

Obrázek 1: Vliv vybraných optimalizačních metod na dobu výpočtu parametrizace MFCC 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.

Robert Krejčí

49

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í funkce Součá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ů: S

b j o t =∏ s=1

N o ;  ,  =

[∑

s

]

Ms

m=1

c jsm N o st ;  jsm ,  jsm 

1

−

n

2  ∣∣

1 2

;

−1

e o−'  o− ,

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í: M s ≈40

ln b j o t  = ln

∑

m=1



K ≈51

2

exp X jm ∑ [ otk − jmk ⋅y jmk ] k=1



3.2.1. Optimaizovaná funkce – TVAR 1 – MAC V 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: p k = o 2k , z jmk = y 2k , v jmk = −2  jmk y 2jmk , X M



ln b j o t  = ln ∑ exp X i=1

K≈ 51

 ∑ z jmk⋅p jmk v jmk⋅o jmk

jm1

k =1

K≈51

2 2 jm1 = X jm  ∑  y jmk⋅ jmk  k=1



Zů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-MAC Naš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

50

Robert Krejčí

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

2 k , m ,i

51

, t k , m ,i = −2 k ,m ,i , U k , m = X k , m∑  y k ,m ,i⋅k , m , i  2

i =1

M ≈40

ln b k o  = ln

2

∑

m=1



K≈51



exp U k , m ∑ z k , m ,i⋅o i⋅oit k , m , i  i=1

Tato 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.

Obrázek 2: Blokové schéma modifikované funkce b(o) - TVAR 2 - ADD-MAC 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ím V 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.

Robert Krejčí

51

4. Závěry Tento č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: . [TMS06] Texas Instruments: TMS320C6000 Programmer's Guide [online]. 2006 [cit. 2010-09-07]. URL: . [TMS10] Allred Daniel: C6Run DSP Software Development Tool [online]. 2010 [cit. 2010-11-24]. URL: . [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: .

52

Ondřej Kučera

Zpětnovazební regulace mikromechanických experimentů Ondřej Kučera České vysoké učení technické v Praze, Fakulta elektrotechnická Ústav fotoniky a elektroniky, Akademie věd ČR [email protected] Abstrakt: Scanning Probe Microscopy, a large family of microscopies, is a sensitive technique with unprecedented resolution. However, the accuracy of the measurement is strongly influenced by the setting and calibration of the apparatus. In this article, we present system analysis approach to the problem of 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 velmi malou sondou (typicky s poloměrem křivosti v řádu desítek, nebo jednotek nanometrů) a vzorkem. 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í sondou je také široce používaná jako nezobrazovací technika umožňující měření lokální vlastností vzorku, jako je Youngův modul, permitivita, hustota náboje, přilnavost, atd. Pro udržení požadované přesnosti měření se používá zpětná vazba. V tomto článku uvádíme analýzu problému nastavení zpětné vazby a přizpůsobení tuhosti sondy.

Ondřej Kučera

2.

53

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-vzorek následující rovnicí pohybu d2 z dz (1) + (k + κ)z = κzo , +β 2 dt dt kde m je efektivní hmotnost, k je tuhost, a β je tlumení sondy, zo je topografie nezatíženého vzorku, κ je tuhost povrchu vzorku a z je vzdálenost sondy od povrchu. Chceme-li být přesní, musíme konstatovat, že κ obsahuje i tuhost interakce, ale to není v lineárním režimu důležité. m

Cantilever

β

k

m

z

κ Sample

zo

F

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 pracovního bodu a zvýšení citlivosti v případě vysokého poměru k/κ. Můžeme ji do rovnice 1 zavést následujícím způsobem   

t

  

∫   d2 z dz  m 2 +β + (k + κ)z = κ zo − P z − I z(τ )dτ  , dt dt   

|

0

{z

zp.vazba

(2)

}

kde P a I jsou proporcionální a integrální zisk zpětné vazby. Přenos zpětnovazebně řízeného systému sonda-vzorek můžeme psát jako sκ , (3) H(s) = 3 s m + s2 β + s(k + κ + P κ) + Iκ kde s je komplexní Laplaceův parametr. V ideálním případě L−1 H → 0 a informace o topografii je obsažená v signálu zpětné vazby s přenosem Hh (s) =

s3 m

+

s2 β

κ(sP + I) , + s(k + κ + P κ) + Iκ

(4)

a odpovídající amplitudovou frekvenční charakteristikou v u u |Hh (jω)| = t

kde ω je úhlová frekvence.

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

(5)

54

Nelineární model

F(z)

3.

Ondřej Kučera

0

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), tuhost interakce je pak ( ) 12

4ε ζ κ(z) = κo ∥ ∇  z z

−

( )6 

ζ z

,

(6)

kde ε je konstanta (hloubka potenciálové jámy), ζ představuje konečnou (!) vzdálenost nulového potenciálu a κo je tuhost samotného vzorku. Existují samozřejmě i další nelinearity v systému, ale lze je mnohdy kompenzovat, či zanedbat. Výsledné rovnice je tř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íku namí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řesto může být citlivost naopak významně zhoršená špatným nastavením zpětné vazby, a to zejména pro vysoký poměr κ/k. Hodnota |Hh (jω)| (rov. 5) vyjádřená v procentech představuje vertikální přesnost měření. Cílem nastavení zpětné vazby je upravit P a I zisky způ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/13 GA ČVUT.

Ondřej Kučera

55

Transfer H [−]

1 0.75 0.5 0.25 0 0.01

0.1

1

10

100

Relative sample stiffness κ/k [−]

Obrázek 3: Efekt nastavení zpětné vazby na přenos v propustném pásmu v závislosti na pomě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é mezi 0,1 a 100. Červená křivka značí případ bez použití zpětné vazby.

Transfer H [dB]

50

0

−50

−100

0

10

2

10

4

10

Frequency f [kHz]

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), respectively. Jsou ukázány tři odlišné hodnoty poměru tuhostí sondy a vzorku: k = κ = 0.05 N/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 and medicine. Journal of Molecular Recognition 20, 6 (2007), 418–431.

56

Ondřej Kučera

[5] Poggi, M.; Gadsby, E.; Bottomley, L.; King, W.; Oroudjev, E.; Hansma, H. Scanning probe microscopy. Anal. Chem 76, 12 (2004), 3429–3444. [6] Salapaka, S.; Salapaka, M. Scanning probe microscopy. IEEE Control Systems Magazine 28, 2 (2008), 65–83. [7] Stark, R.; Drobek, T.; Heckl, W. Thermomechanical noise of a free v-shaped cantilever for atomic-force microscopy. Ultramicroscopy 86, 1-2 (2001), 207–216.

Tomáš Lustyk

57

Hodnocení koktavosti a experimenty s adaptivním prahem u bayesovského spektrálního detektoru Tomáš Lustyk České vysoké učení v Praze, Fakulta elektrotechnická [email protected] 22. listopadu 2010

Abstrakt: V příspěvku je uvedeno stručné shrnutí výsledků výzkumu automatického hodnocení neplynulosti řeči. K hodnocení jsou užity parametry jak v časové 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 koeficientu 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 pomocí 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.

58

Tomáš Lustyk

Č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 foniatry, kteří tíži poruchy popsali pomocí 5-ti stupňové klasifikace na základě relativní četnosti neplynulý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.

Tomáš Lustyk

3.2

59

Parametr ve frekvenční oblasti

Parametr zkoumající neplynulou řeč ve frekvenční oblasti, se snaží zohlednit přítomnost prolongací. 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 obsahuje významné i méně významné změny. Pro jejich odlišení je třeba zavést práh. Práh je odvozen jako zlomek vybraného maxima ve výstupu detektoru (např. pátého nejvyššího). Parametr 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 oblast zkoumají další parametry ve frekvenční oblasti jako např. směrodatná odchylka z intervalů 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. 3

BD

2 1

0

0

10

20

30

40

50

BD

3 2 1 0 11

12

13

14

15 c a s [s ]

16

17

18

19

Obrázek 2: Ukázka výstupu bayesovského detektoru a pevného prahu. Nahoře: Výstup detektoru 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 technickou 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.

60

Tomáš Lustyk

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

0,297 0,446 0,454 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) Maxima 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 detektoru 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. 1

0. 8

0. 6

0. 4

0. 2

0 0

0. 2

0.4

0.6

0. 8

1

1.2

1.4

1. 6

1.8

2

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 zkreslení 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.

Tomáš Lustyk

61

Obrázek 4: Shora dolů: Výstup Bayesova detektoru a práh s regulací pomocí směrodatné odchylky (č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 nahrá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 pomocí 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 nejvyšší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í koeficient -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 katedř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é parametry 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.

62

Tomáš Lustyk

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 kombinují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í adaptivní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ŠMT MSM6840770012 „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 Studentské 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 elektrotechnická, Č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.

Martina Nejepsová

63

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.

30

20

10

0

2

3

4

5

6

7

8

9

10

11

12

13

Věk [rok] Obrázek 1: Počet promluv v jednotlivých věkových kategoriích

64

Martina Nejepsová

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 37 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.

Martina Nejepsová

65

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. 1-I

číslo hodnotitele ø odchylka [rok] pořadí

2-I

2-II

2-III

2-IV

1,49 -0,90 -0,16 -0,56 -0,22 7 5 1 3 2 Tabulka 3: Hodnocení průměrné odchylky v TEST 1

číslo hodnotitele ø kladná odchylka [rok] ø záporná odchylka [rok]

3-I

3-II

-1,24 6

-0,88 4

1-I

2-I

2-II

2-III

2-IV

3-I

3-II

1,83 -0,85

1,27 -1,38

1,18 -1,31

1,17 -1,10

0,91 -1,00

0,72 -0,85

0,58 -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 číslo hodnotitele

1-I

2-I

2-II

2-III

2-IV

3-I

3-II

ø absolutní odchylka [rok] 1,71 1,26 1,06 0,90 0,74 1,35 pořadí 7 5 4 2 1 6 Tabulka 5: Hodnocení absolutní průměrné odchylky v TEST 1

1,02 3

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).

66

Martina Nejepsová

skupiny hodnotitelů

1

2

3

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

1

2

3

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

ø odchylka [rok] 0,00 -0,37 XXX XXX XXX ø absolutní odchylka [rok] 0,66 XXX XXX XXX 0,00 Tabulka 9: Odchylka hodnocením v TEST 1 a v TEST 2

3-I

3-II

-0,08 0,23

0,10 0,21

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] [2] [3] [4] [5] [6] [7]

Otčenášek, Z.: O subjektivním hodnocení zvuku, AMU, Praha 2008 Melka, A: Základy experimentální akustiky, AMU, Praha 2005 Janda, J.: Studie věkově závislých akustických parametrů v dětské řeči, k odborné rozpravě, ČVUT, Praha 2010 http://www.cs.waikato.ac.nz/ml/weka/ http://www.fon.hum.uva.nl/praat/ http://www.speech.kth.se/wavesurfer/ http://www.mathworks.com/

Studie

Václav Procházka

67

Příprava a analýza Českého Web 1T 5-gram korpusu pro použití v jazykovém modelu Václav Procházka České vysoké učení technické v Praze, Fakulta elektrotechnická [email protected] Abstrakt: V této práci je popsán postup analýzy českého Web 1T 5-gram korpusu. Korpus byl analyzován a byly vyhodnoceny jeho základní charakteristiky před a v průběhu zpracování. Při zpracování byl slovník korpusu filtrován různými metodami, tak aby pokud možno obsahoval pouze smysluplná slova. Z pročištěného korpusu byly vygenerovány jazykové modely pro Large Vocabulary Continuous Speech Recognition (LVCSR) a spočítána jejich perplexita. Pro srovnání stejnými filtrovacími postupy byl také zpracovaný 5gramový 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 (Large Vocabulary Continuous Speech Recognition (LVCSR) LVCSR). Tento článek dále popisuje analý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 zdlouhavé, 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ůstem počtu webových stránek neustále zvětšuje, a tím také roste atraktivita internetu jako jejich zdroje, 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 firmou Google z webových stránek vydaný v Linguistic Data Consortium (LDC) [4]. Použitím tohoto korpusu je možné ušetřit mnoho náročné práce se sběrem korpusových dat z WEBových stránek. Velké množství dat (∼136 miliard slov) je již předzpracováno do podoby n-gramů s následujícími vlastnostmi [6]: • kódování n-gramů je sjednocené (UTF8), • obsahuje n-gramy 1 až 5 řádu,

68

Václav Procházka

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

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

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

• speciální token <S> resp. 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 pouze na konci n-gramu,

• čísla o více než 16 číslicích a slova delší než 32 znaků jsou nahrazena tokenem , • slova se znaky z neevropských jazyků, neplatnými UTF8 znaky a kontrolními ASCII znaky jsou také nahrazena tokenem , • slova vyskytující se méně než 40 krát jsou nahrazena tokenem

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

Prvotní náhled do slovníku ukázal, že přes výše uvedené rozsáhlé předzpracování se v tomto korpusu vyskytuje velké množství nevhodných nebo dokonce nesmyslných slov jako 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. Absence n-gramů jdoucích přez hranice vět a ořezání nejméně častých n-gramů, v některých pří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 synchronní 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 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 korpus ve větší míře pouze číselné výrazy, jiná nevhodná slova jako například náhodné shluky písmen jsou méně časté.

Václav Procházka

69

původní korpus unigramy trigramy

9 786 424 117 264 988

pročištěný korpus po prvním kole po druhém kole 1 804 682 954 771 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 unigramy trigramy

2 554 028 189 152 100

pročištěný korpus po prvním kole po druhém kole 1 799 0057 859 403 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 :

– všechna slova, která obsahovala jiné znaky než písmena než ty z české abecedy a 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 druhém kole), • interpunkce byla vynechána bez náhrady (jen Web 1T 5-gram, SYN2006PUB 5gram interpunkci neobsahuje), • nahrazení tokenu tokeny a <s> (jen SYN2006PUB 5-gram , Web 1T 5-gram již oba tokeny obsahuje), • totožné n-gramy vzniklé v předchozích krocích byly sjednoceny a jejich četnosti sečteny. Operace vynechání interpunkce ve svém důsledku způsobuje, že řád zpracovávaného ngramu 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 zahozen zcela. 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.

70

Václav Procházka

slovník 60 000 120 000 180 000 240 000

OOV rate cut off 17.76% 40 12.14% 40 9.44% 40 7.63% 40

bigramový model perplexity bigramy 49 507 11 437 940 126 001 14 904 654 153 682 16 522 260 206 831 17 666 327

trigramový model perplexity trigramy 162 509 35 166 960 255 144 42 125 506 664 866 43 882 665 935 077 45 443 605

Tabulka 3: Parametry jazykových modelů natrénovaných z pročištěného českého Web 1T 5-gram korpusu po prvním kole. cut off 1 slovník 60 000 120 000 180 000 240 000

OOV rate cut off 10.93% 1 6.36% 1 4.18% 1 3.06% 1

bigramový model perplexity bigramy 932 17 587 483 1 226 29 928 179 1 432 35 987 677 1 571 39 963 558

trigramový model perplexity trigramy 928 65 783 406 1 197 96 243 374 1 400 108 421 863 1 521 115 381 912

Tabulka 4: Parametry jazykových modelů natrénovaných z pročištěného SYN2006PUB 5-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 konec věty, byl tento token nahrazen dvojicí tokenů a <s>. Tím došlo ke zvýšení řádu n-gramu resp. k jeho následnému rozdělení na dva n-gramy se stejným řádem a počtem vý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 odstranění krátkých slov, ve slovníku Web 1T 5-gram korpusu stále zbylo mnoho špatných nebo nevhodných slov, například cizí slova, slova s překlepy, slova s odstraněnou diakritikou nebo i jen pouze nesmyslné shluky písmen (kód výrobku, jméno firmy). Tyto problémy jsou ř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ích n-gramů jsou shrnuté v tabulce 1. Perplexity na jazykových modelech postavených z trigramů každého korpusu jsou shrnuty v tabulkách 3 a 4 převzatých z [9]. Perplexita byla počítána na podmnožině přepisů z České SPEECON databáze [8], [7]. Tato podmnožina obsahuje 148 557 slov v 14 914 větách. Perplexity napočítané na jazykových modelech z 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é vykazují 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í diakritikou. Typické vlastnosti internetových stránek mohou být dalšími důvody. Například je mož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ů. Odstranění prvních dvou zmíněných problémů, cizích slov a slov s chybějící diakritikou bylo provedeno v druhém kole zpracování.

Václav Procházka

slovník 60 000 120 000 180 000 240 000

71

OOV rate cut off 13.02% 40 8.14% 40 5.75% 40 4.48% 40

bigramový model perplexity bigramy 337 226 1 512 233 570 781 1 797 817 525 619 1 920 150 1 184 050 1 985 640

trigramový model perplexity trigramy 120 476 1 765 425 249 281 1 991 280 380 179 2 079 416 507 689 2 123 533

Tabulka 5: Parametry jazykových modelů natrénovaných z pročištěného českého Web 1T 5-gram korpusu po druhém kole. cut off 1 slovník 60 000 120 000 180 000 240 000

OOV rate cut off 10.77% 1 6,71% 1 4,78% 1 3,77% 1

bigramový model perplexity bigramy 1 727 30 205 968 2 215 38 181 993 2 578 41 981 434 2 848 44 146 507

Tabulka 6: Parametry jazykových modelů natrénovaných z pročištěného SYN2006PUB 5-gram korpusu po druhém kole. V druhém kole zpracování byly kromě operací nových zachovány i všechny operace předchozího kola. Na rozdíl od prvního kola, všechny pročišťovací operace, které ze své podstaty 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é jsou v n-gramech nahrazeny speciálním tokenem . Zpracování a čištění slovníku v první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ím jsou v tabulkách 5 a 6. U modelů se slovníkem o velikosti 60000 slov došlo podle očekávání k poklesu OOV. U Web 1T 5-gram modelů byl pokles výrazný, u SYN2006PUB 5-gram modelů jen nepatrnému. Pro slovníky ostatních velikostí OOV modelů Web 1T 5-gram korpusu rovnoměrně klesalo, oproti tomu u modelů ze SYN2006PUB 5-gram OOV mírně vzrostlo, důvodem je pravděpodobně příliš agresivní filtrování, které odstranilo i správná slova. Perplexity trigramových modelů postavených z Web 1T n-gramů se trochu zlepš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ěpodobně způsobeno vlastnostmi korpusu: např. ořezáním, obsahem interpunkce a krokem vypuštění interpunkce, kdy pak 5-gramy zredukované na bigramy již dostatečné dobře nereprezentují původní data. Drobný vzrůst perplexit bigramů ze SYN2006PUB 5-gram korpusu 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íčky nástrojů Hidden Markov Model Toolkit (HTK) a Stanford Research Institute Language Modeling Toolkit (SRILM).

72

Václav Procházka

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ším zaměř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 LVCSR nebo samostatně. Část nástrojů pro práci s jazykovými modely byla použita především v 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ším v oblastech rozpoznávání řeči, statistického značkování a segmentace či strojového překladu. 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žít i v HTK. Tyto nástroje byly požity pro tvorbu a analýzu jazykových modelů v druhém kroku.

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-gram korpusu, 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 daty z 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 korpusu již vyřešené. • Perplexita modelů vytvořených z n-gramů českého Web 1T 5-grams po první fáze byla velmi vysoká, více než 105 v porovnání s výsledky perplexit SYN2006PUB 5gram korpusu které byly mezi 900 a 1600. Takto vysoké perplexity českého Web 1T korpusu 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 automatickou kontrolu pravopisu. Při tomto filtrování došlo k významnému poklesu slov ve slovníku. OOV a perplexity trigramových modelů z českého Web 1T 5-grams korpusu 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š agresivnímu automatickému filtrování. • Další pokračování těchto pokusů s čištěním n-gramových korpusů a tvorby jazykový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č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í”.

Václav Procházka

Reference [1] GNU Aspell. accesed: 17.11.2010, http://aspell.net/. [2] Český národní korpus (Czech National Corpus) - SYN2006PUB. Institute of the Czech 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. Linguistic Data 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 validation. 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 Comparison with Czech National Corpus Data. In LNAI 6231, (TSD 2010) 2010, P. Sojka et al., Ed., Springer-Verlag Berlin Heidelberg, pp. 181–188. [10] Stolcke, A. Srilm—an extensible language modeling toolkit. In In Proceedings of the 7th 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.

73

74

Barbora Richtrová

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.

Barbora Richtrová

2.

75

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.

76

Barbora Richtrová

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.).

Barbora Richtrová

3.

77

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

78

Jan Rusz

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

Jan Rusz

79

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 Jitter:local (%) Jitter:RAP (%)

Description

NHR (-) HNR (dB) RPDE (-) DFA (-)

Average absolute difference between consecutive periods, divided by the average period. Relative Average Perturbation, the average absolute difference between a period and the average of it and its two neighbours, divided by the average period. 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. Average absolute difference between consecutive differences between consecutive periods, divided by the average period. Average absolute difference between the amplitudes of consecutive periods, divided by the average amplitude. Average absolute base-10 logarithm of the difference between the amplitudes of consecutive periods, multiplied by 20. 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. 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. 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. Average absolute difference between consecutive differences between the amplitudes of consecutive period. Noise-to-Harmonics-Ratio, the amplitude of noise relative to tonal components. Harmonics-to-Noise-Ratio, the amplitude of tonal relative to noise components. Recurrence Period Density Entropy, the extension of the concept of periodicity. Detrended Fluctuation Analysis, the stochastic self-similarity of the turbulent noise.

PPE (-)

Pitch Period Entropy, the inefficiency of voice frequency control.

Jitter:PPQ5 (%) Jitter:DDP (%) Shimmer:local (%) Shimmer:local (dB) Shimmer:APQ3 (%)

Shimmer:APQ5 (%)

Shimmer:APQ11 (%)

Shimmer:DDA (%)

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

80

Jan Rusz

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

Jan Rusz

81

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. PD 0.5

HC (a)

signal amplitude

0.5

−0.5

−0.5

pitch period

6.7

6.72

6.74

pitch period

1.42 1.425 1.43 1.435

t (s) H norm = 0.42

t (s) 0.15

(c)

0.1

P(T)

P(T)

0.15

0.05 0 0

100

0.05 0 0

200

100

1

(e)

0.5

0 2

3

0

4

P(rp)

P(rp)

1

0.5

0 −2

0

rp (semitones)

(f)

2

3

4

log L (g)

PPE = 0.88

α norm = 0.56

0.5

log L 1

200

T (samples)

log F(L)

log F(L)

α norm = 0.64

(d)

H norm = 0.32

0.1

T (samples) 1

(b)

0

x

x

0

signal amplitude

2

(h)

PPE = 0.23

0.5

0 −2

0

2

rp (semitones)

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.  

82

Jan Rusz

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

Jitter:local (%) Jitter:RAP (%) Jitter:PPQ5 (%) Jitter:DDP (%) Shimmer:local (%) Shimmer:local (dB) Shimmer:APQ3 (%) Shimmer:APQ5 (%) Shimmer:APQ11 (%) Shimmer:DDA (%) NHR (-) HNR (dB) RPDE (-) DFA (-)

Mean 1.77 1.04 0.89 3.11 8.05 0.76 4.07 4.56 6.41 12.22 0.21 14.73 0.34 0.61

(SD) 1.78 1.08 0.83 3.24 5.03 0.44 2.56 2.98 3.86 7.69 0.28 6.99 0.09 0.02

Mean 0.62 0.34 0.31 1.02 3.43 0.32 1.73 1.91 2.78 5.19 0.03 22.37 0.27 0.61

(SD) 0.58 0.39 0.25 1.16 2.58 0.23 1.43 1.44 1.78 4.29 0.04 5.21 0.07 0.02

Differences between groups P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .001 P < .01 P = .086

PPE (-)

0.48

0.28

0.30

0.12

P < .01

Measurement

PD

HC

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

0.02

1

2

0 0

3

Jitter (%)

0.05 0 0

2

0.5

0.03 0.02 0.01 0 0

1

NHR (−)

20

40

HNR (dB)

0.04 0.1

P(DFA)

P(RPDE)

0.1

Shimmer (dB)

0.03 0.02 0.01 0 0

1

0.2

0.4

RPDE (−)

0.6

P(PPE)

0 0

0.05

P(HNR)

P(Shimmer)

P(Jitter)

0.04

0.04

P(NHR)

0.1

0.06

0.05

0 0.5

0.6

DFA (−)

0.7

0.03 0.02

PD

0.01

HC

0 0

0.5

1

1.5

PPE (−)

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.

Jan Rusz

83

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 postlevodopa 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).

84

Jan Rusz

[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 harmonicsto-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, 197202, (1995).

Jan Sova

85

Detekce náhlých změn v intrakraniálním EEG pomocí vlastních čísel korelační matice Jan Sova České vysoké učení technické v Praze, Fakulta elektrotechnická [email protected] Abstrakt: Práce ukazuje možnosti využití korelační matice a vlastních čísel korelač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 epileptologie, 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 v populaci 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 obrazu záchvatů. V lékařské veřejnosti se uvádí klasifikace dle oblastí mozku, kde záchvaty vznikají (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, generalizované záchvaty bez křečí, generalizované záchvaty s křečemi) apod. Nebezpečné jsou tzv. farmakorezistentní epilepsie, tj. epilepsie, u nichž se pomocí medikamentů nedaří dosá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áce je navrhnout takové analýzy signálů intrakraniálního2 EEG3 , které ve výsledku pomohou lékaři při rozhodování o typu a závažnosti onemocnění a volbě následné léčby. Motivace pro 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 s farmakorezistentní 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 diagnostikovaných pacientů. 1

padoucnice, padoucí nemoc, latinsky eufemicky také morbus sacer nebo morbus divinus vnitrolebečního 3 Nutno dodat, že intrakraniální EEG je jen jednou z řady diagnostických nástrojů epiletologa. 2

86

Jan Sova

2.

Motivace

Jako motivaci pro svou práci mohu uvést dva problémy, s kterými se setkávají epileptologové v klinické praxi [4] a které tak mohou být i konkretizovanými požadavky na výsledky analýzy EEG signálů. První problém můžeme označit jako tzv. pseudotemporální epilepsie. V případě pseudotemporální epilesie se jedná o možnou záměnu extratemporální epilepsie za temporální na základě vyšetření pomocí intrakraniálního (interiktálního a někdy i 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 epileptogenní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ých záchvatů vede k chybné diagnóze epilepsie temporálního5 laloku (TLE). MRI nález u postiž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 EEG vyšetření. Méně charakteristické jsou už výsledky dalších vyšetření, jako neuropsychologie, interiktálního PET a/nebo iktálního SPECT vyšetření. Ve skutečnosti mohou být interiktá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ích vě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ě ovlivň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 generována v různých částech rozsáhlého neuronálního okruhu v závislosti na chování jeho další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 synchronizač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: rozpoznávání typu záchvatů (parciální, generalizované, sekundárně generalizované), segmentaci 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 (viz dále). 4 5

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

Jan Sova

87

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 ρ21 ρ22 .. .. . . ρn1 ρn2

 · · · ρ1n · · · ρ2n   , . . . ..  .  · · · ρnn

(1)

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

E((X − µx )(Y − µy )) cov(i, j) = . σi σj σ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ým a 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žováno znaménko. Pro zjednodušení bude také dále uvažovano, že matice R neobsahuje nulové 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í vektor náležící vlastnímu číslu λ. Rovnice (8) se obvykle zapisuje jako homogenní soustava rovnic

88

Jan Sova

(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 čísel n kořenů. Rovnici A(λ) = 0 nazýváme charakteristická rovnice matice A a jejími koř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 vlastních vektorů pro případ kladné reálné symetrické matice (což je právě případ korelační matice R). Věta 1. Je-li A reálná symetrická matice řadu n, potom jsou všechny kořeny charakteristické rovnice reálné, tj. všechna vlastní čísla reálné symetrické matice jsou reálná. Věta 2. Dva vlastní vektory odpovídající různým vlastním číslům matice A jsou navzájem ortogonální. Věta 3. Ke každé symetrické matici A existuje ortonormální matice P taková, že PBPT = B je diagonální matice - prkvy bii na diagonále matice B jsou všechna vlastní čísla λi matice A (počítána i s jejich násobností) a sloupcové vektory matice P jsou jednotkové vzájemně ortogonální6 vlastní vektory matice A příslušné k vlastním číslům λi . Definujeme-li dále stopu matice tr(A) =

n X

aii ,

(8)

i=1

pak můžeme s využitím věty 37 ukázat, že pro symetrickou maitici A platí tr(A) = tr(PT BP) = tr(APPT ) = tr(B).

(9)

Protože tr(R) =

n X i=1

ρii =

n X

1 = tr(B) =

i=1

n X

λi ,

(10)

i=1

dostáváme n X

λi = n.

i=1

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

6

(11)

Jan Sova

89

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 Pnn tedy vždy dostaneme n reálných kladných vlastních čísel, přičemž bude platit i=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 λi P (λi ) = Pn

j=1

4.

Výsledky

λj

=

λi . n

(12)

Dále budu diskutovat výsledky aplikace výše popsané metody na signál s epileptickou aktivitou, který byl poskytnut z Kliniky dětské neurologie 2. LF UK FN Motol. Při výpočtu korelač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áchvatovou EEG aktivitu je korelační matice bez významných shluků lokální synchronizace nebo desynchronizace (viz obrázek 2a)), objevuji se prakticky bezprostředně po vzniku záchvatu mí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 mozku prakticky 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 je výkon rozdělen mezi přibližně 10% vlastních čísel, dochází během vzniku a šíření k přesunu vý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] na velké 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 diskutováno její předpokládané využití v epileptologii (včetně klasifikace záchvatu). Do budoucna bych 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 s další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í intrakraniálního EEG se většinou u pacientů s dobře lokalizovaným epileptogenním ložiskem a s negeneralizovanými (parciálními) záchvaty nedělá.

90

Jan Sova

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

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] Arhuis, M.; Valton, L.; Régis, J. Impaired consciousness during temporal lobe seizure 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 and subdural eeg studies. Epilepsia (2001), 491–503.

Jan Sova

91

[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-depend multuvariate 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 experimental seizure. 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 by analysing 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 emerging field. 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 fast intracerebral 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.

92

Adam Stráník

Klasifikace mezi /s/ a /š/ na základě parametrizace vstupního signálu pomocí LSF Adam Stráník České vysoké učení technické v Praze, Fakulta elektrotechnická [email protected] 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 o ostrý, jasně ohraničený vysokofrekvenční zvuk, který vzniká nevhodným přiblížením mluvidel při tvorbě /s/. Oproti tomu se často vyskytuje jev, kdy jsou sykavky tzv. tupé a mohou přecházet až do hlásky /š/. V příspěvku je představena 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 mluvidel. 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ých okrajů 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řitisknutím okrajů jazyka po stranách k horní dásni, hlavní zúžení je ovšem posunuté více dozadu, 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ři tvorbě 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 parametrizace signá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áci budou 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é jsou do sebe zasunuté. Pokud bude takovými válci proudit vzduch, bude docházet k různým rezonancím v závislosti na tom, zda bude tento model na konci otevřený nebo uzavřený

Adam Stráník

93

(otevření a uzavření hlasivek, úst atp.), přičemž počet rezonančních frekvencí závisí na počtu válců, kterými je hlasový trakt modelován, tedy na řádu modelu. Šířku a pozici těchto rezonančních frekvencí je možné popsat právě pomocí LSF, které přímo souvisí s LSP. 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 X

m=1

am x[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á je dá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 dostaneme až ke vztahu      R(0) R(1) . . . R(M − 1) a1 R(1)  R(1)     R(2) . . . R(M − 2)    a2  =  R(2)  , (3)      . . ... . . .  R(M − 1) R(M − 2) . . . R(0) aM R(M )

kde R je autokorelační funkce a am jsou váhy vzorků. Matici koeficientů am pak jednoduše dopočí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íky ně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 PARCOR (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 + a1 z + a2 z 2 + . . . + aM z M ,

(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 a mají k polynomu A(z) následující vztah A(z) = 1 2

P (z) + Q(z) , 2

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

(5)

94

Adam Stráník

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 ),

Q(z) = A(z) + z 3.3.

−(M +1)

(6)

−1

A(z ).

(7)

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čítat následovně   <Θk −1 ωk = tan . (8) =Θk

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 si odpoví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í minimum energie mezi nimi je (3. a 4. LSF).

2

10

2. LSF

1. LSF

6. LSF

5. LSF 1

H(f) [dB]

10

0

10

odhad LPC spektra

3. LSF

4. LSF

−1

10

500

1000

1500 2000 2500 → frequency [Hz]

3000

3500

4000

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

Adam Stráník

4.

95

Experimenty

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

s with stridence

Frequency

Frequency

1

x 10

2

2

1.5

1.5

2 1.5

4

s

x 10

Frequency

4

x 10

1

0

0.5

1 1.5 Time

0

2

(a)

1 0.5

0.5

0.5

sh

1

2 Time

0

3

(b)

0.2 0.4 0.6 0.8 Time

1

1.2

(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 /š/.

s with stridence

1

1

0

10

1

H(f) [dB]

10

H(f) [dB]

H(f) [dB]

sh

s

10

0

10

0

10

10 −1

10

−1

0.5 1 1.5 2 → frequency [Hz] x 104

0.5 1 1.5 2 → frequency [Hz] x 104

(a)

(b)

10

0.5 1 1.5 2 → frequency [Hz] x 104

(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 a b = i. Tento parametr jednoduchým způsobem popisuje, jak je významný nárůst energie ve frekvenční oblasti mezi těmito frekvenčními páry – čím je diference menší, tím je nárůst významnější a naopak. Pokud je diference příliš velká, může se jednat o lokální minimum energie.

96

Adam Stráník

4.2.

Průměr LSF – mLSFa,b

Parametr mLSFa,b je dán vztahem mLSFa,b =

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

(10)

kde ω je frekvence vypočtená ze vztahu (8) a i je přirozené kladné sudé číslo, a = i − 1 a b = i. Tento parametr zhruba odpovídá pozici lokálního maxima případně minima energie mezi LSP – 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) je dá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šších frekvencích. s with stridence

s

sh

s with stridence

silence

5000

4000

4000 LSF diffs [Hz]

LSF diffs [Hz]

5000

3000 2000

sh

silence

3000 2000 1000

1000 0 2000

s

3rd and 4th LSF

1st and 2nd LSF

3000

4000 5000 6000 LSF means [Hz]

7000

0 6000

8000

7000

8000

(a)

9000 10000 11000 12000 13000 LSF means [Hz]

(b) s with stridence

s

sh

silence

5th and 6th LSF 6000

LSF diffs [Hz]

5000 4000 3000 2000 1000 0 1.35

1.4

1.45

1.5 1.55 LSF means [Hz]

1.6

1.65

1.7 4

x 10

(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:

Adam Stráník

• • • •

97

/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]. Vstupní data byla rozdělena na 10 skupin se stejným poměrovým zastoupením každé třídy a hledané klasifikátory byly testovány pomocí crossvalidace. Výsledný klasifikátor je založený na rozhodovacím stromu, který je zobrazen na obr. 5. V tab. 1 je konfúzní matice tohoto 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 267 2 0 0

S Sh silence ← klasifikováno jako 27 0 0 Sstridence 384 0 0 S 1 1525 1 Sh 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 spektru u těchto hlásek. Detekce hvízdání je možná při popisu vstupního signálu LPC modelem 6. řádu, ze kterého jsou následně vypočteny diference LSF. Pokud je hvízdání v signálu pří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. 5 a jeho konfúzní matice je v tab. 1.

98

Adam Stráník

Poděkování Tato práce je podporována z: GACR102/08/H008 Analýza a modelování biomedicínský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 in fricatives 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. 448467. Dostupný také z WWW: <www.sciencedirect.com>. [4] BÄCKSTRÖM, Tom; MAGI, Carlo. Properties of line spectrum pair polynomials: 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 Explorations, 2009, Volume 11, Issue 1.

Daniel Špulák

99

Vliv průměrování na možnosti odhadu krevního tlaku z EKG a PPG Daniel Špulák České vysoké učení v Praze, Fakulta elektrotechnická [email protected] 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

100

Daniel Špulák

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. Datum

Místo

Označení

18.3.2009 23.4.2009 23.4.2009 23.4.2009 16.3.2010 16.3.2010 13.4.2010 13.4.2010 13.4.2010 13.4.2010

FNM FNM FNM FNM NH NH NH NH NH NH

2 4 5 6 16 17 25 26 27 28

Délka Rozsah tlaku záznamu/s 180 velký 224 velký 211 střední 210 malý 306 malý 319 střední 248 střední 304 malý 257 střední 345 malý

Tabulka 1: Soupis náměrů; FNM – Fakultní nemocnice Motol, NH – nemocnice Na Homolce 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. Korelační koeficienty (abs. h.): systolický tlak vs. čas do minima PPG Subjekt Počet R-R intervalů pro 2 16 17 25 26 27 28 průměrování 1 0,25 0,47 0,11 0,14 0,56 0,00 0,03 3 0,15 0,54 0,25 0,56 0,72 0,14 0,06 6 0,16 0,51 0,27 0,53 0,76 0,27 0,16 11 0,13 0,30 0,29 0,52 0,80 0,37 0,26 16 0,16 0,10 0,33 0,59 0,86 0,43 0,32

Tabulka 2: Absolutní hodnoty korelačních koeficientů mezi systolickým tlakem a PAT

Daniel Špulák

101

Korelační koeficienty (abs. h.): diastolický tlak vs. čas do minima PPG Subjekt Počet R-R intervalů pro 2 16 17 25 26 27 28 průměrování 1 0,20 0,45 0,08 0,29 0,53 0,28 0,00 3 0,01 0,56 0,16 0,53 0,69 0,07 0,02 6 0,13 0,57 0,15 0,53 0,72 0,12 0,13 11 0,34 0,45 0,14 0,51 0,76 0,17 0,23 16 0,46 0,35 0,17 0,57 0,81 0,22 0,28

Tabulka 3: Absolutní hodnoty korelačních koeficientů mezi diastolickým tlakem a PAT

Absolutní hodnota korelačního koeficientu Systolický tlak vs. čas do minima PPG

Korel. koef. (abs. hodn.)

1,00 0,80

2 16 17

0,60

25 26

0,40

27 28

0,20 0,00 1

3

6

11

16

Počet R-R intervalů pro průměrování

Graf 1: Závislost korelačních koeficientů na průměrování pro jednotlivé subjekty; systolický tlak

Absolutní hodnota korelačního koeficientu Diastolický tlak vs. čas do minima PPG

Korel. koef. (abs. hodn.)

1 0,8

2 16 17

0,6

25 26

0,4

27 28

0,2 0 1

3

6

11

16

Počet R-R intervalů pro průměrování

Graf 2: Závislost korelačních koeficientů na průměrování pro jednotlivé subjekty; diastolický tlak

102

Daniel Špulák

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.

Daniel Špulák

103

[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

104

Václav Turoň

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), HilbertHuangova 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 zexp(i 2πlt ) = zcos(2πlt ) + i zsin(2πlt ) l

l

m=0

m =0

= ∑ a 2′ m cos(2πmt ) + i ∑ b2′ m sin (2πmt ), l ∈ Z ,

(1)

kde a2′ m , b2′ m jsou koeficienty Zolotarevových polynomů a l je jejich řád. Báze FT může být zapsána pomocí komplexní exponenciály jako zexp(i 2πlt ) = cos(2πlt ) + i sin(2πlt ), l ∈ Z ,

(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.

Václav Turoň

105

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

 2πn  , n = 0K N − 1, N = 2000 s (n ) = cos f f S  

(3)

Frekvence signálu f je změněna z 600Hz na 3100Hz podle podmínky

 600, n ∈ 0, N / 2 f = . 3100, n ∈ N / 2 + 1, N − 1 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í

(4)

106

Václav Turoň

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.

Václav Turoň

107

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ů

108

Václav Turoň

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