Repositorio Institucional de Documentos

2012 14 Mariano Llamedo Soria Signal processing for automatic heartbeat classification and patient adaptation in the electrocardiogram Departamento Director/es Instituto de Investigación en Ingeniería [I3A] Martínez Cortés, Juan Pablo Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Mariano Llamedo Soria SIGNAL PROCESSING FOR AUTOMATIC HEARTBEAT CLASSIFICATION AND PATIENT ADAPTATION IN THE ELECTROCARDIOGRAM Director/es Instituto de Investigación en Ingeniería [I3A] Martínez Cortés, Juan Pablo Tesis Doctoral Autor 2012 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA en Instituto Universitario de Investigación en Ingeniería de Aragón PhD Thesis Signal Processing for Automatic Heartbeat Classification and Patient Adaptation in the Electrocardiogram Mariano Llamedo Soria Advisor Juan Pablo Martínez Cortés, PhD Zaragoza, June 2012 ii Abstract Cardiovascular diseases are currently the biggest single cause of death in developed countries, so the development of better diagnostic methodologies could improve the health of many people. Arrhythmias are related to the sudden cardiac death, one of the challenges for the modern cardiology. On the other hand, the classification of heartbeats on the electrocardiogram (ECG) is an important analysis previous to the study of arrhythmias. The automation of heartbeat classification could improve the diagnostic quality of arrhythmias, specially in Holter or long-term recordings. The objective of this thesis is the study of the methodologies for the classification of heartbeats on the ECG. First we developed and validated a simple heartbeat classifier based on features selected with the focus on an improved generalization capability. We considered features from the RR interval (distance between two consecutive heartbeats) series, as well as features computed from the ECG samples and from scales of the wavelet transform, at both available leads. The classification performance and generalization were studied using publicly available databases: the MIT-BIH Arrhythmia, the MIT-BIH Supraventricular Arrhythmia and the St. Petersburg Institute of Cardiological Technics (INCART) databases. The Association for the Advancement of Medical Instrumentation (AAMI) recommendations for class labeling and results presentation were followed. A floating feature selection algorithm was used to obtain the best performing and generalizing models in the training and validation sets for different search configurations. The best model found comprehends 8 features, was trained in a partition of the MIT-BIH Arrhythmia, and was evaluated in a completely disjoint partition of the same database. The results obtained were: global accuracy (A) of 93%; for normal beats, sensitivity (S) 95%, positive predictive value (P+) 98%; for supraventricular beats, S77%, P+39%; for ventricular beats S81%, P+87%. In order to test the generalization capability, performance was also evaluated in the INCART, with results comparable to those obtained in the test set. This classifier model has fewer features and performs better than other state of the art methods with results suggesting better generalization capability. With an automatic classifier developed and validated, we evaluated two improvements. One, to adapt the classifier to ECG recordings of an arbitrary number of leads, or multilead extension. The second improvement was to improve the classifier with a nonlinear multilayer perceptron (MLP). For the multilead extension, we studied the improvement in heartbeat classification achieved by including information from multilead ECG recordiii iv ings in the previously developed and validated classification model. This model includes features from the RR interval series and morphology descriptors for each lead calculated from the wavelet transform. The experiments were carried out in the INCART database, available in Physionet, and the generalization was corroborated in private and public databases. In all databases the AAMI recommendations for class labeling and results presentation were followed. Different strategies to integrate the additional information available in the 12-leads were studied. The best performing strategy consisted in performing principal components analysis to the wavelet transform of the available ECG leads. The performance indices obtained for normal beats were: S98%, P+93%; for supraventricular beats, S86%, P+91%; and for ventricular beats S90%, P+90%. The generalization capability of the chosen strategy was confirmed by applying the classifier to other databases with different number of leads with comparable results. In conclusion, the performance of the reference two-lead classifier was improved by taking into account additional information from the 12-leads. The improvement of the linear classifier classifier by means of a MLP was developed with a methodology similar to the one presented above. The results obtained were: Aof 89%; for normal beats, S90%, P+99%; for supraventricular beats, S83%, P+34%; for ventricular beats S87%, P+76%. Finally we studied an algorithm based on the methodologies previously described, but able to improve its performance by means of expert assistance. We presented a patientadaptable algorithm for ECG heartbeat classification, based on a previously developed automatic classifier and a clustering algorithm. Both classifier and clustering algorithms include features from the RR interval series and morphology descriptors calculated from the wavelet transform. Integrating the decisions of both classifiers, the presented algorithm can work either automatically or with several degrees of assistance. The algorithm was comprehensively evaluated in several ECG databases for comparison purposes. Even in the fully automatic mode, the algorithm slightly improved the performance figures of the original automatic classifier; just with less than 2 manually annotated heartbeats (MAHB) per recording, the algorithm obtained a mean improvement for all databases of 6.9% in A, of 6.5% in Sand of 8.9% in P+. An assistance of just 12 MAHB per recording resulted in a mean improvement of 13.1% in A, of 13.9% in Sand of 36.1% in P+. For the assisted mode the algorithm outperformed other state-of-the-art classifiers with less expert annotation effort. The results presented in this thesis represent an improvement in the field of automatic and patient-adaptable heartbeats classification on the ECG. Resumen Las enfermedades cardiovasculares son en la actualidad la mayor causa de muerte individual en los países desarrollados, por lo tanto cualquier avance en las metodologías para el diagnóstico podrían mejorar la salud de muchas personas. Dentro de las enfermedades cardiovasculares, la muerte súbita cardíaca es una de las causas de muerte más importantes, por su número y por el impacto social que provoca. Sin lugar a duda se trata uno de los grandes desafíos de la cardiología moderna. Hay evidencias para relacionar las arritmias con la muerte súbita cardíaca. Por otro lado, la clasificación de latidos en el electrocardiograma (ECG) es un análisis previo para el estudio de las arritmias. El análisis del ECG proporciona una técnica no invasiva para el estudio de la actividad del corazón en sus distintas condiciones. Particularmente los algoritmos automáticos de clasificación se focalizan en el análisis del ritmo y la morfología del ECG, y específicamente en las variaciones respecto a la normalidad. Justamente, las variaciones en el ritmo, regularidad, lugar de origen y forma de conducción de los impulsos cardíacos, se denominan arritmias. Mientras que algunas arritmias representan una amenaza inminente (Ej. fibrilación ventricular), existen otras más sutiles que pueden ser una amenaza a largo plazo sin el tratamiento adecuado. Es en estos últimos casos, que registros ECG de larga duración requieren una inspección cuidadosa, donde los algoritmos automáticos de clasificación representan una ayuda significativa en el diagnóstico. En la última década se han desarrollado algunos algoritmos de clasificación de ECG, pero solo unos pocos tienen metodologías y resultados comparables, a pesar de las recomendaciones de la AAMI para facilitar la resolución de estos problemas. De dichos métodos, algunos funcionan de manera completamente automática, mientras que otros pueden aprovechar la asistencia de un experto para mejorar su desempeño. La base de datos utilizada en todos estos trabajos ha sido la MIT-BIH de arritmias. En cuanto a las características utilizadas, los intervalos RR fueron usados por casi todos los grupos. También se utilizaron muestras del complejo QRS diezmado, o transformado mediante polinomios de Hermite, transformada de Fourier o la descomposición wavelet. Otros grupos usaron características que integran la información presente en ambas derivaciones, como el máximo del vectocardiograma del complejo QRS, o el ángulo formado en dicho punto. El objetivo de esta tesis ha sido estudiar algunas metodologías para la clasificación de latidos en el ECG. En primer lugar se estudiaron metodologías automáticas, con capacidad v vi para contemplar el análisis de un número arbitrario de derivaciones. Luego se estudió la adaptación al paciente y la posibilidad de incorporar la asistencia de un experto para mejorar el rendimiento del clasificador automático. En principio se desarrolló y validó un clasificador de latidos sencillo, que utiliza características seleccionadas en base a una buena capacidad de generalización. Se han considerado características de la serie de intervalos RR (distancia entre dos latidos consecutivos), como también otras calculadas a partir de ambas derivaciones de la señal de ECG, y escalas de su transformada wavelet. Tanto el desempeño en la clasificación como la capacidad de generalización han sido evaluados en bases de datos públicas: la MIT-BIH de arritmias, la MIT-BIH de arritmias supraventriculares y la del Instituto de Técnicas Cardiológicas de San Petersburgo (INCART). Se han seguido las recomendaciones de la Asociación para el Avance de la Instrumentación Médica (AAMI) tanto para el etiquetado de clases como para la presentación de los resultados. Para la búsqueda de características se adoptó un algoritmo de búsqueda secuencial flotante, utilizando diferentes criterios de búsqueda, para luego elegir el modelo con mejor rendimiento y capacidad de generalización en los sets de entrenamiento y validación. El mejor modelo encontrado incluye 8 características y ha sido entrenado y evaluado en particiones disjuntas de la MIT-BIH de arritmias. Todas las carácterísticas del modelo corresponden a mediciones de intervalos temporales. Esto puede explicarse debido a que los registros utilizados en los experimentos no siempre contienen las mismas derivaciones, y por lo tanto la capacidad de clasificación de aquellas características basadas en amplitudes se ve seriamente disminuida. Las primeras 4 características del modelo están claramente relacionadas a la evolución del ritmo cardíaco, mientras que las otras cuatro pueden interpretarse como mediciones alternativas de la anchura del complejo QRS, y por lo tanto morfológicas. Como resultado, el modelo obtenido tiene la ventaja evidente de un menor tamaño, lo que redunda tanto en un ahorro computacional como en una mejor estimación de los parámetros del modelo durante el entrenamiento. Como ventaja adicional, este modelo depende exclusivamente de la detección de cada latido, haciendo este clasificador especialmente útil en aquellos casos donde la delineación de las ondas del ECG no puede realizarse de manera confiable. Los resultados obtenidos en el set de evaluación han sido: exactitud global (A) de 93 %; para latidos normales, sensibilidad (S) 95 %, valor predictivo positivo (P+) 98 %; para latidos supraventriculares, S77 %, P+39 %; para latidos ventriculares S81 %, P+87 %. Para comprobar la capacidad de generalización, se evaluó el rendimiento en la INCART obteniéndose resultados comparables a los del set de evaluación. El modelo de clasificación obtenido utiliza menos características, y adicionalmente presentó mejor rendimiento y capacidad de generalización que otros representativos del estado del arte. Luego se han estudiado dos mejoras para el clasificador desarrollado en el párrafo anterior. La primera fue adaptarlo a registros ECG de un número arbitrario de derivaciones, o extensión multiderivacional. En la segunda mejora se buscó cambiar el clasificador lineal por un perceptrón multicapa no lineal (MLP). Para la extensión multiderivacional Contents Title Page i Abstract iii Resumen........................................ v Conclusiones ..................................... x Contents xiii 1 Introduction 1 1.1 Motivation.................................... 1 1.2 Background ................................... 2 1.2.1 Theheart ................................ 2 1.2.2 From the action potentials to the electrocardiogram . . . . . . . . . 4 1.2.3 Arrhythmias .............................. 12 1.2.4 Manifestation of arrhythmias on the ECG . . . . . . . . . . . . . . 17 1.3 Previousworks ................................. 22 1.4 Objective .................................... 25 1.5 OutlineoftheThesis.............................. 26 2 Materials and Methods 29 2.1 ECGDatabases................................. 29 2.1.1 AAMI class labeling recommendations . . . . . . . . . . . . . . . . 30 2.1.2 MIT-BIH Arrhythmia Database (MITBIH-AR) . . . . . . . . . . . 30 2.1.3 MIT-BIH Supraventricular Arrhythmia Database (MITBIH-SUP) . 34 2.1.4 St. Petersburg Institute of Cardiological Technics (INCART) 12lead Arrhythmia Database . . . . . . . . . . . . . . . . . . . . . . . 34 2.1.5 European ST-T Database (ESTTDB) . . . . . . . . . . . . . . . . . 35 2.1.6 The MIT-BIH ST Change Database (MITBIH-ST) . . . . . . . . . 35 2.1.7 The Long-Term ST Database (LTSTDB) . . . . . . . . . . . . . . . 36 2.1.8 American Heart Association (AHA) ECG Database . . . . . . . . . 37 2.2 Supercomputing Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3 SignalProcessing ................................ 39 xiii xiv CONTENTS 2.3.1 ECGpreprocessing........................... 39 2.3.2 WaveletTransform ........................... 41 2.3.3 PrototypeWavelet ........................... 43 2.4 Heartbeat classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.1 Classification Features . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.2 Discriminant Functions . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.4.3 Domain Handling for some Features . . . . . . . . . . . . . . . . . 55 2.4.4 OutlierRemoval............................. 58 2.4.5 Performance evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.4.6 Model Selection and Dimensionality Reduction . . . . . . . . . . . . 65 3 Automatic ECG Heartbeat Classification 69 3.1 Introduction................................... 69 3.2 Methodology .................................. 70 3.2.1 ECGDatabases............................. 70 3.2.2 ECGpreprocessing........................... 70 3.2.3 Features and Classifiers . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.2.4 ExperimentSetup............................ 72 3.3 Results...................................... 73 3.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.A DetailedResults................................. 80 4 Extensions to the Automatic Classifier 83 4.1 Introduction................................... 83 4.2 Multilead classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2.1 Material and methods . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2.1.1 Robust Covariance Matrix Computation . . . . . . . . . . 89 4.2.2 Results.................................. 91 4.2.3 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 Neural network classifier . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.3.1 FeatureSets............................... 95 4.3.2 FeatureSelection ............................ 95 4.3.3 Multi-Layer Perceptron . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.3.4 Classifier Combination . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.3.5 Results.................................. 98 4.3.6 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . 98 4.A DetailedResults.................................100 5 Patient-Adapted ECG Heartbeat Classification 107 5.1 Introduction...................................107 5.2 Methodology ..................................108 CONTENTS xv 5.2.1 ECGdatabases .............................108 5.2.2 Heartbeats classification . . . . . . . . . . . . . . . . . . . . . . . . 109 5.2.3 Automatic classifier . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.2.4 Clustering algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.2.5 Feature selection for clustering . . . . . . . . . . . . . . . . . . . . . 113 5.2.6 Performance evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.3 Results......................................116 5.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.A DetailedResults.................................124 6 Conclusions and Future Work 137 6.1 Summary ....................................137 6.2 Conclusions ...................................138 6.3 Futurework...................................140 Scientific Contributions 143 A Matlab Implementation 145 A.1 Introduction...................................145 A.2 Features .....................................145 A.3 Installation and Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 A.3.1 The power of the command-line . . . . . . . . . . . . . . . . . . . . 148 A.3.2 The power of a high performance computing cluster . . . . . . . . . 150 A.4 Acknowledgments................................152 Acronyms 153 Figures 157 Tables 163 Bibliography 167 xvi CONTENTS Chapter 1 Introduction 1.1 Motivation The World Health Organization places cardiovascular diseases (CVD) as the first single cause of death globally in the present, and forecasts the same ranking up to 2030 [World Health Organization, 2012]. These diseases affect in a higher degree to lowand middleincome countries, but in the same proportion to women and men. Specifically in Argentina and Spain, more than 30% of the deaths are caused by CVD and is by far, the first single cause of death according to the official agencies [Dirección de Estadísticas e Información en Salud, 2012, Instituto Nacional de Estadística, 2012]. A great part of the deaths caused by CVD occur suddenly, starting with a ventricular fibrillation which leads to a cardiac arrest [Bayés de Luna, 2010]. This situation is known as sudden cardiac death (SCD) and is probably the most important challenge of the modern cardiology. This disease is unusual up to the age of 35, but from there the risk of SCD increases specially during the chronic and acute phases of myocardial infarction, or other cardiopathy related to heart failure. The identification or prediction of SCD has been studied more thoroughly for those risk groups with a previous cardiac condition (cardiac arrest, genetic defects, heart failure, heart attack) than for the people in which SCD is the first manifestation. The importance of the last group is that it represents more than the 50% of people who suffer SCD. However, up to the moment, an exhaustive screening of the population is unfeasible from the technical and economical point of view. The improvement of cost-effective methodologies for the prediction of SCD received lot of attention from the scientific community in the last decades. It was studied in several works that arrhythmias are responsible of most of the cases of SCD [Bayés de Luna, 2010]. One important advance in the study of arrhythmias was the use of long-term (or Holter) recordings and the software to aid the cardiologist in the detection and diagnostic of abnormalities in the electrocardiogram (ECG). The study of arrhythmias by means of the computerized analysis of the ECG signal, is in the present a cost-effective and well established tool to analyze the heart function. The improvement of the methodologies used 1 2CHAPTER 1. INTRODUCTION in the study of arrhythmias is likely to aid cardiologists in the diagnostic and screening of SCD. In this thesis we developed and analyzed new algorithms for the classification of ECG heartbeats, which is an important analysis previous to the study of arrhythmias. 1.2 Background As this thesis is entirely focused on the analysis of the ECG signal, a brief description of its origin is included, as well as the basic concepts of cardiac electrophysiology. We will start with a selection of anatomy and physiology concepts, to subsequently inspect some mechanisms at the cellular level and their manifestation on the ECG. Our objective in the following chapters will be the design of a computer algorithm capable of classifying the concepts explained in this section. This section is based on the books [Bayés de Luna, 2010, Natale and Wazni, 2007, Guyton and Hall, 2006, Sörnmo and Laguna, 2005, Malmivuo and Plonsey, 1995], where the reader is referred for further details and references. 1.2.1 The heart The heart is an electromechanical pulsatile pump. From the anatomic point of view, as can be seen in Figure 1.1, there are two separate pumps: one at the right that pumps blood through the lungs, and one at the left that pumps blood through the peripheral organs. Each half includes a two-chamber pump composed of an atrium and a ventricle. The atrium pumps blood for the ventricle, and then the ventricles supply the main pumping force either through the pulmonary circulation, by the right ventricle, or through the peripheral circulation by the left ventricle. There are four valves to force the direction of the blood, as is shown in Figure 1.2, two located between the atria and the ventricles, and two between the ventricles and the arteries. As a periodic electromechanical pump, an electrical impulse is responsible of the mechanical activation of the muscle. Each cycle is initiated by spontaneous generation of an action potential (AP) in the sinus (or sinoatrial in Figure 1.2) node. This node is located in the superior lateral wall of the right atrium near the opening of the superior vena cava. The impulse, or AP, travels through both atria reaching the atrio-ventricular (A-V) bundle, where is delayed about 0.1 seconds. This delay allows the atria to pump blood into the ventricles. After this, the ventricles are filled and ready to be activated. This is done by a special conduction system (SCS), the right and left bundle branches of Purkinje fibers. This system propagates the impulse from the A-V node to the whole ventricular muscle very fast, allowing a synchronized activation and consequently an effective pump of the blood. This cycle is repeated up to the death of the heart. Now we will try to relate the electrical and mechanical behavior of the heart described above. The activation of the cardiac muscle composed of two phases, contraction and re- 1.2. BACKGROUND 3 Brachiocephalic artery Superior vena cava Right pulmonary arteries Brachiocephalic veins Right atrium Atrioventricular (tricuspid) valve Chordae tendineae Right ventricle Inferior vena cava Septum Left ventricle Atrioventricular (mitral) valve Semilunar valves Left atrium Left pulmonary veins Left pulmonary arteries Aorta Left subclavian artery Left common carotid artery Figure 1.1: Structure of the heart, and course of blood flow through the heart chambers and heart valves. Diagrams based on image http://en.wikipedia ... -en.svg under license CS-BY-SA. Left posterior bundle Right bundle His bundle Purkinje fibers Sinoatrial node Atrioventricular node Bachmann's bundle Figure 1.2: Course of the blood flow through the heart, and the electrical conduction system of the heart. Diagrams based on image http://commons.wikimedia ... Heart.svg under license CS-BY-SA. 4CHAPTER 1. INTRODUCTION 50 90 130 0 40 80 120 Aortic pressure Atrial pressure P Q R S TP Q R S T -0.25 0 0.5 Voltage (mV) Electrocardiogram A-V valve closes Aortic valve opens Aortic valve closes A-V valve opens Isovolumic contraction Isovolumic relaxation Rapid inflow Diastasis Atrial systole Ejection Ejection Ventricular pressure Volume (mL) Pressure (mm Hg) Ventricular volume Mechanical part Electrical part Figure 1.3: Wiggers diagram. Events of the cardiac cycle for left ventricular function, showing changes in left atrial pressure, left ventricular pressure, aortic pressure, ventricular volume, and the electrocardiogram. laxation, or in electrical terms as depolarization and repolarization. As the heart function produces an electrical field, the voltage generated can be recorded by the electrocardiograph from the surface of the body. The first wave, called with the letter P, is caused by spread of depolarization through the atria. After the electrical activation, follows the atrial contraction which causes a slight rise in the atrial pressure. About 0.16 seconds after the onset of the P wave, the QRS waves appear as a result of electrical depolarization of the ventricles. This initiates the contraction of the ventricles and causes the ventricular pressure to begin rising. Finally, the ventricular T wave in the electrocardiogram represents the stage of repolarization of the ventricles when the ventricular muscle fibers begin to relax. As can be noted in Figure 1.3, the electrical depolarization is preceded by the corresponding mechanical contraction. 1.2.2 From the action potentials to the electrocardiogram In general heart cells can be grouped in two types: the ones from the SCS and the contractile cells. The first are responsible of the generation of the electrical impulse (rhythmicity) and its conduction to the contractile cells, while the contractile cells are responsible of the pumping or mechanical function. Both cell types are responsible of the electromechanical link. In Figure 1.4 it is showed the waveforms of the voltage, or action potential, and currents measured in the cellular membrane of a contractile cell. Following the depolar- 1.2. BACKGROUND 5 Na K Cl K K KNa ATPase pump Phase 4Phase 3Phase 2Phase 1Phase 0 Ca Ito2Ito Intra cell Extra cell 1-2 03 4 4 0 3 4 1 2 15 mV 0 mV -90 mV -40 mV 0 mV -90 mV Figure 1.4: Reproduced from [Natale and Wazni, 2007]. Top panel: on left, the action potential in contractile cells, and on the right in SCS cell. Bottom panel: predominant currents during the different phases of Na-channel-dependent action potential. ization phases in the same Figure, note that when a cell receives depolarizing current, Na channels are activated resulting in a net inward current manifested as phase 0 of the AP. Phase 1 starts with the opening of a rapid outward potassium current. Phase 2 or the plateau phase of the AP is the result of an L-type Ca current that counteracts the outward K currents. With time, L-type Ca channels are inactivated and the plateau subsides. At the same time, the increase in calcium concentration acts as a trigger for release of more Ca stored in the sarcoplasmic reticulum, which in turn provides a contraction signal to the myocyte contractile elements, producing the contraction of the cell. Phase 3 is due to ‘delayed rectifier’ outward K currents. Phase 4 constitutes a steady, stable, polarized membrane due to voltage-regulated inward rectifiers. Compared to atrial action potential, ventricular AP has a longer duration, a higher phase 2, a shorter phase 3, and more negative phase 4. On the other hand, the SCS cells have the ability to generate a spontaneous action potential using T-type Ca and K rectifier currents. These currents confer the unstable electrical property of phase 4, causing these cells to develop rhythmic spontaneous slow diastolic depolarization. Once AP reaches –40 mV, L-type Ca channels are activated, generating the slow upstroke of the action potential in these types of cells (phase 0). There are three types of SCS cells: 1. P cells, found mostly in the sinus node are responsible of automaticity. 2. The Purkinje cells, are found in the His bundle branches and are responsible of the fast transmission of electrical impulses through the ventricles. 3. The transitional cells, with slow conduction velocity, are typically found between 6CHAPTER 1. INTRODUCTION 15 mV 0 mV -90 mV ARP RRP TRP Normal AP Aberrated AP Not propagated AP Figure 1.5: Based on Figure 2.20 from [Bayés de Luna, 2010]. Refractory period of ventricular cells. During absolute refractory period (ARP) depolarization is not possible. During the relative refractory period (RRP), an increased activation is necessary to depolarize the cell. After the total refractory period, the cell is able to produce a normal AP upon activation. the P, Purkinje and contractile cells. Once any cell is depolarized it takes certain time until it can be normally depolarized again. This time is known as total refractory period (TRP). Also there is a period of time where the cell can not be depolarized, and is known as absolute refractory period (ARP). If the time of arrival of a new activation is greater than ARP, the cell can produce an aberrated AP if the stimulus is big enough. This is known as relative refractory period (RRP). There is a small time window, between RRP and ARP in Figure 1.5, where the cell reacts to an increased activation, but the activation can not be propagated. Automaticity is an intrinsic property of all myocardial cells. In addition to the sinus node, cells with pacemaking capability in the normal heart are located in some parts of the atria and ventricles. However, the occurrence of spontaneous activity is prevented by the natural hierarchy of pacemaker function, causing these sites to be latent or subsidiary pacemakers. The spontaneous discharge rate of the sinus node normally exceeds that of all other subsidiary pacemakers. Therefore, the impulse initiated by the sinus node depolarizes and keeps the activity of subsidiary pacemaker sites depressed before they can spontaneously reach threshold. However, slowly depolarizing and previously suppressed pacemakers in the atrium, A-V node, or ventricle can become active and assume pacemaker control of the cardiac rhythm if the sinus node pacemaker becomes slow or unable to generate an impulse (e.g., secondary to depressed sinus node automaticity) or if impulses generated by the sinus node are unable to activate the subsidiary pacemaker sites (e.g., sinoatrial exit block, or A-V block). The emergence of subsidiary or latent pacemakers under such circumstances is an appropriate fail-safe mechanism, which ensures that ventricular activation is maintained. Once introduced the types of AP of the heart cells, it is possible to imagine that the electrical field which produces the ECG in the body surface, results from the integration 1.2. BACKGROUND 13 Frontal Plane P Q R S T aVR P Q R T P Q R S T aVF P Q R ST aVL Lead II Lead III P Q R S T P R T Lead I QS aVR aVL aVF P T I II III AVR AVL AVF V1 V2 V3 V4 V5 V6 Figure 1.11: Normal Vectocardiogram and the projection to the 12-lead ECG. erally slight variations. However, under normal conditions and particularly in children, it may present slight to moderate changes dependent on the phases of respiration, with the heart rate increasing with inspiration. In adults at rest the rate of the normal sinus rhythm ranges from 60 to 100 beats per minute (bpm). Thus, sinus rhythms over 100 bpm (sinus tachycardia) and those under 60 bpm (sinus bradycardia) may be considered arrhythmias. However, it should be taken into account that sinus rhythm varies throughout a 24-h period and sinus tachycardia and sinus bradycardia usually are a physiologic response to certain sympathetic (exercise, stress) or vagal (rest, sleep) stimuli. Under such circumstances, the presence of these heart rates should be considered normal. The term arrhythmia does not mean rhythm irregularity, as regular arrhythmias can occur often with absolute stability (flutter, paroxysmal tachycardia, etc.), sometimes presenting heart rates in the normal range. On the other hand, some irregular rhythms should not be considered arrhythmias (mild to moderate irregularity in the sinus discharge, particularly when linked to respiration). Moreover, a diagnosis of arrhythmia in itself does not mean evident pathology. In fact, in healthy subjects, the sporadic presence of certain arrhythmias both active (premature complexes) and passive (escape complexes, certain degree of A-V block, evident sinus arrhythmia, etc.) is frequently observed. There are different ways to classify cardiac arrhythmias: •According to the site of origin: arrhythmias are divided into supraventricular (including those having their origin in the sinus node, the atria, and the AV junction) and ventricular arrhythmias. •According to the underlying mechanism: arrhythmias may be explained by: 1) abnormal formation of impulses, which includes increased heart automaticity (extra systolic or parasystolic mechanism) and triggered electrical activity, 2) reentry of different types, and 3) decreased automaticity and/or disturbances of conduction. 14 CHAPTER 1. INTRODUCTION •From the clinical point of view: arrhythmias may be paroxysmal, incessant or permanent. In reference to tachyarrhythmias (an example of an active arrhythmia), paroxysmal tachyarrhythmias occur suddenly and usually disappear spontaneously (i.e. A-V junctional reentrant paroxysmal tachycardia). Permanent tachyarrhythmias are always present (i.e. chronic atrial fibrillation), and incessant tachyarrhythmias are characterized by short and repetitive runs of supraventricular or ventricular tachycardia. •Finally, from an electrocardiographic point of view, arrhythmias may be divided into two different types: active and passive. –Active arrhythmias, due to increased automaticity, reentry, or triggered electrical activity (these mechanisms are explained below), generate isolated or repetitive premature complexes on the ECG, which occur before the cadence of the regular sinus rhythm. The isolated premature complexes may be originated in a parasystolic or extrasystolic ectopic focus. The extra systolic mechanism presents a fixed coupling interval, whereas the para systolic presents a varied coupling interval. Premature complexes of supraventricular origin are generally followed by a narrow QRS complex, although they may be wide if conducted with aberrancy. The ectopic P wave is often not easily seen as it may be hidden in the preceding T wave. In other cases the premature atrial impulse remains blocked in the AV junction, initiating a pause instead of a premature QRS complex. The premature complexes of ventricular origin are not preceded by an ectopic P wave, and the QRS complex is always wide (>120 ms), unless they originate in the upper part of the intraventricular SCS (ISCS). Premature and repetitive complexes include all types of supraventricular or ventricular tachyarrhythmias (tachycardias, fibrillation, flutter). In active cardiac arrhythmias due to reentrant mechanisms, a unidirectional block exists in some part of the circuit. –Passive arrhythmias occur when cardiac stimuli formation and/or conduction are below the range of normality due to a depression of the automatism and/or a stimulus conduction block in the atria, the AV junction, or the ISCS. From an electrocardiographic point of view, many passive cardiac arrhythmias present isolated late complexes (escape complexes) and, if repetitive, slower than expected heart rate (bradyarrhythmia). Even in the absence of bradyarrhythmia, some type of conduction delay or block in some place of the SCS may exist, for example, first-degree or some second-degree sinoatrial or A-V blocks, or atrial or ventricular (bundle branch) blocks. The latter encompasses the aberrant conduction phenomenon. Thus, the electrocardiographic diagnosis of passive cardiac arrhythmia can be made because it may be demonstrated that the ECG 1.2. BACKGROUND 15 changes are due to a depression of automatism and/or conduction in some part of the SCS, without this manifesting in the ECG as a premature complex, as it does in reentry (see Figure 1.12). The mechanisms of cardiac arrhythmias are often the results of many factors including fluctuation in intracellular concentration of Ca, after depolarization currents, refractory period shortening or lengthening, autonomic nervous system innervation, repolarization dispersion, and changes in excitability and conduction. For example, bradyarrhythmia is often caused by abnormalities in excitability. This could be caused by dysfunction in the Na channels or by ischemia-induced elevation in extracellular K concentration. Furthermore, inherent or metabolically induced abnormalities in Na channels, Ca channels, or connexin have been shown to play a role in conduction diseases. Mechanisms of tachyarrhythmias can be grouped into three categories: re-entry, triggered activity and automaticity. •Re-entry is a depolarizing wave traveling through a closed path. There are three prerequisites for re-entry: 1) At least two pathways: slow and fast AV nodal pathways, accessory pathway or the presence of barrier (anatomic: tricuspid valve; pathologic: incisional scars, myocardial infarction, and functional scar). 2) Unidirectional block: This block can be physiologic: caused by a premature complex, or increased heart rate; or pathologic: caused by changes in repolarization gradients. 3) Slow conduction to prevent collision of the head and the tail of the depolarizing wave. In functional re-entry, unidirectional block can be due to dispersion of refractoriness (repolarization) or dispersion of conduction velocity (anisotropic re-entry). See Figure 1.12 for an example of this concept. •Triggered activities are caused by after depolarization currents. They are classified as early (EAD occurring inside AP: phases 2 and 3) or delayed (DAD: phase 4). These currents can in turn be responsible for both focal and reentrant arrhythmias. The former is caused by eliciting an excitatory response exceeding the activation threshold and the latter can be developed when these currents cause prolongation in action potential which facilitates the development of a unidirectional block due to dispersion of refractoriness. •Automaticity is driven by spontaneous phase 4 depolarization. Automatic depolarizations in the atria and ventricles are not manifested normally due to overdrive suppression by the faster depolarization caused by the sinus node. However, during excess catecholaminergic states, phase 4 depolarization may exceed sinus node depolarization, causing depolarization to be driven by the abnormal tissue. Ventricular tachycardias during the acute ischemic and reperfusion phases are good examples of 16 CHAPTER 1. INTRODUCTION S L O W Ectopic focus 5 A zone of slow conduction will also increase the circuit time and allow re-entry Wave front Refractory tissue Excitable tissue Unidirectional conduction only Conduction barrier 1. Intra-atrial re-entry tends to occur around conduction barriers, especially if part of the surrounding tissue conducts in only one direction(clockwise in this example) 4. If the circuit size is larger, the circuit time increases and re-entry can occur 2. In healthy atria the depolarisation wave is likely to encounter refractory tissue when it has travelled one complete circuit 3. If the atrial refractory period is shorter than the circuit time, re-entry can occur 5. A zone of slow conduction will also increase the circuit time and allow re-entry Figure 1.12: Electrical reentry, the mechanism responsible for initiating and maintaining atrial fibrillation. Reproduced from [Grubb and Furniss, 2001]. After depolarization Plateau EAD Late EAD DAD Phase 2 Phase 3 Phase 4 Figure 1.13: Types of after depolarization currents. EAD, early after depolarization; DAD, delayed after depolarization. Reproduced from [Natale and Wazni, 2007]. 1.2. BACKGROUND 17 Sinus Tachycardia Rate 122 Normal Sinus Rhythm Rate 85 Sinus Bradycardia Rate 48 Sinus Arrhythmia V1 Figure 1.14: Several examples of sinus rhythms. automaticity. They are often originated from the border zone between normal and ischemic cells. As described above, the mechanisms that originates arrhythmias are diverse, and therefore the manifestation in the ECG. In the following section we will show the most important mechanisms as they appear in the ECG. 1.2.4 Manifestation of arrhythmias on the ECG In this subsection several examples of the mechanisms enumerated above are shown in the ECG. Normal sinus rhythm is characterized by a regular cardiac rate with normal QRS complexes whose duration must be less than 120 milliseconds, as can be seen in Figure 1.14. The P-waves are normal in shape, and are synchronized with the QRS complexes. The PR interval must be less than 0.2 seconds. Heart rates may range from 60-100 bpm. There are a number of variant types of sinus rhythm, sinus arrhythmia is a normal rhythm in which heart rate varies periodically, usually with the respiratory cycle. There is an acceleration of rate during inspiration, and a slowing of rate during expiration. Escape beats arise from lower (normally latent) pacemakers outside of the sinus node that fire because of either depressed sinus node function or blocked conduction of sinus impulses. Escape beats may originate at any pacemaker site below the sinus node. If the 18 CHAPTER 1. INTRODUCTION Carotid Pre ssure S inus Pause Atrial Escape Be at Atrial Escape Beat Figure 1.15: Example of an atrial escape beat. sinus node slows sufficiently (perhaps due to vagal tone), other latent pacemaker sites in the atrium may emerge to establish heart rate. The P-wave resulting from these beats is usually different in shape from the normal, and in many cases is inverted in polarity. This reflects the fact that the beats originate low in the atrium. Such beats are sometimes referred to as low atrial or coronary sinus beats. A-V nodal escape beats often terminate prolonged sinus pauses. The QRS complex is normal because the impulse is conducted normally to the ventricles. The P-wave is either not visible at all, or may be found just prior to or immediately following the QRS. In general the P wave is abnormal in shape since it is retrogradely conducted. If the P-wave immediately precedes the QRS complex, the beat is referred to as a fast conducted beat. Conversely, if the P-wave follows the QRS, the beat is called a slow conducted beat. Ventricular escape beats protect the heart against asystole in the event of AV block (either fixed or transitory). They are characterized by a wide and usually bizarre QRS complex. The cardiac impulse originates in the ventricular Purkinje system. It is generally conducted with a slow propagation speed (0.5 meter/second) through the myocardium, thus leading to a wide QRS complex (usually greater than 120 ms). Ventricular escape rhythms (idioventricular rhythms) are common in cases of complete heart block, and have rates of about 40 per minute. Ectopic beats could arise from pacemakers outside the sinus node as a result of an abnormal increase in rhythmicity in the ventricular Purkinje system. Atrial premature beats (APB) are seen frequently in normal individuals and have little clinical significance. They are also seen in heart disease, and when frequent, may be an early sign of atrial irritability which may progress to more serious atrial dysrhythmias. In APBs the QRS complexes are normal since they propagate normally through the ventricles via the conduction system. The P-waves are generally slightly abnormal since they originate from an abnormal focus, and propagate in an abnormal pattern. The impulse generally invades the area of the SA node and resets the sinus pacemaker. APBs occurring quite early following the previous beat may be aberrantly conducted, frequently with a right bundle branch block configuration. Aberrant conduction is particularly likely when the APB follows a long RR interval (the Ashman phenomenon). If an APB is extremely early it may run into refractory tissue in the AV node and be non-conducted. Ventricular ectopic beats (VPB) originate from somewhere in the ventricles. The QRS complex is wide (greater than 0.12 seconds) and bizarre. VPBs may exhibit fixed coupling 1.2. BACKGROUND 19 Nodal Escape Beats The last two beats are nodal escape beats which appear as sinus pacemaker slows. Nodal Rhythm in Complete AV Block 1 Ra te 2 3 4 SN Atria A-V ISCS Ventricles SN Atria A-V ISCS Ventricles Slow conducted P-wave rhythm Figure 1.16: Examples of A-V nodal escape beats. 20 CHAPTER 1. INTRODUCTION SN Atria A-V ISCS Ventricles Ventricular Escape Beat Figure 1.17: Example of a ventricular escape beat. Atrial Premature Contractions Aberrantly Conducted APBs (Ashman Phenomenon) Non-conducted (Blocked) APBs Figure 1.18: Examples of atrial premature beats. The blue triangles indicate the premature beats in the top panel, and the non-conducted beats in the bottom. 1.2. BACKGROUND 21 to previous normal beats. They may occur early or late in the cycle. The mechanism for PVCs may be reentry or triggered activity as discussed previously. Some VPBs appear to show no fixed coupling to preceding normal beats. If they show a regular rhythm of their own, they may result from a parasystolic focus. Note that some parasystolic depolarizations experience “exit block” and do not result in ventricular excitation. Parasystolic ventricular ectopic beats are usually considered relatively benign. Most VPBs are followed by a pause. The pause is usually compensatory, meaning that the coupling interval to the preceding normal beat plus the pause following the VPB comprise an interval equal to twice the normal R-R interval. An interpolated VPB is one which is sandwiched between two normal QRS complexes which arrive on time with the sinus normal activation. VPBs are often found in otherwise normal individuals and probably have little significance if they are infrequent. In heart disease, VPBs may be a risk factor for increased incidence of more serious ventricular arrhythmias and sudden death. VPBs may occur singly or in groups and the following ordering of increasing severity of ventricular ectopic activity has been proposed: 1. Occasional: less than 30 per hour VPBs of the same morphology. 2. Frequent: greater than 30 per hour uniform VPBs or bigeminy where every other beat is a VPB 3. Multiform PVCs: different QRS morphologies 4. Couplets: pairs of consecutive VPBs 5. Ventricular Tachycardia: runs of three or more VPBs 6. Ventricular Flutter: rapid ventricular tachycardia with a sinusoidal configuration caused by merging of QRSs and Ts 7. Ventricular Fibrillation chaotic electrical activity without definite QRS complexes VPBs which occur very early in the cardiac cycle such that they fall on the T-wave of the previous beat are considered particularly dangerous. At the time corresponding to the peak of the T wave, the ventricular myocardium is just beginning to repolarize. Some cells may be in the relatively refractory period, while others may be more fully recovered, and still others quite refractory. The electrical properties of the myocardium are thus quite varied, and conditions favoring reentrant loops are likely. Thus, an extra stimulus in the form of an isolated VPB which is very early-cycle may trigger a repetitive ventricular ectopic rhythm such as ventricular tachycardia or ventricular fibrillation. (The period near the T-wave peak is often referred to as the vulnerable period). Proper characterization of ventricular ectopic activity requires long-term (24-hour) ECG monitoring. The classification of heartbeats on the ECG as can be seen, is an important task for the automatic analysis of arrhythmias. This is the first task performed by a cardiologist 22 CHAPTER 1. INTRODUCTION when inspecting a recording, and as shown above, it is a very demanding task. In the next section we will review the state of the art regarding heartbeat classification algorithms. 1.3 Previous works Many algorithms for ECG heartbeats classification were developed in the last decades. Some of the most relevant before the beginning of this thesis are [Hu et al., 1997, Lagerholm et al., 2000, de Chazal et al., 2004, Inan et al., 2006, Christov et al., 2006, de Chazal and Reilly, 2006], while others were published in the last few years [Llamedo and Martínez, 2007, Jiang and Kong, 2007, Park et al., 2008, Ince et al., 2009]. However, due to the lack of standardization in the development and evaluation criteria, comparison of results across most of these works could not be performed fairly or is impossible. In order to overcome this problem, some methodological aspects in the development and evaluation of heartbeat classifiers were followed in recent works [de Chazal et al., 2004, Jiang and Kong, 2007, Ince et al., 2009, Llamedo and Martínez, 2011a]. The most relevant key-points are: •Use of public and standard databases, as the ones available in Physionet [Goldberger et al., 2000]. •Fulfillment of AAMI recommendations for class labeling and results presentation [AAMI-EC57, 1998–2008]. •Patient-oriented data division into training and testing sets, as described in [de Chazal et al., 2004]. Another aspect suggested in recent works is the analysis of the capability of the classifier to retain its performance in other databases not considered during the development [Llamedo and Martínez, 2011a]. We refer to this property of a classifier as generalization capability, and its analysis provides a broader idea of the performance achieved. Up to the writing of this thesis, only few of the reviewed works used more than one database either for the development [Watrous and Towell, 1995, Kiranyaz et al., 2011] or for a generalization assessment [Chudácek et al., 2009, Krasteva and Jekova, 2007, Syed et al., 2007]. The AAMI EC57 recommendations [AAMI-EC57, 1998–2008] for class labeling and results presentation are at the present time broadly accepted [de Chazal et al., 2004, Inan et al., 2006, Llamedo and Martínez, 2007, Jiang and Kong, 2007, Park et al., 2008, Ince et al., 2009]. As any classification problem, the goal is to learn a function that divides in C regions (or classes) a (hyper) space defined by the features, extracted from the ECG, and then make predictions with this function. In other words, this means assigning a label to an unknown heartbeat as a function of the value of some features. It is not difficult to realize that the lesser the amount of classes (small C), the simpler the partition function. Since cardiologists can group heartbeats into a number of classes that is easily higher than 10, the AAMI EC57 recommendations simplifies the problem into Chapter 2 Materials and Methods In this chapter we describe the materials and methods used in the following chapters. The materials are mainly two, ECG databases and computing resources. The methods used are many more, but can be grouped in signal processing and classification. The first group contains the methods to calculate the feature vectors from the ECG signal, while the second group contains the methods to classify the feature vectors into heartbeat classes. 2.1 ECG Databases All experiments performed in this thesis were carried out in several public databases freely available on Physionet [Goldberger et al., 2000], the well known American Heart Association database [American Heart Association], and a database developed at Biosigna GmbH [Fischer et al., 2008]; their relevant details are summarized in Table 2.3. For all databases the AAMI recommendations for class-labeling were adopted. In the next subsection this recommendation is explained. The AAMI Q class (unclassified and paced heartbeats) was discarded since it is marginally represented in all databases. This limitation occurs to a lesser extent with the fusion (F) AAMI class, but instead of discarding the heartbeats of this class, we adopted an alternative labeling scheme already used in [Llamedo and Martínez, 2011a]. It consists in merging the fusion (of normal and ventricular beats) and ventricular classes, as the same ventricular class (V’ in Table 2.3). We will refer to this modification as AAMI2 labeling. This labeling does not compromise the comparability with other AAMI compliant works, since according to AAMI recommendation [AAMIEC57, 1998–2008], errors involving F and Q classes either should not be accounted for the required performance measurements, or are already accounted by considering the V’ class. The databases used include different types of ECG recordings: some of them were recorded during routine ambulatory practice, but others were selected to include less common ventricular, junctional or supraventricular arrhythmias, or baseline ST segment displacement or other ECG abnormalities. As a result, we use in this work a dataset 29 30 CHAPTER 2. MATERIALS AND METHODS with a broad range of normal and pathological ECG recordings to evaluate the algorithm performance. Further details of each database can be found on Physionet1[Goldberger et al., 2000]. We will refer as the development dataset to the union of the MITBIH-SUP database and the 22 recordings included in the DS1 subset of MITBIH-AR defined in [de Chazal et al., 2004], while the evaluation dataset includes the rest of databases described in Table 2.3. The reason of this division was, first for ensuring a fair comparison with de Chazal et al. [de Chazal et al., 2004] results, and second because MIT Arrhythmia databases have heartbeat annotations thoroughly reviewed by several experts, and are therefore more reliable than the rest of databases. 2.1.1 AAMI class labeling recommendations According to Section 4.2 of [AAMI-EC57, 1998–2008], within annotation files beat labels are defined as follows: Nany beat that does not fall into the S, V, F, or Q categories described below (a normal beat or a bundle branch block beat) Sa supraventricular ectopic beat (SVEB): an atrial or nodal (junctional) premature or escape beat, or an aberrated atrial premature beat Va ventricular ectopic beat (VEB): a ventricular premature beat, an R-on-T ventricular premature beat, or a ventricular escape beat Fa fusion of a ventricular and a normal beat Qa paced beat, a fusion of a paced and a normal beat, or a beat that cannot be classified The description and conversion matrices for the formats used in the databases described below are presented in tables 2.1 and 2.2. 2.1.2 MIT-BIH Arrhythmia Database (MITBIH-AR) The MIT-BIH Arrhythmia Database contains 48 half-hour excerpts of two-channel ambulatory ECG recordings, obtained from 47 subjects studied by the BIH Arrhythmia Laboratory between 1975 and 1979 [Moody and Mark, 2001]. The source of the ECGs included in the MIT-BIH Arrhythmia Database is a set of over 4000 long-term Holter recordings that were obtained by the Beth Israel Hospital Arrhythmia Laboratory between 1975 and 1979. Approximately 60% of these recordings were obtained from inpatients. The first group is intended to serve as a representative sample of the variety of waveforms and artifact that an arrhythmia detector might encounter in routine clinical use. A table of random numbers was used to select tapes, and 1www.physionet.org 2.1. ECG DATABASES 31 Table 2.1: Original annotation format used in the databases. MIT format AHA format Sym Description Sym Description N Normal beat E Ventricular Escape ·Normal beat F Fusion Beat L Left bundle branch block beat N Beat of Non-Ventricular Origin R Right bundle branch block beat P Paced Beat A Atrial premature beat Q Questionable Beat - Indeterminate Origin a Aberrated atrial premature beat R R-on-T Beat J Nodal (junctional) premature beat U Unreadable S Supraventricular premature beat V Premature Ventricular Contraction V Premature ventricular contraction [ Beginning of Ventricular Fibrillation or Flutter F Fusion of ventricular and normal beat ] End of Ventricular Fibrillation or Flutter [ Start of ventricular flutter/fibrillation ! Ventricular flutter wave ] End of ventricular flutter/fibrillation e Atrial escape beat j Nodal (junctional) escape beat E Ventricular escape beat / Paced beat f Fusion of paced and normal beat x Non-conducted P-wave (blocked APB) p Non-conducted P-wave (blocked APB) Q Unclassifiable beat | Isolated QRS-like artifact ? Beat not classified during learning + Rhythm change r R-on-T premature ventricular contraction s ST segment change B Bundle branch block beat (unspecified) 32 CHAPTER 2. MATERIALS AND METHODS Table 2.2: AAMI class conversion matrices for the formats used. MIT format AHA format Original AAMI AAMI2 Original AAMI AAMI2 N N N E V V ·N N F F V L N N N N N R N N Q Q Q A S S R V V a S S V V V J S S S S S V V V F V V HES format [ Q Q Original AAMI AAMI2 ! Q Q Q Q Q ] Q Q N N N e N N N N N j N N S S S E V V Q Q Q / F V V V V f F V S S S x Q Q p Q Q Q F V | Q Q ? Q Q AAMI labels + Q Q N Normal r V V S Supraventricular s Q Q V Ventricular B N N F Fusion n N N Q Unknown 2.1. ECG DATABASES 33 then to select half-hour segments of them. Segments selected in this way were excluded only if neither of the two ECG signals was of adequate quality for analysis by human experts. Records in the second group were chosen to include complex ventricular, junctional, and supraventricular arrhythmias and conduction abnormalities. Several of these records were selected because features of the rhythm, QRS morphology variation, or signal quality may be expected to present significant difficulty to arrhythmia detectors; these records have gained considerable notoriety among database users. The subjects were 25 men aged 32 to 89 years, and 22 women aged 23 to 89 years. (Records 201 and 202 came from the same male subject.) In most records, one signal is a modified limb lead II (MLII), obtained by placing the electrodes on the chest. While the other is usually a modified lead V1 (occasionally V2 or V5, and in one instance V4); as for the first signal, the electrodes are also placed on the chest. Normal QRS complexes are usually prominent in the first signal. The lead axis for the second signal may be nearly orthogonal to the mean cardiac electrical axis, however (i.e., normal beats are usually biphasic and may be nearly isoelectric). Thus normal beats are frequently difficult to discern in the second signal, although ectopic beats will often be more prominent. A notable exception is record 114, for which the signals were reversed. Since this happens occasionally in clinical practice, arrhythmia detectors should be equipped to deal with this situation. In records 102 and 104, it was not possible to use modified lead II because of surgical dressings on the patients; modified lead V5 was used for the first signal in these records. The recordings were digitized at 360 samples per second per channel with 11-bit resolution over a 10 mV range. An initial set of beat labels was produced by a simple slope-sensitive QRS detector, which marked each detected event as a normal beat. Two identical 150-foot chart recordings were printed for each 30-minute record, with these initial beat labels in the margin. For each record, the two charts were given to two cardiologists, who worked on them independently. The cardiologists added additional beat labels where the detector missed beats, deleted false detections as necessary, and changed the labels for all abnormal beats. They also added rhythm labels, signal quality labels, and comments. The annotations were transcribed from the paper chart recordings. Once both sets of cardiologists’ annotations for a given record had been transcribed and verified, they were automatically compared beat-by-beat, and another chart recording was printed. This chart showed the cardiologists’ annotations in the margin, with all discrepancies highlighted. Each discrepancy was reviewed and resolved by consensus. The corrections were transcribed, and the annotations were then analyzed by an auditing program, which checked them for consistency and which located the ten longest and shortest R-R intervals in each record (to identify possible missing or falsely detected beats). The annotations provided with the database were used for training and testing purposes, following the recommendations and class-labeling of AAMI (See tables 2.1 and 2.2). 34 CHAPTER 2. MATERIALS AND METHODS We adopted the training (DS1) and test (DS2) set division scheme used in [de Chazal et al., 2004] for comparative purposes of the results. The division scheme is summarized in Table 2.3. 2.1.3 MIT-BIH Supraventricular Arrhythmia Database (MITBIHSUP) The database consists of 78 two-lead recordings of approximately 30 minutes and sampled at 128 Hz. The recordings were chosen to supplement the examples of supraventricular arrhythmias in the MIT-BIH Arrhythmia Database. The annotations of the recordings were performed automatically first, by the Marquette Electronics 8000 Holter scanner and later reviewed and corrected by a medical student [Greenwald, 1990]. These reference annotations are considered a “silver standard” because the granularity of the beat labels only discriminates among normal, ventricular, supraventricular and fusion beats. The original labeling was also adapted to the AAMI recommendations and to the AAMI2 modification. This database will be considered for validation and model selection purposes. The class distribution is shown in Table 2.3. 2.1.4 St. Petersburg Institute of Cardiological Technics (INCART) 12-lead Arrhythmia Database This database consists of 75 annotated recordings extracted from 32 Holter records. Each record is 30 minutes long and contains 12 standard leads, each sampled at 257 Hz. The reference annotation files contain over 175000 beat annotations in all. The original records were collected from patients undergoing tests for coronary artery disease (17 men and 15 women, aged 18-80; mean age: 58). None of the patients had pacemakers; most had ventricular ectopic beats. In selecting records to be included in the database, preference was given to subjects with ECG’s consistent with ischemia, coronary artery disease, conduction abnormalities, and arrhythmias. These diagnoses were confirmed by enzyme assays, coronary angiography, electrophysiological study, and pressure monitoring where necessary. For each record it was included the patient’s age, sex, diagnoses, and a summary of features of the ECG. The annotations were produced by an automatic algorithm and then corrected manually, following the standard PhysioBank beat annotation definitions. The algorithm generally places beat annotations in the middle of the QRS complex (as determined from all 12 leads); the locations have not been manually corrected, however, there may be occasional misaligned annotations as a result. This database will be considered only for testing purposes. More details about the database are shown in Table 2.3. 2.1. ECG DATABASES 35 2.1.5 European ST-T Database (ESTTDB) This database consists of 90 annotated excerpts of ambulatory ECG recordings from 79 subjects. Myocardial ischemia was diagnosed or suspected for each subject, additional selection criteria were established in order to obtain a representative selection of ECG abnormalities in the database, including baseline ST segment displacement resulting from conditions such as hypertension, ventricular dyskinesia, and effects of medication. Each record is two hours in duration and contains two signals, each sampled at 250 Hz with 12-bit resolution over a nominal 20 millivolt input range. The European ST-T Database is intended to be used for evaluation of algorithms for analysis of ST and T-wave changes. This database consists of 90 annotated excerpts of ambulatory ECG recordings from 79 subjects. The subjects were 70 men aged 30 to 84, and 8 women aged 55 to 71. Myocardial ischemia was diagnosed or suspected for each subject; additional selection criteria were established in order to obtain a representative selection of ECG abnormalities in the database, including baseline ST segment displacement resulting from conditions such as hypertension, ventricular dyskinesia, and effects of medication. The database includes 367 episodes of ST segment change, and 401 episodes of T-wave change, with durations ranging from 30 seconds to several minutes, and peak displacements ranging from 100 microvolts to more than one millivolt. In addition, 11 episodes of axis shift resulting in apparent ST change, and 10 episodes of axis shift resulting in apparent T-wave change, have been marked. Compact clinical reports document each record. These reports summarize pathology, medications, electrolyte imbalance, and technical information about each recording. Each record is two hours in duration and contains two signals, each sampled at 250 samples per second with 12-bit resolution over a nominal 20 millivolt input range. The sample values were rescaled after digitization with reference to calibration signals in the original analog recordings, in order to obtain a uniform scale of 200 ADC units per millivolt for all signals. The header files include information about the leads used, the patient’s age, sex, and medications, the clinical findings, and the recording equipment. A complete description of the database can be found in [Taddei et al., 1992]. 2.1.6 The MIT-BIH ST Change Database (MITBIH-ST) This database includes 28 ECG recordings of varying lengths, most of which were recorded during exercise stress tests and which exhibit transient ST depression. We selected the two-lead recordings, resulting in 18 useful recordings. The recordings were sampled at 360 Hz and 12-bit resolution. 36 CHAPTER 2. MATERIALS AND METHODS 2.1.7 The Long-Term ST Database (LTSTDB) The Long-Term ST Database contains 86 lengthy ECG recordings of 80 human subjects, chosen to exhibit a variety of events of ST segment changes, including ischemic ST episodes, axis-related non-ischemic ST episodes, episodes of slow ST level drift, and episodes containing mixtures of these phenomena. The database was created to support development and evaluation of algorithms capable of accurate differentiation of ischemic and non-ischemic ST events, as well as basic research into mechanisms and dynamics of myocardial ischemia. Detailed clinical notes and ST deviation trend plots are provided for all 86 records. The entire Long-Term ST Database is also available from its original home page at the Laboratory for Biomedical Computer Systems and Imaging at the University of Ljubljana, Slovenia The individual recordings of the Long-Term ST Database are between 21 and 24 hours in duration, and contain two or three ECG signals. Each ECG signal has been digitized at 250 samples per second with 12-bit resolution over a range of ±10 millivolts. Each record includes a set of meticulously verified ST episode and signal quality annotations, together with additional beat-by-beat QRS annotations and ST level measurements. Several sources contributed recordings to the Long-Term ST Database: •Eleven of the recordings included in the Long-Term ST Database are from the initial Long-Term ST Database developed under a joint U.S.-Slovenian research project between 1995 and 1998. •Ten additional recordings of the Long-Term ST Database are from the collection originally gathered by the Pisa group for the European ST-T Database, which contains two-hour excerpts of some of these same recordings. The original analog recordings were redigitized for the Long-Term ST Database; since the signals have been rescaled as a result, direct comparison of the annotations in the European ST-T Database records with those for the corresponding portions of the Long-Term ST Database records is not possible. The inclusion of these recordings in the LongTerm ST Database allows study of the dynamics of ischemic ST changes over a much longer period in these previously well-studied subjects. Among the samples available here, record s20021 includes the two-hour segment that was previously digitized to produce record e0113 of the European ST-T Database. •Another 18 of the LTSTDB recordings, those containing recordings with three ECG signals, were contributed to the project by Zymed, Inc. The annotation of the Long-Term ST Database was performed using SEMIA, a program written by the group in Ljubljana for this purpose. Each recording was reviewed independently by expert annotators using SEMIA at each of the three sites (Ljubljana, Pisa, and Cambridge). Participants met several times annually to obtain the consensus reference annotations. For further details, see [Jager et al., 2003]. 2.2. SUPERCOMPUTING RESOURCES 37 2.1.8 American Heart Association (AHA) ECG Database The American Heart Association (AHA) sponsored the development of the AHA Database for Evaluation of Ventricular Arrhythmia Detectors during the late 1970s and early 1980s at Washington University (St. Louis). The first portions of the AHA Database were released in 1982, and it was completed in 1985. Until recently, the only available portion of the AHA database consisted of 80 two-channel excerpts of analog ambulatory ECG recordings, digitized at 250 Hz per channel with 12-bit resolution over a 10 mV range [American Heart Association]. These recordings, designated as the development set, are divided into eight classes of ten recordings each, according to the highest level of ventricular ectopy present: •no ventricular ectopy (records 1001 through 1010) •isolated unifocal PVCs (records 2001 through 2010) •isolated multifocal PVCs (records 3001 through 3010) •ventricular biand trigeminy (records 4001 through 4010) •R-on-T PVCs (records 5001 through 5010) •ventricular couplets (records 6001 through 6010) •ventricular tachycardia (records 7001 through 7010) •ventricular flutter/fibrillation (records 8001 through 8010) The final thirty minutes of each recording are annotated beat-by-beat, although supraventricular ectopic beats are not distinguished from normal sinus beats. At the time the AHA Database was created, a second set of 75 recordings (designated as the test set) was constructed according to the same criteria as the development set (only 5 recordings in the R-on-T PVC class were included in the test set). The test set was intended for evaluations without any possibility that the detectors might have been tuned (optimized) for the test data; for this reason, the test set was unavailable until recently. 2.2 Supercomputing Resources The development and evaluation of the algorithms presented in this thesis involves a lot of computation. This kind of tasks are unfeasible in ordinary computers, since the time required for an adequate evaluation can easily reach several days. All the results presented in this thesis were calculated using the resources of the Instituto de Investigación en Ingeniería de Aragón (I3A). HERMES is the I3A’s high throughput computing cluster, and this is a brief overview of its main features: 38 CHAPTER 2. MATERIALS AND METHODS Table 2.3: Databases and datasets used in this thesis with its class representation. Database Purpose Length Leads Fs (Hz) N S V’ Q #Rec V F Development DS1of MITBIH Arrhythmia (MITBIH-AR)1,2A 30’ 2 360 45673 929 3755 412 6 22 MITBIH Supr. Arrhythmia (MITBIH-SUP)1A 30’ 2 128 162271 12195 9940 23 79 78 Evaluation DS2of MITBIH Arrhythmia (MITBIH-AR)1,2A 30’ 2 360 44053 1833 3202 388 7 22 American Heart Association (AHA)3A 30’ 2 250 319125 0 32745 1266 0 155 European ST-T (ESTTDB)1ST 2h 2 250 784568 1095 4467 354 0 90 MITBIH ST Change (MITBIH-ST)1ST <1h 2 360 46215 798 319 0 0 18 MITBIH long-term database (MITBIH-LT)1LT 14-22h 2 128 594672 1499 63584 2785 0 7 Long-Term ST (LTSTDB)1LT/ST 21-24h 2/3 250 6422003 30905 37894 0 86 Biosigna4A 1h 12 500 286246 1326 2541 0 0 56 St. Petersburg Inst. of Card. Tech. (INCART)1A 30’ 12 257 153651 1959 20005 219 6 75 Total 8858840 52556 178502 5449 100 609 A: Arrhythmia/Heartbeat classification; ST: Ischemia detection; LT: Long-term detection/classification Heart beats classes are N: normal, S: supraventricular, V: ventricular, F: fusion, Q: unknown, and V’:AAMI2 Ventricular. 1[Goldberger et al., 2000]; 2[de Chazal et al., 2004]; 3[American Heart Association]; 4[Fischer et al., 2008] Dataset MITBIH-AR recording names DS1101, 106, 108, 109, 112, 114, 115, 116, 118, 119, 122, 124, 201, 203, 205, 207, 208, 209, 215, 220, 223, 230 DS2100, 103, 105, 111, 113, 117, 121, 123, 200, 202, 210, 212, 213, 214, 219, 221, 222, 228, 231, 232, 233, 234 2.3. SIGNAL PROCESSING 45 å(t) 1.5 ×(t) Time (s) -1.5 -0.75 0 0.75 3 -3 0 Wavelet prototype and smoothing functions Figure 2.4: Wavelet prototype used in this thesis. A quadratic spline that matches the derivative of the convolution of four rectangular pulses. and the low-pass and high-pass FIR filters have transfer functions [Li et al., 1995] Hejω=ejω/2cos ω 23 , Gejω= 4 j ejω/2sin ω 2, (2.11) with associated impulse responses h[k]=1/8·{δ[k+ 2]+3δ[k+ 1]+3δ[k]+δ[k−1]}, g[k]=2·{δ[k+ 1]−δ[k]}.(2.12) This prototype wavelet has been applied to detection and delineation of ECG signals with good results [Li et al., 1995, Martínez et al., 2004]. As the analysis filters in equation (2.12) have linear phase [Li et al., 1995], the outputs of the filters can be realigned in order to present the same delay with respect to the original signal s(t). The equivalent frequency responses Qm(ejω)for this prototype wavelet using the algorithme à trous can be calculated from equations (2.8) and (2.11). The frequency responses are plotted in Figure 2.5 for the first five scales a= 2m|m=1,2,...,5, considering a sampling frequency Fs of 250 Hz. In the rest of the thesis, the wavelet scale a= 2mwill be referred as the m-th scale. For a value of Fsdifferent from 250 Hz the bands in Figure 2.5 would appear scaled in frequency. It is important to keep the scale fitting to the ECG features on the algorithms independent of Fs: for that purpose a new set of filters, having equivalent frequency responses as close as possible to the ones of Figure 2.5 are constructed for each Fs. The new filters are obtained by resampling adequately the equivalent filter impulse responses 46 CHAPTER 2. MATERIALS AND METHODS Frequency (Hz) Magnitude (dB) Magnitude Response (dB) 0 20 40 60 80 100 120 -90 -70 -50 -30 -10 10 1 2 3 4 5 Scales Figure 2.5: Transfer functions of the filter-bank used to calculate the DWT up to the 5th scale. at 250 Hz. Such a procedure is required to construct a system able to handle equivalently ECG signals with different sampling frequencies. As a final precision of Figure 2.5, the DWT filter-bank response is 250 Hz, but the Fsfor this implementation is 360 Hz, which is the sampling rate of the recordings included in the MITBIH-AR. Following the conclusions of [Martínez et al., 2004], the resulting DWT framework allows an analysis robust to the typical interferences present in routine ECG recordings, so the features derived from the DWT are expected to inherit this desirable property. 2.4 Heartbeat classification This section includes much of the work carried out in this thesis. In the first two subsections are described the features and classifiers used, which ultimately are the methodologies that will be used for the classification of heartbeats. The other subsections are important during the development of the classification algorithm. 2.4.1 Classification Features In this subsection all the features that were used in this thesis are defined. However, only a small subset of them are selected to be implemented in the final classification model, as will be described below in section 2.4.6. Many of these features were initially designed to 2.4. HEARTBEAT CLASSIFICATION 47 mean(RR ) 20min mean(RR ) 1min current beat RR(n) RR(n-1) RR(n+1) Figure 2.6: Transfer function of the low-pass filter used for ECG preprocessing. be used in two-lead recordings, as those included in the MITBIH-AR and MITBIH-SUP databases. Later in chapter 4 we present a strategy to adapt these features to recordings with an arbitrary amount of leads. Following the conclusions of previous works [Hu et al., 1997, de Chazal et al., 2004], we included in our model both rhythm and morphological features. The rhythm features are designed to model the manner in which the heartbeats succeed each other. They are responsible of modeling the sinus node activity, an increased automaticity in any part of the heart or a complete A-V block. As rhythm features we used features from the RR interval (RR) sequence, being RR[i] = Q[i]−Q[i−1] (2.13) the current RR interval the difference between the current (Q[i]) and previous (Q[i−1]) QRS fiducial points (FP). Then we used RR[i−1],RR[i]and RR[i+ 1] to describe the local time evolution of the heart rhythm. In order to assess the local variation of the heart rhythm, the feature RRV[i] = 1 X j=−1|dRR[i−j]|,(2.14) being dRR[i] = RR[i]−RR[i−1], characterizes the amount of RR variation in the surrounding heartbeats. The prematurity of a heartbeat, defined as PRR[i] = RR[i] Pi+1 k=i−1RR[k],(2.15) measures how anticipated is a heartbeat respect to the previous and next RR interval. We also included estimates of the local and global rhythm by the mean RR interval in the last 1, 5, 10 and 20 minutes (RRPbeing P∈ {1,5,10,20},the interval in minutes of aggregation). Examples of the rhythm features used are shown in Figure 2.6. The morphological features are designed to model the manner in which heartbeats are conducted through the heart. The features that we used can be grouped in three cate- 48 CHAPTER 2. MATERIALS AND METHODS QRSw normal QRSw ventricular V1 V1 Normal beat Ventricular beat QRSoff QRSon QRSoff Figure 2.7: QRS width measured for a normal and a ventricular heartbeat. gories depending on whether they were calculated in the ECG signal, the two-dimensional vectocardiogram (VCG) loop formed by both available leads or in the DWT of the ECG signal. 1) The QRS width is obtained from the delineation of the ECG, as a simple difference QRSW=QRSoff −QRSon of the offset and offset of the QRS complex. The QRS location is assumed to be known, for this purpose we used the annotations included in the databases. Following the QRS complex detection positions, the delineation of each heartbeat was performed with the delineator described in [Martínez et al., 2004]. It is well known that not properly conducted heartbeats, as the ventricular class, have widened QRS complexes as is shown in Figure 2.7. From the clinical point of view, a QRSW>120 ms is considered to be widened. 2) From the 2-D VCG loop constructed with the two available leads (see Figure 2.8) we calculated two features: the maximal vector of the QRS loop (V CGM) and the angle of this vector (V CGφ). As any problem in the conduction of the impulse through the heart should change the loop morphology, and therefore this two features, as shown in Figure 2.8 for a normal and a ventricular heartbeat. 3) Regarding the features calculated from the DWT of the ECG, four types can be defined: 3.a) The first type includes 7 features (per lead) that were calculated from peak amplitudes and positions from the fourth scale of the DWT (Ws 4(k)), since this scale (between 12.25–22.5 Hz) has good projection of the ECG information. These 7 features are the 2 greatest absolute values of the QRS complex, the 2 greatest absolute values of the T 2.4. HEARTBEAT CLASSIFICATION 49 MIT-BIH Rec. 202 Normal Ventricular MLII V1 V1 MLII VCG ¿ VCG ¿ VCGM VCGM Normal Ventricular Figure 2.8: Illustration of the features calculated from the VCG loop computed with the two available leads, for a normal (continuous line) and ventricular (dotted line) beats. The maximum value of the loop and the angle at this point are shown. 50 CHAPTER 2. MATERIALS AND METHODS Normal Ventricular Scale 3 Scale 4 Scale 5 Scale 6 MLII MIT-BIH Rec. 202 S =3.9 QRS S =4.8 QRS d1 d3 d2 d1 d3 d2 Figure 2.9: Illustration of the features calculated from the wavelet transform for the same normal and ventricular beat in Fig. 2.8. The two most important peaks from the QRS complex and T wave are indicated with an asterisk, and the relative distances (di) to the most important peak in the fourth scale. Also the scale where the QRS complex is centered (SL QRS) is shown for both types of heartbeats used for its calculation (only for one lead). wave, and their 3 relative positions (to the position of the greatest peak in the heartbeat, see Fig. 2.9). 3.b) The second type is also calculated from the fourth scale of the DWT, but from the autocorrelation sequence r(4) x(k) = N X i=0 Wx 4(i)·Wx 4(i−k)(2.16) of the xlead in a window of Nsamples. This sequence of 2N−1samples is symmetric, as can be seen in Figure 2.10. The cross-correlation is very similar r(4) xy (k) = N X i=0 Wx 4(i)·Wy 4(i−k)(2.17) but involves the two available leads, xand y. In contrast, the cross-correlation sequence 2.4. HEARTBEAT CLASSIFICATION 51 is not symmetric as shown in the bottom of Figure 2.10. The autocorrelation sequence for both leads (r(4) x(k)and r(4) y(k)) and the inter-lead cross-correlation sequence (r(4) xy (k)) were calculated within a time window which starts 130 ms before the FP and ends 200 ms after. One remarkable aspect is that features calculated from the correlation signals will essentially be synchronized in time, even if the FP is not accurately determined. We calculated for the 3 signals the location and value of the absolute maximum, and for r(4) xand r(4) ythe location of the first zero-crossing, as shown in Fig. 2.10. 3.c) The third type of feature is the wavelet scale where the QRS complex is centered for each lead. It is known that fast evolving signals (like a normal beat) tend to have their energy located in lower wavelet scales (higher frequency content). The QRS center scale for each lead (SLead QRS )is calculated as the weighted sum SL QRS =P6 s=1 AL s.s P6 s=1 AL s (2.18) where AL sis the mean absolute amplitude of the QRS peaks at scale sof the DWT, and lead L AL s=1 D D X d=1 WL ss(ld), s = 1,2, . . . 6(2.19) being Dthe number of detected peaks (1 or 2) and ldthe positions of the peaks. In Figure 2.9 the peaks are marked with asterisks, and an arrow points the SQRS feature for a normal and ventricular heartbeats. 3.d) The last type of features computed in the DWT were energies measured at several scales and locations of the heartbeat. Specifically during the occurrence of the P wave, and the QRS complex. For that purpose we defined a window of length 400 ms ending at the QRSon. This window was divided in 4 sections of 100 ms, and the energy of the WRMS s(k) signal was computed in each section in scales 3 to 5. This amounts 12 features per beat (4 sections ×3 scales). The same concept was used to study the QRS complex morphology with a window starting at 100 ms before the FP, or the QRS complex location, and ending 120 ms after, for scales 2 to 5, as can be seen in Figure 2.11. In total 16 features per beat. 2.4.2 Discriminant Functions In this thesis we focused our efforts in developing a classification model which includes features that provide most of the classification power. To encourage this property we used in our experiments very simple classifiers as linear and quadratic discriminant classifiers (or functions). All these classifiers assume that the data is distributed following a Gaussian distribution, which is rarely true. In spite of this evident limitation, these classifiers obtain an acceptable trade-off between the sub-optimal performance achieved, and the simplicity to understand them during operation, and to interpret the effect of their features in the 52 CHAPTER 2. MATERIALS AND METHODS Normal Ventricular 130 ms 200 ms 130 ms 200 ms MLII Scale 4 V1 Scale 4 k MLII Scale 4 V1 Scale 4 MIT-BIH Rec. 202 k M x k Z x k M y k Z y k M xy k M xy k k M y k Z y k M x k Z x FP FP r (k) x (4) r ( ) xk M x (4) r (k) x (4) r ( ) xk M x (4) r (k) y (4) r ( ) xk M y (4) r (k) xy (4) r (k) xy (4) r ( ) k M xy xy (4) r ( ) k M xy xy (4) r ( ) xk M y (4) r (k) y (4) Figure 2.10: Illustration of the features calculated from the wavelet correlation signals for the same normal and ventricular beats. The autocorrelation signal of the QRS complex at scale 4 is shown for both leads (rxand ry) as well as the cross-correlation signal (rxy) at the bottom. The zero-crossings and peaks of interest are indicated with an asterisk. 2.4. HEARTBEAT CLASSIFICATION 53 x(k) QRS 220 ms 400 ms MLII V1 W (k) 22 x W (k) 23 x W (k) 24 x W (k) 25 x y(k) W (k) 22 y W (k) 23 y W (k) 24 y W (k) 25 y w W (k) 22 RMS W (k) 23 RMS W (k) 24 RMS W (k) 25 RMS Figure 2.11: Excerpt of record 201 of MIT-BIH database. Normal (N) and ventricular (V) AAMI class heart beats. In the top Figures both ECG leads are shown with their corresponding wavelet decomposition (scales 2-5). The lower panel depicts the RMS composition of both leads wavelet transform (WRMS s(k)). Some features measured in the WRMS s(k)signal are also shown. 54 CHAPTER 2. MATERIALS AND METHODS discriminant functions. This later property is not feasible in other classifiers as support vector machines and neural networks for example. Under the assumption of normally distributed data, the maximum a posteriori classification criterion (MAP) leads to quadratic discriminant functions, broadly used for classification purposes [van der Heijden et al., 2005]. In the general case, the quadratic discriminant function of the i-th class and feature vector x, can be written as gi(x) = −1 2xTΣ−1 ix+µT iΣ−1 ix−1 2µT iΣ−1 iµi−1 2log(|Σi|) + log(P(ωi)),(2.20) being µi,Σiand P(ωi)the mean vector, covariance matrix and prior probability of the i-th class. The classification rule assigns xto the class iwhich results in the maximum posterior probability P(i|x) = egi(x) PC j=1 egj(x).(2.21) As can be seen in (2.21), the denominator is equal for all classes and finding the maximum posterior probability is equivalent to find the maximum gi(x). The values of µiand Σi were computed from the training data with the sample mean and covariance matrix expressions µi=1 Mi Mi X m=1 xm(2.22) Σi=1 Mi−1 Mi X m=1 (xm−µi)·(xm−µi)T(2.23) being Mithe number of examples (xm) of the i-th class. The values for the prior probabilities P(ωi)were considered the same for all classes.In the case that the covariance matrix Σis considered to be the same for all classes (Σi=Σj=Σ,∀i6=j), the quadratic discriminant classifier (QDC) becomes linear in xleading to the linear discriminant classifier (LDC) g, i(x) = µT iΣ−1x−1 2µT iΣ−1µi+ log(P(ωi)),(2.24) where Σcan be estimated as the weighted sample covariance Σ=1 PC i=1 wi C X i=1 wiPMi m=1(xm−µi)·(xm−µi)T Mi ,(2.25) being Cthe total amount of classes and withe class-weighting coefficients. This classweighting possibility is of much interest due to the heavy imbalance of the class-sizes inherent to this application, where the normal class is in general one order of magnitude (at least) more represented than other classes. We refer as LDC to the linear classifier where wi=wj,∀i6=j, any other weight scheme will be referred as compensated linear 2.4. HEARTBEAT CLASSIFICATION 61 Feature 1 Feature 2 Dataset A semi-wrapped 0 Þ Þ 2 Þ 2 - Þ - 0ÞÞ 2 -Þ 2 Þ - Feature 1 Feature 2 Dataset A semi-wrapped 0 Þ Þ 2 Þ 2 - Þ - 0ÞÞ 2 -Þ 2 Þ - Feature 1 Feature 2 Dataset A semi-wrapped 0 Þ Þ 2 Þ 2 - Þ - 0ÞÞ 2 -Þ 2 Þ - Feature 1 Feature 2 Dataset A semi-wrapped Figure 2.16: In the top panels, semi-wrapped dataset using the approximated wrapped Gaussian distribution and the linear Gaussian distribution respectively. In the bottom panels the decision region for both classes and both distributions is showed. Region filled with red color is for Class 1, corresponding to the blue crosses. Note the huge difference on classification performance. 62 CHAPTER 2. MATERIALS AND METHODS With outliers -0.5 0 0.5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 ln(RR(i)) ln(RR(i+1)) -0.5 0 0.5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 ln(RR(i)) ln(RR(i+1)) Without outliers -0.5 0 0.5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 ln(RR(i)) ln(RR(i+1)) -0.5 0 0.5 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 ln(RR(i)) ln(RR(i+1)) Figure 2.17: Example of the outliers removal when estimating the parameters for both types of classifiers used. In both examples the features used are RR[i−1] and RR[i]. Note the different shape of the classification regions as a consequence of the the outliers removal in the lower panels. outliers. With this assumption of slightly contaminated data, we set an operating point for the trade-off between discarding useful data and allowing the presence of outliers in the parameter estimation process. In Figure 2.17 the use of the algorithm for outliers removal is evidenced in the classification regions, as can be further seen for the QDC classifier. 2.4.5 Performance evaluation We use three approaches to evaluate the performance of a classification experiment. The first consists in estimating the parameters of the classifier in a training dataset, then with these parameters, evaluate their predictive performance in a disjoint test dataset. The data division should be performed by subject, which means that feature-vectors (or 2.4. HEARTBEAT CLASSIFICATION 63 subject 1 subject 2 subject N Whole database heartbeat 1 heartbeat 2 heartbeat M subject 1 subject 2 subject P subject 1 subject 2 subject Q Train Test Params. Features Model Trained Model Results Posteriors subject 2 subject 3 subject N Train 1 a) Normal evaluation b) Crossvalidation c) Optimistically biased subject 1 Test 1 FM TM R1 subject 1 subject 3 subject N Train 2 subject 2 Test 2 FM TM R2 subject 1 subject 2 subject N-1 Train N subject N Test N FM TM RN subject 1 subject 2 subject N Train - Test Features Model Params. Trained Model Results Posteriors Results P+Q = N subject N FM: Features Model TM: Trained Model Rx: Results X Figure 2.18: Scheme showing the three data division methods used. heartbeats) of the same subject can not be included in both the training and testing sets. In [de Chazal et al., 2004] it was shown that the heartbeat oriented data division leads to an optimistic bias in the performance estimation calculated in the test dataset. The second approach used is known as cross-validation, and is adopted when there are not many examples to build the train and test datasets. It consists in dividing a dataset into kdisjoint folds (we used 10 folds), and use them as test sets, obtaining therefore k performance measures, one for each fold. The division is again performed by patients, to avoid the presence of heartbeats of the same patient in both training and test datasets. Note that each cross-validation fold implies training in k−1 /kof the database patients, and testing in the remaining 1 /k. The resulting performance is the mean of the kperformances. In the last approach, although training and testing in the same dataset leads to an optimistically biased performance estimation, this measure could be useful to have an idea of how good can perform a given classifier in the most favorable situation. This optimistically biased performance serves as an upper bound, and represents the performance of the classifier if the probability distributions of examples in both training and test datasets were identical. In Figure 2.18 a scheme shows the data division methods used. Now we focus in how the performance is evaluated for any of the approaches described above. According to Section 4.3 of [AAMI-EC57, 1998–2008], the end product of a beatby-beat comparison between the algorithm and reference labels is a matrix in which each element is a count of the number of beat label pairs of the appropriate type: 64 CHAPTER 2. MATERIALS AND METHODS Algorithm label N s v f q Reference label N Nn Ns Nv Nf Nq S Sn Ss Sv Sf Sq V Vn Vs Vv Vf Vq F Fn Fs Fv Ff Fq Q Qn Qs Qv Qf Qq Moreover, from Subsection 2.4.2, we showed that the outcome of the LDC, LDC-C or QDC classifiers (eq. 2.20 and 2.24), were the posterior probabilities of the classes for each of the heartbeats presented to a trained classifier (see also Figure 2.18). As said before, for each heartbeat the classification rule selects the class with greater posterior. Therefore, given the ground truth or true classes for each example, the number of possible outcomes for a multiclass classification problem are C×C. The amount of events for each of the possible outcomes are represented in a square matrix of dimension C, known as confusion matrix. This is the same matrix with AAMI classes, but for any multiclass problem. Estimated classes True classes 1. . . i . . . C 1 . . . i . . . C             nT 11 . . . nF 1i. . . nF 1C . . ..... . . nF i1. . . nT ii . . . nF iC . . .. . .... nF C1. . . nF Ci . . . nT CC             N1 . . . Ni . . . NC P1. . . Pi. . . PCNT. For the i-th class nT ii is the number of correctly classified examples and nF ij is the number of examples of class iclassified as class j;Niis the total number of examples for class i,Piis the number of examples classified as class iand NTis the total number of examples in the dataset. Ni=nT ii +X m6=i nF im Pi=nT ii +X m6=i nF mi NT= C X i=1 Ni= C X i=1 Pi The performance is calculated from this matrix, for each class, in terms of the class sensitivity (Si) and class positive predictive value (P+ i) Si=nT ii Ni (2.36) 2.4. HEARTBEAT CLASSIFICATION 65 P+ i=nT ii Pi ,(2.37) or globally, with the global accuracy (A), sensitivity (S) and positive predictive value (P+) A=1 NT C X i=1 nT ii = C X i=1 Ni NT Si(2.38) S=1 C C X i=1 Si(2.39) P+=1 C C X i=1 P+ i(2.40) as suggested in [AAMI-EC57, 1998–2008]. From these equations it is clear that any imbalance in the class representation affects P+,P+ iand Acalculation, but not Sand Si. Although the AAMI recommendation does not suggest any measure to deal with the strong class size imbalance (see Table 2.3), we considered weighting the classes previous to the calculation of P+ iand Ain order not to neglect the performance of the less represented classes. The balancing approach used in this work consists in multiplying each row of the confusion matrix by a constant such that the sum of each row Niis equal for all classes, or Ni=Nj,∀i6=j. This is equivalent to repeat examples of the less represented classes, in order to balance the class presence. 2.4.6 Model Selection and Dimensionality Reduction The amount of features included in a classification model should be limited, mainly, because the training set is rarely able to fully describe the class distributions. Then, one of the greatest challenges during the training of a classifier is to avoid overfitting. Overfitting happens when the trained classifier is adapted in excess to the training examples. It is somehow the opposite of generalization in classification. In section 2.4.2 we presented our first strategy to limit the complexity of our classifier, by using linear or quadratic discriminant functions. In this section we present our approach to control the complexity of the features model. It is well known that low dimensional models generalize better to examples not presented during the training phase, resulting in a more robust and realistic classifier. There are two approaches typically used to achieve this, feature transformation (FT) or feature selection (FS). As can be anticipated, both have pros and cons. FT can be defined as a mapping y=W(z)(2.41) of the input space of dimension Ninto another of dimension D, given that DN. One of the advantages of FT is that take into account all the features in the input space, which can also be thought as a disadvantage in case that some features are useless 66 CHAPTER 2. MATERIALS AND METHODS for classification. There are many criteria to the design of the mapping function W, which imposes assumptions over the data [van der Heijden et al., 2005]. The same design requirements in term of training data, overfitting and generalization should be taken into account for the correct development of W. This means that the design of Wis as important as the training process itself. On the other hand, the wrapper approach consists in reducing the feature space by selecting the Dmost important features given a criterion function J(·). We can comment similar pros and cons of FT for FS. However, a property of FS is that the feature model selected is tailored for a specific classifier. If the classifier changes, there is no guarantee that the same feature model maximizes the criterion J(·). In general, the FS problem consists in evaluating J(·)across several feature models. The exhaustive and optimal evaluation across all model combinations requires N!evaluations, which is unfeasible even for moderate N. However, many algorithms exist to overcome this, obtaining a suboptimal model. In this thesis a sequential floating feature selection algorithm (SFFS) was used [Pudil et al., 1994] to obtain the smallest and best performing model. The SFFS algorithm can be briefly explained as the combination of two simpler steps, a sequential forward selection (SFS) algorithm followed by a sequential backwards selection (SBS) algorithm. Following Figure 2.19, the SFFS iterates for all model sizes, starting from a model of one feature, and registering all the best performances found for each model size. Each iteration starts with an SFS step, and from a model size greater than two features after each SFS step, an SBS step is repeated until the performance of the model found is not greater than the registered for this smaller model size. This way the algorithms goes forward and backwards (like floating) searching at each step for the path of maximum performance. The algorithm ends when the specified greater model size is reached. The result of the algorithm is the model found with maximum performance. The interested reader is referred to [Pudil et al., 1994] for a detailed description and to [Duin et al., 2008] for an implementation of the SFFS algorithm. Two examples of the results obtained from this algorithm are shown in Figure 3.2. The performance metrics used by the FS algorithm were a weighted class Se and PP, calculated as JS=PC i=1 πi.Si PC i=1 πi (2.42) JP+=PC i=1 πi.P+ i PC i=1 πi (2.43) with C classes and being Siand P+ ithe class sensitivities and positive predictive defined in the previous subsection. The class weights πiallows the possibility of directing the search to specific class performances, specially in the case of the P+which is directly influenced by the class size unbalance. Then, as a result of using the JScriterion, the 2.4. HEARTBEAT CLASSIFICATION 67 Current best subset of size k Apply a step of the SFS algorithm Start k 0 k k +1 k = d ? Yes No Stop Conditionally excude one feature applying a step of the SBS algorithm Return the conditionally excluded feature back to the current solution k k - 1 Remove the conditionally excluded feature from the current solution Yes No Is this subset better than the current best subset of size k - 1 Figure 2.19: Flow diagram of the sequential floating feature selection (SFFS) algorithm used for the feature selection among dfeatures. SFFS is likely to find sensitive feature models, in the other hand, the JP+criterion will yield specific feature models. 68 CHAPTER 2. MATERIALS AND METHODS Chapter 3 Automatic ECG Heartbeat Classification 3.1 Introduction In this chapter we describe the development of a simple heartbeat classifier, based on the methodologies described in the previous chapter, that will be enhanced in the following chapters. In section 1.3 we described other approaches published by other groups, but in our opinion the work of de Chazal et al. [de Chazal et al., 2004] is the most representative of the state of the art, and will be used as a reference. The main novelties with respect to [de Chazal et al., 2004] are: 1) the morphological features used are based on the wavelet transform (Section 2.4.1), 2) the use of a feature selection algorithm to select those features with generalization capability (Section 2.4.6) and finally 3) the evaluation of the generalization of the classifier outside the MITBIH Arrhythmia database (MITBIHAR). The results of this chapter were published in [Llamedo and Martínez, 2011a]. In summary, the objective pursued in this chapter is to develop and evaluate a heartbeat classification algorithm according to the following conditions: •Perform automatic ECG classification. •Follow AAMI [AAMI-EC57, 1998–2008] recommendations for class labeling and results presentation. •Use a simple classifier (as linear or quadratic discriminant functions) to ensure that the classification performance is due to the features selected. •The features used should have a physiological interpretation, being simple to compute and robust to the typical kind of noise present in the ECG. •Preference for features that can be used for ECG delineation, as those from the discrete wavelet transform (DWT), since the heartbeat classifier is intended to be used after a DWT based ECG delineator [Martínez et al., 2004]. 69 70 CHAPTER 3. AUTOMATIC ECG HEARTBEAT CLASSIFICATION Table 3.1: Class distribution of the databases used and division of the MITBIH-AR database into training (DS1) and testing (DS2) sets. Recordings with paced beats were excluded. MITBIH-AR Dataset Purpose N S V F Q #Rec DS1train 45784 940 3783 413 8 22 DS2test 44188 1835 3218 388 7 22 Totals 89972 2775 7001 801 15 44 Heartbeat classes are N: normal, S: supraventricular, V: ventricular and F: fusion Dataset MITBIH-AR recordings DS1101, 106, 108, 109, 112, 114, 115, 116, 118, 119, 122, 124, 201, 203, 205, 207, 208, 209, 215, 220, 223, 230 DS2100, 103, 105, 111, 113, 117, 121, 123, 200, 202, 210, 212, 213, 214, 219, 221, 222, 228, 231, 232, 233, 234 Other Databases Database Purpose N S V F Q #Rec MITBIH-SUP validation 161902 12083 9897 193 78 78 INCART test 153517 1958 19991 219 5 75 Heartbeat classes are N: normal, S: supraventricular, V: ventricular and F: fusion •Use a multidatabase validation approach for feature selection to ensure better generalization properties of the selected feature set. 3.2 Methodology 3.2.1 ECG Databases In this chapter we used the well-known MITBIH-AR [Moody and Mark, 2001] for training and testing purposes. Additionally, the MITBIH Supraventricular Arrhythmia database (MITBIH-SUP) [Mark et al., 1990] and the St. Petersburg Institute of Cardiological Technics (INCART) database were used for evaluation and testing purposes, in order to assess the generalization achieved by the classification models developed in the MITBIHAR. All databases are freely available on Physionet [Goldberger et al., 2000] and their details were summarized in Table 3.1 and Section 2.1. 3.2.2 ECG preprocessing The ECG recordings of the MITBIH-SUP and INCART databases were first resampled to 360 Hz, which is the sampling frequency of the MITBIH-AR. This was performed with a tenth order lowpass FIR filter without observing any notorious distortion (resample 3.4. DISCUSSION AND CONCLUSIONS 77 Table 3.5: Performance comparison between the model selected in Table 3.2 and the reference classifier [de Chazal et al., 2004] separating all AAMI classes in DS2of MITBIHAR. Truth de Chazal et al. [de Chazal et al., 2004] Algorithm f n s v Total F 347 33 1 7 388 N 3509 38444 1904 303 44160 S 16 173 1395 252 1836 V 176 117 321 2504 3118 Total 4048 38767 3621 3066 49502 Truth Model selected in Table 3.2 Algorithm f n s v Total F 370 11 2 5 388 N 8031 34270 1807 80 44188 S 28 124 1403 280 1835 V 321 46 182 2669 3218 Total 8750 34451 3394 3034 49629 Performance Fusion Normal Suprav. Ventr. Total calculation mode Classifier S P+S P+S P+S P+A S P+ Imbalanced 1 95 4 78 99 76 41 83 88 78 83 58 2 89 9 87 99 75 39 80 81 86 83 57 Balanced 1 95 76 78 88 76 88 83 83 83 83 84 2 89 86 87 80 76 84 80 83 83 83 83 1: [Llamedo and Martínez, 2011a] 2: [de Chazal et al., 2004] Both models were trained in DS1of the same database. First the confusion matrices for both models are shown, and below the classes and total performances are summarized. The performances are expressed in percentages for both, balanced and imbalanced class presence in the dataset. 78 CHAPTER 3. AUTOMATIC ECG HEARTBEAT CLASSIFICATION Table 3.6: Confusion matrix as a result of separating all AAMI2 classes in the INCART database. Truth Algorithm n s v’ Total N 140983 10576 1958 153517 S 84 1660 214 1958 V’ 644 3007 16559 20210 Total 141711 15243 18731 175685 Performance Normal Suprav. Ventr. Total calculation mode Dataset S P+S P+S P+A S P+ Imbalanced DS2MITBIH-AR 95 98 77 39 81 87 93 84 75 INCART 92 99 85 11 82 88 91 86 66 Balanced DS2MITBIH-AR 95 79 77 88 81 88 84 84 85 INCART 92 92 85 80 82 87 86 86 86 By recording DS2MITBIH-AR 95 83 61 73 75 86 94 80 82 INCART 93 90 64 66 71 86 91 79 85 The model used is the selected from Table 3.2, trained in DS1of the MITBIH-AR database. The performances are expressed in percentages, and grouped by performance calculation mode for easy comparison with the results obtained from DS2of MITBIH-AR (Table 3.3). kx Mand ky M; which are described in Table 3.4. As it can be noted, the selected features are computed without exception from time interval measurements. This could be explained given that the used databases do not always include the same pair of ECG leads in each recording. Therefore the classification performances of features which are calculated from amplitudes are heavily degraded. The directional features (like the V CGφ) were also probably affected by this fact, even if the clinical importance of this kind of features is well-known by cardiologists [Taylor, 2002]. In contrast, intervals seem to retain the classification ability with independence of the pair of leads chosen. The first four features in the model are clearly connected to the evolution of heart rhythm, while the other four can be understood as surrogate measurements of the QRS width, and therefore the QRS morphology. As a result, the model found has the evident advantage of a lower size, which results in a computational saving and lower error in the parameter estimation during the training phase. In addition, it only relies on the QRS fiducial point detection, making the classifier model robust to degraded signals where the delineation of the ECG waves is not reliable. It is worth noting that the performance achieved by the reference classifier [de Chazal et al., 2004] in the union of training and validation dataset (Table 3.2) is lower for all classes than the obtained in the final performance reported in Table 3.3. The same phenomenon happens with the suggested model in a smaller degree, with the exception of the supraventricular performance. This phenomenon was also reported in [de Chazal et al., 3.4. DISCUSSION AND CONCLUSIONS 79 2004], obtaining better performance in the test set than in the training set. These results suggest that DS2dataset may not be a good data sample to measure the actual performance of a classifier. To avoid this bias in the actual performance, it may be convenient in future works that the final performance estimation is computed applying other methodologies or redefining the test dataset. One reason that could be biasing the results in DS2 is the different amount of examples by recording for the supraventricular class. As can be seen in Table 3.7, recordings 232 and 222 concentrate the majority of the examples for the supraventricular class, which means that failing in these recordings impacts considerably to the S class performance. For this reason, the average performances presented in Table 3.7 could also be of importance since each recording or subject is equally weighted in the average. The results presented in [de Chazal and Reilly, 2006], where the automatic classifier of [de Chazal et al., 2004] is assisted by a local expert to improve its performance, are also compared in Table 3.3. This suggests that a similar approach of combining the knowledge of a local expert with our model, could also lead to a comparable improvement in the baseline performance. An additional assessment of the suggested model classifying the four AAMI (N, S, V and F) classes is presented in Table 3.5. The results verify the validity of the model achieving slightly lower performance than the results presented in [de Chazal et al., 2004]. It must be noted that the model presented in this work was optimized for the AAMI2 labeling (N, S and V’), and the classifier is mainly misclassifying normal heartbeats as fusion, as shown in Table 3.5. The results in Table 3.6 suggest that the selected features have good generalization capability when evaluating the performance in heartbeats not considered during the development phase, as the ones from the INCART database. The imbalanced performance is comparable for all classes except the supraventricular where a decrease in the P+occurred. This could be explained by an increased class imbalance in the INCART database which is about 75-to-1, while in MITBIH-AR is 22-to-1 approximately. This is confirmed by the balanced results (equivalent to a class balance of 1-to-1) in the same table, where the performance figures are very similar. The validity of the generalization capability of the proposed model, is somehow restricted to the available data, and should be corroborated in future works by including new databases in the analysis or other methodologies. Despite this limitation, the degree of generalization of the suggested model is expected to be better than models obtained considering only the MITBIH-AR database. One limitation of the presented approach is the Gaussian assumption of the data imposed by the classifier, since many features were observed not to fulfill this requirement. Despite this evident limitation, the linear decision regions in the feature space defined by the LDC-C allowed us to select those features which inherently provide better classification performance. Considering the proposed classifier and feature model as a reference for future improvements, the effect of the lack of Gaussianity can be mitigated using more 80 CHAPTER 3. AUTOMATIC ECG HEARTBEAT CLASSIFICATION Table 3.7: Detailed results grouping by recording (or subject). Number of beats Normal Supraventr. Ventricular Totals Rec N S V’ S P +S P+S P+A S P+ 100 2235 33 1 100% 77% 70% 100% 100% 100% 100% 90% 92% 103 2079 2 0 99% 50% 0% 0% – – 99% 50% 25% 105 2524 0 41 97% 100% – – 51% 94% 96% 74% 97% 111 2120 0 1 99% 100% – – 100% 99% 99% 100% 100% 113 1785 6 0 99% 100% 100% 99% – – 99% 100% 100% 117 1531 1 0 100% 100% 100% 100% – – 100% 100% 100% 121 1858 1 1 99% 100% 100% 99% 100% 100% 99% 100% 100% 123 1512 0 3 100% 100% – – 0% 0% 100% 50% 50% 200 1733 30 826 96% 58% 27% 81% 92% 91% 94% 72% 77% 202 2059 55 20 67% 87% 93% 56% 50% 87% 68% 70% 77% 210 2419 22 205 94% 86% 91% 81% 69% 87% 92% 85% 85% 212 2745 0 0 100% 100% – – – – 100% 100% 100% 213 2637 28 582 100% 63% 46% 100% 44% 47% 89% 63% 70% 214 2000 0 257 97% 100% – – 94% 98% 97% 96% 99% 219 2080 7 65 86% 46% 0% 0% 82% 100% 86% 56% 49% 221 2029 0 396 93% 99% – – 99% 100% 94% 96% 100% 222 2271 208 0 72% 92% 89% 76% – – 73% 81% 84% 228 1685 3 362 100% 60% 33% 84% 93% 100% 99% 75% 81% 231 1565 1 2 98% 49% 0% 0% 50% 100% 97% 49% 50% 232 397 1381 0 100% 90% 78% 100% – – 83% 89% 95% 233 2227 7 841 100% 92% 71% 90% 83% 74% 95% 85% 85% 234 2697 50 3 100% 78% 72% 100% 100% 100% 99% 91% 93% Average 44188 1835 3606 95% 83% 61% 73% 75% 86% 94% 80% 82% Gross 95% 98% 77% 39% 81% 87% 93% 84% 75% For the model selected in Table 3.2 separating all AAMI2 classes in DS2of MITBIH-AR, following AAMI recommended performance measures. Average statistics gives each subject equal weight. Gross statistics weight each heartbeat equal. complex classifiers, like ANN’s or mixture of Gaussians. These classifiers allow more complex decision regions in the feature space, retaining details of the training data which may improve the classification performance. Despite the improved results presented in this work, there is still room for improvement in the field since the Sand P+for the supraventricular class are of 77% and 39%, and for the ventricular class (though better) are of 81% and 87%. These results suggest that other features, classifiers or meta-classifier strategies (like local expert assistance) should be developed in order to improve the performance, specially in the supraventricular class. 3.A Detailed Results In this section we present the performances by recording that the AAMI EC57 standard suggests for performance comparison. 3.A. DETAILED RESULTS 81 Table 3.8: Detailed results grouped by recordings in the INCART database Number of beats Normal Supraventr. Ventricular Totals Rec N S V S P+S P+S P+A S P+ I01 2410 0 344 95% 100% – – 100% 95% 95% 98% 98% I02 2442 0 229 94% 100% – – 44% 96% 90% 69% 98% I03 2322 2 125 80% 61% 0% 0% 99% 59% 81% 60% 40% I04 2267 16 138 64% 55% 56% 55% 62% 77% 64% 61% 62% I05 1517 0 256 70% 74% – – 27% 48% 64% 49% 61% I06 2434 48 9 100% 84% 81% 100% 100% 100% 99% 94% 95% I07 2637 65 1 100% 94% 94% 100% 100% 100% 100% 98% 98% I08 1775 2 350 94% 100% 100% 65% 53% 100% 87% 82% 88% I09 2953 0 41 80% 97% – – 88% 100% 81% 84% 99% I10 3596 0 83 95% 96% – – 96% 96% 95% 96% 96% I11 2099 0 4 95% 100% – – 75% 100% 95% 85% 100% I12 2797 1 8 76% 100% 100% 50% 25% 100% 76% 67% 83% I13 1791 0 230 87% 100% – – 3% 100% 77% 45% 100% I14 1799 0 64 94% 100% – – 0% 0% 91% 47% 50% I15 2629 0 3 96% 100% – – 100% 100% 96% 98% 100% I16 1518 0 2 100% 100% – – 0% 0% 100% 50% 50% I17 1643 0 27 98% 100% – – 0% 0% 96% 49% 50% I18 2660 0 422 76% 90% – – 89% 100% 77% 83% 95% I19 1212 0 849 95% 100% – – 97% 100% 95% 96% 100% I20 2358 179 111 86% 94% 99% 70% 67% 100% 86% 84% 88% I21 2070 104 8 95% 98% 98% 77% 75% 100% 95% 89% 92% I22 2814 124 186 76% 93% 98% 69% 76% 100% 77% 83% 87% I23 2190 0 13 100% 100% – – 46% 100% 99% 73% 100% I24 2562 0 6 99% 100% – – 83% 100% 99% 91% 100% I25 1702 2 5 94% 100% 100% 49% 0% 0% 94% 65% 50% I26 1496 7 4 99% 57% 100% 99% 25% 100% 99% 75% 85% I27 1883 0 719 99% 100% – – 100% 100% 99% 100% 100% I28 1710 0 4 96% 100% – – 50% 100% 96% 73% 100% I29 1833 0 783 82% 100% – – 57% 99% 75% 70% 100% I30 1703 0 755 95% 100% – – 51% 100% 82% 73% 100% I31 1843 0 1363 92% 100% – – 79% 98% 86% 86% 99% I32 1559 0 57 97% 100% – – 0% 0% 94% 49% 50% I33 1244 589 1 98% 97% 87% 98% 100% 92% 95% 95% 96% I34 1426 536 0 97% 97% 96% 97% – – 97% 97% 97% I35 3198 0 475 72% 86% – – 79% 100% 73% 76% 93% I36 3445 0 462 81% 86% – – 87% 99% 81% 84% 93% I37 2007 0 452 98% 99% – – 99% 100% 98% 99% 100% I38 2153 0 542 100% 99% – – 90% 100% 98% 95% 100% I39 1459 0 313 100% 100% – – 96% 100% 99% 98% 100% I40 2566 6 92 82% 83% 67% 95% 98% 75% 83% 82% 84% I41 1622 5 1 99% 62% 40% 28% 0% 0% 99% 46% 30% I42 1544 0 1561 94% 100% – – 85% 100% 90% 90% 100% I43 1084 0 1121 98% 100% – – 94% 100% 96% 96% 100% I44 1801 8 683 100% 89% 88% 99% 99% 100% 100% 96% 96% continues on the next page For the model selected in Table 3.2 separating all AAMI2 classes, following AAMI recommended performance measures. 82 CHAPTER 3. AUTOMATIC ECG HEARTBEAT CLASSIFICATION concluded from previous page Number of beats Normal Supraventr. Ventricular Totals Rec N S V S P+S P+S P+A S P+ I45 1434 0 491 99% 100% – – 99% 100% 99% 99% 100% I46 2230 1 425 86% 97% 100% 87% 96% 100% 88% 94% 95% I47 1857 1 92 93% 48% 0% 0% 47% 100% 91% 47% 49% I48 2117 2 236 99% 100% 100% 99% 100% 100% 99% 100% 100% I49 2117 0 27 89% 100% – – 100% 100% 89% 95% 100% I50 2992 0 4 95% 100% – – 100% 100% 95% 98% 100% I51 1968 3 802 73% 52% 33% 55% 100% 100% 81% 69% 69% I52 1608 0 137 100% 100% – – 100% 100% 100% 100% 100% I53 1149 0 1109 99% 100% – – 100% 100% 100% 100% 100% I54 2338 1 22 95% 91% 0% 0% 77% 43% 94% 57% 45% I55 2145 1 17 100% 100% 0% 0% 94% 48% 100% 65% 49% I56 1669 26 7 97% 93% 31% 64% 86% 58% 96% 71% 72% I57 2839 3 24 91% 69% 67% 89% 92% 100% 91% 83% 86% I58 2310 0 12 99% 92% – – 75% 100% 99% 87% 96% I59 2064 0 81 77% 100% – – 5% 99% 75% 41% 100% I60 2472 0 0 87% 100% – – – – 87% 87% 100% I61 1450 1 0 100% 100% 100% 100% – – 100% 100% 100% I62 1451 9 807 96% 71% 67% 60% 43% 80% 77% 69% 70% I63 1844 1 146 100% 74% 100% 86% 49% 100% 96% 83% 87% I64 1883 0 26 100% 72% – – 31% 100% 99% 66% 86% I65 2271 5 386 86% 59% 40% 61% 89% 100% 86% 72% 73% I66 2136 1 200 95% 49% 0% 0% 56% 100% 92% 50% 50% I67 2435 5 532 93% 54% 20% 73% 99% 100% 94% 71% 76% I68 2479 2 161 100% 67% 50% 100% 100% 100% 100% 83% 89% I69 1997 1 168 100% 99% 0% 0% 99% 50% 100% 66% 50% I70 1537 126 0 100% 100% 1% 100% – – 92% 51% 100% I71 1631 35 0 90% 97% 94% 91% – – 90% 92% 94% I72 1872 8 386 99% 86% 88% 94% 92% 100% 98% 93% 93% I73 1888 32 70 100% 96% 97% 72% 61% 100% 98% 86% 89% I74 2079 0 322 100% 79% – – 73% 100% 96% 87% 90% I75 1482 0 618 100% 95% – – 95% 100% 98% 98% 98% Average 153517 1958 20210 93% 90% 64% 66% 71% 86% 91% 79% 85% Gross 92% 99% 85% 11% 82% 88% 91% 86% 66% For the model selected in Table 3.2 separating all AAMI2 classes, following AAMI recommended performance measures. Average statistics gives each subject equal weight. Gross statistics weight each heartbeat equal. Chapter 4 Extensions to the Automatic Classifier 4.1 Introduction In this chapter we study two improvements to the original classifier developed in the previous chapter. The first one allows the classification of recordings of an arbitrary number of leads, while the second one explores the utility of more complex classifiers, as neural networks. 4.2 Multilead classification The room for improvement in the field of heartbeat classification, together with the availability of 3and 12-lead Holter devices makes necessary the development of algorithms capable of exploiting the increase of recorded information. Moreover, the St. Petersburg Institute of Cardiological Technics 12-lead Arrhythmia Database (INCART) has become recently freely available on Physionet [Goldberger et al., 2000], making possible the evaluation of multilead heartbeat classifiers in a comparable way. The objective of this study is to find an effective way to include morphologic information present in multilead ECG signals. For that purpose, we compare several multilead classification strategies against the reference two-lead classifier that we developed in [Llamedo and Martínez, 2011a]. We assess the improvement in classification performance as well as the generalization capability to other databases not considered during the development. The main novelty presented in this work is the generalization of the model developed in [Llamedo and Martínez, 2011a] to an arbitrary number of leads. 4.2.1 Material and methods In this study we used the well-known MITBIH Arrhythmia database (MITBIH-AR) [Moody and Mark, 2001] and other public databases already described in Section 2.1. 83 84 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER Table 4.1: Databases used in this work. Heart beats classes are N: normal, S: supraventricular, V: ventricular, F: fusion, and Q: unknown. Database N S V F Q #Rec INCART 153651 1959 20005 219 6 75 Biosigna 286246 1326 2541 0 0 56 MITBIH-AR 90089 2779 7007 802 15 44 MITBIH-SUP 162271 12195 9940 23 79 78 Totals 692257 18259 39493 1044 100 253 All public databases are available on Physionet [Goldberger et al., 2000] and their details were summarized in Table 4.1. For all databases the AAMI recommendations for class-labeling were adopted (Section 4.2 in [AAMI-EC57, 1998–2008]). Additionally, we used a private database called Biosigna. This database was developed at Biosigna GmbH, and consists of 56 recordings containing a broad set of pathologies. Each recording is one hour length, sampled at 500 Hz with an amplitude resolution of 12-bit over a range of 10mV. The recordings were manually annotated by experienced annotators. More detailed information about this database can be found in [Fischer et al., 2008]. For the preprocessing of the ECG recordings, we used the same methodologies described in Section 2.3.1 of the previous chapter. We follow the results obtained in Chapter 3 where we developed a heartbeat classifier with good generalization capability, using rhythm and morphological features together with a linear classifier compensated for class-imbalance, called LDC-C. This classifier has the possibility of weighting the class contribution to the covariance matrix (equation (2.25)) during the training, to deal with the class imbalance seen in Table 4.1. The equations of the LDC-C were described in Section 2.4.2. The features obtained from the sequential floating feature selection (SFFS) in [Llamedo and Martínez, 2011a] are shown in Table 4.2. As the rhythm features of the model do not depend on the number of available leads, the first four features in Table 4.2 remain the same. Therefore we will focus the analysis on those features describing heartbeat morphology, which are the ones that can be improved by the addition of new leads. The morphology features used in the model are the first zero-crossing (kL Z) and maximum position (kL M) of the autocorrelation sequence of the ECG DWT at scale 4 for each lead L(see Section 2.4.1 for details). Both features were calculated in four sets of leads to study the most suitable way of integrating the information from all leads: 1. The first strategy consists of just including kL Zand kL Mfrom all available ECG leads and is referred as 12L (or 3L when only 3 leads are available), resulting in two morphology features per lead. 2. The second strategy computes kL Zand kL Mfrom the three vectocardiogram (VCG) 4.2. MULTILEAD CLASSIFICATION 85 Table 4.2: Features used in the model obtained in [Llamedo and Martínez, 2011a] only for two-lead recordings. Feature Description ln(RR[i]) Current RR interval1 ln(RR[i+ 1]) Next RR interval1 ln(RR1)Average RR interval in the last minute1 ln(RR20)Average RR interval in the last 20 minutes1 kx ZZero-cross position of the WT autocorrelation signal in lead 12 ky ZZero-cross position of the WT autocorrelation signal in lead 22 kx MMaximum position of the WT autocorrelation signal in lead 12 ky MMaximum position of the WT autocorrelation signal in lead 22 1See Figure 2.6 2See Figure 2.10 leads X, Y and Z, transformed from 12L by the Dower matrix. This strategy can only be performed in 12-leads recordings. 3. In the third strategy, referred as ECG-PCA, we apply principal component analysis (PCA) to the available ECG leads, then the morphology features are computed from the two most important components. 4. Finally for the fourth strategy, called WT-PCA, the PCA is applied not to the ECG, but to the fourth scale of the DWT (W4s(k)), and the two morphological features kL Zand kL Mwere calculated from the principal components. The last two strategies are the result of projecting all ECG leads or its fourth scale DWT (Ws 4(k)), into the two most important basis of a principal component analysis (PCA). The PCA consists of finding an orthogonal linear transformation such that the first component of the transformation comes to lie in the direction of the greatest variance of the data. For our multilead ECG signal S, with each lead of Nsamples as a column resulting in a N×Lmatrix, the transformation is defined by R=S·P,(4.1) where Pis the matrix which defines the linear transformation and Ris the transformed ECG. The matrix Pis the result of computing the eigenvectors of the sample covariance matrix C, calculated as described in the next subsection, P−1CP =Q,(4.2) being C=STS.(4.3) 86 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER 0 200 400 -200 0 0 200 400 600 -100 0 100 II V1 Lead II Lead V1 PCA1 PCA2 PCA2 PCA1 Signal interpretation Geometrical interpretation First component (PCA1) Second component (PCA2) Rotation Rotation Figure 4.1: Toy example where a two-lead ECG excerpt is also interpreted geometrically, and PCA transformation is performed. Note the rotation involved in the PCA transformation. The associated diagonal matrix of eigenvalues Q, is related to the importance of each eigenvector of Pin the sense of the variance of the original data S. That is, if we extract only the first two components of R, we are retaining a fraction of the original QRS complex variance which is typically above the 90% for normal sinus rhythm heartbeats. In the rest of this work we will refer to perform PCA to calculate the two most important projections of a multilead signal. In the toy example presented in Figure 4.1, a two-lead ECG excerpt is interpreted as a multilead signal (typical interpretation), and as a bivariate collection of unrelated samples (geometrical interpretation). The latter interpretation neglects the time relation between each sample of the sequence. This interpretation is useful to visualize the main directions of variation of the signal, and the rotation involved in a PCA transformation. It is useful to think the columns of the Pmatrix, as weights of a linear transformation. Each element of the column vectors, are related to the contribution of each lead, as can be seen in the center of Figures 4.2 and 4.3. The variance shown in the right part of the figures, is the ratio between the first two eigenvalues of Qand the sum of all eigenvalues. Note the similar weight patterns for the normal and supraventricular classes, which are both similarly conducted through the ventricles. These patterns depend on the heartbeat morphology as can be seen for the ventricular and fusion classes in Figure 4.3. The PCA is performed for each heartbeat in a 160 ms window centered at the QRS complex detection sample, or fiducial point (PCA window in Figure 4.4). Then, the multilead signal is projected into the PCA basis in a wider window starting 130 ms before and ending 200 ms after the QRS fiducial point. As it is known from previous works 4.2. MULTILEAD CLASSIFICATION 93 model to other databases for different number of leads. For all databases available for a given number of leads, we assessed the performance using all possible pairs of different databases as train and test sets. We also evaluated the crossvalidated performance within each database. To have an upper-bound reference, we additionally assessed the performance of the model when trained and tested in the same database. This optimistically biased performance serves as an upper bound, and represents the performance of the model if the distributions of the examples in both training and test datasets were identical. These results, grouped by test database, are presented in Table 4.5 for databases with 12, 3 and 2 leads. Results show that the reference model extended with the selected WT-PCA multilead strategy presents good generalization properties for 3 and 12 leads, while certain degradation is observed when using only two leads. 4.2.3 Discussion and conclusions In this work we have adapted and improved a two-lead heartbeat classifier by including the additional morphology information present in multilead recordings, like those of INCART database. We followed the concept of the morphology features assessed in [Llamedo and Martínez, 2011a], but calculating these features in sets of leads obtained by following different lead transformation strategies. The simplest strategies consisted in computing the features kL Zand kL Min all available leads (12L/3L), and in the derived orthogonal leads (VCG). The other two strategies apply PCA to the ECG or its WT previously to the morphology feature computation. The results suggest that strategies using PCA performed better. Moreover the WT-PCA strategy obtains the best improvements with respect to the two-lead classifier obtained in [Llamedo and Martínez, 2011a], either in recordings of 12 or 3 leads. This can be explained because in WT-PCA, the PCA is calculated in Ws 4(k)where most of the noise and other components not related to the QRS have been filtered out [Martínez et al., 2004], and therefore PCA provides a better representation of the multilead evolution of the QRS complex. It must be remembered that although both PCA and WT are linear transformations, the eigenvector calculation is not linear, and therefore sets ECG-PCA and WT-PCA differ in how the multilead signal is projected into the principal components. Another important improvement achieved with the WT-PCA strategy is the robustness against lead misplacement or recordings with undocumented leads. Tables 4.3 and 4.4 show that only the WT-PCA strategy showed the largest observed performance improvement with respect to the two-lead reference model developed in [Llamedo and Martínez, 2011a]. Results in Table 4.5 confirm the generalization capability of the model using the selected WT-PCA strategy to the rest of the studied databases. For the case of 12-lead recordings, as those included in the INCART and Biosigna databases, results show very good generalization for both databases since the performance is slightly lower than the biased performance when training in the other database. In the same table, almost the 94 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER same figures can be seen for the case of 3-lead recordings. The last results presented in Table 4.5 are for databases of two leads. It can be noted that the inter-database dispersion in the performance increased, probably because of the heterogeneity of the databases considered. The results suggest that databases INCART and MITBIH-AR share similar distributions in the feature space, since both obtained the maximal reciprocal performance. Other interesting aspect regarding the MITBIH databases is the lower performance obtained even for the biased case. This fact evidences the diversity of patients and ECG contamination (noise, lead disconnection and misplacement) included in these databases; and therefore the limitation of our classifier to model the data and achieve higher performances. Certainly the biased performance can be thought as a metric of how difficult to classify is a database by a given classification model, the closer to 100%, the easier. This last result reinforces the importance of evaluating the performance of a classifier in several databases. One advantage of the proposed approach is that it can be used for an arbitrary number of leads, because after the PCA we only retain the two most important components for the morphological feature calculation. These components are calculated specifically for the QRS complex, and in the Ws 4(k)signal (with a band between 11.25 and 22.5 Hz), typically where the ECG presents high SNR. However in case of a large-scale artifact during the QRS complex (as a lead disconnection), the PCA calculation would be corrupted, being this the main limitation found for this approach. This problem is addressed by the robust MCD algorithm to compute the covariance matrix. The performance improvement with respect to [Llamedo and Martínez, 2011a] is however moderate, probably because the automatic classification approach is close to the performance limit achievable with the current classification model. The worst aspect of performance remains classification of supraventricular ectopic beats, where further study is needed. Regarding the ventricular class, techniques of patient adaptation, as described in [de Chazal and Reilly, 2006], will be presented in Chapter 5. These results represent an improvement in performance with respect of the previous two-lead classification model, concluding that the adequate addition of multilead information allows the performance improvement of a heartbeat classifier. 4.3 Neural network classifier In this section, the objective is to improve the classification method used. It is not difficult to understand that a simple classifier as the LDC, which is capable of dividing the hyperspace with hyperplanes, is a suboptimal solution for a complex problem such as the classification of heartbeats. The design of an automatic classifier based on neural networks shares the same methodology described in Chapter 3, but were used in a different set of features. The results presented in this section were published in [Mar et al., 2011], 4.3. NEURAL NETWORK CLASSIFIER 95 in collaboration with the Institut für Biomedizinische Technik, at the Dresden University of Technology, Germany. 4.3.1 Feature Sets We used a set of 71 features, divided (as shown in Table 4.6) into the following categories: •Temporal features, which already proved, in different studies, to be the most relevant [Lannoy, 2010]. This category includes heart-rate features, which were the only features computed just once for both channels, and features obtained from the segmentation information yielded by ecgpuwave [Laguna et al., 1994]. •Morphological features, also previously assessed as being of great relevance, made the bulk of the feature set. Direct samples from the ECG signal, and computed measurements such as area, power or extrema, were included. •Statistical features completed the feature set, including different order momentbased indexes and histogram variance. Unlike in temporal features, which were acquired from time domain signals exclusively, the DWT of the ECG signal was used to obtain some of the features belonging to the morphological and statistical categories. DWT features were based on the same heartbeat intervals as the features obtained from the time domain, but using the scales 2, 3, 4 and 5 of the DWT ECG signal. Details about the DWT can be found in Section 2.3.2. Finally, all features were individually normalized, by computing the necessary scaling to make the features from DS1 signals be mean 0 and variance 1, and normalizing the corresponding features from DS2 with the obtained scaling factors. 4.3.2 Feature Selection The SFFS procedure is the same described in Section 2.4.6, but the optimization criterion was different. We propose a new performance measure which tackles specifically the problem of providing in a single value complete information about how good an ECG classification has been. To this end, this new index was chosen to be a combination of two values defined in Section 2.4.5: the jindex, which specifically evaluates the discrimination of the most important ECG arrhythmias (S and V beats), j=SS+SV+P+ S+P+ V(4.6) and the Kappa (κ) index, which globally evaluates the confusion matrix [Cohen, 1960]. This index, despite having been proposed as an evaluation coefficient several decades ago, and its potential convenience, had never been applied before in the heartbeat classification. 96 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER Table 4.6: Features distributed by categories. Features Time Domain DWT Domain* Temporal Previous RR, current RR, average RR, average RR of the last 10 beats, QRS duration, T-Wave duration, P-Wave flag. Morphological Downsampled (10 samples) QRS, downsampled (9 samples) T-Wave, QRS area, QRS power, QRS max, QRS min, QRS Max-Min ratio, peak width at 70% Max, peak slope, beat area, beat power, beat max, beat min, beat Max-Min ratio. Max(3,4,5), Min(3,4,5), difference between Max and Min (3,4,5), distance (in samples) between Max and Min (3,4), power(2,3,4,5), power ratio(3-2,4-3,5-4). Statistical QRS variance, QRS skewness, QRS kurtosis, QRS histogram (20 slots) variance, beat mean, beat variance, beat skewness, beat kurtosis, beat histogram (20 slots) variance. Mean (3,4), standard deviation(3,4), skewness (3,4). *Numbers between parenthesis represent the scales at which the feature was obtained. From its definition, κ=PC i=1 nT ii −PC i=1 Di NT−PC i=1 Di (4.7) where Di=(NiPi) NT (4.8) is known as the weighted detections, it can be seen that it evaluates the global quality of the classification: like the multiway accuracy, it also represents a complete evaluation of the confusion matrix (in a single value and weighting each beat equally), but it is much less influenced by the class imbalance. The resulting combined index, which we named jκ index (Ijκ), takes into account the misclassification and the imbalance present between all the considered classes, thanks to the included κindex, and at the same time emphasizes the discrimination of the most important arrhythmias (S and V), thanks to the jindex (Ij). Ijκ =w1κ+w2Ij(4.9) As jtakes values in the 0-4 range and κin the 0-1 range, w1was set to 1 /2and w2to 1 /8, so that both factors influence equally the overall result. Consequently Ijκ takes values between 0 and 1, where 1 indicates perfect classification. 4.3. NEURAL NETWORK CLASSIFIER 97 Table 4.7: DS1 Division Scheme for MLP Evaluation Dataset N S V F Total Eval1 DS1 23379 442 1936 40 25797 Eval2 DS1 22283 501 1842 373 24999 Total (DS1) 45662 943 3778 413 50796 Eval1 DS1 comprises data from the recordings 109, 114, 118, 119, 124, 201, 203, 205, 215, 220 and 223 Eval2 DS1 comprises data from the recordings 101, 106, 108, 112, 115, 116, 122, 207, 208, 209 and 230. 4.3.3 Multi-Layer Perceptron The multilayer perceptron (MLP) belongs to the class of supervised learning networks, on which the discriminative power is gained through a preliminary learning phase, where labeled examples are presented to the network. The most common training strategy, also used in the present study, is the backpropagation (BP) algorithm [Rumelhart et al., 1986]. It works by computing the error between the returned and the known, desired output, employing it to adjust the MLP weights. Although the training process requires a rather long time, the implementation and execution of a trained MLP is very simple, making this paradigm very suited too for classification on ambulatory settings. On the other hand, its characteristics make this paradigm very inadequate to guide the feature selection (FS) process: The random initialization makes MLPs’ results not constant, which renders the FS procedure unreliable if only one MLP is evaluated for each tested subset. The unreliability could be overcome by training many MLPs for each tested subset, and performing statistical analysis to obtain a result that would lead to the next step in the FS process. Unfortunately, due to the many subsets tested by the SFFS procedure, plus the relatively long time that training each MLP requires, the time and resources that a reliable MLP-SFFS procedure would take are beyond our computing capabilities. Therefore, in the present study the MLP paradigm was only applied to classify ECG arrhythmias with the feature set obtained from the SFFS. The values of the different parameters governing the MLP were determined by applying 2-fold cross-validation on DS1, training with one fold the MLP parameterized with the desired combination and evaluating with the remaining one, and vice versa, averaging the results. Again, this sub-division was performed inter-patiently in such a way that all heartbeat classes were similarly represented in each of the folds, as shown on Table 4.7. MLPs with a single hidden layer of 25 neurons were used, trained with a learn rate of 0.25 and a momentum of 0.03 to avoid getting stuck into local minima. The number of training cycles was chosen to be the one for which the mean results from the 2-fold evaluation began to get worse, i.e. when symptoms of over-learning appeared. 98 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER 4.3.4 Classifier Combination As mentioned above, information from both ECG leads was considered throughout the whole study. Except for the heart-rate ones, all features were obtained separately for each lead, and the two resulting feature sets applied independently to perform classification. The posterior probabilities obtained after classification with each feature set were then combined using the Bayesian product integration scheme [de Chazal et al., 2004] P(i|x) = QL l=1 Pl(i|x) PC c=1 QL l=1 Pl(c|x),(4.10) for L= 2 leads and C classes. Finally, each heartbeat was labeled with the class with higher posterior probability after the combination. 4.3.5 Results In the present study aiming at reducing the classifiers complexity as much as possible, but without making its performance worse, we selected as the most suited for our purposes the classifier with the smallest number of features which achieved at least the performance obtained with the original feature set on DS1. The reduced subset accomplishing this criterion contains 9 features which includes Previous RR, current RR, RR average, beat min, beat max, QRS max-min ratio, peak slope, max-min Difference on DWT scale 3. After evaluating the performance of the SFFS procedure with the matched classifier on DS1, the original feature subset and the selected from the SFFS were tested on DS2 to carry out the final evaluation of the classifier model. Additionally, these subsets were also tested on DS2 with the MLP classifier, in order to analyze the generalizing capability of the SFFS procedure in the case where the criterion function and the classifier paradigm do not match. At the same time, this analysis also tackles the suitability of the MLP for heartbeat classification, in direct comparison to the LDC classifier. In Table 4.8, complete classification description is displayed in the form of the confusion matrices for the results obtained by applying the reduced feature set with either classifier paradigm. These matrices provide insight into the beat-by-beat performance and ease future comparison attempts by other authors. For the rest of studied feature sets, results both for LDA and MLP classification are given through the considered performance measures on Table 4.8. 4.3.6 Discussion and conclusions In addition to the analysis of the SFFS process itself, it is also interesting to identify the most relevant features obtained. Inspecting the selected subset we can observe the previous RR,current RR and RR average features (RR[i−1],RR[i]and RRAin Figure 2.6) were present. This indicates the uttermost importance that heart-rate features have on ECG classification. This fact was also corroborated the model found in Chapter 3 and 4.3. NEURAL NETWORK CLASSIFIER 99 Table 4.8: Classifier Performances on DS2 Obtained for the most Relevant Studied Classifier Models Truth LDC Algorithm n s v f Total N 37384 2726 691 3260 44061 S 60 1517 237 16 1830 V 45 225 2782 156 3208 F 137 1 50 200 388 Total 37626 4469 3760 3632 49487 Truth MLP Algorithm n s v f Total N 39497 2778 771 1015 44061 S 122 1523 93 92 1830 V 104 235 2783 86 3208 F 125 6 20 237 388 Total 39848 4542 3667 1430 49487 Normal Suprav. Ventr. Fusion Total Classifier S P+S P+S P+S P+A S P+jκ Ijκ LDC 85 99 83 34 87 74 52 6 85 77 53 2.78 0.51 0.60 MLP 90 99 83 34 87 76 61 17 89 80 57 2.79 0.60 0.65 Table 4.9: Relevant Indices for AAMI standard and Inter-Patient Division Conform Studies, Including Present Study’s ones. Total Classifier A S jκ jκ LDC [Mar et al., 2011] 85 77 2.78 0.51 0.60 MLP [Mar et al., 2011] 89 80 2.79 0.60 0.65 LDC [de Chazal et al., 2004] 86 83 2.764 0.532 0.612 LDC* [de Chazal and Reilly, 2006] 94 88 3.234 0.754 0.781 LDC [Llamedo and Martínez, 2007] 80 80 2.372 0.421 0.507 SVM [Park et al., 2008] 86 76 – – – SVM [Lannoy, 2010] – 83 – – – LDC** [Llamedo and Martínez, 2011a] 78 83 2.887 0.412 0.567 *Patient adapting: Requires expert intervention. ** Feature set optimized for classification of N, S and V’ classes, where V’ class included V and F beats. 100 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER presented in Table 4.2. Regarding the results for the feature set obtained, which contained a much larger number of features, including many statistical ones, it is remarkable to observe that the 6 non-heart-rate features are all morphological features too, and even more noteworthy, that 5 of them represent extrema. This fact suggests that extrema values possess the highest discriminative power among all morphological features. In spite of the large number of studies in which the MLP classifier paradigm has been applied for ECG classification, none among them could be found in which the results were evaluated in conformance with the AAMI standard and inter-patient dataset distribution. Yet, results show that, when applying reduced feature sets, the MLP can clearly outperform LDC in the task of heartbeat classification. Comparing the results of both paradigms, an improvement in the range of the 4% can be observed in the global accuracy (A)and in the global sensitivity (S). Nevertheless, it should be noted that these numbers are just orientative of the possible improvement, as, due to their random initialization, successive evaluations of the MLP with the same feature set could yield different results. Comparing the achieved performances with those of previous studies, it provides further insight on the suitability of the proposed techniques. As mentioned, this comparison can only be objectively done among those studies following the same constrains. Thus, the performances obtained in the present study on both classifier paradigms have been compared with the results of the other AAMI conforming studies which followed the inter-patient division scheme (see Table 4.9). However, the ones achieved with the MLP outperform all previous non-adapting proposed methods. These results evince that the non-linear classification capabilities of this type of MLP are extremely suitable to perform heartbeat classification, which is intrinsically non-linear too. Moreover, they also show, in the context of this study, a greater generalization capability of the MLP when compared to algorithmic methods such as LDC, suggesting that they may be a more appropriate tool for ECG heartbeat classification. 4.A Detailed Results In this section we present the confusion matrices for ease the comparison of future works. The summarized performances presented in the previous sections were calculated from these numbers, according to the methodologies described in Section 2.4.5. 4.A. DETAILED RESULTS 101 Table 4.10: Confusion matrices of the results presented in Table 4.5 for 12-leads databases. Truth biased INCART Algorithm n s v Total N 151649 1553 449 153651 S 69 1741 149 1959 V 481 1361 18382 20224 Total 152199 4655 18980 175834 Truth biased Biosigna Algorithm n s v Total N 280874 5031 341 286246 S 77 1222 27 1326 V 260 337 1944 2541 Total 281211 6590 2312 290113 Truth crossval INCART Algorithm n s v Total N 151077 2055 519 153651 S 91 1685 183 1959 V 526 1674 18024 20224 Total 151694 5414 18726 175834 Truth crossval Biosigna Algorithm n s v Total N 280555 5206 485 286246 S 83 1212 31 1326 V 266 357 1918 2541 Total 280904 6775 2434 290113 Truth Biosigna - INCART Algorithm n s v Total N 151488 1399 764 153651 S 88 1685 186 1959 V 828 1267 18129 20224 Total 152404 4351 19079 175834 Truth INCART - Biosigna Algorithm n s v Total N 280163 6586 805 287554 S 109 1189 37 1335 V 213 366 1993 2572 Total 280485 8141 2835 291461 102 CHAPTER 4. EXTENSIONS TO THE AUTOMATIC CLASSIFIER Table 4.11: Confusion matrices of the results presented in Table 4.5 for 3-leads databases. Truth biased INCART Algorithm n s v Total N 150937 1857 857 153651 S 70 1811 78 1959 V 565 1606 18053 20224 Total 151572 5274 18988 175834 Truth biased Biosigna Algorithm n s v Total N 280840 5176 230 286246 S 74 1227 25 1326 V 289 391 1861 2541 Total 281203 6794 2116 290113 Truth crossval INCART Algorithm n s v Total N 150307 2278 1066 153651 S 88 1705 166 1959 V 557 1841 17826 20224 Total 150952 5824 19058 175834 Truth crossval Biosigna Algorithm n s v Total N 279365 5567 1314 286246 S 83 1219 24 1326 V 298 400 1843 2541 Total 279746 7186 3181 290113 Truth Biosigna - INCART Algorithm n s v Total N 151516 1393 742 153651 S 89 1653 217 1959 V 1114 1778 17332 20224 Total 152719 4824 18291 175834 Truth INCART - Biosigna Algorithm n s v Total N 275435 5401 5410 286246 S 84 1181 61 1326 V 231 389 1921 2541 Total 275750 6971 7392 290113 5.2. METHODOLOGY 109 5.2.2 Heartbeats classification Following the scheme presented in Figure 5.1, the patient-adaptable algorithm includes a linear discriminant classifier (LDC) and an expectation-maximization clustering algorithm (EMC). Both LDC and EMC work independently and each performs a preliminary classification/clustering task in different feature spaces. The LDC was developed and trained as described in [Llamedo and Martínez, 2011a], while the EMC development will be described later. Finally, the heartbeat and cluster labels provided by the LDC and EMC respectively, are integrated with a voting scheme into a final heartbeat label. Three modes of operation are proposed, depending on the degree of expert assistance available in the application scenario: 1) automatic, 2) slightly assisted and 3) assisted. The algorithm performs the following procedures: a) cluster and centroid identification, b) LDC automatic classification and c) expert assistance. For the automatic mode, in each record, Kclusters and centroids are identified, corresponding to groups of similar heartbeats, while at the same time, the LDC computes the labels for each heartbeat. Then for each cluster, the algorithm tests if any label obtains a qualified majority, meaning that the most represented label exceeds the αpercent of the cluster population. In case this label exists, it is assigned to the whole cluster, superseding the LDC labels. If the qualified majority is not reached, the uncertainty is considered to be too high to change the labels, an thus the LDC labels remain unchanged. The slightly assisted mode is similar to the automatic, with the exception that in case of not finding a class with qualified majority, expert assistance is required to label the cluster centroid and propagate it to the whole cluster, ignoring LDC labels. The procedure of expert assistance is simulated by inspecting the true labels provided with each database. Finally in the assisted mode,Kclusters and centroids are identified, then the expert is required to label each centroid. The algorithm concludes assigning these labels to the rest of heartbeats in each cluster. To better understand the three modes of operation, a toy example can be found in the center of Figure 5.1. The LDC by itself makes 4 errors in cluster 1 and 3 in cluster 2. For the automatic mode, clusters 2 and 3 have majority of V and N classes respectively. Then, votation propagates centroid labels to the rest of examples within both clusters, occurring 1 mistake in cluster 2. In cluster 1 there is no qualified majority for α= 50%, so the LDC labels remain unchanged and 5 mistakes happen. The automatic mode made 1 mistake less than the LDC. For the slightly-assisted mode, only cluster 1 would be modified, by propagating the true label of the centroid S, resulting 4 errors for this cluster. Finally for the assisted mode, the result is the same as in the previous mode, four errors in cluster 1, one error in cluster 2 and no errors in cluster 3. In summary, 7 errors for the LDC, 6 for the automatic mode and 5 for the slightly and assisted modes. As it was shown, the algorithms rely heavily in the ability of the EMC to cluster the heartbeats adequately. 110 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION nvn LDC Feature vector LDC EMC LDC Feature vector LDC Feature vector EMC Feature vector EMC Feature vector EMC Feature vector s nv n nn n n n nn n n v v v v vv v v vvv v v n n n nnn s sss s Voting Expert 1) Automatic 2) Slightly assisted 3) Assisted Expert labels all centroids Heartbeats labels Heartbeats labels Heartbeats labels clusters centroids labeled centroids labels labels ECG Feature extraction Feature extraction Nn Nn Nn Nn Nn Nn Nn Nn Nn Vn Vn Vv Vv Vv Vv Vv Vv Vv Vv Vv Vv Nn Nn Nn Nn Nn Nn Nn Vn Nn Ns Sn Ss Ss Ss Ss Ss Sn Nv Sv Vv Cluster 1 (8s, 3n, 1v) Cluster 2 (1n, 12v) Cluster 3 (15n) True labels N: Normal LDC labels centroids 1 2 3 LDC when no qualified majority Expert assistance when no qualified majority S: Suprav. V: Ventricular n: Normal s: Suprav. v: Ventricular Figure 5.1: Overview of the proposed algorithm. There is a graphical description in the center of the scheme about the task carried out by each block. The toy-example in the middle is also commented in the text to better understand the three modes of operation. 5.2. METHODOLOGY 111 Table 5.1: feature model used by the automatic classifier for recordings of 2 or more leads . Feature Description ln(RR[i]) Current RR interval1 ln(RR[i+ 1]) Next RR interval1 ln(RR1)Average RR interval in the last minute1 ln(RR20)Average RR interval in the last 20 minutes1 k1 ZZero-cross position of the WT autocorrelation signal in lead 12 k2 ZZero-cross position of the WT autocorrelation signal in lead 22 k1 MMaximum position of the WT autocorrelation signal in lead 12 k2 MMaximum position of the WT autocorrelation signal in lead 22 1See Figure 2.6 2See Figure 2.10 5.2.3 Automatic classifier We follow a scheme similar to the one in [Llamedo and Martínez, 2011a, 2012a], where we developed a multilead heartbeat classifier with good generalization capability. We used a linear classifier compensated for the class-imbalance, while as feature model we adopted rhythm and morphological features computed in a multilead manner. Regarding to the classifier used, we found that linear discriminant functions were suitable for the heartbeat classification task in terms of performance and generalization capability. The details and equations of this classification model can be found in Section 2.4.2. The features used by the automatic classifier are described in Table 5.1. The morphology features kL Zand kL Mfor lead L, are calculated in the two principal ECG leads after integrating the multilead information with a principal component analysis (PCA). In Chapter 4 it was shown that WT-PCA strategy was a good strategy to include multilead morphology information. Therefore these features account for a multilead morphological description of the QRS complex. For a detailed description of the features and the multilead strategy used see Chapters 3 and 4. 5.2.4 Clustering algorithm The EMC algorithm used in this chapter is based on the mixture of Gaussians model [van der Heijden et al., 2005]. It consists of estimating the parameters of a density function p(x|Ψ) = K X k=1 πk·f(x;µk,Σk)(5.1) f(x;µk,Σk) = K X k=1 πk 1 q(2π)m|Σk|exp−1 2(x−µk)TΣ−1 k(x−µk),(5.2) 112 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION where the m-dimensional vector xis modeled by KGaussians with mixing coefficients πk, in order to retain a more realistic structure of the data. The parameter set Ψ = {πk,µk,Σk|k= 1, . . . , K}is estimated by maximum likelihood criterion. We maximize the log likelihood L(X|Ψ) = ln N Y n=1 p(xn|Ψ),(5.3) for the Nheartbeats in each recording named X={x1,..., xN}. Since there is not a closed form solution for Ψby maximizing L(X|Ψ), the well-known expectation-maximization algorithm (EM) is used to obtain the estimation equations of the parameters Ψ, which are the mixing coefficient for each cluster ˆπk=1 N N X m=1 ˆ βm,k,(5.4) the cluster mean ˆµk=1 Nˆπk N X m=1 ˆ βm,kxm(5.5) and cluster covariance matrix ˆ Σk=1 Nˆπk N X m=1 ˆ βm,k(xm−ˆµk)·(xm−ˆµk)T.(5.6) Where ˆ βm,k is known as the ownership variable, which indicates the probability of sample xmto have been generated by the k-th component ˆ βm,k =ˆπk·f(xm;ˆµk,ˆ Σk) PK j=1 ˆπj·f(xm;ˆµj,ˆ Σj).(5.7) The EM algorithm iteratively computes the weight, location and dispersion for each of the Kclusters (Eq. (5.4)-(5.6) respectively), until ˆ βm,k does not change significantly, which is equivalent to obtaining stable clusters. The interested reader is referred to [van der Heijden et al., 2005, Duin et al., 2008] for details, equations and the implementation used in this chapter. The EM algorithm guarantees the convergence, at least, to a local optimum. The mathematical demonstration of this can be found in [van der Heijden et al., 2005]. However, in Figure 5.2 an example of this is shown. The toy example is build from three Gaussian distributions with different parameters. Even for the case where the algorithm tried to find less components (K= 2), the algorithm converges. For the rest of the cases, the redundancy of the components is notorious. This leads to one of the most critical aspects of this kind of clustering algorithms, how to estimate the complexity of the data to cluster. As there is not a reliable methodology to answer this, one recommendation 5.2. METHODOLOGY 113 Feature 1 Feature 2 0.5 1 1.5 0 0.5 1 1.5 Feature 1 Feature 2 0.5 1 1.5 0 0.5 1 1.5 Feature 1 Feature 2 0.5 1 1.5 0 0.5 1 1.5 Feature 1 Feature 2 0.5 1 1.5 0 0.5 1 1.5 2 clusters 4 clusters 8 clusters 6 clusters Initial Middle Last Figure 5.2: Toy example of a non-Gaussian distribution with several amount of clusters to be found. The elliptic equiprobable contour shows the estimated Gaussian distribution component at each step. Three situations of the EM algorithm are shown: initial, middle and last iteration. Note how the components are adapted to the data. is to follow the a priori knowledge of the problem. As will be explained below in the results section, as our classification problem involves 3 classes, we will allow between 3 and 4 clusters per class, which means 9 or 12 clusters in total. An example with real data is shown in Figure 5.3. In this example two clusters are shown where the EMC grouped normal and ventricular heartbeats. In the top-right of the same figure, among the more distant heartbeats, some misclassified examples can be seen. Note the presence in the morphology details of some widened QRS complexes. 5.2.5 Feature selection for clustering Regarding the feature model used with the EMC, we followed the same feature selection procedure described in Section 2.4.6, by means of a sequential floating feature selection algorithm (SFFS) [Pudil et al., 1994, van der Heijden et al., 2005]. The complete pool 114 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION 1.6" 0.8" 0.0" 0.8" 1.6" 0 0.5 1 1.5 0.5 0.55 0.6 0.65 0.7 Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails Centroid this heartbeat 1.6" 0.8" 0.0" 0.8" 1.6" 1 0 1 2 0.4 0.6 0.8 1 Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails 10 examples far from the centroid this heartbeat 1.6" 0.8" 0.0" 0.8" 1.6" 1 0.5 0 0.5 1 1.5 2 0.4 0.5 0.6 0.7 0.8 Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails Centroid this heartbeat 1.6" 0.8" 0.0" 0.8" 1.6" 2 1 0 1 2 Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails 10 examples far from the centroid this heartbeat 0.4 0.6 0.8 1 1.6" 0.8" 0.0" 0.8" 1.6" 1 0 1 2 Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails 10 examples close to the centroid this heartbeat 0.5 0.6 0.7 0.8 1.6" 0.8" 0.0" 0.8" 1.6" Local RR inte rval evolution 0.2" 0.0" 0.2" 0.4" Morphology de tails 10 examples close to the centroid this heartbeat 1 0 1 2 0.4 0.6 0.8 1 Ventricular cluster Normal cluster Figure 5.3: Clustering algorithm applied to the recording 208 of MITBIH-AR database. Only two clusters are shown for simplicity. In the top panel one with normal heartbeats, and below one with ventricular heartbeats. Within each panel, some heartbeats sampled from each cluster, from left to right, the centroid heartbeat and the 10 closer and farther heartbeats to the centroid. The red dashed lines indicates the heartbeat position. The rhythm evolution and the morphology details are also shown below. Note some bad-clustered heartbeats in the farther examples. 5.2. METHODOLOGY 115 Table 5.2: Features used with the EMC algorithm. Feature Description ln(RR[i]) Current RR interval1 ln(RR[i−1]) Previous RR interval1 ln(PRR)Prematurity of the heartbeat1 ln(dRRL)Local RR interval variation1 ln(RR20)Mean RR interval within the last 20 minutes1 ln(S1 QRS)QRS mean wavelet scale in the first principal component2 k1 M First minimum position of the WT autocorrelation sequence in the 1st principal component3 rQRST(kM) Value of the first maximum in the QRST complex crosscorrelation sequence between WT scale 3 of the3first two principal components 1See Figure 2.6 2See Figure 2.9 3See Figure 2.10 of features consisted of 61 features, described in Section 2.4.1. For the case of clustering, instead of looking for features with generalization capability or interpatient separability, we looked for those with high intrapatient separability. This criterion was achieved by modifying the SFFS’ optimization criterion explained in Section 2.4.6, in order to find a feature model that provides as much intrapatient class separability as possible, facilitating the clusters identification. The first modification consisted in evaluating our clustering algorithm in a patient by patient fashion since this is how this algorithm will be used in practice. The second is that the performance will be evaluated in an optimistically biased fashion, described in Section 2.4.5, assuming that we know a priori the true labels of the heartbeats. The feature selection experiments were carried out in a dataset formed by the union of MITBIH-SUP with DS1 subset of MITBIH-AR [de Chazal et al., 2004]. As the SFFS performs thousands of model evaluations, this task is very demanding in processing power, specially for the random and iterative nature of the EMC. For this reason we replaced the EMC, only for the feature selection task, for a classifier based on mixture of Gaussians (MoG), which uses the same algorithm used for cluster discovery. The classifier based on MoG (MoGC) models each AAMI class with KGaussian distributions, in contrast with the LDC-C that models each class with a mean vector and a pooled covariance matrix for all classes (Eq. (2.22) and (2.23) respectively). The MoGC uses during training the EM algorithm for the estimation of the Gaussian components. This modification results, first, in moving through deterministic paths through the performance surface evaluated with the SFFS, and second in easing the EM iteration since the heartbeat labels are known a priori. As a result, a model of 8 features was obtained. This model also includes a description of the rhythm and morphology of heartbeats as shown in Table 5.2. Among the rhythm features used in the model, some of them have already been used in previous works. The prematurity of a heartbeat 116 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION PRR[i] = RR[i] Pi+1 k=i−1RR[k],(5.8) measures how anticipated is a heartbeat respect to the previous and next RR interval. The local RR interval variation is defined as dRRL[i] = Pi+1 k=i−1|dRR[k]|, where dRR[i] = RR[i]−RR[i−1]. One of the morphology related features is the wavelet scale where the QRS complex is mostly projected. It is known that fast evolving signals, as a normal heartbeat, tend to be projected in lower wavelet scales or contains higher frequency components. The QRS center scale for each lead (SLead QRS )is calculated as the weighted sum SL QRS =P6 s=1 AL s.s P6 s=1 AL s (5.9) where AL sis the mean absolute amplitude of the QRS peaks at scale sof the DWT, and lead L AL s=1 D D X d=1 WL ss(ld), s = 1,2, . . . 6(5.10) being Dthe number of detected peaks (1 or 2) and ldthe positions of the peaks. The last morphology feature used is the maximum of the autocorrelation sequence of the ECG WT at scale 3 (rQRST(kM)), which describes the QRST complex similarity between PCA leads at scale 3 of the WT. This feature is related to changes in the multilead morphology and the depolarization axis of the QRST complex. See Figure 2.10 for details about the calculation of all the morphology features used. 5.2.6 Performance evaluation The performance is calculated from the confusion matrix after performing a classification experiment, in terms of the class sensitivity (Si), class positive predictive value (P+ i), global accuracy (A), global sensitivity (S) and global positive predictive value (P+) as suggested in [AAMI-EC57, 1998–2008] and described in [Llamedo and Martínez, 2011a, 2012a]. As the initialization of the EMC is random, the results of the clustering algorithm are not deterministic. Then each experiment is repeated 30 times to evaluate the mean and standard deviation of the performance estimates. The amount of expert assistance required in the patient-adaptable modes of operation will be also accounted for each experiment. 5.3 Results We performed two experiments, in the first one we studied the values of the algorithm parameters, that will be used in the second to evaluate its performance. The objective of the first experiment was to set up the number of clusters (K) and the qualified majority percentage used in votations (α), both parameters used in automatic and slightly 5.3. RESULTS 117 Table 5.3: Performance obtained in the development dataset for the election of Kand α parameters. Normal Supraventricular Ventricular Total Operation mode α K S P+S P+S P+A S P+ Slightly assisted 50 996 ±0 97 ±0 58 ±2 53 ±1 79 ±1 67 ±1 93 ±0 77 ±0 72 ±0 12 96 ±0 97 ±0 56 ±2 51 ±1 78 ±1 66 ±1 92 ±0 77 ±1 72 ±0 75 9 97±0 98±0 60±1 61±1 84±1 71±1 94±0 80±1 77±1 12 97 ±0 98 ±0 61 ±1 60 ±1 84 ±1 71 ±1 94 ±0 80 ±0 76 ±0 Automatic 0996 ±0 97 ±0 52 ±3 49 ±1 78 ±1 65 ±1 92 ±0 75 ±1 71 ±1 12 95 ±0 97 ±0 50 ±2 49 ±1 77 ±1 63 ±1 92 ±0 74 ±1 70 ±1 50 9 96±0 97±0 53±2 51±2 77±1 65±1 92±0 75±1 71±0 12 95 ±0 97 ±0 52 ±1 48 ±1 77 ±1 64 ±0 92 ±0 75 ±0 70 ±0 75 995 ±0 97 ±0 50 ±1 47 ±0 77 ±0 61 ±0 92 ±0 74 ±0 68 ±0 12 95 ±0 97 ±0 50 ±1 47 ±0 77 ±0 61 ±0 92 ±0 74 ±0 68 ±0 K: number of clusters; α: majority threshold (in percent) assisted modes of operation. These parameters were assessed in the development dataset (MITBIH-SUP and DS1 subset of MITBIH-AR), and then used for the final performance evaluation in the remaining datasets. Table 5.3 shows the results of this experiment for two values of the evaluated parameters. As a result of this experiment we adopted K= 9 and α= 50% for the automatic mode, and K= 9 and α= 75% for the slightly assisted mode. The final evaluation of the algorithm was performed in a broad set of databases in order to obtain a realistic estimation of its performance, as done in [Chudácek et al., 2009, Llamedo and Martínez, 2012a]. The three modes of operation were evaluated for each database with the parameter values obtained in the first experiment. The results of this experiment are presented in Tables 5.4 and 5.5 grouped by dataset. Comparison with the most relevant algorithms found in the literature are presented separately in Table 5.4. In Table 5.5 the performance obtained for all databases are presented. For each database we present the performance of our previous classifier [Llamedo and Martínez, 2011a] at the bottom for comparison, and a biased performance estimation on top as an upper bound. This biased performance is obtained when a quadratic classifier [Llamedo and Martínez, 2011a, van der Heijden et al., 2005, Duin et al., 2008] and the feature model presented in Table 5.2 are trained and tested in the same patient, for each patient in a database. This optimistically biased performance serves as an upper bound, and represents the performance of the model if it could be re-trained for each patient. From the results presented in Table 5.4, the proposed algorithm outperforms almost all reviewed algorithms, except the algorithms of Jiang [Jiang and Kong, 2007] and Ince [Ince et al., 2009], both in a small subset of MITBIH-AR, and the algorithm of Kiranyaz [Kiranyaz et al., 2011] in the MITBIH-LT. Finally the results showed in Table 5.5 evidence that the algorithm improves the baseline performances obtained by the LDC. 118 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.4: Performance comparison with reference algorithms. Normal Supraventricular Ventricular Total Dataset Observation #MAHB S P+S P+S P+A S P+ MITBIH-AR (DS2)• [de Chazal and Reilly, 2006] 500 94 99 88 47 95 82 94 92 76 [Llamedo and Martínez, 2012b] 12 100±0 99±0 92±1 90±3 92±1 97±1 99±0 94±1 96±1 [de Chazal et al., 2004] 0 87 99 76 39 87 43 87 83 60 [Mar et al., 2011] 0 90 99 83 34 87 61 89 87 65 [Llamedo and Martínez, 2012b] 0 93±0 99±0 77±0 39±0 82±0 70±0 92±0 84±0 69±0 MITBIH-AR (DS1-DS2)• [Ince et al., 2009] ≈300†98 98 64 54 85 87 96 82 80 [Llamedo and Martínez, 2012b] 12 100±0 99±0 89±2 88±3 90±1 97±0 98±0 93±1 95±1 MITBIH-AR (24 rec.)5 [Jiang and Kong, 2007] ≈300†99 96 51 68 85 93 95 78 86 [Llamedo and Martínez, 2012b] 12 99±0 98±0 91±1 88±2 90±1 96±1 98±0 93±1 94±1 MITBIH-AR (11 rec.)◦ [Hu et al., 1997] ≈300†– – – – 79 76 – – – [de Chazal and Reilly, 2006] 500 – – 76 39 76 91 – – – [Jiang and Kong, 2007] ≈300†– – 75 79 94 96 – – – [Ince et al., 2009] ≈300†– – 82 63 90 92 – – – [Llamedo and Martínez, 2012b] 12 99±0 99±0 89±2 88±3 93±1 97±1 98±0 85±1 85±2 MITBIH-LT [Kiranyaz et al., 2011] ≤900?99 100 40 17 98 99 99 79 72 [Llamedo and Martínez, 2012b] 20 99±1 99±0 16±16 0±0 94±0 91±4 98±0 70±5 74±9 ◦Comparison against results presented in Table II of [Ince et al., 2009]. 524 recordings from 200 to 234 used in [Jiang and Kong, 2007]. •DS1 and DS2 datasets defined in [de Chazal et al., 2004]. †Heartbeats in the first 5 min of each recording. ?Heartbeats in approx. 15 min of each recording. #MAHB: manually annotated heartbeats per recording. 5.A. DETAILED RESULTS 125 Table 5.7: Confusion matrices of the results presented in Table 5.5 for AHA database. Truth biased Algorithm n s v Total N 265807 – 1175 266982 S – – – – V 177 – 33833 34010 Total 265984 – 35008 300992 Truth Slightly assisted Algorithm n s v u Total N 305655±423 5548±361 6515±251 18±10 317735 S – – – – – V 1398±164 3399±284 29157±280 33±18 33987 U 171±21 7±5 118±25 276±27 572 Total 307224±478 8954±429 35789±365 327±44 352294 Truth Assisted 12 MAHB Algorithm n s v u Total N 316325±1054 – 642±112 31±20 317735 S – – – – – V 809±80 – 32885±334 61±40 33987 U 131±35 – 93±14 348±40 572 Total 317266±1068 – 33619±364 440±83 352294 Truth Automatic Algorithm n s v u Total N 296739±736 10589±363 10406±647 – 317735 S – – – – – V 1660±131 6653±231 25674±250 – 33987 U 201±23 65±11 306±24 – 572 Total 298600±752 17308±469 36387±818 – 352294 Truth Assisted 9 MAHB Algorithm n s v u Total N 316865±151 – 846±149 24±12 317735 S – – – – – V 975±140 – 32965±146 47±19 33987 U 145±26 – 114±22 314±28 572 Total 317984±267 – 33925±265 385±41 352294 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 289116 15746 12873 – 317735 S – – – – – V 1616 7003 25368 – 33987 U 253 77 242 – 572 Total 290985 22826 38483 – 352294 126 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.8: Confusion matrices of the results presented in Table 5.5 for ESTTDB database. Truth biased Algorithm n s v Total N 534551 1486 979 537016 S 12 1056 8 1076 V 96 49 4673 4818 Total 534659 2591 5660 542910 Truth Slightly assisted Algorithm n s v u Total N 759079±2090 20811±1122 4748±1184 0 784638 S 619±113 429±110 47±14 0 1095 V 585±61 359±111 3875±104 2±5 4821 U 8±3 2±2 1±1 1±1 11 Total 760291±2066 21601±1081 8671±1272 3±6 790565 Truth Assisted 12 MAHB Algorithm n s v u Total N 783911±394 177±102 532±300 18±27 784638 S 565±53 478±53 52±2 0 1095 V 418±127 18±9 4375±132 11±3 4821 U 10±1 0 0 1±1 11 Total 784904±509 673±149 4959±375 30±40 790565 Truth Automatic Algorithm n s v u Total N 724671±1956 52226±2271 7742±1147 – 784638 S 687±173 289±147 119±34 – 1095 V 608±95 493±107 3720±106 – 4821 U 7±2 1±1 1±1 – 11 Total 725973±1810 53011±2289 11581±1155 – 790565 Truth Assisted 9 MAHB Algorithm n s v u Total N 781204±7044 199±147 908±347 0 784638 S 632±87 413±94 50±11 0 1095 V 583±106 61±71 4175±120 2±5 4821 U 9±2 1±2 1±1 1±1 11 Total 782428±7048 673±295 5133±411 3±6 790565 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 700021 76527 8090 – 784638 S 433 517 145 – 1095 V 297 748 3776 – 4821 U 3 8 0 – 11 Total 700754 77800 12011 – 790565 5.A. DETAILED RESULTS 127 Table 5.9: Confusion matrices of the results presented in Table 5.5 for INCART database. Truth biased Algorithm n s v Total N 146861 886 458 148205 S 31 1917 10 1958 V 167 25 20030 20222 Total 147059 2828 20498 170385 Truth Slightly assisted Algorithm n s v u Total N 147498±904 5787±918 418±99 – 153703 S 313±85 1445±89 203±52 – 1960 V 512±139 944±85 18774±166 – 20230 U 3±1 0 3±1 – 6 Total 148325±925 8176±938 19398±281 – 175899 Truth Assisted 12 MAHB Algorithm n s v u Total N 153088±231 255±238 359±78 2±5 153703 S 266±56 1667±55 27±20 0 1960 V 405±85 20±18 19805±89 0 20230 U 3±1 0 3±1 0±1 6 Total 153761±278 1941±265 20195±144 2±6 175899 Truth Automatic Algorithm n s v u Total N 137391±936 15876±1052 436±241 – 153703 S 200±62 1464±105 296±85 – 1960 V 623±76 1806±206 17800±200 – 20230 U 3±1 1±1 3±1 – 6 Total 138217±931 19147±1082 18535±383 – 175899 Truth Assisted 9 MAHB Algorithm n s v u Total N 153046±169 184±133 472±114 2±5 153703 S 359±88 1550±91 51±50 0 1960 V 496±114 16±13 19718±118 0 20230 U 2±1 0 3±1 0±1 6 Total 153903±266 1750±144 20244±244 2±6 175899 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 135069 18013 621 – 153703 S 65 1496 399 – 1960 V 607 2561 17062 – 20230 U 3 0 3 – 6 Total 135744 22070 18085 – 175899 128 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.10: Confusion matrices of the results presented in Table 5.5 for LTSTDB database. Truth biased Algorithm n s v Total N 7651109 63997 12862 7727968 S 2071 63461 644 66176 V 437 894 72174 73505 Total 7653617 128352 85680 7867649 Truth Slightly assisted Algorithm n s v u Total N 8692959±1939 33376±8330 41710±6390 – 8768045 S 20152±448 19881±573 5588±125 – 45620 V 17462±480 2426±239 67070±719 – 86958 U 371±2 8±2 91±4 – 469 Total 8730943±2869 55691±8661 114458±5792 – 8901092 Truth Assisted 12 MAHB Algorithm n s v u Total N 8746770±123 13167±787 5698±910 – 8768045 S 18748±2845 24273±2032 2600±813 – 45620 V 14777±1051 6150±3415 66031±2364 – 86958 U 367±18 15±7 87±11 – 469 Total 8783073±1653 43605±603 74415±2256 – 8901092 Truth Automatic Algorithm n s v u Total N 8412810±54249 239484±9806 96661±5438 – 8768045 S 18640±968 22590±712 4371±253 – 45620 V 13328±2069 16229±1542 57388±527 – 86958 U 315±19 126±10 27±8 – 469 Total 8445092±54044 278429±8976 158446±5156 – 8901092 Truth Assisted 9 MAHB Algorithm n s v u Total N 8743249±6069 14798±3169 9999±2901 – 8768045 S 18225±301 23718±285 3678±585 – 45620 V 16059±1896 2643±28 68257±1924 – 86958 U 353±20 10±0 107±20 – 469 Total 8777885±8285 41168±2856 82039±5429 – 8901092 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 8048986 576535 142524 – 8768045 S 10563 25940 9117 – 45620 V 8663 17024 61271 – 86958 U 264 152 53 – 469 Total 8068476 619651 212965 – 8901092 5.A. DETAILED RESULTS 129 Table 5.11: Confusion matrices of the results presented in Table 5.5 for MITBIH-AR database. Truth biased Algorithm n s v Total N 67839 519 604 68962 S 38 2727 8 2773 V 123 32 7649 7804 Total 68000 3278 8261 79539 Truth Slightly assisted Algorithm n s v u Total N 88554±163 1286±184 287±91 0±1 90127 S 307±59 2225±144 248±118 0 2780 V 995±185 269±48 6547±188 0 7811 U 12±2 0 2±1 0±1 15 Total 89868±301 3780±202 7084±268 1±2 100733 Truth Assisted 12 MAHB Algorithm n s v u Total N 89681±121 226±101 219±70 1±1 90127 S 277±67 2477±70 26±12 0 2780 V 651±148 105±45 7055±136 0 7811 U 12±1 0 2±1 1±1 15 Total 90621±248 2807±185 7303±180 2±2 100733 Truth Automatic Algorithm n s v u Total N 86716±350 2375±256 1036±358 – 90127 S 303±72 2117±74 360±28 – 2780 V 1112±212 409±58 6290±212 – 7811 U 12±2 1±1 3±1 – 15 Total 88142±444 4901±233 7689±472 – 100733 Truth Assisted 9 MAHB Algorithm n s v u Total N 89555±122 318±103 254±75 0±1 90127 S 362±51 2358±52 34±13 0 2780 V 785±181 103±36 6923±189 0 7811 U 12±1 0 2±1 0±1 15 Total 90714±246 2806±122 7212±226 1±2 100733 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 85465 3220 1442 – 90127 S 284 2113 383 – 2780 V 1061 783 5967 – 7811 U 10 0 5 – 15 Total 86820 6116 7797 – 100733 130 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.12: Confusion matrices of the results presented in Table 5.5 for DS2 of MITBIH-AR database. Truth biased Algorithm n s v Total N 31871 288 168 32327 S 19 1812 1 1832 V 42 26 3536 3604 Total 31932 2126 3705 37763 Truth Slightly assisted Algorithm n s v u Total N 43007±217 1106±237 146±97 – 44259 S 74±28 1585±119 178±120 – 1837 V 324±132 77±31 3208±127 – 3609 U 6±1 0 1±1 – 7 Total 43411±291 2768±235 3534±222 – 49712 Truth Assisted 12 MAHB Algorithm n s v u Total N 44032±96 140±83 88±37 – 44259 S 133±33 1688±34 16±9 – 1837 V 246±71 66±41 3298±73 – 3609 U 6±1 0 1±1 – 7 Total 44416±131 1893±122 3402±84 – 49712 Truth Automatic Algorithm n s v u Total N 41902±335 1562±225 795±354 – 44259 S 7483±29 1476±21 277±26 – 1837 V 478±177 83±32 3048±173 – 3609 U 5±1 0 2±1 – 7 Total 42469±421 3121±218 4122±421 – 49712 Truth Assisted 9 MAHB Algorithm n s v u Total N 43942±109 208±91 110±57 – 44259 S 134±35 1684±36 20±10 – 1837 V 289±110 69±33 3251±105 – 3609 U 6±1 0 1±1 – 7 Total 44369±190 1961±101 3382±148 – 49712 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 41209 2045 1005 – 44259 S 141 1413 283 – 1837 V 434 211 2964 – 3609 U 4 0 3 – 7 Total 41788 3699 4255 – 49712 5.A. DETAILED RESULTS 131 Table 5.13: Confusion matrices of the results presented in Table 5.5 for MITBIH-LT database. Truth biased Algorithm n s v Total N 579074 9217 6381 594672 S 56 1431 12 1499 V 2489 2872 61008 66369 Total 581619 13520 67401 662540 Truth Slightly assisted Algorithm n s v u Total N 570716±4038 – 29516±4038 – 600232 S 1287±31 – 213±31 – 1500 V 11197±991 – 55806±991 – 67003 U – – – – – Total 583200±4998 – 85535±4998 – 668735 Truth Assisted 20 MAHB Algorithm n s v u Total N 593776±3263 148±149 6309±3114 – 600232 S 1029±291 242±244 230±48 – 1500 V 4118±303 0 62885±303 – 67003 U – – – – – Total 598923±3252 390±393 69423±2859 – 668735 Truth Automatic Algorithm n s v u Total N 529319±1917 39867±2631 31047±714 – 600232 S 1313±36 24±3 163±33 – 1500 V 26086±877 16092±252 24826±1129 – 67003 U – – – – – Total 556718±1076 55983±2886 56035±1810 – 668735 Truth Assisted 9 MAHB Algorithm n s v u Total N 594826±1502 121±122 5286±1624 – 600232 S 1288±30 11±11 202±19 – 1500 V 5055±1711 2061±2078 59888±3789 – 67003 U – – – – – Total 601168±3183 2192±2210 65376±5394 – 668735 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 513227 55080 31925 – 600232 S 873 493 134 – 1500 V 26256 15288 25459 – 67003 U – – – – – Total 540356 70861 57518 – 668735 132 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.14: Confusion matrices of the results presented in Table 5.5 for MITBIH-ST database. Truth biased Algorithm n s v Total N 13772 100 5 13877 S 17 779 1 797 V 0 0 314 314 Total 13789 879 320 14988 Truth Slightly assisted Algorithm n s v u Total N 33065±794 12323±690 839±171 – 46228 S 194±68 604±68 0±1 – 798 V 13±7 25±362 281±34 – 319 U – – – – – Total 33272±817 12952±706 1121±194 – 47345 Truth Assisted 12 MAHB Algorithm n s v u Total N 46118±46 107±46 3±2 – 46228 S 96±21 702±21 0±1 – 798 V 13±6 4±2 302±6 – 319 U – – – – – Total 46227±54 813±52 306±7 – 47345 Truth Automatic Algorithm n s v u Total N 30781±341 14351±366 1096±229 – 46228 S 316±72 470±67 12±11 – 798 V 9±5 164±30 146±30 – 319 U – – – – – Total 31106±320 14984±342 1254±222 – 47345 Truth Assisted 9 MAHB Algorithm n s v u Total N 46076±75 149±75 3±1 – 46228 S 128±44 669±44 0±1 – 798 V 16±8 4±2 299±8 – 319 U – – – – – Total 46220±103 822±101 303±8 – 47345 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 30168 14339 1721 – 46228 S 266 491 41 – 798 V 1 149 169 – 319 U – – – – – Total 30435 14979 1931 – 47345 5.A. DETAILED RESULTS 133 Table 5.15: Confusion matrices of the results presented in Table 5.5 for MITBIH-SUP database. Truth biased Algorithm n s v Total N 156674 4810 787 162271 S 359 11558 278 12195 V 105 515 9343 9963 Total 157138 16883 10408 184429 Truth Slightly assisted Algorithm n s v u Total N 155136±485 4000±452 3202±165 2±6 162340 S 3311±268 7530±228 1355±155 2±4 12198 V 595±170 930±169 8439±209 3±5 9966 U 23±7 27±5 28±6 2±2 79 Total 159065±785 12486±663 13024±450 8±15 184583 Truth Assisted 12 MAHB Algorithm n s v u Total N 159991±300 1869±327 455±132 25±37 162340 S 2597±223 9081±206 507±106 12±20 12198 V 606±108 600±103 8754±112 7±11 9966 U 33±7 11±3 26±5 6±7 79 Total 163230±458 11561±459 9742±253 50±56 184583 Truth Automatic Algorithm n s v u Total N 153420±415 4496±301 4424±312 – 162340 S 3993±266 5751±420 2454±304 – 12198 V 518±71 1305±55 8143±76 – 9966 U 20±5 27±3 32±4 – 79 Total 157951±513 11579±550 15053±326 – 184583 Truth Assisted 9 MAHB Algorithm n s v u Total N 159842±510 2008±517 468±132 22±32 162340 S 2820±321 8871±313 498±113 9±15 12198 V 687±213 620±162 8656±231 4±6 9966 U 38±7 10±4 26±7 5±5 79 Total 163386±876 11509±853 9648±439 40±53 184583 Truth Automatic [Llamedo and Martínez, 2012a] Algorithm n s v u Total N 151246 6040 5054 – 162340 S 2735 6185 3278 – 12198 V 704 1365 7897 – 9966 U 12 37 30 – 79 Total 154697 13627 16259 – 184583 134 CHAPTER 5. PATIENT-ADAPTED ECG HEARTBEAT CLASSIFICATION Table 5.16: Confusion matrices of the results presented in Table 5.6 for MITBIH-AR database. Truth Assisted 12 MAHB Algorithm n s v f u Total N 89681±121 226±101 174±67 45±30 1±1 90127 S 277±67 2477±70 26±12 0 0 2780 V 417±133 96±42 6373±148 121±71 0 7008 F 234±48 9±5 159±100 402±116 0 803 U 12±1 0 2±1 1±1 1±1 15 Total 90621±248 2807±185 6734±215 569±184 2±2 100733 Truth Assisted 9 MAHB Algorithm n s v f u Total N 89555±122 318±103 202±68 51±41 0±1 90127 S 362±51 2358±52 33±13 0 0 2780 V 540±157 93±34 6178±185 197±105 0 7008 F 245±71 11±5 146±110 402±86 0 803 U 12±1 0 1±1 1±1 0±1 15 Total 90714±246 2806±122 6561±250 651±190 1±2 100733 Truth Slightly assisted Algorithm n s v f u Total N 88602±206 1251±206 202±68 44±35 28±39 90127 S 311±57 2222±138 246±112 0 1±2 2780 V 707±158 111±41 5995±189 195±106 1±1 7008 F 290±82 11±5 145±110 358±130 0 803 U 12±1 0 2±1 1±1 0±1 15 Total 89922±338 3594±225 6589±324 597±236 30±41 100733 Truth Automatic Algorithm n s v f u Total N 86958±277 2359±248 344±119 374±160 91±59 90127 S 273±65 2056±121 449±99 0±1 2±6 2780 V 652±111 398±60 5868±108 86±56 4±2 7008 F 511±49 6±3 95±60 191±63 0 803 U 10±2 1±1 3±1 0 1±1 15 Total 88405±275 4819±268 6760±207 651±156 95±59 100733