UNIVERSITAT LINZJOHANNES KEPLER JKU

Technisch-NaturwissenschaftlicheFakultat

Automatisierte Analyseverfahren furElektrokardiogramme

BACHELORARBEIT

zur Erlangung des akademischen Grades

Bachelor of Science

im Bachelorstudium

Mechatronik

Eingereicht von:

Hannes Browa

Angefertigt am:

Institut fur Signalverarbeitung

Beurteilung:

Univ.-Prof. Dr. Mario Huemer

Betreuung:

DI Alexander MelzerDI Oliver LangDI (FH) Christoph Mahringer

Linz, 9. Februar 2016

Inhaltsverzeichnis

1 Einfuhrung 31.1 Allgemeines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Ziele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Uberblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Das Elektrokardiogramm (EKG ) . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4.2 Nomenklatur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Grundlagen der Wavelet-Transformation 52.1 Motivation der Wavelet-Transformation . . . . . . . . . . . . . . . . . . . . . 52.2 Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2.1 Das Haar-Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.2 Das Mexican-Hat-Wavelet . . . . . . . . . . . . . . . . . . . . . . . . 62.2.3 Das Morlet-Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.4 Daubechies-Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3 Wichtige Eigenschaften von Wavelets . . . . . . . . . . . . . . . . . . . . . . 82.4 Die kontinuierliche Wavelet-Transformation . . . . . . . . . . . . . . . . . . . 9

2.4.1 Variation der CWT-Parameter . . . . . . . . . . . . . . . . . . . . . . 92.5 Die diskrete Wavelet-Transformation . . . . . . . . . . . . . . . . . . . . . . 102.6 Skalierungsfunktion und Multiskalenanalyse . . . . . . . . . . . . . . . . . . . 112.7 Beispiele der Wavelet-Transformation . . . . . . . . . . . . . . . . . . . . . . 12

2.7.1 Einfaches Testsignal: Rechteck- und e-Funktion . . . . . . . . . . . . . 122.7.2 Analyse der Momentanfrequenz . . . . . . . . . . . . . . . . . . . . . 142.7.3 Ein EKG-Signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Implementierung der EKG-Zeitbereichsanalyse 173.1 Detektion der QRS-Komplexe . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.1.1 Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.1.2 Signale der Einzelstufen . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2 Abschatzung der Herzfrequenz . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.3.1 QRS-Unterdruckung . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3.2 Weitere Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.4 Lokalisierung der Einzelwellen (Kovarianzfilter) . . . . . . . . . . . . . . . . . 223.5 Approximation durch Gaußkurven-Segmente . . . . . . . . . . . . . . . . . . 223.6 Vergleich mit Referenz-Markern . . . . . . . . . . . . . . . . . . . . . . . . . 253.7 Statistische Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Grafische Benutzeroberflache (GUI) 284.1 Idee und Funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.2 Umsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5 Performance-Bewertung 305.1 Synthetischer Datensatz (rauschfrei) . . . . . . . . . . . . . . . . . . . . . . 30

5.1.1 Herzfrequenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305.1.2 Detektionsrate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305.1.3 Standardabweichung . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.1.4 Mittelwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5.2 Synthetischer Datensatz (stark verrauscht) . . . . . . . . . . . . . . . . . . . 325.2.1 Herzfrequenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.2.2 Detektionsrate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.2.3 Standardabweichung . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.2.4 Mittelwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.3 Gegenuberstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345.4 PhysioNet QT-Datenbank . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

6 Fazit 37

7 Ausblick 37

Abbildungsverzeichnis 39

Tabellenverzeichnis 39

1

Abstract

The aim of this bachelor’s thesis is to investigate methods of automated analysis ofelectrocardiogram (ECG) data with focus on time-domain based algorithms.As a first part the thesis provides an introduction to electrocardiogram signals as well astheory and examples of the wavelet analysis, which has established a reputation as state-of-the-art method in large fields of digital signal processing. The thesis then illustratesa plain ECG analysis algorithm based on simple time domain filters and further outlinesadvantages and drawbacks of the implementation. It’s main objective is to perform thedetection and approximation of particular signal features as accurately and reliably aspossible. The reliability of the algorithm is then examined by applying it on a series ofsynthetically generated reference ECG data sets and genuine patient data from a publicmedical database.In conclusion, this thesis provides a comprehensive overview of digital electrocardiogramsignal processing and provides a basis for the step to a wavelet-transform algorithm,which is planned for future research.

2

1 Einfuhrung

1.1 Allgemeines

In der modernen Medizin hat sich die Messung von Elektrokardiogrammen, kurz EKG, als wich-tiges Diagnoseverfahren der Herzaktivitat etabliert. Die Analyse der großen Datenmengen, diedabei generiert werden, bedarf jedoch eines immensen Zeitaufwandes speziell geschulten Fach-personals, weshalb es unabdingbar ist, zumindest einen Teil der Auswertung auf automatisierte(rechnergestutze) Methoden auszulagern. Diese Bachelorarbeit befasst sich mit der digitalisier-ten Verarbeitung von EKG-Signalen, wobei der Fokus auf einfache und verstandliche Methodenim Zeitbereich gelegt wurde. Oberste Prioritat hat die hochstmogliche Zuverlassigkeit der be-schriebenen Algorithmen sowie eine intuitive Implementierung. Damit ist der Grundstein furVergleiche mit anderen Methoden wie z. B. der Wavelet-Analyse gelegt, die fur diese Anwen-dung als leistungsstarkes Instrument der digitalen Signalverarbeitung in Frage kommt.

1.2 Ziele

Mit dieser Bachelorarbeit soll eine Grundlage fur die automatisierte Analyse von Elektrokardio-grammen geschaffen werden. Hierzu werden zunachst robuste und relativ einfache Zeitbereichs-methoden prasentiert, welche im Zuge zukunftiger Forschung durch andere Methoden ersetztwerden konnen. Zukunftiges Hauptziel ist hier die Anwendung der Wavelet-Transformation,deren Grundlagen diese Arbeit daher ebenfalls umreißt. Außerdem wird die Performance derbeschriebenen Routinen bewertet und mit Ergebnissen anderer Arbeiten verglichen.

1.3 Uberblick

Eine kurze Erlauterung des EKG-Signals ist in Kapitel 1.4 angefuhrt, diese ist lediglich alsUberblick ohne detailliertes medizinisches Fachwissen anzusehen. Kapitel 2 gibt einen Einblickin die Wavelet-Transformation, deren mathematische Grundlagen sowie einige Anwendungs-beispiele. Der implementierte Signalverarbeitungs-Algorithmus wird in Kapitel 3 ausfuhrlichbeschrieben, die dazugehorige grafische Benutzeroberflache wird in Kapitel 4 erlautert. Kapi-tel 5 beinhaltet die Analyse der Ergebnisse.

1.4 Das Elektrokardiogramm (EKG)

1.4.1 Definition

Das Elektrokardiogramm, kurz EKG (engl. ECG ), dient der Analyse der Herzaktivitat. Dazuwerden Elektroden nach standardisierten Mustern1 am Korper angebracht. Diese erfassen dieelektrischen Potentiale, welche durch die Zellaktivitaten im Herzen entstehen. Dadurch ist esnicht nur moglich, die Herzfrequenz exakt zu registrieren; anhand der Breite einzelner Signal-merkmale konnen auch Ruckschlusse auf die Aktivitaten verschiedener Abschnitte des Herzengezogen werden. Somit ist eine Fruherkennung sowie eine genaue Uberwachung typischerKrankheitsbilder moglich. Mit dem großen Vorteil, dass lediglich Elektroden auf der Hautaufgebracht werden mussen, und keine weiteren kostspieligen oder gefahrlichen medizinischenEingriffe notwendig sind, hat sich das EKG als essenzielles kardiologisches Diagnoseverfahrendurchgesetzt [2]. Typische Einsatzgebiete sind Ruhe-EKG, Belastungs-EKG und Langzeit-EKG.

1Die Platzierung und Anzahl der Elektroden hangt davon ab, wieviele und welche EKG-Kanale, so-genannte Ableitungen, aufgezeichnet werden sollen.

3

1.4.2 Nomenklatur

Den einzelnen Teilwellen (Zacken) eines Herzzyklus sind zur Unterscheidung mit den Buch-staben P bis T (manchmal auch U) gekennzeichnet. Man spricht von P-Welle, T-Welle sowieQRS-Komplex, da bei letzterem die drei Teilwellen unmittelbar ineinander ubergehen. Auchdie Zeitdauer zwischen An- oder Abklingen zweier Wellen kann von Interesse sein - relevantsind vor allem PQ-Strecke, ST-Strecke, PR- und QT-Intervall [3] (siehe Abb. 1.1).Die P-Welle entspricht der Vorhoferregung, der QRS-Komplex der Kammererregung und dieT-Welle der Erregungsruckbildung. Genaueres medizinisches Verstandnis ist fur die hier be-schriebenen Methoden der Signalverarbeitung nicht von Noten, weshalb an dieser Stelle keineweiteren Details erlautert werden.

P− QRS− T−Welle Komplex Welle

PR− QT−Intervall Intervall

Abbildung 1.1: Nomenklatur des EKG-Signals2

2Synthetisches EKG-Signal des Generators von McSharry und Clifford [11]

4

2 Grundlagen der Wavelet-Transformation

2.1 Motivation der Wavelet-Transformation

Als Standard der Signalanalyse im Frequenzbereich hat sich in vielen Anwendungsgebietendie Fourier-Transformation etabliert. Einerseits ist das Frequenzspektrum intuitiv “lesbar”, an-dererseits kann z. B. die schnelle Fourier-Transformation (engl. fast Fourier transform, kurzFFT ) sehr recheneffizient implementiert werden. Aufgrund der unendlichen Ausdehnung desFourier-Integrals konnen jedoch nur uber das gesamte Signal gemittelte Spektralinformationengewonnen werden, wodurch es fur die Beschreibung transienter Signalformen zu einem spe-zifizierten Zeitpunkt sehr schlecht geeignet ist. Um die zeitliche Veranderung des Spektrumsdarstellen zu konnen, kann die Kurzzeit-Fourier-Transformation (engl. short-time Fourier trans-form, kurz STFT ) angewendet werden, bei der ein bewegtes Zeitfenster von festgelegter Langezum Einsatz kommt.Eine wesentlich flexiblere Moglichkeit der lokalen Spektralanalyse bietet die Wavelet-Transfor-mation, welche Fenster mit variabler Breite nutzt. Somit konnen gleichzeitig langsame (nie-derfrequente) und schnelle (hochfrequente) Signalanteile erfasst und zeitlich genau lokalisiertwerden. Dadurch eignet sich diese Methode hervorragend fur die Analyse transienter bzw.aperiodischer Signale, bei denen gewisse Frequenzanteile zu einem bestimmten Zeitpunkt vonInteresse sind. Ein weiterer Vorteil der Wavelet-Transformation liegt darin, dass eine Vielzahlmoglicher Tragerfunktionen (Wavelets) zur Verfugung steht. Je nach Signalform kann dar-aus die passendste ausgewahlt werden, im Gegensatz zur Fourier-Transformation, die auf denSinus-Trager beschrankt ist. Unter Verwendung einer orthogonalen Basis eignet sich die diskre-te Wavelet-Transformation ebenfalls zur Codierung bzw. Kompression, wie etwa die diskreteCosinus-Transformation.

2.2 Wavelets

Als Grundlage der Wavelet-Transformation sei im Folgenden eine Kurvenschar ψa,b definiert,die mit den zwei Parametern a und b eine Wavelet-Funktion ψ(t) skaliert und verschiebt

ψa,b(t) =1√|a|

ψ

(t− ba

), a 6= 0, b ∈ R. (2.1)

Dabei ist a der Dilatations- bzw. Dehnungsfaktor und b der Translations- bzw. Verschiebungs-parameter. Die unveranderte Funktion ψ(t) bezeichnet man auch als Mutter-Wavelet (engl.mother wavelet).

2.2.1 Das Haar-Wavelet

Das alteste bekannte Wavelet wurde 1909 von Alfred Haar eingefuhrt [4]. Es ist zugleichdie einfachste bekannte Wavelet-Funktion, da sie sich lediglich aus zwei Rechteckfunktionengemaß

ψ(t) =

1 0 ≤ t <

1

2

−11

2≤ t < 1

0 sonst.

(2.2)

zusammensetzt. Damit ist die Implementierbarkeit extrem einfach - allerdings ist die Funktionunstetig und damit nicht differenzierbar (Abb. 2.1a).

5

−1.5 −1 −0.5 0 0.5 1 1.5

−1

−0.5

0

0.5

1

t

ψ(t

)

Haar−Wavelet

−4 −2 0 2 4

−0.5

0

0.5

1

t

ψ(t

)

Mexican−Hat−WaveletStandardnormalverteilung

Abbildung 2.1: (a) Haar-Wavelet, (b) Ricker-Wavelet (“Mexican hat”)

2.2.2 Das Mexican-Hat-Wavelet

Als Beispiel fur eine kontinuierliche Wavelet-Funktion sei hier das gebrauchliche Ricker-Waveletangefuhrt - aufgrund seiner Sombrero-artigen Form wird es in der englischsprachigen Literaturmeistens als “Mexican hat” bezeichnet, was auch im deutschsprachigen Raum ubernommenwurde. Es wird durch

ψ(t) = (1− t2) e−t2

2 (2.3)

beschrieben, was der negativen zweiten Ableitung der Gauß’schen Standardnormalverteilung3

entspricht. Streng genommen musste die Ableitung ebenso deren Normierungsfaktor 1/√

2πaufweisen, der aber in der verbreiteten Definition des Ricker-Wavelets vernachlassigt wird. DasWavelet sowie die zugrunde liegende Normalverteilung sind in Abb. 2.1b dargestellt.Grundsatzlich liefern alle Ableitungen der Gauß-Funktion gultige Wavelets, die Auswahl hangtvon der Anwendung ab. Am gebrauchlichsten ist die zweite Ableitung, manchmal kommt auchdie erste zum Einsatz. Hohere Ableitungen findet man in der Praxis eher selten. [7]

2.2.3 Das Morlet-Wavelet

Es existieren auch komplexwertige Wavelets. Eines der gebrauchlichsten ist das Morlet-Wavelet.Es besteht aus einer komplexen Tragerschwingung mit der Frequenz ω0 multipliziert mit einereinhullenden Standardnormalverteilungsfunktion. Da D. Gabor bereits 1946 postulierte, dassdies der beste Kompromiss aus Frequenz- und zeitlicher Auflosung sei [5], wird die Funktionauch als Gabor-Wavelet bezeichnet. Das Gabor-Wavelet ist gemaß

ψ(t) =14√π

(eiω0t − eω

20t)

e−t2

2 (2.4)

definiert. Diese Definition beinhaltet wegen der Admissibilitatsbedingung (siehe Kapitel 2.3)einen Korrekturterm eω

20t. Dieser kann fur ω0 > 5 vernachlassigt werden was zum vereinfachten

Morlet-Wavelet fuhrt

ψ(t) =14√π

eiω0t e−t2

2 (2.5)

. Abb. 2.2 zeigt das Wavelet in reeller und komplexer Darstellung.

3Manche Autoren geben die Definition auch allgemeiner fur beliebige Werte von σ an.

6

−4 −2 0 2 4

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

t

ψ(t

)

Reelles Morlet−Wavelet ω0=5

−5

0

5

−1

0

1−1

−0.5

0

0.5

1

tRe(ψ(t))

Im(ψ

(t))

Komplexes Morlet−Wavelet ω0=5

Abbildung 2.2: Das Morlet-Wavelet in reeller und komplexer Darstellung

2.2.4 Daubechies-Wavelets

1988 fuhrte Ingrid Daubechies erstmals stetig differenzierbare orthogonale Wavelet-Funktionenmit kompakten Tragern ein. Diese weisen eine maximale Anzahl N verschwindender Momen-te fur eine gegebene Tragerlange auf [6]. Bei der Implementierung ist diese Eigenschaft vongroßer Bedeutung, da die Filterordnung 2N die Lange des Tragers bestimmt. Die Daubechies-Wavelets werden ublicherweise nach ihrer (meist geradzahligen) Ordnung N als “dbN” be-zeichnet, z.B. db2, db4, db20. Eine andere gebrauchliche Form ist die Bezeichnung D mit derFilterordnung 2N , so ist z. B. “D4” gleichbedeutend mit “db2”. Einen Sonderfall stellt dieniedrigste Ordnung dar: db1 gleicht dem Haar -Wavelet. Fur hohere Ordnungen existiert keinexpliziter Ausdruck. Statt einer ausfuhrlichen Herleitung sei an dieser Stelle auf die MATLAB-Dokumentation [1] verwiesen.

Eigenschaften der dbN -Wavelets:

� Ordnung N ∈ N+

� Die Wavelet-Funktion ψ hat N verschwindende Momente

� Die Lange des Tragers von Wavelet ψ und Skalierungsfunktion φ betragt 2N − 1

� Die Filterlange betragt 2N

� Daubechies-Wavelets sind hochgradig unsymmetrisch

7

0 1 2 3

−1

−0.5

0

0.5

1

1.5

db20 2 4 6

−0.5

0

0.5

1

db40 5 10

−1

−0.5

0

0.5

1

db6

0 5 10 15

−1

−0.5

0

0.5

db8

0 10 20

−0.5

0

0.5

db12

0 10 20 30

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

db20

Abbildung 2.3: Daubechies-Wavelets

2.3 Wichtige Eigenschaften von Wavelets

Um als Wavelet klassifiziert zu werden, muss eine Funktion folgende mathematische Kriterienerfullen:

1. Beschranktheit

E =

∫ ∞−∞|ψ(t)|2 dt = 〈ψ,ψ〉 <∞ (2.6)

Die Energie (L2-Norm) des Wavelets muss endlich sein.

2. AdmissibilitatEs sei Ψ(ω) = F(ψ(t)) die Fourier-Transformierte der Funktion ψ(t), dann gilt

Cψ =

∫ ∞−∞

|Ψ(ω)|2

ωdω <∞. (2.7)

Damit diese Funktion bei w=0 definiert ist, muss gelten

Ψ(0) =

∫ ∞−∞

ψ(t) dt = 0. (2.8)

Daraus folgt, dass ein Wavelet mittelwertfrei sein muss, d. h. der GleichspannungsanteilΨ(0) ist null (Glg. 2.8). Cψ bezeichnet man als Admissibilitatskonstante, die auch furdie Rucktransformation benotigt wird. Man spricht von einem normierten Wavelet, wennCψ = 1 gilt.

3. Die Fourier-Transformierte von komplexen (oder analytischen) Wavelets muss reell seinund darf keine negativen Frequenzanteile enthalten. [7]

Die Einhaltung der angefuhrten Bedingungen stellt sicher, dass die Transformation und ihreInverse berechenbar sind.

8

2.4 Die kontinuierliche Wavelet-Transformation

Die kontinuierliche Wavelet-Transformation (engl. continuous wavelet transform, CWT ) isteine Zeit-Frequenz-Analysemethode, die im Unterschied zur herkommlichen STFT eine nahe-zu beliebig genaue4 zeitliche Lokalisierung hochfrequenter Signalanteile ermoglicht. Das liegtan der variablen Fensterbreite der CWT, die mit dem Beobachtungsmaßstab skaliert [8]. DerBegriff kontinuierlich bezieht sich dabei lediglich auf die Variation der Transformationspa-rameter, wie spater in diesem Kapitel beschrieben. Die erlauterten Berechnungsvorschriftenfur kontinuierliche Signale gelten jedoch i. d. R. ohne Einschrankungen genauso fur diskrete(abgetastete) Signale, wie sie ublicherweise in der digitalen Signalverarbeitung vorliegen. DieQualitat der Approximation5 ist hierbei maßgeblich fur den resultierenden Fehler.Die Berechnung der kontinuierlichen Wavelet-Transformation eines beliebigen Zeitsignals x(t)erfolgt durch

T (a, b) =1√a

∫ ∞−∞

x(t)ψ∗a,b dt ∼= 〈x,ψa,b〉 , (2.9)

was dem Skalarprodukt des Zeitsignals x(t) mit dem um b zeitlich verschobenen und mit demFaktor a skalieren Wavelet entspricht. Durch Einsetzen von (2.1) erhalt man

T (a, b) =1√a

∫ ∞−∞

x(t)ψ∗(t− ba

)dt, (2.10)

wobei ψ∗ die konjugiert komplexe Funktion des Wavelets ψ bezeichnet. Aus der Linearitat desSkalarproduktes folgt auch die Linearitat der Wavelet-Transformation.Wie bei der Fourier-Transformation kann das ursprungliche Zeitsignal durch eine Rucktransfor-mation rekonstruiert werden. Die inverse kontinuierliche Wavelet-Transformation (engl. inversecontinuous wavelet transform (ICWT ) ist definiert als

x(t) =1

Cψ

∫ ∞−∞

∫ ∞0

T (a, b)ψa,b(t)da db

a2. (2.11)

Die in (2.7) ermittelte Admissibilitatskonstante Cψ skaliert die Rucktransformation, sofernkeine normierte Wavelet-Funktion verwendet wird.

2.4.1 Variation der CWT-Parameter

Aus der Definition der CWT ergibt sich, dass die Parameter a und b bei kontinuierlichen Signa-len ebenfalls beliebig kontinuierlich variiert werden konnen. Dabei gilt fur die Zeitverschiebung−∞ < b <∞ und fur die Skalierung 0 < a <∞.In der Praxis wird bei der CWT eines abgetasteten Zeitsignals die Verschiebung b mit demAbtastintervall variiert und die Skalierung a im zu untersuchenden Frequenzband logarith-misch diskretisiert [8]. Ein einfaches Anwendungsbeispiel der CWT auf einen diskreten Di-racimpuls bei t = 0 fur das Mexican-Hat- und das Daubechies-db6-Wavelet ist in Abb. 2.4illustriert. In diesem Beispiel variiert die Verschiebung b in 1-Sample-Schritten im Bereich−256 ≤ b ≤ 256, also insgesamt Nb = 513 Zeitpunkte. Die Skalierung a steigt logarithmischvon a = 1 bis a = 64, wobei insgesamt Na = 301 Maßstabe ausgewertet werden. Das ProduktNa · Nb ergibt die gesamte Anzahl der CWT-Auswertungspunkte NCWT , in diesem BeispielNa ·Nb = 301 · 513 = 154413 = NCWT .

4Ahnlich der Heisenbergschen Unscharferelation gibt es ein unteres Limit fur das Produkt der Stan-dardabweichungen eines Wavelets in Zeit- und Frequenzbasis σtσf ≥ 1/4π. [7]

5Auflosung der diskretisierten Signale, Uberfuhrung der Integrale in Summation etc.

9

t

Deh

nung

aMexican Hat

−256 −128 0 128 256

1

2

4

8

16

32

64

t

Deh

nung

a

Daubechies db6

−256 −128 0 128 256

1

2

4

8

16

32

64

Abbildung 2.4: Beispiele der CWT fur Einheitsimpuls δ(t)

2.5 Die diskrete Wavelet-Transformation

Wahrend der Skalierungskoeffizient a und der Verschiebungskoeffizient b bei der CWT kon-tinuierlich von 0 bis ∞ bzw. von −∞ bis ∞ variiert werden, so wird die diskrete Wavelet-Transformation (engl. discrete wavelet transform, DWT ) nur fur diskrete Wertepaare (a,b)evaluiert. Dazu wird der Zeitdehnungsparameter a in diskreten logarithmischen Schritten va-riiert und entsprechend mit dem Verschiebungsparameter b verknupft

ψm,n(t) =1√am0

ψ

(t− nb0am0

am0

), m ∈ N, n ∈ N. (2.12)

Somit konnen Zeitdehnung und -verschiebung mit den ganzzahligen Parametern m und ngesteuert und uber die Konstanten a0 > 1 und b0 > 0 skaliert werden. Ublicherweise werdendie Stutzstellen in einem dyadischen Gitter gewahlt, d. h. in ganzzahligen Zweierpotenzen.Das fuhrt zu orthonormalen Wavelet-Basisfunktionen, wodurch die DWT irredundant undverlustfrei umkehrbar wird. Deshalb ist eine ubliche Wahl der Parameter a0 = 2 und b0 = 1,was im Folgenden als Standardfall angenommen wird, sofern nicht anders angegeben. Einsetzender Werte in (2.12) fuhrt zur kompakteren Schreibweise der dyadischen Wavelet-Funktion

ψm,n(t) = 2−m/2 ψ(2−m t− n

). (2.13)

Die Wavelets sind genau dann orthonormal, also orthogonal und auf Energie 1 normiert, wenngilt ∫ ∞

−∞ψm,n(t)ψm′,n′(t) dt =

{1 m = m′, n = n′

0 sonst.(2.14)

Die Berechnung der DWT eines beliebigen Zeitsignals x(t) ist als

Tm,n =

∫ ∞−∞

x(t)ψm,n(t) dt (2.15)

definiert, was der Bestimmung der Wavelet- oder Detailkoeffizienten Tm,n der DWT mit denSkalierungs- und Zeitindizes m und n unter Verwendung des dyadischen Wavelet-Gitters (2.13)entspricht.Die Wahl einer orthonormalen Wavelet-Basis ψm,n(t) ermoglicht die Rekonstruktion des Zeit-signals aus den Wavelet-Koeffizienten Tm,n in Form einer Summation uber alle Indizes m undn gemaß

x(t) =∞∑

m=−∞

∞∑n=−∞

Tm,n ψm,n(t), (2.16)

10

was als inverse diskrete Wavelet-Transformation (IDWT ) bezeichnet wird. Außerdem gilt indiesem Fall Energiegleichheit,∫ ∞

−∞|x(t)|2 dt =

∞∑m=−∞

∞∑n=−∞

|Tm,n|2 , (2.17)

ahnlich dem Satz von Parseval fur die Fourier-Transformation. Aufgrund dieser Vorteile hatsich die dyadische orthonormale Basis der DWT in der Praxis etabliert und wird ublicherweiseals Standardfall vorausgesetzt.

2.6 Skalierungsfunktion und Multiskalenanalyse

Ahnlich zur schnellen Fourier-Transformation existiert auch eine schnelle (diskrete) Wavelet-Transformation (engl. fast wavelet transform, FWT ). Fur deren Herleitung werden die or-thonormalen dyadischen Wavelet-Funktionen ψm,n(t) mit sogenannten Skalierungsfunktionenverknupft. Eine Skalierungsfunktion beschreibt eine skalierte und verschobene Version einesZeitsignals φ(t) gemaß

φm,n(t) = 2−m/2 φ(2−m t− n

), (2.18)

φ0,0(t) = φ(t) wird auch als Vater-Wavelet bezeichnet. Anders als beim mittelwertfreienMutter-Wavelet gilt fur das Vater-Wavelet∫ ∞

−∞φ0,0(t) dt = 1. (2.19)

Eine Skalierungsfunktion ist orthogonal zu ihren Verschiebungen, aber nicht zu Dilatationen.Die Faltung der Skalierungsfunktion mit einem Signal x(t) liefert die sogenannten Approxima-tionskoeffizienten als eine gewichtete Mittelwertbildung

Sm,n =

∫ ∞−∞

x(t)φm,n(t) dt, (2.20)

Sm,n ist somit ein Maß dafur wie stark die Funktion φm,n(t) in x(t) vorhanden ist.Als nachstes wird mit xm(t) eine geglattete Version des Signals x(t) in Abhangigkeit desSkalierungsmaßstabs m beschrieben. Diese ergibt sich zu

xm(t) =∞∑

n=−∞

Sm,nφm,n(t). (2.21)

Fur sehr kleine Maßstabe (m → −∞) nahert sich die kontinuierliche Approximation xm(t)dem ursprunglichen Zeitsignal x(t). Bei einer geeigneten Wahl von Wavelet- und Skalierungs-funktion ψ(t), φ(t) lasst sich das Signal durch eine kombinierte Reihe der Approximationsko-effizienten und der DWT-Detailkoeffizienten

x(t) =∞∑

n=−∞

Sm0,nφm0,n(t) +

m0∑m=−∞

∞∑n=−∞

Tm,n ψm,n(t) (2.22)

fur einen beliebigen Skalierungsindex m0 reprasentieren. Fur die mathematische Herleitungdes erforderlichen Zusammenhangs zwischen φ(t) und ψ(t) sei hier auf Fachliteratur verwie-sen [7], [8], [9].Das Zeitsignal kann also durch seine Approximation fur den Maßstab m0 sowie die Summeder Details mit den Skalierungen m0 bis −∞ beschrieben werden. Als Bezeichung der konti-nuierlichen Signaldetails fur den Maßstab m sei

11

dm(t) =∞∑

n=−∞

Tm,n ψm,n(t), (2.23)

dann ergibt sich (2.22) durch Einsetzen von (2.21) und (2.23) zu

x(t) = xm0(t) +

m0∑m=−∞

dm(t). (2.24)

Da x−∞ → x(t) gilt, fuhrt Einsetzen von dm(t) = xm−1(t) − xm(t) zur sogenannten Multi-skalenanalyse (engl. multiresolution analysis, MRA oder multiscale approximation, MSA)

xm−1(t) = xm(t) + dm(t). (2.25)

Durch Addition von Approximation und Detail bei einem bestimmten Maßstab m erhalt mandie Approximation fur die nachsthohere Auflosung (m− 1) [9]. Dieser Zusammenhang ist dieGrundlage der schnellen Wavelet-Transformation.

2.7 Beispiele der Wavelet-Transformation

Zum besseren Verstandnis der kontinuierlichen und diskreten Wavelet-Transformation sind imFolgenden einige einfache Anwendungsbeispiele angefuhrt. Diese beschranken sich auf die qua-litative Aussage, die hier im Vordergrund steht. Auf Berechnungen und quantitative Angabenwird daher verzichtet. Das gilt auch fur die Visualisierung, fur beste Erkennbarkeit wurdendie Ergebnisse mit dem hochstmoglichen Kontrast zwischen Minimal- und Maximalwert dar-gestellt.

2.7.1 Einfaches Testsignal: Rechteck- und e-Funktion

In Kapitel 2.4 wurde bereits die Anwendung der CWT auf einen Einheitsimpuls illustriert(Abb. 2.4). Wie aus der Signaltheorie bekannt, bleibt die Waveletfunktion dadurch unverandert,d. h. das Ergebnis der CWT stellt das Wavelet in Abhangigkeit des Dehnungsmaßstabes dar.Analog dazu entspricht das Verhalten an einer Sprungfunktion dem mit der Sprunghohe ska-lierten Integral uber das Wavelet, wobei nicht vergessen werden darf, dass dieses mittelwertfreiist. Das bedeutet, dass der Korrelationswert wieder null wird, sobald sich die gesamte Wave-letfunktion mit einem neuen konstanten Signalverlauf uberlappt. Dies kann bei der Analysedes Rechtecksignals besonders gut beobachtet werden, wie in Abb. 2.5 dargestellt. Hier wur-de fur positives und negatives Rechteck jeweils eine Breite von 32 Samples gewahlt. DasMexican-Hat-Wavelet liefert hier z. B. bei einer Dehnung von 16 am starksten variierende Kor-relationswerte. Kleinere Maßstabe “verschwinden” in den konstanten Signalteilen, bei großerenMaßstaben ist die Anderung entsprechend geringer und uber eine langere Dauer “verschmiert”.Die charakteristischen Spitzen der CWT, die genauso bei den Flanken der e-Funktionen auf-treten, weisen darauf hin, dass bei allen analysierten “Frequenzen” ein Sprung mit der selbenPhase auftritt. Bei der DWT lasst sich beobachten, dass fur die Darstellung des relativ nie-derfrequenten Rechtecksignals nur Koeffizienten bis m = −4 erforderlich sind. Die feinerenBander enthalten hier also keine Energieanteile. Bei den Exponentialfunktionen kann man gutbeobachten, wie die DWT-Koeffizienten an den steilsten Stellen auch die hochfrequentestenAnteile aufweisen.

12

0 100 200 300 400 500 600 700 800 900 1000−1

−0.5

0

0.5

1

Testsignal

Am

plitu

de

CWT (Mexican Hat)

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128

CWT (Daubechies db6)

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128

DWT (Haar)

Samples

Ska

lieru

ngsi

ndex

m

100 200 300 400 500 600 700 800 900 1000

−9−8−7−6−5−4−3−2−1

0

min max

Abbildung 2.5: CWT und DWT von Rechteck- und e-Funktion

13

2.7.2 Analyse der Momentanfrequenz

Grundsatzlich eignen sich alle Wavelet-Typen aufgrund ihrer Skalierbarkeit zur Analyse vonFrequenzkomponenten in beliebigen Bandern. Allerdings wird die Qualitat der Korrelation vonder Wahl der Wavelets in Abhangigkeit der erwarteten Signalform beeinflusst. In diesem Beispielwird daher ein Sinussignal mit steigender Frequenz (“chirp”) mit dem ebenfalls sinusformigenkomplexen Morlet-Wavelet transformiert (Abb. 2.6). Die betragsmaßig starksten Auspragungender Wavelet-Transformierten kennzeichnen also die entsprechende “Momentanfrequenz” desTestsignals. Auch die DWT-Koeffizienten folgen dem steigenden Frequenzverlauf.

2.7.3 Ein EKG-Signal

Abbildung 2.7 zeigt die diskrete und kontinuierliche Wavelet-Transformation eines EKG-Signalsmit dem Haar- bzw. Mexican-Hat-Wavelet. Die DWT-Koeffizienten verdeutlichen, dass imBereich der QRS-Komplexe die starksten Signalanteile in mittel- bis hochfrequenten Bandernliegen. Die CWT spiegelt die Periodizitat des Signals wider und hebt die QRS-Spitzen deutlichhervor. Um die Nulldurchgange zu verdeutlichen, ist das CWT-Ergebnis außerdem mit loga-rithmischer Skalierung dargestellt. Dadurch heben sich positive bzw. negative Bereiche starkervom Nullwert ab.

14

0 100 200 300 400 500 600 700 800 900 1000−1

−0.5

0

0.5

1

TestsignalA

mpl

itude

CWT (komplexes Morlet−Wavelet, Realteil)

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128256512 min

max

CWT (komplexes Morlet−Wavelet, Betrag)

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128256512 max

min

DWT (Haar)

Samples

Ska

lieru

ngsi

ndex

m

100 200 300 400 500 600 700 800 900 1000

−9−8−7−6−5−4−3−2−1

0 min

max

Abbildung 2.6: CWT und DWT eines Chirp-Signals

15

0 100 200 300 400 500 600 700 800 900 1000

0

1

2

Testsignal (qtdb/sel103)

Am

plitu

de

DWT (Haar)

Ska

lieru

ngsi

ndex

m

100 200 300 400 500 600 700 800 900 1000

−9−8−7−6−5−4−3−2−1

0

CWT (Mexican Hat)

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128

CWT logarithmisch

Samples

Deh

nung

a

100 200 300 400 500 600 700 800 900 1000

1248

163264

128

min max

Abbildung 2.7: DWT und CWT eines EKG-Signals im Vergleich

16

3 Implementierung der EKG-Zeitbereichsanalyse

Um ein EKG-Signal detailliert zu analysieren, wurden zunachst Zeitbereichsmethoden gewahlt,da diese intuitiv am einfachsten erscheinen. Als langfristigere Zielsetzung sollen im Zuge wei-terer Arbeiten andere Methoden hinsichtlich Zuverlassigkeit und Genauigkeit mit den hierangefuhrten verglichen werden.Der gesamte Algorithmus besteht aus mehreren Arbeitsschritten, die in den folgenden Kapitelnerortert werden.

1. Detektion der QRS-Komplexe

2. Abschatzung der Herzfrequenz

3. Preprocessing

4. Lokalisierung der Einzelwellen (Kovarianzfilter)

5. Approximation durch Gaußkurven-Segmente

6. Vergleich mit Referenz-Markern

7. Statistische Auswertung

3.1 Detektion der QRS-Komplexe

Als solide Grundlage der EKG-Analyse sollte eine bewahrte Methode der QRS-Erkennungherangezogen werden. Die Wahl fiel auf einen Algorithmus von Jiapu Pan und Willis J. Tomp-kins [10], der bereits als fertige Matlab-Funktion zum Download zur Verfugung steht [15].Abb. 3.1 zeigt ein Beispiel fur die QRS-Erkennung mit der gewahlten Methode.

0 5 10 15−0.5

0

0.5

1

1.5

2

Zeit/s

Am

plitu

de/m

V

Pan−Tompkins−QRS−DetektionR detektiert

Abbildung 3.1: Ausgabe des Pan-Tompkins-Algorithmus (MATLAB-Code [15])

3.1.1 Algorithmus

Die Grundidee ist, den relativ amplitudenstarken und hochfrequenten QRS-Komplex durch Dif-ferenzieren hervorzuheben. Da positive und negative Flanken auftreten, wird das Ergebnis nochquadriert, womit automatisch große Amplituden nochmals eine hohere Verstarkung erfahren.Zusatzlich kommen zur Qualitatsverbesserung verschiedene Filterfunktionen zum Einsatz.

17

Vollstandiger Algorithmus im Detail:

1. Bandpassfilterung - Das Filter mit einem Durchlassbereich von 5 Hz bis 11 Hz (−3 dB)6

unterdruckt Storungen wie Rauschen, Grundliniendrift und Netzbrumm.

2. Differentiation - Ein digitales Differentiationsfilter (naherungsweise ideal von 0 Hz bis30 Hz) liefert die Flanken des QRS-Komplexes.

3. Quadrierung - Durch punktweises Quadrieren werden alle Datenpunkte positiv gemachtund nichtlinear verstarkt.

4. Integration - Die moving-average-Integration (uber ein bewegtes Zeitfenster) fasst Q-,R- und S-Wellen zusammen, wobei die Fensterbreite fur optimale Ergebnisse ungefahrgleich der erwarteten Dauer eines QRS-Komplex gewahlt wird (150 ms).

5. QRS-Erkennung - Die steigende Flanke des Integrals entspricht der Position eines QRS-Komplexes.

6. Plausibilitatsprufung - Anhand verschiedener Parameter wie z. B. Signal- bzw. Rauscham-plitude sowie zuvor detektierter (regelmaßiger) R-R-Intervalle wird das Ergebnis bewertetund ggf. nachkorrigiert7.

3.1.2 Signale der Einzelstufen

Die Signalformen der jeweiligen Filterstufen sind in Abb. 3.2 dargestellt. Dazu wurde einsynthetisches Testsignal verwendet (Generator von McSharry & Clifford [11], Ausschnitt in derLange von etwas mehr als einem Herzzyklus, Parametersatz siehe Tab. 3.1).

Parametersatz fur synthetisches EKG (rauschfrei)Name Wert Dim. Bedeutung

sfecg 512 Hz Sampling-FrequenzN 16 - Anzahl von HerzzyklenAnoise 0 mV Additives gleichverteiltes Rauschenhrmean 60 bpm Durchschnittliche Herzfrequenzhrstd 0 bpm Standardabweichung der Herzfrequenzlfhfratio 0.5 - Verhaltnis niedrigster zu hochster Frequenzsfint 256 Hz Interne Sampling-Frequenzti - rad Extremwert-Winkel (default)ai - - Extremwert-Amplituden (default)bi - - Zackenbreiten (default)

Tabelle 3.1: Parametersatz fur ecgsyn.m

6Optimal ware ein Durchlassbereich von 5 Hz bis 15 Hz [12], hier wurde ein Kompromiss zugunsteneffizienter Implementierbarkeit auf einem Mikroprozessor eingegangen [10].

7Wenn es z. B. scheint, als ware ein QRS-Komplex “ubersehen” worden, wird die QRS-Erkennungim entsprechenden Bereich mit einem niedrigeren Schwellwert erneut durchgefuhrt. Ebenso werden mut-maßliche falsche Positiverkennungen durch P-Wellen verworfen, wenn sie innerhalb von 360 ms nach derletzten gultigen QRS-Erkennung auftreten und ihre Integralflanken weniger als halb so steil sind.

18

0 100 200 300 400 500 600−0.5

0

0.5

1

Synthetisches EKG

0 100 200 300 400 500 600

−0.5

0

0.5

1

Bandpass−gefiltert

0 100 200 300 400 500 600−1

0

1

differenziert

0 100 200 300 400 500 6000

0.5

1

quadriert

0 100 200 300 400 500 6000

0.2

0.4

integriert

0 100 200 300 400 500 600−0.5

0

0.5

1

1.5

Samples

OriginalsignalQRS detektiert

Abbildung 3.2: Einzelstufen des Pan-Tompkins-Algorithmus

19

3.2 Abschatzung der Herzfrequenz

Fur viele der angewendeten Zeitbereichsfilter ist es notwendig, die ungefahre Periodendauerder Herzzyklen zu kennen. Da eine Naherung ausreichend ist, wird die Durchschnittszeit durchdas Verhaltnis der gesamten Signaldauer zur Anzahl detektierter QRS-Komplexe abgeschatzt(Glg. 3.1). Der Kehrwert liefert die Pulsfrequenz, welche ublicherweise in Schlagen pro Mi-nute (bpm) angegeben wird (Glg. 3.2). Diese kann in einem sehr großen Bereich variieren(ca. 30 bis uber 200 bpm), welcher nicht bzw. nur mit erheblichen Einbußen von einem ein-zigen “Universalfilter” abgedeckt werden kann. Es wurden zwar auch Versuche mit parallelenFilterbanken fur verschiedene Frequenzbereiche unternommen, diese Methode wurde jedochwieder verworfen, da der Gewinn an Stabilitat und Signalqualitat nicht zufriedenstellend aus-fiel. Die grobe Anpassung der Filter auf die momentane Herzfrequenz hingegen ist intuitiv,einfach zu implementieren und lieferte sehr stabile Resultate.

Tav ≈Tges

NQRS − 1=tend − tstartNQRS − 1

(3.1)

HR ≈ 1

Tav

60 s

1 min[HR] = min−1 (3.2)

Naturlich mussen mindestens 2 Herzschlage detektiert werden, um eine Schatzung der Fre-quenz zu ermoglichen. Sollte dies nicht der Fall sein, terminiert das Programm an dieser Stellemit einer Fehlermeldung.

3.3 Preprocessing

3.3.1 QRS-Unterdruckung

Nach der Erkennung der QRS-Komplexe sollen in weiterer Folge P- und T-Welle genau lokali-siert werden. Zur Verbesserung von Genauigkeit und Zuverlassigkeit der entsprechenden Filterwerden dazu die Signalsegmente in einem definierten Zeitintervall rund um jede detektierteR-Zacke unterdruckt (Abb. 3.3). Die Interpolationsfunktion ahnelt einem Zero-Order-Hold-Glied, nur dass anstelle des Funktionswertes zum Startzeitpunkt der Mittelwert von beidenIntervallgrenzen gewahlt wurde. Gleichung 3.3 beschreibt die Berechnungsvorschrift fur dasQRS-unterdruckte Signal s aus dem Originalsignal s, wobei ∆t die geschatzte halbe Dauer ei-nes QRS-Komplexes zum Zeitpunkt TR definiert. Als Standardwert wurde ∆t = 75 ms aus [10]ubernommen.

s(t) =

s (TR −∆t) + s (TR + ∆t)

2TR −∆t < t < TR + ∆t

s sonst(3.3)

3.3.2 Weitere Filter

Durch die beschriebene Unterdruckung der QRS-Komplexe konnen gegebenenfalls storendeUnstetigkeiten im Signal entstehen. Um diese zu eliminieren, wurde ein Tiefpassfilter vorgese-hen, wobei die Wahl auf ein finite impulse response (FIR)-Filter fiel, um Verfalschungen des Si-gnals durch Laufzeitverzerrungen zu vermeiden. Auch dem Originalsignal kann moglicherweisehochfrequentes Rauschen uberlagert sein, welches ebenfalls mit dem FIR-Tiefpass unterdrucktwird. Praxistests8 zeigten jedoch, dass der Einfluss hochfrequenter Storungen außerst geringist, weshalb die Zusatzfilter standardmaßig deaktiviert werden.

8Verschiedene mehr oder weniger stark verrauschte synthetische Signale wurden mit und ohne Tiefpass-filter analysiert. Dabei wurde eine merkliche Qualitatsverbesserung durch den Tiefpass erst fur ubermaßig

20

3.5 4 4.5 5 5.5 6 6.5 7 7.5 8−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

Synthetisches EKG

3.5 4 4.5 5 5.5 6 6.5 7 7.5 8−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Zeit/s

Am

plitu

de/m

V

QRS unterdrückt

Abbildung 3.3: Preprocessing: Unterdruckung der QRS-Komplexe

21

3.4 Lokalisierung der Einzelwellen (Kovarianzfilter)

Die Position der Maxima der einzelnen EKG-Teilwellen soll moglichst genau bestimmt wer-den. Da lokale Maxima durch Storungen wie Rauschen oder Nulliniendrift verfalscht werdenkonnen, kommt hierfur ein Kovarianzfilter zum Einsatz. Aufgrund der Ahnlichkeit mit den EKG-Signalteilen wurde als Filter die Gauß’sche Normalverteilungsfunktion (Abb. 3.5) gewahlt, diedurch ihr Bandpassverhalten bereits einen Großteil unerwunschter Signalanteile unterdruckt.Naturgemaß sind die Ergebnisse der Kovarianzfunktion genau dann am aussagekraftigsten,wenn als Filterfunktion die gesuchte Signalform verwendet wird (vgl. “Matched Filter”). Daherist es naheliegend, bei den Filterfunktionen fur die Einzelwellen jeweils verschiedene a prioridefinierte Standardabweichungen zu verwenden. Die verwendeten Werte wurden durch em-pirische Tests mit synthetischen EKG-Signalen ermittelt und sind in Tabelle 3.2 angefuhrt,wobei 100 % der durchschnittlichen Dauer eines Herzzyklus entspricht. Abbildung 3.4 zeigt alsBeispiel die Ergebnisse der Filter fur P- und S-Wellen.Vorerst wurde keine Plausibilitatsuberprufung der detektierten Maxima vorgesehen, um auszu-schließen, dass ursprunglich korrekt erkannte Wellen als mutmaßliche Fehlerkennungen wiederverworfen werden konnten. Somit wird vor und nach jedem QRS-Komplex eine P- und ei-ne T-Welle detektiert, unabhangig von der Qualitat der Resultate. An dieser Stelle sollte mitmedizinischem Fachwissen angesetzt werden, um moglichst zuverlassige Entscheidungskriterienzu finden.Die Kreuzkovarianz is laut MATLAB-Dokumentation [1] definiert als

cxy(m) =

M−|m|−1∑

n=0

(x(n+m)− 1

N

N−1∑i=0

xi

)(y∗n −

1

N

N−1∑i=0

y∗i

)m ≥ 0

c∗yx(−m) m < 0.

(3.4)

Welle σFilter

P 4,0 %Q 0,5 %R 1,0 %S 0,5 %T 3,5 %

Tabelle 3.2: Relative Standardabweichungen der Gauss-Kovarianzfilter

3.5 Approximation durch Gaußkurven-Segmente

Zur Charakterisierung der EKG-Wellen sind neben Position und Amplitude der Maxima auchBreite, Zeitdifferenz sowie Anfangs- und Endpunkte der einzelnen Zacken von Bedeutung.All diese Informationen konnen durch eine Approximation der Teilwellen durch Segmente vonGauß’schen Glockenkurven gewonnen werden. Mit einer Least-Squares-Naherung erhalt manden Mittelwert µ und die Standardabweichung σ, womit die Position des Maximums sowiedie “Breite” der Welle bestimmt werden. Außerdem konnen damit Start- und Endzeitpunktals µ± 1, 8σ definiert werden. Dieser Skalierungsfaktor geht aus empirischen Vergleichen mitvon Experten manuell gesetzten Markern hervor [13], da es aufgrund von Nulliniendrift etc.

verrauschte Testsignale beobachtet, bei denen die eigentliche Signalform unter dem Rauschen nicht mehrerkennbar war. Dieser Extremfall wurde aber als praktisch irrelevant eingestuft.

22

5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

Synthetisches EKGP ReferenzQ ReferenzR ReferenzS ReferenzT Referenz

5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5−0.2

−0.1

0

0.1

0.2

0.3

Zeit/s

Am

plitu

de/m

V

Kovarianz PP detektiert

5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

Zeit/s

Am

plitu

de/m

V

Kovarianz SS detektiert

Abbildung 3.4: Ergebnis der Kovarianzfilterung fur P- und S-Wellen

23

0 10 20 30 40 50 60 70 80 900

0.01

0.02

0.03

0.04

0.05

Samples n

Impu

lse

Res

pons

e h P

[n]

Abbildung 3.5: Gaußfilter fur P-Wellen (fS = 256 Hz, σ = 0,04 s, T = 0,35 s)

praktisch nicht moglich ist, eindeutige Anfangs- bzw. Endpunkte der Wellen zu definieren.Abb. 3.6 zeigt einen Ausschnitt eines synthetischen EKG-Signals, dessen Teilwellen durchGaußkurven approximiert wurden. Die hervorgehobenen Segmente decken jeweils den Bereichvon µ− 1, 8σ ≤ t ≤ µ+ 1, 8σ ab. Das Approximationsproblem

s(t) ∼= A · 1

σ√

2π· e−

12( t−µσ )

2

+ d

∣∣∣∣tstart≤t≤tend

(3.5)

fuhrt zu Amplitudenskalierungsfaktor A, Standardabweichung σ, Mittelwert µ sowie Nullinien-Offset d als Schatzwerte der Least-Squares-Optimierung. Als Zeitintervall fur die Optimierungwird rund um die zuvor bestimmten Maxima ein Bereich anhand vordefinierter Werte fur die zuerwartenden Standardabweichungen gewahlt. Diese wurden durch empirische Tests ermitteltund sind in Tab. 3.3 angefuhrt, die Prozentwerte beziehen sich auf die durchschnittliche Dauereines Herzzyklus. Es zeigte sich ebenfalls, dass die Approximation umso schlechtere Ergebnisseliefert, je starker die Teilwellen asymmetrisch um den Maximalwert sind. Da aber aus medizini-scher Sicht zur Charakterisierung von Krankheitsbildern oft die moglichst exakte Bestimmungvon Anfangs- oder Endpunkt der Wellen relevant ist, z. B. QT-Intervall [13], kann in diesemFall die Approximation nur auf der steigenden oder fallenden Seite angewendet werden. Dadie Ergebnisse stark vom gewahlten Intervall abhangen, wurden Mittelpunkt und 1,8-facheStandardabweichung der vorherigen Naherungslosung fur die Intervallgrenzen herangezogen.

Welle σ

P 5,0 %Q 0,5 %R 1,2 %S 0,4 %T 4,0 %

Tabelle 3.3: Standardabweichungen fur Gaußkurven-Approximation

24

11.7 11.8 11.9 12 12.1 12.2 12.3 12.4−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

Abbildung 3.6: Reprasentation der Teilwellen durch Gaußkurven-SegmenteHier ist zu erkennen, dass aufgrund schlechter Symmetrie von Q- und S-Wellen keine

gute Approximation durch die symmetrischen Gaußkurven moglich ist (dunne Linien).Abhilfe schafft das getrennte Approximieren von linkem und rechtem Ast; somit sind

wesentlich prazisere Naherungen moglich (dicke Linien).

3.6 Vergleich mit Referenz-Markern

Ein wichtiges Kriterium fur die automatisierte EKG-Datenverarbeitung ist die Genauigkeit undZuverlassigkeit bei der Erkennung der jeweiligen Extremwerte9. Um diese beurteilen zu konnen,werden die Resultate des Algorithmus mit einer bekannten Referenz verglichen. Diese soge-nannten Marker werden bei synthetischen EKG-Signalen automatisch generiert; bei “echten”EKG-Messungen wurden die Marker bei einigen Datensatzen manuell von einer Gruppe vonExperten gesetzt (z.B. PhysioNet-Datenbank [16]).Bei der implementierten Losung werden die Referenzdaten fur P- bis T-Wellen mit den de-tektierten Daten verglichen. Dazu werden zuerst die nachsten Nachbarn ermittelt, welcheeinander anschließend eindeutig zugeordnet werden, d. h. nur die besten Ubereinstimmungenwerden als gultig angesehen. Beiderseits werden Marker ohne Ubereinstimmung hervorgeho-ben, um Fehlerkennungen zu verdeutlichen. An dieser Stelle konnten zusatzliche Plausibi-litats- bzw. Qualitatsbedingungen hinzugefugt werden, z. B. limitierter Zeitversatz fur gultigeUbereinstimmungen, Verwerfen von Ausreißern mit geringen Amplituden. Um die Bewertungder gewonnenen Resultate nicht zu verfalschen, wurde darauf zunachst verzichtet.Der Vergleich der Marker ist in Abb. 3.7 in einem Beispiel dargestellt, bei dem ein unerkann-ter QRS-Komplex auftritt11. Da fur die Analyse stets der Bereich zwischen zwei R-Zackenherangezogen wird, gibt es somit fur die ersten funf Wellen keine Ubereinstimmungen, wasdurch die fett gedruckten Marker hervorgehoben wird. Auf diese Weise wurden auch falschepositive Erkennungen markiert werden, die allerdings in diesem Beispiel nicht erscheinen. Un-ter der Voraussetzung, dass die erste und letzte R-Zacke des Signals nicht detektierbar sind,alle anderen hingegen korrekt detektiert werden, hinge die Erkennungsrate ausschließlich vonder gesamten Signallange ab. Aus diesem Grund wurden unvollstandige Herzzyklen an denSignalrandern von der statistischen Auswertung ausgenommen.

9Fur manche Anwendungen sind außerdem Anfangs- und Endpunkt der Teilwellen10 relevant.11Da sich der erste QRS-Komplex exakt am Anfang des Datensatzes befindet, kann er mit der gewahlten

Methode nicht mit ausreichender Sicherheit detektiert werden.

25

0 0.5 1 1.5 2 2.5 3 3.5−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

ReferenzmarkerP ReferenzQ ReferenzR ReferenzS ReferenzT ReferenzP FehlerQ FehlerR FehlerS FehlerT Fehler

0 0.5 1 1.5 2 2.5 3 3.5−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

Detektierte MarkerP detektiertQ detektiertR detektiertS detektiertT detektiert

Abbildung 3.7: Referenz- und detektierte Marker im VergleichDer unvollstandige QRS-Komplex am Rand des Signals wird nicht erkannt, weshalb es

fur die Referenzmarker des ersten Herzzyklus keine Ubereinstimmungen mit detektiertenMarkern gibt.

26

3.7 Statistische Auswertung

Um die Qualitat der EKG-Analyse beurteilen zu konnen, werden ublicherweise die Gute unddie Zuverlassigkeit der Wellenerkennung verglichen. Dazu werden fur die jeweiligen Wellen dieAbweichungen von allen erkannten Markern an den Zeitpunkten Ti zu den Referenzzeiten Tiberechnet. Aus diesen werden der Mittelwert

µ =1

N

N∑i=1

(Ti − Ti

), (3.6)

und die Standardabweichung

σ =

√√√√ 1

N − 1

N∑i=1

(Ti − Ti

)2(3.7)

bestimmt, wobei N fur die Gesamtzahl korrekt detektierter Wellen steht. Die Abweichung deskleinsten Schatzwertes

∆min = min∀i

(Ti − Ti

)(3.8)

und die des großten

∆max = max∀i

(Ti − Ti

)(3.9)

werden ebenfalls ermittelt, um Symmetrie und systematische Fehler beurteilen zu konnen. DieGaußkurven-Approximationen werden mit der gleichen Methode bewertet, was einen direktenQualitatsvergleich ermoglicht. Außerdem wird die jeweilige Fehleranzahl fur nicht detektier-te oder falsch positiv detektierte Zacken anhand der Ubereinstimmungen mit den Referenz-Markern ausgegeben.Kapitel 5 enthalt einige Beispiele fur die Ausgabe der statistischen Daten.

27

4 Grafische Benutzeroberflache (GUI)

4.1 Idee und Funktion

Mit dem Matlab-Tool “guide” wurde eine Grafische Benutzeroberflache (engl. Graphical UserInterface, kurz GUI) zum Auswerten der EKG-Signale erstellt, da diese zahlreiche Filter- undVerarbeitungsstufen durchlaufen. Um die Auswirkungen der einzelnen Prozesse analysierenzu konnen, sind zahlreiche Visualisierungen der Zwischenschritte erforderlich, was besondersbei großen Datenmengen schnell zu außerst unubersichtlichen Resultaten fuhrt. Der erstellte“Waveform-Viewer” steigert die Benutzerfreundlichkeit in dieser Hinsicht mit der Funktion,zwei beliebige Signale direkt miteinander vergleichen zu konnen. Dazu stehen eine komfortableZoom-Funktion sowie weitere Visualisierungsoptionen zur Verfugung.

4.2 Umsetzung

Das gesamte GUI ist in Abb. 4.1 dargestellt. In den Plotfenstern (oben) werden zwei Si-gnale dargestellt, deren Quellen aus den jeweiligen Listen (rechts) unabhangig voneinanderausgewahlt werden konnen. Sofern der Datensatz Markerinformationen enthalt, konnen dieseoptional den Diagrammen hinzugefugt werden. Der untere Plot dient als Ubersicht fur das Ge-samtsignal, wobei der darunter gelegene Schieberegler gleichzeitig den Zoombereich anzeigt.Signalquelle des Ubersichtsfensters und Zoomstufe werden uber die Steuerelemente (rechtsunten) eingestellt. Die effiziente Dokumentation der Ergebnisse erfolgt uber die Schaltflache“Save Screenshot”, wodurch eine .png-Datei mit Datum und Uhrzeit der Erstellung als Da-teiname im Unterordner “fig/” gespeichert wird.

28

4.5 5 5.5 6 6.5 7 7.5 8−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

4.5 5 5.5 6 6.5 7 7.5 8−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

0 5 10 15−0.5

00.5

11.5

Synthetisches EKGPQRST

Detektierte MarkerPQRST

Abbildung 4.1: GUI: Waveform-Viewer

29

5 Performance-Bewertung

Ursprunglich wurde eine Auswertung anhand echter EKG-Daten mit von Experten gesetztenMarkern angestrebt. Dazu bot Fachpersonal des AKH Linz freundlicherweise seine Kooperati-on an. Aufgrund von Schwierigkeiten beim Generieren der extrem umfangreichen Datensatzewurde letztendlich doch entschieden, das Programm mit synthetischen Signalen zu testen. Lei-der liefert der EKG-Generator von McSharry & Clifford [11] als Referenz nur die Zentren derEinzelwellen, weshalb die Analyse der Wellenbreite bzw. Start- und Endpunkte hiermit nichtbewertet werden kann.

5.1 Synthetischer Datensatz (rauschfrei)

Zunachst soll der Optimalfall bewertet werden, also ein rauschfreies Signal mit einer geringenVarianz der Herzfrequenz. Ein Zeitfenster von 5 Minuten deckt dabei ca. 300 Herzzyklen mit 60bpm ab, was aufgrund der nahezu idealen Periodizitat als ausreichend angesehen werden kann.Die Parameter fur den EKG-Generator und die Ausgabe der Statistikdaten sind in Tab. 5.1und Tab. 5.2 angefuhrt.

Parametersatz fur synthetisches EKG (rauschfrei)Name Wert Dim. Bedeutung

sfecg 256 Hz Sampling-FrequenzN 300 - Anzahl von HerzzyklenAnoise 0 mV Additives gleichverteiltes Rauschenhrmean 60 bpm Durchschnittliche Herzfrequenzhrstd 1 bpm Standardabweichung der Herzfrequenzlfhfratio 0.5 - Verhaltnis niedrigster zu hochster Frequenzsfint 256 Hz Interne Sampling-Frequenzti - rad Extremwert-Winkel (default)ai - - Extremwert-Amplituden (default)bi - - Zackenbreiten (default)

Tabelle 5.1: Parametersatz 1 fur ecgsyn.m

5.1.1 Herzfrequenz

Die geschatzte Herzfrequenz HR = 60, 0 bpm stimmt erwartungsgemaß exakt mit der Vor-gabe uberein. Das Ergebnis uberrascht nicht, da die Frequenz des Generators naherungsweisekonstant ist und die Schatzung auf der Anzahl der detektierten QRS-Komplexe bezogen aufdie Gesamtdauer basiert (Glg. 3.2). Sofern also alle QRS-Komplexe korrekt erkannt werden,stimmt der Durchschnittswert der Herzfrequenz sehr genau.

5.1.2 Detektionsrate

Die Detektionsrate ist ein wichtiges Maß fur die Zuverlassigkeit des Algorithmus. Auffalligist, dass bei diesem Datensatz von jedem Wellentyp jeweils zwei Vorkommen nicht detektiertwurden. Die Ursache liegt darin, dass der erste und letzte QRS-Komplex des Signals unvoll-standig sind, und daher nicht als gultig erkannt werden (Abb. 3.7). Die restlichen Teilwellenwerden jeweils zwischen zwei aufeinanderfolgenden R-Zacken bestimmt, weshalb der Bereichvor der ersten bzw. nach der letzten QRS-Erkennung keine weiteren positiven Resultate liefern

30

Geschatzte Herzfrequenz

HR = 60, 0 bpm

Marker-Schatzung durch Kovarianzfilterµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P -0.2981 1.4266 -3.9063 3.9063 2 0Q 7.2392 1.4051 3.9063 11.7188 2 0R -0.5341 1.3433 -3.9063 0 2 0S -5.6415 1.9428 -7.8125 -3.9063 2 0T 0.8638 1.6227 0 3.9063 2 0

Marker-Schatzung durch Gaußkurven-Approximationµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P -0.2832 1.1748 -2.8715 2.1229 2 0Q -0.1513 1.5641 -2.9350 11.7188 2 0R -0.3082 1.1622 -2.4195 1.7743 2 0S 0.8030 1.3425 -1.2968 2.7245 2 0T -0.4091 1.2071 -3.7083 1.9352 2 0

Tabelle 5.2: Analyse des rauschfreien synthetischen EKGMan beachte, dass die Abweichungen ∆min und ∆max die “negativste” und “positivste”Differenz zwischen Schatzwert und wahrem Wert beschreiben (3.8), (3.9), und nicht etwaden betragsmaßig kleinsten/großten Fehler. Dies ist von Noten, da ∆min und ∆max dasselbe Vorzeichen aufweisen konnen, also alle Abweichungen positiv bzw. negativ sind.

kann. Diese Rander zu berucksichtigen ware also sinnlos, da in diesem Fall die Detektionsratelediglich von der willkurlichen Wahl der Signallange abhinge. Beschrankt man die Bewertungauf das Signal zwischen erster und letzter detektierbarer R-Zacke, liegt die Erkennungsrate furalle Teilwellen bei 100%.Falsche Positiverkennungen treten nicht auf. Im rauschfreien Testsignal treten keine Artefakteauf, die irrtumlich als QRS-Komplexe identifiziert werden.

5.1.3 Standardabweichung

Um die Qualitat der Schatzung der Zeitpunkte der Wellenzentren zu beurteilen, sollte dieermittelte Standardabweichung auf die Sampling-Frequenz bezogen werden. Fur die Abtast-zeit ergibt sich TS = 1/fS = 1/256Hz ≈ 4 ms. Die Standardabweichungen liegen also in derGroßenordnung von 1/3 bis 1/3 der Abtastzeit, was bedeutet, dass im Mittel uber 95% derdetektierten Marker genauer als ±1 Sample mit den Referenzen ubereinstimmen. Durch Ap-proximation mit Gaußkurvensegmenten kann das Ergebnis sogar noch geringfugig verbessertwerden.

5.1.4 Mittelwerte

Die Mittelwerte der geschatzten Marker weichen fur P-, R- und T-Wellen nur einen Bruchteilder Abtastzeit ab, was hinsichtlich der Genauigkeit bereits mehr als zufriedenstellend ist. Beiden unsymmetrischen Q- und S-Zacken liefert die einfache Kovarianzfilterung mit Gaußkurvenhingegen merklich schlechtere Resultate. Auch diese konnen aber durch die Kurvenapproxima-tionen erheblich verbessert werden.

31

5.2 Synthetischer Datensatz (stark verrauscht)

Im Folgenden wird ein synthetisches EKG ausgewertet, dem gleichverteiltes Rauschen (engl.uniform noise) uberlagert ist. Um den Vergleich mit dem rauschfreien Datensatz nicht zuverfalschen, wurden alle anderen Parameter unverandert belassen. Das Rauschmodell ist wahr-scheinlich keine realistische Approximation der zu erwartenden Storungen bei der EKG-Messung.Da es aber als Parameter beim McSharry-Clifford-Generator definiert werden kann, ermoglichtes zumindest reproduzierbare Performancetests der Programmroutinen. Tabelle 5.3 zeigt denvollstandigen Parametersatz des EKG-Generators. Das Resultat der Analyse ist in Tabelle 5.4angefuhrt.

Parametersatz fur synthetisches EKG (verrauscht)Name Wert Dim. Bedeutung

sfecg 256 Hz Sampling-FrequenzN 300 - Anzahl von HerzzyklenAnoise 0.25 mV Additives gleichverteiltes Rauschenhrmean 60 bpm Durchschnittliche Herzfrequenzhrstd 1 bpm Standardabweichung der Herzfrequenzlfhfratio 0.5 - Verhaltnis niedrigster zu hochster Frequenzsfint 256 Hz Interne Sampling-Frequenzti - rad Extremwert-Winkel (default)ai - - Extremwert-Amplituden (default)bi - - Zackenbreiten (default)

Tabelle 5.3: Parametersatz 2 fur ecgsyn.m

Geschatzte Herzfrequenz

HR = 60, 0 bpm

Marker-Schatzung durch Kovarianzfilterµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P -0.4587 8.7054 -27.3438 42.9688 2 0Q 6.2225 5.0941 -11.7188 19.5313 2 0R -0.3891 1.8117 -3.9063 3.9063 2 0S -4.9459 4.5627 -15.6250 7.8125 2 0T 1.0091 11.7467 -35.1563 35.1563 2 0

Marker-Schatzung durch Gaußkurven-Approximationµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P -0.6140 7.9593 -24.6264 29.1798 2 0Q 4.9017 5.6156 -14.2210 17.7959 2 0R 0.0345 2.1586 -5.5295 6.8570 2 0S -3.4853 5.1654 -16.9644 21.2372 2 0T -0.2881 6.7430 -22.4528 23.4375 2 0

Tabelle 5.4: Analyse des verrauschten synthetischen EKG

32

5.2.1 Herzfrequenz

Trotz des Rauschens werden alle QRS-Komplexe des Testsignals korrekt detektiert, womit diegeschatzte Herzfrequenz HR = 60, 0 bpm wiederum exakt stimmt.

5.2.2 Detektionsrate

Wie im rauschfreien Fall bleiben auch hier lediglich Zacken an den Signalrandern unerkannt.Vernachlassigt man diese, liegt die Detektionsrate fur alle Teilwellen bei 100%. Es gibt keinVorkommen falscher Positiverkennungen.

5.2.3 Standardabweichung

Durch das Hinzufugen von Rauschen steigt die Standardabweichung der geschatzten Zeitpunk-te erwartungsgemaß erheblich. Bei diesem Datensatz ergibt sich der großte Wert zu ST ≈ 8 ms.Bezogen auf die Abtastzeit TS ≈ 4 ms betragt die schlechteste Standardabweichung also ca.2 Samples. In Anbetracht der schwierig zu identifizierenden Wellenberge (siehe Abb. 5.1) istdieses Resultat immer noch als zufriedenstellend zu werten. Vergleichswerte aus [14] liegen furverschiedene aktuelle Methoden z. B. fur σR im Bereich von 5 bis 20 ms, der hier erzielte Wertliegt bei σR ≈ 2,2 ms.Die Approximation der Teilwellen durch Gaußkurven bringt fur die schmalen Q-, R- und S-Zacken keinen Gewinn an Genauigkeit, auch fur die P-Wellen ist die Verbesserung nicht nen-nenswert. Die Standardabweichung der breiteren T-Wellen kann hingegen fast halbiert werden.

129 129.5 130 130.5 131 131.5−1

−0.5

0

0.5

1

1.5

Zeit/s

Am

plitu

de/m

V

Abbildung 5.1: Wellenerkennung des verrauschten Testsignals

5.2.4 Mittelwerte

Die Schatzung der Mittelwerte wird durch das uberlagerte Rauschen nicht nennenswert ver-schlechtert. Auf den umfangreichen mathematischen Beweis dieser Beobachtung wird verzich-tet, da die Erwartungstreue der Schatzer nicht im Mittelpunkt der Arbeit steht - die intuitiveErklarung ist durch die Eigenschaften des Rauschsignals gegeben: Aufgrund der gleichverteil-ten Amplitude ist der Erwartungswert des Rauschanteils E(Anoise) = 0. Die Storungen imBereich der Wellenberge konnen also dazu fuhren, dass diese zu fruh oder zu spat lokalisiertwerden, wodurch deren Standardabweichung zunimmt - dieser Fehler verschwindet aber imstatistischen Mittel.

33

5.3 Gegenuberstellung

Zusatzlich zu den in Kapitel 5.1 und 5.2 ausfuhrlich beschriebenen Ergebnissen fur ein ganzlichrauschfreies sowie ein stark verrauschtes Signal wurde im Sinne statistischer Signifikanz ei-ne umfangreiche Testreihe durchgefuhrt, die in Tabelle 5.5 aufgelistet ist. Dabei wurden dieParameter der durchschnittlichen Herzfrequenz, deren Standardabweichung sowie der Ampli-tude des uberlagerten Rauschens willkurlich variiert, wobei zunachst eine Abtastfrequenz vonfS = 256 Hz sowie eine ungefahre Signallange von N ≈ 300 Herzzyklen als globale Parame-ter festgelegt wurden. Eine zweite Testreihe mit den globalen Parametern fS = 512 Hz undN ≈ 100 zeigt, dass diese (in einer sinnvollen Großenordnung gewahlt) nahezu keinen Einflussauf die Ergebnisse haben.Vergleicht man die Schatzwerte der P- bis T-Wellen, fallt als erstes die insgesamt ausgezeich-nete Qualitat der Mittelwerte auf. Die mittleren Abweichungen der geschatzten Zeitpunktebetragen in den meisten Fallen deutlich weniger als 1 ms und liegen damit klar unter derAbtastzeit von TS1 ≈ 4 ms bzw. TS2 ≈ 2 ms. Lediglich bei den Q-Wellen zeichnet sich einetwa doppelt so großer Fehler ab, der damit aber noch immer in einer sehr zufriedenstel-lenden Großenordnung liegt. Außerdem ist deutlich erkennbar, dass alle Abweichungen einensystematischen Fehler12 aufweisen. Dieser ist durch die Art der Approximation zu erklaren,da die gewahlten Gaußkurven-Segmente naturlich nicht perfekt mit asymetrischen Wellenubereinstimmen. Durch entsprechende Korrekturen konnten also diese systematischen Feh-ler reduziert werden, was erwartungsgemaß zu einer weiteren Verbesserung der Genauigkeitfuhren wurde.Die Standardabweichungen der Schatzwerte liegen bei einer Abtastfrequenz von fS1 = 256 Hzbei geringen Rauschamplituden in der Großenordnung von 1 ms, was ebenfalls als außerst gutzu bewerten ist. Außerdem ist eindeutig ersichtlich, dass das uberlagerte Rauschen direktenEinfluss auf die Unsicherheit hat, was sich entsprechend besonders stark auf die schmalen Q-und S-Wellen auswirkt. Das Verdoppeln der Abtastfrequenz auf fS1 = 512 Hz hat in etwa eineHalbierung der Standardabweichungen zur Folge, auch die Mittelwerte bei Q- und S-Wellenverbessern sich merklich. Das verdeutlicht, dass die erzielbare Genauigkeit unmittelbar durchdie Wahl der Abtastrate beeinflusst werden kann.

12Die uber 100 Ergebnisse fur die Abweichungen der Mittelwerte weisen bis auf eine Ausnahme stetsdas selbe Vorzeichen fur die jeweiligen Teilwellen auf.

34

Parameter Schatzwerte (µ± σ)/ms Det.rate

HR/Hz σHR/Hz AN/µV P Q R S T DR

Abtastfrequenz: fS = 256 Hz, Anzahl der Herzzyklen: N ≈ 300

60 1,0 50 −0, 45± 1, 9 −0, 35± 4, 2 −0, 44± 1, 2 0, 16± 2, 8 −0, 80± 1, 7 100%60 1,0 50 −0, 43± 2, 1 0, 52± 4, 8 −0, 35± 1, 1 0, 28± 2, 6 −1, 03± 1, 9 100%60 1,5 20 −0, 37± 1, 4 −1, 58± 2, 0 −0, 36± 1, 1 0, 63± 1, 7 −1, 04± 1, 3 100%60 1,5 20 −0, 44± 1, 4 −1, 69± 1, 6 −0, 39± 1, 1 0, 65± 1, 8 −0, 94± 1, 3 100%80 2,0 10 −0, 30± 1, 2 −1, 74± 1, 1 −0, 24± 1, 1 0, 76± 1, 0 −0, 68± 1, 3 100%80 2,0 10 −0, 33± 1, 2 −1, 71± 1, 2 −0, 33± 1, 2 0, 84± 1, 0 −0, 72± 1, 2 100%120 2,5 10 −0, 46± 1, 2 −1, 58± 1, 2 −0, 23± 1, 1 1, 32± 1, 2 −0, 67± 1, 1 100%120 2,5 10 −0, 48± 1, 2 −1, 62± 1, 1 −0, 20± 1, 1 1, 24± 1, 2 −0, 64± 1, 2 100%150 3,0 20 −0, 81± 1, 3 −1, 74± 1, 1 −0, 18± 1, 2 1, 19± 1, 2 −0, 58± 1, 3 100%150 3,0 20 −0, 65± 1, 4 −1, 72± 1, 2 −0, 27± 1, 1 1, 25± 1, 2 −0, 61± 1, 2 100%180 5,0 50 −0, 45± 2, 4 −1, 44± 1, 8 −0, 30± 1, 2 1, 07± 1, 6 −0, 79± 2, 2 100%180 5,0 50 −0, 40± 2, 7 −1, 60± 1, 8 −0, 27± 1, 1 0, 96± 1, 5 −0, 74± 2, 2 100%200 5,0 0 −0, 33± 1, 1 −1, 41± 1, 2 −0, 20± 1, 1 1, 04± 1, 2 −0, 43± 1, 2 100%200 5,0 0 −0, 38± 1, 2 −1, 43± 1, 2 −0, 14± 1, 1 1, 02± 1, 1 −0, 44± 1, 2 100%

Abtastfrequenz: fS = 512 Hz, Anzahl der Herzzyklen: N ≈ 100

60 2,0 20 −0, 39± 1, 0 −0, 70± 0, 8 −0, 43± 0, 6 0, 12± 0, 7 −0, 98± 0, 9 100%80 4,0 100 −0, 39± 2, 3 −0, 43± 2, 7 −0, 45± 0, 6 0, 24± 2, 0 −0, 77± 2, 0 100%100 5,0 10 −0, 36± 0, 7 −0, 32± 0, 6 −0, 32± 0, 6 0, 42± 0, 6 −0, 65± 0, 7 100%120 2,0 0 −0, 29± 0, 6 −1, 01± 0, 6 −0, 19± 0, 6 0, 56± 0, 6 −0, 27± 0, 5 100%150 5,0 0 −0, 42± 0, 6 −0, 78± 0, 6 −0, 22± 0, 6 0, 35± 0, 5 −0, 19± 0, 6 100%180 3,0 50 −0, 53± 1, 5 −1, 14± 1, 9 −0, 29± 0, 5 0, 44± 1, 2 −0, 76± 1, 2 100%200 1,0 20 −0, 59± 0, 8 −0, 53± 0, 9 −0, 27± 0, 6 0, 42± 0, 7 −0, 96± 0, 7 100%

Tabelle 5.5: Statistische Auswertung mehrerer synthetischer Datensatze

35

5.4 PhysioNet QT-Datenbank

Um die Schatzung der Q-/T- Start- und Endpunkte ebenfalls bewerten zu konnen, wurdedie QT-Datenbank von PhysioNet [16] herangezogen. Tabelle 5.6 zeigt ein Beispiel fur dieStatistik-Ausgaben anhand der Auswertung des Datensatzes qtdb/sel103.

Geschatzte Herzfrequenz

HR = 70, 7 bpm

Marker-Schatzung durch Kovarianzfilterµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P 13.8213 2.2664 8.0000 20.0000 1 0Q 7.2340 2.0298 4.0000 16.0000 1 0R 3.8386 0.9450 0 8.0000 0 0S -7.0043 6.9057 -16.0000 4.0000 1 0T 0.7319 1.8122 -4.0000 4.0000 0 0

Marker-Schatzung durch Gaußkurven-Approximationµ/ms σ/ms ∆min/ms ∆max/ms Nicht det. Falsch pos.

P 8.0485 2.2260 -0.2760 15.3980 1 0Q 6.6030 1.7265 2.5732 15.1103 1 0R 3.9887 1.2020 1.9766 6.1585 0 0S -6.9514 6.4075 -15.3047 5.4379 1 0T 3.0784 1.4643 -1.2305 7.4452 0 0

P/T Start- und Endzeitpunkte (onset, offset)µ/ms σ/ms ∆min/ms ∆max/ms µref/ms [13] σref/ms [13]

Pon -1.1438 8.4517 -86.6240 11.2929 1.4 28.3Poff 25.1398 5.9522 -21.6699 38.4630 10.0 25.2Ton 20.4661 4.2026 10.7243 45.9978 -1.4 45.8Toff -12.2829 3.6949 -23.4218 -2.1074 2.8 34.9

Tabelle 5.6: Analyse von qtdb/sel103 (Samples 1 bis 100000)

Die schlechteste Genauigkeit der Mittelwerte erreicht hier µPoff ≈ 25 ms, die Standardabwei-chungen liegen im Bereich von SPon ≈ 8 ms oder besser. Zum Vergleich werden die Werteaus [13] fur die Analyse der selben Datenbank herangezogen - diese betragen µPoff = 10 msund SPon ≈ 28 ms. Der hier beschriebene Algorithmus liefert also wesentlich bessere Stan-dardabweichungen, die Mittelwerte zeigen hingegen starkere Fehler. Durch ausfuhrliche Testsin Form von manueller Auswertung zahlreicher Problemfalle konnte beobachtet werden, dassdiese Ungenauigkeit sehr stark von der Form der Wellen abhangt, die hier durch Gaußkurvenapproximiert werden. Daher kann das Ergebnis durch Anpassen einfacher Skalierungsparameterfur eine typische Signalform optimiert werden, was jedoch fur andere Signale zu schlechterenResultaten fuhren kann. Aus diesem Grund wurde auf diese Art der Optimierung verzichtet, dasie als “Ergebniskosmetik” angesehen wurde. Sobald echte EKG-Daten von vielen verschiede-nen Personen mit gesetzten Markern zur Verfugung stehen, kann ein Parametersatz ermitteltwerden, der im Durchschnitt die besten Resultate liefert. Alternativ konnte die Optimierungfur die jeweilige Wellenform durch manuelles Setzen einiger Marker durchgefuhrt werden.

36

6 Fazit

Die beschriebene Methode der EKG-Signalanalyse liefert verhaltnismaßig geringe Standardab-weichungen fur die ermittelten “Marker”, welche die Teilwellen sowohl anhand der Positionender Wellenberge als auch durch Start- und Endpunkte charakterisiert. Auch die Mittelwerte derSchatzungen weichen in einer Großenordnung ab, die verglichen mit anderen Veroffentlichungenals “state of the art” bezeichnet werden kann.Somit konnte gezeigt werden, dass die Detektion der EKG-Wellenberge durch relativ einfacheZeitbereichsfilter sehr zuverlassig und genau bewerkstelligt werden kann, deren echtzeitfahigeImplementierung auf beinahe jedem handelsublichen Signalprozessor moglich sein sollte. DieLeast-Squares-Approximierung von Gaußkurvensegmenten erweist sich hingegen fur langereSignalteile als verhaltnismaßig rechenintensiv und fur sehr kurze Abschnitte als etwas unge-nau und weniger zuverlassig. Bei der A-posteriori-Auswertung synthetischer “Laborsignale” aufeinem rechenstarken Prozessor stellt dies kein Problem dar - die Eignung fur die Echtzeittaug-lichkeit wird sich bei Tests mit echten Messdaten zeigen.Insgesamt fiel die Auswertung des implementierten Programms sehr zufriedenstellend aus,allerdings fehlt noch ein umfangreicher Test mit echten Patientendaten. Hiermit ist eine viel-versprechende Basis fur den Ubergang auf Wavelet-basierte Methoden geschaffen.

7 Ausblick

Bei der Programmierung lag der Fokus zunachst auf einer moglichst intuitiven und ein-fach nachvollziehbaren Losung, weshalb noch keine Algorithmen fur Lernfahigkeit, automa-tische Optimierung oder zur ruckwirkenden Fehlerkorrektur fur schlechte Ergebnisse vorgese-hen wurden. Die in Kapitel 5 beobachteten systematischen Fehler werden dzt. nicht kom-pensiert. Außerdem wurden viele a-priori-Parameter nur anhand einiger empirischer Testsgeschatzt oder aus anderen Veroffentlichungen ubernommen, jedoch nicht durch eine umfang-reichere Testreihe mit echten EKG-Daten optimiert. Grund dafur ist, dass diese Programm-routine hauptsachlich als Basis bzw. Referenz fur andere Methoden wie z. B. der Wavelet-Transformation vorgesehen war. Aufgrund der sehr genauen Analyseergebnisse ist jedoch auchnicht ausgeschlossen, das Programm tatsachlich fur die automatisierte Auswertung von EKG-Signalen heranzuziehen, besonders in Anbetracht des angefuhrten Verbesserungspotentials.Um auch die medizinisch relevanten Aspekte einzubinden, ware dafur jedenfalls umfangreicheKooperation mit kardiologisch geschultem Fachpersonal von Noten.Eine vielversprechende Verbesserung ware die Erweiterung der implementierten Losung aufein zwei-Kanal-EKG. Mit nur einem Kanal (einer “Ableitung”) ist die Gefahr groß, dass dieSignalform stark verfalscht wird, z. B. aufgrund von Winkelversatz. In diesem Fall funktioniertder Algorithmus nur noch sehr unzuverlassig. Auch die Anfalligkeit gegen starken Nullinien-drift ist sehr hoch. Vermutlich konnte die Erweiterung auf zwei (oder mehr) Ableitungen dieZuverlassigkeit erheblich steigern und evtl. sogar einen Gewinn an Genauigkeit mit sich brin-gen. Besonders interessant wird es, wenn Routinen auf Grundlage der Wavelet-Transformationeingebunden werden, da man z. B. mit einer komplexen Wavelet-Funktion zwei (komplexe)Ableitungen gleichzeitig auswerten kann.

37

Literatur

[1] MATLAB and Statistics Toolbox Release 2010b, The MathWorks, Inc., Natick, Massachu-setts, United States

[2] Daniel E. Becker, “Fundamentals of Electrocardiography Interpretation”, Anesthesia Pro-gress 53.2 (2006): 53-64

[3] John R. Hampton, “The ECG Made Easy”, ISBN 0443072523, Churchill Livingstone, May2003

[4] Alfred Haar, “Zur Theorie der orthogonalen Funktionensysteme”, Mathematische AnnalenBand 69, 1910

[5] Alexandre Bernardino and Jose Santos-Victor, “A Real-Time Gabor Primal Sketch forVisual Attention”, IBPRIA - 2nd Iberian Conference on Pattern Recognition and ImageAnalysis, June 2005

[6] Ingrid Daubechies, “Orthonormal Bases of Compactly Supported Wavelets”, Communica-tions on Pure and Applied Mathematics 41, 1988

[7] Paul S Addison, “The Illustrated Wavelet Transform Handbook: Introductory Theory andApplications in Science Engineering, Medicine and Finance”, Bristol: Institute of PhysicsPublishing, 2002

[8] Paul S Addison, “Wavelet transforms and the ECG: a review”, IOP Publishing Ltd, Phy-siological Measurement 26, August 2005

[9] Stephane G. Mallat, “A theory for multiresolution signal decomposition: the wavelet repre-sentation”, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTEL-LIGENCE, VOL. II, NO. 7. JULY 1989

[10] Jiapu Pan and Willis J. Tompkins, “A Real-Time QRS-Detection Algorithm”, IEEE,TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-32, NO. 3, MARCH1985

[11] Patrick E. McSharry, Gari D. Clifford, Lionel Tarassenko, Leonard A. Smith, “A DynamicalModel for Generating Synthetic Electrocardiogram Signals”, IEEE TRANSACTIONS ONBIOMEDICAL ENGINEERING, VOL. 50, NO. 3, MARCH 2003

[12] N. V. Thakor, J. G. Webster, and W. J. Tompkins, “Optimal QRS detector”, MED BIOLENG COMPUT, vol. 21, no. 3, pp. 343-350, 1983

[13] D. Hayn, A. Kollmann, G. Schreier, “Automated QT Interval Measurement from MultileadECG Signals”, ISSN 0276-6547, Computers in Cardiology 2006;33:93-96

[14] A. G. Ramakrishnan, “Threshold-Independent QRS Detection Using the Dynamic PlosionIndex”, IEEE SIGNAL PROCESSING LETTERS, VOL. 21, NO. 5, MAY 2014

[15] Hooman Sedghamiz, Linkoping University, hooman.sedghamiz@medel.com

[16] http://www.physionet.org/

38

Abbildungsverzeichnis

1.1 Nomenklatur des EKG-Signals . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Haar- und Mexican-Hat-Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Morlet-Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3 Daubechies-Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.4 Beispiele der CWT fur Einheitsimpuls δ(t) . . . . . . . . . . . . . . . . . . . 102.5 CWT und DWT von Rechteck- und e-Funktion . . . . . . . . . . . . . . . . . 132.6 CWT und DWT eines Chirp-Signals . . . . . . . . . . . . . . . . . . . . . . . 152.7 DWT und CWT eines EKG-Signals im Vergleich . . . . . . . . . . . . . . . . 163.1 Ausgabe des Pan-Tompkins-Algorithmus (MATLAB-Code [15]) . . . . . . . . 173.2 Einzelstufen des Pan-Tompkins-Algorithmus . . . . . . . . . . . . . . . . . . 193.3 Preprocessing: Unterdruckung der QRS-Komplexe . . . . . . . . . . . . . . . 213.4 Ergebnis der Kovarianzfilterung fur P- und S-Wellen . . . . . . . . . . . . . . 233.5 Gaußfilter fur P-Wellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.6 Reprasentation der Teilwellen durch Gaußkurven-Segmente . . . . . . . . . . . 253.7 Referenz- und detektierte Marker im Vergleich . . . . . . . . . . . . . . . . . 264.1 GUI: Waveform-Viewer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295.1 Wellenerkennung des verrauschten Testsignals . . . . . . . . . . . . . . . . . 33

Tabellenverzeichnis

3.1 Parametersatz fur ecgsyn.m . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2 Relative Standardabweichungen der Gauss-Kovarianzfilter . . . . . . . . . . . 223.3 Standardabweichungen fur Gaußkurven-Approximation . . . . . . . . . . . . . 245.1 Parametersatz 1 fur ecgsyn.m . . . . . . . . . . . . . . . . . . . . . . . . . 305.2 Analyse des rauschfreien synthetischen EKG . . . . . . . . . . . . . . . . . . 315.3 Parametersatz 2 fur ecgsyn.m . . . . . . . . . . . . . . . . . . . . . . . . . 325.4 Analyse des verrauschten synthetischen EKG . . . . . . . . . . . . . . . . . . 325.5 Statistische Auswertung mehrerer synthetischer Datensatze . . . . . . . . . . 355.6 Analyse von qtdb/sel103 (Samples 1 bis 100000) . . . . . . . . . . . . . . . 36

39