826

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO.

9.

SEPTEMBER

1990

An Approach to Cardiac Arrhythmia Analysis Using Hidden Markov Models

Abstract-This paper describes a new approach to ECG arrhythmia analysis based on “hidden Markov modeling” (HMM), a technique successfully used since the mid-1970’s to model speech waveforms for automatic speech recognition. Many ventricular arrhythmias can be classified by detecting and analyzing QRS complexes and determining R-R intervals. Classification of supraventricular arrhythmias, bowever, often requires detection of the P wave in addition to the QRS complex. The bidden Markov modeling approach combines structural and statistical knowledge of the ECG signal in a single parametric model. Model parameters are estimated from training data using an iterative, maximum likelihood reestimation algorithm. Initial results suggest that this approach may provide improved supraventricular arrhythmia analysis through accurate representation of the entire beat including the P wave.

I. INTRODUCTION ARDIAC arrhythmias are alterations of cardiac rhythm that disrupt the normal synchronized contraction sequence of the heart and reduce pumping efficiency. Causes include rate variations of the cardiac pacemaker, ectopic pacemaker sites, and abnormal propagation of pacing impulses through the specialized cardiac conduction system. Type and frequency of occurrence of arrhythmias provide an important indication of the electrical stability of the heart. In particular, certain ventricular arrhythmias are thought to indicate susceptibility to life-threatening conditions. Since arrhythmias can be suppressed by anti-arrhythmic drugs, early recognition is important. Arrhythmias are identifiable in the electrocardiogram (ECG) recorded from the surface of the chest. Long-term ambulatory ECG monitoring (‘ ‘Holter” monitoring) is used for detection of infrequent arrhythmias, evaluation of patients at risk for sudden cardiac-death, monitoring patients after myocardial infarction, documenting the effectiveness of anti-arrhythmic drug therapy, and assessing the functioning of implanted cardiac pacemakers [ 11. The high volume of data and the need for timely results create the requirement for automated analysis of contin-

C

Manuscript received March 22, 1989; revised November 21, 1989. This work was supported in part by Allegheny-Singer Research Institute and by the National Institutes of Health under Grant 1 R 0 1 HL40525-01. D. A. Coast, G. G. Cano. and S. A. Briller are with Allegheny-Singer Research Institute, Pittsburgh, PA 152 12. R. M. Stern is with the Department of Electrical and Computer Engineering, Biomedical Engineering Program, Carnegie Mellon University, Pittsburgh, PA 15213. IEEE Log Number 9036838.

uous ECG’s. A typical Holter monitor records over 100 000 heartbeats in each of two channels in 24 h. Arrhythmia detection and classification is complicated by artifact from skeletal muscle, electrode movement, and power line interference. Automated systems for arrhythmia analysis have proliferated since the early 1960’s and many are used clinically [2]. These systems use only the QRS complex and the R-R interval to group arrhythmias by origin into ventricular or supraventricular categories and to further analyze ventricular arrhythmias. The QRS complex is detected with a specialized algorithm and then classified by either template-matching based on cross-correlation measures or decision theoretic techniques. No attempt is made to detect P waves, preventing complete analysis of supraventricular arrhythmias [2], [3]. One system which does detect P waves is described by Jenkins et al. [4]. An esophageal electrode with attached lead wire is swallowed by the patient and then positioned near the atria to record high-amplitude P waves. Although this approach provides improved analysis (especially for supraventricular arrhythmias), its use has been limited, probably due to the perceived discomfort of the esophageal electrode and lead wire. Another system which incorporates P-wave information for analysis is described by Gustafson and Wang [ 5 ] . This approach analyzes R-R and P-R intervals using statistical signal analysis techniques, and is an extension of an earlier system for R-R interval analysis described by Gustafson et al. [6], [7]. The method used for P-R and R-R interval determination is not specified, but several P waves were not detected and had to be manually inserted. Motivated by the need for improved performance, model-based systems employing syntactic [8] or structural [9], [lo] pattern recognition have been reported. These systems attempt to represent the entire ECG signal including P waves but difficulty in accurately extracting primitive structural waveform components has yet to be overcome. This paper considers the application of hidden Markov models to arrhythmia analysis. This technique has been successfully applied to similar problems in automated speech recognition systems since the mid- 1970’s. Hidden Markov modeling is a statistical modeling technique that characterizes an observed data sequence by a probability density function which varies according to the state of an

0018-9294/90/0900-0826$01.OO O 1990 IEEE

~~

~~

3

-

~~

--

I 827

COAST cf d . : CARDIAC ARRHYTHMIA ANALYSIS USING MARKOV MODELS

underlying Markov chain. Its principal advantage is that the Markov chain topology preserves structural characteristics while state parameters account for the probabilistic nature of the observed data. The hidden Markov modeling approach is quite different from approaches using Markov chains to model R-R interval sequences such as those described by Gersch et al. [ l l ] and Vinke et al. [12]. These approaches characterize rhythms by the probability that a given R-R interval (quantized as short, normal, or long) will follow any other R-R interval. A Markov chain with three states (representing short, normal, and long R-R intervals) is used with state-transition probabilities estimated for each type of rhythm. In contrast, hidden Markov modeling uses a Markov chain to model the sequence and duration of waveforms and intervals within each beat. The goal of this research is to accurately identify each beat by its wavefront components so that complete arrhythmia analysis can be achieved. This paper is structured as follows: Section I1 introduces the concept of hidden Markov modeling, Section I11 describes the specification and training of hidden Markov models to represent ECG signals, Section IV describes the application of these models to classify beats and presents results from evaluation of a standardized database, and Section V summarizes performance of the technique and suggests some issues for future research. 11. HIDDENMARKOVMODELING Hidden Markov models have been applied to signal analysis problems such as speech processing since the mid 1970’s. Popularity of these models is primarily due to an automatic parameter estimation algorithm discovered by Baum and Eagon [13] in the late 1960’s and subsequently applied to speech processing (e.g., [14], [15]). Refinements in the theory and implementation have greatly enhanced this technique so that it has become the dominant approach to automatic speech recognition (e.g., [ 161, ~ 7 1 [181). , This section introduces the terminology and summarizes major results of hidden Markov modeling theory. Refer to Rabiner [19] for a tutorial on hidden Markov models including detailed descriptions of the algorithms and application examples. The problem addressed by the hidden Markov modeling technique is how to infer the state sequence S = [ sf 1) t = 0, , T - 1 1 of a “hidden” random process from values of a related process, the observation sequence 0 = [ o , 11 t = 0, * . , T - 1 1 . The sequence S represents the desired information. In arrhythmia analysis, S represents the electrical activation of the heart, and 0 is the measured ECG signal. A hidden Markov model M for continuous data representation is characterized by a 4-tuple, ( Q , U , A, B) , consisting of a set of states Q = [ q , 11 i = 0, * * , N 1 1 , a vector of initial probabilities U = [a,11 i = 0, * . . , N - 1 1 , a matrix of transition probabilities A = [a,, I(i = 0, , N - 1 ; j = 0, * * , N - 1 1 , and a

-

---

vector of probabilistic output functions B = [ b, ( x ) I( i = 0, * . . , N - 1 ] where x is an observed signal value. a, = Pr ( 4 , at t = 0)

a, = Pr ( q5 at t

6, (x)

=

Pr

(0,=

+ I I q, at t)

x J q ,at t ) .

(1) (2) (3)

To apply this modeling technique, one considers the data sequence to be analyzed as the observation sequence of a hidden Markov process. It is assumed that this sequence was generated in the following manner. At time t = 0 the model starts in one of the states q, with probability a,.A random observation value x is produced with probability b, ( x ) . The model then jumps to state q, with probability a,,, produces a new output value, and so on until Toutputs have been generated. In this paper, the output functions B are assumed to be univariate Gaussian probability density functions. This assumption is based on amplitude histograms computed for each ECG waveform component using hand-segmented data. Each output function is specified by the values of its mean and standard deviation parameters, p, and U,2

(4) Thus, for each hidden Markov model, the required parameters are N means, N standard deviations, N initial state probabilities, and N 2 transition probabilities. In this paper, a left-to-right model topology is used wherein only 2N of the transition probabilities are nonzero and one state is designated as the initial state (initial probability of 1 .O). Three problems are considered fundamental in hidden Markov modeling applications: 1) Estimation of hidden Markov model parameters from a set of representative training data (parameters include state transition probabilities, output function mean vectors, and covariance matrices). 2) Efficient calculation of Pr ( 01 M )-the probability that a given observation sequence was generated by a particular hidden Markov model M. 3) Determination of S the most likely underlying state sequence corresponding to a given observation sequence 0 such that Pr ( 0I S I, M ) is maximum. The importance of solving the first problem is obvious; model parameters must first be estimated before the models can be used for classification purposes. The solution to the second problem is often used as the basis for a classification system. By computing Pr ( 0 I Mi ) for each of i possible models and choosing the most likely model, classification can be inferred. An alternative classification approach uses the solution of thc third problem to find the single best state sequence which maximizes Pr ( 0 I S ’, M ) . Classification can then be inferred by choosing the model with the most likely best state sequence, which requires less computation than determining the most likely model. Another benefit from solving the third probI,

828

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37. NO. 9. SEPTEMBER 1990

lem is that the signal undergoing analysis is automatically segmented when the best state sequence is found. The solution to the second problem is the most obvious and will be considered first.

A. Calculating Pr ( 0I M ) The probability Pr ( 0 1 M ) can be calculated directly by enumerating all possible state sequences S and summing the probability of generating the observation sequence from each of the sequences. Although this is a straightforward approach, it is not practical since the number of state sequences for an N-state model increases exponentially with the length T as NT. Instead, this problem can be solved with computational complexity linearly proportional to T by defining forward probabilities a,( i ) where

complexity) increases exponentially with the number of states. Instead, an iterative procedure known as the forward-backward algorithm [ 151, [ 171 is employed which limits computational complexity to a linear increase with the length of the sequence. This algorithm uses the recursively defined forward probabilities a, ( i ) and defines backward probabilities 0,( i ) as well. By storing these intermediate values in two N x T tables, significant computational savings are realized in calculating the expected values of statistics required to estimate model parameters. For each iteration, the forward and backward probabilities are computed and the results are used to compute new parameter estimates. A detailed description of the forward-backward algorithm can be found in Rabiner [ 191.

C. Determining the Most Likely State Sequence

The state sequence with the highest probability of hav( 5 ) ing generated a particular observation sequence could be determined by computing the probability for each possible state sequence and choosing the maximum. However, N- I computation increases exponentially with the length of the ~ ~ t ( =i ) ,E a , - l ( j ) ajibi(o,) J=o sequence. Using the Viterbi algorithm [21], [22], this fort=1,2;-.,T-l. ( 6 ) computation increases only linearly with the sequence length. This algorithm finds the best path through a stateThe total probability Pr ( 0 I M ) is easily computed as time trellis representing all possible model states for every N- I element of the observation sequence. For a complete dePr(OlM) = j C = O ~ ~ - ~ ( j ) . (7) scription of this algorithm, refer to Rabiner [ 191.

-

, 0,and s, = qi I M ) and then computing them using the recursion a , ( i )= Pr (oo,

Direct computation of Pr ( 0 I M ) requires on the order of 2TN calculations whereas the forward computation uses only on the order of N 2 T calculations (assuming all elements of A are nonzero).

B. Estimating Model Parameters Model parameters could be directly estimated using maximum likelihood formulas if the underlying state sequence S which generated the given training data 0 were known. The transition probability estimates a,j would be the ratio of the number of times the pair s,sJ occurred in S to the number of times s, appeared. The output function parameters would be estimated using sample mean and sample covariance formulas [20] applied to the observation values generated from each state. Unfortunately, the underlying state sequence S which is assumed to have occurred in generating the training data is not known. Using knowledge of the observation sequence 0, however, and assuming an initial set of parameter values, the expected values of statistics required for parameter estimation formulas can be computed. The model parameter estimates are then computed by substituting the expected values for the unknown values in the maximum likelihood estimation formulas. This approach, first described by Baum and Eagon [ 131 and often called Baum-Welch reestimation, produces new parameter estimates which have an equal or greater likelihood of having generated the training data. Direct implementation of this approach is not practical since the number of state sequences (and computational

D. Application to Arrhythmia Analysis This research formulates the arrhythmia analysis problem as one of decoding the observation sequence of a hidden Markov process (Fig. 1). The objective is to recover an unobservable Markov state sequence representing the electrical activation of the heart from the observed ECG signal. Each waveform or interval in the ECG signal is assumed to correspond to a state of this Markov process. The observation sequence elements are digital samples of the ECG signal or derived from digital ECG samples. In Section 111, procedures necessary for creating and training hidden Markov models are discussed, and Section IV describes the application of these models for ECG analysis. 111. MODELING ECG WAVEFORMS This section describes the method by which ECG signals are represented in terms of hidden Markov models. Preliminary investigations [23] established a suitable model structure and revealed the need for evaluation of alternate observation sequence representations and for development of a training procedure to obtain accurate parameter estimates. A. Model Structure

Each class of beat is represented by a separate model such as the one shown in Fig. 2. One state is assigned to each waveform or interval of the ECG waveforms; states are connected in a left-to-right order. This topology re-

COAST e!

U/.:

829

CARDIAC ARRHYTHMIA ANALYSIS USING MARKOV MODELS

Models Input ECG

-HQE Processing Signal

Decoder

Recognized Beats

Fig. 1 . Hidden Markov modeling applied to ECG analysis.

P

f,O

PR

R

s

ST

T

f40

PR

Fig. 2 . Hidden Markov modeling of isolated normal beats. Left panel: Markov chain topology. Right panel: the top trace shows the beat analyzed, the bottom trace shows the best path through the state-time trellis derived using the Viterbi algorithm. The best path of segmentation is represented as a step-like waveform. with each level representing a different model state.

f

Q.6.Q.Q.Q

\\

Supraventricular

ti

(si\

I

Ventricular cycle (v)

V

l

t

i

I

s

ti

I

1

/ I

normal cycle (H)

\

N

l

/

>tares

T

Fig. 3. Combined hidden Markov model (three branches) for continuous ECG analysis. Left panel: Markov chain topology. Right panel: continuous ECG data with normal beats ( N ), a ventricular arrhythmia ( V ) and a supraventricular arrhythmia (S) and the corresponding best state-time path.

flects the sequential activity of the cardiac conduction system. Using the hidden Markov modeling approach, one observation sequence element (digitized sample of the ECG signal) is generated for each state transition so there is a one-to-one registration between the observation sequence and the underlying state sequence. Transitions from each state back to itself allow repetitions of model state symbols so that the variable duration of each waveform segment can be represented. The right panel in Fig. 2 illustrates the representation accuracy of the hidden Markov modeling analysis by comparing the beat (top trace) with the automatically derived state sequence represented as a step-like waveform (bottom trace)-the levels of the waveform steps have no significance other than to indicate different model states. This sequence corresponds to the best path through the state-time trellis representing possible correspondence be-

-

1'

-___

tween the model states and the ECG data. Note that all of the waveform events ( P wave, R wave, etc.) are correctly identified and accurately located. A combined model is used to detect and classify several beat categories and is a parallel combination of individual models such as the one shown in Fig. 2. Analysis requires decoding an observed sequence of ECG values to recover the most likely path through state-time trellis of the combined model, inferring classification. Fig. 3 shows a combined hidden Markov model designed to detect beats in any of three classes. A circular transition from the common final state back to the common initial state enables analysis of continuous ECG data. For example, the ECG data shown in Fig. 3 consist of normal beats with one ventricular premature beat and one supraventricular premature beat. The lower trace of Fig. 3 shows a step-like waveform which represents the actual

830

IEEE TRANSACTIONS ON B I O M E D I C A L E N G I N E E R I N G . VOL. 37. NO. Y. SEPTEMBER I Y ~ O

state sequence derived automatically using the Viterbi algorithm. Note that all beats are accurately segmented and that absence of a P wave for the supraventricular beat is recognized.

The filtering operation is implemented using only integer arithmetic (shifts, addition, subtraction) and avoids multiplications or divisions. The normalized frequency response is shown in Fig. 4 .

B. Observation Sequence Selection

C. Training Procedure Development Ideally, an arrhythmia analysis system should be capable of automatic analysis of any new patient’s ECG recording. Such a system can be called “patient-indepenConversely, the term “patient-dependent’’ dent. describes a system requiring supervised training to analyze ECG recordings from each new patient. The system described in this paper was implemented in a patient-dependent manner to test the applicability of hidden Markov models to ECG analysis. Further research is underway to extend this approach to a patient-independent analysis system. Once the model structure is established, model parameters must be estimated. The process of estimating these involves first providing an initial estimate and then refining the estimate from a training set with the forward-backward algorithm. Characteristics of the training procedure were investigated [25] with results as summarized below. 1 ) Size of Training Set: The training set for each beat category consists of several example beats isolated from the current patient’s ECG recording by marking beginning and end points. In an experiment to determine the amount of training data necessary, accuracy of estimated parameters was compared to true parameters by computing the probabilistic model distance proposed by Juang and Rabiner [26]. This experiment showed that the amount of training data required for accurate parameter estimates was only one or two beats for very regular noise-free data. Using the forward-backward algorithm, parameter estimates converge rapidly (within five iterations) provided the initial estimates of output means are accurate. However, inaccurate initial estimates prevent the forward-backward algorithm from converging to correct parameter values, regardless of the number of iterations. This phenomenon is discussed next. 2) Accuracy of Initial Estimates: Rabiner et al. [27] have shown that good initial parameter estimates are essential for the forward-backward algorithm to reach the globally optimum parameter estimates in a speech recognition application. In preliminary experiments with hidden Markov models applied to arrhythmia analysis, the forward-backward algorithm did not produce parameter estimates which resulted in correct segmentation of the training data unless accurate initial estimates for output function means were used. An experiment to determine the severity of this problem indicated that initial estimates of output means must be within one standard deviation of the true value for the training algorithm to accurately estimate the model parameters [ 2 5 ] . A training procedure which satisfied the initialization requirements was devised. First, three examples of each

In speech recognition applications, the signal characteristics of short frames of data are frequently represented with Fourier transform coefficients, linear predictive coding (LPC) coefficients, or LPC cepstral coefficients [ 2 4 ] . Thus the observation sequence is vector-valued. These representations are inadequate for ECG analysis since important diagnostic information is contained in waveform morphology and interval timing. In ECG signal analysis, the observation sequence can be scalar-valued as in samples of a single ECG channel, or vector-valued as in samples from multiple ECG channels. Each ECG channel can be separately filtered to reduce artifact before combining the signal values to create the observation vector. The results presented in this paper are based on an observation sequence derived from a single ECG channel. The observation sequence used for the examples shown above was simply the sequence of digitized ECG sample values from a single channel. A significant problem with using the ECG samples directly is that ambulatory recordings are subject to low-frequency artifact, causing signal characteristics (especially the mean value) associated with each waveform segment to vary greatly from beat to beat. Other problematic sources of artifact include muscle noise and power-line interference. In preliminary experiments, several digital filters designed to minimize artifact were compared. These filters were applied first to a single ECG channel to form scalar observation sequences and then to both ECG channels to form vector observations sequences. Based on these initial experiments, a relatively simple digital filter applied to a single ECG channel was selected [ 2 5 ] .The filter consists of two stages: 1) a derivative approximation; and 2) a digital low-pass filter. The derivative approximation is implemented as a two-point central difference ( x [ n ] x [ n - M 1 , M = 6 ) . The low-pass filter is created by cascading two linear-phase moving average filters ( N = 9 ) , each with the transfer function N- I

H(z) =

c z-“

n=O

=

~

1 -2-N 1 - z-I‘

(8)

These filters are implemented in the recursive fashion implied by the second equality. The difference equation of the derivative approximation combined with the low-pass filter is y [ n ] = 2 y [ n - 1 1 - y [ n - 21 + x [ n ] - x [ n - M ]

-2(x[n

-

N]

-

x[n - N -

+ x [ n - 2 N ] - x [ n - 2N

-

MI) MI.

(9)

To compensate for the gain of (9), the resulting values are divided by 32 (implemented by right-shifting five bits).

”

83 I

COAST et al.: CARDIAC ARRHYTHMIA ANALYSIS USING MARKOV MODELS

I

10

I

Fig. 4. Frequency response of the ECG preprocessing filter. The preprocessing consists of a two-point central-difference (TPCD) with a separation of + / - 3 points followed by two cascaded nine-point moving average (MA) filters.

beat type are manually segmented; then the output mean values are initialized as the sample means of the associated segments; finally, at least five iterations of the forward-backward algorithm are used to obtain the final model parameter estimates. Although it may seem that there is no need to use the forward-backward algorithm if accurate initial estimates are obtained ahead of time, the final parameter estimates do improve beat recognition accuracy by approximately 5 % relative to the initial estimates. IV. CLASSIFICATION OF ECG WAVEFORMS The previous section described the development of a training procedure for estimating the parameters of hidden Markov models. In this section, a procedure for analyzing ECG signals usicg hidden Markov models is developed. The development entails repeated evaluation of performance and modifications of the basic Viterbi decoding algorithm to achieve satisfactory arrhythmia detection accuracy. At this point, arrhythmia detection performance on an independent set of test data is reported, followed by an evaluation of P-wave detection accuracy.

detection, T P ) . Otherwise, the candidate beat is declared a missed detection (false negative detection, FN) or an extra detection (false positive detection, F P ) . Two statistics are defined to measure algorithm performance for both QRS detection and ventricular ectopic beat (VEB) detection: 1) Sensitivity ( S e ) is the fraction of real events that are correctly detected: TD

Se =

11

(TP

+ FN)

+

2) Positive predictivity ( P ) is the fraction of detections that are real events:

+P =

TP (TP

+ FP)’

Together, these provide a means to comparatively assess accuracy. The typical measures of sensitivity and bias from classical signal detection theory are not applicable because the number of true negative detections ( T N ) required to determine these measures cannot be meaningfully tabulated for QRS detections.

A . Performance Evaluation E . Analysis Procedure Development A quantitative performance evaluation requires a dataAn analysis procedure specifically adapted to the rebase of ECG recordings with hand-annotated beat classiquirements of ECG arrhythmia analysis was developed fications. The American Heart Association (AHA) venusing six ECG recordings from the AHA ventricular artricular arrhythmia database, developed for standardized evaluation of ventricular arrhythmia detectors, was used rhythmia database. Hidden Markov models representing the six database [28]. This database consists of 80 1 /2-hour 2-channel ECG recordings digitized at 250 Hz. Each beat is identi- files were created, initialized, and trained following the fied in an annotation file prepared by expert cardiologist- procedure outlined in Section 11. The model structure for each is illustrated in Fig. 5 ; individual branch models are annotators. A procedure recommended by AAMI is used to evalu- represented by boxed classification labels. The labels are ate HMM arrhythmia detection performance [29]. This N for normal beats, V for ventricular ectopic beats, and F procedure requires a beat-by-beat matching (within 150 for “fusion” beats. msec) between beat labels from an annotated database and ‘Initially, only a sample set was available which consisted of 1201-1207, those produced by the algorithm under test (true positive excluding 1205.

’

832

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. 9. SEPTEMBER 1990

Fig. 5 . Hidden Markov model topologies used to represent the six AHA database files. Each model branch, indicated by the boxed symbols, represents a different beat morphology. The number of branches combined to form each model depends on the number of different beat morphologies present in the database file. Note that 2201 and 4201 have identical topologies as do 6201 and 7201 (but not identical model parameters).

The Viterbi decoding algorithm was used for ECG signal analysis. The 35-min database files were analyzed in increments of 20 s due to memory limitations. Circular hidden Markov models were used to analyze the continuous sequence of beats within each buffer. Performance results were computed for each file in the subset and the combined results are listed as the first row in Table I . While QRS Se scores were above 99%, QRS P scores were rather poor at 80.85 % . Closer examination revealed that the poor QRS + P scores were due to an excessive number of false positive ( F P ) beat detections. Performance requirements necessitated incremental modifications to improve state duration modeling, to provide parameter adaptation, and to implement isolated beat analysis. 1 ) Modeling State Duration: In a hidden Markov model, the duration of a waveform or interval is modeled as a state with one transition back to itself and one or more transitions to other states. This topology leads to an exponential distribution for expected duration. For example, let aii be the probability of taking the transition from state i to state i . Then the probability that the process remains in state i for n state transitions is

+

Pr(d = n ) = ( 1 - ajj)a:-’

n

> 0.

(12)

Examination of histograms of ECG segment duration shows that an exponential duration distribution is not accurate for modeling ECG waveform durations. Truncated Gaussian probability density functions provide a better match to the observed duration histograms. Rabiner et al. [27] examine the duration modeling problem for a speech recognition task (isolated word recognition) and they propose two approaches to correct the modeling inaccuracy. In the first approach, the Viterbi algorithm recursion is modified to include an internal duration model, while the second approach applies an a posteriori weighting of the total representation probability according to the actual state durations determined. When applied to arrhythmia analysis, an internal duration model reduces false positive beat detections resulting from in-

correct state durations but requires excessive computation since each step in the recursion now includes a search over all possible durations of the previous state. The second approach requires minimal additional computation, but it does nothing to correct erroneous state sequences representing false positive detections. To circumvent these problems, a new duration modeling approach was devised which specifies minimum and maximum durations (average duration plus or minus two standard deviations) for each state. Rather than searching over all possible durations of the previous state, the Viterbi recursion tracks the current duration for each state. When choosing the best path into a state, the Current duagainst the minration of the originatingstate is and maximum duration for that state. If the duration is outside the range, the likelihood score is attenuated by multiplying the log likelihood by a constant factor of 1.5 (experimentally selected). This “max/min” duration model increases computation time by less than five percent over the basic Viterbi algorithm. The max/min duration model improved QRS + P but reduced VEB P . It eliminates many false positive ( F P ) detections caused by poor duration modeling, but new FP errors result from inaccurate duration limits. If the duration of an ECG segment falls below the minimum limit or above the maximum limit, the model state sequence becomes misaligned with the data sequence and a FP or FN detection can occur. These misalignment errors are most prevalent in files which exhibit noticeable signal variations over time which implies that adaptation of model parameters is needed. 2) Temporal Parameter Adaptation: Parameters of a hidden Markov model trained to represent a particular Holter ECG recording must adapt to the changing signal characteristics in order to maintain accurate modeling. ECG signal amplitude often varies due to changes in the electrode connection or changes in body position. Waveform intervals and durations change with variations in heart rate. A recursive adaptation scheme was devised to adjust the model parameters every 20 beats. The model parameters after adaptation period k , represented by ek, are a combination of model parameters of the previous period, 6 k - I , and cuEent estimates computed from the most recent 20 beats e k :

+

ek =

o.%& I + 0.18k

(13)

The current estimates for model parameters are computed from the decoded model state sequence provided by the Viterbi algorithm. Performance results computed after incorporating parameter adaptation demonstrated significant improvement in + P measures, but still were significantly worse than results reported for commercial Holter analysis systems. 3) Isolated-Beat Viterbi Algorithm Evaluation: Although parameter adaptation improves the accuracy of state duration limits, misalignment errors are still responsible for the majority of remaining FP errors. Errors of

I 833

COAST er al.: CARDIAC ARRHYTHMIA ANALYSIS USING MARKOV MODELS

TABLE 1 COMBINED AHA DATABASE SUBSET EVALUATION RESULTS ( I NTERMS OF SENSITIVITY. Se A N D POSITIVE + P ) FOR FOURHMM-BASEDSYSTEM IMPLEMENTATIONS(N-NORMAL.VEBPREDICTIVITY, THE BASE VENTRICULAR ECTOPIC, QRS-ALL BEATS,BOTHN A N D VEB). THEFIRSTROW REPRESENTS INCREMENTAL IMPROVEMENTS TO T H E Row IMMEDIATELY SYSTEM (BASE).SUBSEQUENT ROWSREPRESENT ABOVE.(MAXIMIN = MAXIMUM/MINIMUM DURATION MODEL.ADAPT.= PARAMETER ADAPTATION. IBA = ISOLATED-BEAT ANALYSIS). AHA Database Subset Evaluation: Incremental Improvement Analysis

N Se

N + P

VEB Se

VEB + P

QRSSe

QRS + P

Base +Max/Min + Adapt. +IBA

95.27% 75.69% 83.86% 98.79%

86.05% 97.47% 98.17% 99.84%

98.10% 98.89% 97.70% 98.42%

40.05% 25.06% 35.22% 87.96%

99.98% 98.50% 99.75% 99.87%

80.85% 91.93% 97.01 % 99.78%

this type can propagate for several beats until the model and the ECG data realign. To solve this problem, isolated-beat analysis (IBA) was implemented. This approach analyzes only a length of ECG data equal to the maximum possible beat length with the circular transition from final state to initial state removed. Since the initial and final state locations are now fixed at the beat boundaries, only one beat detection is possible which eliminates any false positive detections. The actual beat length is automatically determined by searching the final state row of the Viterbi state-time trellis for the node yielding the maximum likelihood per length. Using this position as the initial state location for the following beat reduces the chance of misalignment between the model and the ECG data. Table I lists performance results for each successive algorithm modification. Note that each algorithm modification resulted in a significant increase in QRS P , while VEB + P increases only when using the isolated-beat algorithm. Improved duration modeling and parameter adaptation prevent FP beat detections, but only isolated-beat analysis significantly improves discrimination of normal beats from VEB's.

+

C. Performance of the HMM-Based Analysis System The arrhythmia detection results reported in Table I may not be representative of the true performance of this HMM analysis system since the same data were used for both development and testing. A test set of six additional AHA tapes2 not used in developing the training or analysis procedures were analyzed. These results are presented in Table 11. The combined performance measures for QRS detection and VEB detection shows in Table I1 are only slightly lower than the IBA results in Table I. The individual file results show that 3202 is largely responsible for the decreased performance. These errors are due to inadequate adaptation of model parameters to changing heart rate in this file. A larger development set is necessary to design an adaptation scheme capable of compensating for all types of signal variability. '1202-7202, excluding 2202 because of pacemaker artifact.

TABLE I1 AHA DATABASE TESTSETDETECTION PERFORMANCE AHA Database Test-Set Evaluation QRS Detection Results

File

Total

TP

FP

FN

Se

+P

1202 3020 4202 5202 6202 7202

2596 2943 2377 2352 1954 2128

2596 2926 2376 2350 1954 2122

8 44 4 28 51 24

0 17 1 2 0 0

100.00% 99.42% 99.96% 99.91% 100.00% 99.72%

99.69% 98.52% 99.83% 99.82% 97.46% 98.88%

Combined

14350

14324

159

26

99.82%

98.90%

VEB Detection Results File

Total

TP

FP

FN

Se

+P

1202 3202 4202 5202 6202 7202

0 59 120 164 237 219

0 59 113 158 235 7

0 92 8 2 21 7

0 0 7 6 2 7

100.00% 94.17% 96.34% 99.16% 96.80%

39.07% 93.39% 98.75% 91.80% 96.80%

Combined

799

777

130

22

97.25%

85.67%

The performance results for the hidden Markov modeling analysis system can be compared to those of commercial Holter analysis systems for which AHA arrhythmia database detection statistics have been published. The QRS Se and + P measures were similar for all systems compared, so only the VEB measures will be reported. The Marquette Electronics Inc. Series 8000 system [30] had a VEB Se score of 87.9% and VEB + P score of 98.1 %, while the CardioData Corp. Mk3 system [31] had scores of 97.1 % for VEB Se and 96.8% for VEB +P . Referring to the HMM results from the last row of Table 11, the score of 97.25% for VEB Se compares favorably while the VEB + P score of 85.67% is somewhat lower. Note that the HMM results were computed using only six AHA database tapes while results for the commercial systems were computed using 59 and 60 tapes, respectively. The AHA arrhythmia database is useful for computing

834

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. 9. SEPTEMBER 1990

TABLE 111 AHA DATABASE P-WAVEDETECTION PERFORMANCE: (MEAN= MEANABSOLUTE ERROR I N DETECTED IN MS. SD = STANDARD DEVIATION OF ABSOLUTE ERROR IN DETECTED P-WAVE P-WAVE LOCATION I N MS) LOCATION AHA Database Subset P-wave Detection Evaluation File

Total

TP

FP

FN

Se

+P

Mean

SD

1201 2201 320 I 420 I 620 I 720 1

292 455 389 254 449 428

2 84 449 330 254 447 425

7 5 42 1 3 2

8 6 59 0 2 3

97.26% 98.68% 84.83% 100.00% 99.55% 99.30%

97.59% 98.90% 88.71% 99.61% 99.33% 99.53%

8 I1 13 9 13 IO

7 13 13 IO

Combined

2267

2189

60

78

96.56%

97.33%

11

IO

beat detection statistics only. For comprehensive evaluation of P-wave detection, a database which includes handlabeled P-wave locations is needed. All P waves in the first 5 min of the six AHA database subset ECG segments were annotated and their positions recorded after analysis of the data using the hidden Markov modeling system was completed. Sensitivity and positive predictivity measures are used to quantify the number of P waves detected and missed. Standard deviation of the detected P-wave location error is computed to verify accuracy. For reference, P-wave location errors for three different annotators were compared and the average standard deviation between annotators was found to be 8 ms. Table I11 presents the results of P-wave detection evaluation for the AHA database subset. The individual sensitivity and positive predictivity scores for each file are all above 97 % with the exception of 3201. The mean absolute P-wave location error ranges from 8 to 13 ms with a combined average of 11 ms. Since a typical PR interval measures 120 to 200 ms, this implies at most a 10% error in average detected PR interval. The standard deviation of the absolute P-wave location error ranges from 7 to 13 ms with a combined average of 10 ms. This performance nearly matches results for human annotators who had an 8 ms standard deviation in P-wave location error. These P-wave detection results demonstrate that hidden Markov modeling is capable of accurately detecting the occurrence and location of P waves. V. DISCUSSION Initial results for an arrhythmia analysis system based on hidden Markov modeling show ventricular arrhythmia detection performance approaching that of commercial analysis systems. Beyond ventricular arrhythmias, hidden Markov modeling addresses the problem of detecting lowamplitude P waves in typical ambulatory recordings. P-wave detection results from a subset of the AHA database support continued development of the hidden Markov modeling approach as a means of improving supraventricular arrhythmia detection. However, other problems, such as that of coincident waveform components and noncoupled P waves, have not

10

9

been addressed. Since many supraventricular arrhythmias are characterized by P waves coincident with other waveform components, means of representing this situation within the hidden Markov modeling framework are being investigated. One possibility is to model atrial and ventricular activity separately, using parallel hidden Markov models. The system described in this paper is patient-dependent: an operator is required to hand-segment a few examples of each new type of beat before and during analysis. The results obtained have shown that hidden Markov modeling is effective for arrhythmia analysis and helped identify characteristics of the technique which must be considered in applying the technique to patient-independent arrhythmia analysis. Future research will investigate patient-independent implementations (requiring no supervision or operator interaction) of the hidden Markov modeling technique. Hidden Markov modeling allows a great deal of flexibility in the choice of observation sequence to be modeled. In this paper, a sequence of observations derived from a single ECG channel is used. This approach can be easily extended for multiple ECG channels by simply changing to a vector-valued observation sequence. No changes are necessary in the training and analysis algorithms other than replacing the univariate output functions with multivariate output functions. The computation required to implement Markov modeling is considerable, even with the use of the efficient Viterbi algorithm. For the models discussed in this paper, the actual computation time for analysis was approximately five times real time or about 2.5 h for each 35-min AHA database tape. The popularity of hidden Markov modeling in speech recognition applications has led to the development of special-purpose hardware implementations for Viterbi algorithm computation [32]-[34]. Using such hardware, an arrhythmia analysis system based on hidden Markov modeling could be operated at approximately 60 times faster than the recording speed. VI. CONCLUSION Hidden Markov modeling has been applied to cardiac arrhythmia analysis with promising results. Ventricular

I

COAST et u l . : CARDIAC ARRHYTHMIA ANALYSIS USING MARKOV MODELS

arrhythmia detection statistics approaching those of commercial Holter analysis systems have been presented. More significantly, P-wave detection results show accurate detection of the low amplitude P wave in ambulatory ECG recordings. No currently available commercial system provides such detection of P waves.

REFERENCES [ I ] R. G . Mark and K. L. Ripley, “Ambulatory ECG monitoring: Realtime analysis versus tape scanning systems,’’ MD Computing, vol. 2, no. I , pp. 38-50, 1985. 121 C. L. Feldman, “Computer detection of cardiac arrhythmias: Historical review.” Amer. Rev. Diag., vol. 2, no. 5 , pp. 138-145, 1983. [3] R. G . Mark, “Arrhythmia monitoring: Current status and future challenges,” in Patient Monitoring and Data Management, pp. 7-1 I . AAMI, 1985. AAMI Technology Analysis and Review: TAR No. 11-85. [4] J . M. Jenkins, D. Wu, and R. C . Arzbaecher, “Computer diagnosis of supraventricular and ventricular arrhythmias,” Circdation, vol. 60, no. 5 , pp. 977-985, 1979. 151 D. E. Gustafson and J.-Y. Wang. “Cardiac rhythm interpretation using statistical P and R wave analysis,” in Frontiers ofEngineering in Health Cure. New York: IEEE, 1981, pp. 267-272. 16) D. E. Gustafson, A. W. Willsky, J.-Y. Wang, M. C . Lancaster, and J . H. Triebwasser, “ECGIVCG rhythm diagnosis using statistical signal analysis-I. Identification of persistent rhythms.” IEEE Truns. Biomed. Eng.. vol. BME-25, no. 4 , pp. 344-353, 19781 [7] -- “ECGIVCG rhythm diagnosis using statistical signal analysis11. Identification of transient rhythms,” IEEE Trans. Biomed. Eng., vol. BME-25, no. 4 , pp. 353-361, 1978. 181 E. Skordalakis, “Syntactic ECG processing: A review,” Pattern Recogn., vol. 19, no. 4, pp. 305-313, 1986. [9] G . C . Stockman and L. N. Kanal, “Problem reduction representation for the linguistic analysis of waveforms,” IEEE Trans. Pattern A n d . Much. Intell., vol. PAMI-5, pp. 287-298, May 1983. [IO] F. L. Xiong, B. A. Lambird, and L. N. Kanal, “An experiment in recognition of ekctrocard,iogrdms using a structural analysis algorithm,” in Proc. IEEE Con$ Cybernet. Soc., Bombay and New Delhi, India, IEEE SMC Society, 1984. [ 111 W. Gersch. P . Lilly. and E. Dong, J r . . “PVC detection by the heartbeat interval data-Markov chain approach,” Compur. Biomed. Res., vol. 8 , pp. 370-378. 1975. [ 121 U .V. H. Vinke et a l . , “Classification of cardiac rhythms using theory of Markov chains,” in Computers in Cardiology, Vol. 6, New York: IEEE, 1979, pp. 255-258. 1131 L. E. Baum and J . A. Eagon, “An inequality with applications to statistical prediction for functions of Markov processes and to a model for ecology,” Bull. Amer. Math. Soc., vol. 73. pp. 360-363, 1967. [ 141 I . K. Baker, Stochastic Modeling f o r Automatic Speech Understanding. New York: Academic. pp. 521-542. 1975. [ 151 F. Jelinek, “Continuous speech recognition by statistical methods,” Proc. IEEE, vol. 64, pp. 532-556, Apr. 1976. 1161 L. R. Bahl, F. Jelinek, and U. L. Mercer, “A maximum likelihood approach to continuous speech recognition,” IEEE Trans. Partern Anal. Mach. Intell., vol. PAMI-5, no. 2, pp. 179-190, 1983. 117) L. R. Rabiner and B.-H. Juang, “An introduction to hidden Markov models,” IEEE Acoust., Speech, Signal Process. Mag., pp, 4-16, Jan. 1986. 1181 K.-F. Lee. “Large-vocabulary speaker-independent continuous speech recognition: the SPHINX system,” Ph.D. thesis. Dep. Comput. Sci., Carnegie Mellon Univ., Pittsburgh, PA, 1988. [I91 L. R. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition,” Proc. IEEE, vol. 77, pp. 257286, Feb. 1989. 1201 R. 0. Duda and P. E. Hart, Puttern Classificution ond Scene Antrlysis. New York: Wiley. 1973. 1211 A. J. Viterbi, “Error bounds for convolutional codes and an asymptotically optimal decoding algorithm,” IEEE Trans. Inform. Theory, vol. IT-13, pp. 260-269, Apr. 1967. 1221 G . D. Forney. Jr.. “The Viterbi algorithm.” Proc. IEEE. vol. 61. pp. 268-278, Mar. 1973. 1231 D. A. Coast, R. M. Stern, G . G . Cano, and S. A. Briller, “Cardiac

835

arrhythmia analysis using hidden Markov models,” in Proc. Ninth Annu. Con$., IEEE Eng. Med. Biol. Soc., Nov. 1987, pp. 901-903. 1241 L. R. Rabiner, K. C. Pan. and K. K. Soong, “On the performance of isolated word recognizers using vector quantization and temporal energy contours,” AT&TBell Lab. Tech. J . , vol. 63, pp. 1245-1260, Sept. 1984. [25] D. A. Coast, “Cardiac arrhythmia analysis using hidden Markov models,” Ph.D. thesis, Dep. Elect. Comput. Eng., Carnegie Mellon Univ., Pittsburgh, PA, 1988. [26] B.-H. h a n g and L. U. Rabiner. “A probabilistic distance measure for hidden Markov models,” AT&T Tech. J . , vol. 64, no. 2, pp. 391408, 1985. [27] L. R. Rabiner, B.-H. Juang, S. E. Levinson, and M. M. Sondhi. “Some properties of continuous hidden Markov model representations,” AT&T Tech. J . , vol. 64, pp. 1251-1269, July-Aug. 1985. [28] R. E. Hermes, D. B . Geselowitz, and G . C. Oliver, “Development distribution and use of the AHA data base for ventricular arrhythmia detector evaluation,” in Computers in Curdiology. New York: IEEE, 1980, pp. 263-266. [29] Association for the Advancement of Medical Instrumentation (AAMI). “Recommended practice for testing and reporting performance results of ventricular arrhythmia detection algorithms,” Tech. Rep. AAMI ECAR-1987, AAMI, Arlington, VA, 1987. 1301 P. Schluter and D. Clapham et al., “The design and evaluation of a computer based system for Holter tape analysis.” in Computers in Cardiology. New York: IEEE, 1984, pp. 193-196. [3 I ] C. L. Feldman, Computer-Based Holter Scanners-Current Stutus. Northboro, MA: Cardiodata. 1984. [32] T. S . Anantharaman and R. Bisiani, “A hardware accelerator for speech recognition algorithms,” in 13th Internat. Symp. Conipur. Architect., Tokyo, Japan, 1986. [33] S . C. Glinski, T. M. Lalumia, D. R. Cassiday, T . Koh, C. Gerveshi. G . A. Wilson, and J. Kumar. “The graph search machine (GSM): A VLSI architecture for connected speech recognition and other applications,” Proc. IEEE. vol. 75, pp. 1172-1184, Sept. 1987. [34] R. Bisiani, T. Anantharaman, and L. Butcher. “BEAM: An accelerator for speech recognition,” in Proc. Internut. Conf. Acoustics, Speech. and Signul Process., Glasgow, Scotland, pp. 217-220. IEEE. May 1989, pp. 217-220.

Douglas A. Coast (S’81-M’87) was born in Harrisville, PA, in 1960. He received the B.S., M . S . , and Ph.D. degrees in electrical-biomedical engineering from Carnegie Mellon University, Pittsburgh. PA in 1982. 1984. and 1988 respectively. He joined Allegheny Singer Research Institute, Pittsburgh, PA, in 1987 where he now holds the position of Associate Research Scientist. His research interests include biomedical signal processing, pattern recognition. and bioinstrumentation.

Richard M. Stern (M’77) was born on July 5, 1948, in New York City. He received the S . B . degree from the Massachusetts Institute of Technology. in 1970, the M.S. degree from the University of California, Berkeley, in 1972, and the Ph.D. degree from MIT in 1977, all in electrical engineering. He has been on the faculty of Carnegie Mellon University since 1977. where he is currently an Associate Professor of Electrical Engineering. In addition to biomedical signal analysis, he is actively involved in research in automatic speech recognition and auditory perception. His speech research is concerned with the development of algorithms to improve the robustness of speech recognition systems in the presence of noise and reverberation, and he has previously,wwked on problems in sentence parsing, speaker adaptation, and statistical classification procedures. His research in auditory perception has been primarily con-

836

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 31. NO. 9, SEPTEMBER 1990

cerned with the development of models to describe and predict the results of psychoacoustical studies of binaural lateralization and related phenomena from the response of the auditory periphery to sound. Dr. Stem is a member of the Acoustical Society of America and the Audio Engineering Society.

Gerald G. Can0 (M’85) was born in Pittsburgh, PA in 1946. He received the B.S.E.E. and M.S.E.E. degrees from the University of Notre Dame, Notre Dame, IN, in 1968 and 1970, respectively, and the Ph.D. E.E.-B.M.E. degrees from Carnegie Mellon University, Pittsburgh, PA, in 1980. From 1980 to 1983 he was a bioengineer with the Heart Station, Allegheny General Hospital, Pittsburgh. Currently he is Senior Scientist and Director, Cardiology-Medical Electronics Lab,

Allegheny-Singer Research Institute, Pittsburgh. He is also Senior Lecturer, Biomedical Engineering Program, Carnegie-Melion University. His research interests include modeling, signal processing and microprocessorbased medical instrumentation.

Stanley A. Briller received the B.S. degree from Yale University, New Haven, CT, in 1945, the M.D. degree from New York University College of Medicine in 1947. and the M.A. (Hon) degree from the University of Pennsylvania, Philadelphia, in 1971. He is currently Professor of Medicine at the Medical College of Pennsylvania, Professor of Engineering in Medicine at Carnegie-Mellon University, Clinical Professor of Medicine at the University of Pittsburgh and directs the Heart Station at Allegheny General Hospital, Pittsburgh, PA. His research interests are centered about fundamental electrocardiography.