ČASOVĚ-FREKVENČNÍ ANALÝZA ELEKTROGRAMŮ

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY

FAKULTA ELEKTROTECHNIKY A KOMUNIKAČNÍCH TECHNOLOGIÍ ÚSTAV BIOMEDICÍNSKÉHO INŽENÝRSTVÍ FACULTY OF ELECTRICAL ENGINEERING AND COMMUNICATION DEPARTMENT OF BIOMEDICAL ENGINEERING

ČASOVĚ-FREKVENČNÍ ANALÝZA ELEKTROGRAMŮ TIME-FREQUENCY ANALYSIS OF ELECTROGRAMS

DIPLOMOVÁ PRÁCE MASTER'S THESIS

AUTOR PRÁCE

Bc. PETR DOLEŽAL

AUTHOR

VEDOUCÍ PRÁCE SUPERVISOR

BRNO 2015

doc. Ing. JANA KOLÁŘOVÁ, Ph.D.

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ Fakulta elektrotechniky a komunikačních technologií Ústav biomedicínského inženýrství

Diplomová práce magisterský navazující studijní obor Biomedicínské inženýrství a bioinformatika Student: Ročník:

Bc. Petr Doležal 2

ID: 136466 Akademický rok: 2014/2015

NÁZEV TÉMATU:

Časově-frekvenční analýza elektrogramů POKYNY PRO VYPRACOVÁNÍ: 1) Seznamte se s elektrofyziologií srdce a základními studiemi z oblasti experimentální kardiologie. 2) Prostudujete problematiku snímání elektrické aktivity z izolovaných srdcí. Seznamte se databází signálů EG dostupné na UBMI. 3) Seznamte se s metodou časově-frekvenční analýzy využívající algoritmus Matching Pursuit a vytvořte základní programovou funkci. 4) Navrhněte programovou aplikaci pro práci s daty. 5) Analyzujte elektrogramy pomocí metody Matching Pursuit a spojité vlnové transformace a srovnejte výsledky. 6) Výsledky vhodně prezentujte a zhodnoťte možnosti použití metody Matching Pursuit. DOPORUČENÁ LITERATURA: [1] MALLAT, S. G. and ZHANG, Z. Matching Pursuits with Time-Frequency Dictionaries, IEEE Transactions on Signal Processing, December 1993, pp. 3397–3415. [2] SORNMO, L., LAGUNA, P. Bioelectrical Signal Processing in Cardiac and Neurological Applications. Elsevier Academic Press, 2005. Termín zadání:

9.2.2015

Termín odevzdání:

22.5.2015

Vedoucí práce: doc. Ing. Jana Kolářová, Ph.D. Konzultanti diplomové práce:

prof. Ing. Ivo Provazník, Ph.D. Předseda oborové rady UPOZORNĚNÍ: Autor diplomové práce nesmí při vytváření diplomové práce porušit autorská práva třetích osob, zejména nesmí zasahovat nedovoleným způsobem do cizích autorských práv osobnostních a musí si být plně vědom následků porušení ustanovení § 11 a následujících autorského zákona č. 121/2000 Sb., včetně možných trestněprávních důsledků vyplývajících z ustanovení části druhé, hlavy VI. díl 4 Trestního zákoníku č.40/2009 Sb.

Abstrakt Tato práce se zabývá časově-frekvenční analýzou elektrogramů měřených na izolovaných srdcích morčat perfundovaných podle Langendorffa. Časově-frekvenční analýza je založena na algoritmech Matching Pursuit a Wigner-Villovy distribuce. V teoretické části práce jsou popsány základy elektrokardiografie, měření na izolovaných srdcích, teorie aproximační metody Matching Pursuit a její kombinace s Wigner-Villovou distribucí zobrazující spektrum hustoty energie signálu. Také jsou představeny další běžné přístupy časově-frekvenční analýzy včetně teorie spojité vlnkové transformace. Uvedené algoritmy byly otestovány na souboru elektrogramů, u kterých byla v rámci měření vyvolána ischemie a následně reperfuze. Navržená metoda umožňuje rychlou detekci ischemie bez jakékoli apriorní znalosti signálu, a taktéž slouží jako nástroj k rozměření významných bodů a intervalů elektrogramu. V závěru práce byla prezentována úspěšnost metody a diskutovány možnosti jejího využití.

Klíčová slova Elektrofyziologie, EKG, izolované srdce podle Langendorffa, Matching Pursuit, Gaborovy atomy, Wigner-Villova distribuce, Spojitá vlnková transformace, časově-frekvenční analýza, detekce ischemie.

Abstract This thesis deals with time-frequency analysis of electrograms measured on isolated guinea pig hearts perfused according to Langendorff. Time-frequency analysis is based on algorithms Matching Pursuit and Wigner-Ville Distribution. The theoretical part describes the basics of electrocardiography, measurement on isolated hearts, the theory of approximation method Matching Pursuit and its combination with the Wigner-Ville distribution spectrum showing the energy density of the signal. Also other common approaches of time-frequency analysis are presented including the theory of continuous wavelet transform. The presented algorithms were tested on a set of electrograms, on which were induced ischemia within measurement followed by reperfusion. The proposed method allows for the fast detection of ischemia without any a priori knowledge of the signal, and also serves as a tool for measurement of EG important points and intervals. In the conclusion efficacy of the method was presented and its possible uses has been discussed.

Key words Electrophysiology, ECG, isolated heart perfused by Langendorff, Matching Pursuit, Gabor’s atoms, Wigner-Ville Distribution, Continuous Wavelet Transform, time-frequency analysis, ischemia detection. I

Bibliografická citace: DOLEŽAL, P. Časově-frekvenční analýza elektrogramů. Brno: Vysoké učení technické v Brně, Fakulta elektrotechniky a komunikačních technologií, 2015. 70 s. Vedoucí diplomové práce doc. Ing. Jana Kolářová, Ph.D.

Prohlášení Prohlašuji, že svou diplomovou práci na téma Časově-frekvenční analýza elektrogramů jsem vypracoval samostatně pod vedením vedoucího diplomové práce a s použitím odborné literatury a dalších informačních zdrojů, které jsou všechny citovány v práci a uvedeny v seznamu literatury na konci práce. Jako autor uvedené diplomové práce dále prohlašuji, že v souvislosti s vytvořením tohoto projektu jsem neporušil autorská práva třetích osob, zejména jsem nezasáhl nedovoleným způsobem do cizích autorských práv osobnostních a jsem si plně vědom následků porušení ustanovení § 11 a následujících autorského zákona č. 121/2000 Sb., včetně možných trestněprávních důsledků vyplývajících z ustanovení § 152 trestního zákona č. 140/1961 Sb.

V Brně dne ..............................

......................................... Petr Doležal

Poděkování Děkuji vedoucí mé diplomové práce doc. Ing. Janě Kolářové, Ph.D. za odborné vedení a cenné rady při zpracování mé práce. Také děkuji za cenné konzultace Ing. Martinu Vítkovi, Ph.D. a Ing. Lukáši Smitalovi, Ph.D. Velký dík patří též mé rodině a přítelkyni za neustálou podporu.

V Brně dne ..............................

......................................... Petr Doležal

II

Obsah: ÚVOD

1

1.

2

ELEKTROFYZIOLOGIE SRDCE

1.1.

Fyziologický popis srdce

2

1.2.

Převodní systém srdeční

3

1.3.

Vznik a šíření elektrického signálu v myokardu

3

1.4.

Měření EKG

5

1.5.

Časové a spektrální vlastnosti EKG signálu

6

1.6.

Patologická elektrofyziologie

8

1.7.

Vznik a projevy srdeční ischemie

9

1.8.

Detekce QRS a rozměření signálu EKG

2.

EXPERIMENTÁLNÍ KARDIOLOGIE

11 12

2.1.

Studium farmakologických látek

12

2.2.

Optické snímání elektrické aktivity srdce

12

3.

SNÍMÁNÍ ELEKTRICKÉ AKTIVITY Z IZOLOVANÝCH SRDCÍ

13

3.1.

Izolované srdce podle Langendorffa

13

3.2.

Snímání elektrogramů z izolovaných srdcí

14

3.3.

Vlastnosti EG signálů morčat

16

4.

DATABÁZE ELEKTROGRAMŮ

17

5.

ČASOVĚ-FREKVENČNÍ ANALÝZA SIGNÁLU

18

6.

VLNKOVÁ TRANSFORMACE

19

6.1.

Spojitá vlnková transformace (CWT)

19

III

6.3.

Mateřské vlnky

21

6.4.

Výběr vlnky

22

6.5.

CWT a analýza diskrétních signálů

22

7.

ALGORITMUS MATCHING PURSUIT A WIGNER-VILLOVY DISTRIBUCE

23

7.1.

Wigner-Ville Distribution

24

7.2.

Matching Pursuit

27

7.3.

Odhad hustoty energie signálu pomocí Matching Pursuit a Wigner-Villovy distribuce 31

8.

REALIZACE ALGORITMU

32

8.1.

Detekce QRS

34

8.2.

Vznik časově-frekvenční mapy

38

8.3.

Zobrazení detailů a energie signálu v časově-frekvenční mapě

40

8.4.

Srovnání s jinými metodami časově-frekvenční analýzy

44

8.5.

Automatická detekce srdeční ischemie

45

9.

DISKUSE VÝSLEDKŮ

53

ZÁVĚR

55

SEZNAM LITERATURY

57

SEZNAM ZKRATEK

60

PŘÍLOHY

61

IV

Seznam obrázků Obr 1 Srdce [2] ........................................................................................................................... 2 Obr 2 Převodní systém srdeční [2] ............................................................................................. 3 Obr 3 Akční potenciál myokardu [6] ......................................................................................... 4 Obr 4 Zobrazení unipolárních svodů (Einthovenův trojúhelník) [2] ......................................... 5 Obr 5 Významné vlny EKG křivky a rozměření typických intervalů [7] ................................. 7 Obr 6 Amplitudové spektrum, převzato z [8], upraveno ............................................................ 7 Obr 7 Srdeční ischemie [10]....................................................................................................... 9 Obr 8 Elevace a deprese ST segmentu [13] ............................................................................. 10 Obr 9 Schéma obecně platného postupu detekce QRS ............................................................ 11 Obr 10 Snímání akčního napětí na povrchu srdce pomocí optické metody [21] ..................... 12 Obr 11 Izolované srdce připojené ke kanyle [26] .................................................................... 13 Obr 12 Langendorffův perfuzní set [22] .................................................................................. 14 Obr 13 Zapojení při měření komorového tlaku [27] ................................................................ 15 Obr 14 Rozměření významných bodů EG zvířecího modelu [18] ........................................... 16 Obr 15 Vliv změny posunu a měřítka ...................................................................................... 20 Obr 16 Vlnka Mexican hat, Morletova vlnka, Meyerova vlnka a Haarova vlnka ................... 21 Obr 17 Časově-frekvenční mapa algoritmu Wigner-Ville Distribution ................................... 26 Obr 18 Gaborovy atomy [32] ................................................................................................... 28 Obr 19 Matching Pursuit aproximace signálu a použité Gaborovy atomy .............................. 29 Obr 20 Gaborovy atomy použité pro aproximaci MP signálu na Obr 19 ................................ 30 Obr 21 Výsledek metody Matching Pursuit společně s WVD ................................................. 31 Obr 22 Celý záznam jednoho EG (morče č.2, svod I) ............................................................. 32 Obr 23 Průměrný EG v každé minutě měření .......................................................................... 33 Obr 24 Srovnání průběhů stabilizace, ischemie a reperfuze .................................................... 33 Obr 25 Detekce R vln ............................................................................................................... 36 Obr 26 Počátečních 5 min ischemie ......................................................................................... 37 Obr 27 Průběh vzniku časově-frekvenční mapy ...................................................................... 38 Obr 28 Průběh vzniku časově-frekvenční mapy ...................................................................... 39 Obr 29 Průběh vzniku časově-frekvenční mapy ...................................................................... 39 Obr 30 Průběh vzniku časově-frekvenční mapy ...................................................................... 40 Obr 31 Zobrazení energie signálu v časově-frekvenčním měřítku .......................................... 41 Obr 32 Rozměření bodů P, Q a T na základě MPWVD-S ....................................................... 42

V

Obr 33 Porovnání časově-frekvenčních map (stabilizace, ischemie, reperfuze) ..................... 43 Obr 34 Srovnání metod Matching Pursuit, Spojité vnlnkové transformace a STFT ............... 44 Obr 35 Průběh vytváření průměrné časově-frekvenční mapy pomocí MPWVD-E................. 45 Obr 36 Průběh vytváření průměrné časově-frekvenční mapy pomocí CWT ........................... 45 Obr 37 Srovnání histogramů pro stabilizaci, ischemii a reperfuzi ........................................... 46 Obr 38 Detail změn frekvenčního zastoupení použitých atomů při MPWVD-E ..................... 47 Obr 39 Histogramy stabilizačních fází měření. Zavedení průměr. referenčního histogramu. . 48 Obr 40 Zjednodušený vývojový diagram metody automatické detekce ischemie ................... 49 Obr 41 Průběh detekce ischemie pomocí MPWVD................................................................. 50 Obr 42 Detekovaná ischemie ................................................................................................... 50 Obr 43 Histogramy CWT ( stabilizace, ischemie a reperfuze) ............................................... 51 Obr 44 Průběh detekce ischemie pomocí CWT ....................................................................... 51 Obr 45 Pozitivní detekce ischemie pomocí CWT .................................................................... 52

Seznam tabulek Tabulka 1 Přehled běžných hodnot významných intervalů a velikost vln u EKG morčete ..... 16 Tabulka 2 Přehledová tabulka účinnosti detektoru QRS ......................................................... 35 Tabulka 3 Přehled časů detekce ischemie metodou MPWVD-E a CWT ................................ 52

VI

Úvod Analýza elektrogramů je významnou oblastí experimentální medicíny, neboť je důležitým nástrojem diagnostiky srdečních onemocnění, které jsou hlavní příčinou mortality v rozvinutých zemích. S rostoucí potřebou rychlého vyhodnocení elektrokardiografických dat roste i význam jejich automatického rozměření či detekce významných abnormalit. Mezi ty patří především detekce srdeční ischemie, která se jakožto příčina infarktu myokardu podílí velkou měrou na úmrtnosti jak v České republice, tak ve světě. Při vzniku ischemického stavu rozhoduje především čas o rozsahu a závažnosti poškození tkáně, které může vést až k srdeční dysfunkci či akutnímu srdečnímu selhání. Proto má včasná detekce tak velký význam. Cílem této práce je na základě časově-frekvenční analýzy rozpoznat projevy ischemie uměle vyvolané na izolovaných srdcích a nalézt souvislosti prokazující výskyt ischemie dříve, než je to možné určit z časového průběhu elektrogramu. V první části práce je podrobně vysvětlena elektrofyziologie srdce, vznik a šíření elektrického impulsu myokardem a možnosti jeho záznamu. Elektrogram měřený ex vivo na zvířecím modelu je svým průběhem a vlastnostmi velmi věrně podobný elektrokardiogramu lidského jedince. Proto jsou nejprve popsány vlastnosti EKG křivky, uvedeny typické hodnoty a časy vztahující se k významným bodům tohoto signálu. Část práce je věnována možným patologickým stavům a jejich projevům na EKG křivce. Následně je prezentována oblast experimentální kardiologie se zaměřením na využití izolovaných srdcí perfundovaných podle Langendorffa. Je popsán průběh experimentu a princip záznamu elektrogramů a jsou uvedeny běžné hodnoty elektrogramu srdce morčete. Dále jsou představeny možnosti detekce QRS komplexu, která úzce souvisí s jakýmkoli hodnocením srdeční činnosti. Po představení časově-frekvenční analýzy v obecném pojetí je prezentována metoda spojité vlnkové transformace, a poté metoda Matching Pursuit společně s Wigner-Villovou distribucí. Právě metoda Matching Pursuit publikovaná poprvé Mallatem a Zhangem [1] je stěžejní metodou této práce a v praktické části je testována na souboru elektrogramů a srovnávána s vlnkovou transformací, která doposud vévodila časově-frekvenční analýze signálů. V rámci praktické části je nejprve realizována vlastní modifikace detektoru QRS. Následně je představen princip získání časově-frekvenční mapy. Dále jsou představeny dvě modifikace metody Matching Pursuit kombinované s Wigner-Villovou distribucí (MPWVD) a na měřených datech je prezentována jejich účinnost a využití. V závěru práce je zhodnocen přínos navržené formy časově-frekvenční analýzy a popsány její inovativní prvky, díky nimž může metoda nalézt široké praktické využití.

1

1.

Elektrofyziologie srdce

Srdce zajišťuje nepřetržité čerpání krve nezbytné pro fungování všech orgánů v těle. Funkční diagnostika srdce je jedním z nejzákladnějších vyšetření jedinců ve zdravém, emergentním i urgentním stavu. Základní analýza EKG signálu je dnes běžně součástí mnoha různých diagnostických postupů, neboť změna tvaru EKG signálu, jiná doba trvání některých intervalů atd. téměř vždy souvisí s poruchou elektrické aktivity srdce. Dnes dostává mapování srdeční činnosti nové formy vizualizace (jako např. intrakardiální 3D mapování).

1.1.

Fyziologický popis srdce

Srdce funguje jako pumpa s dvěma objemovými jednotkami, které pracují v sérii. Každá polovina sestává ze síně a komory oddělené chlopněmi. Do pravé síně přitéká neokysličená krev žílami, přes chlopně se dostává do pravé komory a odtud je vypuzována tlakem přes chlopně do malého krevního oběhu, do plic, odkud se okysličená vrací do srdce levou síní, a poté je z levé komory hnána do aorty a dále do velkého krevního oběhu do celého těla.

Obr 1 Srdce [2]

Činnost srdce je přísně regulována, a to jednak autoregulačními mechanismy, a také vegetativním nervovým systémem, tj. sympatikem a parasympatikem. Sympatický nervový systém moduluje činnost srdce prostřednictvím β1 adrenergních receptorů, parasympatický pomocí M2 cholinergních receptorů. Modulována je především rychlost geneze akčního potenciálu pacemakerových buněk sinoatriálního uzlu a vedení atrioventrikulární junkcí. Zjednodušeně můžeme říci, že sympatikus činnost srdce povzbuzuje a urychluje, zatímco parasympatikus tlumí. Činnost srdce je mechanická a elektrická, přičemž mezi těmito složkami existuje velmi těsný a důležitý vztah. Srdeční tep je normálně řízen autonomně, avšak při změně nároků na srdeční tempo je to (z větší části) na eferentních srdečních nervech. [3], [4], [5]

2

1.2.

Převodní systém srdeční

Část srdce specializovaná na převod elektrických impulsů se nazývá převodní systém srdeční. Ten pak spouští mechanickou práci srdce. Převodní systém srdeční tvoří sinoatriální uzel (SA), preferenční síňové dráhy, atrioventrikulární uzel (AV), Hisův svazek, pravé a levé Tawarovo raménko a Purkyňovy vlákna. Tento systém je tvořen dvěma druhy buněk, první typ najdeme v SA a AV uzlu a je specializován na vytváření spontánních vzruchů (zdroj srdeční automacie, pacemaker), a druhý typ buněk je specializován na rychlé vedení elektrických impulsů srdeční tkání. Další buňky myokardu jsou pracovní, tvoří zbytek síní a komor a jejich koordinovaným stahem je zajištěna mechanická práce srdce. [3]

Obr 2 Převodní systém srdeční [2]

1.3.

Vznik a šíření elektrického signálu v myokardu

Buňky síňového a komorového myokardu jsou navzájem propojeny pomocí skupinového spojení (gap junctions). To znamená, že podnět, který vznikl v síních či komorách, vyvolá úplnou kontrakci obou komor i síní, tzv. zákon „vše nebo nic“. Normálně podnět vzniká v SA uzlu na základě pacemakerového potenciálu. Tento buněčný potenciál (napětí) v SA uzlu nikdy není konstantní. Po repolarizaci, kdy dosahuje -70mV, nastává pomalu narůst depolarizace až do dosáhnutí prahového potenciálu (okolo -40mV), pak je znovu vyvolán akční potenciál. To vše se děje díky iontovým proudům, které vznikají změnami iontové vodivosti plazmatické membrány. Při -70mV nastává nahrnutí kationtů dovnitř buňky, což způsobuje onu pomalou depolarizaci.

3

Po dosažení hodnoty prahového potenciálu se zvýší vodivost pro vápenaté ionty, zvýší se proud Ca2+ dovnitř a nastane strmý nástup potenciálu odpovídající rychlému vzestupu depolarizace. Ve chvíli, kdy se potenciál dostane ze záporných hodnot do kladných, začne docházet k vytékání K+ ven, což opět repolarizuje pacemakerovou buňku zpět na -70mV, tedy na maximální diastolický potenciál. Jeden takový vzestup potenciálu myokardu odpovídá jednomu srdečnímu stahu, tepová frekvence je tedy přímo určena frekvencí impulzů pacemakeru. Standardně je pacemakerem SA uzel, stejnou schopnost spontánní depolarizace mají sice i další části převodního systému srdečního, avšak jejich vlastní frekvence je pomalejší než u SA uzlu. Slouží tedy jako záložní pacemakery v případě poruchy spontánní depolarizace v SA uzlu. [5]

Obr 3 Akční potenciál myokardu [6]

Klidný úsek mezi dvěma akčními potenciály buněčných membrán myokardu se nazývá plató. Fáze plató je velmi důležitá, neboť díky ní jsou části myokardu, které byly podrážděné jako první, v refrakterní fázi až do chvíle podráždění posledních částí myokardu, což zabraňuje vzniku krouživých vzruchů nazývaných reentry. Vlivem akčního potenciálu jsou otevřeny napěťově řízené kanály pro Ca2+, který díky tomu začne vtékat z extraceluárního prostoru do buňky. To má za následek vyšší koncentraci Ca2+ v buněčném cytosolu, což spustí (jako tzv. triggereffect) otevření ligandů řízených kanály sarkoplazmatického retikula (které slouží jako zásobárna Ca2+). Odtud začne Ca2+ vytékat do cytosolu, to vyvolá elektromechanické spřažení, tedy kontrakci srdce. [5]

4

1.4.

Měření EKG

Elektrokardiogram (EKG) zobrazuje potenciálové rozdíly mezi podrážděnými a nepodrážděnými částmi srdce. Vůbec nepodrážděný či naopak zcela podrážděný myokard nevykazuje na EKG žádnou změnu (žádná změna potenciálu). Při depolarizaci vznikají různé potenciály reprezentované jednotlivými vektory různé velikosti a směru, jejich sumací vznikne sumační (integrační) vektor. EKG zobrazuje právě okamžitý sumační vektor pomocí končetinových a hrudních svodů. Rozlišujeme dva typy svodových systémů: standartní 12 svodový systém a ortogonální svodové systémy. Mezi nejvíce rozšířené ortogonální systémy patří Frankův, McFeeův a SVEC III. Vytvářený obraz elektrické aktivity srdce vzniká pomocí průmětu elektrického vektoru do tří rovin ortogonálního systému. Běžně používaný je však systém 12 svodový, který skládá signál z 6 hrudních a 4 končetinových elektrod. Končetinové elektroda R je na pravé ruce, L na levé ruce, F na levé noze a referenční N na pravé noze. První tři elektrody nacházející se ve frontální rovině tvoří dohromady Einthovenův trojúhelník. Rozdílem potenciálů mezi těmito elektrodami získáme bipolární Eintovenovy končetinové svody I, II a III. Aritmetický průměr ze tří končetinových potenciálů tvoří referenční bod nazývaný Wilsonova svorka. Augmentové bipolární svody podle Goldbergera se vyjádří jako rozdíl potenciálu jedné elektrody vůči průměrnému potenciálu protilehlé strany. Značí se aVR, aVL a aVF. Unipolární hrudní svody V1-V6 jsou umístěné v rovině horizontální a nabízí detailnější pohled na srdce.

Obr 4 Zobrazení unipolárních svodů (Einthovenův trojúhelník) [2]

5

1.5.

Časové a spektrální vlastnosti EKG signálu

EKG křivka odráží důležité fáze činnosti myokardu. Tvar jednotlivých vln a kmitů spolu s časovými intervaly mezi nimi jsou nezbytným základem jakékoli další analýzy signálu. Vlna P reprezentuje depolarizaci síní. Repolarizace síní není vidět, neboť je překryta kmity Q, R a S, tvořící komplex QRS, které jsou projevem depolarizace komor. Vlna T značí repolarizaci komor. Kmitem R je označen první pozitivní kmit v komplexu QRS, kmity R tedy nejsou ve všech svodech synchronní. Vlna T má stejný směr jako kmit R, což je většinou pozitivní. Vlna P – Vzruch vychází z SA uzlu a depolarizace se šíří svalovinou pravé a levé síně. Vlna P má většinou hladký monofázický tvar, její velikost nepřekračuje 300 μV a netrvá déle jak 120ms. Často je těžké určit její přesný začátek a konec kvůli její malé velikosti. Interval PR – Jakmile dojde depolarizace do AV uzlu, dojde k náhlému zpomalení šíření vzruchu (AV uzel vede vzruch nejpomaleji z celého myokardu). Tento úsek je důležitý pro určení izoelektrické linie EKG signálu. Komplex QRS – Vlna Q odpovídá depolarizaci mezikomorového septa ve směru od levé komory k pravé, následuje depolarizace srdečního hrotu značená kmitem R a nakonec depolarizace stěn komor od endokardu k epikardu reprezentovaná vlnou S. Celý komplex QRS trvá 70-110 ms a je nejvýraznější částí EKG křivky díky výšce dosahující až 2-3 mV. Interval ST – Po depolarizaci komor je na chvíli elektrická aktivita srdce nulová, protože srdeční vlákna nacházející se ve fázi plató mají stejný elektrický náboj a nikde netečou žádné proudy. Vlna T – Po fázi plató nastává repolarizace komor probíhající od epikardu k endokardu. Vlna se nachází zpravidla 300 ms za QRS komplexem, ale její pozice je ovlivňována tepovou frekvencí (při vyšších tepových frekvencích je blíže komplexu QRS). Po vlně T někdy následuje vlna U, jejíž původ není zatím zcela objasněn. Interval RR – Vyjadřuje dobu trvání jednoho srdečního cyklu měřenou mezi dvěma po sobě následujícími R vlnami. RR interval je základním parametrem každé analýzy EKG. Interval PQ – Vyjadřuje dobu šíření vzruchu myokardem od SA uzlu do komor, měří se od počátku depolarizace síní po začátek depolarizace komor. Interval QT – Odpovídá době trvání celkové depolarizace a repolarizace komor. Délka trvání intervalu QT závisí na tepové frekvenci. Zkracování intervalu odpovídá vyšší tepové frekvenci. Prodlužování intervalu je spojeno s různými poruchami srdce s rizikem náhlé smrti. [5]

6

Obr 5 Významné vlny EKG křivky a rozměření typických intervalů [7]

Spektrem signálu dle Fourierovy transformace je signál rozložený na harmonické složky vyjádřený pomocí amplitudy a fáze. Tímto převodem z časové oblasti do oblasti frekvenční získáme zastoupení jednotlivých složek EKG signálu, díky čemuž můžeme signál snáze klasifikovat. K tomu nám slouží zpravidla amplitudové spektrum signálu, viz Obr 6.

Obr 6 Amplitudové spektrum, převzato z [8], upraveno

V běžném EKG nepřekračuje frekvenční rozsah vln P a T hodnotu 10Hz. Nejvyšších spektrálních hodnot nabývá QRS komplex s maximem mezi 10-20 Hz. Většina energie spektrální funkce je do 35-45 Hz.

7

1.6.

Patologická elektrofyziologie

Patofyziologie srdce může být různá, mohou nastat poruchy srdeční frekvence, srdečního rytmu či šíření vzruchu. Nastane-li hypoxie, typotermie či změny koncentrace K+ a Ca2+ iontů, mohou se stát zdrojem depolarizace i jiné buňky než obvykle a nastává porucha srdečního rytmu. Arytmie dělíme dle charakteristického klinického projevu na bradykardie a tachykardie, podle místa vzniku na supraventrikulární (oblast síní) a ventrikulární (oblast AV uzlu a komor), a podle patogenetického mechanismu rozeznáváme poruchu tvorby vzruchu, poruchu vedení vzruchu a poruchu šíření depolarizační vlny. Při poruše tvorby vzruchu v SA uzlu nastává sinusová brady/tachykardie, je změněna frekvence, ale rytmus zůstává pravidelný a jak komplex QRS tak i vlna P má normální tvar. Nebo může nastat nemoc chorého sinu (Sic Sinus syndrom), kdy se střídají fáze rychlé a pomalé srdeční akce. V případě, že se stane aktivním zdrojem vzruchů vedoucím k depolarizaci a kontrakcím myokardu jiné místo než SA uzel, dojde k předčasnému srdečnímu stahu, extrasystole. Podle místa předčasného vzruchu dělíme extrasystoly na síňové, junkční a komorové. [9] Spuštěná aktivita (triggered aktivity) je cyklický děj, kdy se myokard spontánně opakovaně depolarizuje (bez impulzu z převodního systému srdečního). Elektrofyziologickým projevem je tachyarytmie s typickým kolísáním amplitudy kmitu R (torsade de pointes). V případě kolísání mezi úplnou repolarizací a úplnou depolarizací nastává pozdní následná depolarizace, elektrofyziologickým projevem je tachyarytmie. Dalším projevem spuštěné aktivity je reentry. V myokardu vzniká místo, ze kterého se šíří depolarizace, po vlastní repolarizaci je znovu depolarizováno vzruchem, který se do tohoto místa vrátil díky aberantnímu šíření, a tak vzniká patologická smyčka vracejícího se vzruchu. Smyčka je rychlejší než původní pacemaker, proto jej nahradí a vznikají různé druhy tachyarytmií. [9] Fibrilace síní se vyznačuje chaotickým šířením depolarizace myokardem bez zřejmého pacemakeru. Rytmus je zcela nepravidelný, vlny P jsou nahrazeny vlnami f (nepravidelné vlnky s frekvencí až 600/min). Pro fibrilaci komor je charakteristický sled nepravidelných vlnek o frekvenci 250-500/min. Jedná se o bezprostředně život ohrožující stav, neboť nastává zástava oběhu, do 10 s ztráta vědomí, a po 4-5minutách je poškození mozku již ireverzibilní. [9] Flutter síní je tachyarytmie s vysokou frekvencí stahů srdečních síní, obvykle spojenou také s rychlou frekvencí stahů komor, což se projeví rychlým pulzem (pravidelný rychlý rytmus střídání depolarizace a repolarizace). Tachykardie je rychlý rytmus srdeční vyvolaný následným časným/pozdním spuštěním rytmu. Dalšími poruchami mohou být AV blok I, II III, blokáda levého a pravého Tawarova raménka, asystolie, bigeminie, couplet, triplet, komorová salva či fenomén R na T. Více informací nalezneme například v [9].

8

1.7.

Vznik a projevy srdeční ischemie

Mezi velmi častý a závažný patologický stav patří bezesporu srdeční ischemie. Při měření elektrogramů analyzovaných v této práci byla ischemie způsobována na izolovaných srdcích záměrně, aby mohly být sledovány změny v naměřeném signálu a testovány metody včasné detekce této srdeční poruchy. Ischemie je zapříčiněna nedostatečným prokrvením tkáně. To nastává nejčastěji v souvislosti s trombózou, kdy je ucpána jedna či více koronárních tepen. Důsledkem špatného či žádného prokrvení myokardu je hypoxie, tedy nedostatek kyslíku, který má za následek nedostatek živin a hromadění odpadních látek. Takto postižená tkáň přestává normálně fungovat a nekrotizuje. Jak nekróza, tak jizva je elektricky němá, což se projeví na EKG patologickými kmity Q. Tyto změny na EKG se dělí na morfologické a nemorfologické. Ischemická choroba srdeční má různá stadia, od akutního až po chronické.

Obr 7 Srdeční ischemie [10]

Morfologické změny Morfologické změny se hodnotí přímo z tvaru a rozměření EKG křivky a jejich hodnocení je velmi rozšířené v klinické praxi. Automatická detekce změn je však v tomto případě náročná a velmi náchylná na šum v signálu. K nejvýraznějším změnám dochází v segmentu ST odpovídající zóně poškození a vlně T vypovídající o samotné zóně ischemie. Někdy dochází i k výrazným změnám v QRS, například nastane-li transmurální infarkt (procházející celou stěnou díky trombu v koronární tepně). Transmurální infarkt je spjat s nálezem patologického kmitu Q, jinak se jedná o netransmurální infarkt myokardu. [11],[12]

9

T vlna Vlna T je nejcitlivější částí EKG cyklu kvůli mimořádným nárokům na množství ATP, potřebném pro repolarizaci myokardu. Na vlně T je ischemie dobře viditelná, projevy však mohou být různé. Vlna T může být zvýšená, invertovaná, roztažená/zploštělá či bifazická. Avšak inverze vlny T ve svodech III, aVR a V1 může být zcela fyziologická. Vysoký hrotnatý tvar T vlny na hrudních svodech je jedním z příznaků akutního infarktu myokardu (IM). Vysoké T vlny ve svodech V1 – V3 jsou zpravidla důsledkem ischemie zadní stěny levé komory. Normální velikost vlny T by měla být do 1mV a v poměru k vlně R dosahovat 1/8 až 2/3 vlny R. ST segment Nejnápadnějším změnou ST segmentu je jeho deprese či elevace. Deprese ST je způsobena postižením subendokardu a elevace ST postižením subepikardu. Deprese ST segmentu projevující se při hypertrofii levé komory je nejvíce viditelná ve svodech V4 – V6. Horizontální deprese ST je známkou srdeční nedostatečnosti. Elevace ST je patrná při vazospastické angině pectoris. Elevaci ST segmentu často provází právě inverze vlny T.

Obr 8 Elevace a deprese ST segmentu [13]

Nemorfologické změny Pro nalezení nemorfologických změn není potřeba provádět rozměření EKG signálu, které bývá náročné a složité. Proto jsou takovéto změny vhodnější pro automatickou detekci ischemie. Pro získání nemorfologických parametrů můžeme signál zpracovat pomocí vlnkové transformace, analýzy hlavních komponent, analýzou srdečního rytmu či jen ve frekvenční oblasti pomocí výkonového spektra signálu. Výstupem takovýchto metod je soubor příznaků, který je rozdílný pro ischemický a neischemický EKG signál. [11], [12]

10

1.8.

Detekce QRS a rozměření signálu EKG

Komplex QRS je nejvýraznějším útvarem v EKG signálu, proto jej lze také nejsnáze detekovat, a je tak vhodný jako orientační bod k rozměření signálu EKG na dílčí úseky. Již ze samotné frekvence QRS lze odvodit tepovou frekvenci srdce či variabilitu srdeční činnosti. Určení polohy QRS komplexu usnadní určení méně výrazných vln P a T. Detekce QRS je základním stavebním kamenem jakékoli analýzy EKG. [14] Při detekci QRS je kladen důraz především na co nejvyšší účinnost. Ta je ztížena především variabilitou EKG signálu, přítomností artefaktů a šumu, nebo stimulačními impulzy kardiostimulátoru. Mezi běžné přístupy detekce můžeme řadit zejména algoritmy využívající první derivace či digitální filtraci, algoritmy založené na vlnkové transformaci, neuronových sítích či genetických algoritmech. Dále například detekci pomocí skrytých Markovových řetězců, přizpůsobené filtraci, průchodu nulou a další. [15] Jednotlivé přístupy zde popsány nebudou, ale můžeme uvést to, co mají společné, a sice předzpracování signálu a rozhodovací algoritmus (viz Obr 9). Při předzpracování se snažíme co nejvíce zvýraznit QRS komplex a zároveň potlačit nežádoucí složky šumu a různých artefaktů. Rozhodovací algoritmus pak porovnává transformovaný signál s prahovou hodnotou a dle toho indikuje přítomnost QRS.

Obr 9 Schéma obecně platného postupu detekce QRS

11

2.

Experimentální kardiologie

Laboratoře experimentální kardiologie se zabývají fyziologií a patofyziologií srdeční činnosti. Na různých experimentálních zvířecích modelech je analyzována srdeční činnost ve více úrovních biologické komplexity – od molekulární úrovně, měření membránových iontových proudů, po měření (na srdcích) in/ex vivo. Elektrická aktivita srdce in vivo je označována jako elektrokardiogram (EKG), zatímco elektrická aktivita izolovaného srdce je označována jako elektrogram (EG).

2.1.

Studium farmakologických látek

Izolovaná perfundovaná srdce slouží ve velké míře ke studiu a testování farmakologických látek. Hodnotí se například funkční a morfologické změny, aktivace kináz, reakce na vystavení různým farmakům a jejich vliv na post-ischemickou fázi. Zkoumána může být také například změna růstového faktoru při vystavení záření způsobující akutní radiační poškození [20]. Můžeme také hodnotit afinitu konkrétního ligandu k nějakému receptoru, aj.

2.2.

Optické snímání elektrické aktivity srdce

Optické snímání elektrické aktivity srdce je založeno na využití fluorescenčního jevu, který se projevuje po excitaci navázaných molekul na membránu srdečních buněk napěťově citlivých barviv. Napěťově závislá fluorescence se odvíjí od intra- a extracelulárních změn v závislosti na napěťovém gradientu. Barviva můžeme dělit na rychlá a pomalá dle jejich reakčního času. V elektrofyziologii se využívají rychlá barviva s reakčním časem na úrovni mikrosekund. Tato reakční doba je dostatečná pro sledování změn depolarizace buněk, která je řádově v milisekundách. Princip lze zjednodušeně vyjádřit tak, že po excitaci při kladné změně potenciálu nastává posun emisního spektra barviva směrem k vyšším vlnovým délkám. Výsledný snímaný signál je úměrný intenzitě emisního světla, který odpovídá průběhu akčního napětí. [21]

Obr 10 Snímání akčního napětí na povrchu srdce pomocí optické metody [21]

Jednou z nejrozšířenějších oblastí experimentální kardiologie je snímání elektrické aktivity izolovaných srdcí. Této problematice je věnována následující kapitola. 12

3.

Snímání elektrické aktivity z izolovaných srdcí

Izolované perfundované srdce je funkční preparát srdce ex vivo. Na aortu vystupující ze srdce je napojena kanyla, přes níž je srdce prostřednictvím koronárního řečiště vyživováno. Existuje několik variant napojení izolovaného srdce. Měření a analýzy na izolovaném srdci mají široké spektrum využití a patří mezi nejrozšířenější techniky experimentální kardiologie.

3.1.

Izolované srdce podle Langendorffa

Metoda perfuze izolovaného savčího srdce byla zavedena v roce 1895 prof. Oscarem Langendorffem [24]. Metoda izolovaného srdce umožňuje měření pro fyziologii, farmakologii, toxikologii. Výběr zvířete je prováděn s maximální snahou získat podobné výsledky jako u člověka, nejčastěji se používají srdce králíků, morčat a potkanů. Při práci s laboratorními zvířaty jsou vždy dodržována všechna etická pravidla a koncepce 3R (Reduction – nejméně zvířat, Refinement – nejméně stresu a bolesti, Replacement – nahrazení modely). Provede se celková anestezie zvířete, následuje provedení sternotomie, oddělení srdce od velkých cév a vyjmutí srdce. Kvůli rychlému nástupu ischemie je nutné, aby tato operace zabrala co nejméně času. Srdce se nejprve vloží do studeného roztoku (ledová tříšť fyziologického roztoku NaCl), aby se nejvíce zpomalil srdeční metabolismus. Srdce zbavíme perikardu a připevníme jej aortou ke kanyle perfuzního přístroje a ponoříme do termostaticky vyhřívané lázně na 37°C. Lázeň je tvořena vyživovacím roztokem shodným s perfuzátem. Fyziologickým roztokem ještě nejprve propláchneme levou komoru z důvodu případné krevní sraženiny. Vyhřívání perfuzního systému je umožněno cirkulující destilovanou vodou ohřívanou pomocí termostatu. Kanylou je do srdce přiváděn perfuzát, který srdci zajišťuje přísun živin a kyslíku do myokardu. Díky perfuzátu srdce po chvíli samo obnoví svou činnost. Perfuzní roztok je do srdce hnán pod takovým tlakem (70-80 mmHg dle hodnoty systolického tlaku), aby došlo k uzavření chlopní retrográdním tokem roztoku. Perfuzát je veden do koronárního systému a srdeční dutiny jsou v podstatě v průběhu experimentu prázdné. Koronární průtok je přibližně 6-8 ml/min. Jako perfuzát je nejčastěji používán KrebsůvHenseleitův roztok (118mM NaCl, 24mM NaHCO3, 4.2mM KCl, 1.2 mM KH2PO4, 1.2 mM MgCl2, 1.2 mM CaCl2, 5.5mM glukóza a 10mM Taurin). Roztok je provzdušňován směsí 95% O2 a 5% CO2. [25]

Obr 11 Izolované srdce připojené ke kanyle [26]

13

3.2.

Snímání elektrogramů z izolovaných srdcí

Elektrickými projevy izolovaného srdce jsou bifázické akční potenciály, tzv elektrogramy (EG). Elektrody snímacího zařízení (zpravidla Ag-AgCl) jsou umístěny v lázni, kde je ponořeno srdce. Elektrické signály jsou zesilovačem až tisícinásobně zesíleny, filtrovány anti-aliasingovým filtrem zabraňujícím překryvu spekter, a následně digitalizovány. Elektrody se zapojí do svodů, které přibližně odpovídají končetinovým svodům EKG u člověka (pravá síň, levá síň a levá komora). Takto napojené srdce je schopné pracovat podobně jako srdce živé díky schopnosti srdeční automacie (vlastní srdeční činnosti). Parametry digitalizace jsou určeny vlastnostmi signálu a požadovanými parametry pro zpracování a analýzu signálu. Hlavními vlastnostmi signálu jsou dynamický a frekvenční rozsah (0-100 Hz), vzorkovací frekvence (250, 500, 1000, 2000 Hz), bitová hloubka (8, 12 či 16 bitů). [21]

Obr 12 Langendorffův perfuzní set [22] (A – dvouplášťová nádoba, B,C – snímání teploty, D – okysličovaný roztok, E – kanyla s odstraňovačem bublin, F – systém tlakové kontroly, G – přepadový odtok, H – stříbro-chloridové elektrody)

14

Perfuze začíná v Langendorffově aparátu s průtokovou rychlostí 6-8 ml/min, aby byl umožněn bazální perfuzní tlak 70 mmHg. Rychlost průtoku během experimentu kolísá vlivem teploty perfuzátu, dle tepové frekvence a též díky výskytu extrasystol. Koronární perfuzní tlak je měřen skrz laterální spojení v perfuzní kanyle a levo-ventrikulární tlak je měřen s použitím latexového balónku nafouknutého na diastolický tlak 5-10 mmHg. Měření tlaku bývá prováděno například pomocí tlakových převodníků Statham, viz Obr 13. [27]

Obr 13 Zapojení při měření komorového tlaku [27]

Během měření EG analyzovaných v rámci této práce byl nejprve měřen na EG normální stav po dobu 30 min. Poté byl přívod perfuzátu zastaven a následovala fáze ischemie trvající 15 min. Poté byl přítok perfuzátu obnoven a následovalo 15 min reperfuze, během níž se srdce znovu zásobuje kyslíkem a ionty.

15

3.3.

Vlastnosti EG signálů morčat

Záznam elektrické aktivity srdce morčete se zásadně od lidského nikterak neliší. Nejprve se objevuje vlna P, následovaná komplexem QRS a poté vlnou T. Přece jen však zde nepatrné rozdíly jsou. Už jenom díky rozdílu mezi velikostí malého srdce morčete a oproti tomu podstatně většího srdce lidského lze očekávat kratší doby trvání jednotlivých srdečních časových úseků. Kratší časové průběhy jsou dány šířením akčního napětí po kratší dráze. Změny lze pozorovat též ve velikosti výchylky jednotlivých vln a také ve spektru signálu. Tepová frekvence se pohybuje v rozmezí 150-300 tepů za minutu. Zaznamenané experimentální signály elektrické aktivity srdce morčat tak nemůžeme považovat za ryze fyziologické. Taktéž proto, že samotné srdce je izolováno a na zaznamenávané elektrické aktivitě srdce se tak nepodílí celý kardiovaskulární systém, nehovoříme o elektrokardiografu (EKG), ale elektrografu (EG) [16], [17]. Názorné rozměření můžeme vidět na Obr 14. Tabulka 1 Přehled běžných hodnot významných intervalů a velikost vln u EKG morčete [16], [18]

Parametr

doba trvání, výchylka

P vlna, délka trvání

15 - 30 ms

P vlna, amplituda

0,01 mV

PQ interval

25 - 50 ms

RR interval

180 ms

PR interval

60 ms

QRS komplex

8 - 20 ms 1,1 – 1,9 mV

R vlna ST segment

55 ms

QT

110 ms

T vlna

60 ms

Obr 14 Rozměření významných bodů EG zvířecího modelu [18]

16

4.

Databáze elektrogramů

Data analyzovaná v rámci této práce pocházejí z databáze signálu Ústavu biomedicínského inženýrství (UBMI). Jedná se o sérii EG signálů morčat měřených na pracovišti fyziologického ústavu lékařské fakulty MU. Měřeno bylo sedm sérií orthogonálních EG signálů (EG I, EG II – horizontální, EG III - vertikální) na izolovaném srdci podle Landendorffa. Signály jsou zesíleny pomocí předzesilovače DAM-50 (World Precision Instruments, Inc.) a digitalizovány sběrnou kartou s 16bitovým dynamickým rozsahem a vzorkovací frekvencí 2000 Hz. Signály jsou získány za pomocí LabView kompatibilní datově akviziční multifunkční karty PCI-6250 (National Instruments, USA). Pro všechna měřená data je společný průběh měření: přibližně 30 min stabilizace, následně 15 min průtokové ischemie a 15 min reperfuze. Perfuzní tlak byl vždy 80 mmHg, teplota perfuzátu 37°C. Všechna měřená morčata byly samičky, hmotnost morčat se pohybovala v rozmezí 375-405g, hmotnost samotného srdce pak v rozmezí 1,6 – 2,6g. Anestezie byla provedena pomocí éteru a cervikální dislokace. V následující části práce bude představeno téma časově-frekvenční analýzy a jejich konkrétních metod, které byly aplikovány na měřenou databázi elektrogramů. Analýza se zaměří především na projevy ischemie v EG signálu a na za základě srovnání časověfrekvenčních spekter se snaží o automatickou detekci ischemie. Měřená data tak poslouží k testování nové metody, která by mohla najít své uplatnění v klinické praxi.

17

5.

Časově-frekvenční analýza signálu

Převážná většina reálných technických signálů patří do oblasti nestacionárních signálů. Fourierova transformace a její modifikace jsou techniky vhodné především ke zpracování stacionárních signálů. Můžeme je využít pro analýzu nestacionárních signálů v případě, že nás zajímají frekvenční komponenty obsažené v celém signálu, ale nedávají nám přehled o časovém výskytu frekvenčních složek. Proto pro určení časové lokalizace frekvenčních komponent je nutné použít časově-frekvenčních transformací. Metody časově-frekvenční analýzy dělíme na lineární (krátkodobá Fourierova transformace, vlnkové transformace) a nelineární (Cohenovy, Afinní, Hyperbolické a další). Velkou výhodou lineárních transformací je především rychlost výpočtu a poměrně uspokojivé časově-frekvenční rozlišení. Hlavní nevýhodou je však limitace lineárních transformací tzv. Heisenbergovým principem neurčitosti, který říká, že nelze přesně určit, jaká frekvence se bude vyskytovat v daném časovém okamžiku. To má za následek, že nemůžeme přesně určit lokalizaci dané složky signálu v časově-frekvenčním prostoru, ale je možné určit pouze pozici uvnitř obdélníka Δt ⋅ Δω v dané časově-frekvenční oblasti (Δt je minimální časový interval a Δω je minimální frekvenční interval). [28] Časově-frekvenční analýza podává informaci o změně frekvence v čase. Nalézt však přístup umožňující přesnou časovou lokalizaci dané části spektra při dostatečném časovém a frekvenčním rozlišení není tak jednoduché. Mezi základní časově-frekvenční zobrazení patří spektrogram spočítaný pomocí krátkodobé Fourierovy transformace (STFT – Short-Time Fourier Transform). Ten však nepřináší pro podrobnou analýzu dostatečné rozlišení kvůli omezenosti pevnou velikostí okna. Metodou vycházející z STFT je vlnková transformace (WT, z angl. Wavelet Transform). Ta umožňuje změnu okna vedoucí k optimálního poměru rozlišení v čase a frekvenci. Základní model vlnkové transformace bude představen v následující kapitole. Další metodou časově-frekvenční analýzy je Wigner-Villova transformace, která dosahuje velmi dobrých výsledků. Objevují se v ní sice artefakty interference, ale ty lze úspěšně potlačit využitím metody Matching Pursuit. Metoda Matching Pursuit kombinovaná WignerVillovou transformací je základem této práce. Její teoretické uvedení i praktické provedení je představeno v další části práce, spolu s porovnáním vůči metodě spojité vlnkové transformace.

18

6.

Vlnková transformace

Vlnková transformace (WT – z angl. Wavelet Transform) nachází uplatnění zejména pro analýzu signálů nestacionárních, aperiodických či zašumělých. Transformace je založena na lineární okénkové operaci se dvěma parametry, měřítkem a posunem. Transformaci si lze představit jako korelaci signálu s bázovými funkcemi (vlnkami). Bázové vlnky jsou odvozeny od vlnky mateřské. Právě změnou měřítka (dilatací) mateřské vlnky a jejím posunem získáme časově-měřítkovou reprezentaci, která představuje lokalizaci různých frekvenčních pásem v čase. Jde tedy v podstatě o časově-frekvenční analýzu (v různých rozlišovacích oblastech), která se dá využít ke sledování nespojitosti, špiček a trendů v signálu. Pokud budeme měnit měřítko a posun spojitě, jedná se o takzvanou spojitou vlnkovou transformaci (CWT – z angl. Continuous Wavelet Transform). Pokud budeme tyto parametry měnit s nějakým diskrétním krokem, jedná se o diskrétní vlnkovou transformaci (DWT – z angl. Discrete Wavelet Transform) a pokud máme signál diskrétní v čase, jedná se o takzvanou vlnkovou transformaci s diskrétním časem (DTWT – z angl. Discrete-Time Wavelet transform). Vzhledem k zadání této práce zde bude představena pouze CWT, která je následně aplikována na reálná data v rámci praktické části a její výsledky jsou porovnávány s Matching Pursuit, primární metodou této práce. [28], [30]

6.1.

Spojitá vlnková transformace (CWT)

Spojitá vlnková transformace je metoda podobná krátkodobé Fourierově transformaci (STFT – z angl. Short-Time Fourier Transform). U obou metod se provádí výpočet vždy pro danou část signálu vybranou okénkem. Okno je však u STFT vždy stejné, což má za následek problém s rozlišením ve frekvenční a časové oblasti. Spojitá vlnková transformace používá proměnné okénko dle frekvence. Díky tomu lze změnou okna (jak šířky, tak tvaru) dosáhnout optimálního poměru rozlišení v čase a frekvenci. To umožňuje pro vyšší frekvence dobré časové rozlišení a horší frekvenční rozlišení, a u nižších frekvencí dobré frekvenční rozlišení a horší časové. CWT tedy umožňuje analýzu s násobným rozlišením (označení MRA – z angl. Multi Resolution Analysis). Taková vlastnost je užitečná pro analýzu signálů, které obsahují vysokofrekvenční komponenty po krátkou dobu a zároveň nízkofrekvenční komponenty po dlouhou dobu. Takových signálů je však veliké množství a lze mezi ně zařadit i signály EKG. Spojitá vlnková transformace je definována vztahem 𝑊𝑎 𝑥(𝑏) =

1 √𝑎

+∞

∫ −∞

𝑡−𝑏 𝑥(𝑡)Ψ∗ ( ) 𝑑𝑡, 𝑎

𝑎>0

kde x(t) je transformovaný signál, funkce je popsána dvěma spojitými parametry: a (měřítko – dilatace vlnky) a b (posunutí vlnky po časové ose). Dilatace určuje frekvenční spektrum 19

příslušné vlnky (konstanta normalizuje energii jednotlivých vlnek). Časový posun umožňuje postupně pokrýt vlnkami určitého konečného trvání celý časový rozsah signálu. Dohromady tvoří tzv. měřítkovou funkci (scaling function). Funkce Ψ∗ označuje základní neboli mateřskou vlnku (mother wavelet). Symbol * značí komplexně sdruženou funkci.

Transformace je tedy v podstatě časově-frekvenční rozklad na základě korelace signálu 𝑥(𝑡) s bázovými vlnkami. Čím větší je parametr 𝑎, tím více se zúží spektrum vlnky a dojde k posunu k nižším frekvencím. Stlačením vlnky tedy získáme zvýšení časové rozlišitelnosti a snížení frekvenční rozlišitelnosti, a naopak. Násobení konstantou

1 √𝑎

slouží k normalizaci

energie vlnky při změnách měřítka. [28]

Obr 15 Vliv změny posunu a měřítka

6.2.

Vlastnosti CWT

Linearita CWT vyplývá přímo z definičního vztahu. Invariance v čase znamená, že posun analyzované funkce po časové ose způsobí stejný posun vlnkových koeficientů po ose polohy. Lze to odvodit z představy, že CWT může být interpretována pomocí množiny lineárních časově-invariantních filtrů. Diskrétní WT vlastnost invariance v čase nemá (díky decimacím), což představuje velkou nevýhodu. Lokalizační vlastnost říká, že odezvou na Diracův impuls v 𝑡0 je vlnka se středem v 𝑡0 . Dilatace popisuje závislost mezi CWT originální funkce a její zúženou či roztaženou podobou. Ve vlnkových koeficientech dochází k roztažení v ose polohy a k posunu v ose měřítka.

20

Mateřské vlnky

6.3.

Mateřská vlnka musí mít nulovou střední hodnotu a musí být nenulová jen na konečném časovém intervalu (nebo se zanedbatelnými hodnotami mimo tento interval). Roztažení mateřské vlnky vede k zúžení spektra a posunu spektra k nižším kmitočtům. Následující přehled uvádí některé rozšířené vlnky, které se používají pro CWT. [31] Vlnka Mexican hat Vlnka Mexican hat (Mexický klobouk, Marrova vlnka) má tvar druhé derivace průběhu hustoty pravděpodobnosti Gaussova rozdělení. Je symetrická, nemá kompaktní nosič a není ortogonální (proto nelze použít pro DWT). Vhodná pro CWT. Morletova vlnka Morletova vlnka má tvar komplexní sinusovky modulované Gaussovským oknem. Nabízí kompromis mezi polohovou lokalizací jednorázových dějů (Mexican hat) a frekvenčním rozlišením (Fourier transform). Vlnka je symetrická, komplexní, nemá kompaktní nosič, není ortogonální a je vhodná pro CWT. Meyerova vlnka Meyerova vlnka je definována ve frekvenční oblasti, v časové oblasti ji nelze vzorcem explicitně vyjádřit. V originálním tvaru ji nelze realizovat FIR filtry, proto se používá její diskrétní aproximace (konvolučními) filtry. Je symetrická, nemá kompaktní nosič, ortogonální, vhodná jak pro CWT, tak DWT. Haarova vlnka Haarova vlnka má velmi jednoduchý „obdélníkový“ tvar, díky tomu však neumožňuje hladkou rekonstrukci signálu. Vlnka je symetrická, má kompaktní nosič, je ortogonální, vhodná jak pro CWT, tak DWT, a díky své jednoduchosti je jednoduchá a efektivní i její implementace. Velkou nevýhodou pro aplikovatelnost vlnky je její nespojitost. Vlnky Daubechies Vlnky Daubechies jsou vlnky různé řádu 𝑁 ≥ 1, nemají explicitní vyjádření, jsou asymetrické, kompaktní nosič, vhodné pro CWT i DWT, ortogonální. vlnka Mexican Hat

Morletova vlnka

1

Haarova vlnka

Meyerova vlnka

1

1.5

1.5

0.8

0.8

0.6

0.6

1

1

0.4

amplituda

0.2

amplituda

0.5

0.2

0.4

amplituda

amplituda

0.5

0

0

0

-0.2

-0.5

-0.4

0

-0.5

-0.6

-1

-0.2 -0.8

-0.4 -8

-6

-4

-2

0 cas [s]

2

4

6

8

-1 -8

-6

-4

-2

0 cas [s]

2

4

6

-1 -8

8

-6

-4

-2

0 cas [s]

2

4

6

8

-1.5 -0.5

Obr 16 Vlnka Mexican hat, Morletova vlnka, Meyerova vlnka a Haarova vlnka

21

0

0.5 cas [s]

1

1.5

6.4.

Výběr vlnky

Ze studií zabývajících se aplikací vlnkové transformace vyplývá, že nejčastěji je výběr vhodné vlnky prováděn heuristicky nebo intuitivně. Lze však popsat několik souvislostí mezi charakterem analyzovaného signálu a vhodnou vlnkou: 

Komplexní vlnky detekují dobře oscilace a špatně detekují osamocené singularity. Také se hodí k současné detekci amplitudy a fáze.



Reálné vlnky s málo oscilacemi detekují dobře špičky a osamocené singularity.



Antisymetrické vlnky detekují dobře změny gradientu.



Symetrické vlnky nezpůsobují fázový posun mezi špičkou v signálu a projevem vlnkových koeficientů.

6.5.

CWT a analýza diskrétních signálů

Jelikož v praxi pracujeme nejčastěji s již diskrétními signály, nepočítáme spojitou WT, ale realizuje se aproximace CWT pro diskrétní signály. Před vlastní analýzou pomocí vlnkové transformace je potřeba vybrat nejprve typ transformace. Můžeme volit mezi CWT, paketovou DTWT, dyadickou DTWT a dalšími. V rámci této práce je dle zadání volena CWT a její srovnání s metodou Matching Pursuit, na níž je práce zaměřena primárně. Dále je potřeba před analýzou zvolit typ mateřské vlnky a její parametry, vhodná měřítka frekvenčních pásem pro dobrou viditelnost detailů, a způsob zobrazení, který může být například vrstevnicový či konturový. Realizace CWT v rámci této práce byla provedena pomocí funkce cwt.m zahrnuté ve Wavelet Toolbox. Bylo testováno několik vlnek, nejvíce se osvědčilo použití vlnky Meyerovy. Časově-frekvenční analýza EKG signálu pomocí CWT byla provedena kvůli porovnání výsledků a schopnosti frekvenčního a časového rozlišení vůči metodě Matching Pursuit.

22

7.

Algoritmus Matching Pursuit a Wigner-Villovy distribuce

Pro analýzu signálu metodou Matching Pursuit je nezbytný soubor známých funkcí, který nazýváme slovníkem. Abychom určili, které z funkcí slovníku nejlépe vyjadřují zkoumaný signál, potřebujeme nějak určit a kvantifikovat podobnost dvou signálů. Při hodnocení shodnosti hraje hlavní roli skalární součin, ortogonalita, frekvence a fáze. Nejrozšířenější metodou měření shodnosti signálů je jejich skalární součin. Princip určení hodnoty skalárního součinu je prostý, hodnoty jednotlivých bodů prvního signálu vynásobíme korespondujícími hodnotami signálu druhého, poté získané násobky sečteme, a získáme tak jediné číslo, jehož hodnota vypovídá o shodnosti signálů. Čím vyšší je tato hodnota, tím jsou si porovnávané signály podobnější. Nejvyšší hodnotu tedy logicky mají signály shodné, nicméně neplatí to vždy. Stačí změnit fázi, signály tedy vzájemně posunout, maxima jednoho signálu teď mohou korespondovat s minimy signálu druhého a hodnota skalárního součinu rapidně klesne. Signály, u nichž je hodnota skalárního součinu rovna nula, nazýváme ortogonální. A právě slovníky, jejichž funkce jsou vzájemně ortogonální, jsou nejvhodnější pro analýzu signálů. Nejoblíbenějším slovníkem je slovník odvozený od jedné sinusové vlnky, kdy pomocí násobků původní frekvence vznikne množství rozmanitých funkcí, které jsou však vzájemně ortogonální. Velkou výhodou sinusového signálu je také to, že prochází jakýmkoli lineárním časově-invariantním systémem (LTI) nezměněn, což je výhodné v mnoha situacích. Pro výpočet vezmeme v potaz sinusové vlnky všech frekvencí až do poloviny vzorkovací frekvence signálu (do Nyquistovy frekvence). Pro každou z těchto frekvencí hledáme nejlépe vyhovujjící vlnku zkoušením různé fáze. Vybereme tu, jejíž fáze nejlépe odpovídá originálnímu signálu a určíme skalární součin. Vynesením hodnot skalárního součinu pro všechny frekvence do grafu a jejich spojením spojitou čárou získáme výkonové spektrum, které informuje o energetických poměrech signálu. Každý bod spektra koresponduje se součinem signálu a sinusovky odpovídající frekvence (a optimální fáze). [32] Výpočet spektra můžeme provádět na krátkých epochách signálu. Takovou metodou je krátkodobá Fourierova transformace (STFT z angl. short-time Fourier transform). Vykreslením rozsahu STFT pak získáme spektrogram. Nevykreslujeme však jednotlivá spektra, nýbrž vývoj frekvencí v časově-frekvenčním rovině. Funkčně podobným zobrazením jako spektrogram je zobrazení časově frekvenční roviny na základě Wigner-Villovy distribuce (WVD). Takové zobrazení oproti STFT netrpí na ztrátu informace překryvem okna (leakage effects), zato se zde objevuje problém nežádoucí křížové interference (cross-terms). Řešením může být kombinace WVD a algoritmu Matching Pursuit.

23

7.1.

Wigner-Ville Distribution

Wigner-Villova distribuce uvedená poprvé v roce 1932 nese název po Eugen Wignerovi a Jean-André Villovi. Vychází z rovnice autokorelace signálu, kterou je možno zapsat jako: ∞

∞

𝑅𝑥𝑥 (𝜏) = ∫ 𝑥(𝑡 + 𝜏)𝑥

∗ (𝑡)𝑑𝑡

−∞

= ∫ 𝑥(𝑡)𝑥 ∗ (𝑡 − 𝜏)𝑑𝑡 −∞

Kde 𝜏 značí časové zpoždění, respektive časový posun. Křížová korelace je tak většinou vykreslována jako funkce časového zpoždění. Běžná konvence je vykreslovat časové zpoždění jako pozitivní, zatímco druhá funkce je posouvána doleva vzhledem k funkci první. Drobnou úpravou odpovídající posunu jedné funkce o polovinu časového zpoždění doleva a druhé funkce o polovinu časového zpoždění doprava získáme tvar: ∞ 𝜏 𝜏 𝑅𝑥𝑥 (𝜏) = ∫ 𝑥(𝑡 + )𝑥 ∗ (𝑡 − )𝑑𝑡 2 2 −∞

jehož jádro je základem Wigner-Villovy transformace. Podle Wiener-Kninchinova teorému můžeme druhou mocninu Fourierovy transformace signálu (což odpovídá výkonovému spektru signálu) spočítat pomocí Fourierovy transformace z autokorelace signálu, tedy: ∞

|ℱ(𝑥)|2

𝜏 𝜏 = ℱ (∫ 𝑥(𝑡 + )𝑥 ∗ (𝑡 − )𝑑𝑡) 2 2 −∞

(1)

Kde ℱ značí Fourierovu transformaci signálu počítanou dle vzorce ∞

ℱ(𝑥) = ∫ 𝑒 −𝑖𝜔𝑡 𝑥(𝑡)𝑑𝑡

(2)

−∞

Pokud aplikujeme rovnici (2) na (1), získáme: ∞ ∞ 𝜏 𝜏 |ℱ(𝑥)|2 = ∫ 𝑒 −𝑖𝜔𝜏 (∫ 𝑥(𝑡 + )𝑥 ∗ (𝑡 − )𝑑𝑡) 𝑑𝜏 2 2 −∞ −∞

(3)

Odstraněním prostředního integrálu odpovídajícímu integrování v čase získáme časově závislou spektrální hustotu odpovídající dvourozměrné funkci. Pro analytický (komplexní) signál získáme vztah: ∞ 𝜏 𝜏 𝑊(𝑡, 𝑓) = ∫ 𝑒 −𝑖2𝜋𝑓𝜏 𝑥 (𝑡 + ) 𝑥 ∗ (𝑡 − ) 𝑑𝜏 2 2 −∞

(4)

Což odpovídá Wigner-Villově distribuci signálu. Tato transformace je bezpochyby považována za jednu z fundamentálních transformací pro časově-frekvenční reprezentaci signálu. [32]

24

Pro reálný signál pak rovnici (4) můžeme upravit na: ∞

𝜏 𝜏 𝑊(𝑡, 𝑓) = ∫ 𝑒 −𝑖2𝜋𝑓𝜏 𝑥 (𝑡 + ) 𝑥 (𝑡 − ) 𝑑𝜏 2 2 −∞

(5)

𝜏

Následně můžeme zavést substituci 𝜃 = 2, čímž získáme tvar ∞

𝑊(𝑡, 𝑓) = 2 ∫ 𝑒 −𝑖4𝜋𝑓𝜃 𝑥(𝑡 + 𝜃)𝑥(𝑡 − 𝜃) 𝑑𝜃

(6)

−∞

Což lze také interpretovat jako 𝑊(𝑡, 𝑓) = {2 ∙ ℱ𝜃→2𝑓 (𝑥(𝑡 + 𝜃)𝑥(𝑡 − 𝜃))}

(7)

Kde ℱ𝜃→2𝑓 vyjadřuje zdvojnásobení intervalu vzorkované frekvence, neboť časové zpoždění (nyní značeno 𝜃) se posouvá u každé funkce o celou hodnotu 𝜏. Předpokládejme, že 𝑥(𝑡) je ohraničený signál vzorkovaný Nyquistovým kmitočtem, délky N a počtem vzorků n. Pak můžeme zapsat: 𝑊(𝑡, 𝑓) = 𝑊(𝑛 ∙ ∆𝑡, 𝑘 ∙ ∆𝑓) = 2 ∙ 𝐷ℱ𝑇𝑚→2𝑘 [𝑥(𝑛 + 𝑚)𝑥(𝑛 − 𝑚)]

(8)

kde ∆t odpovídá intervalu vzorkování času, ∆f intervalu vzorkované frekvence a m časovému zpoždění vyjádřeném nyní jako proměnný vektor hodnot zahrnující všechny pozitivní a negativní hodnoty τ včetně nuly pro všechny varianty překryvu signálů při autokorelaci.

Mezi kvadratickými časově-frekvenčními reprezentacemi zaujímá Wigner-Villova distribuce přední místo díky své matematické eleganci a některým základním vlastnostem. I zde se však objevuje běžný problém, Wignerova-Villova distribuce není vždy nulová v časech, kdy je hodnota signálu nula, neboli není vždy nulová pro frekvence nezastoupené ve spektru. Tento jev je označován jako křížová interference (v angl. literatuře cross-terms). Nežádoucí interferenční složky jsou přítomny v různých formách ve všech kvadratických časověfrekvenčních metodách odhadu hustoty energie signálu.

25

Interferenci můžeme názorně vidět na Obr 17, kde je zobrazeno časově-frekvenční

12

12

10

10

8

8

frekvence [Hz]

frekvence [Hz]

spektrum signálu složeného ze dvou úseků harmonických funkcí o frekvencích 2 Hz a 10 Hz. Pro názornost interference jsou tyto harmonické funkce dostatečně odděleny jak časově, tak frekvenčně. Hodnota signálu (frekvence) je uprostřed mezi těmito dvěma úseky nulová, avšak WVD přesto vykazuje v časově-frekvenčním zobrazení energii odpovídající právě interferenci.

6

6

4

4

2

2

0 1.5

1

0.5

0

0

0

2

4

6

8

vykon.spektrum

10

12

14

16

18

16

18

cas [s] vytvoreny signal (f1 = 2 Hz, f2 = 10 Hz)

amplituda

1

0

-1

0

2

4

6

8

10

12

14

20

cas [s]

Obr 17 Časově-frekvenční mapa algoritmu Wigner-Ville Distribution

Na Obr 17 je taktéž patrné, že na okrajích výskytu harmonické funkce je energie málo konvergovaná ve frekvenční rovině. To je způsobeno diskontinuitou začátku a konce harmonické vlny, neboť její ostrý zlom je tvořen složením velkého množství harmonických funkcí. Slabá konvergence hran je nazývána Gibbsův efekt. [32] Interferenci můžeme sice vyhladit, ale bez jakékoli apriorní znalosti signálu není možné vyhladit pouze křížovou interferenci, vyhlazení by ovlivnilo i vlastní interferenci (v angl. literatuře nazývanou auto-terms), která je v obraze žádoucí, neboť zobrazuje skutečné struktury signálu. Využití metody Matching Pursuit pro zobrazení Wigner-Villovy distribuce umožňuje křížovou interferenci úspěšně potlačit.

26

7.2.

Matching Pursuit

Matching Pursuit je algoritmus adaptivní aproximace signálu využívající tzv. redundantního slovníku funkcí. Algoritmus byl poprvé publikován Mallatem a Zhangem v roce 1993 [27]. Signál je rozložen na lineární expanzi atomů vlnek, které jsou vybírány z daného slovníku funkcí tak, aby nejlépe aproximovaly vstupní signál. Běžně je používán Gaborův slovník funkcí (tzv. Gaborovy atomy), který nabízí několik výhod v časově-frekvenční analýze signálu. Tvarem odpovídají Gaborovy atomy Gaussovskému okénku modulovanému kosinovou funkcí. [32], [33] Hodně metod analýzy signálů vychází z myšlenky, že daný neznámý signál 𝑥 lze vyjádřit pomocí součtu známých váhovaných funkcí 𝑔𝑛 které byly vybrány ze slovníku D. Pokud je slovník D složen z ortonormálních bázových vlnek, pak jsou koeficienty 𝑎𝑛 dány pouze skalárním součinem slovníkových funkcí se signálem. 𝑁

𝑥 = ∑ 𝑎𝑛 𝑔𝑛

(9)

𝑛=0

Adaptivní aproximace Matching Pursuit nezíská hned přesné vyjádření původního signálu, ale snaží se postupně nalézt v redundantním slovníku D funkce 𝑔𝛾𝑛 , které by optimálně reprezentovaly daný signál x. Vyjádření původního signálu se tedy stává aproximací za použití jenom některých funkcí 𝑔𝑛 vybraných z redundantního slovníku D. V praxi slovník obsahuje mnohem více potencionálních funkcí 𝑔𝑛 , než je pak skutečný počet M funkcí vybraných pro aproximaci signálu. [32] 𝑀−1

𝑥 ≈ ∑ 𝑎𝑛 𝑔𝛾𝑛

(10)

𝑛=0

Pokud by se měly zkoušet všechny možné kombinace M funkcí ze slovníku, aby byla chyba reprezentace co nejmenší, bylo by to prakticky nereálné kvůli obrovskému počtu kombinací již při malém slovníku D. Matching Pursuit algoritmus proto hledá suboptimální řešení prostřednictvím postupných iterací. [32]

Vlastní MP algoritmus Na začátku algoritmu je residuem signálu samotný vstupní signál x, tedy 𝑅 0 𝑥 = 𝑥. V každé iteraci vybereme ze slovníku funkci (vlnku), která nejlépe aproximuje reziduum signálu,

tj.

taková

vlnka,

která

maximalizuje

𝑔𝛾𝑛 = arg max𝑔𝛾 ∈𝐷 |〈𝑅 𝑛 𝑥, 𝑔𝛾𝑖 〉|. 𝑖

27

skalární

součin

|〈𝑥, 𝑔𝛾 〉|,

tedy

Reziduum je pak postupně zmenšováno vždy po odečtení výsledku předchozí iterace. Postupným skládáním vybraných vlnek dostáváme signál, který s narůstajícím M (počet použitých vlnek) konverguje k originálnímu signálu x [32]:

𝑀−1

𝑥 ≈ ∑ 〈𝑅 𝑛 𝑥, 𝑔𝛾𝑛 〉𝑔𝛾𝑛

(11)

𝑛=0

Slovník Gaborových atomů Pro časově-frekvenční analýzu využíváme slovník Gaborových funkcí, což jsou Gaussovské obálky modulované sinusovou oscilací. Gaussovskou funkci můžeme vyjádřit jako 𝑔𝛾 (𝑡) = 𝐾(𝛾)𝑒 −𝜋(

𝑡−𝑢 2 ) 𝑠

cos(𝜔(𝑡 − 𝑢) + 𝜑),

(12)

kde parametr s značí měřítko, u posun, 𝜔 frekvenci a 𝜑 fázi. Touto rovnicí jsme schopni parametricky popsat širokou škálu tvarů vlnek/funkcí. Gaborův slovník používaný pro algoritmus Matching Pursuit obvykle obsahuje také Diracovy a Fourierovy vlnky. Prakticky se však vždy musíme omezit na diskrétní a konečnou množinu trojrozměrného prostoru parametrů (u, s, 𝜔). [32]

Obr 18 Gaborovy atomy [32]

28

Realizace algoritmu Matching Pursuit Algoritmus MP byl sestaven dle rovnice (11). Hledáme maximální skalární součin signálu a vlnky ze slovníku. Reziduum (v prvním kroce původní signál) o tuto vlnku umenšujeme každou iterací, dokud není splněno ukončovací kritérium. Sečtením použitých vlnek ze slovníku získáme aproximaci původního signálu (viz Obr 19).

EG signal a jeho aproximace 6 EG

5

aproximace MP

napeti [mV]

4 3 2 1 0 -1 -2

0

0.05

0.1

0.15 cas [s]

0.2

0.25

0.3

0.25

0.3

vlnky slovniku pouzite k aproximaci 0.6 0.4

amplituda

0.2 0 -0.2 -0.4 -0.6 -0.8

0

0.05

0.1

0.15 cas [s]

0.2

Obr 19 Matching Pursuit aproximace signálu a použité Gaborovy atomy

29

Důležité jsou použité Gaborovy atomy neboli použité funkce ze slovníku. V případě testování algoritmu (viz Obr 19) bylo ukončovací kritérium nastaveno na 30 iterací. Získali jsme tak celkem 30 atomů (viz Obr 20), jejichž sumace odpovídá aproximaci původního signálu.

Obr 20 Gaborovy atomy použité pro aproximaci MP signálu na Obr 19

Důležitější než dokonalá aproximace původního signálu je však pro nás časověfrekvenční rozložení hustoty energie jednotlivých atomů. Takovouto reprezentaci získáme po aplikaci Wigner-Villovy transformace na tyto vybrané atomy. Časově-frekvenční mapa ukazuje rozložení energie atomů. Mapa je podkladem pro samotnou analýzu elektrogramů.

30

7.3.

Odhad hustoty energie signálu pomocí Matching Pursuit a Wigner-Villovy distribuce Aplikací Wigner-Villovy distribuce (5) přímo na rovnici (10) dostaneme 𝑀−1

𝑊(𝑥) ≈ 𝑊 ( ∑ 𝑎𝑛 𝑔𝛾𝑛 ) =

(13)

𝑛=0

Což odpovídá: 𝑀−1

𝑀−1

𝑀−1

2

∑ 𝑎𝑛 𝑊(𝑔𝛾𝑛 ) + ∑ 𝑛=0

∑ 𝑎𝑛 𝑎𝑘 𝑊(𝑔𝛾𝑛 𝑔𝛾𝑘 )

(14)

𝑛=0 𝑘=1,𝑘≠𝑛

kde 𝑎𝑛 = 〈𝑅 𝑛 𝑥, 𝑔𝛾𝑛 〉. Dvojitá suma zahrnuje křížovou interferenci (cross-terms). Použitím lineární expanze (11), můžeme tuto dvojitou sumu vynechat a vyjádřit časově-frekvenční reprezentaci hustoty energie signálu (zahrnutím pouze auto-terms): 𝑀−1 2

𝜀𝑥 = ∑ |〈𝑅 𝑛 𝑥, 𝑔𝛾𝑛 〉| 𝑊𝑔𝛾𝑛

(15)

𝑛=0

10

10

9

9

8

8

7

7

frekvence [Hz]

frekvence [Hz]

Na Obr 21 můžeme vidět výsledek eliminaci interference. Energie časově-frekvenční mapy přesně odpovídá jak časové, tak frekvenční lokalizaci.

6 5 4

6 5 4

3

3

2

2

1

1

0 0.025

0.02

0.015

0.01

0.005

0

0

0

2

4

6

8

vykon.spektrum

10

12

14

16

18

16

18

cas [s] vytvoreny signal (f1 = 3 Hz, f2 = 6,5 Hz)

amplituda

1

0

-1

0

2

4

6

8

10

12

14

20

cas [s]

Obr 21 Výsledek metody Matching Pursuit společně s WVD

Vzorec (15) je základem programového řešení této práce. Postupně je aplikován na databázi signálů a získané časově-frekvenční mapy reprezentující hustotu energie signálů slouží k následné analýze elektrogramů. 31

8.

Realizace algoritmu

V rámci této práce byla realizována aplikace MPWVD* umožňující časově frekvenční analýzu EG signálů. Realizace řešení byla vytvořena v programovém prostředí MATLAB verze R2012b. Aplikace byla testována na souboru EG dat ze sedmi experimentů, které byly zaměřeny na sledování vlivu ischemie a následné reperfuze na elektrickou aktivitu zvířecího izolovaného srdce (morče). Data byla převedena z formátu *.bin do formátu *.mat a vybrán nejvhodnější ze svodů EG záznamu. Celý záznam měření jednoho experimentu byl vždy v délce cca 7mil vzorků odpovídající přibližně hodině měření. (stabilizace 25 min, ischemie 15 min, reperfuze 15 min). Z důvodu značných artefaktů na počátku každého signálu a taktéž z důvodu sjednocení délky signálu ze všech měření pro další práci s daty, byly signály zkráceny o prvních pět minut signálu (600 000 vzorků) a o redundantní konec měření tak, aby takto získaný signál měl délku 6 000 000 vzorků, tedy celkem 50 min při vzorkovací frekvenci 2000 Hz. Takto připravené EG signály byly použity pro další zpracování a analýzu. Zobrazení jednoho z EG signálu můžeme vidět na Obr 22. Červené linie odpovídají reálné době počátku vyvolané ischemie a počátku reperfuze. Stabilizace

Ischemie

Reperfuze

6

5

4

napeti [mV]

3

2

1

0

-1

-2

-3

0

5

10

15

20

25 cas [min]

30

35

40

45

Obr 22 Celý záznam jednoho EG (morče č.2, svod I)

________________________ * MPWVD – vlastní zkratka pro algoritmus kombinující metody Matching Pursuit a Wigner-Villovy Distribuce

32

50

Pro hrubou orientaci změn v rámci celého naměřeného signálu slouží vykreslení průměrného EG v každé minutě měření (viz Obr 23). Toto zobrazení nelze použít k žádnému dalšímu měření, neboť nezohledňuje nepatrné změny v signálu a náhlé změny není možné vidět okamžitě. Avšak pro zevrubný přehled o vývoji signálu je takovéto zobrazení naprosto dostačující. Červeným puntíkem je označena minuta, ve které byla při experimentu způsobena na izolovaném srdci ischemie. Můžeme si všimnout, že na časovém průběhu se projeví ischemie zřetelnou elevací ST segmentu až o několik minut později. Fialovým puntíkem je označena minuta, kdy byla ischemie ukončena a zahájen proces reperfuze. Zde si můžeme povšimnout poměrně rychlého návratu do takřka stejného stavu jako v průběhu stabilizace. Vliv ischemie a podobnost reperfuze vzhledem k stabilizaci je zobrazen na Obr 24. 1. min

2. min 6 4 2 0 -2

0

20 11. min

6 4 2 0 -2

0

20 12. min

6 4 2 0 -2 0

20

0

20 31. min

6 4 2 0

0

20 41. min

6 4 2 0 -2 20

20

20 33. min

0

20

0

20 43. min

20

0

20

20 16. min

0

0

20

20 35. min

20 17. min

0

20 36. min

20 18. min

0

20 37. min

0

20 38. min

6

0

4

4

4

2

2

2

2

0

0

0

20 44. min

0

20 45. min

6 4 2 0 -2 0

20

0

20 46. min

6 4 2 0 -2 0

20

20 47. min

6 4 2 0 -2 0

20

0

20

20

5

ischemie

4

reperfuze stabilizace

3 2 1 0 -1 -2 0

0.05

0.1

0.15 cas [s]

0.2

0.25

Obr 24 Srovnání průběhů stabilizace, ischemie a reperfuze

33

0.3

20 30. min

0

20 40. min

0

20 50. min

6 4 2 0 -2 0

20 49. min 6 4 2 0 -2

0

Obr 23 Průměrný EG v každé minutě měření

0

20 39. min

6 4 2 0 -2 0

20 20. min

6 4 2 0 -2 0

20 48. min

6 4 2 0 -2 0

20

6 4 2 0 -2

0 0

0

29. min 6 4 2 0 -2

0

4

20 19. min 6 4 2 0 -2

28. min 6 4 2 0 -2

0

0

20

10. min 6 4 2 0 -2

6 4 2 0 -2

27. min 6 4 2 0 -2

0

0

20

9. min 6 4 2 0 -2

6 4 2 0 -2

26. min 6 4 2 0 -2

8. min 6 4 2 0 -2

6 4 2 0 -2

25. min

0

7. min 6 4 2 0 -2

6 4 2 0 -2 0

20 34. min

6 4 2 0 -2 0

20 15. min

6 4 2 0 -2

2 0

6 4 2 0 -2 0

20

6 4

0

0

24. min 6 4 2 0 -2

0

20 42. min

20 14. min

6. min 6 4 2 0 -2

6 4 2 0 -2 0

6 4 2

0

0

23. min

20 32. min

6 4 2 0 -2 0

0 6 4 2 0 -2

6 4 2 0 -2 0

20 13. min

5. min 6 4 2 0 -2

6 4 2 0 -2

22. min 6 4 2 0 -2

0

0

20

4. min 6 4 2 0 -2

6 4 2 0 -2

21. min 6 4 2 0 -2

3. min 6 4 2 0 -2

napeti [mV]

6 4 2 0 -2

20

0

20

8.1.

Detekce QRS

Detekce QRS komplexu hraje důležitou roli při každém rozměření EG signálu. Komplex QRS je nejvýraznějším útvarem celé EG křivky, a díky tomu je i vhodným orientačním bodem. Z počtu QRS v minutě signálu lze určit tepová frekvence (HR z angl. Heart rate), srdeční variabilita, ale především je komplex QRS důležitý pro rozměření typických bodů a segmentů EG křivky. Přístupů detekce QRS je celá řada. Nejrozšířenějšími jsou algoritmy využívající první derivace, algoritmy založené na vlnkové transformaci, genetické algoritmy či neuronové sítě, dále například skryté Markovovské řetězce či přizpůsobená filtrace. Více například v [15] a [34]. Tato práce není primárně zaměřena na navržení nového přístupu detekce QRS, ale spíše na vlastní analýzu EG cyklů (úseků PQRST), a to v časově-frekvenčním měřítku. Možnostem detekce QRS a jejím pokročilým metodám tedy nebude věnováno tolik prostoru. Avšak navržená metoda časově-frekvenční analýzy pracuje s konkrétními cykly EG a ne s celým signálem naráz, proto je zde detekce QRS využita pro lokalizaci a indexaci jednotlivých EG cyklů. Pro specifičnost průběhu analyzovaných dat byl použit vlastní přístup detekce QRS. Princip detekce QRS odpovídá obecnému schématu detekce (viz Obr 9). Díky častému impulsnímu rušení projevujícím se jako úzké píky velké výchylky, se jevilo jako vhodné použít mediánovou filtraci. Ta by takovéto nežádoucí impulsy snižující účinnost detektoru spolehlivě potlačila, avšak nakonec mediánová filtrace použita nebyla, neboť její aplikací se výrazně snížila (a zkreslila) výchylka také R vln, které pak detektor nedokázal rozeznat od vln T. Dalším nemalým problémem byly ořezané špičky QRS v některých svodech způsobené nastavením velkého zesílení předzesilovače během záznamu měření. Proto byla navržena vlastní metoda eliminující falešně pozitivní detekci spočívající ve stanovení rozhodovacího pravidla porovnávající velikosti výchylek vln. Velikost maximální výchylky vlny detekované detektorem jako vlna R musí být v rozmezí 60-140% velikosti předchozí R vlny. Zavedením spodní hranice se podařilo potlačit falešné označení vln T, které v některých naměřených signálech dosahují přibližně polovinu výchylky vlny R. A díky omezení horní hranicí se podařilo většinově vyloučit detekci impulsního rušení. Další předzpracování detekce spočívá v zvýraznění QRS umocněním signálu. Stanovení prahu detekce je automaticky nastaveno pouze na 10% velikosti nejvyšší výchylky signálu. Takto nízká hodnota prahu je volena z toho důvodu, že ischemické QRS některých měření mají mnohem menší výchylky než je tomu ve fázi stabilizace, a tak by vysoký práh způsobil při ischemii nezaznamenání žádného QRS. Byla testována i adaptivní velikost prahu, jehož hodnota se vždy aktualizovala vzhledem k předešlým detekcím, nicméně tento přístup měl problém se začákem reperfuze, kdy dochází ke skokovému navýšení výchylky QRS. Samotná detekce pozice R vlny je realizována hledáním nejprve nadprahové hodnoty a následně určením maxima v plovoucím okně. Tento princip eliminuje chybu detektoru označit několik R vln těsně vedle sebe.

34

Úspěšnost detekce V následující tabulce jsou uvedeny výsledky testování detektoru přímo na experimentálních datech. Míru kvality detektoru můžeme vyjádřit parametry senzitivitou (Se) a pozitivní prediktivitou (P+). Sensitivita detektoru vyjadřuje jeho schopnost správně detekovat QRS komplex. Pozitivní prediktivita vyjadřuje vztah mezi falešnou a reálnou detekcí, zda označená vlna patří skutečně QRS komplexu či nikoli. Vypočítat je můžeme pomocí vztahů:

𝑆𝑒 =

𝑇𝑃 𝑇𝑃 + 𝐹𝑁

(16)

𝑃+ =

𝑇𝑃 𝑇𝑃 + 𝐹𝑃

(17)

Kde TP (true positive) značí počet správně detekovaných událostí, FN (false negative) jsou události, které nejsou detekovány, ale měli být, a FP (false positive) jsou detekované události, které však nenáleží hledanému objektu.

Tabulka 2 Přehledová tabulka účinnosti detektoru QRS

Počet QRS

FN

FP

Se [%]

P+ [%]

morče č. 1

4563

0

0

100,00

100,00

morče č. 2

4810

0

1

100,00

99,98

morče č. 3

5058

0

0

100,00

100,00

morče č. 4

5359

2

0

99,96

100,00

morče č. 5

4683

32

0

99,32

100,00

morče č. 6

5778

2

8

99,97

99,86

morče č. 7

5014

1

3

99,98

99,94

99,90

99,97

senzitivita a pozitivní prediktivita detektoru:

V Tabulka 2 jsou zobrazeny hodnoty senzitivity a pozitivní prediktivní hodnoty pro jednotlivá měření. Detektor byl testován na všech signálech pro fázi stabilizace, kdy byl vytvářen referenční model. Dole je taktéž uvedena hodnota napříč měřeními, čímž byla získána hodnota sensitivity a pozitivní prediktivní hodnoty detektoru.

35

Na Obr 25 můžeme vidět část signálu EG. Červenými trojúhelníky jsou označeny detekované R vlny. Na základě detekce R vln mohly být z celého signálu indexovány jednolité EG cykly (PQRST úseky), na které je následně postupně aplikován algoritmus Matching Pursuit. Zelené trojúhelníky značí počátek a fialové konec jednoho EG cyklu. Začátek a konec indexace jednotlivých úseků je dán dle velikosti RR intervalu. Celkem délka indexovaného EG cyklu odpovídá 120% délky RR intervalu (přesah je z důvodu měnících se intervalů QT apod.) s rozdělením na 40% RR před detekovanou R vlnou (počátek indexace) a 80% RR za ní (konec indexace). Tyto hodnoty byly experimentálně otestovány jako dostačující, aby zahrnuly právě jeden EG cyklus, a zároveň do něj nebyla zahrnuta například vlna T z předchozího, či vlna P z následujícího cyklu. Případně je také možnost změnit tyto hodnoty manuálním nastavením detektoru.

originalni signal detekovane R vlny

5

4

napeti [mV]

3

2

rozmezi indexace jednoho EG cyklu

1

0

-1

-2 1.664

1.6645

1.665

1.6655

1.666

1.6665 vzorky [-]

Obr 25 Detekce R vln

36

1.667

1.6675

1.668

1.6685

1.669 6

x 10

Na Obr 26 je zobrazeno prvních pět sekund z každé minuty na počátku ischemie. V první minutě se ještě vliv ischemie na tvaru EG křivky nikterak neprojevil. Obdobně je tomu v průběhu druhé a třetí minuty.

napeti [mV]

1.minuta ischemie 4 2 0 -2

napeti [mV]

napeti [mV]

napeti [mV]

napeti [mV]

0

0.5

1

1.5

2

2.5 3 cas [s] 2.minuta ischemie

3.5

4

4.5

5

4 2 0 -2 60

60.5

61

61.5

62

62.5 63 cas [s] 3.minuta ischemie

63.5

64

64.5

65

4 2 0 -2 120

120.5

121

121.5

122

123.5

124

124.5

125

4 2 0 -2 180

180.5

181

181.5

182

183.5

184

184.5

185

4 2 0 -2 240

240.5

241

241.5

242

243.5

244

244.5

245

122.5 123 cas [s] 4.minuta ischemie

182.5 183 cas [s] 5.minuta ischemie

242.5 cas [s]

243

Obr 26 Počátečních 5 min ischemie

Výraznější změny jsou patrné až v páté minutě, kde se projeví elevace ST segmentu. Změny v důsledku ischemie se tedy projevují pozvolna a odhadem z tvaru EG křivky lze ischemii určit až relativně pozdě. Proto neustále probíhá vývoj nových přístupů detekce ischemie za použití nejrůznějších příznaků (například [35]). Je obecnou snahou správně automaticky detekovat ischemii co nejdříve. Nejčastěji se používají různé principy na základě podrobného rozměření signálu a zjištění souvislostí změn. Ke změnám na EG signálu však dochází nejen v časové, nýbrž i frekvenční rovině, proto má velký význam právě časověfrekvenční analýza.

37

8.2.

Vznik časově-frekvenční mapy

Princip, kterým je tvořen časově-frekvenční mapa, představuje následující série obrázků. Nalevo je vždy zobrazováno časově-frekvenční spektrum (WVD) jednotlivých atomů, které byly použity k aproximaci MP původního signálu. Algoritmus Matching Pursuit vybírá ze slovníku atomy s postupně klesající úrovní energie. To znamená, že složky signálu s vysokou energií jsou v podobě atomů extrahovány jako první. Napravo můžeme vidět postupnou sumaci této energie, z níž získáme výstupný hledaný obraz. Obrázek je doplněn o výkonová spektra ukazující aktuální výkon signálu, respektive v případě rozkladu na atomy nejvýraznější frekvence signálu. Následuje série obrázků ilustrující vývoj získávání časově-frekvenční mapy. vlnka c.1, fmax = 44 Hz 1000 900

900

900

800

800

800

800

700

700

700

700

600

600

600

600

500

500

500

500

400

400

400

400

300

300

300

300

200

200

200

200

100

100

100

100

0 0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0 5000 0 vykon.spektrum

0

4

4

2

2

napeti [mV]

0 5000 0 vykon.spektrum

frekvence [Hz]

900

napeti [mV]

frekvence [Hz]

1000

0 -2 -4

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0 -2 -4

Obr 27 Průběh vzniku časově-frekvenční mapy

Mezi první atomy (viz Obr 27) rozhodně patří takové, které zastupují QRS komplex, neboť ten je nejvýraznější částí EG signálu. Postupně se objevují také atomy reprezentující vlnu T (viz Obr 28).

38

vlnka c.8, fmax = 44 Hz 1000

900

900

900

900

800

800

800

800

700

700

700

700

600

600

600

600

500

500

500

500

400

400

400

400

300

300

300

300

200

200

200

200

100

100

100

100

0 100 50 0 vykon.spektrum

frekvence [Hz]

frekvence [Hz]

1000

0 0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0

0.35

0 2 1 0 vykon.spektrum 4 x 10

0.4

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

4 napeti [mV]

napeti [mV]

0.05

6

0.2 0 -0.2

2 0

-0.4 -0.6

0

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

-2

0.35

Obr 28 Průběh vzniku časově-frekvenční mapy

Nejvyšší frekvence mají atomy aproximující zákmit Q, špičku R vlny a kmit S. Ty dosahují frekvencí až 350 Hz (viz Obr 28), přestože maximální frekvence celého komplexu QRS měřených EG nabývá hodnot 30-40 Hz. To je dáno tím, že QRS komplex je rozložen atomy různých frekvencí, a kromě „nosných“ funkcí se na rozkladu QRS podílí právě i funkce tak vysokých frekvencí. vlnka c.13, fmax = 350 Hz 1000 900

900

900

800

800

800

800

700

700

700

700

600

600

600

600

500

500

500

500

400

400

400

400

300

300

300

300

200

200

200

200

100

100

100

100

0

0

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0 2 1 0 vykon.spektrum 4 x 10

1

6

0.5

4

napeti [mV]

0 10 5 0 vykon.spektrum

frekvence [Hz]

900

napeti [mV]

frekvence [Hz]

1000

0 -0.5 -1

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

2 0 -2

Obr 29 Průběh vzniku časově-frekvenční mapy

39

vlnka c.30, fmax = 0 Hz 1000

900

900

900

900

800

800

800

800

700

700

700

700

600

600

600

600

500

500

500

500

400

400

400

400

300

300

300

300

200

200

200

200

100

100

100

100

0

0

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0 2 1 0 vykon.spektrum 4 x 10

0.05

6

0

4 napeti [mV]

napeti [mV]

0 10 5 0 vykon.spektrum

frekvence [Hz]

frekvence [Hz]

1000

-0.05 -0.1 -0.15

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0.35

2 0 -2

Obr 30 Průběh vzniku časově-frekvenční mapy

Na Obr 30 vpravo můžeme vidět celkovou časově-frekvenční mapu složenou z 30 časově-frekvenčních spekter jednotlivých 30 atomů vybraných z Gaborova slovníku. Pod touto časově-frekvenční mapou pak vidíme aproximaci původního EG vybranými 30 atomy. Počet atomů, které použijeme k aproximaci signálu, a tím k vytvoření časově-frekvenční mapy není teoreticky limitován (generovaný slovník jich obsahuje několik tisíc), avšak s větším počtem atomů roste samozřejmě i výpočetní náročnost. Výhodou časově frekvenční analýzy pomocí MPWVD je, že ve velké míře postačí malý počet atomů. Čím méně atomů je použito, tím spíše jejich složení nevytvoří naprosto přesnou aproximaci původního signálu, nicméně i malé množství atomů postačí k zobrazení právě těch nejvýznamnějších frekvencí s lokalizací v čase.

8.3.

Zobrazení detailů a energie signálu v časově-frekvenční mapě

Po výběru úseku signálu byla tedy provedena detekce R vln a podle toho byl rozdělen signál na dílčí úseky tak, aby každý úsek obsahoval právě jeden PQRST záznam. Poté je vytvořen slovník Gaborových atomů a zvolen práh maximálního počtu vybraných atomů ze slovníku. Tento práh, který odpovídá počtu iterací rozkladu signálu metodou Matching Pursuit, může být libovolně zvolen. Pro detekci největšího množství energie signálu odpovídající QRS komplexu stačí práh relativně nízký (hodnoty 5-15), zatímco pro zobrazení jednotlivých drobných změn a detailů v signálu je lepší volit práh větší (25-50). Celkové zobrazení časověfrekvenční mapy je realizováno kombinací metody Matching Pursuit a Wigner-Villovy distribuce signálu pro každý atom vybraný ze slovníku. Tím získáme rozložení energie

40

jednotlivých atomů a jejich sumačním zobrazením pak energii signálu v časově-frekvenčním měřítku. Z předchozích obrázků je patrné, že shluky energie odpovídající frekvenčním změnám jsou zároveň přesně lokalizované v čase. Získali jsme tedy dostatečně přesné frekvenční rozlišení při zachování perfektního časového rozlišení. Při sumaci běžných časově-frekvenčních map bylo použito vždy přeškálovaných hodnot spektra každého atomu tak, aby maximální hodnota odpovídala jedné. To umožňuje zobrazení velkých detailů, neboť při sčítání nejsou potlačeny méně výrazné atomy. Takovéto zobrazení je užitečné pro hledání detailů v signálu, jako například jednoduchého a zároveň přesného určení vln P, Q, a T. Pro porovnávání časově-frekvenčních map různých signálů, jako například klidového a ischemického, se však jeví tento princip jako méně vhodný. Neboť při velmi malé změně ve tvaru EG je použit jiný soubor atomů (vybíráno dle maximálního skalárního součinu, viz dříve) a vzhledem k přeškálování hodnot atomů tak vzniká výrazně odlišný obraz. Není tedy prakticky možné stanovit reprezentační soubor atomů, respektive časově-frekvenční mapu, představující například standartní hodnoty při stabilizaci, či ischemii. Pokud ovšem budeme časově-frekvenční mapy jednotlivých atomů přímo sčítat beze změny škály hodnot, získáme časově-frekvenční mapu zobrazující největší energii EG v časově-frekvenčním měřítku. Porovnáme-li tedy první přístup zobrazení MPWVD se škálováním (MPWVD-S) s druhým přístupem MPWVD s důrazem na energii signálu (MPWVD-E), vidíme sice zřetelný rozdíl, kdy energické MPWVD nenabídne tolik detailů (viz Obr 31), avšak právě přístup MPWVD-E se osvědčil při analýze ischemických změn. MPWVD energy 500

500

450

450

450

450

400

400

400

400

350

350

350

350

300

300

300

300

250

250

250

250

200

200

200

200

150

150

150

150

100

100

100

100

50

50

50

50

0

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0 5000 0 vykon.spektrum

0.3

0

6

6

4

4

napeti [mV]

0 5000 0 vykon.spektrum

frekvence [Hz]

500

napeti [mV]

frekvence [Hz]

MPWVD scale 500

2 0 -2

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

0

0.05

0.1

0.15 0.2 cas [s]

0.25

0.3

2 0 -2

0.3

0

Obr 31 Zobrazení energie signálu v časově-frekvenčním měřítku

41

MPWVD-S Na Obr 32 nahoře je zobrazena časově-frekvenční mapa MPWVD-S, která poskytuje dostatek detailů pro rozměřený významných bodů. Na základě testování souboru EG signálů byly vytvořeny kritéria pro automatickou detekci začátku vlny P, Q a konce vlny T. Právě tyto jinak obtížně stanovitelné zlomové body vln lze díky rozkladu signálu na atomy poměrně snadno nalézt. Díky tomu je možné též stanovit délku PQ a QT intervalu (viz Obr 32 dole). MPWVD-S MPWVD-E

500

400

300

200

100

0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

cas [s]

PQ interval = 44.5 [ms], QT interval = 170 [ms] 6 MP aproximace 5

originalni signal

4

zacatek P zacatek Q

3

konec T

2 1 0 -1 -2

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Obr 32 Rozměření bodů P, Q a T na základě MPWVD-S

Toto rozměření bylo testováno ve fázi měření stabilizace. Vzhledem k charakteristice experimentu, který se soustředí na změny po navození ischemie a poté návratu k standartním hodnotám během reperfuze, je tato práce koncipována s větším zaměřením na časověfrekvenční analýzu ischemických dat. Následující část práce se tak bude věnovat možnostem automatické detekce ischemie.

42

MPWVD-E Zobrazíme-li časově-frekvenční mapy MPWVD-E v různých fázích měření, jsou viditelné zřetelné rozdíly jak mezi časovým průběhem, tak v časově-frekvenčních mapách (viz Obr 33). Díky změně EG křivky v návaznosti na elevaci ST segmentu provázejícího ischemii byly použity pro dekompozici EG signálu jiné atomy, což se projevilo i odlišným výskytem hlavního shluku energie EG signálu v časově-frekvenční mapě. V mapě prezentující ischemii lze vidět velké zastoupení nízkých frekvencí, které odráží použití většiny atomů o nízkých frekvencích. Vlivem elevace ST dochází k úplnému potlačení kmitu S a zmírnění ostrosti zákmitu v Q a R. Při použití nízkého počtu atomů (zde 10) jsou při dekompozici pomocí MP použity pouze atomy reprezentující zevrubný tvar křivky, což umocní potlačení ostrosti v kmitech Q a R. Na základě porovnání časově-frekvenční mapy stabilizace a ischemie lze tak snadno určit, která mapa odpovídá ischemii. Ve fázi reperfuze vidíme takřka totožné časověfrekvenční spektrum odrážející skutečnost, že došlo k návratu do klidových hodnot stabilizačního stavu. ischemie

reperfuze

450

450

450

400

400

400

350

350

350

300

300

300

250 200

frekvence [Hz]

500

frekvence [Hz]

500

250 200

250 200

150

150

150

100

100

100

50

50

50

0

5

10

15 20 cas [ms]

25

0

30

0

5

10

15 20 cas [ms]

25

0

30

6

6

6

4

4

4

2 0 -2

amplituda [mV]

amplituda [mV]

0

amplituda [mV]

frekvence [Hz]

stabilizace 500

2 0 -2

0

4

8

12

16 20 cas [ms]

24

28

32

0

5

10

15 20 cas [ms]

25

30

2 0 -2

0

4

8

12

16 20 cas [ms]

24

28

32

0

4

8

Obr 33 Porovnání časově-frekvenčních map (stabilizace, ischemie, reperfuze)

43

12

16 20 cas [ms]

24

28

32

8.4.

Srovnání s jinými metodami časově-frekvenční analýzy

Jako základní časově-frekvenční zobrazení lze považovat Fourierovu Transformaci. Ta však neumožňuje žádnou lokalizaci v čase. Toto omezení do jisté míry řeší modifikace nazývaná krátkodobá Fourierova transformace (STFT). STFT je však omezena pevnou velikostí okna, což může působit problém v dosažení optimálního poměru rozlišení v čase a frekvenci. Výhodou spojité vlnkové transformace (CWT) oproti STFT je proměnné časověfrekvenční rozlišení, čímž dokážeme nalézt optimální rozlišení pro libovolný kmitočet.

Obr 34 Srovnání metod Matching Pursuit, Spojité vnlnkové transformace a spektrogramu STFT

Na Obr 34 lze vidět srovnání časově-frekvenčních map STFT, CWT a MPWVD. Histogram STFT je na tom z uvedených možností jednoznačně nejhůře. Nenabízí dostatečné časové ani frekvenční rozlišení, respektive jsme nuceni si vybrat mezi volbou buď lepšího časového, nebo frekvenčního rozlišení. CWT je na tom mnohem lépe, jsou zřetelné vlny R a T s poměrně přesnou časovou lokalizací. Nejlépe je na tom zobrazení vycházející z algoritmu Matching Pursuit. Její výhodou je získání hodnot výrazných frekvenčních složek s časovým výskytem podle použitých bázových funkcí (atomů). U CWT a STFT to vyžaduje podrobnější analýzu mapy. Metoda MPWVD-S tedy nabízí rozměření významných bodů obdobně, jak to umožňuje spojitá vlnková transformace, která je v tomto ohledu robustnější metodou V následující části práce budou porovnány možnosti využití časově-frekvenční analýzy MPWVD-E a CWT k detekci ischemie.

44

8.5.

Automatická detekce srdeční ischemie

Nejprve byl sledován vývoj změn v časově-frekvenční mapě. V rámci jednotlivých fází měření (stabilizace, ischemie a reperfuze) byly změny minimální. K největším změnám došlo při přechodu ze stabilizační fáze do ischemické. Bylo spočítáno více než 300 časověfrekvenčních map z průběhu stabilizace a ischemie každého měření, a to jak pomocí MPWVDE, tak CWT (viz Obr 35 a Obr 36). A z každého takového souboru byla vytvořena průměrná časově-frekvenční mapa reprezentující stabilizaci a ischemii. Stejným způsobem byla spočítána průměrná časově-frekvenční mapa i pro reperfuzi. Tím byly získány referenční modely pro dané fáze měření. cas: 4min 58s 500 450 400

frekvence [Hz]

350 300 250 200 150 100 50 0

0

5

10

15 cas [ms]

20

16 cas [ms]

20

25

6

napeti [mV]

4 2 0 -2

0

4

8

12

24

28

Obr 35 Průběh vytváření průměrné časově-frekvenční mapy pomocí MPWVD-E

meritko

cas: 4min 58s

100

200

300 cas [ms]

400

500

600

6

napeti [mV]

4 2 0 -2

0

4

8

12

16 cas [ms]

20

24

28

Obr 36 Průběh vytváření průměrné časově-frekvenční mapy pomocí CWT

45

K analýze získaných časově-frekvenčních map byly použity histogramy. Umožní nám zobrazit „četnost“ energie vytvořené v časově-frekvenční mapě atomy vzhledem k dané frekvenci. Tím získáme přehled, jak velká energie (důležitá složka signálu, neboť tam byl použit atom) přísluší konkrétní frekvenci. Tento histogram tak není pouhým vyjádřením množství jednotlivých frekvencí elektrogramu, což by vyjadřovalo výkonové spektrum (PSD), ale reprezentuje zastoupení energie atomů. Pokud zobrazíme histogramy pro různé fáze měření, vidíme, že během stabilizace jsou nejvýraznější frekvenční složky v rozsahu 40-50 Hz, během ischemie dochází k frekvenčním změnám projevujících se výskytem komponent s majoritní energií na frekvencích 2-20 Hz. Během reperfuze dochází k postupným změnám zpět, až po určité době dojde k ustálení na hodnotách velmi blízkých stabilizační fázi.

Obr 37 Srovnání histogramů pro stabilizaci, ischemii a reperfuzi na celém frekvenčním spektru

Na obrázku Obr 37 můžeme vidět relativní zastoupení atomů při stabilizaci, ischemii a reperfuzi. Nárůst počtu atomů při ischemii je tak obrovský, že lze velmi těžko rozeznat signály stabilizace a reperfuze. Přiblížením získáme následující výřez.

46

Obr 38 Detail změn frekvenčního zastoupení použitých atomů při MPWVD-E

Z Obr 38 je patrné, že tvar histogramu reperfuze těsně kopíruje histogram stabilizace, zatímco histogram ischemie má, jak již bylo řečeno, enormně zvýšené hodnoty četnosti atomů patřící nízkým frekvencím. Tento jev byl pozorován u všech měření, a proto se od něj odvíjí metody detekce ischemie. Získáním průměru všech histogramů měření stabilizace byl stanoven referenční model reprezentující stabilizaci. Budeme-li vykreslovat histogram postupně pro všechny EG cykly, vidíme postupnou změnu od histogramu stabilizačního k ischemickému. K větším změnám dochází již relativně brzy, proto se mohou stát podkladem pro automatickou detekci ischemie. Nejprve je třeba stanovit referenční histogram reprezentující stabilizační fázi. Nejjednodušší cesta, a taktéž nejpřesnější, je ze získaného signálu, u kterého známe hranice stabilizace, vzít celý úsek stabilizace, pro každý EG cyklus spočítat MPWVD časověfrekvenční mapu a histogram a z nich průměrem stanovit referenční model. Výpočet pro úplně všechny EG cykly stabilizace však není nutný, neboť během stabilizace se signál výrazně nemění, a tak můžeme vzít v úvahu buď celou délku stabilizace, ale z každé minuty jen několik EG cyklů, nebo naopak zvolit kratší časový interval, z nějž hustěji vybereme EG cykly a stanovíme referenční model stabilizace. Nejlepší volbou je vytváření referenčního modelu z konečné fáze stabilizace odpovídající nejlepšímu stavu srdce. To je však možné pouze při práci s experimentálními daty, u nichž je známo, kdy k ischemii dochází. Pro práci se signálem bez jakékoli apriorní znalosti je třeba referenční model vytvářet postupně, aktualizovat jej v čase, a vytvořit tak model, který bude nejlépe odpovídat průběhu stabilizace srdce.

47

Obr 39 Histogramy stabilizačních fází měření. Zavedení průměrného referenčního histogramu.

Na Obr 39 je uveden celkový histogram a přiblížení detailu průběhu stabilizací. Červené hodnoty značí vytvořený histogram referenčního modelu. Následně byly procházeny jednotlivé signály a počítána postupně MPWVD u vyskytujících se EG cyklů. Histogram každé časověfrekvenční mapy byl porovnáván s histogramem časově-frekvenční mapy referenčního modelu a dle nastaveného prahu zjišťován výskyt ischemie. Vzhledem k malému souboru dat a velkému rozptylu stabilizačních histogramů, byla podmínka detekce ischemie nastavena na dvojnásobek hodnot referenčního modelu. Po překročení těchto hodnot, detektor zaznamenal ischemii. Účinnost byla 100%, ischemii se podařilo detekovat vždy do 3 minut, někdy i dříve. Výsledky by byly daleko přesnější při vytváření referenčního modelu pouze z daného signálu, či při dostatečně velkém souboru měření. Velký vliv má také správná detekce R vln. Byly tedy porovnávány histogramy všech naměřených signálů a z nich stanoven referenční model. Metoda je robustnější, univerzálnější, neměl by být problém ji aplikovat na jakýkoli EG či EKG signál, nicméně kvůli univerzálnosti není tak přesná, neboť musíme nastavovat vyšší hodnoty prahu detekce, abychom dosáhli stoprocentní úspěšnosti detekce, což je žádoucí. To je vykoupeno delší dobou, od vzniku ischemie po její detekci. Případně můžeme práh nastavit přísně dle referenčního modelu, čímž získáme skutečně univerzální detektor s včasnou detekcí, ovšem s horší účinností. Pak bychom jen procházeli signál, pro aktuální detekovanou R vlnu určili celý EG cyklus jí náležící, spočítali MPWVD, histogram a porovnáním s histogramem referenčního modelu určili ischemii. Detekce zohledňuje právě nejvýraznější změny v oblasti 2-10 Hz.

48

Nakonec byl zvolen přístup kombinující obě předchozí metody. Je navržen pro libovolný signál, ale referenční model je určen „online“ přímo adaptivně během procházení signálu. Není tedy nutná žádná apriorní znalost signálu a přitom je referenční model vytvořen pro signál na míru, což přináší vyšší rychlost a spolehlivost detekce. Na Obr 40 je zjednodušený vývojový diagram detekce ischemie pomocí metody MPWVD-E. START

vstupní signál odhad RR určení okna

procházení signálu oknem

nalezena R vlna

NE

posun v signalu o delku okna

ANO určení bodu počátku a konce EG indexace EG cyklu výpočet MPWVD vytvoření aktuálního histogramu

posun k dalšímu RR intervalu

počet histogramů zahrnutých do referenčního modelu <100

NE

ANO přidání hist. do referenčního modelu přepočítání referenčního modelu

pronásobení referenčního modelu konstantou

NE

zastoupení vybraných frekvencích aktuálního histogramu > zastoupení vybraných frekvencích referenčního modelu

ANO detekována ischemie

KONEC

Obr 40 Zjednodušený vývojový diagram metody automatické detekce ischemie na základě MPWVD-E

49

Na Obr 41 je zobrazen průběh detekce ischemie na základě MPWVD. V titulku obrázku křivky EG se ukazuje aktuální čas detekované R vlny a v závorce je čas, kdy byla reálně vyvolána ischemie. Kvalita referenčního modelu znamená, kolik časově-frekvenčních map bylo použito k získání referenčního modelu. Ukázalo se, že 100 map je dostačující počet pro robustní referenční model. Takže po dosáhnutí alespoň 100 iterací algoritmu, což odpovídá 100 časověfrekvenčním mapám, je kvalita referenčního modelu 100%.

Obr 41 Průběh detekce ischemie pomocí MPWVD

Na Obr 42 je stejné okno o několik minut později. Hodnoty aktuálního histogramu na nejnižších frekvencích přesáhly stanovenou mez a byla detekována ischemie. Pozitivní detekce je oznámena velkým červeným textem ISCHEMIE, který se objeví u dané EG křivky.

Obr 42 Detekovaná ischemie

50

Obdobný princip byl použit i pro testování detekce na základě CWT. Byly spočítány koeficienty časově-frekvenční mapy a vykreslen jejich histogram. Průměrováním histogramů jednotlivých fází měření byly získány tři histogramy reprezentující tři fáze měření (viz Obr 43).

Obr 43 Histogramy CWT ( stabilizace, ischemie a reperfuze)

I zde je histogram odpovídající reperfuzi takřka totožný s histogramem stabilizace. Oproti tomu histogram ischemie vykazuje značné výchylky. Tyto změny byly sledovány a dle jejich výskytu byl stanoven práh detekce ischemie. Průběh sledování změn při detekci ischemie je zobrazen na Obr 44. Cas: 14:59 (isch.vyvolana v: 23:23)

aktualni EG cyklus

50

6 prah detekce 5

-50 4 -100 -150

-200

3

0

50

100

150

200

250

meritko

meritko

CWT 248 235 222 209 196 183 170 157 144 131 118 105 92 79 66 53 40 27 14 1

napeti [mV]

suma cwt na danem meritku

histogram CWT 0

2

1

0

-1

-2

100

200

300 cas [s]

400

500

600

-3

0

Obr 44 Průběh detekce ischemie pomocí CWT

51

100

200

300 cas [s]

400

500

600

Metoda detekce pomocí CWT byla otestována na všech signálech měření (pozitivní detekce na Obr 45). V Tabulka 3 je přehled časů, během nichž byla ischemie detekována vzhledem ke skutečnému začátku ischemie. ISCHEMIE! V case: 25:19, realny cas vyvolane ischemie: 23:23, tj. po: 117 s

aktualni EG cyklus

50

6 prah detekce 5

-50 -100

4

-150 3

-200 -250

0

50

100

150

200

250

napeti [mV]

suma cwt na danem meritku

histogram CWT 0

meritko

meritko

CWT 248 235 222 209 196 183 170 157 144 131 118 105 92 79 66 53 40 27 14 1

2

1

0

-1

-2

100

200

300 cas [s]

400

500

-3

600

0

100

200

300 cas [s]

400

500

600

Obr 45 Pozitivní detekce ischemie pomocí CWT Tabulka 3 Přehled časů detekce ischemie metodou MPWVD-E a CWT

měření

MPWVD-E [s]

CWT [s]

morče č. 1

110

117

morče č. 2

152

231

morče č. 3

131

169

morče č. 4

136

195

morče č. 5

281

291

morče č. 6

115

126

morče č. 7

120

172

Teoreticky je možné použít metodu MPWVD přímo na celý úsek signálu, pomocí ní detekovat QRS, signál rozměřit a analyzovat. Metoda však vytváří matici o rozměru NxN, kde N odpovídá délce vybraného úseku signálu. Proto při velkém úseku signálu tak získáme matici velmi velkých rozměrů. To vytváří velké nároky na paměť a náročnost výpočtů vůbec. Proto při vysokofrekvenčních EG byl použit algoritmus pracující s jednotlivými EG cykly určenými dle jednoduché prahové detekce QRS. Další možností by bylo signál podvzorkovat, čímž však ztrácíme některé důležité detaily a dochází ke zkreslení. Případně by mohla být kompromisem úprava algoritmu na nový detektor, s velkým plujícím oknem přes několik EG cyklů, v rámci něhož by se metoda počítala. Zde by bylo potřeba empiricky zjistit vhodnou délku okna pro zachování výpočetní rychlosti při dostatečné informační výtěžnosti. 52

9.

Diskuse výsledků

V Tabulka 2 je uveden souhrn výsledků testování přesnosti detektoru QRS. Elektrogramy vykazují komplikovaný průběh kvůli výskytu impulsního rušení a nepatřičných zákmitů. Navíc některé ze signálů databáze trpí ořezem špičky QRS, které bylo způsobeno nesprávnou akvizicí. To vše se podílí na problematické detekci R vlny. Obecně se ukázalo u EG signálů z izolovaných srdcí nezbytné nastavit práh na nižší hodnotu, neboť při vyšší přísnější hodnotě nebyl vůbec zaznamenán výskyt R vln po dobu ischemie, kde velikost výchylky R vlny radikálně klesá, než je její velikost ve fázi reperfuze znovu obnovena. Zavedení nižšího prahu prospěje snížení FN, avšak na druhou stranu značně navýší FP. Proto bylo stanovena podmínka detektoru stanovující minimální a maximální velikost výchylky. Díky tomu se falešně pozitivní detekce způsobená především vlnami T s velkou úspěšností eliminovala. Tato podmínka byla vytvořena s adaptivní aktualizací mezních hodnot na základě předchozí korektní detekce. Detekce QRS použitá při analýze se tedy ukázala jako vhodná a velmi přesná, což dokazují i hodnoty senzitivity a pozitivní prediktivní hodnoty. Tato metoda detekce je však navržena přímo pro EG signály s typickým průběhem. V případě aplikace metody detekce QRS na reálné EKG signály by musela být metoda rozšířena o vhodnou filtraci (myopotenciálového šumu, brumu, korekci nulové izolinie a další) a upraveny její parametry. Prezentovaný algoritmus Matching Pursuit se prokázal na základě hodnocení kvality aproximace původního signálu jako naprosto vyhovující. Již použitím několika málo atomů ze slovníku funkcí dokázal s dostatečnou přesností nahradit původní signál se zanedbatelným reziduem. Díky jednoduchosti a efektivnosti výběru vhodné vlnky na základě skalárního součinu je algoritmus velmi rychlý a díky malému množství použitých atomů nejsou další operace výpočetně náročné. Kombinací algoritmu Matching Pursuit s Wigner-Villovou distribucí se podařilo odstranit nežádoucí výskyt interference v časově-frekvenční mapě a vlastními modifikacemi metody vznikly dva přístupy (označované jako MPWVD-S a MPWVD-E), z nichž každý našel své uplatnění v časově-frekvenční analýze. Metody se odvíjí od úspěšné detekce QRS, respektive vlny R, která slouží jako orientační bod pro výřez signálu, na němž je aplikován algoritmus MPWVD. Díky modifikaci MPWVD-S provádějící přeškálování hodnot použitých atomů je možné zobrazit podrobnou časově-frekvenční mapu se všemi důležitými detaily. To umožnilo rozměření několika významných bodů EG cyklu. Zjištění hodnot PQ a QT intervalu se prokázalo jako poměrně úspěšné během stabilizační fáze. K problematickému určení došlo až během ischemie v důsledku elevace ST segmentu, kde však se stanovením konce vlny T mají problémy i mnohé přístupy vlnkové transformace. Odladění algoritmu pro úspěšnější detekci vln by vyžadovalo hlubší prozkoumání této problematiky a stanovení citlivějších kritérií hodnocení.

53

Největším přínosem je metoda MPWVD-E, která se zaměřuje na detekci ischemie. Získané histogramy reprezentují četnost jednotlivých vlnek použitých při dekompozici vzhledem k frekvenci. Z těchto histogramů vytvořených pro různé průběhy elektrogramů je patrné, že při změně stabilizační fáze na ischemickou dochází ke změně rozložení energie. Při ischemii dochází k nárůstu zastoupení nízkých frekvencí, vyšších frekvencí ubývá. Ve fázi následné reperfuze se rozložení energie jednotlivých frekvenčních složek vrací do původního stavu. Tato změna rozložení energie byla využita při detekci ischemie. Testování na signálech s apriorní znalostí počátku ischemie přineslo pozitivní výsledky stoprocentní úspěšnosti označení signálu za ischemický. To vedlo k vytvoření referenčního modelu stabilizace, podle nějž by mohla být ischemie detekována v libovolném signále. Tento přístup má však řadu úskalí, neboť vlastnosti referenčního modelu odrážejí specifika měřeného souboru a při snaze zachovat včasnou detekci ischemie ztrácí metoda na úspěšnosti. Nejlepších výsledků bylo dosaženo třetí variací detektoru kombinujícího předchozí přístupy. Tato metoda dokáže bez jakékoli apriorní znalosti signálu vytvářet vlastní referenční model stabilizace a postupně jej adaptivně měnit. Z tabulky hodnot detekovaných časů (viz Tabulka 3) je patrné, že detekce byla ve všech případech úspěšná a velmi rychlá, neboť ve většině případů byla ischemie detekována v čase kolem 2 min, kdy ještě není zřejmá z časového průběhu. Metoda tak může konkurovat nejnovějším technikám automatické detekce ischemie. Časově-frekvenční analýza metodou MPWVD se tak ukázala být velice účinným nástrojem pro analýzu experimentálních dat a výsledky slibují dobré uplatnění i v praxi, ať už při analýze detailů naměřeného signálu offline, či formou online zpracování signálu v reálném čase.

54

Závěr V této diplomové práci byla prezentována elektrofyziologie srdce, od základních poznatků o vzniku a šíření elektrického impulsu srdcem až po možné patologické jevy včetně jejich projevů na elektrokardiografu. Byly popsány časové a spektrální vlastnosti EKG a principy jejich měření. Taktéž byly v rámci uvedení do tématiky experimentální kardiologie prezentovány přístupy získávání EG a zhodnocena podobnost EG z izolovaného srdce vzhledem k reálnému EKG člověka. V další části práci byla podrobně ilustrována problematika snímání elektrické aktivity z izolovaných srdcí a prezentována databáze signálů EG použitých v rámci této práce k testování navržených metod. Navazující část diplomové práce pojednává o časově-frekvenční analýze signálů a postupně představuje její metody. Nejvíce prostoru je věnováno spojité vlnkové transformaci, algoritmu Matching Pursuit a Wigner-Villově distribuci s podrobným matematickým rozborem. Na základě předchozích teoretických poznatků byla navržena programová aplikace pro práci s daty v programovém prostředí MATLAB. Nejprve byla na vzorcích signálu testována úspěšnost aproximace Matching Pursuit a poté přesnost detektoru R vln. Z důvodu specifických vlastností souboru EG byla navržena vlastní modifikace detekce QRS. Na základě testování se detekce projevila jako velmi přesná s hodnotami senzitivity a pozitivní prediktivní hodnoty téměř 100%. Konkrétní specifika a využitelnost detektoru QRS v obecnějším pojetí je diskutována v závěrečné části práce. Následně byly elektrogramy analyzovány metodou Matching Pursuit v kombinaci s Wigner-Villovou transformací a taktéž pomocí spojité vlnkové transformace. Na základě srovnání výsledků přínosu pro časově-frekvenční analýzu byla zhodnocena metoda využívající algoritmus Matching Pursuit jako optimálnější. Byly vytvořeny dvě modifikace této metody aplikující Matching Pursuit na Wigner-Villovu transformaci, z nichž každá našla své uplatnění na poli časově-frekvenční analýzy. První varianta MPWVD-S slouží k detekci vln P, Q a T na křivce PQRST elektrogramu. Druhá varianta MPWVD-E je testována na celém souboru EG signálů se zaměřením na detekci ischemie, a srovnávána s metodou detekce ischemie pomocí spojité vlnkové transformace. Z několika navržených přístupů metody detekce pomocí MPWVD-E je nakonec vyprofilována metoda automatické detekce ischemie fungující se 100% úspěšností bez jakékoli apriorní znalosti signálu. Velkou předností této metody je její včasná detekce, která v některých případech dosahovala rychlosti detekce do dvou minut a v drtivé většině případů do tří minut. Metoda tak může nalézt veliké uplatnění v klinické praxi, neboť dosahuje spolehlivé a mnohem rychlejší automatické detekce, než je nyní v praxi běžné. Například pro detekci infarktu myokardu je sice včasná detekce ischemie v praxi de facto nepoužitelná, neboť pacienti přicházejí na základě symptomů bolesti a srdečního selhání až několik minut či spíše hodin po 55

rozvoji ischémie, avšak nabízí se několik možností jiného využití, kde by včasná detekce mohla být velkým přínosem. Jedním z nich mohou být zátěžové testy v kardiologii, které podstupují jak pacienti s podezřením na ischemickou chorobu srdeční (ICHS), tak pacienti s již prokázanou ICHS. Během zátěžových testů je snaha vyprovokovat ischémii a její přítomnost má zásadní význam pro další terapeutické rozhodování a prognostické zvažování, proto by zde mohl algoritmus MPWVD najít své velké uplatnění. Další možné využití detekce ischemie se nabízí na jednotkách intenzivní péče (JIP), kde je srdeční aktivita pacientů stále monitorována. Zde by v případě automatického detektoru napojeného na snímání EKG mohla být včasná detekce užitečná k oznámení blížícího se akutního stavu. Velký přínos by mohla navržená metoda nabídnout v oblasti telemedicíny, kdy je pacientův EKG monitorován na dálku. U takto monitorovaných pacientů je často očekáván možný výskyt ischemie například v souvislosti s prodělaným infarktem myokardu, a díky včasné detekci ischemie a pohotovému řešení situace by se dalo předejít vážnějším komplikacím a život ohrožujícím stavům.

56

Seznam literatury [1]

[2]

[3] [4] [5] [6] [7] [8]

[9]

[10]

[11]

[12] [13]

G.MALLAT, Stephane a Zhifeng ZHANG. Matching Pursuits With Time-Frequence Dictionaries. IEEE transactions on signal processing. 1993, roč. 41, č. 12, s. 3397– 3415. School of Health Sciences: Practice Learning Resources. The University of Nothingham. The University of Nottingham: Cardiology Teaching Package [online]. 2000 [cit. 2015-04-18]. Dostupné z: http://www.nottingham.ac.uk/nursing/ practice/resources/cardiology/function/ WILHELM, Zdeněk. Stručný přehled fyziologie člověka pro bakalářské studijní programy. 4. vyd. Brno: Masarykova univerzita, 2010, 117 s. ISBN 978-80-210-5283-3. BIENERTOVÁ VAŠKŮ, Julie. Praktikum z patologické fyziologie. Brno: kolektiv pracovníků Ústavu patologické fyziologie LF MU, 2007. ISSN 1802-128X. SILBERNAG, Stefan a Agamemnon DESPOPOULOS. Atlas fyziologie člověka. 6. přeprac. a rozš. vyd. Praha: Grada, 2004, XII, 435 s. ISBN 80-247-0630-X. GANONG, William F. Přehled lékařské fyziologie. 20. vyd. Praha: Galén, 2005, 890 s. ISBN 80-726-2311-7. ROZMAN, J. Elektronické přístroje v lékařství. In Elektronické přístroje v lékařství. Academia. Praha: Academia, 2006. s. 1-432. ISBN: 80-200-1308- 3. WATFORD, Christopher. Understanding ECG Filtering. EMS 12 Lead: Cardiac Rhythm Analysis, 12-Lead ECG Interpretation, Resuscitation [online]. 2014 [cit. 201505-20]. Dostupné z: http://www.ems12lead.com/2014/03/10/understanding-ecgfiltering/ POKORNÝ, Jan. Elektrofyziologie [online]. Praha, 2010, 374 s. [cit. 2015-05-20]. České vysoké učení technické v Praze. Dostupné z: http://fbmi.cvut.cz/files/nodes/657/ public/Elektrofyziologie.pdf Myocardial Ischemia due to Coronary Artery Blockage. Celebrate Life, Unlimited. [online]. 2011 [cit. 2015-05-20]. Dostupné z: https://liveimmensely.wordpress.com/ 2011/08/27/myocardial-ischemia-due-to-coronary-artery-blockage/ CHANNER, K. ABC of clinical electrocardiography: Myocardial ischaemia. BMJ. 2002, roč. 324, č. 7344: s. 1023-1026. DOI: 10.1136/bmj.324.7344.1023. ISSN 09598138. Dostupné také z: http://www.bmj.com/cgi/doi/10.1136/bmj.324.7344.1023 Haman, P., Výukový web EKG [online]. [cit. 2015-05-20]. Dostupné z: http://ekg.kvalitne.cz/ Patofyziologie ischemické srdeční choroby. Ústav Patologické fyziologie [online]. 2012 [cit. 2015-05-20]. Dostupné z: http://pfyziollfup.upol.cz/castwiki2/?p=1037

57

[14]

[15]

[16] [17]

[18] [19] [20]

[21] [22]

[23]

[24]

[25]

SÖRNMO, Leif a Pablo LAGUNA. Bioelectrical signal processing in cardiac and neurological applications. Boston: Elsevier Academic Press, 2005, 668 s. ISBN 01-2437552-9. HAMILTON, P. Open source ECG analysis. Computers in Cardiology [online]. IEEE, 2002, roč. 29, s. 101-104 [cit. 2015-04-18]. DOI: 10.1109/CIC.2002.1166717. Dostupné z: http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1166717 WAGNER, Joseph E a Patrick J MANNING. The Biology of the guinea pig. New York: Academic Press, 1976, 317 s. ISBN 01-273-0050-3. BELL, Robert M., Mihaela M. MOCANU a Derek M. YELLON. Retrograde heart perfusion: The Langendorff technique of isolated heart perfusion. Journal of Molecular and Cellular Cardiology. 2011, roč. 50, č. 6, 940-950 s. DOI: 10.1016/j.yjmcc.2011.02.018. ISSN 00222828. Dostupné také z: http://linkinghub.elsevier.com/retrieve/pii/S0022282811000952 MAYER, Jörg a Thomas M DONNELLY. Clinical veterinary advisor. Missouri: Elsevier, 2013, 752 s. ISBN 978-141-6039-693. REECE, William O. 2011. Fyziologie a funkční anatomie domácích zvířat. 1. české vyd. Praha: Grada, 473 s. ISBN 978-80-247-3282-4. BOERMA, Marjan, Junru WANG, Vijayalakshmi SRIDHARAN, Jean-Marc HERBERT, Martin HAUER-JENSEN a Gayle E. WOLOSCHAK. Pharmacological Induction of Transforming Growth Factor-Beta1 in Rat Models Enhances Radiation Injury in the Intestine and the Heart. PLoS ONE [online]. 2013, roč. 8, č. 7, [cit. 201504-18]. DOI: 10.1371/journal.pone.0070479. Dostupné z: http://dx.plos.org/10.1371/journal.pone.0070479 KOLÁŘOVÁ, J., Optické mapování elektrické aktivity izolovaného srdce. Brno, 2009. Veřejná přednáška konaná v rámci habilitačního řízení. Vysoké učení technické v Brně. NOVÁKOVÁ, Marie, Jiří MOUDR a Pavel BRAVENÝ. A modified perfusion system for pharmacological studies in isolated hearts. In Analysis of Biomedical Signals and Images, Biosignal 2000. 15. vyd. Brno: Vutium Press, Brno University of Technology, 2000. s. 162-164. ISBN 80-214-1610-6. PROVAZNÍK, Ivo. Experimentální metody získávání, zpracování a analýzy biologických signálů: Experimental methods biological signals of acquisition, processing, and analysis : teze přednášky k profesorskému jmenovacímu řízení v oboru Elektronika a sdělovací technika [online]. Brno: VUTIUM, 2005.ISBN 80-214-2933-x. LANGENDORFF, O. Untersuchugen am überlebenden Säugethierherzen, Pflügers Archiv für die gesamte Physiologie des Menschen und der Tiere, roč. 65, č. 1, s. 291332, 1895. NOVÁKOVÁ, Marie. Effects of sigma receptor ligand BD737 in rat isolated hearts. Scripta medica, Brno: Masaryk University, Faculty of Medicine, 2007, roč. 80, č. 6, s. 255-262. ISSN 1211-3395 58

[26]

[27]

[28]

[29]

[30] [31]

[32] [33]

[34]

[35]

SUTHERLAND, FIONA J. a DAVID J. HEARSE. The Isolated Blood And Perfusion fluid Perfused Heart. Pharmacological Research. 2000, roč. 41, č. 6, s. 613-627. DOI: 10.1006/phrs.1999.0653. ISSN 10436618. Dostupné také z: http://linkinghub.elsevier.com/retrieve/pii/S1043661899906539 GRANADO Miriam, Nuria FERNÁNDEZ, Luis MONGE, Juan Carlos FIGUERAS, Gonzalo CARREÑO-TARRAGONA, Sara AMOR, Angel Luis GARCÍA-VILLALÓN a Tamas ZAKAR. Effects of Coronary Ischemia-Reperfusion in a Rat Model of Early Overnutrition. Role of Angiotensin Receptors. PLoS ONE [online]. 2013, roč. 8, č. 2, [cit. 2015-04-18]. DOI: 10.1371/journal.pone.0054984. Dostupné z: http://dx.plos.org/10.1371/journal.pone.0054984 JAKSCH, I. Časově frekvenční analýza signálů [online]. Liberec: Technická univerzita v Liberci, 2010, 27 s. [cit. 2015-05-20]. Dostupné z: http://www.rss.tul.cz/download/tdg/P6_cas_freq_analyza.pdf KOZUMPLÍK, Jiří. Vlnkové transformace a jejich využití pro filtraci signálů EKG: Wavelet transforms and their use for filtering of ECG signals. Brno: VUTIUM, 2005. 27 s. ISBN 80-214-3045-1. Zkrácená verze habilitační práce. Vysoké učení technické v Brně. KOZUMPLÍK, J. Multitaktní systémy. Elektronická skripta. Brno, FEKT VUT v Brně. 2005. p. 1 - 57. ŠMÍD, Radislav. Úvod do vlnkové transformace [online]. ČVUT FEL, 2001 [cit. 201504-18]. Dostupné z: http://measure.feld.cvut.cz/groups/diag/download/Waveletintro8859.pdf DURKA, P. Matching pursuit and unification i EEG analysis. Boston: Artech House, 2007, 184 s. ISBN 978-1-58053-304-1. MAUTNER, P., MOUČEK, R. Využití metody matching pursuit pro detekci ERP vln. In Kognice a umělý život. Opava: Slezská univerzita, 2009. s. 201-205. ISBN: 978-807248-516-1 CLIFFORD, Gari D, Francisco AZUAJE a Patrick MCSHARRY. Advanced methods and tools for ECG data analysis. Boston: Artech House, 2006, 384 s. Artech House engineering in medicine. ISBN 1580539661. MOCHIZUKI, Seibu (ed.). The Ischemic heart: Progress in Experimental Cardiology. S.l.: Springer-Verlag New York, 2013. ISBN 978-1475783001.

59

SEZNAM ZKRATEK ATP

adenosintrifosfát (z angl. adenosine triphosphate)

AV

atrio-ventrikulární

CWT

spojitá vlnková transformace (z angl. Continuous Wavelet Transform)

DTWT

vlnková transformace s diskrétním časem (z angl. Discrete-Time Wavelet Transform)

DWT

diskrétní vlnková transformace (z angl. Discrete Wavelet Transform)

EG

elektrogram

EKG

elektrokardiogram

FN

počet falešně negativních detekcí (z angl. False Negative)

FP

počet falešně pozitivních detekcí (z angl. False Positive)

HR

srdeční rytmus (z angl. Heart Rate)

IM

infarkt myokardu

LTI

lineární časově invariantní (z angl. Linear Time-Invariant)

MP

Matching Pursuit

MPWVD

Matching Pursuit a Wigner-Ville Distribution algoritmus

MRA

analýza s násobným rozlišením (z angl. Multi Resolution Analysis)

P+

pozitivní prediktivita (z angl. Positive Predictivity)

PQ

interval mezi začátkem depolarizace síňí a začátkem depolarizace komor

PSD

výkonové spektrum signálu (z angl. Power Spectral Density)

QT

interval mezi začátkem depolarizace a koncem repolarizace komor

RR

interval mezi dvěma po sobě jdoucími vlnami R

SA

sino-atriální

Se

senzitivita detekce (z angl. Sensitivity)

STFT

krátkodobá Fourierova transformace (z angl. Short-Time Fourier Transform)

TP

počet pravdivě pozitivních detekcí (z angl. True Positive)

UBMI

Ústav biomedicínského inženýrství

WT

vlnková transformace (z angl. Wavelet Transform)

WVD

Wigner-Ville Distribution 60

PŘÍLOHY Popis vytvořených funkcí a informace ke spuštění algoritmu Všechny implementované funkce byly realizovány v programovém prostředí MATLAB verze R2012b. Z důvodu plné kompatibility je doporučeno spouštět algoritmus ve shodné verzi. Funkce lze spouštět přímo z příkazového řádku, avšak pro lepší orientaci a plynulou návaznost byl vytvořen skript DP.m, kde jsou všechny funkce uvedeny spolu s předvolenými hodnotami vstupů. Hodnoty je tak možné jednoduše měnit a buď spouštět postupně jednotlivé bloky či celý skript najednou. U několika funkcí je volba proměnných realizována tak, že se po spuštění funkce objeví v příkazovém řádku výzva k vložení hodnoty proměnné spolu s doporučeným rozmezím. Také je možné využít nápovědy naprogramovaných funkcí, která se zobrazí pomocí příkazu help a jména funkce. Následuje výčet všech funkcí s krátkým popisným komentářem.

nacteniEG.m

Načte signál ve formátu *.bin ze zvolené složky a uloží do proměnné ve formátu *.mat. Je zde také volena vzorkovací frekvence a případně body počátku ischemie a počátku reperfuze.

vykresliSignal.m

Vykreslí celý načtený signál s označením hranic ischemie a reperfuze.

vyberEG.m

Provádí detekci QRS a dle detekované R vlny indexuje ze signálu jednotlivé PQRST úseky, které ukládá do matice EG

EGzaMinutu.m

Zobrazí průměrné EG v každé minutě měření

zobraz5min.m

Zobrazí prvních 5 sekund na začátku každé minuty z prvních pěti minut ischemie

slovnik_img.m

Vykresluje jednotlivé atomy Gaborova slovníku

prubeh_stab_isch_rep.m

Funkce prezentuje rozdílnost průběhů stabilizace ischemie a reperfuze na vybraném úseku EG signálu

________________________________________________________________________________ MatchingPursuit.m

Provádí základní algoritmus Matching Pursuit a prezentuje úspěšnost jeho aproximace na části signálu.

WignerVilleDistribution.m

Počítá časově-frekvenční mapu signálu pomocí Wigner-Villovy distribuce.

MP_WVD.m

Kombinuje metodu Matching Pursuit a Wigner-Villovu distribuci a počítá časově-frekvenční mapy pro jednotlivé atomy získané dekompozicí signálu. Jejich sloučením vzniká výsledná časověfrekvenční mapa daného EG cyklu (PQRST úseku).

MP_WVD_komplet.m

Aplikuje funkci MP_WVD.m na dlouhý úsek signálu a ukládá jednotlivé časově-frekvenční mapy do 3D matice

61

krivkaAnimace.m

Plynule vykreslí jednotlivé indexované EG úseky v rychlém sledu za sebou a vytvoří z nich video, které je pak možné přehrát a případně pomalým krokováním sledovat změny v analyzovaném signálu.

ukazka_MP.m

Prezentuje názorně princip skládání časově-frekvenční mapy z příspěvků jednotlivých atomů.

rozmereniPQT.m

Provádí rozměření signálu na základě MPWVD-S algoritmu. Určuje začátek P a Q vlny a konce vlny T.

MPhistogram.m

Realizuje výpočet histogramu ze získané časově-frekvenční mapy.

st_is_re.m

Zobrazuje srovnání časově-frekvenčních map jednotlivých fází měření (stabilizace, ischemie a reperfuze)

srovnej1.m

Zobrazuje vývoj změn časově-frekvenční mapy a signálu v rámci měření.

srovnej2.m

Zobrazuje srovnání dvou časově-frekvenčních map.

srovnej3.m

Zobrazuje srovnání tří časově-frekvenčních map. Například porovnává výsledky různých přístupů časově-frekvenční analýzy (Matching-Pursuit, spojité vlnkové transformace a krátkodobé Fourierovy transformace).

interpolHist.m

Interpoluje hodnoty histogramu na požadovanou délku pro možné porovnání hodnot různě dlouhých histogramů.

________________________________________________________________________________ WVDimage.m

Podpůrná funkce zobrazující časově-frekvenční mapu.

WVDi.m

Podpůrná funkce zobrazující tři různé časově-frekvenční mapy.

WVDi3.m

Podpůrná funkce zobrazující tři různé časově-frekvenční mapy různých metod časově-frekvenční analýzy.

WVDplay.m

Podpůrná funkce zobrazující vývoj změn časově-frekvenční mapy a signálu v rámci měření.

________________________________________________________________________________ detekce_isch_referHist.m

Provádí automatickou detekci ischemie na základě hodnocení histogramu MPWVD a jeho srovnání s předem známým referenčním modelem.

detekce_isch_MPWVD.m

Provádí automatickou detekci ischemie na základě histogramů MPWVD avšak funkce si sama vytváří adaptivní referenční model.

detekce_isch_CWT.m

Provádí automatickou detekci ischemie na základě spojité vlnkové transformace a hodnocení jejích histogramů.

62