Č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: <http://focus.ti.com/dsp/docs/dsphome.tsp?
sectionId=46&DCMP=TIHomeTracking&HQS=Other+OT+home_p_dsp>.
[TMS06]
Texas Instruments: TMS320C6000 Programmer's Guide [online]. 2006 [cit.
2010-09-07]. URL: <http://www.ti.com/litv/pdf/spru198i>.
[TMS10]
Allred Daniel: C6Run DSP Software Development Tool [online]. 2010 [cit.
2010-11-24]. URL: <http://www.ti.com/litv/pdf/spry143>.
[UHL07]
UHLÍŘ, Jan, et al. Technologie hlasových komunikací. Praha: Nakladatelství
ČVUT, 2007. 276 s. ISBN 978-80-01-03888-8.
[YOU09]
Young Steve. et al., The HTK Book. cambridge University Engineering
Department, 2006. [online] URL: <http://htk.eng.cam.ac.uk/ftp/software/htkbook.pdf.zip>.
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 <UNK> pro neznámá nebo neplatná slova,
• speciální token <S> resp. </S> pro začátek resp. konec věty,
• neobsahuje n-gramy jdoucí přez hranici věty, tzn. <S> se vyskytuje pouze na začátku a </S> pouze na konci n-gramu,
• čísla o více než 16 číslicích a slova delší než 32 znaků jsou nahrazena tokenem
<UNK>,
• slova se znaky z neevropských jazyků, neplatnými UTF8 znaky a kontrolními ASCII
znaky jsou také nahrazena tokenem <UNK>,
• slova vyskytující se méně než 40 krát jsou nahrazena tokenem <UNK>
• a všechny n-gramy které po předchozích operacích měly výskyt menší než 40 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 </s> označující hranici věty,
• obsahuje n-gramy jdoucí přes hranice vět
• a neobsahuje interpunkční znaménka.
Z nevhodných slov vyskytujících se v korpusu Web 1T 5-gram obsahuje tento 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 <UNK>:
– 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 </s> tokeny </s> 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ů </s> 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 <UNK>. 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
Download

analýza a zpracování řečových a biologických signálů sborník prací