the Hilbert transform to build a Holter Monitoring System is proposed. Using the Parks-McClellan algorithm, the Remez approach was present, and a digi...

0 downloads 31 Views 1MB Size

No documents

Loading...

A Thesis Submitted to the College of Graduate Studies and Research in Partial Fulfillment of the Requirements For the Degree of Master of Science in the Department of Electrical Engineering University of Saskatchewan Saskatoon, Saskatchewan

by Xiangling Wang

© Copyright Xiangling Wang, September 2006, All Rights Reserved

PERMISSION TO USE I agree that the Library, University of Saskatchewan, may make this thesis freely available for inspection. I further agree that permission for copying of this thesis for scholarly purpose may be granted to the professor or professors who supervised the thesis work recorded herein or, in their absence, by the Head of the Department or the Dean of the College in which the thesis work was done. It is understood that due recognition will be given to me and to the University of Saskatchewan in any use of the material in this thesis. Copying or publication or any other use of this thesis for financial gain without approval by the University of Saskatchewan and my written permission is prohibited. Request for permission to copy or to make any other use of the material in this thesis in whole or part should be addressed to:

Head of the Department of Electrical Engineering 57 Campus Drive University of Saskatchewan Saskatoon, Saskatchewan Canada S7N 5A9

i

ACKNOWLEDGEMENTS I wish to express my gratitude to the following people who not only made this thesis possible but also an enjoyable experience: Dr. Ronald J. Bolton: my supervisor, for his valuable guidance, criticisms and consistent encouragement throughout the course of this research work. My husband, Zhanghai Wang: for his love and encouragement. My parents, my sister Xiangrong Wang and brother-in-law Xianggang Yu: for the support they provided to me. All my friends in Saskatoon, Sha Li, Song Hu, Jing Yin, Ying Cui, Quan Wan, Yajun Wang, Xin Xu and Yanan Xing: for making me feel welcome. They will always be special friends in my life. The Department of Electrical Engineering: for supplying the opportunity to study in Canada and the necessary facilities with which to work.

ii

ABSTRACT Many people have abnormal heartbeats from time to time. A Holter monitor is a device used to record the electrical impulses of the heart when people do ordinary activities. Holter monitoring systems that can record heart rate and rhythm when you feel chest pain or symptoms of an irregular heartbeat (called an arrhythmia) and automatically perform electrocardiogram (ECG) signal analysis are desirable. The use of the Hilbert transform (HT) in the area of electrocardiogram analysis is investigated. A property of the Hilbert transform, i.e., to form the analytic signal, was used in this thesis. Subsequently pattern recognition can be used to analyse the ECG data and lossless compression techniques can be used to reduce the ECG data for storage. The thesis discusses one part of the Holter Monitoring System, Input processing. Four different approaches, including the Time-Domain approach, the Frequency-Domain approach, the Boche approach and the Remez filter approach for calculating the Hilbert transform of an ECG wave are discussed in this thesis. By comparing them from the running time and the ease of software and hardware implementations, an efficient approach (the Remez approach) for use in calculating the Hilbert transform to build a Holter Monitoring System is proposed. Using the Parks-McClellan algorithm, the Remez approach was present, and a digital filter was developed to filter the data sequence.

iii

Accurate determination of the QRS complex, in particular, accurate detection of the R wave peak, is important in ECG analysis and is another task in this thesis. A program was developed to detect the R wave peak in an ECG wave. The whole algorithm is implemented using Altera’s Nios SOPC (system on a program chip) Builder system development tool. The performance of the algorithm was tested using the standard ECG waveform records from the MIT-BIH Arrhythmia database. The results will be used in pattern recognition to judge whether the ECG wave is normal or abnormal.

iv

Table of Contents PERMISSION TO USE ............................................................................................. i ACKNOWLEDGEMENTS ....................................................................................... i ABSTRACT ............................................................................................................. iii Table of Contents ...................................................................................................... v List of Tables........................................................................................................... vii List of Figures......................................................................................................... viii List of Abbreviations................................................................................................ xi Chapter 1 Introduction............................................................................................... 1 1.1 Research Motivation........................................................................................ 1 1.1.1 Electrocadiogram...................................................................................... 2 1.1.2 Advantanges of Holter Monitoring System.............................................. 6 1.1.3 Structure of the Holter Monitoring System .............................................. 8 1.2 Objectives ...................................................................................................... 10 1.3 Outline of Thesis ........................................................................................... 11 Chapter 2 Background............................................................................................. 13 2.1 Holter Monitoring Review ............................................................................ 13 2.2 Hilbert Transform Review............................................................................. 16 2.2.1 Definition................................................................................................ 17 2.2.2 Frequency Response of the Hilbert Transform ...................................... 18 2.3 Hilbert Transform Properties......................................................................... 20 2.4 Hilbert Transform Applications .................................................................... 23 2.4.1 Analytic Signal ....................................................................................... 23 2.4.2 Analytic Signal Applied in Pattern Recognition [11]............................. 24 Chapter 3 Computation of the Hilbert Transform ................................................... 26 3.1 Time-Domain Approach................................................................................ 26 3.2 Frequency-Domain Approach ....................................................................... 36 3.3 Boche Approach ............................................................................................ 41 3.4 Remez Approach ........................................................................................... 46

v

3.5 Comparison.................................................................................................... 51 Chapter 4 Implementation ....................................................................................... 53 4.1 Nios Embedded Processor Overview ............................................................ 53 4.2 Digital Filter .................................................................................................. 59 4.3 Implementation.............................................................................................. 63 4.3.1 Filter Order ............................................................................................. 63 4.3.2 Filter Coefficients................................................................................... 69 4.3.3 Digital Filter ........................................................................................... 77 4.3.4 Detector for R Wave Peak...................................................................... 78 4.3.5 R Wave Peak Detection Test................................................................. 79 4.4 Nios Implementation .................................................................................... 84 Chapter 5 Results..................................................................................................... 91 5.1 Experimental Results..................................................................................... 91 5.2 Complete ECG Testing................................................................................ 104 Chapter 6 Summary and Conclusion..................................................................... 112 6.1 Summary…………………………………………………………………...112 6.2 Conclusion…………………………………………………………………115 References ............................................................................................................. 116 Appendix A ........................................................................................................... 119 The MIT-BIH Arrhythmia Database ..................................................................... 119 A.1 Introduction ................................................................................................ 119 A.2 File Structure .............................................................................................. 120 A.3 Notational and Other Conventions ............................................................. 121 A.4 File Format Specifications.......................................................................... 122 A.5 Annotation Specifications........................................................................... 124 Appendix B............................................................................................................ 128 B.1 The Parks-McClellan Algorithm ................................................................ 128

vi

List of Tables Table 1.1 Comparison general Holter Monitor with new system in this thesis ........ 7 Table 3.1 The comparison of the four methods for computing Hilbert transform. . 51 Table 4.1 Four types of the linear phase FIR filter.................................................. 62 Table 4.2 Filter order comparison ........................................................................... 67 Table 4.3 XW_1 extreme points and values............................................................ 81 Table 4.4 XW_2 extreme points and values............................................................ 84 Table 4.5 XW_1 R wave points and their Hilbert transform points........................ 90 Table 4.6 XW_2 R wave points and their Hilbert transform points........................ 90 Table 5.1 mit212from100 extreme points and values ............................................. 94 Table 5.2 mit212_1from47200 extreme points and values ..................................... 95 Table 5.3 mit213_1from67000 extreme points and values ..................................... 98 Table 5.4 mit213_1from1450 extreme points and values ..................................... 100 Table 5.5 mit223_1from47700 extreme points and values ................................... 102 Table 5.6 mit223_1from7900 extreme points and values ..................................... 103 Table 5.7 R wave detection performance .............................................................. 110

vii

List of Figures Figure 1.1 ECG signal [1] ......................................................................................... 2 Figure 1.2 The different components of the QRS complex [2]. ................................ 3 Figure 1.3 A normal ECG wave [3]. ......................................................................... 5 Figure 1.4 Bradycardia [3]. ....................................................................................... 5 Figure 1.5 Tachycardia [3]. ....................................................................................... 5 Figure 1.6 An irregular heartbeat wave [3]. .............................................................. 5 Figure 1.7 Three abnormal ECG waveforms [3]....................................................... 6 Figure 1.8 The structure of the Holter Monitoring System. ...................................... 8 Figure 2.1 A man with the Holter monitor [5]. ....................................................... 15 Figure 2.2 The Hilbert transform of a square wave………………………………..18 Figure 3.1 Input waveform: sin(2π * 0.02 * 500) ..................................................... 32

Figure 3.2 (a) Output waveform: The Hilbert transform of sin(2π * 0.02 * 500) . .. 32 Figure 3.2 (b) Output waveform using the equation from [12]……………………33 Figure 3.3 A normal ECG wave. ............................................................................. 34 Figure 3.4 The Hilbert transform of the ECG wave................................................ 35 Figure 3.5 Input wave: sin(2π * 0.02 * 500) . .......................................................... 38 Figure 3.6 The Hilbert transform of sin(2π * 0.02 * 500) . ...................................... 38 Figure 3.7 Input ECG wave..................................................................................... 39 Figure 3.8 Output: the Hilbert transform of the ECG wave. ................................... 40 Figure 3.9 (a) f (t )

(b) f 26 (c) f 26 (t ) − f (t ) (d) t i . ......................................... 44

Figure 3.10 (a) fˆ (t ) (b) fˆ26 (t ) (c) fˆ26 (t ) − fˆ (t ) .................................................... 45

Figure 3.11 The input wave: y = sin(2π * 0.02 * 500) . ........................................... 48 Figure 3.12 The output wave: HT of y = sin(2π * 0.02 * 500) . .............................. 48 Figure 3.13 An ECG wave. ..................................................................................... 50 Figure 3.14 The Hilbert transform of an ECG wave. .............................................. 50 Figure 4.1 Nios Embedded Processor System......................................................... 54

viii

Figure 4.2 Hardware/Software development flow for a Nios processor system [15] ................................................................................................................................. 56 Figure 4.3 Nios SDK Shell (bash)........................................................................... 57 Figure 4.4 Nios Development Board Components [17].......................................... 58 Figure 4.5 Filter. ...................................................................................................... 60 Figure 4.6 Digital filter transposed structure........................................................... 60 Figure 4.7 Illustration of four types of impulse response symmetry....................... 62 Figure 4.8 Input waveform: sin(2 * π * 0.02 * 500) .................................................. 64 Figure 4.9 Output waveform, M = 51 . .................................................................... 64 Figure 4.10 Output waveform, M = 71 . .................................................................. 65 Figure 4.11 Output waveform, M = 91 . .................................................................. 65 Figure 4.12 Output waveform, M = 101 . ................................................................. 66 Figure 4.13 Output waveform, M = 201 .................................................................. 66 Figure 4.14 Frequency response, M = 100 ............................................................. 68 Figure 4.15 Frequency response, M = 101 . ............................................................ 68 Figure 4.16 Frequency response of the ideal and Remez design filter.................... 70 Figure 4.17 Coefficients of the filter when M = 100 . ............................................ 71 Figure 4.18 Frequency response of Remez filter..................................................... 72 Figure 4.19 Initial parameter of remez.c ................................................................ 75 Figure 4.20 Coefficients from the C program. ........................................................ 76 Figure 4.21 Frequency response for a M = 100 Hilbert transform filter................ 77 Figure 4.22 Input wavefom: XW_1......................................................................... 80 Figure 4.23 Output waveform for XW_1. ............................................................... 81 Figure 4.24 Input waveform: XW_2. ...................................................................... 82 Figure 4.25 Output waveform for XW_2. ............................................................... 83 Figure 4.26 Nios SDK Shell Prompt. ...................................................................... 85 Figure 4.27 Nios-Build messages............................................................................ 86 Figure 4.28 Nios SDK shell prompt. ....................................................................... 87 Figure 4.29 Input ECG waveform: XW_1. ............................................................. 88 Figure 4.30 The Hilbert transform of XW_1........................................................... 88 Figure 4.31 Input ECG waveform: XW_2. ............................................................. 89 ix

Figure 4.32 The Hilbert transform of XW_2........................................................... 89 Figure 5.1 Input waveforms: mit212from100.txt. ................................................... 93 Figure 5.2 Output waveform: outMIT212from100.txt............................................ 93 Figure 5.3 Input waveform: mit212_1from47200.txt.............................................. 94 Figure 5.4 Output waveform: outMIT212from47200.txt........................................ 95 Figure 5.5 Input waveform: mit213_1from67000.txt.............................................. 97 Figure 5.6 Output waveform: outMIT213_1from67000.txt.................................... 98 Figure 5.7 Input waveform: mit213_1from1450.txt................................................ 99 Figure 5.8 Output waveform: outMIT213_1from1450.txt...................................... 99 Figure 5.9 Input waveform: mit223_1from47700.txt............................................ 101 Figure 5.10 Output waveform: outMIT223_1_from47700.txt.............................. 101 Figure 5.11 Input waveform: mit223_1from7900.txt............................................ 102 Figure 5.12 Output waveform: outMIT223_1_from7900.txt................................ 103 Figure 5.13 The portion of results of MIT212. ..................................................... 106 Figure 5.14 The portion of results of MIT213. ..................................................... 107 Figure 5.15 The portion of results of MIT213. ..................................................... 108 Figure 5.16 The portion of results of MIT223. ..................................................... 109 Figure 5.17 The portion of results of MIT223. ..................................................... 110 Figure B.1 Flowchart of Parks-McClellan algorithm………………………….....134

x

List of Abbreviations APBs

atrial premature beats

BPM

beats per minute

DHCP

dynamic host configuration protocol

ECG

electrocardiogram

FIR

finite impulse response

FT

Fourier transform

FFT

fast Fourier transform

FPGA

field programmable gate array

HW

hardware

HT

Hilbert transform

IFFT

inverse fast Fourier transform

ISA

instruction set architecture

LED

light emitting diode

MIT-BIH

Massachusetts Institute of Technology- Beth Israel Hospital

PCG

polarcardiogram

PVC

premature ventricular contraction

PHY/MAC

physical layer/media access control

RISC

reduced instruction set computer

RMS

root mean square

SDK

software development kit

SOPC

system on a program chip

xi

SW

software

VCG

vectorcardiogram

VPB

ventricular premature beat

xii

Chapter 1 Introduction

Since the development of medical science, many instruments for improving people’s health have been developed. The electrocardiogram (ECG) monitoring system is one of them. The most common type of ECG monitoring is called Holter monitoring. Holt monitoring is a portable recording tool and can help doctors make a precise diagnosis.

1.1 Research Motivation Many people have irregular heartbeats from time to time. Some heart problems occur during certain activities, such as eating, exercise or even sleeping. Sometimes the irregular heartbeats don’t influence life style and are usually harmless in normal hearts. But it is also possible that these irregular heartbeats with pre-existing illness can cause heart attacks that lead to death. A device that can record the activities of the heart is very useful in preventing heart attacks. The Holter monitoring system was developed for this objective in an ambulatory situation. A Holter monitoring system is a small recording instrument that is used to capture ECG data of the heart’s electrical activities over a period of time. The

1

patient can carry it in a pocket or in a small pouch. The monitor is battery operated. The electrocardiogram is saved in a memory card. The electrocardiographer can analyse the recordings visually by means of a computer.

1.1.1 Electrocadiogram The standard ECG is a representation of the heart’s electrical activities recorded from electrodes that are placed on different parts of patient’s body. The electrocardiogram is composed of complexes and waves. In normal sinus rhythm, waves and complexes are the P wave, QRS complex, ST Segment, T wave and U wave. Measurements are PR interval, QRS duration, QT interval, RR interval and PP interval. Figure 1.1 illustrates a typical waveform of normal heartbeats and intervals as well as standard time and voltage measures.

Figure 1.1 ECG signal [1]

2

Different parts of the ECG waves are caused differently. Detailed information will be given below to explain each part of the ECG waveform. • The P wave is due to the electrical activation (depolarization) of the heart

(atria). It is usually positive, low amplitude and smooth. In normal situation, the time of the P wave should be smaller than 0.12 seconds. • The QRS complex represents right and left ventricular depolarization. It is

high amplitude in normal situations. The shape of the QRS complex will be changed if the electrodes are placed on different parts of the body. It also changes when abnormal heartbeats occur. A QRS complex can have positive (upwards) or negative (downwards) deflections. The figure below summarizes the nomenclature used to define the different components of the QRS complex.

Figure 1.2 The different components of the QRS complex [2].

• The ST segment represents the time following the QRS it takes for

depolarization of the ventricles before repolarization. Repolarization of the atria is a low amplitude signal that occurs during the time of the high amplitude QRS and consequently it can’t be seen on a standard ECG.

3

• The T wave is caused by the repolarization of the ventricles. Usually it is

positive and rounded. • The reason that causes the U wave is not that clear, “afterdepolarizations” in

the ventricles maybe is the answer. • The PR interval is the time interval from the beginning of the P wave to the

beginning of the QRS complex. In normal situation the PR interval should be 0.120.2 s, Short PR < 0.12 s, Prolonged PR >0.2 s. • The QRS duration is the time of ventricular depolarization. Normal: 0.06 s-

0.1 s, Prolonged QRS duration: >0.1s. • The QT interval represents the duration of ventricular depolarization and

repolarization. It is between the onset of the QRS complex and the end of the T wave. It normally depends on heart rate. • The RR interval is the duration of ventricular cardiac cycle. The value of

the RR interval indicates the ventricular rate. • The PP interval is the duration of atrial cycle. It indicates the atrial rate.

The normal adult heart beats regularly between 60 to 100 beats per minute. Bradycardia occurs once the heart rate is slower than 60 beats per minute. The waveform is similar to the normal ECG wave, but the RR interval is longer. A rate of above 100 beats per minute is called tachycardia, in this case the RR interval is shorter and the waveform is also similar to the regular sinus rhythm. Each P wave is following by a QRS complex. A waveform of a regular ECG wave is shown in

4

Figure 1.3. The wave of the bradycardia is shown in Figure 1.4. Figure 1.5 illustrates the tachycardia. Figure 1.7 shows some abnormal ECG waves.

Figure 1.3 A normal ECG wave [3].

Figure 1.4 Bradycardia [3].

Figure 1.5 Tachycardia [3].

Figure 1.6 An irregular heartbeat wave [3].

5

Figure 1.7 Three abnormal ECG waveforms [3].

The importance of irregular heartbeats depends on the type of pattern they produce, how often they occur, how long they last, and whether they occur at the same time the patient had symptoms.

1.1.2 Advantanges of Holter Monitoring System As discussed previously, the Holter monitoring system records the electrical activity of heart during usual daily activities. A recording is much more likely to detect any abnormal heartbeats that occur during these activities. During the late 1960s, computerized ECG's came into use in many of the larger

hospitals.

General

Holter

monitoring

system

records

continuous

electrocardiographic measurements of the heart’s rhythm. Usually the recording

6

time is around 24 to 48 hours. That means even when the heartbeat is normal, the Holter monitor also works as well. The system discussed in this thesis automatically records the ECG wave when the user is not feeling good or the heartbeats are not regular. The recording algorithm is not continuous any more. It also can record the heartbeats manually; the wearers can record the heartbeat if wanted when the heart rhythms are ordinary. The differences between the general Holter monitor and the system developed in this thesis are shown in the Table 1.1. Table 1.1 Comparison general Holter Monitor with new system in this thesis General Holter Monitor

New system in this thesis

Continuity

Continuous

Intermittent

Saving time

24-48 hours

More than 48 hours

Operation

Manually operated

Automatic/Manually operated

The advantages in the system discussed in this thesis are: • Record the ECG wave automatically when wearer does not feel good. • Save memory space and extend the recording time. • Record activities manually when the wearer wants. • Normally record for a few hours or for a few seconds.

7

1.1.3 Structure of the Holter Monitoring System The structure of the Holter monitoring system discussed in this thesis is shown in Figure 1.8.

Raw ECG Data

Input Processing (Hilbert Transform)

Compression Nios

Nios

Storage Compact Flash Card (Removable)

User Interface Pattern Recognition Normal/Abnormal

Abnormal

Nios

Figure 1.8 The structure of the Holter Monitoring System.

From Figure 1.8 it can be seen that the system includes four sub-systems: the Input processing sub-system, the Pattern Recognition sub-system, the Compression sub-system and the Storage sub-system. The input data is the raw ECG data. These data record the activities of the heart. Every heartbeat is caused by a section of the heart generating an electrical signal that then conducts through specialized pathways to all parts of the heart.

8

These electrical signals also get transmitted through the chest to the skin where they can be recorded [4]. 1.1.3.1 Input Processing

The main objective of this sub-system is to implement the Hilbert transform of the input ECG data on the Nios embedded processor. Then the zero crossing points corresponding to the input R wave peaks are found and the results information is sent to the Pattern Recognition sub-system. 1.1.3.2 Pattern Recognition

Obtain the results of the input processing, and by using properties of the Hilbert transform, create the analytic signal and use it to assist in doing pattern recognition to determine if the ECG wave is normal or abnormal. 1.1.3.3 Compression

Compress all the abnormal data obtained from the “Input processing” subsystem. 1.1.3.4 Storage

Save all the compressed data in compact flash card. This card is removable and can be inserted into the computer to read the information recorded in this card. Electrocardiographer can analyze the data or by means of some software to draw the waveform to see what happened visually.

9

Holter monitoring gives doctors the record of patient’s heart rate and rhythm over a period of time. The Holter monitor can record heart rate and rhythm when the patient feels chest pain or symptoms of an irregular heartbeat (an arrhythmia). The doctor can then look at the time when the patient noticed their symptoms. Reading this printout will give the doctor an idea about the nature of the heart problem. A Holter monitor provides the physician a better opportunity to capture any abnormal heartbeats or rhythms that may be causing the patient’s symptoms. It’s necessary to develop an efficient (i.e., long time durations and automatic pattern recognition analysis) system to record these irregular heartbeats so that the doctor can know what had happened when patient had those symptoms. This system is a very good tool to prevent people from fatal heart problems. Briefly, the reasons for using a Holter monitor may include: to detect problems missed in a regular ECG; to check activity after an arrhythmia; to see if a new pacemaker works and to see if drug therapy is working. The work of this thesis is to implement the input processing sub-system. That is to develop an efficient algorithm to compute the Hilbert transform of the input ECG waves and determine R-R intervals.

1.2 Objectives The objectives of this thesis are to address the aforementioned problems and to propose solutions to the following problems. z To use different methods to compute the Hilbert transform of an input signal.

10

z To compare these algorithms with each other from running time and

hardware/software implementation. z To develop an efficient algorithm for calculating the Hilbert transform to build

a Holter ECG Monitoring System. z To develop a detector to find all of the zero crossing points that correspond to

the R wave peaks in the output wave, i.e., the Hilbert transform sequence of an ECG wave.

1.3 Outline of Thesis Chapter 2 gives the background materials of the Holter monitor and the Hilbert transform including the past history, notation and definition. Chapter 3 is devoted to four different methods used to compute the Hilbert transform of an input signal. The four approaches include the Time-Domain approach, the Frequency-Domain approach, the Boche approach and the Remez filter approach. In this chapter, the examples and the results are given for using different approaches to calculate the Hilbert transform. By comparing four methods with each other in running time and the ease of software and hardware implementations, an efficient algorithm for the Hilbert transform to build a Holter ECG Monitoring System will be present. Chapter 4 describes how to design and implement the Remez filter approach that was mentioned in Chapter 3, and also describes how to design and implement a filter with an optimal fit. The information about Alter’s Nios SOPC Builder system development tool is introduced as well. 11

Chapter 5 concentrates on analyzing the results and compares them to the MIT-BIH Arrhythmia database to make sure the results are correct. The thesis is summarized in Chapter 6. At the end of the thesis, an appendix is given.

12

Chapter 2 Background

In this chapter, a brief review of Holter Monitoring is given. The definition, the properties and the applications of the Hilbert transform are also contained in this chapter.

2.1 Holter Monitoring Review As discussed in the previous chapter, an electrocardiogram (ECG) is a record produced by an electrocardiography, which indicates the electrical voltage in the heart. An ECG provides information on the condition and performance of the heart. It is one of the simplest and fastest procedures used to evaluate the heart. Because an arrhythmia can occur irregularly, it will be difficult to record when the patient is in the doctor’s office. In 1949, American physician Norman Jeff Holter (1914-1983) developed a 75 pound backpack that can record the ECG of the wearer. This portable monitoring device is called the Holter monitor, named after its inventor. The Holter monitor is battery-powered and can continuously record the electrical activities of the heart over a specified period of time, normally 24 to 48 hours. Usually the patient will undergo Holter monitoring as an outpatient, meaning 13

that the monitor will be placed on the body of the patient by a technician in a cardiologist’s office. Then the patient will go home and do normal activities. With the development of technology, the Holter monitor is greatly reduced in size. It is now very compact and combined with digital recording and used to record ECGs. The Holter monitor can be easily carried without interfering with the patient’s activities. At the end of the recording period, the patient will go back to the doctor’s office to remove the Holter monitor. The data saved in the Holter monitor will be analyzed by an electrocardiographer and a computer. The analysis results will provide the information about the patient’s heart rhythm, the frequency of the beats and the irregularities. This portable monitor can be an effective and powerful diagnostic tool that can directly determine how the physician treats the patient’s condition. As Figure 2.1 shows, the Holter monitor is a small-size recording device. The monitor has wires called leads. The leads attach to metal disks called electrodes, which the user wears on his chest. These electrodes are very sensitive, and they can pick up the electrical impulses of the heart. The impulses are recorded by the Holter monitor record the heart’s electrical activity.

14

Figure 2.1 A man with the Holter monitor [5].

If necessary, the ECG can be transmitted by telephone to a computer at the hospital or doctor's office for an immediate reading as soon as symptoms occur. The use of the effective home care monitor of the heart patients will decrease the incidence of the readmissions and lower the costs of the hearth care. Advanced Holter monitors have been developed that employ digital electrocardiographic recordings, extended memory for more than 24 hours recording, pacemaker pulse detection and analysis, software for analysis of digital ECG recordings that are downloaded and stored on a computer, and capability of transmission of results over the internet [6]. The system discussed in this thesis is to record the rhythm of the heartbeat automatically when the symptoms occur. The monitor does not continuously record

15

but records the heart rate and rhythm when the patients feel symptoms of an irregular heartbeat or when abnormal heart beats or rhythms occur. That means the device automatically capture the arrhythmias when they occur. It also can be activated manually by the patients when chest pain is felt during a symptomatic event. So it can record for a long time. Up to now, the Holter monitor in the market can record over 24-48 hours; the longest one is not much more than 72 hours. The system discussed in this thesis will much improve the recording time.

The

advantages of the system have been discussed in Chapter 1.

2.2 Hilbert Transform Review In 1893, the physicist Arthur E. Kennelly (1861-1939) and the scientist Charles P. Steinmetz (1865-1923) first used the Euler formula e jz = cos( z ) + j sin( z )

(2.1)

which was derived by a famous Swiss mathematician Lenonard Euler (1707-1783) to introduce the complex notation of harmonic wave forms in electrical engineering, that is: e jwt = cos(ωt ) + j sin(ωt ) ,

(2.2)

where j is the imaginary unit. In the beginning of the 20th century, the German scientist David Hilbert (1862-1943) proved that the Hilbert transform of the function cos(ωt ) is sin(ωt ) . This is the one of properties of the Hilbert transform, i.e., basic ±

16

π 2

phase-shift.

2.2.1 Definition In mathematics and in signal processing, the Hilbert transform xˆ (t ) of a real time function x(t ) is defined as: xˆ (t ) = H [ x(t )] =

1

π∫

∞

−∞

x(τ ) dτ t −τ

(2.3)

when the integral exists. It can be seen from the Equation (2.3) that the independent variable is not changed as result of this transformation, so the output xˆ (t ) is also a time dependent function. It is normally not possible to calculate the Hilbert transform as an ordinary improper integral because of the possible singularity atτ = t . The integral is to be considered as a Cauchy principal value. In mathematics, the Cauchy principal value of certain is defined as b −ε

lim ⎡ ∫ ε →0+ ⎢ ⎣a

x(t )dt + ∫

c

b +ε

x(t )dt ⎤ , ⎥⎦

(2.4)

where b is a point at which the behaviour of the function x(t ) is such that

∫

b

a

and

∫

c

b

x(t )dt = ±∞

for any a < b

x(t )dt = m ∞ for any c > b . [7]

So when the Hilbert transform exists, it is written as presented at Equation (2.3). Other forms for H ( x(t )) can be obtained by change of variable, that is

17

H [ x(t )] =

H [ x(t )] =

π∫

∞

x(t − τ )

−∞

τ

1

∞

x(t + τ )

−∞

τ

1

π∫

dτ

(2.5)

dτ .

(2.6)

A Hilbert transform of a square wave is shown below:

Figure 2.2 The Hilbert transform of a square wave.

2.2.2 Frequency Response of the Hilbert Transform From the Equation (2.3), (2.4) and (2.5), it can be seen that Hilbert transform is a convolution: H [ x(t )] = xˆ (t ) =

1 * x(t ) πt

(2.7)

Equation (2.7) shows that xˆ (t ) is a linear function of x(t ) . It is obtained from x(t ) applying convolution with 1 /(πt ) .

18

According to the convolution theorem (the Fourier transform of a convolution of two functions is the product of their Fourier Transforms.), it can be seen that the spectrum of H [ x(t )] is related to that of x(t ) . Applied the Fourier transform to the Equation (2.7), that is F {xˆ (t )} =

⎧1⎫ F ⎨ ⎬ F{x(t )} π ⎩t ⎭ 1

(2.8)

where F is the Fourier transform. Since ∞ 1 ⎧1⎫ F ⎨ ⎬ = ∫ e − j 2π * fx dx = − jπ sgn( f ) ⎩ t ⎭ −∞ x

(2.9)

where

sgn f is +1 for f > 0, 0 for f = 0 and -1 for f < 0 . Rewriting Equation (2.9), that is: ⎧− j ⎛1⎞ F ⎜ ⎟ = − j sgn( f ) = ⎨ ⎝ πt ⎠ ⎩+ j Therefore, the Fourier transform of

for f > 0 . for f < 0

1 is − j sgn( f ) , which is equal to − j πt

for positive frequency and + j for negative frequency. Hence the Hilbert transform is equivalent to a kind of filter, in which the amplitudes of the spectral components are left unchanged, but their phases are altered by

π 2

, positively or negatively

according to the sign of frequency [8]. Therefore, the Fourier transform of the Hilbert transform of f (t ) is given by F {( xˆ )} = − j sgn fF {x(t )} .

19

(2.10)

The time domain result can be obtained performing an inverse Fourier transform.

2.3 Hilbert Transform Properties In this part, some properties of the Hilbert transform will be discussed.

(1) The Hilbert transform of a real function is linear. As discussed in the section 1.2, the Hilbert transform of a function f (t ) is defined as H [ f (t )] =

1

f (τ ) dτ −∞ t − τ

π∫

∞

.

Because of the possible singularity at t = τ , the integral is to be considered as a Cauchy principal value. It is expressed on the form H [ f (t )] = lim ε →0

1

π∫

x − t >ε

f (τ ) dτ . t −τ

(2.11)

Suppose f (t ) = c1 f1 (t ) + c 2 f 2 (t ) , then the Hilbert transform of f (t ) should be H [ f (t )] = H [c1 f1 (t ) + c 2 f 2 (t )]

= lim ε →0

1

π∫

= c1 lim ε →0

c1 f 1 (τ ) + c 2 f 2 (τ ) dτ x − t >ε t −τ

1

π∫

x − t >ε

f1 (τ ) 1 dτ + c 2 lim ε →0 π t −τ

= c1 H [ f1 (t )] + c 2 H [ f 2 (t )].

20

∫

x − t >ε

f 2 (τ ) dτ t −τ

(2.12)

Equation (2.12) shows the linearity of the Hilbert transform.

(2) The Hilbert transform of a Hilbert transform is the negative of the original function.

[

]

ˆ H [ fˆ (t )] = fˆ (t ) = F −1 − j sgn(ω ) F [ fˆ (t )]

= F −1 [− j sgn(ω )[− j sgn(ω ) F ( jω )]] = F −1 [− F ( jω )]

= − f (t )

(2.13)

(3) The Hilbert transform of the derivative of a function is equivalent to the derivative of the Hilbert transform of a function. [9] H [ f (t )] =

f (τ ) 1 dτ = ∫ − ∞ t −τ π π 1

∞

d {H [ f (t )]} = d ⎧⎨ 1 dt dt ⎩π

So

=

=

1

π∫

∞

1

∞

−∞

π∫

−∞

∫

∞

−∞

∫

∞

−∞

f (t − s ) ds s

f (t − s ) ⎫ ds ⎬ s ⎭

f ′(t − s ) ds s f ′(τ ) dτ t −τ

= H [ f ′(t )] . (4) A function and its Hilbert transform are orthogonal

∫

∞

−∞

∞ ⎡ 1 f (t ) fˆ (t )dt = ∫ f (t ) ⎢ −∞ ⎣ 2π

∫

∞

−∞

⎤ − j sgn(ω ) F ( jω )e jωt dω ⎥dt ⎦

21

(2.14)

=

1 2π

=

2 −j ∞ sgn( ω ) F ( j ω ) dω 2π ∫−∞

∞ − j sgn(ω ) F ( jω ) ⎡ ∫ f (t )e jωt dt ⎤ dω ⎢⎣ −∞ ⎥⎦ −∞

∫

∞

=0.

(2.15) 2

Since integrand sgn(ω ) F ( jω ) is an odd function which is integrated over

symmetric limits, the result is 0. Equation (2.15) proves that a real function and its Hilbert transform are orthogonal. This property can be used in energy and power signals.

(5) The energy in a real function and its Hilbert transform are equal. The signal and its Hilbert transform have identical energy because a phase shift does not change the energy of the signal only amplitude changes can do that. The energy of f (t ) and F (ω ) is defined as Ef = ∫

2

∞

f (t ) dt =

−∞

1 2π

∫

∞

−∞

2

F (ω ) dω .

(2.16)

So the energy of the Hilbert transform of f (t ) can be computed as E fˆ = ∫

2

∞

fˆ (t ) dt

−∞

=

1 2π

=

1 2π

∫

∞

−∞

∫

∞

−∞

2

− j sgn(ω ) F (ω ) dω 2

F (ω ) dω .

From Equation (2.16) and Equation (2.17), it shows that E fˆ = E f . 22

(2.17)

2.4 Hilbert Transform Applications The Hilbert transform is a very useful tool for the analysis of problems in various research areas. The Hilbert transform has a variety of applications, such as in the field of radio and signal processing, communication and power area.

2.4.1 Analytic Signal In digital signal processing, it is often needed to look at the relationships between the real part and imaginary part of a complex signal. The relationships are usually described by Hilbert transforms. Hilbert transform is also used to create special signals called analytic signals which are especially important in simulation. An analytic signal is a complex function created by taking a signal and then adding in quadrature its Hilbert transform [10]. An analytic signal is defined as z (t ) = f (t ) + jfˆ (t ) = A(t )e jθ ( t ) ,

(2.18)

where f (t ) is the input signal. fˆ (t ) is the Hilbert transform of the input signal.

z (t ) is the analytic signal constructed from f (t ) and its Hilbert transform. It is called the pre-envelope of f (t ) . The real and imaginary parts can be expressed in polar coordinates as: z (t ) = A(t )e jθ (t ) . where

23

(2.19)

A(t ) is the “envelope” or amplitude of the analytic signal given as A(t ) =

f 2 (t ) + fˆ 2 (t ) .

(2.20)

θ (t ) is the phase of the analytic signal given as ⎛ fˆ (t ) ⎞ ⎟. ⎟ f ( t ) ⎝ ⎠

θ (t ) = arctan⎜⎜

(2.21)

A(t ) and f (t ) have common tangents and the same values at the points where fˆ (t ) = 0 , i.e., the envelop determined using Equation (2.20) will have the same slope and magnitude of the original signal f (t ) at its maxima.

2.4.2 Analytic Signal Applied in Pattern Recognition [11] In this thesis, given a real function f (t ) , such as an ECG wave, it is possible to compute the Hilbert transform, fˆ (t ) . This allows the calculation of the envelope of f (t ) and also the phase of the pre-envelope of f (t ) and fˆ (t ) . If the two functions are then plotted in polar form (polar plot), the result is a waveform display very similar to a Vectorcardiogram (VCG) or a Polar-cardiogram (PCG). Thus the resulting magnitude versus angle plot is used for further analysis. The temporal dependence of the ECG data is removed. In effect the data has been shifted from a magnitude versus time system into a magnitude versus angle system. At the same time that because using a sampled-data waveform sampled at a fixed frequency (usually 360 Hz), the time information is still implicitly available to the

24

user. A major disadvantage of the time normalisation is that it implicitly assumes linear distortion of the ECG waveform over the length of the normalised segment. As mentioned previously, while the Hilbert transform display (polar plot) is reasonably familiar to a vectorcardiographic display, the data display is subject to different interpretation. Depending on the abnormality occurring in the ECG data, different displays are presented to the user. The main use of this display format is to monitor the data being placed in the pattern recognition so that different waveform segments (P,Q,R,S,T) and different time locations (before QRS ,after QRS) could be found. According to these information, the input ECG wave is normal or abnormal can be judged

25

Chapter 3 Computation of the Hilbert Transform In this chapter, four methods of implementing the Hilbert transform are given. It includes the Time-domain approach, the Frequency-Domain approach, the Boche approach and Remez filter approach.

3.1 Time-Domain Approach The Hilbert transform of a signal y (t ) at time t is given by yˆ (t ) =

1

y (τ ) dτ . −∞ t − τ

π∫

+∞

(3.1)

Assume that the signal y (t ) has been sampled every Δt second to give the sequence y k = y (kΔt ), k = 1,2,3,K, N and that the sampled Hilbert transform signal yˆ (k ) is to be computed. If the signal y (t ) is assumed to vary linearly during the sampling interval [12], for time from Δt to NΔt , the Hilbert transform at time Δt is

yˆ k = yˆ (kΔt ) = =

1

π∫

NΔt

Δt

1

π

(∫

2 Δt

Δt

y (τ ) dτ kΔt − τ kΔt ( k −1) Δt y (τ ) y (τ ) y (τ ) dτ + K + ∫ dτ + ∫ dτ ( k − 2 ) Δt kΔt − τ ( k −1) Δt kΔt − τ kΔt − τ

26

+∫

( k +1) Δt

1

k −2

kΔt

=

π

NΔt ( k + 2 ) Δt y (τ ) y (τ ) y (τ ) +∫ dτ + K + ∫ dτ ) ( N −1) Δt kΔt − τ kΔt − τ ( k +1) Δt kΔt − τ

(∑ I i( f ) + I k( c ) + i =1

N

∑I

i =k +2

(r ) i

),

(3.2)

where I i( f ) ≡ ∫

( i +1) Δt

iΔt

y (τ ) dτ kΔt − τ

( k +1) Δt y (τ ) y (τ ) dτ + ∫ dτ ( k −1) Δt kΔt − τ kΔt kΔt − τ

I k( c ) ≡ ∫

kΔt

I i( r ) ≡ ∫

iΔt

( i −1) Δt

y (τ ) dτ . kΔt − τ

When y (t ) is linear during the sampling period, y (t ) = y i + ( y i +1 − y i )(t − iΔt ) / Δt for iΔt ≤ t ≤ (i + 1)Δt y (t ) = y i + ( y i −1 − y i )(t − iΔt ) /(− Δt ) for (i − 1)Δt ≤ t ≤ iΔt . So I i( f ) ≡ ∫

( i +1) Δt

iΔt

=∫

( i +1) Δt

=∫

Δt

iΔt

0

=∫

Δt

0

y (τ ) dτ kΔt − τ y i + ( y i +1 − y i )(τ − iΔt ) / Δt dτ kΔt − τ

y i + ( y i +1 − y i ) x / Δt dx ( x ≡ τ − iΔt ) kΔt − iΔt − x

y i + ( y i +1 − y i )τ / Δt dτ (k − i )Δt − τ

(3.3)

and y (τ ) dτ ( i −1) Δt kΔt − τ

I i( r ) ≡ ∫

iΔt

27

=∫

iΔt

( i −1) Δt

y i + ( y i −1 − y i )(τ − iΔt ) /(− Δt ) dτ kΔt − τ

y i + ( y i −1 − y i ) x / Δt dx ( x ≡ −(τ − iΔt )) kΔt − iΔt + x

= −∫

0

= −∫

Δt

Δt

0

y i + ( yi −1 − y i )τ / Δt dτ (i − k )Δt − τ

(3.4)

and ( k +1) Δt y (τ ) y (τ ) dτ + ∫ dτ ( k −1) Δt kΔt − τ kΔt kΔt − τ

I k( c ) ≡ ∫

kΔt

( k +1) Δt y + ( y y k + ( y k −1 − y k )(τ − kΔt ) /(− Δt ) k k +1 − y k )(τ − kΔt ) / Δt dτ + ∫ dτ ( k −1) Δt kΔt kΔt − τ kΔt − τ Δt y + ( y 0 y + (y − y k k k −1 ) x / Δt k +1 − y k )ω / Δt =∫ dx + ∫ k dw − Δt 0 −w −x

=∫

kΔt

( x ≡ τ − kΔt , w ≡ τ − kΔt ) =∫

0

− Δt

Δt y + ( y y k + ( y k − y k −1 )τ / Δt k +1 − y k )τ / Δt dτ + ∫ k dτ . 0 −τ −τ

(3.5)

Further calculation gives I i( f ) = ∫

y i + ( y i +1 − y i )τ / Δt dτ (k − i )Δt − τ

Δt

0

= −∫

Δt

0

y i + ( y i +1 − y i )(τ − (k − i )Δt ) / Δt + ( y i +1 − y i )(k − i ) dτ τ − (k − i)Δt

Δt

Δt

0

0

= − ∫ ( y i +1 − y i ) / Δtdτ − ∫

yi + ( y i +1 − y i )(k − i ) dτ τ − (k − i)Δt Δt

= −( y i +1 − y i ) − ( y i + ( y i +1 − y i )(k − i )) ln τ − (k − i )Δt 0 = y i ln

k −i k −i + ( y i +1 − y i )(−1 + (k − i ) ln ) k − i −1 k − i −1

and

28

(3.6)

I i( r ) = − ∫

y i + ( y i −1 − y i )τ / Δt dτ (i − k )Δt − τ

Δt

0

=∫

y i + ( yi −1 − y i )(τ − (i − k )Δt ) / Δt + ( yi −1 − y i )(i − k ) dτ τ − (i − k )Δt

Δt

0

Δt

Δt

0

0

= ∫ ( y i −1 − y i ) / Δtdτ + ∫

y i + ( y i −1 − yi )(i − k ) dτ τ − (i − k )Δt Δt

= ( y i −1 − y i ) + ( y i + ( y i −1 − y i )(i − k )) ln τ − (i − k )Δt 0 i − k −1 i − k −1 + ( y i −1 − y i )(1 + (i − k ) ln ) i−k i−k

= y i ln

(3.7)

and I k( c ) = ∫

Δt y + ( y y k + ( y k − y k −1 )τ / Δt k +1 − y k )τ / Δt dτ + ∫ k dτ 0 −τ −τ

0

− Δt

= −∫

0

Δt

=∫

Δt

Δt y + ( y y k + ( y k − y k −1 )(− x) / Δt k +1 − y k )τ / Δt dx − ∫ k dτ 0 x τ

y k + ( y k − y k −1 )(−τ ) / Δt

τ

0

dτ − ∫

Δt

0

y k + ( y k +1 − y k )τ / Δt

τ

dτ

Δt

= ∫ ( y k −1 − y k +1 ) / Δtdτ . 0

= y k −1 − y k +1 .

(3.8)

So the Hilbert transform of y k is yˆ k =

1

π =

k −2

(∑ I i =1

1

π

k −2

(f) i

+I

(∑ ( y i ln i −1

(c) k

+

N

∑I

i =k + 2

(r ) i

)

k −i k −i + ( y i +1 − y i )(−1 + (k − i ) ln )) k − i −1 k − i −1

+ ( y k −1 − y k +1 )

29

+

N

∑ (y

i =k +2

i

ln

i − k −1 i − k −1 + ( y i −1 − y i )(1 + (i − k ) ln ))) . i−k i−k

(3.9)

The results given on the reference [5] is Ii

(f)

Ik

Ii

(c)

(r )

= y i ln

k −i k −i + ( y i +1 − y i )(−1 + (k − i ) ln ) k − i −1 k − i −1

= y k −1 − y k +1

= y i ln

i−k i−k + ( y i −1 − y i )(−1 + (i − k ) ln ). i − k −1 i − k −1

Noticed that I i same. But I i

(r )

(f)

and Equation (3.6), I k

(c)

and Equation (3.8) are the exactly

is different from Equation (3.7). The equation does influence the

result of the Hilbert transform of a signal. So it is important to prove which one (Equation (3.10) or the one given on the reference [12]) is correct, This can be proved from the example below. Here an example is given. Consider the signal f ( x) = sin( x) . From Equation (3.1), its Hilbert transform is H [sin( x)] =

1

sin(τ ) dτ . −∞ x − τ

π∫

∞

Letting s = x − τ , get H (sin x) =

1

sin( x − s ) (−ds ) −∞ s

π∫

=−

=−

∞

1

sin x cos s − cos x sin s ds −∞ s

π∫

∞

∞ cos x sin s 1 ⎡ ∞ sin x cos s ⎤ ds − ∫ ds ⎥ ∫ ⎢ − ∞ − ∞ π⎣ s s ⎦

30

(3.10)

=−

It is well known that

So

1

∞ cos s ∞ sin s 1⎡ ⎤ − sin x ds cos x ds ⎥ ∫ ∫ ⎢ − ∞ − ∞ π⎣ s s ⎦

cos s 1 ds = 0 and −∞ π s

π∫

∞

H (sin x) = cos x

1

π∫

∞

−∞

sin s ds = 1 . −∞ s

∫

∞

sin s ds s

= cos x .

(3.11)

In order to test that the Equation (3.9) is right, here an example is given. Assuming the input function is y = sin(2π * f * N ) , letting f = 0.02 Hz, N = 500 . From the Equation (3.11), it can be seen that the Hilbert transform of y = sin(2π * 0.02 * 500) should be cos(2π * 0.02 * 500) . Writing a MATLAB program for this algorithm, Figure 3.1 and Figure 3.2 were obtained. Figure 3.1 shows the input waveform and Figure 3.2 (a) illustrates the Hilbert transform waveform of this input wave using Equation (3.9). Figure 3.2 (b) shows the results using the equation from reference [12].

31

Figure 3.1 Input waveform: sin(2π * 0.02 * 500) .

Figure 3.2 (a) Output waveform: The Hilbert transform of sin(2π * 0.02 * 500) .

32

Figure 3.2 (b) Output waveform using the equation from [12].

As shown in Figure 3.2 (a), the output waveform is the Hilbert transform of the input sine wave. It is a cosine wave. It also illustrates that the Hilbert transform of a real function does not change the amplitude of the signal but only changes its phase by

π 2

rad/s.

The waveform shown in Figure 3.2 (b) is not a cosine wave. From the waveform it also can be seen that I i

(r )

from [12] is not correct and Equation (3. 9)

is correct.

As discussed in previous chapters, the input wave in the Holter monitoring system is an ECG wave. Figure 3.3 is the wave obtained from the MIT-BIH

33

(Massachusetts Institute of Technology-Beth Israel Hospital) arrhythmia database [see Appendix A] MIT213. According to the Time-Domain Approach, using Equation (3.9), the Hilbert transform of this ECG wave can be computed. The input ECG waveform is shown in Figure 3.3. According to the explanation about the ECG waves in Chapter 1, from the value of the P wave, the QRS complex, PR interval, QRS duration, RR interval and PP interval, it can be seen that Figure 3.3 is a normal ECG wave. Because each R wave stands for a beat, Figure 3.3 describes 4 beats of a heart.

Figure 3.3 A normal ECG wave.

34

Figure 3.4 The Hilbert transform of the ECG wave.

Figure 3.4 shows the results of the Hilbert transform of this ECG wave. The Hilbert transform waveform of the ECG wave should oscillate from negative to positive or from positive to negative around the X-axis. The points corresponding to peak values of R wave should be zero in the output waveform. But from Figure 3.4, it can be seen that the output waveform is distorted. It’s not the correct waveform, so this method of computing the Hilbert transform may not be suitable for the ECG wave.

35

3.2 Frequency-Domain Approach To further investigate the Hilbert transform, the frequency domain analysis is very useful. The second method to compute the Hilbert transform of a function is Frequency-Domain approach. As shown before, the Hilbert transform of the function y (t ) is: 1 ) y (t ) =

π∫

∞

−∞

y (τ ) dτ . t −τ

Because the usual time domain definition based on the Cauchy principal value of an integral is usually not easy to calculate, the Hilbert transform in the frequency ) domain is defined. Suppose Y ( f ) and Y ( f ) are the Fourier transform of y (t ) and ) yˆ (t ) . Y ( f ) and Y ( f ) are defined as ∞

Y ( f ) = ∫ y (t )e − j 2πft dt −∞

) Y ( f ) = − j sgn( f )Y ( f ) .

(3.12)

Applying the Fourier transform to the convolution defined in the equation above, can obtain

) ) ⎡1⎤ Y ( jω ) = FT [ y (t )] = FT ⎢ ⎥Y ( jω ) = − j sgn(ω )Y ( jω ) [5], ⎣ πt ⎦

(3.13)

where FT [] represents Fourier transform. This equation indicates that the Hilbert transform can be interpreted in the frequency domain.

36

) So given a sampled signal y k , the sequence y k can then be computed using fast Fourier transform (FFT) techniques as

) y k = IFFT [− j sgn(ω n ) FFT [ y k ]] ,

(3.14)

where FFT [] represents the fast Fourier Transform operation. IFFT [] represents inverse fast Fourier transform operations.

ω n represents the n th frequency of the discrete Fourier transform. sgn is the sign function. This formula can be used to calculate the Hilbert transform, by first taking the Fourier transform of y k , multiplying it by − j sgn(ω n ) , then taking the inverse

π ) Fourier transform, thus obtaining y k . Thus the Hilbert transform is a − phase 2 ) shifter when viewed as a linear system whose input is y k and output is y k . Here an example is also given. The input is the sine wave, sin(2π * 0.02 * 500) , used previously, and the output should be the Hilbert transform of the input, i.e., a cosine wave. Figure 3.5 shows the input waveform sin(2π * 0.02 * 500) and Figure 3.6 shows the output results using FrequencyDomain approach to compute the Hilbert transform of a function.

37

Figure 3.5 Input wave: sin(2π * 0.02 * 500) .

Figure 3.6 The Hilbert transform of sin(2π * 0.02 * 500) .

38

Figure 3.7 shows the same section of ECG wave taken from the MIT-BIH database, MIT213, used previously (Figure 3.3). Using the Frequency-Domain approach to obtain the Hilbert transform, the output waveform is shown in Figure 3.8.

Figure 3.7 Input ECG wave

39

Figure 3.8 Output: the Hilbert transform of the ECG wave.

The Hilbert transform waveform shown in Figure 3.8 is good. The wave oscillates around the X-axis. The zero crossing points corresponding to the R peak wave are right. From the examples given above, the figures show that the Hilbert transform of the sine wave is cosine wave. The output wave oscillates around zero as it should when the input wave is the ECG wave. Even though this algorithm works, it’s also not suitable for the system developed in this thesis. The reason will be given later.

40

3.3 Boche Approach The Boche approach [13] presents a new algorithm by reconstructing a bandlimited function from samples to calculate the Hilbert transform. The algorithm can be described as follows: A set of discrete instants {t i } are given with the corresponding values { y i } where exists a function f (t i ) = y i . A statement can be made for approximating function as Equation (3.15): n

f n (t ) = ∑ bk ,n k =1

n

f n (t i ) = ∑ bk ,n k =1

sin π (t − t k ) π (t − t k )

(3.15)

sin π (t i − t k ) . π (t i − t k )

(3.16)

To calculate the coefficients bk ,n of the Equation (3.16), a system of n linear equations, Equation (3.17), have to be solved: ⎡ f n (t1 ) ⎤ ⎡ a11 ⎢ M ⎥=⎢ M ⎢ ⎥ ⎢ ⎢⎣ f n (t n )⎥⎦ ⎢⎣a n1

where aik =

sin π (t i − t k ) π (t i − t k )

L aik L

a1n ⎤ ⎡b1n ⎤ M ⎥⎥ ⎢⎢ M ⎥⎥ , a nn ⎥⎦ ⎢⎣bnn ⎥⎦

(3.17)

t i ≠ t k for i ≠ k .

Using the iteration method to solve the linear equations, the coefficients bk ,n were obtained. The Hilbert transform of the Equation (3.17) can be derived as follows: It is well known that sin(π (t − t k )) dt = 1 −∞ π (t − t k )

∫

∞

41

and cos(π (t − t k )) dt = 0 . −∞ π (t − t k )

∫

∞

According to the definition of the Hilbert transform of a real function, the sin(π (t − t k )) Hilbert transform fˆ (t ) of f (t ) = is derived as follows: π (t − t k ) 1 fˆ (t ) =

π∫

=

1

∞

−∞

sin(π (t − t k )) 1 dt π (t − t k ) λ − t

sin(π (t − t k )) 1 1 ( + )dt −∞ π (λ − t k ) λ − t t − t k

π∫

∞

=

∞ sin(π (t − t )) sin(π (t − λ + λ − t k )) 1 k ( − )dt ∫ π (λ − t k ) −∞ π (t − t k ) π (t − λ )

=

∞ sin(π (t − t )) sin(π (t − λ )) cos(π (λ − t k )) + cos(π (t − λ )) sin(π (λ − t k )) 1 k ( − )dt ∫ − ∞ π (λ − t k ) π (t − t k ) π (t − λ )

=

1 (1 − cos(π (λ − t k ))) π (λ − t k )

=

1 − cos(π (λ − t k )) . π (λ − t k )

(3.18)

It is shown that the Hilbert transform of the series n

f n (t ) = ∑ bk ,n k =1

sin π (t − t k ) π (t − t k )

is n

H { f n (t )} = ∑ bk ,n k =1

1 − cos(π (t − t k )) , π (t − t k )

where bk ,n are the coefficients in Equation (3.17).

42

(3.19)

Here an example is given to demonstrate the algorithm for the sampled signal as well as its Hilbert transform. Assuming the sampled function g is given by f (t ) =

sin(2t ) 3 sin(3(t − 5)) 1 sin( 2.5(t + 6)) − + . 2t 4 3(t − 5) 2 2.5(t + 6)

The samples were taken in the interval − 10 ≤ t ≤ 10 with a sampling interval of 0.25, thus yielding 81 sampling points. Based on the Equation (3. 19), i.e., 1 − cos(π (λ − t k )) fˆ (λ ) = π (λ − t k )

and the linear property of the Hilbert transform that discussed in Chapter 2, the Hilbert transform, fˆ (t ) , of the function , f (t ) , can be written as follows : 1 − cos(2t ) 3 (1 − cos(3(t − 5))) 1 (1 − cos(2.5(t + 6))) fˆ (t ) = − + . 2t 4 3(t − 5) 2 2.5(t + 6) Running a MATLAB program, the results are shown in Figure 3.9 and Figure 3.10. Figure 3.9 (a) shows the waveform of original function f (t ) , Figure3.9 (b) is the 26th approximation f 26 (t ) after 26 iterations. The error function is obtained using f 26 (t ) minus f (t ) , and the result is shown in Figure 3.9 (c). The error function is smaller than 0.0001 shown in Figure 3.9 (c). The sample sequence t i is given in Figure 3.9 (d).

43

Figure 3.9 (a) f (t )

(b) f 26 (c) f 26 (t ) − f (t ) (d) t i .

Figure 3.10 (a) shows the Hilbert transform fˆ (t ) of the original function f (t ) , and the 26th approximation fˆ26 (t ) is shown in Figure 3.10 (b), and then the error function of the Hilbert transform is illustrated in Figure 3.10 (c).

44

Figure 3.10 (a) fˆ (t ) (b) fˆ26 (t ) (c) fˆ26 (t ) − fˆ (t )

Boche algorithm permits reconstruct the bandlimited function from samples and recovery of the Hilbert transform of this function. Compared with other known solutions for computing Hilbert transform of a function, this algorithm does not need to calculate integrals. However a set of linear equations has to be solved in each iteration step. Here the algorithm is just given as a reference. In this thesis, the difficulty of computing the solution to a variable set of linear equations is not easier than calculating an integral. So this method is not suitable for the Holter monitoring system in this thesis.

45

3.4 Remez Approach This section mainly talks about using “ remez ” function in MATLAB to obtain the Hilbert transform of a signal. The Remez Exchange FIR filter design approach (also called the ParksMcClellan or Optimal method) is a popular technique used to design FIR filters. The well known Parks-McClellan algorithm uses this approach and Chebyshev approximation theory to generate filters with an optimal fit between the desired frequency responses and actual frequency responses. Filters designed in this way illustrate an equiripple wave in the frequency response. By implementing the ParksMcClellan algorithm [see Appendix B], the Remez approach designs a linear-phase FIR filter. The syntax of the remez function can be written the following way: b = remez(n, f , a) b = remez(n, f , a) returns a row vector b including the n + 1 coefficients of

the order n FIR filter. where “ n ” represents the order of the filter. “ f ” represents a vector of pairs of normalized frequency points. The frequencies are specified in the range between 0 and 1, where 1 corresponds to the Nyquist frequency. The frequencies must be in increasing order. “ a ” represents a vector containing the desired amplitudes at the points specified in f . f and a are the same length. The length is an even number.

46

The

output

coefficients

in

b

satisfied

the

symmetry

relation

b(k ) = b(n + 2 − k ), k = 1,..., n + 1 ;

Remez function can specify the different filter type: b = remez(n, f , a, ' ftype' ) . b = remez(n, f , a, w, ' ftype' ) is used when the special filter type is needed.

where ' ftype'

represents the filter type parameter. It includes three types:

Multiband, Differentiator and Hilbert transform. The one used in this thesis is 'hilbert', that is the Hilbert transformer, for linear-phase filters with odd symmetry. The output filter coefficients in b satisfies b(k ) = −b(n + 2 − k ) , k = 1,2,K , n + 1 . The Hilbert transformer has the desired amplitude of 1 across the entire band.

Here an example is given: h = remez (100, [0.05 0.95], [1 1], ' hilbert ')

designs an approximate FIR Hilbert transformer of length 100. The frequencies are specified from 0.05 to 0.95 and their corresponding amplitudes are 1. The amplitude will be 0 at other frequencies.

47

Figure 3.11 The input wave: y = sin(2π * 0.02 * 500) .

Figure 3.12 The output wave: HT of y = sin(2π * 0.02 * 500) .

48

As before, the input signal is the same sine wave y = sin(2π * 0.02 * 500) , used previously. The remez function discussed above was used to generate a Hilbert transformer and obtain the Hilbert transform of the input sine signal. The input waveform is shown in Figure 3.11. Figure

3.12

shows

the

Hilbert

transform

of

the

input

signal

y = sin(2π * 0.02 * 500) . Because the order of the filter used here is M = 100 , the

filter phase delay should be 100 / 2 = 50 . From the output waveform, it can be seen clearly that the filter has a phase delay for N from 0 to 50. The output results corresponding to the input signal should be calculated from N = 51 . Figure 3.12 illustrates that the Hilbert transform of a sine function is a cosine function. Another example is given here. The input wave is still the normal ECG wave from the MIT-BITH Arrhythmia database MIT213 used previously (see Figure 3.3).

49

Figure 3.13 An ECG wave.

Figure 3.14 The Hilbert transform of an ECG wave.

50

Figure 3.14 shows the output, i.e., the Hilbert transform of the input ECG signal. For the same reason, the filter phase delay, the output wave lags by 50 samples.

3.5 Comparison Four methods for computing Hilbert transform are discussed above. From this section, the comparison of these four methods is given. Table 3.1 shows the different running time and the easy level of hardware (HW) and software (SW) implementation. All of programs are MATLAB programs and timing was done using the tic/toc functions. Table 3.1 The comparison of the four methods for computing Hilbert transform. Time Domain

FrequencyDomain

Approach

Approach

Sine

ECG

(N=500) (N=800) Running

2.7740s

7.1110s

Sine

ECG

Remez Filter Boche Approach

0.010s

Sine

ECG

(N=500) (N=800)

(N=500) (N=800) 0.010s

Approach

0.6110

0.0140s

0.030s

Time

HW/SW

HW&SW

HW

HW&SW

Problem

Problem

Problem

OK

For Time-Domain approach, even though the equations for computing the Hilbert transform have been derived and it is not needed to calculate the integral anymore, it still needs to compute the “ ln ” function. It is hard to implement in the Nios system that will be used in this thesis because there is no hardware “ ln ”

51

function in the Nios processor. The running time was measured for sine wave and ECG wave individually. Note that the sample number is 500 for sine wave and 800 for ECG wave. The results are shown in the Table 3.1. Using this method, the output waveform that corresponding to an input ECG signal is distorted to some extent. For Frequency-Domain approach, the running time is very fast, but it is inconvenient to be implemented on the hardware because it needs to compute the FFT and IFFT for this method. This is not that easy to implement on the hardware. The Boche approach supplied a simple way to calculate the Hilbert transform of the bandlimited function. Even though it does not need to calculate the integral, a set of linear equations has to be solved. The size of the equations is variable. It’s not that easy to implement in the Nios processor used in this thesis. Compared with the first three methods discussed above, the Remez filter approach is a better choice to compute Hilbert transform of a function in this thesis. Its running time is shorter and can be implemented on both hardware and software.

The computer performance that used in this thesis: Operation System: Microsoft Window XP CPU: Intel® Pentium®4 2.66GHz RAM: 512MB.

52

Chapter 4 Implementation

In the previous chapters, four methods to compute the Hilbert transform applied to the input data have been discussed. After comparing them in terms of running time and the complexities of the software and hardware implementation, the Remez approach was selected to be the best method for the ECG Holter Monitoring System. The problem now, is how to implement the Hilbert transform algorithm to build the Holter monitoring system on the Nios system. In the following sections, this will be discussed in detail.

4.1 Nios Embedded Processor Overview The Nios embedded processor is a user-configurable, 16-bit ISA (Instruction Set Architecture), general-purpose RISC (Reduced Instruction Set Computer) embedded processor that was designed to be a very flexible and powerful processor solution [14]. The Nios embedded processor has become a commonly used embedded processor because of its ease-of-use and flexibility. The Nios embedded processor system structure is shown in Figure 4.1.

53

The Quartus II software, the SOPC (System on a Programmable Chip) Builder system development tool, is used to build and evaluate custom processorbased systems.

Designers can use SOPC Builder to integrate one or more

configurable Nios CPUs with any number of standard peripherals, gluing the system together [14]. Using SOPC Builder, a user can combine the Nios processor with user logic and program it into a FPGA (Field Programmable Gate Array) easily.

C/C++ VHDL/Verilog

System Component

Develop Design SOPC Builder [email protected] II Software Embedded Software Development Tools Third-Party Tools

Development Boards & Kits

CPU

Nios Development Kit, Stratix Edition

Avalon Switch Fabric

Nios Development Kit, Cyclone Edition

Peripherals

Nios Development Kit Stratix Professional Edition

On-Chip Debugging

Third-Party Boards&Kits

Figure 4.1 Nios Embedded Processor System.

54

In this thesis, the development tool is SOPC Builder, using a Nios development kit, Stratix Professional Edition. It is a complete embedded systems development kit for the Nios embedded processor. There are a number of necessary steps to create a Nios system on the Nios development board. The procedure is shown in Figure 4.2. The flow includes both the hardware and software design tasks required to create a working system. The right side illustrates the software development flow and the left side illustrates the hardware design flow. Based on the system requirements, the hardware design begins with the SOPC builder system integration software. At this point, the designer can begin writing device-independent C/C++ software. After the hardware designer defines the customer Nios processor hardware system using SOPC Builder, SOPC Builder generates a custom software development kit (SDK) that forms the foundation for the software development flow [16]. With the SDK, the designer can begin writing software that interacts at the low level with hardware components. The Nios SDK Shell provides an UNIX bash shell environment board on a PC platform. It is a very useful utility. Figure 4.3 shows the Nios SDK shell (bash environment).

55

Step1: Predesign Activity z Analyze System Requirements (Performance &Throughout) z Defne Nios Processor Subsystem (CPU, Peripherals, Memory Structure, etc.)

Software Development

Hardware Development Standard System Components

Step2: Begin C/C++ Development

Define Nios Processor System Module with SOPC Builder

Custom SDK

Step3: Develop Drivers& Routines For Custom Hardware User-Defined Components

Software Libraries

OS Kernel

Assign Device, Layout Pins &Compile Hardware with the Quartus II Softwar Hardware Prolotype on development Board

Create Custom Acceleration Hardware

No

Step5: Software Prototype on Development Board

Does System

Step4: Compile & Link, Targeting Custom Hardware Platform

No

Meet Goals?

Yes Step6 Successful Prototype of Nios System Module

Figure 4.2 Hardware/Software development flow for a Nios processor system [15]

56

Drivers &Routines

Figure 4.3 Nios SDK Shell (bash).

This “bash” environment can be used for all related development work for the Nios system and communicate with the Nios development board. The Nios development board is shown in Figure 4.4. The Nios development kit includes many Nios-specific utilities that can run in the Nios SDK Shell to generate and debug software. The Nios SDK Shell also can be used to run test programs on the Nios development board.

57

Figure 4.4 Nios Development Board Components [17]

58

Figure 4.4 shows the Nios development board components. It includes [18] • Stratix EP1S40F780 device • MAX EPM7128AE CPLD configuration control logic • SRAM (1 Mbyte in two banks of 512 Kbytes, 16-bit wide) • SDR SDRAM (16 Mbytes, 32-bit wide) • Flash (8 Mbytes) • CompactFlash connector header for Type 1 CompactFlash cards • 10/100 Ethernet physical layer/media access control (PHY/MAC) • Ethernet connector (RJ-45) • Two serial connectors (RS-232 DB9 port) • Two 5-V-tolerant expansion/prototype headers • Two JTAG connectors • 50-MHz crystal (socket), external clock input • Mictor connector for debugging • Four user-defined push-button switches • Eight user-defined LEDs • Dual 7-segment LED display • Power-on reset circuitry.

Hardware designers can use the Nios development board as a platform to prototype complex embedded systems. Software developers can use the Nios reference design pre-programmed on the development board to begin prototyping software immediately.

4.2 Digital Filter In digital signal processing, an important function of a filter is to remove unwanted parts of the signal, such as random noise, or to extract useful parts of the signal, such as the components lying within a certain frequency range [19]. Figure 4.5 illustrates the basic concept.

59

raw (unfiltered) signal

filtered signal FILTER Figure 4.5 Filter.

The filter function is implemented as a direct form II transposed structure as shown in Figure 4.6.

Figure 4.6 Digital filter transposed structure.

60

For a linear time-invariant system (Figure 4.6), its input and output satisfy the following equation N

M

k =1

k =0

y[n] − ∑ a k y[n − k ] = ∑ bk x[n − k ]

(4.1)

with the corresponding rational system function M

H ( z) =

∑b z k =0 N

−k

k

1 − ∑ ak z

[20],

(4.2)

−k

k =1

where: x[n − k ] is the previous input. y[n] is the output. y[n − k ] is the previous output. a k and bk are the filter coefficients. H ( z ) is the filter’s Z transform.

For causal FIR (Finite Impulse Response) system, the system function has only zeros (except for poles at z=0), FIR filter does not depend on the past values of the output. FIR filters are therefore non-recursive. Since the coefficients a k are all zero, the Equation (4.1) reduces to M

y[n] = ∑ bk x[n − k ] .

(4.3)

k =0

When the filter sequence (impulse response) of FIR filter is either symmetric or anti-symmetric, the filter is of linear phase. Such filters do not distort the phase

61

of the input signal. It is well known that FIR filters can always be designed such that they exhibit the desirable characteristic. If the impulse response of the FIR filter satisfies the condition h[ M − n] = h[n]

for n = 0,1,K , M ,

it is called symmetric. If the impulse response of the FIR filter satisfies the condition h[ M − n] = − h[n]

for n = 0,1, K, M ,

it is called anti-symmetric. Table 4.1 categorizes linear phase filters according to their symmetry and length. Table 4.1 Four types of the linear phase FIR filter Type

Impulse Response

1

symmetric

Length (M+1) is odd

2

symmetric

Length (M+1) is even

3

anti-symmetric

Length (M+1) is odd

4

anti-symmetric

Length (M+1) is even

Examples of the four types of impulse response sequences are shown in Figure 4.7.

Figure 4.7 Illustration of four types of impulse response symmetry.

62

A Type 1 filter may be used to implement any desired bandpass frequency response. A Type 2 filter may not be used to define a highpass filter since the symmetry condition requires H (π ) = 0 . It can be used instead of Type 1 in cases where an even length filter is preferable [21]. Antisymmetric filters can be used to design FIR differentiators and Hilbert transformers. Differentiators are antisymmetric FIR filters with approximately linear magnitude responses. Hilbert transformers are anti-symmetric FIR filters with approximately constant magnitude.

4.3 Implementation In this section, the implementation for performing the Hilbert transform of an input signal will be discussed in detail.

4.3.1 Filter Order As discussed in previous chapter, the order of the filter should be determined to meet certain filter specifications including the passband ripple, stopband attenuation and the transition bandwidth. MATLAB is a perfect tool for this purpose. A sine wave signal was used as the test signal. Let M represent the filter order. The Hilbert transform experiments were conducted with the MATLAB program and 501 samples of the test signal sin(2 * π * 0.02 * n) with n = 0, K,500 . The input signal is shown in Figure 4.8. The Hilbert transform outputs are shown in Figures 4.9-4.13 for the filter order M = 51,71,91,101 and 201. 63

Figure 4.8 Input waveform: sin(2 * π * 0.02 * 500) .

Figure 4.9 Output waveform, M = 51 .

64

Figure 4.10 Output waveform, M = 71 .

Figure 4.11 Output waveform, M = 91 .

65

Figure 4.12 Output waveform, M = 101 .

Figure 4.13 Output waveform, M = 201 .

66

The group delay caused by the FIR filter with the anti-symmetric coefficients and order

M

is

M . 2

As result, the phase delay when the input

sin(2 * π * 0.02 * n)

signal

passes

through

the

filter

is

M π * 0.02 * 2π + = (0.02 * M + 0.5)π . 2 2

The output maximum and the minimum values were checked. The amplitude error and the phase delay for different order filters were also illustrated in Table 4.2. From the Table 4.2, it can be seen that the higher the order is, the less the error is when M is changed from 51 to 101. In order to not take long time for compute the Hilbert transform of the ECG wave, the filter order should be chosen properly.

Table 4.2 Filter order comparison M

Maximum

Minimum

Error

Phase Delay

51

0.934020

-0.933397

<6.63%

1.52π

71

0.975858

-0.967608

<2.83%

1.92π

91

1.000212

-0.983314

<0.8%

2.32π

101

1.008626

-0.987704

<0.2%

2.52π

201

1.046485

-0.997695

<2.44%

4.52π

The frequency responses are shown in Figure 4.14 and Figure 4.15 when the order of the filter is 100 or 101.

67

Figure 4.14 Frequency response, M = 100 .

Figure 4.15 Frequency response, M = 101 .

68

Based on the delay analysis above, the phase response is − 360 *

M f * − 90 2 2

degree for the normalized sampling frequency of 2. The slopes of the phase responses in Figure 4.14 and 4.15 are -9090 degree/Hz and -9180 degree/Hz, respectively. From the Figure 4.9, Figure 4.10, Figure 4.11 and Figure 4.12 and Figure 4.13 it can be seen that the order of the filter must be bigger than 91 so that the amplitude of the output wave is ≈ 1 . The phase shift depends on the filter order as shown above. In addition, in this thesis, the error requirement is to be smaller than 0.5%. From the Figure 4.14 and Figure 4.15, when the order is odd, the frequency response is not symmetrical (see the right hand side of Figure 4.14 and Figure 4.15). So the order of the filter in this thesis is determined as 100.

4.3.2 Filter Coefficients 1. When using MATLAB, the remez function: b = remez(n, f , a, ' filtertype' )

can be used to get the coefficients of the filter directly. That is

b = remez(100, [0.05 0.95], [1 1], ' h' ) , where n = 100 is the filter order.

f = [0.05 0.95] is a vector of pairs of normalized frequency point, specified in the range 0 to 1, where 1 corresponds to the Nyquist frequency. This frequency range is determined according to a few tests.

69

a = [1 1] is a vector containing the desired amplitudes 1 at the points specified

in f = [0.05 0.95]. It has the same length as f . ' filtertype' =' h' specifies that the filter is Hilbert transformer. This parameter allows specifying one of the following filters: Multiband, Differentiator, and Hilbert transformer. Running a MATLAB program, Figure 4.16 shows that the calculated remez filter is an equi-ripple bandpass filter with a symmetrical magnitude response around f = 0.5 . The coefficient b was obtained as shown in Figure 4.17.

Figure 4.16 Frequency response of the ideal and Remez design filter.

70

Figure 4.17 Coefficients of the filter when M = 100 .

71

From the output coefficients shown in the Figure 4.17, it can be seen that the Hilbert transformer has negative symmetry.

Figure 4.18 Frequency response of Remez filter.

Figure 4.18 demonstrates that the calculated filter has a symmetrical magnitude response around f = 0.5 and a linear phase response with the slope of 9090 degree/Hz.

2. Using C program remez.c to calculate the coefficients of the filter This program uses the Remez exchange algorithm to design linear phase FIR digital filters with minimum weighted Chebyshev error in approximating a desired

72

ideal frequency response [22]. The program has a special built-in section for specifying the more common ideal filter types such as multi-band, bandpass filters, Hilbert transform filters, and differentiators [23]. remez.c calculates the optimal FIR filter impulse response for a set of given

band edges, the desired response and the weight on those bands. It includes a main program that handles the input, sets up the appropriate approximation problem and handles the output of the optimal filter coefficients. Function remez (h[], numtaps, numband , bands, desired [], weight[], type)

has input values numtaps , numband , bands , desired[] , weight[] , type and output value h[] , where

h[] is the impulse response of the filter, i.e. the coefficients of the filter. numtaps is an integer. Specifying number of the filter coefficients. It should

be M + 1 , M is the order of the filter. numband is an integer, specifying number of bands in filter specification. bands is a double variable, specifying user_specified band edges, using upper

and lower cutoff frequencies. The bands array specifys the set F to be the form F = UBi where each frequency band Bi is a closed subinterval of the frequency

axis [0, 1/2]. The number of bands should be 2 * numband .

desired[] is an array, which is the user_specified band responses, the desired frequency response in each band. The number of desired[] should equal to the number of bands . 73

weight[] is an array, which is the user_specified error weights, a positive

weight function in each band. The number of weight[] equals to the number of bands. The array desired [] and weight[] specify the ideal response and weight function in each band. type is the type of the filter. It includes:

(a) Multi-band filter; (b) Bandpass filters; (c) Hilbert transform filters; (d) Differentiators. In this thesis, it has been explained that the order of the filter is M = 100 , so the numtaps should be M + 1 = 101 . The numband is 3. The bands , the desired [] and the weight[] is defined as shown in Figure 4.19. The type parameter is ' HILBERT ' . The output result h[] is saved in a text file “coefile.txt”. To test program for the remez() function, the appropriate arguments to remez () was used to generate a filter. The initial parameters are shown in Figure

4.19. The resulting coefficients are shown as the Figure 4.20.

74

Figure 4.19 Initial parameter of remez.c .

A C program was run in Visual C++ and the filter coefficients are obtained. They are shown in Figure 4.20. The frequency response of the filter is illustrated in Figure 4.21.

It can be noticed that the coefficients obtained from the C code are not exactly the same as those from MATLAB. Since the source code in MATLAB is not available, it is not possible to check the differences in the specific calculations. The frequency responses in Figure 4.18 and Figure 4.21 also show the differences. It looks like the response from C code is closer to the ideal filter response (Figure 4.16).

75

Figure 4.20 Coefficients from the C program.

76

Figure 4.21 Frequency response for a M = 100 Hilbert transform filter.

This filter should also have a linear phase response as the one in Figure 4.18. However the phase response of the plot does not have the same slope in the stopband as one in the passband. It is caused by the calculation error because the transfer function values in the stopband are very small.

4.3.3 Digital Filter According

to

the

FIR

filter

definition,

the

function

fir _ filter ( M , Coef , Xinput ) was written in C program that is used to apply filtering operations on the data sequence in vector Xinput , where

77

M is an integer. It is the number of coefficients of the filter.

Coef is the coefficients of the filter. Xinput is the input data.

The following codes illustrate how to use this function.

where “xin” contains the input data, it could be the ECG data. “infile” is the file name which saved the input data. It is named by users. “h” is the coefficients obtained from the remez.c . “outfile” is the file name which saved the results filtered from the filter. It is named by user. In this case, the outfile saved the Hilbert transform of the input ECG data.

4.3.4 Detector for R Wave Peak Accurate determination of the QRS complex, in particular accurate detection of the R wave peak, is essential in computer-based ECG analysis [24]. As described in the previous chapter, one of the properties of the Hilbert transform is that it is an odd function. That is to say that it will cross zero on the x -axis every time that there is an inflection point in the original waveform [24]. Similarly a crossing of the zero between consecutive positive and negative inflection points in the Hilbert transformed conjugate will be represented as a peak in its original

78

waveform. Using this characteristic, a detector for determining the R wave peak in the input ECG waveform was developed. The Hilbert transform of the ECG wave was obtained in section 4.3.3. The peaks in the Hilbert transform sequence h(n) represent regions of high probability of finding R wave peaks. An adaptive threshold is used to locate the peaks in the h( n) sequence. For finding the R wave peak accurately, a moving 1000 points

window is used to subdivide the Hilbert transform h(n) sequence. The RMS (Root Mean Square) value and the maximum amplitude in the present window are then calculated. The threshold in this window was determined according to the criteria below: 1. If (the RMS value) ≥ (18%*maximum amplitude) in the Hilbert transform sequence, the threshold is set up at (39%*maximum amplitude). 2. If (the present maximum amplitude) ≥ (2* previous maximum amplitude), the threshold will be (39%*previous maximum amplitude). 3. If (the RMS value) < (18%* maximum amplitude) in the Hilbert transform sequence, the threshold will be (1.6*RMS value). 4. If the two peaks in the h(n) sequence are two close together, only one of them is the real R peak.

4.3.5 R Wave Peak Detection Test In this section, two ECG waveforms taken from the MIT-BIH Arrhythmia database are used as the test signals.

79

The MIT-BIH Arrhythmia database consists of 48 records, each containing 30 minutes of two-channel ECG with heartbeat and rhythm annotations. The recordings were digitized at 360 samples per second per channel with 11-bit resolution over a 10 mV range. All samples are represented as positive numbers. The entire 30-minute record is annotated. XW_1 and XW_2 are the normal ECG waves. They are the small sections of MIT213 from the MIT-BIH database. Figures 4.22 and 4.23 illustrate the input data XW_1 and its filtered version, i.e., the Hilbert transform of the input signal.

Figure 4.22 Input wavefom: XW_1.

As discussed in the previous chapter, the standard ECG is a representation of the heart electrical activity recorded from electrodes on the body face. Figure 4.22

80

shows a small section beats of a normal heart. There are 4 heartbeats since each R wave represents a heartbeat.

Figure 4.23 Output waveform for XW_1. Table 4.3 XW_1 extreme points and values

Input (Figure 4.22)

Output (Figure 4.23)

N

Extreme Value

N

Value

( N output − N Input − 50)

104

1459

155

≈0

1

329

1419

380

≈0

1

524

1404

575

≈0

1

724

1435

774

≈0

0

81

Error

Figure 4.23 shows the Hilbert transform of the input wave, i.e. XW_1. For those points in Figure 4.22 where the slope changed from positive to negative, their outputs should be the zero crossing points with the extreme values changing from negative to positive. Because of the filter delay, the output results lag the input by 50 ( M / 2) . Table 4.3 shows N and extreme values corresponding to every heartbeat, i.e., the R wave peak. Compared the results with the MIT-BIH database annotation file, it can be seen that the output values are close. Figure 4.24 and Figure 4.25 are for data XW_2.

Figure 4.24 Input waveform: XW_2.

82

Figure 4.25 Output waveform for XW_2.

Figure 4.24 shows another channel waveform for the same person. The reason why Figure 4.24 is different from the Figure 4.22 is the electrode is placed at a different location. XW_1 is obtained when the electrode is placed at the front of the body while XW_2 is obtained when the electrode is placed on the back of the body. For this case, for those points which slope changed from negative to positive, the output point should be the zero crossing point with the extreme values changing from positive to negative.

83

Table 4.4 XW_2 extreme points and values

Input (Figure 4.24)

Output (Figure 4.25)

N

Extreme Value

N

Value

( N output − N Input − 50)

103

518

152

≈0

-1

327

466

377

≈0

0

522

436

572

≈0

0

722

416

771

≈0

-1

Error

Table 4.4 shows the results, for the same reason with the previous wave, filter phase delay, the output wave lag the input by 50 ( M / 2 ). Comparing the output results with the MIT-BIH database annotation file, the error between the output results and the annotation results is acceptable.

4.4 Nios Implementation In order to use the Nios development board [Figure 4.4], connect the host PC and open a Nios SDK shell and type $ nios − run − t − r This establishes a simple terminal connection with the development board. Press the SafeConfig button on the Nios development board to reset the Nios development board and reconfigure the Stratix FPGA. The reference design emits a text message to the serial port when the Nios processor boots. After the LEDs begin

84

to blink and the LED displays network-initialization status messages, press SW3 to abort DHCP network configuration. Text will display in the Nios SDK Shell window. Press the Enter key on the PC several times to provide stimulus to the reference design. The interface is shown in Figure 4.26.

Figure 4.26 Nios SDK Shell Prompt.

85

If the activity in the Nios SDK shell looks like in Figure 4.26, then the PC is communicating correctly with the Nios development board. Press Ctrl + C to exit the terminal program and return to the bash shell. To compile all programs at the at the bash prompt, use command $ nios − build The GNU C/C++ compiler and linker will be invoked. Several intermediate files and an executable ( .srec ) file will be produced. The messages are shown in Figure 4.27.

Figure 4.27 Nios-Build messages.

86

To download and run the code compiled in Figure 4.27, download the srec file to the Nios development board. Redirect stdout to a data file outputdata1.txt . The message is shown in Figure 4.28.

Figure 4.28 Nios SDK shell prompt.

The input waveforms XW_1 and XW_2 are a section wave of MIT 213 taken from the MIT_BIH Arrhythmia database (Figure 4.29 and Figure 4.31). The output file record the data of the Hilbert transform of the ECG wave. The output waveforms are shown in the Figure 4.30 and Figure 4.32.

87

Figure 4.29 Input ECG waveform: XW_1.

Figure 4.30 The Hilbert transform of XW_1.

88

Figure 4.31 Input ECG waveform: XW_2.

Figure 4.32 The Hilbert transform of XW_2.

89

From the output waveforms in Figure 4.30 and Figure 4.32, it can be seen that the waveforms are almost the same as the results obtained from Visual C++.

Table 4.5 XW_1 R wave points and their Hilbert transform points

Error

Input (Figure 4.29)

Output (Figure 4.30)

N

Extreme Value

N

Value

( N output − N Input − 50)

104

1459

155

≈0

1

329

1419

380

≈0

1

524

1404

575

≈0

1

724

1435

774

≈0

0

Table 4.6 XW_2 R wave points and their Hilbert transform points

Error

Input (Figure 4.31)

Output (Figure 4.32)

N

Extreme Value

N

Value

( N output − N Input − 50)

103

518

152

≈0

-1

327

466

377

≈0

0

522

436

572

≈0

0

722

416

771

≈0

-1

Table 4.5 and Table 4.6 show the N values of the R wave points and their Hilbert transform points. It’s the same as the one discussed in the section 4.2. The results illustrates that the programs can work correctly on the Nios embedded processor.

90

Chapter 5 Results

In the previous chapter, implementation of the algorithm for the Hilbert transform of an ECG wave has been discussed. In order to prove the program is working well in the Nios processor, some ECG waves taking from the MIT-BIH Arrhythmia database were used as verifications. The results will be compared with the annotation files recorded in the MIT-BIH database.

5.1 Experimental Results Further examples of some ECG waves and their Hilbert transformed output waveforms are shown in this part. The file mit 212 from100.txt is a small section from MIT212 recorded in the MIT-BIH arrhythmia database. It includes 1000 data samples from 100 to1099 in MIT212. Figure 5.1 shows the original ECG waveform. The output, i.e., the Hilbert transform

of

mit 212 from100.txt

is

shown

in

Figure

5.2.

The

file

mit 212 _ 1 from47200.txt is another section from the MIT212. It also includes 1000

91

data samples. The original waveform and the Hilbert transform of this ECG wave are shown in Figure 5.3 and Figure 5.4. As discussed preciously, accurate detection of the R wave peaks is important in ECG analysis. The zero crossing points corresponding to the true R wave peak in the Hilbert transformed data of the original ECG waveform are shown in Table 5.1 and Table 5.2. The actual position of R wave points in the MIT-BIH arrhythmia database is also given. The error between the located R and actual R is calculated using Equation 5.1: Error = N output − N MIT − BIH − 50 ,

(5.1)

where N ouput is the number of the zero crossing point in the Hilbert transform

sequence of the original ECG data located by detector. N MIT − BIH is the actual R wave peak location recorded in the MIT-BIH

annotation file. Since the order of the filter used in this thesis is 100, the output waveform lagged the input waveform 50.

92

Figure 5.1 Input waveforms: mit212from100.txt.

Figure 5.2 Output waveform: outMIT212from100.txt.

93

Table 5.1 mit212from100 extreme points and values Input

Output

MIT-BIH

Error

N

Extreme Value

N

Value

Value

116

1287

165

≈0

114N

1

353

1294

402

≈0

351N

1

597

1265

647

≈0

596N

1

825

1218

875

≈0

824N

1

Figure 5.3 Input waveform: mit212_1from47200.txt.

94

Figure 5.4 Output waveform: outMIT212from47200.txt.

Table 5.2 mit212_1from47200 extreme points and values Input

Output

MIT-BIH

Error

N

Extreme Value

N

Value

Value

48

1255

97

≈0

46N

1

255

1169

304

≈0

254N

0

477

1294

527

≈0

476N

1

714

1324

763

≈0

713N

0

963

1330

95

The ECG waveforms (Figure 5.1 and Figure 5.3) are normal because the P waves, QRS complexes, ST segments, T waves and U waves are normal. Each R wave peak presents a heartbeat. Table 5.1 and Table 5.2 show the positions of the zero crossing points in the output waveform which corresponding to the R peaks in the input ECG wave. Because of the filter delay, the output results lag the input by 50 ( M / 2) . From Table 5.1 and Table 5.2, it can be seen that the N in the output waveform is lagging about 50. The error between the zero crossing points obtained in the output waves with the value recorded in the MIT-BIH annotation file is very small. It proves that the algorithm for calculating the Hilbert transform of the original ECG wave is correct. Note in the Table 5.1 and Table 5.2, the meaning of the letter in the “MIT-BIH value” is the type of the beat recorded in the MIT-BIH annotation file. N means normal QRS, V means premature ventricular contraction (PVC). F means fusion PVC. The details about the MIT-BIH annotation code are shown in the Appendix A. The files mit 213 _ 1 from67000.txt and mit 213 _ 1 from1450.txt are small sections from another record, MIT213, in the MIT-BIH arrhythmia database. Figure 5.5 and Figure 5.7 show the ECG waveform. The outputs, i.e., the Hilbert transform of the input wave, were illustrated in Figure 5.6 and Figure 5.8. It can be seen there is an irregular heartbeat between N = 600 and 700 in Figure 5.5. This kind of beat is called a PVC (Premature Ventricular Contraction). PVCs are premature heartbeats originating from the ventricals of heart. PVCs are early or extra heartbeats that commonly occur and are usually harmless in normal hearts, but can

96

cause problems in hearts with pre-exiting disease. A person with PVCs may or may not feel the irregular heartbeat, usually as a skip heartbeat [25]. The characteristic of PVCs is that there is no P wave and PR interval and the QRS complex is greater than 0.12s.

Figure 5.5 Input waveform: mit213_1from67000.txt

97

Figure 5.6 Output waveform: outMIT213_1from67000.txt.

Table 5.3 mit213_1from67000 extreme points and values Input

Value

Output

Error

N

ExtremeValue

N

Value

(MIT-BIH)

51

1451

101

≈0

50 N

1

244

1549

293

≈0

243 F

0

439

1463

489

≈0

438 N

1

637

1582

689

≈0

635 V

4

834

1476

884

≈0

833 N

1

98

Figure 5.7 Input waveform: mit213_1from1450.txt.

Figure 5.8 Output waveform: outMIT213_1from1450.txt.

99

Table 5.4 mit213_1from1450 extreme points and values Input

Output

Value

Error

N

ExtremeValue

N

Value

(MIT-BIH)

17

1388

67

≈0

16 N

1

206

1388

255

≈0

205 N

0

400

1396

450

≈0

399 N

1

598

1415

647

≈0

597 N

0

From the results obtained in Table 5.3 and Table 5.4, it can be seen that the error for normal heartbeat and for the PVC beat are small and are acceptable.

Here other examples are given, the files mit 223 _ 1 from47700.txt and

mit 223 _ 1 from7900.txt are sections from the MIT 223 in the MIT-BIH database. The input waves are shown in Figure 5.9 and Figure 5.11, the Hilbert transforms of the input ECG waves are shown in Figure 5.10 and Figure 5.12. In this case, there are a few irregular heartbeats, such as PVCs. In Figure 5.9, the PVCs happened between N = 0 to 100 and N = 800 to 900. Table 5.5 and Table 5.6 show the compared results.

100

Figure 5.9 Input waveform: mit223_1from47700.txt.

Figure 5.10 Output waveform: outMIT223_1_from47700.txt.

101

Table 5.5 mit223_1from47700 extreme points and values Input

Output

Value

Error

N

Extreme Value

N

Value

(MIT-BIH)

60

1376

110

≈0

60 V

0

393

1327

442

≈0

391 N

1

667

1321

716

≈0

664 N

2

862

1367

912

≈0

861 V

1

Figure 5.11 Input waveform: mit223_1from7900.txt.

102

Figure 5.12 Output waveform: outMIT223_1_from7900.txt.

Table 5.6 mit223_1from7900 extreme points and values Input

Output

Value

Error

N

ExtremeValue

N

Value

(MIT-BIH)

79

1303

128

≈0

76 N

2

352

1328

401

≈0

349 N

2

591

1143

641

≈0

590 V

1

879

1298

928

≈0

876 N

2

103

In Figure 5.11, the PVC happened between N = 500 to 600. The R wave peak is not that clear, there are two peaks and the extreme values are close together. Table 5.5 and Table 5.6 record the R wave peaks in the input and output waveform. From Table 5.5 and Table 5.6, it can be seen that the output results lag the input data by 50 because the filter delay is 50 ( M / 2) . As discussed in the Chapter 4, the points which are of interest are those zero crossing points with the extreme value changing from negative to positive or from positive to negative. Comparing the results shown in Table 5.1-Table 5.6 with the data recorded in the MIT-BIH arrhythmia database, the error is small and is acceptable. It illustrated that the algorithm and the program are correct and working properly.

5.2 Complete ECG Testing In the previous section, just small sections from the record of the MIT-BIH Arrhythmia database are tested. As mentioned previously, each MIT-BIH excerpt contains 30 minutes of ECG with heartbeat recording. In this section, the whole excerpt will be used to test. Every excerpt includes 649999 data samples. MIT212 is the heart recording of an adult (31 years old). The predominant rhythm is normal sinus at rates of 75-90 BPM.

104

MIT213 is the record of a 61year_old person’s heart; heart rate is 100-110 BPM. Predominant is normal but there are occasional APBs (Atrial Premature Beats). MIT223 is the record of an 84 year_old person’s heart. APBs are present throughout. There is high-grade ventricular ectopic activity with frequent multifocal VPBs, couplets, and runs of Vtach. The entire excerpt from the MIT-BIH Arrhythmia database was tested by the system discussed in this thesis. The Hilbert transform of this excerpt is obtained and the R wave peaks are detected. Since the data is too long, just portion of the results of MIT212, MIT213 and MIT223 were shown in Figure 5.13, Figure 14-15 and Figure 5.16-17. Beat by beat comparison was performed according to the annotation file. The results are shown in Table 5.7.

105

Figure 5.13 The portion of results of MIT212.

106

Figure 5.14 The portion of results of MIT213.

107

Figure 5.15 The portion of results of MIT213.

108

Figure 5.16 The portion of results of MIT223.

109

Figure 5.17 The portion of results of MIT223. Table 5.7 R wave detection performance MIT-BIH record

Failed detection number 10

Detection error rate

Average error

MIT212

Actual number of beats in record 2748

0.00364

1.11

MIT213

3251

9

0.00277

1.26

MIT223

2605

10

0.00384

2.43

110

The detection error rate and average error were calculated using the following equations respectively: DetectionErrorRate =

FailedDetectionNumber AcutralNumberOfBeats

⎛ k ⎞ ⎜ ∑ locatedR − actualR ⎟ ⎠, AverageError = ⎝ i =1 K

(5.2)

(5.3)

where

K is the total number of R correctly located by the detector.

From the detection error rate and the average error, it can be seen that the results are closed to the data recorded in the MIT-BIH annotation file. The average error in MIT 223 is bigger since MIT 223 records the heart beats of an 83 year old person and there are a lot of PVCs, the beat by beat error is larger. That is why the average error is larger than others. But the error is still acceptable (< 3 samples). The results of the test illustrate the algorithm for calculating the Hilbert transform and detecting the R wave peaks is effective.

111

Chapter 6 Summary and Conclusion In the previous chapters, the information about the Holter monitoring system and the process for calculating the Hilbert transform of an ECG signal are discussed in detail. This final chapter summarises the results of the research and the contents of the thesis.

6.1 Summary The Holter ECG Monitoring System mainly consists of four parts: Input Processing, Pattern Recognition, Compression and Storage. The main objective of the input processing is to calculate the Hilbert transform of the input ECG data. Pattern recognition uses vectorcardiograph and polarcardiograph representations and the concepts of pre-envelope and envelope of a real waveform given by the Hilbert transform to judge whether the ECG wave is normal or abnormal. Lastly data is compressed from the abnormal ECG wave and saved in the flash card. In this thesis, only the first part of the Holter ECG monitoring system, i.e., Input Processing is discussed.

112

The background of the Holter ECG Monitoring System was presented. The Hilbert transform applications and the basic mathematics and properties of a Hilbert transformer are also presented. The Hilbert transform is a mathematical method for analysing signal waveforms, and has been widely used in the areas of communication systems analysis. The resulting display of the Hilbert transformed data is similar to that obtained from conventional vectorcardiographic systems. It allows easy visual indication of the different classes of normal and abnormal morphologies. It allows quick, precise segmentation of the incoming ECG into individual heartbeats and also allows the detection of Q-, R-, S-, and T-wave complexes in the data. In this thesis, four approaches to compute the Hilbert transform including the Time-Domain approach, the Frequency-Domain approach, the Boche approach and the Remez Filter approach have been discussed in detail. The algorithms are deduced and examples are also given for every approach. After comparing them in running time and the ease of implementation, the Remez Filter approach which implements the Parks-McClellan algorithm to design and apply a linear-phase filter was determined to be the best and is used in computing the Hilbert transform of the ECG wave in this thesis. The results for every approach were shown in Chapter 3 and the comparisons were also given. In the software implementation part, as a first step, the filter order and the frequency range are determined by analyzing a set of test results using MATLAB. The remez function was used in MATLAB to test the Hilbert transform of the ECG data.

113

In additional, a C program ( remez.c ) was developed to implement the algorithm of the Remez filter approach. remez.c uses the Parks-McClellan algorithm, i.e., uses the Remez exchange algorithm and Chebyshev approximation theory to design a filter with the optimal fit between the desired frequency response and actual frequency response. A main program that handles the input, sets up the appropriate approximation problem and handles the output of the optimal filter coefficients was included. Specifically remez.c has a build-in section for the more common ideal filter types such as multi-band, bandpass filters, differentiators and the Hilbert transform filters. All results including the order, frequency range and the coefficients of the filter were shown in Chapter 4. A digital filter was developed to apply filtering operations on the data sequence. Using the coefficients achieved from remez.c and the input data, i.e., the ECG data sequence, the Hilbert transform of the input ECG data was obtained. All results and waveforms are also presented in Chapter 4. Accurate determination of the QRS complex, in particular, accurate detection of the R wave peak, is essential in ECG analysis and is another task in this thesis. In the system discussed in this thesis, the method of ECG waveform analysis uses vectorcardiograph and polarcardiograph representations and examines the concepts of pre-envelope and envelope of a real waveform given by the Hilbert transform. A prototype two stage QRS detector was used based on the determination of a zero crossing in the Hilbert transformed data of the original ECG waveform. The positions of the zero crossing points that correspond to the R wave peaks are useful in judging whether the heartbeat is normal or abnormal.

114

A C program was developed to detect the QRS complex, in particular, to detect the R wave peak used in ECG analysis. All the results were shown in Chapter 4 and Chapter 5. The simulation results are also presented in Chapter 5. The Nios embedded processor was introduced. SOPC Builder was used to be development tool. Using SOPC builder, users can combine the Nios processor with user logic and program it into an FPGA easily. All programs were run on the Nios embedded processor. The results, i.e., the Hilbert transform of the ECG data sequence, are almost the same as the results obtained in the C program except the first 50 data samples. The first 50 data samples in the output sequence can be ignored since the phase delay of the filter is 50. The results prove that the program can work properly in the Nios embedded processor.

6.2 Conclusion The data used in test were from standard ECG waveform records in the MITBIH arrhythmia database. The performance of the chosen algorithm was tested. The test results were compared with the annotation files recorded in the MIT-BIH arrhythmia database. The detection error rate is smaller than 0.005. The average error is smaller than 3 samples. The results were given in Chapter 5. The error is acceptable. The results illustrated that the algorithm performed effectively with accurate R wave detection.

115

References [1] Alan E. Lindsay, “ECG Learning Center in Cyberspace”, Available at: http://medlib.med.utah.edu/kw/ecg/ecg_outline/Lesson1/index.html Accessed February 2006. [2] Richard E. Klabunde, “Cardiovascular Physiology Concepts”. Available at: http://www.cvphysiology.com/Arrhythmias/A009.htm Accessed March 2006. [3] Molson Medical Informatics Student Project, “Cardiac Arrythmias”, Available at: http://sprojects.mmi.mcgill.ca/cardiophysio/sinustachycardia.htm Accessed March 2006. [4] “Electrocardiograms”, Available at: http://www.healthyhearts.com/ecg.htm Accessed March 2006. [5] Merck& Co., Inc., “The Merck Manual of Health & Aging”, Available at: http://www.merck.com/pubs/mmanual_ha/figures/fg47_1.html

Accessed

March 2006. [6] Aetna, Inc. Available at: http://www.aetna.com/cpb/data/CPBA0019.html Accessed March 2006. [7] Wikipedia, the free encyclopedia, “Cauchy principal value”, Available at: http://en.wikipedia.org/wiki/Cauchy_principal_value Accessed April 2006.

116

[8] Kunio Takaya, “Digital Signal Processing Electrical Engineering EE-880”, University of Saskatchewan, 2004. [9] Mathias Johansson, M.Sc. thesis, “The Hilbert transform”, VaxJo University, 1999. [10] Easy Fourier Analysis, “Signal processing & simulation newsletter”, Available at: http://www.complextoreal.com/tcomplex.htm Accessed April 2006. [11] Ronald

J.

Bolton,

Ph.D.

thesis,

“Hilbert

transform

processing

of

Electrocardiograms”, University of Queensland, 1983. [12] Jose A. Inaudi and James M. Kelly, “Linear hysteretic damping and the Hilbert transform”, Journal of Engineering Mechanics, Vol.121, Issue 5, pp.626-632, May 1995. [13] Holger Boche and Marcus Protzmann, “A new algorithm for the reconstruction of bandlimited functions and their Hilbert transform”, IEEE Transactions on Instrumentation and Mesurement, Vol. 46. No.2, pp. 442-444, April 1997. [14] Altera

website,

http://www.altera.com/products/ip/processors/nios2/ni2-

index.html Accessed April 2006. [15] Altera

Corporation,

“Nios

Hardware

Development

Tutorial”,

TU-

Software

Development

Tutorial”,

TU-

NIOSHWDV-1.2, January 2004. [16] Altera

Corporation,

“Nios

NIOSSFTWR-1.3, July 2003. [17] Altera Corporation, “Nios II Development Board Reference Manual, Stratix Professional Edition”, MNL-N2DEVLSTXPRO-1.1, September 2004.

117

[18] Altera Corporation, “Nios II Development Kit Getting Started User Guide”, UG-NIOSIIDEVKIT-2.2, May 2006. [19] A. Robin, “Digital Filters: An Introduction”, TechOnLine, Available at: http://www.techonline.com/community/ed_resource/feature_article/20365 Accessed May 2006. [20] Alan V. Oppenheim and Ronald W. Schafer, “Discrete-Time Signal Processing”, Prentice Hall, Englewood Cliffs, NJ, 1989. [21] Wolfram Research, “Digital Image Processing Documentation”, Available at: http://documents.wolfram.com/applications/digitalimage/UsersGuide/FilterDes ign/ImageProcessing9.2.html , Accessed May 2006. [22] “Parks-McClellan

algorithm

for

FIR

filter

design”,

Available

http://search.cpan.org/src/MLEHMANN/PDL-Audio-1.1/remez.c

at:

Accessed

May 2005. [23] The digital signal processing committee, “Programs for Digital Signal Processing”, IEEE Acoustics, Speech, and Signal Processing Society, The Institute of Electrical and Electronics Engineers INC., pp5.1-1 to pp.5.1-13, 1979. [24] D. Benitez, P.A. Gaydecki, A. Zaidi, A.P. Fitzpatrick “The use of the Hilbert transform in ECG signal analysis”, Computers in Biology and Medicine , Vol.31, pp.399-406, 2001. [25] Encyclopedia, “Premature Ventricular Contractions”, Komo 1000 news, Available at: http://ww3.komotv.com/Global/story.asp?S=1230267 Accessed July 2006.

118

Appendix A The MIT-BIH Arrhythmia Database

A.1 Introduction The database consists of 48 records, each containing 30 minutes of twochannel ECG with beat and rhythm annotations. Each digital record has been copied from an analog recording made with an Avionics 445 two-channel recorder. Annotations have been made by two independent cardiologists with consultation to resolve disagreements. The data base is recorded on twelve 2400-foot (730m) ANSI standard 9-track tapes at 800 bpi, odd parity, with NRZI recording. A detailed catalog of the contents of each tape, with illustrations, is included with the database. The format is that which will be used for the American Heart Association (AHA) Database for Ventricular Arrhythmia Detectors, with these differences: Each 9-track tape contains four 30-minute records. Sampling frequency is 360 Hz per channel. Sampling precision is 11 bits, and all samples are represented as positive numbers. The entire 30-minute record is annotated. Annotations are referenced to samples (rather than milliseconds).

119

An additional “0” annotation has been added to the AHA set to specify nonbeat annotations (e. g., rhythms, artifact).The “R” (R-on-T PVC) annotation is not used. Space allocated but unused in the AHA format annotation blocks is used to specify rhythms and beat types more precisely than is allowed using the AHA annotation codes alone. Atrial ectopic beats and conduction defects are among the items specified in this way.

A.2 File Structure Each tape contains 16 files separated by ANSI standard end-of-file (EOF) marks. The last file is terminated with two EOFs to indicate the end of the tape. Each record, corresponding to 30 minutes of ECG and annotations, is comprised of four files: an ID block file, a sample data file, an annotation file, and a second ID block file. The order of files is: (Record 1)

(Record 2)

ID Block (EOF) Sample Data Blocks (EOF) Annotation Blocks (EOF) ID Block (EOF) ID Block (EOF) Sample Data Blocks (EOF) Annotation Blocks (EOF) ID Block (EOF)

120

(Record3)

(Record4)

Sample Data Blocks (EOF) Annotation Blocks (EOF) ID Block (EOF) Sample Data Blocks (EOF) Annotation Blocks (EOF) ID Block (EOF)

A.3 Notational and Other Conventions Multiple-byte Numbers In this specification, the least significant 8-bit byte of a multiple-byte number is referred to as byte i, the next most significant byte as byte ii, and so on. (The first byte read from the tape in a given block is called byte 1.) The AHA format specifies that: 16-bit numbers are stored in the order: byte i, byte ii. 32-bit numbers are stored in the order: byte iii, byte iv, byte i, byte ii.

ASCII text In the ID block and in certain annotation labels (see below) brief comments are present. These are coded as ASCII characters, and should be read from the tape in byte-sequential order.

TOCs “TOC” means “time of occurrence”, TOCs are always represented as 32-bit numbers. TOCs in the AHA database are given as the number of MILLISECONDS

121

from the beginning of the annotated segment of the record. In the MIT-BIH database, TOCs are given as the number of SAMPLE COUNTS from the beginning of the record. To convert sample counts to milliseconds, multiply sample counts by 1000/360 (=2.777….).

A.4 File Format Specifications ID block file The first and fourth files in each record each consist of a single 512 byte ID block, The AHA specification for the ID block is: Bytes

Use

1-8 9-10 11-16 17-20 21-24 25-26 27-32 33-36 37-40 41-42 43-512

record ID (8 ASCII characters) number of annotations unused time of first sample in the annotated segment of the record time of last sample in the annotated segment of the record number of bytes of sample data, divided by 512 and rounded upward unused TOC-first annotation, relative to the beginning of the annotated segment TOC-last annotation, relative to the beginning of the annotated segment numbers of bytes of annotation data, divided by 512 and rounded upward unused

In each record in the MIT-BIH database, the entire record has been annotated; thus the time of the first sample is always zero. Each tape has exactly 649999 samples, and an end-of-sample-data mark, per channel (30 minutes and 5.444 seconds), and the annotated segment is considered to end after the end-of-sampledata mark, so that the time of the last sample in the annotated segment is always 650000, and the number of bytes of sample data divided by 512 is always 5079

122

(649999 samples per channel, times 2 channels, times 2 bytes per sample, divided by 512, rounded up). NOTE THAT THE UNITS OF TIME ARE SAMPLE COUNTS, NOT MILLISECONDS

Sample data file The second file in each record is the sample data file, which consists of exactly 2540 blocks, each 1024 bytes long. Each block contains 256 2-byte samples from each channel. Samples are stored alternately in the block: Byte

Use

1 2 3 4 5 M 1021 1022 1023 1024

Channel 1, sample 1, byte i Channel 1, sample 1, byte ii Channel 2, sample 1, byte i Channel 2, sample 1, byte ii Channel 1, sample 2, byte i Channel 1, sample 256, bytei Channel 1, sample 256, byteii Channel 2, sample 256, bytei Channel 2, sample 256, byteii

The AHA database has been recorded using a 12-bit A/D converter with a range of -10V to +10V, and preamp gain adjusted so that a QRS complex is nominally 1V peak-to-peak, or about 200 ADC units. The MIT-BIH database has been recorded using an 11-bit A/D converter with a range of -5mV to +5mV, and the unamplified QRS complexes are nominally 1 mV, or about 200 ADC units; thus the scales are the same though the ranges differ. Both positive and negative (two’s complement, with sign extension to 16 bits) samples are recorded in the AHA database; in the MIT-BIH database, all samples are positive (in the range of 0 to 2047). 123

Block 2540 contains the last fifteen samples for each channel. The end of the sample data is marked in the last sample block by two consecutive sample values of 10000 (base 8) following the last samples. The remainder of the last sample block is padded with zeroes.

Annotation file The third file in each record is the annotation file, which consists of a variable number (typically 20 to 50) of blocks, each 1024 bytes long. Each block contains 64 annotation labels, each 16 bytes long. Annotation labels are stored in strict chronologic order. The AHA format leaves a number of unused bytes in each annotation label, some of which are used in the MIT-BIH database. The last annotation block is padded with all ones (177 base 8) following the last annotation. If there is no room following the last annotation, an entire block of 177s is written.

A.5 Annotation Specifications Annotation labels Byte

AHA format

MIT-BIH format

1 2 3-6 7-8

unused AHA annotation code TOC (milliseconds) Annotation label serial number unused unused unused

unused(0) AHA annotation code TOC (sample counts) Annotation label serial number unused (0) MIT-BIH annotation code ASCII text *

9 10 11-16

124

*The ASCII text field is filled with zero bytes unless the MIT-BIH annotation code is 22 or 28 (see next page following).

AHA annotation codes The AHA annotation codes are ASCII characters: Character

Value (base 8) Meaning

N V E F R P Q U [ ]

116 126 105 106 122 120 121 125 133 135

supraventricular beat premature ventricular contraction (PVC) ventricular escape beat fusion PVC R-on-T PVC paced beat beat of indeterminate origin data unreadable between preceding and following beat labels beginning of ventricular flutter or fibrillation end of ventricular flutter/fibrillation

The “R” code does not appear in the MIT-BIH database. An additional code, “O” (117 base 8), has been defined to permit inclusion of rhythm labels, artifact labels, and comments. “O” labels are never QRS labels, and may be ignored for the purpose of counting beats.

MIT-BIH annotation codes The MIT-BIH annotation codes are not ASCII characters, but numbers between 1 and 37:

Code 1 2 3

AHA equivalent Meaning N N N

normal QRS left bundle branch block beat right bundle branch block beat 125

4 5 6 7 8 9 10 11 12 13 14 15 16 17-21 22 23-24 25 26 27 28 29-30 31 32 33 34 35 36 37 38 *

N V F N N N E N P Q O,U O O O O O N O O O O O [ ] N N O O O

aberrantly conducted beat premature ventricular contraction (PVC) fusion PVC *** nodal premature beat atrial premature bat (APB) nodal or atrial premature beat ventricular escape beat nodal escape beat paced beat beat of indeterminate origin beginning of noise * end of noise * single QRS-like artifact reserved for future use ** comment (text) annotation *** reserved for future use ** left or right bundle branch block beat non-captured pacemaker spike axis shift rhythm onset (text) annotation *** reserved for future use ** ventricular flutter wave onset of ventricular flutter or ventricular fibrillation end of ventricular flutter/fibrillation atrial ectopic beat atrial ectopic beat nodal ectopic beat missed beat blocked APB reserved for future use **

Annotation codes 14 and 15 are used in pairs. If AHA code corresponding to

the code 14 “U”, no beats are labeled until the next code 15; otherwise, all beats are labeled. **

Where codes designated “reserved for future use” appear, they should be

ignored. *** “Text” annotations use the last six bytes of the annotation label for an ASCII string (terminated by a zero byte if less than six characters).

126

**** In the context of paced rhythm (tapes 102,104,107,217) annotation code 6 is used for pacemaker fusion beats.

Rhythm onset annoations Rhythm onset annotation (MIT-BIH annotation code 28) include an ASCII string which begins with a “(“: String

Meaning

(AB (AFIB (AFL (B (BI (BII (BIII (IVR (N (NOD (PAT (PREX (SVTA (T (VFIB (VFL (VT

atrail bigeminy atrial fibrillation atrial flutter ventricular bigeminy first degree heart block second degree heart block third degree heart block idioventricular rhythm normal sinus rhythm normal (A-V junctional) rhythm paroxsysmal atrial tachycardia pre-excitation (WPW) supraventricular tachyarrhythmia ventricular trigeminy ventricular fibrillation ventricular flutter ventricular tachycardia

Comment annotations Sparse comment annotations exist on a few records. They are: PSE TS

pause tape slippage

127

Appendix B

B.1 The Parks-McClellan Algorithm Consider a particularly effective and widely used algorithm procedure for the design of FIR filters with generalized linear phase. Although only type I filters are considered in detail, where appropriate results that apply to types II, III, and IV generalized linear phase filters are indicated. In designing a causal type I linear phase FIR filter, it is convenient first to consider the design of a zero-phase filter, i.e., one for which he [n] = he [−n],

(1)

and then to insert sufficient delay to make it causal. Consequently, consider he [n] satisfying the condition of Eq. (1). The corresponding frequency response is given by Ae (e jω ) =

L

∑ h [n]e

n=− L

e

− jωm

,

(2)

with L = M / 2 an integer or, because of Eq. (1), L

Ae (e jω ) = he [0] + ∑ 2he [n] cos(ωn).

(3)

n =1

Note that Ae (e jω ) is a real, even, and periodic function of ω . A causal system can be obtained from he [ n] by delaying it by L = M / 2 samples. The resulting system has impulse response

128

h[n] = he [n − M / 2] = h[ M − n]

(4)

H (e jω ) = Ae (e jω )e − jωM / 2 .

(5)

and frequency response

The Parks-McClellan algorithm is based on reformulating the filter design problem as a problem in polynomial approximation. Specially, the terms cos(ωn) in Eq. (3) can be expressed as a sum of powers of cos ω in the form cos(ωn) =T n(cos ω ),

(6)

where Tn ( x) is an n th-order polynomial. Consequently, Eq. (3) can be rewritten as an L th-order polynomial in cos ω . Specially, L

Ae (e jω ) = ∑ a k (cos ω ) k ,

(7)

k =0

where the a k ’s are constants that are related to he [ n] , the values of the impulse response. With the substitution x = cos ω , Eq. (7) can be expressed as Ae (e jω ) = P ( x)

x = cos ω

,

(8)

where P( x) is the L th-order polynomial L

P ( x) = ∑ a k x k .

(9)

k =0

It is not necessary to know the relationship between the a k ’s and he [n] ; it is enough to know that Ae (e jw ) can be expressed as the L th-degree trigonometric polynomial of Eq.(7).

129

The key to gaining control over ω p and ω s is to fix them at their desired values and let δ 1 and δ 2 vary. To formalize the approximation problem, define an approximation error function E (ω ) = W (ω )[ H d (e jω ) − Ae (e jω )],

(10)

where the weighting function, W (ω ) , incorporates the approximation error parameters into the design process. In this design method, the error function E (ω ) , the weighting function W (ω ) , and the desired frequency response H d (e jω ) are defined only over closed subintervals of 0 ≤ ω ≤ π . Parks and McClellan applied the following Alternation Theorem of approximation theory to the filter design problem. Let FP denote the closed subset consisting of the disjoint union of closed subsets of the real axis x . P( x) denotes an r th-order polynomial r

P ( x) = ∑ a k x k . k =0

Also, DP (x) denotes a given desired function of x that is continuous on FP , and E P (x) denotes the weighted error E P ( x) = WP ( x)[ DP ( x) − P ( x)].

The maximum error E is defined as E = max E P ( x) . x∈FP

A necessary and sufficient condition that P( x) is the unique r th-order polynomial that minimizes E is that E P (x) exhibit at least (r + 2) alternations, i.e., there

130

must exist at least (r + 2) values xi in FP such that x1 < x 2 < L < x r + 2 and such that E P ( xi ) = − E P ( xi +1 ) = ± E for i = 1,2,K, (r + 1). The alternation theorem states necessary and sufficient conditions on the error for optimality in the Chebyshev sense. Although the theorem does not state explicitly how to find the optimum filter, the condition is phrased in terms of type I lowpass filters, the algorithm easily generalizes. From the alternation theorem, the optimum filter Ae (e jω ) will satisfy the following set of equations: W (ω i )[ H d (e jωi ) − Ae (e jωi )] = (−1) i +1 δ ,

i = 1,2,K, ( L + 2),

(11)

where δ is the optimum error and Ae (e jω ) is given by either Eq. (3) or Eq. (7). Using Eq. (.7) for Ae (e jω ) , these equations can be written as

⎡ ⎢1 x1 ⎢ ⎢1 x 2 ⎢ ⎢M M ⎢ ⎢1 x L+2 ⎢⎣

x12

L

x1L

x 22

L

x 2L

M

M

x L2+ 2 L x LL+ 2

1 ⎤ W (ω1 ) ⎥ ⎥ −1 ⎥ W (ω 2 ) ⎥ ⎥ M L+2 ⎥ (−1) ⎥ W (ω L + 2 ) ⎥⎦

jω ⎡ a 0 ⎤ ⎡ H d (e 1 ) ⎤ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ jω 2 ⎢ ⎢ a 1 ⎥ H d (e ) ⎥ ⎥, ⎢ ⎥ =⎢ M ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢M ⎥ ⎢ ⎥ ⎢ jω ⎣⎢δ ⎦⎥ ⎢⎣ H d (e L + 2 )⎥⎦

(12)

where xi = cos ω i . This set of equation serves as the basis for an iterative algorithm for finding the optimum Ae (e jω ). The procedure begins by guessing a set of alternation frequencies ω i , i = 1,2,K , ( L + 2) . Note that ω P and ω s are fixed and are members of the set of alternation frequencies. Specifically if ω l = ω p , then

ω l +1 = ω s . The set of equations (12) could be solved for the set of coefficients a k and δ . However, a more efficient alternative is to use polynomial interpolation.

131

Specifically, Parks and McClellan found that for the given set of the extremal frequencies, δ is given by the formula L+2

δ=

∑b H k =1 L+2

k

d

( e jω k ) ,

bk (−1) k +1 ∑ k =1 W (ω k )

(13)

where L+2

bk = C i =1 i≠k

1 ( x k − xi )

(14)

and, as above, xi = cos ω i . That is, if Ae (e jω ) is determined by the set of coefficients a k that satisfy Eq. (12), with δ given by Eq. (13), then the error function goes through ± δ at the ( L + 2) frequencies ω i or, equivalently, Ae (e jω ) has values 1 ± Kδ if 0 ≤ ω i ≤ ω p and ± δ if ω s ≤ ω i ≤ π . Now since Ae (e jω ) is known to be an L th-order trigonometric polynomial, can interpolate a

trigonometric polynomial through ( L + 1) of the ( L + 2) known values E (ω i ) (or equivalently Ae (e jω ) ). Parks and McClellan used the Lagrange interpolation formula to obtain L +1

jω

Ae (e ) = P(cos ω ) =

∑ [d k =1 L +1

k

∑ [d k =1

/( x − x k )]C k , k

(15)

/( x − x k )]

where x = cos ω , xi = cos ω i , C k = H d (e jωk ) −

132

(−1) k +1 δ , W (ω k )

(16)

and L +1

dk = ∏ i =1 i≠k

bk 1 . = ( x k − xi ) ( x k − x L + 2 )

(17)

Now, Ae (e jω ) is available at any desired frequency without the need to slove the set of Eq. (12) for the coefficients a k . The polynomial of Eq. (15) can be used to evaluate Ae (e jω ) and also E (ω ) on a dense set of frequencies in the passband and stopband. If E (ω ) ≤ δ for all ω in the passband and stopband, then the optimum approximation has been found. Otherwise a new set of extremal frequencies must be found. In this algorithm all the impulse response values he [n] are implicitly varied on each iteration to obtain the desired optimum approximation, but the values of he [n] are never explicitly computed.

Figure B.1 shows the flowchart of Parks-McClellan algorithm.

133

Initial guess of (L+2) extremal frequencies

Calculate the optimum δ on extremal set

Interpolate through (L+1) points to obtain Ae (e jω )

Calculate error E (ω ) and find local maxima where E (ω ) ≥ δ

More than (L+2) extrema?

Yes

Retain (L+2) largest extrema

No

Changed

Check whether the extremal points changed

Unchanged Best approximation

Figure B.1 Flowchart of Parks-McClellan algorithm.

134