Estrategias estáticas y dinámicas para el agrupamiento de latidos mediante acumulación de evidencia

TESIS DE DOCTORADO ESTRATEGIAS ESTÁTICAS Y DINÁMICAS PARA EL AGRUPAMIENTO DE LATIDOS MEDIANTE ACUMULACIÓN DE EVIDENCIA David González Márquez ESCUELA DE DOCTORADO INTERNACIONAL PROGRAMA DE DOCTORADO EN INVESTIGACIÓN EN TECNOLOXÍAS DA INFORMACIÓN SANTIAGO DE COMPOSTELA AÑO 2017 Life is like an ECG diagram. If it goes smoothly, you are dead. Internet Mi corazón, mi corazón es un músculo sano pero necesita acción. Dame paz y dame guerra, y un dulce colocón y yo te entregaré lo mejor. Calamaro Agradecimientos Me gustaría comenzar agradeciendo a mis directores Paulo Félix y Abraham Otero todo su trabajo, apoyo y dedicación. Han sido ellos los que me han acompañado, guiado y arrastrado (cuando hacía falta) y sin ellos nada de esto sería posible. Agradecer especialmente su paciencia y trato conmigo, sé que en este camino los tropiezos han sido míos y siempre habéis estado ahí para intentar apagar los fuegos. Gracias a Paulo por acogerme en su casa (el CiTIUS) y abrirme las puertas de par en par, poniendo todos los medios necesarios para que esta tesis haya podido seguir adelante. Gracias a Abraham por introducirme en el mundo de la investigación y llevarme de la mano todos estos años (y los que nos quedan). Todo agradecimiento queda corto. Debo dar también las gracias a la profesora Ana Fred, por acogerme hasta en dos ocasiones en Lisboa y por sus consejos y dirección, sin ella esta tesis tampoco habría sido posible. Igualmente agradecer al Instituto Superior Técnico de Lisboa y al grupo “Pattern and Image Analysis – Lx” del Instituto de Telecomunicações que me recibieron durante mis dos estancias en Lisboa. También debo tener en cuenta a todos los compañeros y amigos que encontré en Lisboa, que hicieron de esos meses de estancia una experiencia para recordar (tanto que tuve que volver). Dar las gracias al CiTIUS, que ha sido mi casa durante estos años y me ha permitido formarme y desarrollar mi tesis en un entorno óptimo. Entorno que no sería posible sin todas las personas que lo forman, desde los conserjes hasta el servicio técnico (especial mención a Jorge, todo un santo). También merecen una mención aparte los compañeros de laboratorio que me acogieron como uno más (a pesar de que el gallego no sea lo mío) y sin los cuales probablemente estaría interno en algún psiquiátrico. Agradecimiento extraordinario se merece Tino que por si no fuera poco con andar por este camino que compartimos tiene que ir arreglando los desastres que le dejo a mi paso. XII Índice general 2.2.3. Selección de características . . . . . . . . . . . . . . . . . . . . . . . 39 2.3. Coste computacional de la representación de Hermite . . . . . . . . . . . . . 41 2.4. Discusión .................................... 42 3 Cálculo de la representación de Hermite utilizando GPUs 47 3.1. Programación utilizando GPUs . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2. Implementación optimizada en C . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3. Implementación paralela . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1. Precomputación de las funciones de Hermite . . . . . . . . . . . . . 53 3.3.2. Caracterización del complejo QRS mediante Hermite . . . . . . . . . 53 3.3.3. Optimización de la transferencia de datos . . . . . . . . . . . . . . . 54 3.4. ResultadosyDiscusión............................. 56 3.4.1. Test A: Procesado en diferido (offline) de registros cortos . . . . . . . 57 3.4.2. Test B: Procesado en diferido (offline) de registros largos . . . . . . . 63 3.4.3. Test C: Procesado en tiempo real (online) . . . . . . . . . . . . . . . 65 4 Agrupamiento de latidos mediante acumulación de evidencia 69 4.1. Agrupamiento mediante acumulación de evidencia positiva y negativa (PNEAC) ...................................... 72 4.1.1. Generación de las particiones del ensemble .............. 73 4.1.2. Combinación de las particiones de datos . . . . . . . . . . . . . . . . 74 4.1.3. Extracción de la partición final de los datos . . . . . . . . . . . . . . 76 4.1.4. Análisis de complejidad computacional de PN-EAC . . . . . . . . . 78 4.2. Agrupamiento de latidos con PN-EAC . . . . . . . . . . . . . . . . . . . . . 79 4.2.1. Estrategias para la generación de particiones . . . . . . . . . . . . . 81 4.2.2. Resultados ............................... 86 4.2.3. Discusión................................ 89 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC . . . . . . 93 4.3.1. Estrategias para la generación de particiones . . . . . . . . . . . . . 95 4.3.2. Resultados ............................... 97 4.3.3. Discusión................................ 100 5 Agrupamiento dinámico de latidos mediante acumulación de evidencia 105 Índice general XIII 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa(EPN-EAC) .............................. 106 5.1.1. Descripción de la técnica EPN-EAC . . . . . . . . . . . . . . . . . . 107 5.1.2. Inicialización del algoritmo . . . . . . . . . . . . . . . . . . . . . . 109 5.1.3. Búsqueda del par de objetos más similares en W........... 110 5.1.4. Fusión de dos objetos e inclusión del nuevo objeto xz......... 112 5.1.5. Recolección de evidencia . . . . . . . . . . . . . . . . . . . . . . . . 113 5.1.6. Extracción de la partición final . . . . . . . . . . . . . . . . . . . . . 114 5.1.7. Ejemplo ilustrativo del funcionamiento del algoritmo . . . . . . . . . 115 5.1.8. Análisis de complejidad computacional . . . . . . . . . . . . . . . . 121 5.2. Aplicación al agrupamiento de latidos . . . . . . . . . . . . . . . . . . . . . 122 5.3. Resultados.................................... 125 5.4. Discusión .................................... 128 Conclusiones y trabajo futuro 135 Conclusions and future work 141 Bibliografía 147 Índice de figuras 161 Índice de tablas 167 CAPÍTULO 1 INTRODUCCIÓN El campo de la electrocardiografía tiene su origen en los trabajos de Galvani [48] que observó las contracciones de unas ancas de rana al estimularlas con electricidad. Esta mal llamada en su momento “electricidad animal” marcó el inicio de un periodo de investigación sobre la interacción entre la electricidad y los sistemas biológicos presentes en el cuerpo humano. Serían posteriormente Matteucci, Kolliker y Muller, entre otros, a los que correspondería iniciar este nuevo campo y demostrar la relación entre el ciclo cardíaco y la actividad eléctrica del miocardio [98]. Este desarrollo teórico fue acompañado de avances en los instrumentos de medida de la mano de Waller [142] y Einthoven [118] (véase Figura 1.1). Desde que por primera vez Pipberger digitalizó el electrocardiograma (ECG) en 1959 y diseñó los primeros programas para su análisis [114], el interés en esta señal fisiológica no ha parado de crecer, siendo en la actualidad una prueba fundamental en la rutina clínica. Su bajo coste y su sencillez lo convierten en una herramienta excepcional para el diagnóstico y seguimiento de enfermedades cardiacas, particularmente en países en vías de desarrollo. Además, su carácter no invasivo y la simplicidad de la instrumentación necesaria para su adquisición, hacen de él un candidato ideal para realizar una monitorización de larga duración durante la actividad rutinaria del paciente. Este interés en el ECG se ve reforzado por la alta mortalidad de las enfermedades cardiovasculares, que actualmente son la primera causa de muerte por enfermedad en el mundo [106]. Lejos de resolverse esta circunstancia, en el futuro cercano se prevé que continuará incrementando (véase Figura 1.2), especialmente en los países en vías de desarrollo, por los cambios en la dieta y estilo de vida derivados del mayor poder adquisitivo de sus ciudadanos [95]. 2Capítulo 1. Introducción Powered by TCPDF (www.tcpdf.org) Figura 1.1: Antiguo electrocardiógrafo comercial construido en 1911 por la Cambridge Scientific Instrument Company. (Fuente: [24]) Cardiopatía isquémica Accidente cerebrovascu... Infecciones vías respirato... EPOC Diarrea Tuberculosis VIH/sida Prematuridad Cánceres tráquea, bron... Asfixia al nacer 0246810 Millones de defunciones Causa de defunción Las 10 principales causas de defunción en el mundo en el año 2000 Cardiopatía isquémica Accidente cerebrovascu... Infecciones vías respirato... EPOC Cánceres tráquea, bron... Diabetes mellitus Alzheimer y otras deme... Diarrea Tuberculosis Accidentes de tránsito 0246810 Millones de defunciones Causa de defunción Las 10 principales causas de defunción en el mundo en 2015 Figura 1.2: Principales causas de defunción en el mundo en 2000 y 2015. (Fuente: [106]) 3 El electrocardiograma, además de servir como instrumento al servicio del diagnóstico de la patología cardiaca, es una fuente inagotable de información [92]. Gracias a su estudio es posible comprender mejor la compleja interacción entre los distintos procesos fisiopatológicos que concurren en las alteraciones del impulso eléctrico en el miocardio. Su análisis es objeto de innumerables trabajos científicos y la aparición de nuevas aplicaciones es continua: la estimación de la salud fetal en obstetricia [90], el seguimiento de pacientes crónicos como en el caso de la diabetes [25], la enfermedad pulmonar obstructiva crónica [64], la apneahipopnea del sueño [110], o incluso el diseño de nuevos fármacos [144], son solo algunos ejemplos. En este contexto es indudable el interés en desarrollar nuevos procedimientos y métodos para el análisis computacional del ECG. Cualquier mejora en este sentido redundará en un gran beneficio por su aplicación en un amplio conjunto de escenarios que hacen uso de esta prueba médica. Estas mejoras podrían permitir una monitorización más precisa, ampliar el conocimiento sobre el miocardio y su patología y en última instancia podrían resultar en una mejora del tratamiento y la consecuente reducción en la pérdida de vidas humanas y mejora de calidad de vida. Sin embargo, parece que en los últimos años los avances en el campo de la electrocardiografía se han ralentizado [4]. La instrumentación de adquisición, excepto pequeñas mejoras incrementales, no ha experimentado avances significativos a la hora de permitir obtener una mayor información acerca del funcionamiento electrofisiológico del miocardio. Por ello, es de esperar que los futuros avances no surgirán de un cambio en la tecnología de adquisición, sino de mejoras en el procesamiento del ECG, concretamente del procesado computacional. Son de especial importancia en el análisis computacional del ECG las tareas que se sitúan en la primera fase de detección e identificación del latido cardíaco. A partir de estas tareas se caracterizará el ciclo cardiaco y se construirá el resto del análisis. Errores en esta primera fase pueden invalidar el resto de la interpretación realizada. Por ello estas tareas requieren una atención especial por parte del cardiólogo y en muchas ocasiones siguen siendo realizadas mediante inspección visual. Los registros Holter, habitualmente utilizados para el diagnóstico, pueden durar varios días; por tanto, su inspección visual resulta tediosa y el resultado de esta depende excesivamente del cardiólogo que la realiza. Además, debemos tener en cuenta la gran cantidad de pruebas de electrocardiograma que se realizan cada año (aproximadamente 200 millones en todo el mundo) [117]. 4Capítulo 1. Introducción En la siguiente sección realizaremos una introducción al funcionamiento electrofisiológico del corazón que servirá de base para comprender los problemas que serán abordados en la presente tesis doctoral. 1.1. Introducción a la electrocardiografía El electrocardiograma es la representación gráfica de la actividad eléctrica del corazón en función del tiempo. Para su estudio es necesario tener una noción básica del funcionamiento del corazón, prestando especial atención al aspecto eléctrico. En esto consistirá la primera parte de esta introducción al campo de la electrocardiografía. 1.1.1. Funcionamiento eléctrico del corazón El corazón es un músculo hueco que mantiene la circulación sanguínea por medio de contracciones y dilataciones rítmicas. Está compuesto por cuatro cavidades, dos ventrículos y dos aurículas, situadas en las mitades superior e inferior del corazón, respectivamente. En la Figura 1.3 se muestran dichas cavidades junto con el nódulo sinoauricular (SA), el nódulo auriculoventricular (AV) y el Haz de His (las fibras de Purkinje no se muestran en la figura por claridad, pero aparecerían a continuación de las ramas del Haz de His alrededor de los ventrículos). Como veremos a continuación, estas estructuras juegan un papel fundamental a la hora de comprender la propagación del estímulo eléctrico. La aurícula derecha recibe la sangre venosa del cuerpo y la envía al ventrículo derecho, el cual la envía a su vez a los pulmones, lugar en el que la sangre se oxigena. Simultáneamente la aurícula izquierda recibe sangre oxigenada de los pulmones, que envía al ventrículo izquierdo para que la distribuya por todo el cuerpo. Para que este ciclo funcione de forma adecuada las contracciones se deben realizar de forma cíclica y ordenada. Para ello existe un sistema de estimulación y conducción eléctrica compuesto por fibras del músculo cardíaco especializadas en la generación y transmisión de impulsos eléctricos. El correcto funcionamiento de este sistema da lugar a un registro normal de la actividad eléctrica del corazón, como el que se observa en la Figura 1.4). La acción de bombeo del corazón proviene de un sistema intrínseco de conducción eléctrica. El nódulo sinoauricular (también conocido como el “marcapasos del corazón”) está compuesto por células especializadas cuya velocidad de despolarización espontánea es superior al resto de las células cardíacas. El impulso eléctrico normalmente es generado en 1.1. Introducción a la electrocardiografía 5 Powered by TCPDF (www.tcpdf.org) Figura 1.3: Representación del corazón con sus principales componentes. (Fuente: [19]) Figura 1.4: Fragmento de un ECG con latidos normales. (Fuente: MIT-BIH Arrhythmia Database, registro 100, entre 0:0:0 y 0:0:10) este nódulo y propagado por el miocardio, causando la contracción ordenada de cada una de sus cavidades. Este nódulo genera un impulso eléctrico de 60 a 100 veces por minuto en condiciones normales. Dicho impulso se propagará por las aurículas provocando su contracción, y alcanzará el nódulo auriculoventricular. La velocidad de conducción eléctrica de las células de este nódulo es inferior a la asociada al tejido muscular en aurículas y ventrículos, lo que constituye una barrera en situaciones normales. Esto evita que los ventrículos desarrollen contracciones rápidas, desincronizadas e ineficientes como ocurre, por ejemplo, en la fibrilación ventricular (véase Figura 1.5). Cuando el impulso eléctrico alcanza este nódulo se produce un pequeño retraso (aproximadamente 0.1 segundos). Si el impulso se detiene en este punto se podría desembocar en una parada cardíaca, a menos que algunas células se despolaricen independientemente de dicho nódulo, originando una sucesión de latidos denominados “latidos de escape” (véase Figura 1.6). En los ventrículos existen haces de células en torno a las cuales la velocidad de conducción del impulso eléctrico es muy superior a la asociada a las paredes musculares, lo que permite la contracción simultánea de toda la cavidad, fundamental para un bombeo adecuado. Estos 6Capítulo 1. Introducción Figura 1.5: Fragmento de un ECG en el que se aprecia fibrilación ventricular. (Fuente: MITBIH Arrhythmia Database, registro 207, entre 0:0:40 y 0:0:50) Figura 1.6: Fragmento de un ECG en el que aparecen una serie de latidos normales hasta que, aproximadamente en la mitad del fragmento, las distancias entre latidos se incrementan y aparecen latidos de escape. (Fuente: MIT-BIH Arrhythmia Database, registro 222, entre 0:12:36 y 0:12:46) Figura 1.7: Fragmento de un ECG en el que aparecen latidos con bloqueo de rama derecha. (Fuente: MIT-BIH Arrhythmia Database, registro 212, entre 0:0:0 y 0:0:10) haces integran el haz de His, que desciende desde el nódulo AV a lo largo de los ventrículos dividiéndose en dos ramas, derecha e izquierda. El bloqueo de estos haces origina situaciones de bloqueo de rama derecha o izquierda (véase Figura 1.7). A causa de la alta velocidad de conducción, en condiciones normales, el impulso se propaga rápidamente por el haz de His hasta alcanzar las fibras de Purkinje, generando la contracción ventricular. Este proceso se corresponde con la etapa de sístole que en situaciones de normalidad tiene una duración inferior a 100 ms. Tras la despolarización se produce la repolarización ventricular durante la cual las fibras musculares de los ventrículos se preparan para el próximo latido. Dado que los procesos de despolarización y repolarización se producen de forma coordinada en miles de células, las corrientes resultantes son suficientemente grandes como para generar diferencias de potencial de varios mV entre los distintos puntos de la superficie del cuerpo humano, las cuales pueden medirse con electrodos superficiales para obtener la señal de ECG. 1.1. Introducción a la electrocardiografía 7 P Q R S T Intervalo QT Intervalo PR Complejo QRS Segmento ST Segmento PR Figura 1.8: Imagen de un latido con las principales ondas y segmentos señalados. (Adaptado de [128]) 1.1.2. Trazo eléctrico del latido cardíaco Los elementos principales del trazado de un latido en el ECG son la onda P, el complejo QRS y la onda T (véase Figura 1.8). Estas ondas y los segmentos que las conectan tienen una relación directa con los eventos eléctricos que se producen en el corazón. La onda P corresponde a la despolarización auricular, como se aprecia en la Figura 1.9. La repolarización de las aurículas quedará eclipsada por la despolarización ventricular (complejo QRS) por lo que en un ECG normal no será posible visualizarla. El complejo QRS corresponde a la despolarización ventricular, que provoca la contracción de los ventrículos derecho e izquierdo (véase Figura 1.10). Finalmente, la onda T representa la repolarización de los ventrículos que recuperan su estado de reposo hasta la siguiente contracción (véase Figura 1.11). El intervalo QT mide la distancia entre el inicio del complejo QRS y el final de la onda T y se corresponde al tiempo entre la despolarización y la repolarización ventricular. Un complejo QRS normal tiene una duración de entre aproximadamente 60 y 120 ms [71] y generalmente es la parte central y más visible del ECG. La masa cardíaca de los ventrículos es mucho mayor que la de las aurículas, ya que deben bombear la sangre hacia los pulmones 14 Capítulo 1. Introducción de cara a la interpretación del ECG. La interpretación requiere una gran cantidad de tiempo, especialmente en registros de larga duración, necesarios en aquellos casos en que los síntomas aparecen intermitentemente [73]. A mayor duración del registro, mayor el tiempo necesario para la interpretación y mayor la necesidad de soporte computacional. Aunque la detección de arritmias es uno de los primeros pasos en el análisis del ECG, aportando información fundamental para el diagnóstico de una multitud de enfermedades cardíacas, en la práctica totalidad de las propuestas de la bibliografía se observa un comportamiento excelente sobre las bases de datos de referencia y un comportamiento mediocre cuando han de abordar el análisis de registros correspondientes a nuevos pacientes [63]. Estamos, por tanto, ante un problema abierto, y aún lejos de proporcionar una respuesta lo suficientemente satisfactoria. Entre las dificultades que se identifican en la resolución de este problema se encuentran: 1) la variabilidad de los procesos fisiológicos entre distintos pacientes, y aún en el mismo paciente; 2) el carácter estocástico de estos procesos; 3) la posibilidad de que múltiples procesos fisiológicos concurran simultáneamente, y sus múltiples interacciones; 4) la presencia de ruido y artefactos que enmascaran los procesos fisiopatológicos de interés; 5) la ausencia de modelos de funcionamiento del miocardio suficientemente completos y fiables, dejando en manos de la experiencia del cardiólogo la resolución del problema; y 6) el conocimiento tácito, subjetivo y difícilmente formalizable que constituye la experiencia del cardiólogo. Dado que las arritmias se deben a cambios en el punto o el instante de la activación del impulso eléctrico y/o alteraciones en el camino de propagación, estas se reflejan visualmente en el ECG, afectando a la morfología del latido (véanse Figuras 1.5 y 1.7) o a la separación entre latidos (véase Figura 1.6). Es por ello que la identificación automática mediante técnicas computacionales de las distintas morfologías de latido presentes en un registro sería de gran ayuda al cardiólogo como paso previo a su interpretación. En la bibliografía aparecen distintas técnicas computacionales para esta labor, pudiendo distinguirse dos grandes tipos: técnicas de clasificación y técnicas de agrupamiento [21]. 1.2.1. Clasificación de latidos En la bibliografía la tarea de clasificación de latidos consiste generalmente en asignar a cada latido una etiqueta identificando el origen del latido. Idealmente dicha etiqueta sería la misma que asignaría un cardiólogo al latido. En [37] la Association for the Advancement of Medical Instrumentation (AAMI), entre otras recomendaciones para la evaluación del 1.2. Identificación morfológica de latidos 15 rendimiento de algoritmos de análisis del ritmo cardíaco, propone una división de los latidos en 5 tipos: latido normal (N), supraventricular (S), ventricular (V), fusión (F) e indeterminado (Q). Esta clasificación ha sido aceptada y utilizada por la mayor parte de los algoritmos de clasificación, convirtiéndose en un estándar de facto. La mayoría de los clasificadores requieren una etapa de entrenamiento para la cual es necesario contar con un conjunto de datos suficientemente representativo y etiquetado, lo que en el caso que nos atañe supondría contar con un conjunto amplio de latidos ya preinterpretados por el cardiólogo. En la bibliografía aparece una multitud de técnicas de aprendizaje automático aplicadas a esta tarea [87], como el Análisis Discriminante Lineal de Fisher [30], los Modelos Ocultos de Markov [8], Filtros Bayesianos [122] o las Máquinas de Vectores Soporte [107]; y técnicas que proceden de la inteligencia artificial, como las Redes de Neuronas Artificiales [40], los Algoritmos Evolutivos [54], o la Teoría de la Resonancia Adaptativa [23]. El principal punto débil de las técnicas de clasificación es la fuerte dependencia del resultado con el conjunto de entrenamiento y la diversidad presente en dicho conjunto. Las diferencias entre pacientes hacen que no se pueda asumir que un clasificador entrenado en un conjunto de datos producirá un resultado válido en nuevos pacientes [85] (véase Figura 1.13). Por ello, en ocasiones se opta por entrenar un clasificador sobre una amplia base de datos de latidos obtenidos de múltiples pacientes, incorporando a posteriori otra etapa de entrenamiento específica para cada paciente sobre un conjunto de latidos anotados. Las propuestas de la bibliografía con mejores resultados siguen este método, acercándose más a la idea de un clasificador asistido que a un clasificador completamente automático [31][63][72]. Incluso en este escenario, las diferencias observadas a lo largo del tiempo en la morfología de los latidos de un mismo paciente hacen que no se pueda asumir que la morfología de un tipo de latido no sufrirá alteraciones (véase Figura 1.14), pudiendo incluso aparecer nuevos tipos de latidos, no presentes en el conjunto de entrenamiento [67]. Otra dificultad para los clasificadores es que latidos con diferencias significativas a nivel morfológico pueden corresponderse con un mismo tipo de latido para el cardiólogo (véase Figura 1.13) [143]. Esto complica el diseño de un clasificador que use las mismas clases que el cardiólogo. Por otro lado, para mantener el número de dimensiones del espacio de características usado en la representación del latido en unos márgenes razonables, la gran mayoría de los trabajos de la literatura emplean una o dos derivaciones como fuente de dichas características [30][68][77]. Si en la fase de entrenamiento del clasificador se usaron 16 Capítulo 1. Introducción latidos provenientes siempre de las mismas derivaciones, la generalidad del clasificador al aplicarse sobre derivaciones diferentes no está garantizada, por lo que su uso se ve restringido a aquellos escenarios donde se dispone de las mismas derivaciones sobre las cuales se realizó el entrenamiento. Por contra, usar en el entrenamiento latidos provenientes de derivaciones diferentes resulta un reto significativo por la gran variabilidad morfológica de un mismo latido registrado en diferentes derivaciones (véase Figura 1.17). 1.2.2. Agrupamiento de latidos La tarea del agrupamiento consiste en particionar un conjunto de datos en grupos, de tal modo que dicha partición refleje la estructura subyacente a los datos. En este caso no es necesario un conjunto de entrenamiento con latidos etiquetados, ni ningún tipo de entrenamiento específico para cada paciente. Al no existir una etapa previa de entrenamiento, tampoco existe dependencia con ninguna derivación concreta a la hora de realizar el análisis del ECG. Un agrupamiento exitoso permite resumir de forma eficaz registros de ECG, presentando al cardiólogo una serie de grupos donde los latidos pertenecientes a un mismo grupo presentan características similares; esto es, un mismo origen y camino de propagación del impulso eléctrico por el miocardio. El agrupamiento permite analizar las distintas familias morfológicas, algunas de ellas claramente asociadas a comportamientos patológicos, establecer las familias de normalidad, e identificar cuándo se producen cambios morfológicos, para luego proceder, mediante inspección del registro, a analizar su contexto temporal. En este escenario el hecho de que latidos con morfologías diferentes se correspondan con un mismo tipo de latido no supone un reto, ya que el algoritmo puede encontrar grupos diferentes para cada morfología del mismo tipo de latido, y el cardiólogo puede establecer la correspondencia de varios grupos con un mismo tipo. El uso de técnicas de agrupamiento para la identificación morfológica de latidos resuelve buena parte de las desventajas del uso de clasificadores, a costa de requerir la intervención del cardiólogo para realizar la interpretación final de los grupos. Mientras que con un algoritmo de clasificación es teóricamente posible proporcionar directamente una interpretación de cada latido cardiaco (implícita en la etiqueta asignada), en el agrupamiento será la interacción con el cardiólogo la que permita llegar a dicha interpretación [93]. Si tenemos en cuenta que en la actualidad en la rutina clínica siempre se lleva a cabo un proceso de inspección visual del ECG, el escenario en el cual en vez de analizar directamente el ECG simplemente se valoran algunos representantes de un número reducido de grupos de latidos constituye una mejora 1.2. Identificación morfológica de latidos 17 notable. Además, desde nuestro punto de vista, una aproximación donde la decisión final está en manos del cardiólogo, como es el caso del agrupamiento, es más susceptible de ser incorporada a la rutina clínica que una solución que trate de automatizar todo el proceso. Es por ello que en esta tesis doctoral se hace una apuesta por las técnicas de agrupamiento como una herramienta de apoyo al cardiólogo en la identificación morfológica de latidos. En el agrupamiento de latidos es habitual reducir la representación del latido a su rasgo más significativo: el complejo QRS. El motivo principal es la difícil identificación y caracterización de todos los constituyentes del latido: en particular, la forma y tamaño de las ondas P y T hacen que en ocasiones aparezcan enmascaradas por el ruido de la señal de ECG. Incluir estas ondas podría hacer que en el agrupamiento se encontraran diferencias derivadas de inexactitudes en la representación y no de distintas morfologías. Limitando la representación del latido al complejo QRS se gana en sencillez y en uniformidad en la representación, además de reducir el tamaño del espacio de características. En la bibliografía podemos encontrar múltiples trabajos que realizan agrupamiento de latidos, de los que aquí revisamos los más relevantes. El trabajo más referenciado en la bibliografía es [77]. En este trabajo se realiza un agrupamiento basado en la morfología del complejo QRS y en la distancia entre latidos. Se utilizó la base de datos MIT-BIH Arrhythmia Database al completo, que fue remuestreada a 1000 Hz. La detección de latidos se realizó mediante un algoritmo propio que detecta la posición de la onda R de cada latido. Posteriormente, se extrajo una ventana fija de 200 ms de señal alrededor de la posición dada por el anotador para representar cada complejo QRS. El agrupamiento posterior se realiza por medio de mapas auto-organizados (SOM). Los mapas auto-organizados están compuestos por varias neuronas que evolucionan mediante aprendizaje por competición, utilizando una función de vecindad para preservar las propiedades topológicas del espacio de representación. En [77] se utilizaron 25 neuronas, dando lugar a 25 grupos por registro. Los latidos fueron divididos en 16 tipos diferentes, excluyendo los latidos que solo contienen onda P, obteniendo en el mejor de los casos un error total cercano al 1.5%. Otra técnica orientada a reducir el número de latidos que tiene que revisar el cardiólogo se presenta en [27]. En este caso la técnica consta de dos etapas: en la primera se intenta reducir el número de latidos a agrupar calculando la similitud entre latidos, en la segunda se realiza el agrupamiento propiamente dicho, con dos alternativas: el algoritmo K-means y el algoritmo Max-Min. La técnica fue probada sobre 27 registros de ECG de la MIT-BIH Arrhythmia Database, que fueron elegidos intentando abarcar el mayor número de patologías 18 Capítulo 1. Introducción posible. Para extraer la información del ECG se utiliza una ventana de señal variable, alrededor de la anotación proporcionada por un detector de latidos. Posteriormente los latidos son normalizados y anotados manualmente por cardiólogos, lo que permite la verificación de los resultados. En la primera etapa el algoritmo redujo el número de latidos de 27412 a 1026. Estos latidos fueron posteriormente agrupados en la segunda etapa, logrando en el mejor de los casos un 7% de error con el algoritmo K-means. Otras aproximaciones más recientes también utilizan el algoritmo Max-Min para realizar el agrupamiento [134]. Este trabajo presenta una técnica novedosa para analizar la información del ECG a partir de un análisis simbólico. Para ello en una primera fase los latidos son detectados y agrupados. Posteriormente a cada grupo distinto de latidos se le asignará un símbolo identificativo. Con ello se consigue aumentar la abstracción del análisis, pasando de analizar la señal de ECG a analizar una secuencia de símbolos. Esto reduce la cantidad de datos y la dimensionalidad, facilitando el descubrimiento de patrones. Mediante el análisis simbólico es posible detectar cambios de ritmo y patrones temporales, entre otras anormalidades, y ver su relación con el estado del paciente y los fenómenos fisiopatológicos que concurren en el miocardio. El algoritmo de agrupamiento obtiene un error sobre la base de datos MIT-BIH Arrhythmia Database del 1.37%, utilizando 6 tipos de latidos y un número variable de grupos por registros (mediana de 22 grupos por registro). En [26] se propone una técnica de agrupamiento específica para la detección y el análisis de latidos ventriculares prematuros. La técnica utiliza una única derivación del ECG, que es filtrada para eliminar el ruido de alta frecuencia y la deriva de la línea base. Se aplica un detector de complejos QRS para localizar la onda R. Alrededor de esta posición se extrae una ventana variable de señal, dependiente de la distancia entre latidos. Para permitir comparar latidos de distinta longitud se utiliza DTW (Dynamic Time Warping) y el agrupamiento se basa en una combinación de algoritmos jerárquicos y algoritmos basados en centroides. En la técnica el número de grupos es ajustado dinámicamente para cada registro. La técnica fue aplicada a los registros de la base de datos MIT-BIH Arrhythmia Database, pero únicamente a los latidos etiquetados como normales o ventriculares prematuros. Para cada registro se utilizaron entre 2 y 30 grupos, con un porcentaje de error por registro menor al 1% en la mayor parte de los registros. En [75] se propone una técnica que combina agrupamiento y clasificación de complejos QRS, utilizando para el agrupamiento algoritmos basados en optimización mediante colonia de hormigas y para la clasificación redes neuronales y el método de los k vecinos más 1.3. Bases de datos de electrocardiografía 19 cercanos. Se utiliza un algoritmo adaptativo para detectar los complejos QRS, que son posteriormente normalizados. Después se realiza el agrupamiento y basándose en los resultados del agrupamiento se lleva a cabo la clasificación. También se produce una retroalimentación de la clasificación al agrupamiento por lo que en la etapa de clasificación se puede decidir volver a recalcular el agrupamiento. La técnica fue aplicada a un subconjunto de la base de datos MIT-BIH Arrhythmia Database obteniendo, en el mejor de los casos, una sensibilidad media del 94.4%. Finalmente, en [20] se propone un algoritmo dinámico de agrupamiento de latidos basado en plantillas que representan los grupos que se crean, se unen, se eliminan y evolucionan con el tiempo. En este trabajo se utilizó la base de datos MIT-BIH Arrhythmia Database completa. Para extraer los complejos QRS se utiliza una ventana de tamaño fijo de 200 ms. Al realizar comparaciones entre latidos se utiliza DTW para alinearlos. En el agrupamiento cada grupo se representa por una plantilla. El algoritmo se ejecuta dinámicamente, por lo que para cada nuevo latido se debe decidir si unirlo a un grupo ya existente o crear un grupo nuevo. Posteriormente, se produce una actualización de los grupos afectados y puede producirse alguna fusión de grupos. Las comparaciones entre los latidos y las plantillas se realizan basándose en el análisis y comparación de puntos relevantes en el trazado del ECG. Adicionalmente, además de comparar la morfología se puede añadir una etapa al agrupador en la que se compara también características derivadas de la distancia entre latidos. El algoritmo propuesto tiene una latencia mínima, siendo el primero de la bibliografía en proponer una solución completa que puede funcionar en tiempo real. Mediante la aplicación de este algoritmo a la base de datos MIT-BIH Arrhythmia Database completa se obtiene un error del 1.44% utilizando los 17 tipos de la base de datos. 1.3. Bases de datos de electrocardiografía En esta sección introduciremos las bases de datos electrocardiográficas que se usarán para validar las técnicas desarrolladas en esta tesis. 1.3.1. MIT-BIH Arrhythmia Database Se utilizará la base de datos más referenciada en la literatura de identificación de arritmias: la MIT-BIH Arrhythmia Database [96]. La gran variedad de pacientes, los distintos tipos de 20 Capítulo 1. Introducción latidos y la gran cantidad de anotaciones, han convertido esta base de datos en un “goldstandard” [16][30][31][77][108][112][152]. Esta base de datos está compuesta de 48 registros de electrocardiograma obtenidos de 47 pacientes distintos. Cada registro consta de dos derivaciones entre las siguientes: MLII, V1, V2, V3, V4 y V5. Los registros están digitalizados a una frecuencia de muestreo de 360 Hz con una resolución de 11 bits. La base de datos no se debe considerar como una muestra representativa de la población ya que los registros fueron seleccionados cuidadosamente para intentar abarcar la mayor variedad de trastornos cardiacos posible. Cada latido fue anotado por al menos dos cardiólogos, siendo aproximadamente el 68% de ellos considerados normales mientras que el resto se dividieron entre 16 tipos de latidos anormales (ver Tabla 1.1). Los latidos en los que solo aparece la onda P (p) no serán utilizados en este trabajo al no tener complejo QRS. La base de datos MIT-BIH Arrhythmia Database es anterior a la recomendación de la AAMI para etiquetar latidos [37]. Por ello utiliza sus propias etiquetas, dividiendo los latidos en 17 tipos distintos. La correspondencia entre las etiquetas de la base de datos MIT-BIH Arrhythmia Database y las etiquetas de la AAMI se muestra en la Tabla 1.2 (para evitar confusiones se mantienen los nombres originales en inglés de los tipos de latidos). En [30] se propone una correspondencia ligeramente distinta, que ha sido adoptada y utilizada por varios autores [33][76][85]. Sin embargo, siguiendo las recomendaciones de la AAMI, la correspondencia propuesta aquí es la correcta. Latidos auriculares de escape (e) y latidos nodales de escape (j) deberían ser asignados al tipo S de la AAMI en vez de al tipo N como se hace en [30]. 1.3.2. St.-Petersburg Institute of Cardiological Technics 12-lead Arrhythmia Database La base de datos “St.-Petersburg Institute of Cardiological Technics 12-lead Arrhythmia Database” (INCARTDB) consta de 75 registros muestreados a una frecuencia de 257 Hz provenientes de 32 pacientes distintos. La variedad de trastornos cardiacos no es tan amplia como en la MIT-BIH Arrhythmia Database, pero a cambio tiene un número mayor de registros y más derivaciones. Es una de las pocas bases de datos de la bibliografía con 12 derivaciones y anotada latido a latido, por lo que su uso está muy extendido [7][55][84][89]. Los registros fueron obtenidos de pruebas Holter realizadas para detectar la arterioesclerosis coronaria, por lo que en muchos de ellos se observan latidos ventriculares 1.3. Bases de datos de electrocardiografía 21 Tabla 1.1: Número de latidos clasificados por tipo en la base de datos MIT-BIH para cada registro. (Fuente: [104]) Record . L R A a J S V F ! e j E P f p Q 100 2239 - - 33 - - - 1 - - - - - - - - - 101 1860 - - 3 - - - - - - - - - - - - 2 102 99 - - - - - - 4 - - - - - 2028 56 - - 103 2082 - - 2 - - - - - - - - - - - - - 104 163 - - - - - - 2 - - - - - 1380 666 - 18 105 2526 - - - - - - 41 - - - - - - - - 5 106 1507 - - - - - - 520 - - - - - - - - - 107 -- - - - - - 59 - - - - - 2078 - - - 108 1739 - - 4 - - - 17 2 - - 1 - - - 11 - 109 -2492 - - - - - 38 2 - - - - - - - - 111 -2123 - - - - - 1 - - - - - - - - - 112 2537 - - 2 - - - - - - - - - - - - - 113 1789 - - - 6 - - - - - - - - - - - - 114 1820 - - 10 - 2 - 43 4 - - - - - - - - 115 1953 - - - - - - - - - - - - - - - - 116 2302 - - 1 - - - 109 - - - - - - - - - 117 1534 - - 1 - - - - - - - - - - - - - 118 -- 2166 96 - - - 16 - - - - - - - 10 - 119 1543 - - - - - - 444 - - - - - - - - - 121 1861 - - 1 - - - 1 - - - - - - - - - 122 2476 - - - - - - - - - - - - - - - - 123 1515 - - - - - - 3 - - - - - - - - - 124 -- 1531 2 - 29 - 47 5 - - 5 - - - - - 200 1743 - - 30 - - - 826 2 - - - - - - - - 201 1625 - - 30 97 1 - 198 2 - - 10 - - - 37 - 202 2061 - - 36 19 - - 19 1 - - - - - - - - 203 2529 - - - 2 - - 444 1 - - - - - - - 4 205 2571 - - 3 - - - 71 11 - - - - - - - - 207 -1457 86 107 - - - 105 - 472 - - 105 - - - - 208 1586 - - - - - 2 992 373 - - - - - - - 2 209 2621 - - 383 - - - 1 - - - - - - - - - 210 2423 - - - 22 - - 194 10 - - - 1 - - - - 212 923 - 1825 - - - - - - - - - - - - - - 213 2641 - - 25 3 - - 220 362 - - - - - - - - 214 -2003 - - - - - 256 1 - - - - - - - 2 215 3195 - - 3 - - - 164 1 - - - - - - - - 217 244 - - - - - - 162 - - - - - 1542 260 - - 219 2082 - - 7 - - - 64 1 - - - - - - 133 - 220 1954 - - 94 - - - - - - - - - - - - - 221 2031 - - - - - - 396 - - - - - - - - - 222 2062 - - 208 - 1 - - - - - 212 - - - - - 223 2029 - - 72 1 - - 473 14 - 16 - - - - - - 228 1688 - - 3 - - - 362 - - - - - - - - - 230 2255 - - - - - - 1 - - - - - - - - - 231 314 - 1254 1 - - - 2 - - - - - - - 2 - 232 -- 397 1382 - - - - - - - 1 - - - - - 233 2230 - - 7 - - - 831 11 - - - - - - - - 234 2700 - - - - 50 - 3 - - - - - - - - - 22 Capítulo 1. Introducción Tabla 1.2: Correspondencia entre las anotaciones de la MIT-BIH Arrhythmia Database y las recomendadas por la AAMI. Anotación AAMI Tipo de latido MIT-BIH (código) N Normal beat (N) Left bundle branch block beat (L) Right bundle branch block beat (R) S Aberrated atrial premature beat (a) Supraventricular premature beat (S) Atrial premature beat (A) Nodal (junctional) premature beat (J) Nodal (junctional) escape beat (j) Atrial escape beat (e) V Ventricular flutter wave (!) Ventricular escape beat (E) Premature ventricular contraction (V) F Fusion of ventricular and normal beat (F) Q Paced beat (/) Unclassifiable beat(Q) ectópicos. Ninguno de los pacientes tiene marcapasos y la edad media de los pacientes es de 58 años. Cada registro es un extracto de 30 minutos de duración y contiene las 12 derivaciones estándar. Todos los registros fueron anotados por un algoritmo automático y posteriormente estas anotaciones fueron corregidas manualmente. Las anotaciones de esta base de datos siguen la clasificación de latidos propuesta por la AAMI. CAPÍTULO 2 REPRESENTACIÓN DE COMPLEJOS QRS UTILIZANDO FUNCIONES DE HERMITE La representación computacional de los datos que forman parte de un problema tiene un impacto considerable a la hora de resolver dicho problema. La elección de una representación inadecuada puede complicar e incluso impedir su resolución. En el caso que nos atañe, el agrupamiento morfológico de latidos, debemos elegir una representación para los latidos que, siendo lo más sencilla y compacta posible, permita resolver satisfactoriamente el problema del agrupamiento. Este primer paso, previo al análisis, influirá en el resto del proceso. Errores en esta fase pueden invalidar las subsecuentes fases y por tanto el resultado del análisis. Si la representación elegida no captura suficientes características del latido para permitir su correcto agrupamiento no habrá forma de corregir la elección en etapas posteriores. Si la representación captura matices irrelevantes o innecesarios su almacenamiento y procesado se complicará, y estos matices pueden añadir ruido que dificulte el agrupamiento correcto. En la bibliografía podemos encontrar diversas propuestas de representación del latido, pero la mayoría se pueden agrupar en tres aproximaciones: 1. Señal. Para representar el latido se utiliza un fragmento del electrocardiograma en forma de señal digital, es decir, una señal discreta en tiempo discreto. Esta señal puede haber sufrido algún tipo de procesamiento previo con el objetivo de eliminar el ruido de la red eléctrica, suprimir la deriva de línea base o eliminar artefactos de alta frecuencia [62][131]. En algunos casos también se aplica una etapa de submuestreo 30 Capítulo 2. Representación de complejos QRS utilizando funciones de Hermite más alejado del valor de la media dentro de dicha ventana, que generalmente será el pico de la onda R. Esta será la nueva posición de la anotación del latido. El algoritmo de corrección de las anotaciones se puede aplicar en una derivación del ECG y emplear las posiciones obtenidas para las restantes, o puede aplicarse sobre cada derivación por separado. Si se aplica únicamente a una derivación, se asume que la posición del pico de la onda R es igual en el resto de derivaciones (lo que no es necesariamente cierto en la práctica). De aplicarse el algoritmo a cada derivación, la localización en el tiempo del latido puede ser ligeramente distinta para cada derivación. En este capítulo probaremos tres estrategias: usar las anotaciones originales de la base de datos; corregir la posición sobre una única derivación y usar dicha posición sobre el resto de derivaciones; y corregir la posición sobre cada derivación independientemente. 2.1.3. Cálculo de la representación de Hermite El complejo QRS es la parte central del latido cardiaco y su característica más importante. Para el problema que nos ocupa, el del agrupamiento de latidos, la onda T proporciona poca información adicional si ya tenemos en cuenta la proporcionada por el complejo QRS. La onda P sí que proporciona información adicional que ayuda a distinguir entre distintas arritmias, como las contracciones auricular y auriculoventricular prematuras, y los latidos de escape auriculares y auriculoventriculares. El problema en este caso reside en identificar y extraer la onda P de forma precisa y estable [34]. Este hecho lleva a que habitualmente se intente obtener información equivalente de otra forma más robusta, siendo la opción más habitual emplear información derivada de medir la distancia entre latidos consecutivos [30][77][139]. Por tanto, en esta sección nos centraremos en la representación del complejo QRS mediante las funciones de Hermite, asumiendo que esta representación se complementará con información extraída de la distancia entre latidos. Para el cómputo de la representación se parte de las anotaciones que sitúan el latido en un instante de tiempo. Alrededor de estas anotaciones se extrae una ventana de 200 ms para cada complejo. El tamaño de esta ventana es suficientemente grande para contener la anchura del complejo QRS de un latido ventricular, pero al mismo tiempo suficientemente estrecha como para dejar fuera las ondas P y T. Este valor para la anchura de esta ventana es uno de los más comunes en la bibliografía [77]. Las funciones de Hermite convergen a cero en ±∞. Para asegurar esta convergencia se añaden 100 ms de ceros en ambos extremos de la ventana de 200 ms que contiene el complejo 2.1. Material y métodos 31 QRS. Se obtiene por tanto una ventana de 400 ms, x[l], representada como: x[l] = N−1 ∑ n=0 cn(σ)φn[l,σ)+ e[l], l=−W·fs 2,−W·fs 2+1,...,W·fs 2, (2.2) siendo Nel número de funciones de Hermite utilizadas, Wel tamaño de la ventana en segundos y fsla frecuencia de muestreo de la señal. Los símbolos b c significan que el valor es redondeado al entero menor más cercano. φn[l,σ)es la nfunción discreta de Hermite obtenida discretizando la función continua φn(t,σ)a una frecuencia fs;cnson los coeficientes de la combinación lineal que representa nuestro latido; e[l]es el error entre x[l]y la representación de Hermite; y σes un parámetro de elongación que controla la anchura de la función de Hermite, permitiendo ajustarla a la anchura del complejo QRS. Las funciones de Hermite φn[l,σ), 0 ≤n<N, son definidas como: φn[l,σ) = 1 pσ2nn!√πe−(l·Ts)2/2σ2Hn(l·Ts/σ)(2.3) siendo Tsel inverso de la frecuencia de muestreo. De esta forma, cada complejo QRS se representa mediante los Ncoeficientes cn(σ), 0 ≤ n<N, y el valor de σ. Podemos ver en la Figura 2.5 cómo varía la representación utilizando distintos valores de N. Generalmente, usar más funciones implica una representación más precisa. Sin embargo, utilizar un mayor número de funciones también conlleva el riesgo de sobreajustar y modelar el ruido de la señal y no la forma real del complejo QRS. En la misma figura se puede ver este comportamiento cuando se utilizan 15 funciones. Para un cierto valor de σ, las funciones de Hermite forman una base ortonormal: Z∞ −∞ φn(σ)φm(σ) = δmn.(2.4) Esto permite un cálculo eficiente de cn(σ). Sin una ventana infinita, o en el caso de funciones discretas, (2.4) no se cumple. Sin embargo, si φn[l,σ)es suficientemente cercano a cero en los extremos y fuera de la ventana, (2.4) resulta una aproximación aceptable. Para los límites de la ventana usaremos el criterio (bastante tolerante) de que φn[l,σ)sea como mucho 1/10 de su máximo valor dentro de la ventana: |φn[−l0,σ)|=|φn[l0,σ)|<1 10 max l∈[−l0,l0]|φn[l,σ)|,(2.5) 32 Capítulo 2. Representación de complejos QRS utilizando funciones de Hermite (a) Latido Original (b) N=3 (c) N=6 (d) N=9 (e) N=12 (f) N=15 Figura 2.5: Latido original y representación de Hermite con N=3, 6 ,9 ,12 y 15 para un mismo valor de σ. donde −l0yl0son respectivamente la primera y última muestra de la ventana. Requeriremos también que el valor de φn[l,σ)fuera de la ventana sea menor que en el límite de la ventana: |φn[l,σ)|≤|φn[l0,σ)| ∀|l|>l0.(2.6) Para un determinado tamaño de ventana y un número fijo de funciones de Hermite, (2.5) y (2.6) imponen un límite máximo en el valor de σ. Por ejemplo, para N=3, 4 y 5, y una ventada de 200 ms de señal, los valores máximos de σque cumplen (2.5) y (2.6) son 62 ms, 55 ms y 51 ms, respectivamente. Una vez fijado el valor de σpodemos calcular los coeficientes cn(σ)minimizando la suma cuadrática del error: ∑ l (e[l])2=∑ l x[l]− N−1 ∑ n=0 cn(σ)φn[l,σ)!2 .(2.7) 2.1. Material y métodos 33 El mínimo de este error cuadrático podemos aproximarlo a partir de la propiedad de ortogonalidad si suponemos que: cn(σ) =~x·~ φn(σ),(2.8) dónde los vectores,~xy~ φn(σ), son definidos como~x={x[l]}y~ φn(σ) = {φn[l,σ)}. Por tanto, podemos obtener el mejor σmediante un proceso iterativo [77]. Para ello se realizará un incremento iterativo de σrecalculando para cada iteración (2.8) y (2.7). En este proceso σcomienza desde 0 y se va incrementando hasta el valor máximo (fijado por (2.5) y (2.6)) en incrementos de fs 1000. Finalmente seleccionaremos el valor que minimiza la suma cuadrática del error. Los valores óptimos de σpara cualquier Nen la base de datos MIT-BIH Arrhythmia Database están mayormente en un rango entre 14 y 21 ms, bastante por debajo que el límite máximo marcado. Terminado este paso ya tenemos los valores que conformarán la representación del complejo QRS: el valor de σy los valores de los coeficientes de las funciones de Hermite c0(σ),c1(σ),c2(σ),...,cN−1(σ). Si en nuestro problema se van a considerar varias derivaciones en el electrocardiograma, cada derivación será procesada independientemente, obteniendo una representación para cada complejo QRS. 2.1.4. Error en la representación Una vez obtenida la representación del complejo QRS es necesario estudiar el error entre la representación y la señal original. Para ello existen múltiples métricas muy variadas. En [77] para calcular el error en la representación de Hermite se utiliza el cociente entre la energía del error cuadrático (2.7), y la energía de la señal original: Err =∑l(e[l])2 ∑l(x[l])2.(2.9) Esta es una métrica del error directa y sencilla pero su valor resulta de difícil interpretación. No obstante, nos servirá para comparar nuestros resultados con los obtenidos en dicho trabajo. Además, calcularemos la raíz del error cuadrático medio normalizada (NRMSD): NRMSD =RMSD xmax −xmin =q∑t(e[l])2 V xmax −xmin ,(2.10) donde Ves el tamaño de la ventana en muestras y xmin yxmax son, respectivamente, el valor mínimo y máximo de la señal en dicha ventana. Esta medida tiene una interpretación más 34 Capítulo 2. Representación de complejos QRS utilizando funciones de Hermite intuitiva que (2.9), siendo el error promedio expresado como un porcentaje del rango de valores de la señal. 2.1.5. Selección de la representación óptima del complejo QRS Añadir más funciones a la representación del latido siempre reduce el error, obteniendo por tanto una representación más exacta, pero a costa de incrementar la dimensionalidad del vector de características que representa el latido y los requerimientos computacionales para su procesamiento. En la bibliografía no hay ningún criterio bien definido para decidir cuál es el número adecuado de funciones a utilizar en la representación, sino que generalmente esta tarea se lleva a cabo mediante una inspección visual del resultado de la reconstrucción de unos pocos latidos. Dado que la representación se basa en un modelo matemático, en esta sección aplicaremos técnicas de la literatura de selección de modelos [17][18][22] que proporcionan un criterio objetivo para la selección del número óptimo de funciones para la representación del latido. En el campo de las técnicas de selección de modelos encontramos principalmente dos aproximaciones: criterios bayesianos y criterios basados en la teoría de la información. Entre las aproximaciones basadas en criterios bayesianos destaca el criterio de información bayesiana (Bayesian Information Criterion, BIC) [124]. Las aproximaciones basadas en teoría de la información son muchas y variadas, pero podemos destacar por su extendido uso el criterio de información de Akaike (Akaike Information Criterion, AIC) [5]. BIC está basado en la función de verosimilitud: L=p(x|ˆ θ,M),(2.11) donde ˆ θson los parámetros, xlos datos observados y Mel modelo. Al añadir parámetros adicionales la verosimilitud siempre tiende a aumentar. Sin embargo, si dejamos que la verosimilitud crezca sin límites estaríamos sobreajustando el modelo. Para resolver esta limitación, BIC añade un término de penalización basado en el número de parámetros del modelo. En el caso que nos atañe cada nuevo parámetro es un nuevo coeficiente, resultado de añadir una función de Hermite a la representación. Por ejemplo, si utilizamos funciones hasta φN−1[l,σ)tendremos N+1 parámetros, los Ncoeficientes c0(σ),c1(σ),c2(σ),...,cN−1(σ) de la combinación lineal y σ. Por tanto, para nuestro problema BIC se formula como: BIC =−2·ln ˆ L+k·ln(m),(2.12) 2.1. Material y métodos 35 dónde mes el número de muestras, kes el número de parámetros (k=N+1 cuando usamos Nfunciones) y ˆ Les el valor máximo de la función de verosimilitud. AIC está basado en la distancia de Kullback-Leibler y en la entropía de la información. Mediante su uso podemos obtener una estimación de la distancia relativa entre un modelo y el proceso desconocido que verdaderamente generó los datos observados. En la práctica se usa para comparar distintos modelos y seleccionar el que presente la menor distancia, pero no permite saber nada acerca de la calidad del modelo en sentido absoluto. Si todos los modelos candidatos se adaptan pobremente a los datos, AIC no lo detectará. Akaike en su trabajo definió AIC como: AIC =−2·ln ˆ L+2·k.(2.13) Esta formulación clásica de AIC es poco adecuada para los casos en que el tamaño muestral (m) no sea varios órdenes de magnitud mayor que el número de parámetros (k), ya que tiende a sobreajustar los datos. Sugiura [133] propuso una variante de segundo orden, AICc, con un término de corrección adicional para paliar este efecto. En el caso de que msea suficientemente grande AICcconverge a AIC y por tanto seleccionarán el mismo modelo. Esta corrección AICcse define como: AICc=AIC+2k(k+1) m−k−1.(2.14) Dado que en nuestro problema el cociente m/kno es grande, haremos uso de esta variante de AIC. Tanto AIC como BIC tienen sus ventajas e inconvenientes. AIC tiende a elegir modelos más complicados de lo estrictamente necesario, especialmente si el número de muestras es pequeño [17]. Pero al estar basado en una minimización del error cuadrático medio se convierte en un candidato idóneo para un problema como el que nos atañe. Por otra parte, BIC debería, en teoría, ser aplicado únicamente para encontrar el modelo real entre un conjunto de modelos candidatos que lo contiene; condición que no se cumple en nuestro caso. Además, BIC muestra una tendencia a seleccionar modelos demasiado simples si el número de muestras es pequeño, justo al contrario que AIC. Aprovechando que ambos criterios abordan el mismo problema desde distintas aproximaciones, aplicaremos tanto AIC como BIC por separado para luego comparar sus resultados [148]. Bajo las hipótesis de normalidad e independencia de los errores, y asumiendo además que la varianza es constante, BIC y AICcpueden ser expresados como [17]: BIC =m·ln(σ2 e)+k·ln(m),(2.15) 36 Capítulo 2. Representación de complejos QRS utilizando funciones de Hermite AICc=m·ln(σ2 e)+2·k+2k(k+1) m−k−1,(2.16) dónde σ2 ees la varianza del error. El estimador insesgado c σ2 eserá usado para estimar σ2 e, que en nuestro caso adopta la siguiente forma: c σ2 e=1 m−1 m ∑ l=1 x[l]− N−1 ∑ n=0 cn(σ)φn[l,σ)!2 =1 m−1 m ∑ l=1 e[l]2.(2.17) Cada complejo QRS de cada derivación se considera un problema de selección de modelos independiente. Para cada uno de ellos calcularemos BIC y AICcusando desde 2 hasta 30 funciones de Hermite para construir la representación. Posteriormente, empleando cada uno de los dos criterios se selecciona el número óptimo de funciones de Hermite a utilizar. 2.1.6. Selección de características Cada complejo QRS es representado por los coeficientes de la combinación lineal de funciones cn(σ)y por el parámetro de anchura σ. Si tenemos en cuenta que generalmente trabajaremos al menos con dos derivaciones del ECG y que podemos llegar a usar hasta 30 funciones para representar el complejo QRS, podríamos tener hasta 62 características para representar cada latido, o incluso más en el caso de utilizar un mayor número de derivaciones (362 para 12 derivaciones en la INCARTDB). En este contexto, podría tener sentido intentar seleccionar un subconjunto de aquellas características que son más relevantes para separar unas familias morfológicas de latidos de otras. Con este fin se aplicaron varias técnicas de selección de características que permiten ordenar las características según la cantidad de información que aportan. Las técnicas utilizadas fueron dos basadas en la información mutua “InformationGain” y “GainRatio” y una basada en el test Chi-cuadrado. Las dos primeras tratan de medir la información que aporta cada nueva característica para tomar decisiones (InformationGain) [115], aplicando GainRatio, una corrección sobre la anterior para penalizar decisiones que den lugar a un gran número de bifurcaciones en la decisión. Por otra parte, el método basado en Chi-cuadrado ordena las características basadas en el valor del estadístico chi-cuadrado con respecto a la clase [91]. Para aplicar estas técnicas se utilizó la base de datos MIT-BIH al completo. La selección de características fue aplicada para cada registro (paciente) de forma independiente y sobre 2.2. Resultados 37 ● ● ● ●● ● ●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ●● ●●●●●●●●●●●●●●●●●●●●●●●● D1 D2 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 0% 1% 2% 3% 4% 5% 6% 7% 8% 0% 1% 2% 3% 4% 5% 6% 7% 8% Número de Funciones Valor del Error NRMSD Anotaciones ●D1&D2 D1 Originales Figura 2.6: Resultados del error NRMSD por derivación para la señal filtrada para las tres estrategias de corrección de posición del latido: sin corrección (Originales), corrección sobre la primera derivación (D1) y corrección sobre las dos derivaciones (D1&D2). toda la base de datos en su conjunto. De esta forma podemos comparar los resultados de ambas aproximaciones. 2.2. Resultados Los algoritmos descritos en este capítulo fueron implementados en el lenguaje Java, con la excepción del filtrado de línea base y el filtrado de alta frecuencia. Dichos filtros fueron implementados en Matlab y desde ese entorno se generó un conjunto completo de los registros de la base de datos MIT-BIH filtrados. Posteriormente se procesaron tanto los registros filtrados como los registros originales con los algoritmos implementados en Java. 38 Capítulo 2. Representación de complejos QRS utilizando funciones de Hermite 2.2.1. Error en la representación Los resultados del error NRMSD (2.10) se muestran en las Figuras 2.6 y 2.7. En dichas figuras se muestra el error de cada derivación; las barras muestran la desviación estándar de dicho error. En la Figura 2.6 se muestran los resultados con la señal filtrada mientras que en la Figura 2.7 se muestran con la señal sin filtrar; correspondiendo las gráficas de la parte superior a la primera derivación y las de la parte inferior a la segunda derivación en ambos casos. Para cada caso aparecen tres valores de error correspondientes con cada una de las tres estrategias de corrección de posición de latido recogidas en la Sección 2.1.2: utilizar las anotaciones originales de la base de datos; corregir la posición solo sobre la primera derivación (D1); y corregir la posición sobre ambas derivaciones por separado (D1&D2). Se eligió la primera derivación al corregir las anotaciones sobre la posición de una única derivación debido a que originalmente fue la derivación utilizada para anotar la base de datos MIT-BIH. Estas tres aproximaciones se muestran en las figuras en distintas series de datos utilizando cuadrados, triángulos y círculos respectivamente. En otra gráfica (véase Figura 2.8) se muestran los resultados del error calculado según la ecuación (2.9); el error mostrado es la media del error en cada una de las dos derivaciones. En este caso se muestra en la misma figura los resultados con la señal filtrada (arriba) y con la señal sin filtrar (abajo) utilizando las mismas tres estrategias de corrección de posiciones de los latidos que en las figuras anteriores. 2.2.2. Representación óptima según AIC y BIC En la Figura 2.9 podemos ver el porcentaje de latidos de la base de datos MIT-BIH cuyos complejos QRS son óptimamente representados, de acuerdo a BIC y AICc, con No menos funciones de Hermite. Por ejemplo, el 57% de los complejos son representados óptimamente con 20 o menos funciones de Hermite, de acuerdo tanto a BIC como a AICc. Ningún complejo QRS se representó óptimamente con un número menor de 8 funciones; algunos complejos fueron representados óptimamente utilizando 9 o 10 funciones, pero su número es muy bajo y por ello no son mostrados en la gráfica. Los resultados mostrados fueron calculados utilizando la señal filtrada, siendo muy similares a los resultados obtenidos sin filtrar. Es necesario especificar que de los latidos que aparecen como representados óptimamente con 30 funciones en la Figura 2.9 pueden en algunos casos requerir un mayor número de funciones para representarse de modo óptimo, por lo que se podría decir que el análisis se 2.2. Resultados 39 ● ● ● ●● ●●●●●●●●●●●●●●●●●●●●●●●● ● ● ● ●● ●●●●●●●●●●●●●●●●●●●●●●●● D1 D2 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 0% 1% 2% 3% 4% 5% 6% 7% 8% 0% 1% 2% 3% 4% 5% 6% 7% 8% Número de Funciones Valor del Error NRMSD Anotaciones ●D1&D2 D1 Originales Figura 2.7: Resultados del error NRMSD por derivación para la señal sin filtrar para las tres estrategias de corrección de posición del latido: sin corrección (Originales), corrección sobre la primera derivación (D1) y corrección sobre las dos derivaciones (D1&D2). realizó únicamente desde 2 hasta 29 funciones, aunque se calcularan también los valores para 30 funciones. Por tanto, para aquellos latidos que hubieran requerido un mayor número de funciones AICcy BIC seleccionarán la representación con 30 funciones como el modelo óptimo. Sin embargo, el 97% de los latidos según BIC y el 99% de los latidos según AICc están ya óptimamente representados con 29 o menos funciones. Por tanto, el número de latidos que se representarían óptimamente con más de 30 funciones sería, sin duda, pequeño. 2.2.3. Selección de características Los resultados obtenidos por los métodos de selección de características no fueron alentadores. Las características más relevantes que fueron seleccionadas por los tres métodos variaban considerablemente entre pacientes. Por otra parte, las características seleccionadas para cada paciente eran distintas a las seleccionadas trabajando con toda la base de datos CAPÍTULO 3 CÁLCULO DE LA REPRESENTACIÓN DE HERMITE UTILIZANDO GPUS En el capítulo anterior se estudió la representación del latido mediante las funciones de Hermite y se analizó el tiempo de ejecución del cálculo de dicha representación. Este tiempo puede llegar a suponer un problema de rendimiento, especialmente si se utiliza el número de funciones sugerido por AICcy BIC (véase Sección 2.3). El paralelismo es una opción habitual en la bibliografía para resolver problemas científicos complejos que requieren procesar grandes cantidades de datos. Entre las muchas aproximaciones disponibles destaca el uso de GPUs (“Graphics Processing Units”), siendo una de las opciones más atractivas y más usadas para la aceleración de procesos computacionales. Las GPUs son relativamente baratas, fáciles de instalar y, en muchos casos proporcionan tanta potencia de cómputo como cientos de procesadores CPU trabajando en paralelo, a la vez que presentan un consumo de energía significativamente menor. Esto facilita su integración en equipos de monitorización electrocardiográfica, que mediante una GPU pueden llegar a realizar, en un tiempo equivalente, los mismos cálculos que un clúster de computación formado por docenas, o incluso cientos, de ordenadores. Su uso presenta ventajas considerables tanto en el ámbito de la monitorización hospitalaria, como en la monitorización domiciliaria, donde los datos podrían enviarse al hospital para su procesamiento. Es fácil encontrar ejemplos de aplicaciones de la tecnología de GPUs en el campo de la biomedicina: reconstrucción MRI [74], simulación tejido cardiaco [49], dinámica molecular [151], bioinformática [102], entre otras muchas. 48 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs Por ello hemos considerado interesante el estudiar cómo acelerar el cálculo de la representación de Hermite mediante el uso de GPUs. Los dispositivos GPU están diseñados para realizar una misma tarea de forma repetitiva sobre una gran cantidad de datos. Si la tarea requiere ejecuciones con flujos condicionales o muestra una gran dependencia entre los datos no se da una situación propicia para el uso de GPUs y el rendimiento total puede no mejorar notablemente, sino incluso disminuir. Previamente, en el capítulo anterior, ya se aplicó una paralelización basada en procesamiento paralelo mediante hilos para calcular la representación de Hermite. Los buenos resultados obtenidos y la idoneidad del problema para su paralelización motivaron que se decidiera utilizar GPUs para optimizar la obtención de la representación de Hermite. En la operación del cálculo de la representación de Hermite cada complejo QRS en cada derivación representa una operación independiente, que no requiere memoria compartida y que además hay que repetir una gran cantidad de veces. Es más, con un manejo hábil de las operaciones resulta una tarea que es posible paralelizar incluso dentro del mismo complejo QRS, a nivel de muestra. Por todo ello, este problema es un escenario perfecto para sacar todo el partido a las GPUs y conseguir una aceleración considerable. 3.1. Programación utilizando GPUs Habitualmente, para programar las GPUs se utilizan lenguajes basados en C. Dichos lenguajes permiten diseñar e implementar los mecanismos de paralelización de forma fácil, permiten una compilación rápida y proporcionan una integración sencilla con otros sistemas. En este trabajo se decidió utilizar un lenguaje basado en C llamado CUDA (“Compute Unified Device Architecture”)[74]. El motivo principal de dicha decisión fue la facilidad de implementación y su extendido uso. El lenguaje CUDA realiza una encapsulación de los detalles del hardware, simplificando la tarea de desarrollo y permitiendo una mayor portabilidad entre distintas plataformas hardware. Un dispositivo GPU está formado por varios multiprocesadores llamados “streaming processors” (SP), cada uno a su vez compuesto por varios núcleos que trabajan en paralelo. La unidad básica de ejecución en la GPU se denomina “kernel”. La GPU ejecuta el mismo kernel en paralelo con conjuntos de datos distintos como entrada. Cada hilo de ejecución ejecuta una instancia concreta del kernel en paralelo con el resto de hilos. Los hilos se agrupan en “warps” y son gestionados por un SP. Todos los hilos de un warp comparten el 3.1. Programación utilizando GPUs 49 mismo código, siguen el mismo camino de ejecución y se espera que se detengan en el mismo punto. El SP se encarga de agrupar a todos los hilos en el mismo punto de ejecución y de que estén siempre sincronizados. Por lo tanto, la presencia de bloques condicionales puede deteriorar considerablemente el rendimiento ya que el SP debe esperar a que todos los hilos alcancen el mismo punto para seguir con la ejecución. En CUDA existen mecanismos que permiten controlar la ejecución y paralelización de los hilos. De esta forma el desarrollador puede controlar cómo se agruparán los hilos; una buena estrategia para diseñar estos grupos es fundamental para lograr aprovechar al máximo el paralelismo del procesador. Estos grupos que forma el desarrollador se denominan bloques y cada uno estará representado por un identificador (ID). De forma similar los bloques se pueden agrupar a su vez en grids, esto es, conjuntos de bloques que también tienen su propio identificador. Finalmente, cada hilo también tendrá un identificador que, combinado con el identificador del bloque y del grid, permite indexar los hilos de forma fácil y sencilla. Combinando estos identificadores cada hilo puede generar localizaciones en memoria únicas a las que acceder de forma privada. Durante la ejecución cada bloque es asignado a un SP que gestiona su ejecución. Respecto a la memoria, las GPUs tienen una jerarquía de memoria piramidal. La memoria de mayor tamaño es una memoria global (DRAM) a la que pueden acceder todos los hilos, seguida de una memoria compartida (SRAM) a la que pueden acceder todos los hilos de un mismo bloque y finalmente una memoria privada, que está formada por los registros privados de cada hilo de ejecución. Como suele ser habitual, cuanto mayor es la memoria mayor es el tiempo de acceso necesario. La memoria global, con una gran capacidad (1-6 GB) es la más lenta, mientras que la memoria compartida es pequeña (16-48 KB) pero mucho más rápida (unos dos órdenes de magnitud más rápida que la memoria global). En la memoria global el acceso es en bloque, ya que las operaciones de lectura y escritura manejan bloques de bits consecutivos (32, 64, 128, etc.). También se dispone de una caché de nivel 2 que optimiza el uso de la memoria global, pero ello no evita que se deba realizar una planificación apropiada del uso de la memoria. Una mala planificación conllevaría un uso ineficiente de la memoria global y por consiguiente un mayor tiempo de procesamiento. El acceso a la memoria compartida se realiza de forma aleatoria, por lo que mientras se observen ciertas restricciones de uso el acceso será rápido. Teniendo todo esto en cuenta, es posible diseñar un algoritmo para calcular una representación de Hermite de forma que aproveche al máximo las capacidades de una GPU. 50 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs Para conseguir esto es necesario realizar una cuidadosa elección de los parámetros y diseñar concienzudamente el flujo de ejecución y de datos del algoritmo. 3.2. Implementación optimizada en C Antes de pasar a diseñar el algoritmo para calcular la representación de Hermite del QRS utilizando CUDA, se realizará una implementación optimizada en C. Con esta implementación ejecutada en una CPU estableceremos una implementación de referencia con la que comparar la aceleración de la implementación en CUDA (CPU vs GPU). Esto nos permite establecer una referencia más justa para la comparación que la implementación previa en Java, que no estaba optimizada. Los cálculos necesarios para realizar la caracterización de los complejos QRS mediante Hermite y el pseudocódigo optimizado se muestran en el Algoritmo 1. Como entrada el algoritmo recibe el registro de ECG y el número de funciones de Hermite a utilizar en la representación. Las salidas son el conjunto de características que representarán al latido: los coeficientes de las funciones de Hermite que mejor representan el latido y la σutilizada. El Algoritmo 1 consta de tres partes principales: i) extracción de los complejos QRS, ii) precomputación de las funciones de Hermite (#Lazo1) y, iii) cálculo de los coeficientes de las funciones y de σ(#Lazo2). Para la extracción de los complejos QRS se usa el mismo procedimiento ya detallado en la Sección 2.1.3, tomando 144 muestras para representar cada complejo QRS (suponiendo frecuencia de muestreo de 360Hz). En el bucle (#Lazo1, líneas 4-6) se calculan los valores de las funciones de Hermite φn[l,σ)para todos los σposibles mediante la Ecuación (2.3). Estas funciones serán utilizadas intensivamente en el siguiente bucle, de ahí el interés en calcularlas previamente y almacenar el resultado, evitando así repetir el cálculo. Para medir el incremento en rendimiento de este cálculo se midió la diferencia entre el tiempo necesario para ejecutar el precomputado de las funciones junto con el segundo bucle (#Lazo2) y el tiempo necesario para ejecutar el segundo bucle (#Lazo2) calculando las funciones dentro del bucle. La alternativa de precalcular los valores consiguió una aceleración de 105x comparada con la otra opción, ejecutando las pruebas en un procesador Intel-i7. Este test demuestra que resulta beneficioso realizar este cálculo previamente. Por otra parte, en el bucle (#Lazo2, líneas 8-21) se realiza el cálculo para obtener el conjunto de características (los coeficientes cn(σ)y la σ) que mejor representan el latido. En 3.2. Implementación optimizada en C 51 Algoritmo 1 Caracterización del complejo QRS sin paralelizar [tbp] Input: Registro de ECG, N(número funciones a utilizar) Output: σycn(σ)que minimizan el error para cada latido 1: for all latido ido 2: Extraer señal del complejo QRS, xi[l], del ECG. 3: end for 4: for all nyσdo #Lazo1 5: Calcular φn[l,σ)#Ecuación (2.3) 6: end for 7: Errmin =∞ 8: for all latido i,xi[l]do #Lazo2 9: for all σdo 10: for all ndo 11: Calcular cn(σ)#Ecuación (2.8) 12: end for 13: Calcular ˆxi[l] = ∑N−1 n=0cn(σ)φn[l,σ)#Ecuación (2.2) 14: Calcular MSE Err=∑l(e[l])2=∑l(xi[l]-ˆxi[l])2#Ecuación (2.7) 15: if Err<Errmin then 16: Sigmabest =σ 17: Cbest =cn(σ) 18: Errmin =Err 19: end if 20: end for 21: end for esta parte tenemos tres bucles anidados. El más externo (línea 8) itera sobre todos los latidos extraídos. El segundo (línea 9) prueba iterativamente todos los posibles valores de σy el más interno calcula los coeficientes de las Nfunciones de Hermite utilizadas (véase Ecuación (2.8)). Tras este cálculo, el segundo bucle selecciona los valores σycn(σ)óptimos (aquellos que minimizan el MSE (2.7)). Analizando el Algoritmo 1 puede verse que tanto (#Lazo1) como (#Lazo2) permiten una paralelización casi completa. En la siguiente sección intentaremos paralelizar ambos bucles adaptando este algoritmo al uso de una GPU de la forma más optimizada posible. 52 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs Algoritmo 2 Caracterización del complejo QRS paralelizada con la GPU [tbp] Input: Registro de ECG, N(número funciones a utilizar) Output: σycn(σ)que minimizan el error para cada latido 1: for all latido ido 2: Extraer señal del complejo QRS, xi[l], del ECG. 3: end for 4: Reservar memoria en la GPU 5: kernel_φ(N) 6: Enviar los complejos xi[l]a la GPU #Escribir en la memoria de la GPU 7: Esperar a que la GPU termine de procesar 8: kernel_Hermite(xi[l],N) #Algoritmo 3 9: Esperar a que la GPU termine de procesar 10: Leer σycn(σ)para cada latido #Escribir en la memoria del ordenador el resultado 3.3. Implementación paralela Partiendo del código del Algoritmo 1, se realizó una paralelización de (#Lazo1) y (#Lazo2) por medio de dos kernels:kernel_φ, encargado de precalcular los valores de las funciones de Hermite, y kernel_Hermite, encargado de calcular el valor de los coeficientes y el valor del error. El nuevo flujo de ejecución se muestra en el Algoritmo 2. En este caso el algoritmo fue diseñado para aprovechar al máximo las capacidades de la GPU. El primer paso en la versión paralela es extraer los complejos QRS y reservar la memoria necesaria para almacenarlos en la GPU (Algoritmo 2, líneas 1-4). Tras esto, se llama al primer kernel (kernel_φ) para generar los valores de las funciones de Hermite (Algoritmo 2, línea 5) y los resultados son almacenados en la GPU. Dicho kernel será ejecutado en la GPU y al mismo tiempo, mientras está siendo ejecutado, los latidos son enviados a la memoria global de la GPU (Algoritmo 2, línea 6). Posteriormente, cuando la GPU termine de procesar el primer kernel se pasará al segundo (kernel_Hermite) (Algoritmo 2, línea 8). Tras terminar esta ejecución los resultados finales son enviados a la memoria del ordenador (Algoritmo 2, líneas 9-10). Este será el esquema general de procesado del algoritmo en la GPU. En las siguientes secciones se concretará en detalle la función de cada uno de los dos kernels, además de presentar la estrategia empleada para optimizar la transferencia de datos para el procesado de latidos en tiempo real. 3.3. Implementación paralela 53 BLOCK(0,0) σ=σ0 TH0 Φ0[0, σ0]TH143 Φ0[143, σ0] ... ... BLOCK(0,N-1) σ=σ0 TH0 ΦN-1[0, σ0]TH143 ΦN-1[143, σ0] ... BLOCK(S-1,N-1) σ=σS-1 TH0 ΦN-1[0, σS-1]TH143 ΦN-1[143, σS-1] ... ... BLOCK(S-1,0) σ=σS-1 TH0 Φ0[0, σS-1]TH143 Φ0[143, σS-1] ... ... ... Figura 3.1: Distribución de hilos en bloques para el kernel_φ. 3.3.1. Precomputación de las funciones de Hermite El kernel_φ, equivalente al (#Lazo1) del Algoritmo 1, es paralelizable directamente sin necesidad de realizar cambios en el código. En el caso de trabajar a 360 Hz, como en la base de datos MIT-BIH Arrhythmia Database, la configuración óptima utiliza bloques de 144 hilos, correspondientes a las 144 muestras que representan cada complejo QRS. Los bloques se distribuirán en un grid de dimensiones S×N, siendo Sel número de posibles valores que toma σyNel número de funciones de Hermite a calcular. Cada bloque se encarga del cálculo de todas las muestras para una función de Hermite φn[l,σ)con un valor concreto de σyn. El bloque tendrá un identificador bidimensional, siendo su primera componente un índice referido al valor de σdentro de S que debe calcular y su segunda componente un índice referido a la función de Hermite. Cada hilo dentro de un bloque se encarga de evaluar la Ecuación (2.3) para una muestra distinta (correspondiente al identificador del hilo). De esta forma se consigue una implementación completamente paralela, calculando simultáneamente el valor de tantas funciones como SP haya en la GPU (véase Figura 3.1). 3.3.2. Caracterización del complejo QRS mediante Hermite En el Algoritmo 1 el núcleo de la caracterización del latido está en el (#Lazo2) donde se obtienen las características que representan al complejo QRS. La paralelización de este 54 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs proceso en el kernel_Hermite requerirá de un análisis más pormenorizado que el realizado para el kernel_φ. En este caso cada bloque trabajará sobre un complejo QRS concreto y tendremos por tanto tantos bloques como latidos en el registro (aproximadamente unos 2000 latidos por registro de la base de datos MIT-BIH Arrhythmia Database). Posteriormente se abordará cómo procesar registros de larga duración como los Holter. Cada bloque deberá contar con tantos hilos como muestras en el complejo QRS (144 por bloque con una frecuencia de muestreo de 360 Hz). Con esta estructura sabremos en todo momento sobre qué latido se están realizando los cálculos (ID del bloque) y sobre qué muestra (ID del hilo). El Algoritmo 3 muestra el pseudocódigo desarrollado para ejecutar el kernel_Hermite. Este mismo código se ejecuta para todos los hilos, por ello lo primero será identificar sobre qué latido y qué muestra se va a trabajar (Algoritmo 3, líneas 1-2). Seguidamente la información del complejo QRS será copiada a la memoria compartida, ya que va a ser utilizada múltiples veces. Esto implica que toda referencia posterior a xi[l]será realizada mediante una lectura rápida de la memoria compartida. De forma análoga al Algoritmo 1, se utiliza un bucle para recorrer todos los posibles valores de σ(Algoritmo 3, líneas 5-23). Para cada valor se calcula el vector de coeficientes cn(σ)(Algoritmo 3, líneas 6-11). Se debe tener en cuenta que la multiplicación (Algoritmo 3, línea 7) se realiza de forma totalmente paralela y que el resultado se guarda en una variable compartida entre los hilos. Posteriormente los resultados parciales obtenidos por cada uno de los hilos se suman para obtener el coeficiente de la función de Hermite; esta suma (Algoritmo 3, línea 10) se realiza mediante la técnica de suma por reducción [74], para evitar realizarla de forma completamente secuencial. Posteriormente, se calcula la representación del latido obtenida a partir de los coeficientes cn(σ)(Algoritmo 3, líneas 12-15). Con esta representación se calculará el error MSE (Algoritmo 3, líneas 16-17). De nuevo la multiplicación es paralela y la suma se realiza mediante la técnica de reducción. Finalmente, solo el hilo 0 actualiza la mejor solución si el MSE calculado es menor que el mínimo actual (Algoritmo 3, líneas 18-22). 3.3.3. Optimización de la transferencia de datos El flujo de datos presentado en el Algoritmo 2 puede ser optimizado si se tiene en cuenta la capacidad de la GPU de realizar al mismo tiempo un cálculo y una transferencia de datos (de la memoria del ordenador a la GPU o de la GPU al ordenador). Esto resulta especialmente interesante al lidiar con el análisis de registros de ECG muy grandes (por ejemplo, registros 3.3. Implementación paralela 55 Algoritmo 3 kernel_Hermite [tbp] Input: Muestras de ECG de un latido xi[l],N(número funciones a utilizar) Output: σycn(σ)que minimizan el error para el latido 1: i = bloque.ID # ID del bloque usado como índice de latido 2: l = hilo.ID # ID del hilo usado como índice de muestra 3: Copiar xi[l]a la memoria compartida # Copia completamente paralela 4: Errmin =∞ 5: for all σdo 6: for all ndo 7: Calcular varTempHilol=xi[l]φn[l,σ)#Ecuación (2.8) 8: end for 9: for all ndo 10: Calcular cn(σ) = (∑lvarTempHilol) #Suma paralela por reducción 11: end for 12: ˆxi[l] = 0 13: for all ndo 14: Calcular ˆxi[l]+ = cn(σ)φn[l,σ)#Ecuación (2.2) 15: end for 16: Calcular e[l]+ = (xi[l]-ˆxi[l])2#Ecuación (2.7) 17: Calcular MSE =∑le[l]# Suma paralela por reducción 18: if l== 0 y MSE <Errmin then # Solo entrará en el condicional el hilo 0 19: Sigmabest =σ 20: Cbest =cn(σ) 21: Errmin =MSE 22: end if 23: end for Send subset (0) kernel_Φkernel_Hermite (0) ... time Send subset (1) Read (0) kernel_Hermite (1) Send subset (2) Read (1) Send subset (3) Read (2) kernel_Hermite (2) Figura 3.2: Esquema del flujo de trabajo optimizado de la GPU. Holter). Dichos registros contienen una gran cantidad de latidos que no caben en la memoria de la GPU y que hacen necesario un procesamiento por lotes en conjuntos de latidos más pequeños. Por otra parte, en el caso de procesamiento en tiempo real también es necesario que los latidos sean procesados en pequeños conjuntos para evitar una gran latencia en la presentación 62 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs 10 100 1000 4400 0 10 20 30 40 50 60 70 80 90 100 No of beats (M) % Computation time distribution (TEST A - N=30) Malloc CPU Malloc GPU Kernel_ 1st transfer Kernel_Hermite Read Figura 3.5: Porcentajes de tiempo de ejecución dedicados a cada tarea en Test A con N=30. 3.4. Resultados y Discusión 63 Tabla 3.2: Resultados para el Test B del tiempo de ejecución para cada uno de los experimentos y aceleración obtenida. N M1 M2 CPU (ms) GPU (ms) Aceleración 6 20 5000 211559 1230 172x 10 21 4762 322829 1735 186x 20 22 4546 599249 4612 130x 30 23 4348 845513 9882 86x el tiempo de transferencia para registros cortos tampoco es desdeñable en este caso. Por otra parte, el tiempo necesario para calcular las funciones de Hermite aumenta al haber aumentado N, pero permanece constante independientemente del número de latidos y por tanto pierde importancia en registros largos. Debe resaltarse que el rendimiento no se incrementa siempre de forma monótona (véase Tabla 3.1). La causa de esto es que para valores de N>20 el número de variables locales en el kernel excede el número de registros en el SP y por tanto es necesario recurrir a la memoria global, incurriendo en el llamado register spilling [101]. Esto supone un impedimento importante para conseguir un rendimiento óptimo en la GPU. 3.4.2. Test B: Procesado en diferido (offline) de registros largos El tiempo de procesado necesario y la aceleración lograda en este escenario se muestran en la Tabla 3.2. En la primera columna de la tabla se indica el número de funciones utilizadas para la representación (N). En la segunda y tercera columnas se indican el número de bloques de latidos procesados (M1) y el número de latidos por bloque (M2). El número total de latidos (M) se puede obtener multiplicando ambos datos (M=M1·M2) y es siempre de aproximadamente 105latidos. La cuarta y quinta columnas contienen, respectivamente, el tiempo de la implementación de referencia (CPU) y la implementación paralelizada (GPU). Para finalizar la última columna muestra la aceleración conseguida (tiempo en la CPU dividido tiempo en la GPU). La aceleración de la GPU es sensiblemente mayor que en el primer test. Esto concuerda con lo apreciado en los resultados del Test A, en los que el mayor rendimiento se lograba al procesar registros largos. Se puede observar cómo de nuevo para valores de Nmayores la aceleración es menor, debido al ya mencionado efecto del register spilling. En la Figura 3.6 se muestra la distribución del tiempo de ejecución entre las distintas tareas para distintos 64 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs 6 10 20 30 0 10 20 30 40 50 60 70 80 90 100 Polynomial order (N) % Computation time distribution (TEST B) Malloc CPU Malloc GPU Kernel_ 1st transfer Kernel_Hermite Read Figura 3.6: Porcentajes de tiempo de ejecución dedicados a cada tarea en el Test B para distintos valores de N. 3.4. Resultados y Discusión 65 Tabla 3.3: Resultados para el Test C del tiempo de ejecución para cada uno de los experimentos y aceleración obtenida. N M2 CPU (ms) GPU (ms) Aceleración 6 1 160652 46882 3.43x 5 9315 17.25x 10 5404 29.73x 100 1494 107.53x 10 1 250691 61864 4.05x 5 12815 19.56x 10 6613 37.90x 100 2267 110.56x 20 1 472305 101674 4.65x 5 20626 22.90x 10 10402 45.40x 100 4965 95.12x 30 1 689779 157276 4.39x 5 29978 23.01x 10 15653 44.07x 100 11359 60.72x valores de N. Como se puede apreciar, la tarea asociada al kernel_Hermite ocupa ahora más del 90% del tiempo de ejecución. En este caso no aparece representado el tiempo necesario para transferir los datos a la GPU ya que este proceso se solapa con la computación del kernel_Hermite, siguiendo la estrategia explicada en la Sección 3.3.3 (véase Figura 3.2). La aceleración obtenida varía entre 86x y 172x, demostrando claramente los beneficios de una GPU para el procesado de registros de larga duración. El tiempo de procesado se reduce de minutos (CPU) a segundos (GPU). Para un registro Holter con 12 canales de 3 días, para obtener la representación de Hermite utilizando 30 funciones, pasamos de necesitar más de 3 horas y media con una CPU a necesitar dos minutos y medio. Esto permite liberar recursos computacionales que podrían ser utilizados para aplicar técnicas de análisis sobre la representación obtenida de los latidos. 3.4.3. Test C: Procesado en tiempo real (online) Los resultados correspondientes a este último escenario se muestran en la Tabla 3.3. En la primera columna se muestra el número de funciones de Hermite utilizadas (N). En la segunda el número de latidos procesados en cada bloque (M2). Cuanto más alto sea M2, mayor será 66 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs 1 5 10 100 103 104 105 106 Beats per block (M2) Time (ms) Computation time: CPU vs GPU (Test C) 6 6 10 10 20 20 30 30 GPU CPU Figura 3.7: Tiempo de ejecución para ambas implementaciones en el Test C. la latencia del sistema. Asumiendo una frecuencia cardiaca ficticia de 1 Hz (un latido por segundo), el tiempo de latencia sería igual a M2 segundos. El valor de M1 no aparece por claridad y fue establecido de tal forma que M=M1·M2≈105. La tercera y la cuarta columnas muestran, respectivamente, el tiempo necesario para la ejecución de la implementación de referencia (CPU) y la implementación paralelizada (GPU). El tiempo de computación en la CPU es constante para un valor dado de Nya que no depende de M2 y el número de latidos Mes siempre el mismo. En la última columna se muestra la aceleración obtenida (tiempo en la CPU dividido tiempo en la GPU). 3.4. Resultados y Discusión 67 La aceleración varía entre 4x y 110x, siendo sensiblemente menor cuando se impone un valor de M2 pequeño. Los resultados obtenidos demuestran que el tiempo de ejecución es menor al usar una GPU incluso cuando se impone una latencia mínima (M2=1). No obstante, la aceleración obtenida en este caso no parece suficiente como para justificar el uso de una GPU. Una CPU con una implementación multihilo más optimizada podría llegar a alcanzar un tiempo similar al de la GPU. Si se aumenta el tamaño de M2 es posible ver cómo el beneficio de utilizar una GPU aumenta rápidamente (véase Figura 3.7). Otro aspecto a tener en cuenta es que, aunque la CPU sea capaz de procesar en tiempo real la caracterización de Hermite (el tiempo de computación es menor que el tiempo que tarda el corazón en generar los latidos), se limitaría en gran medida la capacidad de computo disponible para el procesado posterior sobre esta representación. Sin embargo, si se delega el cálculo de la representación de Hermite en la GPU, es posible usar toda la capacidad de cómputo de la CPU para realizar un análisis adicional del latido. En la Figura 3.8 se puede ver cómo se reparte el tiempo de ejecución en la implementación paralela utilizando 6 funciones de Hermite. Nuevamente, la mayor parte del tiempo lo ocupa la tarea de calcular los coeficientes de Hermite (kernel_Hermite). En este caso, a diferencia de los escenarios anteriores, la tarea de transferir los resultados a la memoria del ordenador ocupa un porcentaje de tiempo visible en el gráfico. Dicho porcentaje es mayor para valores pequeños de M2, lo que implica que para conseguir una latencia menor la GPU va a dedicar un porcentaje considerable de tiempo a esta tarea. El reparto porcentual de tiempo para los otros valores de Nes similar al mostrado. Como se puede apreciar en los resultados, hay un claro compromiso entre la latencia y la aceleración conseguida utilizando una GPU. Esto concuerda con los resultados observados en los escenarios previos, en los que la GPU mostraba mejor rendimiento cuando trabajaba con un gran número de latidos. Por tanto, para caracterizar latidos en tiempo real el uso de una GPU puede ser una alternativa viable para el cálculo de la representación de Hermite incluso cuando se un utilizan 30 funciones, suponiendo que se puedan agrupar los latidos en grupos de 5 o más latidos, lo cual implicaría una latencia de unos 5 segundos. 68 Capítulo 3. Cálculo de la representación de Hermite utilizando GPUs 1 5 10 100 0 10 20 30 40 50 60 70 80 90 100 Block size (M2) % Computation time distribution (TEST C - N=6) Malloc CPU Malloc GPU Kernel_ 1st transfer Kernel_Hermite Read Figura 3.8: Porcentajes de tiempo de ejecución dedicados a cada tarea en el Test C para N=6. CAPÍTULO 4 AGRUPAMIENTO DE LATIDOS MEDIANTE ACUMULACIÓN DE EVIDENCIA Una vez resuelto el problema de la representación del latido abordaremos el problema de su agrupamiento morfológico. Una parte importante del proceso de análisis basado en agrupamiento consiste en seleccionar la técnica apropiada para el conjunto de datos que estamos explorando. Para conseguirlo, generalmente se realizan varias pruebas con diferentes algoritmos y distintas configuraciones, hasta encontrar una que proporcione un resultado satisfactorio, basándose en la información disponible del problema. Sin embargo, este es un proceso lento, típicamente guiado por heurísticas, con una fuerte dependencia del criterio del analista y proclive a error [39]. En la bibliografía de agrupamiento de latidos no existe un consenso con respecto a la superioridad de una técnica de agrupamiento frente a otra, pudiendo obtenerse resultados similares con técnicas muy diversas [20][21][26][27][75][77][93][120] (véase Sección 1.2). Esto resulta coherente con los análisis comparativos en el ámbito del agrupamiento automático [65]; ninguna técnica de agrupamiento ha mostrado hasta ahora una superioridad manifiesta sobre las demás, sino que simplemente presenta una mejor adaptación a las características de algunos tipos específicos de problemas. Tomando inspiración de los trabajos de combinación de clasificadores y de fusión de sensores, en los últimos años se han desarrollado múltiples técnicas de agrupamiento mediante la combinación de particiones obtenidas por distintos algoritmos de agrupamiento [11][36][42][61][138]. Estas nuevas técnicas de agrupamiento por medio de combinación, 70 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia ●●●●●●●●●●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●●●●●● ● ● ● ● ● 5 10 15 20 25 30 5 10 15 20 25 30 V1 a) Datos originales ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● ● ● ● ● ● ● 5 10 15 20 25 30 5 10 15 20 25 30 V1 V2 b) K−means K=3 ● ● ● ●● ●● ●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●●●●●● ●●●●●●● ● ● ● ● ● ● ● ● ● ● 5 10 15 20 25 30 5 10 15 20 25 30 c) K−means K=20 ●●●●●●●●●●●●●●●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ●●●●●● ● ● ● ● ● 5 10 15 20 25 30 5 10 15 20 25 30 V2 d) Ensemble de 100 particiones K=[7−37] Figura 4.1: (a) muestra las tres particiones naturales de los datos; (b) y (c) muestran dos versiones del algoritmo K-means con K=3 y K=20, respectivamente; (d) muestra el resultado de la combinación de 100 particiones de datos distintas creadas por el algoritmo K-means con valores de K aleatorios entre 7 y 37. 71 también llamadas agrupamiento mediante ensembles, se han afianzado como alternativas a las técnicas de agrupamiento tradicionales, ya que a menudo son capaces de mejorar la robustez y la estabilidad de los resultados del agrupamiento. Mediante ellas es posible combinar las particiones generadas por distintos algoritmos conservando en gran parte sus ventajas individuales y contrarrestando sus desventajas. Podemos ver un ejemplo en la Figura 4.1 donde se muestra cómo combinando los resultados de 100 ejecuciones del algoritmo K-means es posible representar las tres particiones originales de los datos, particiones que nunca podrían haber sido obtenidas por una sola aplicación de K-means, dado que este algoritmo siempre obtiene particiones hiperesféricas. Otra ventaja adicional del agrupamiento mediante ensembles es que permite reducir la dependencia entre el algoritmo de agrupamiento y los resultados, ya que posibilita la utilización de múltiples algoritmos con diferentes parámetros de inicialización, combinando posteriormente los resultados obtenidos por cada uno de ellos. El paradigma de acumulación de evidencia (Evidence Acumulation Clustering, EAC) proporciona una de las técnicas que permite realizar agrupamiento mediante ensembles. Propuesto por primera vez en [42], es el resultado de una búsqueda de nuevas técnicas que no impongan una forma o modelo en los grupos. Para ello combina varios resultados de agrupamiento, extrayendo la información individual de cada partición y combinándola para definir los grupos. Esta misma idea ya ha sido extensamente explorada en clasificación, resultando en mejoras en la exactitud y precisión por medio de la combinación de múltiples clasificadores [58][80][125]. En EAC se combinan los resultados individuales de varios agrupamientos para obtener una nueva medida de similitud entre las instancias, que integra y resume la información obtenida por cada uno de los algoritmos utilizados para obtener las particiones. Posteriormente, se obtendrá la partición final de los datos basándose en esta nueva medida de similitud. En la siguiente sección se expondrá una nueva técnica de agrupamiento basada en acumulación de evidencia. Posteriormente, se aplicará dicha técnica al agrupamiento del latido según su origen cardiaco, utilizando para ello la base de datos MIT-BIH Arrhythmia Database. Finalmente se extenderá el método para usarlo con 12 derivaciones y se validará sobre la base de datos INCARTDB. 78 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia 4.1.4. Análisis de complejidad computacional de PN-EAC Para estudiar el rendimiento de PN-EAC se realizó un análisis de la complejidad de cada una de sus partes. La complejidad del algoritmo K-means utilizado para crear las particiones (véase Algoritmo 4, línea 2) ha sido ampliamente estudiada en la bibliografía [52]: O(n·k·d· i), siendo iel número de iteraciones, del número de dimensiones del vector de características ykel número de grupos. El número de dimensiones para la mayor parte de los problemas será fijo y mucho menor que n. Lo mismo ocurre con el número de iteraciones i, que estará limitado por un máximo fijo (típicamente 100). El número de grupos ken nuestro algoritmo se determina según (4.1) por lo que como mucho será √n. Por lo tanto, la complejidad de K-means puede aproximarse por: O(n·k·d·i)≈O(n·√n). Este algoritmo será ejecutado una vez por cada partición (véase Algoritmo 4, líneas 1-3) por lo que la complejidad total de este paso del algoritmo sería O(m·n·√n), siendo mel número de particiones creadas. La extracción de la evidencia de las particiones requiere recorrer la matriz de evidencia, por lo que la complejidad está determinada por el tamaño de la matriz (n×n), siendo O(n2). Dicha matriz se recorre una vez por cada partición generada, por lo que la complejidad computacional de recopilar la evidencia es O(m·n2)(véase Algoritmo 4, líneas 4-6). Finalmente, debemos analizar la complejidad del algoritmo jerárquico utilizado para extraer la partición final. En dichos algoritmos se comienza por calcular la matriz de distancias entre todos los objetos, operación que tiene una complejidad de O(n2). Posteriormente, se realizan niteraciones y en cada una de ellas se encuentran los dos grupos más cercanos, que se unen, y se actualiza la distancia del resto de grupos con el nuevo grupo. La complejidad de cada una de estas iteraciones es de O(n·log(n)), por lo que al ser n iteraciones la complejidad total es O(n2·log(n)), siendo este el término que determina la complejidad del algoritmo jerárquico [28] (véase Algoritmo 4, línea 7). La complejidad del algoritmo jerárquico (O(n2·log(n))) podría llegar a ser superior a la de extraer la evidencia (O(m·n2)), en el caso de que log(n)sea mayor que m. Sin embargo, para valores razonables de nymse cumple que m>> log(n). Por ejemplo, para 300 particiones n tendría que ser mayor que 2300. Por ello, la complejidad de PN-EAC será O(m·n2): O(m·n·√n)+ O(m·n2)+O(n2·log(n)) ≈O(m·n2).(4.5) En cuanto a la complejidad espacial, vendrá dada por el tamaño de la matriz de evidencia y será O(n2), que en una implementación optimizada podría reducirse a O(n2/2), al ser esta matriz simétrica. 4.2. Agrupamiento de latidos con PN-EAC 79 4.2. Agrupamiento de latidos con PN-EAC Una de las mayores dificultades en el agrupamiento de latidos reside en la complejidad morfológica del propio latido, de modo que la variabilidad intrínseca a cualquier familia morfológica, es decir, al conjunto de latidos que comparten un mismo origen de activación y camino de propagación, no se proyecta en un conjunto de características que muestren propiedades de simetría en el espacio de representación. Así es que, sea cual sea la fórmula adoptada para representar el latido, se ha comprobado que las distintas familias morfológicas muestran un conjunto de formas heterogéneas en el espacio de representación. La técnica propuesta en la sección anterior (PN-EAC) propone un agrupamiento en el que la forma de los grupos se define mediante la combinación de resultados individuales, permitiendo para cada grupo una forma diferente y arbitraria. Además, dicha técnica permite combinar de forma natural información de distintas fuentes, lo que permitiría combinar fácilmente información de distintas derivaciones. Para el agrupamiento se utilizará la representación de Hermite del latido, ya presentada en el Capítulo 2. Esta representación solamente captura información sobre el complejo QRS. En ciertos tipos de arritmias, como contracciones prematuras auriculares y atrioventriculares, y latidos de escape auriculares y de unión, la morfología del complejo QRS es similar, y no es suficiente para discriminar y agrupar correctamente los latidos (véanse Figuras 1.6 y 4.3). En estos casos es necesaria información adicional, que idealmente se obtiene mediante la identificación de la onda P. Sin embargo, la dificultad de delinear el latido y extraer la onda P con precisión hace que sea habitual recurrir a obtener dicha información por otras vías [20][31], como en nuestro caso, que se obtendrá mediante información relativa a la distancia entre latidos. Por lo tanto, para permitir al algoritmo distinguir entre ciertos tipos de latidos cuyos complejos QRS tienen una forma muy similar se incluyeron dos características basadas en la distancia entre latidos: R1[i] = R[i]−R[i−1],(4.6) R2[i] = u(α)·α, α= (R1[i+1]−R1[i])−(R1[i]−R1[i−1]),(4.7) 80 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia MLII Figura 4.3: Fragmento de ECG en el que al principio aparecen latidos prematuros atrioventriculares (J) y posteriormente, aproximadamente en la mitad del fragmento, las distancias entre latidos se incrementan y aparecen latidos normales. (Fuente: MIT-BIH Arrhythmia Database Arrhythmia Database, registro 234, entre 0:14:27 y 0:14:37) donde R[i]es el momento de ocurrencia del latido número i, dado por las anotaciones de la base de datos, y u(x)es la función escalón Heaviside: u(x) =    1 si x≥0 0 si x<0(4.8) Esta función permite eliminar los valores negativos poniéndolos a cero y conservar los positivos en la ecuación (4.7). La expresión (4.6) mide la distancia de un latido con el latido previo. Comparando los valores de esta medida es posible identificar si un latido es prematuro. Para ello, la expresión (4.7) compara esta distancia con otras dos distancias: la distancia del latido actual con el siguiente y la distancia del latido previo con el anterior a este. De esta forma se normaliza esta medida para que no dependa del ritmo cardíaco. En consecuencia, un valor alto para (4.7) proporciona evidencia de que el latido puede ser prematuro. Por tanto, nuestro latido estará representado por los coeficientes de las funciones de Hermite para cada una de las derivaciones, los parámetros σy las características dadas por (4.6) y (4.7). Se propone un conjunto de experimentos para aplicar el algoritmo PN-EAC al agrupamiento de latidos. Para estos experimentos se limitará el número de funciones de Hermite utilizadas en la representación a 16. Se utilizará la base de datos MIT-BIH Arrhythmia Database completa (véase Sección 1.3.1), eliminando tanto el ruido de alta frecuencia como la deriva de línea base, mediante los filtros ya presentados en la Sección 2.1.1. Sin embargo, para facilitar la comparación con otros trabajos de la bibliografía, no se aplicará corrección alguna sobre las posiciones de los latidos de la base de datos, utilizando por tanto las posiciones originales. 4.2. Agrupamiento de latidos con PN-EAC 81 Cada latido será representado por 36 características, 16 coeficientes de Hermite por cada derivación de la base de datos, un valor de σpor cada derivación y las dos medidas basadas en la distancia entre latidos, dadas por (4.6) y (4.7). 4.2.1. Estrategias para la generación de particiones A la hora de aplicar el paradigma de acumulación de evidencia al agrupamiento de latidos se emplearán tres estrategias distintas para generar los ensembles y agrupar los latidos. Estrategia 1 La primera estrategia es la aproximación clásica de aprendizaje automático en la que todas las características disponibles de una instancia conforman un único vector. Por tanto, en el vector aparecerán los parámetros de Hermite (coeficientes y valor de σ) para las dos derivaciones y las características dadas por las ecuaciones (4.6) y (4.7): ui= (C1,σ1,C2,σ2,R1,R2),(4.9) donde uies el vector de características que representa el latido xi,Cdson los coeficientes de Hermite que representan el complejo QRS en la derivación d, Cd={c0(σd),c1(σd),...,cN−1(σd)}(en este caso con N=16) y σdes el valor de σ en dicha derivación. Utilizando la representación indicada en la ecuación (4.9) para cada latido se generarán las particiones de datos P1,P2,...,Pm(véase Algoritmo 5, líneas 1-3). Estrategia 2 La segunda estrategia se basa en la hipótesis de que generando particiones con distintas representaciones de datos por separado es posible obtener mejores resultados que agrupando todas las características en el mismo vector. Como ya se vio en la Sección 1.1.3, los cardiólogos suelen apoyarse para el diagnóstico en la interpretación de varias derivaciones de forma paralela. La segunda estrategia imita este modus operandi extrayendo evidencia por separado de la información procedente de cada derivación. También se separa la información de las características dadas por (4.6) y (4.7) al representar una información distinta a la morfológica, con otro marco de referencia. Siguiendo esta estrategia, se dividirán las características en tres representaciones diferentes del latido: una representación para la información morfológica extraída de cada derivación, representada mediante las características de Hermite, y una tercera representación con la información de las características 82 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia Algoritmo 5 Agrupamiento de latidos con PN-EAC (Estrategia 1) Input: Conjunto de latidos X={x1,x2,...,xn} Input: mnúmero de particiones a generar Output: P∗partición final de los nelementos 1: for all latido xido 2: ui←(representación para el latido xisegún la Ecuación (4.9)) 3: end for 4: U={u1,u2,...,un} 5: for j=1 to mdo #Sección 4.1.1 6: k=aleatorio([√n/2,√n]) #Ecuación (4.1) 7: Pj=K−means(U,k) 8: end for 9: G=zeros(n,n)#Inicializar matriz evidencia 10: for j=1 to mdo #Sección 4.1.2 11: for h=1 to ndo 12: for f=1 to ndo 13: if uhyufestán en un mismo grupo en Pjthen 14: G(h,f)=G(h,f)+1 #Ecuación (4.2) 15: end if 16: end for 17: end for 18: end for 19: G∗=G/m#Ecuación (4.2) 20: P∗=enlace.medio(G∗)#Sección 4.1.3 derivadas de la distancia entre los latidos. Por consiguiente, tendremos tres vectores de características: v1 i= (C1,σ1),(4.10) v2 i= (C2,σ2),(4.11) v3 i= (R1,R2),(4.12) dónde v1 i,v2 iyv3 ison los vectores de características de cada una de las tres representaciones del latido xi. Dichas representaciones se utilizarán por separado para generar particiones de las que extraer evidencia positiva y usando esta evidencia se extraerá la partición final (véase Algoritmo 6, líneas 1-5). Estrategia 3 La tercera estrategia es una variante de la segunda, pero en esta ocasión se aplicará el algoritmo PN-EAC, extrayendo evidencia negativa de las características 4.2. Agrupamiento de latidos con PN-EAC 83 Algoritmo 6 Agrupamiento de latidos con PN-EAC (Estrategia 2) Input: Conjunto de latidos X={x1,x2,...,xn} Input: mnúmero de particiones a generar Output: P∗partición final de los nelementos 1: for all latido xido 2: v1 i←(representación para el latido xisegún la Ecuación (4.10)) 3: v2 i←(representación para el latido xisegún la Ecuación (4.11)) 4: v3 i←(representación para el latido xisegún la Ecuación (4.12)) 5: end for 6: V1={v1 1,v1 2,...,v1 n} 7: V2={v2 1,v2 2,...,v2 n} 8: V3={v3 1,v3 2,...,v3 n} 9: for j=1 to mdo #Sección 4.1.1 10: P1 j=K−means(V1,k=aleatorio([√n/2,√n])) 11: P2 j=K−means(V2,k=aleatorio([√n/2,√n])) 12: P3 j=K−means(V3,k=aleatorio([√n/2,√n])) 13: end for 14: G=zeros(n,n)#Inicializar matriz evidencia 15: for j=1 to mdo #Sección 4.1.2 16: for h=1 to ndo 17: for f=1 to ndo 18: if v1 hyv1 festán en un mismo grupo en P1 jthen 19: G(h,f)=G(h,f)+1 # Ecuación 4.2 20: end if 21: if v2 hyv2 festán en un mismo grupo en P2 jthen 22: G(h,f)=G(h,f)+1 # Ecuación 4.2 23: end if 24: if v3 hyv3 festán en un mismo grupo en P3 jthen 25: G(h,f)=G(h,f)+1 # Ecuación 4.2 26: end if 27: end for 28: end for 29: end for 30: G∗=G/(3·m)# Ecuación 4.2 31: P∗=enlace.medio(G∗)#Sección 4.1.3 84 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia N N N R N R R R ML II Figura 4.4: Se muestra un fragmento de ECG con cuatro latidos normales seguidos de cuatro latidos con bloqueo de rama derecha. Resaltar que la distancia entre latidos es similar para los 8, teniendo los latidos patológicos aproximadamente los mismos valores de las características (4.6) y (4.7) que los latidos normales. (Fuente: MIT-BIH Arrhythmia Database, registro 212, entre 0:12:13 y 0:12:18) derivadas de los instantes de ocurrencia de los latidos (v3 i). El hecho de que los valores de (4.6) y (4.7) no cambien, es decir, que no haya cambios en la distancia entre latidos, no implica que no haya un cambio en el tipo de latido; existe una gran variedad de tipos de latido que pueden darse con distancias iguales entre latidos. Por lo tanto, que los valores de estas características sean similares no proporciona información concluyente de si dos latidos pertenecen al mismo tipo (véase Figura 4.4). Sin embargo, dos latidos con diferencias considerables en estas características generalmente pertenecen a distintos tipos de latido (véase Figura 4.5). En esta estrategia para modelar esta información obtenida de las ecuaciones (4.6) y (4.7) se utilizará la evidencia negativa. Se emplearán las mismas tres representaciones del latido cardiaco que en la estrategia anterior (véanse Ecuaciones (4.10), (4.11) y (4.12)), pero en este caso de la representación derivada de la distancia entre latidos (4.12) se extraerá evidencia negativa, mientras que de las otras dos se continuará extrayendo evidencia positiva (véase Algoritmo 7, líneas 1-5). Para la generación de las particiones se utilizará el algoritmo K-means con diferentes inicializaciones y un número de grupos aleatorio en el rango dado por (4.1) (véase Sección 4.1.1). En la primera estrategia, donde todas las características del latido están en un mismo vector ui, se generan 300 particiones de datos. De forma similar, para la segunda y tercera estrategias se generan 100 particiones por cada vector de características (v1 i,v2 iyv3 i), con lo que tenemos en total 300 particiones por cada estrategia. En la segunda estrategia se extrae evidencia positiva de las tres representaciones dadas por las ecuaciones (4.10), (4.11) y (4.12). 4.2. Agrupamiento de latidos con PN-EAC 85 Algoritmo 7 Agrupamiento de latidos con PN-EAC (Estrategia 3) Input: Conjunto de latidos X={x1,x2,...,xn} Input: mnúmero de particiones a generar Output: P∗partición final de los nelementos 1: for all latido xido 2: v1 i←(representación para el latido xisegún la Ecuación (4.10)) 3: v2 i←(representación para el latido xisegún la Ecuación (4.11)) 4: v3 i←(representación para el latido xisegún la Ecuación (4.12)) 5: end for 6: V1={v1 1,v1 2,...,v1 n} 7: V2={v2 1,v2 2,...,v2 n} 8: V3={v3 1,v3 2,...,v3 n} 9: for j=1 to mdo #Sección 4.1.1 10: P1 j=K−means(V1,k=aleatorio([√n/2,√n])) 11: P2 j=K−means(V2,k=aleatorio([√n/2,√n])) 12: P3 j=K−means(V3,k=aleatorio([√n/2,√n])) 13: end for 14: G=zeros(n,n)#Inicializar matriz evidencia 15: G−=zeros(n,n)#Inicializar matriz evidencia negativa 16: for j=1 to mdo #Sección 4.1.2 17: for h=1 to ndo 18: for f=1 to ndo 19: if v1 hyv1 festán en un mismo grupo en P1 jthen 20: G(h,f)=G(h,f)+1 #Ecuación (4.2) 21: end if 22: if v2 hyv2 festán en un mismo grupo en P2 jthen 23: G(h,f)=G(h,f)+1 #Ecuación (4.2) 24: end if 25: if v3 hyv3 fno están en un mismo grupo en P3 jthen 26: G− (h,f)=G(h,f)−1 #Ecuación (4.3) 27: end if 28: end for 29: end for 30: end for 31: G=G/(2·m)#Ecuación (4.2) 32: G−=G−/m#Ecuación (4.3) 33: G∗=G+G−#Ecuación (4.4) 34: P∗=enlace.medio(G∗)#Sección 4.1.3 86 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia N N A N N N N ML II Figura 4.5: En este fragmento de ECG el tercer latido es un latido auricular prematuro, mientras que el resto son latidos normales. Se puede apreciar cómo el latido prematuro es morfológicamente similar a los normales, pero la distancia con el latido anterior y el siguiente cambia. (Fuente: MIT-BIH Arrhythmia Database, registro 223, entre 0:01:02 y 0:01:07) En la tercera estrategia, de la representación derivada de la distancia entre latidos (4.12) se extrae únicamente evidencia negativa (Algoritmo 7, líneas 25-27). En el caso de esta última estrategia, para obtener la matriz final de evidencia se combina la evidencia positiva y la negativa mediante la ecuación (4.4) (Algoritmo 7, línea 33). La partición final de los datos P∗se extrae de la matriz de evidencia aplicando el algoritmo de enlace medio (Algoritmo 5, línea 20; Algoritmo 6, línea 31; Algoritmo 7, línea 34). En nuestro caso, para el agrupamiento de los latidos, esta variante es la que ha mostrado el mejor rendimiento. El número de grupos de esta partición final se determina utilizando el criterio del tiempo de vida. También se obtiene el resultado con un número fijo de 25 grupos por registro. 4.2.2. Resultados Para evaluar el rendimiento del agrupamiento se considerará que cada grupo obtenido pertenece al tipo de latido mayoritario dentro del grupo, de acuerdo con la anotación de la base de datos, considerando por tanto a todos los latidos distintos del tipo mayoritario del grupo como errores. En la rutina clínica, si fuera necesario, la correspondencia entre grupos y tipos de latidos podría obtenerse de un cardiólogo que anotara un latido de cada grupo. En la Tabla 4.1 se muestran los resultados para las tres estrategias, registro a registro. Para cada registro aparece el número de errores con 25 grupos por registro (25C) y el número de grupos y de errores cuando se utiliza el criterio del tiempo de vida (Tiempo de vida). Al final de la tabla aparece el número total de errores y el porcentaje que representan sobre el total de la base de datos (109966 latidos). Al utilizar 25 grupos por registro se obtienen unos porcentajes de error de 2.25%, 1.81% y 1.44% para la primera, segunda y tercera estrategias, 4.2. Agrupamiento de latidos con PN-EAC 87 respectivamente. Al utilizar el criterio del tiempo de vida los porcentajes son de 2.75%, 3.81% y 0.99% con un total de 508, 203 y 6947 grupos, respectivamente. Tabla 4.1. Resultados del agrupamiento para las tres estrategias. Se muestra, para cada estrategia, para cada registro (#), el número de errores (Err) utilizando un número fijo de 25 grupos (25C) y utilizando el criterio del tiempo de vida. En este último caso se muestra también el número de grupos seleccionado por dicho criterio. 1aEstrategia 2aEstrategia 3aEstrategia #25C Tiempo de vida 25C Tiempo de vida 25C Tiempo de vida Err Err Grupos Err Err Grupos Err Err Grupos 100 33 33 5 6 33 2 9 33 5 101 3 3 7 0 3 5 1 3 5 102 7 47 6 13 58 4 28 28 20 103 1 1 15 0 1 2 0 0 81 104 251 257 11 309 351 4 110 106 43 105 11 12 10 5 5 5 5 5 52 106 2 10 9 1 28 7 0 0 85 107 0 1 9 1 1 3 0 0 8 108 11 16 9 9 9 2 6 6 66 109 4 9 14 2 10 2 2 2 39 111 0 0 11 0 0 2 0 0 22 112 2 2 7 1 2 3 1 1 16 113 0 0 11 0 0 2 0 0 49 114 12 16 12 11 16 4 13 6 61 115 0 0 4 0 0 3 0 0 5 116 2 2 13 0 2 2 1 1 11 117 1 1 7 0 0 7 1 1 15 118 96 96 12 58 100 2 34 12 141 119 0 0 7 0 0 2 0 0 33 121 1 1 12 0 1 2 1 1 5 122 0 0 8 0 0 8 0 0 4 123 0 0 6 0 0 2 0 0 16 124 36 43 9 41 41 4 38 36 44 200 129 130 15 117 531 14 84 84 22 94 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia of Cardiological Technics 12-lead Arrhythmia Database (INCARTDB) [103] (véase Sección 1.3.2). Esta base de datos también goza de cierta popularidad en la comunidad científica. Sin embargo, incluso aquellos autores que la utilizan se limitan en su mayoría a seleccionar una o dos derivaciones. En [7] se realiza una clasificación de latidos utilizando árboles de decisión; para ello se utiliza la base de datos MIT-BIH Arrhythmia Database junto con 4000 latidos de la INCARTDB, pero usando una única derivación. En [94] se presenta un clasificador de latidos en dos etapas con un agrupamiento mediante el algoritmo Fuzzy C-means y clasificación utilizando redes neuronales. El trabajo se validó sobre la base de datos INCARTDB, entre otras, pero empleando solamente una derivación. Otro clasificador, en este caso basado en ELM (Extreme Learning Machines) que utiliza una derivación de esta base de datos es [6]. En [82] se utiliza la INCARTDB para evaluar la detección de latidos ventriculares prematuros, pero nuevamente utilizando una sola derivación. En [85] se utilizan dos derivaciones de esta base de datos para desarrollador un clasificador lineal de latidos. El mismo clasificador lineal desarrollado en [85] es utilizado en [83], pero en este caso es mejorado para ser capaz de utilizar las 12 derivaciones. El clasificador se aplica a todas las derivaciones, a un subconjunto de ellas y a la mejor derivación. Finalmente, en [84] se amplía el trabajo publicado en [83] realizando la validación sobre varias bases de datos, utilizando distintas representaciones del latido y con diferentes estrategias para combinar la información de las derivaciones. El principal motivo por el que no se suelen usar más de un par de derivaciones es el incremento en la dimensionalidad del espacio de características empleado en la representación del latido. El uso de 12 derivaciones de modo concurrente significa multiplicar por 12 la dimensionalidad del espacio de características comparado con el uso de una única derivación. Este incremento en un orden de magnitud supone un reto en la aplicación de técnicas de aprendizaje automático. Sin embargo, empleando el paradigma de acumulación de evidencia cada derivación puede emplearse como una fuente complementaria de evidencia a combinar para la creación de la partición final, lo que elimina el problema del incremento en la dimensión del espacio de características y a la vez permite explotar y combinar toda la información de las 12 derivaciones. En esta sección se empleará PN-EAC para crear un algoritmo de agrupamiento de latidos que trabaje de modo simultáneo sobre 12 derivaciones. También se estudiará cómo evolucionan los resultados del agrupamiento según se va incorporando información proveniente de un mayor número de derivaciones. Para ello probaremos a ejecutar el 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC 95 algoritmo sobre 2, 4, 6, 8, 10 y 12 derivaciones y analizaremos los resultados para comprobar si existe una mejoría al incrementar el número de derivaciones. En las pruebas realizadas en este capítulo se utilizaron 71 registros de los 75 de la base de datos INCARTDB (véase Sección 1.3.2). Los registros excluidos fueron el I02, I03, I57 y I58 debido a que en ellos una de las derivaciones está ausente, por lo que se prefirió excluirlos y trabajar solo sobre los registros que contienen las 12 derivaciones. 4.3.1. Estrategias para la generación de particiones Se obtendrá para cada una de las derivaciones la representación de Hermite (coeficientes yσ), utilizando de nuevo 16 funciones. Posteriormente se calcularán características derivadas de la distancia entre latidos dadas por las ecuaciones (4.6) y (4.7) calculadas sobre las anotaciones de la base de datos. En consecuencia, cada latido será caracterizado por doce representaciones de Hermite (16 coeficientes y σ), una por cada derivación, y las dos características derivadas de la distancia entre latidos, en total 206 características por latido. Esta dimensión del espacio de características puede resultar demasiado elevada para la mayoría de los algoritmos de aprendizaje automático, requiriendo generalmente el uso de técnicas de selección de características o de técnicas de reducción de dimensiones [113]. Anteriormente se propusieron tres estrategias para aplicar PN-EAC sobre registros de dos derivaciones (véase Sección 4.2.1). De entre ellas la estrategia basada en obtener evidencia positiva de cada representación de Hermite del latido y evidencia negativa de las características dadas por las ecuaciones (4.6) y (4.7) (Estrategia 3), fue la que obtuvo un mejor resultado (véase Tabla 4.1). Por tanto, nos basamos en esta estrategia para realizar la acumulación de evidencia. Se realizaron pruebas utilizando 2, 4, 6, 8, 10 y 12 derivaciones, extrayendo siempre evidencia de cada derivación de modo independiente. Cada latido xi, con dderivaciones, se representará de la siguiente forma: v1 i= (C1,σ1), v2 i= (C2,σ2), ... vd i= (Cd,σd), vd+1 i= (R1,R2). (4.13) 96 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia Al aplicar el algoritmo sobre 2 derivaciones el procedimiento es análogo al ya presentado en la sección anterior (véase Algoritmo 7); para cada una de las dos derivaciones se generan 100 particiones de las que se extraerá evidencia positiva. Estas particiones, al igual que el resto, se generan mediante el algoritmo K-means. En este caso debemos elegir 2 derivaciones de las 12 disponibles para generar las particiones. Esta elección se realiza de modo aleatorio. Se generan también 100 particiones con las características dadas por (4.6) y (4.7), de las que se extrae evidencia negativa. Todo el proceso se repite 100 veces para obtener una medida del error promedio de usar 2 derivaciones. Para ejecutar PN-EAC con 4 derivaciones se procederá de forma similar. El primer paso es elegir aleatoriamente 4 de las 12 derivaciones disponibles. Para cada una de estas derivaciones se generan 100 particiones con la información de Hermite (coeficientes y σ). Por consiguiente, habrá en total 400 particiones de las que extraer evidencia positiva. Para mantener el mismo peso relativo de evidencia positiva con respecto a la evidencia negativa que el empleado anteriormente (1/2), se generan 200 particiones con las características derivadas de la distancia entre latidos ((4.6) y (4.7)). De este modo se evita que la información procedente de la distancia entre latidos se vaya diluyendo al incrementar el número de derivaciones empleadas para el agrupamiento. Evitar esta dilución es importante ya que la información derivada de las características de Hermite y la información derivada de la distancia entre latidos es complementaria. Para los casos de 6, 8, 10 y 12 derivaciones se procederá de forma análoga, obteniendo 100 particiones por cada derivación para generar evidencia positiva y ajustando el número de particiones generadas con las características dadas por las ecuaciones (4.6) y (4.7) para que el número de particiones que proporcionan evidencia negativa sea siempre 1/2 del número total de particiones que proporcionan evidencia positiva. Para cada número de particiones repetiremos el proceso 100 veces. En cada ocasión las derivaciones utilizadas se eligen aleatoriamente entre las 12 disponibles; excepto en el caso de 12 derivaciones, donde se utilizan todas. La generación de particiones y la extracción de la partición final de la matriz de acumulación de evidencia se realizarán empleando el Algoritmo 7, modificando el número de derivaciones y el peso relativo de la evidencia negativa. En este caso no utilizaremos el criterio del tiempo de vida para seleccionar el número de grupos debido al mal rendimiento que mostró, en especial en lo relativo a la proliferación de grupos; los resultados se obtendrán siempre con 25 grupos. 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC 97 4.3.2. Resultados La Tabla 4.5 muestra el promedio de errores para los distintos números de derivaciones. En todos los casos el número de errores se obtiene utilizando un número fijo de 25 grupos por registro. El porcentaje de error se obtiene de dividir el número total de errores por el número total de latidos de la base de datos (165514 latidos, unos 2300 por registro). Los resultados medios de este porcentaje son 0.6006%, 0.4052%, 0.381%, 0.3579%, 0.3493% y 0.3379% para 2, 4, 6, 8, 10 y 12 derivaciones, respectivamente. Tabla 4.5. Resultados del agrupamiento para 2 (d2), 4 (d4), 6 (d6), 8 (d8), 10 (d10) y 12 (d12) derivaciones. Se muestra para cada registro la media del número de errores en las 100 ejecuciones realizadas. Al final se muestra el promedio del número total de errores y el tanto por ciento de error sobre todos los latidos de la base de datos. d2 d4 d6 d8 d10 d12 I01 0.75 0.07 0 0 0 0 I04 33.79 31.87 30.42 28.58 28.42 27.89 I05 15.73 13.66 13.24 12.73 11.38 10.87 I06 23.35 13.84 10.05 9.44 9.47 9 I07 19.17 5.22 3.72 3.6 3.33 3.07 I08 6.35 3.64 3.76 3.13 3 2.37 I09 9.76 9.51 9.16 9.15 8.67 8 I10 0.85 0.2 0.01 0.01 0 0 I11 24.02 24 24 24 24 24 I12 3.97 2.99 3.11 2.96 3.05 3.15 I13000000 I14000000 I15 0.04 0 0 0 0 0 I16 0.04 0 0 0 0 0 I17 0.19 0.01 0 0 0 0 I18 57.08 52.74 51.52 50.54 50.37 51.7 I19 1.77 1.5 1.43 1.35 1.51 1.82 I20 102.1 60.59 52.78 47.6 47 39.87 I21 51.95 26.31 25.4 24.29 22.59 22.28 I22 63.79 41.17 37.23 36.81 35.93 37.31 I23 0.03 0 0 0 0 0 98 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia d2 d4 d6 d8 d10 d12 I24 0.08 0 0 0 0 0 I25 1.57 1.05 1.14 1 1 1 I26 3.8 2.22 2.22 1.96 1.98 2 I27 0.01 0 0 0 0 0 I28000000 I29 10.55 2.09 1.23 0.61 0.35 0.12 I30 1.9 0.41 0.05 0 0 0 I31 41.17 35.25 45.28 33.29 29.44 19.85 I32 0.06 0 0 0 0 0 I33 90.59 26.28 16.7 15.36 15.07 15.3 I34 63.73 15.93 11.89 11.51 11.85 11.93 I35 59.05 45.68 41.64 37.67 34.94 34.96 I36 47.48 38.37 34.69 36.3 36.67 42.15 I37 1.02 1 1 1 1 1 I38 2.37 0.87 0.64 0.73 0.57 0.33 I39 1.68 0.15 0.03 0.01 0 0 I40 4.26 3.42 3.46 3.23 3.18 3.17 I41 4.96 4.95 4.94 4.82 4.93 5 I42 5.98 5.64 5.19 5.33 5.51 5.18 I43 7.05 5.96 5.07 5.01 5 5 I44 7.04 5.75 5.68 5.48 4.82 4.17 I45 0.02 0 0 0 0 0 I46 7.88 6.02 5.75 5.39 4.97 4.97 I47 2.86 1.12 1.07 1.03 1.02 1 I48 3.15 2.53 2.43 2.66 2.78 3 I49000000 I50 1.91 1.85 1.72 1.82 1.83 1.12 I51 3.34 3.09 2.88 3.01 2.93 3 I52000000 I53 0.1 0.03 0 0 0 0 I54 2.85 0.31 0.16 0.11 0.02 0 I55 1.47 0.45 0.09 0.04 0 0 I56 1.71 0.06 0.02 0 0 0 I59 38.17 37.62 37.83 34.87 33.96 33.05 I60 50.23 48.86 48.58 47.21 45.28 42.41 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC 99 d2 d4 d6 d8 d10 d12 I61 0.88 0.89 0.89 0.93 0.98 1 I62 36.72 26.77 25.6 24.22 23.49 21.92 I63 12.42 10.68 9.72 8.74 9.4 9.48 I64 2.43 1.94 2.2 1.93 1.96 2 I65 11.6 9.47 8.84 8.84 8.94 8.98 I66 1.81 1.77 1.91 1.85 1.81 2 I67 6.15 5.05 5 5 5 5 I68 1.37 1.01 1 1 1 1 I69 0.71 0.37 0.66 0.58 0.76 0.83 I70 0.01 0 0 0 0 0 I71 9.08 6.49 5.97 5.85 5.72 5.93 I72 7.77 7.46 7.56 7.54 7.86 8 I73 3.93 2.02 2.51 2.02 2 2 I74 15.11 12.37 11.66 10.35 11.46 10.04 I75 1.42 0.17 0.03 0 0 0 Total 994.2 670.7 630.8 592.5 578.2 559.2 %0.601 0.405 0.381 0.358 0.349 0.338 Sobre los resultados se aplicaron test estadísticos para verificar la significancia de las diferencias entre utilizar un mayor o menor número de derivaciones [119]. Se rechazó la hipótesis de normalidad de los datos por lo que se aplicaron test no paramétricos. Primeramente, se aplicó el test de Friedman [46]. Dicho test toma como hipótesis nula que los resultados obtenidos con las distintas cantidades de derivaciones son iguales. Esta hipótesis fue rechazada (p−valor <0.001), demostrando por tanto que hay diferencias entre los resultados obtenidos si variamos el número de derivaciones utilizado. Adicionalmente se compararon las distintas opciones de dos en dos, para ver si había diferencias significativas al incrementar el número de derivaciones. Para ello se utilizó el test Wilcoxon [145]. En todos los casos posibles se obtuvo un p-valor menor que 0.001, rechazando siempre la hipótesis nula. Por tanto, hay diferencias significativas en todos los casos al variar el número de derivaciones Posteriormente a la aplicación del test de Friedman, si la hipótesis nula es rechazada, es posible aplicar otros test para obtener un ranking de los números de derivaciones que obtienen 100 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia ● ● ● ● ● ● d2 d4 d6 d8 d10 d12 600 800 1000 1200 1400 Número de derivaciones Número de errores totales Figura 4.6: Resultados del error en toda la base de datos para los distintos números de derivaciones. Los puntos fuera de rango se representan como puntos individuales, las líneas verticales o “bigotes” son proporcionales a la diferencia entre cuartiles y las cajas van desde el valor del primer cuartil Q1 hasta el valor del tercer cuartil Q3, correspondiendo la línea horizontal intermedia con la mediana (Q2). un mejor resultado. Para ello se aplicó el test Nemenyi [100]. Los resultados del test dieron el siguiente orden (de mejor a peor): 12, 10, 8, 6, 4 y 2 derivaciones. También se aplicaron otros test similares (Holm [59] y Finner [41]) obteniendo la misma ordenación. 4.3.3. Discusión La Tabla 4.5 muestra que el error tiende a reducirse cuando se incrementa el número de derivaciones utilizadas. El error se reduce del 0.601% de media usando 2 derivaciones hasta el 0.338% utilizando las 12 derivaciones. Este decremento se puede apreciar en la Figura 4.6, donde se ve que el mayor cambio se produce al pasar de 2 a 4 derivaciones, mientras que al seguir añadiendo derivaciones las mejoras son menos pronunciadas. 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC 101 También se observa que disminuye la variación en los resultados (rango intercuartílico) al añadir más derivaciones. Al revisar las cien ejecuciones realizadas para cada número de derivaciones, se aprecian diferencias en los resultados según las derivaciones concretas que se utilicen (I, II, etc.). Esto es especialmente notable al trabajar con cantidades más pequeñas de derivaciones. Por ejemplo, cuando se utilizan 2 derivaciones, algunas parejas de derivaciones obtienen un error de tan solo 860 latidos en toda la base de datos, mientras que otras parejas llegan hasta los 1400. Esta variabilidad disminuye conforme se aumenta el número de derivaciones utilizadas. Esto se debe a que en algunas de las derivaciones la señal tiene baja calidad (véase Figura 4.7). La evidencia extraída de estas derivaciones es pobre, lo que hace empeorar los resultados del agrupamiento. Cuando se emplea un número elevado de derivaciones el impacto relativo causado por la mala calidad en alguna de las derivaciones es compensado más fácilmente por el resto de las derivaciones. Por ello los resultados son más consistentes y disminuye la variabilidad. Los resultados obtenidos sustentan la hipótesis inicial de que el agrupamiento basado en acumulación de evidencia puede mejorar al utilizar más derivaciones. Lo cual no implica que siempre se obtenga una mejora al utilizar un número superior de derivaciones. Por ejemplo, todas las ejecuciones que emplearon 4 derivaciones mejoran cualquier resultado utilizando 2 derivaciones. Sin embargo, algunas ejecuciones con 10 derivaciones obtienen un mejor resultado que algunas con 12 derivaciones (aunque en promedio se obtienen mejores resultados con 12 derivaciones; véase Figura 4.6). Algunas derivaciones en determinados casos pueden no aportan información nueva e incluso pueden empeorar el resultado añadiendo ruido al agrupamiento (véase Figura 4.7). Sin embargo, no es sencillo seleccionar a priori las mejores derivaciones para agrupar. Esto depende en muchas ocasiones del registro, de los tipos de latidos presentes o de las condiciones de grabación del registro. Por lo que, si no se cuenta con un método fiable para evaluar la calidad de una derivación y decidir si utilizarla o no en el agrupamiento, la mejor opción para aumentar el rendimiento y la robustez del agrupamiento, según los resultados obtenidos, es aumentar el número de derivaciones utilizadas. En la bibliografía de agrupamiento de latidos no hemos encontrado ningún trabajo con el que comparar los resultados obtenidos en esta base de datos. Tampoco hemos encontrado ningún otro trabajo que aplique un algoritmo de agrupamiento a las doce derivaciones del ECG de modo simultáneo. Basándose en una técnica de agrupamiento es posible construir un clasificador asistido [29][109], con la ayuda de un cardiólogo que anote los grupos de latidos 102 Capítulo 4. Agrupamiento de latidos mediante acumulación de evidencia 00:02:05 00:02:15 I II III AVR AVL AVF V1 V2 V3 V4 V5 V6 I II III AVR AVL AVF V1 V2 V3 V4 V5 V6 Figura 4.7: Fragmento de ECG de 12 derivaciones. Obsérvese cómo algunas derivaciones, I o V4, tienen mucho ruido y en ellas es difícil distinguir los latidos y la forma que tienen. (Fuente: St.-Petersburg Institute of Cardiological Technics 12-lead Arrhythmia Database, registro I66, entre 0:02:05 y 0:02:15) (véase Sección 1.2). Los resultados de este clasificador asistido sí posibilitan una comparación con otros algoritmos de clasificación que fueron aplicados sobre la misma base de datos. El clasificador propuesto en [6] obtiene, en el mejor caso, un error del 0.3% sobre la base de datos INCARTDB, pero para obtener dicho error es necesario que se anoten al menos 200 latidos de cada registro. En este mismo trabajo también se muestra el error obtenido sin anotar ningún latido (permitiendo por tanto una comparación más directa con nuestro trabajo), siendo este del 23.6%. En [85] se obtiene un error del 9.38% en la base de datos INCARTDB. En [83], utilizando PCA sobre las 12 derivaciones y wavelets para representar los latidos, se obtuvo un error del 2.88% para la base de datos INCARTDB completa, aunque dividiendo los latidos únicamente en tres clases (en nuestro caso se dividen en las 5 clases recomendadas por la AAMI). 4.3. Agrupamiento de latidos utilizando 12 derivaciones con PN-EAC 103 Los resultados obtenidos avalan la idea de emplear el paradigma de acumulación de evidencia para la integración de información de múltiples derivaciones, así como el interés de usar evidencia negativa a partir de la información extraída de las posiciones de los latidos. Además, el algoritmo es fácilmente adaptable a cualquier número de derivaciones y altamente paralelizable (la evidencia de cada derivación puede obtenerse en paralelo con las demás). La reducción del error al pasar de 2 (0.601%) a 12 (0.338%) derivaciones de casi el 50% hace pensar que si en otras bases de datos, como la MIT-BIH Arrhythmia Database, estuvieran disponibles las 12 derivaciones el rendimiento de PN-EAC algoritmo podría mejorar de forma considerable. 110 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia En la fase de inicialización todos los objetos aparecen en las pparticiones iniciales, por lo que toda la matriz Φse inicializará con este valor, Φ(:,:)=p. Únicamente es necesario almacenar el número de particiones de las que se extrae evidencia positiva ya que será este número el que indique el máximo posible de particiones en las que dos objetos se agruparon juntos. Por último, también es necesario inicializar el vector Y, de tal modo que los oprimeros objetos de Westén vinculados con sus correspondientes objetos de X, (y1=1,y2=2,...,yo=o). En este momento concluye el proceso de inicialización y comienza la ejecución del bucle (véase Algoritmo 8, líneas 2-8) que procesará cada nuevo objeto xzdinámicamente. 5.1.3. Búsqueda del par de objetos más similares en W Al recibir un objeto nuevo es necesario liberar espacio enWpara almacenarlo. Para ello el primer paso es buscar la pareja de objetos (wi,wj)más similares (véase Procedimiento 2). Se comienza calculando la proporción de veces que cada par de objetos wiywjfueron agrupados juntos mediante la división elemento a elemento de las matrices ΩyΦ: Ψ(i,j)=Ω(i,j) Φ(i,j) .(5.4) La matriz Ψproporciona una medida de similitud entre los objetos de W. A mayor valor de Ψ(i,j), mayor la similitud entre los correspondientes objetos wiywj. Por ello, se buscará el máximo de la matriz Ψy por tanto a las parejas de objetos más similares. Debemos identificar todas las celdas de la matriz donde aparece el máximo, ya que las correspondientes parejas de objetos tendrán la misma similitud. Dada la simetría de la matriz únicamente se considerarán las celdas por encima de la diagonal (véase Procedimiento 2, líneas 3-12): I=(i,j)|(i,j) = argmaxl,m(Ψ(l,m)),l<m.(5.5) En el caso de que Isolo contenga una pareja estos serán los objetos a fusionar y se procede al siguiente paso del algoritmo (véase Sección 5.1.4). En otro caso, es necesario determinar cuál de todas las parejas con el mismo valor Ψ(i,j)será fusionada. En el paso anterior se calculó la similitud de dos objetos basándose en si son agrupados juntos. Ahora se medirá de acuerdo a cómo es su agrupamiento con el resto de objetos. El hecho de que un objeto wi tienda a agruparse con los mismos objetos con los que se agrupa wjes un indicio de que wiy wjson similares. 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa (EPN-EAC) 111 Procedimiento 2 Búsqueda pareja de objetos (wi,wj)más similares en W Input: Lista de objetos W=w1,w2,...,wo, Matrices ΦyΩ Output: Pareja de objetos wiywj 1: Function BuscarPareja(W,Φ,Ω) 2: Ψ=Ω/Φ(elemento a elemento) #Ecuación (5.4) 3: Max =−∞,I={} 4: for m=1to o−1do 5: for l=m+1to odo 6: if Ψ(l,m)>Max then 7: I={(wl,wm)},Max =Ψ(l,m)#Ecuación (5.5) 8: else 9: if Ψ(l,m)== Max then I={I,(wl,wm)}end if 10: end if 11: end for 12: end for 13: if size(I) == 1then 14: (wi,wj) = I 15: else 16: MinDistance =∞,J={} 17: for all (wl,wm)∈Ido 18: if distance(Ψ(l,:),Ψ(m,:))<MinDistance then 19: J={(wl,wm)},MinDistance =distance(Ψ(l,:),Ψ(m,:))#Ecuación (5.6) 20: else 21: if distance(Ψ(l,:),Ψ(m,:)) == MinDistance then J={J,(wl,wm)}end if 22: end if 23: end for 24: if size(J) == 1then 25: (wi,wj) = J 26: else 27: MinDistance =∞,K={} 28: for all (wl,wm)∈Jdo 29: if distance(wl,wm)<MinDistance then 30: K={(wl,wm)},MinDistance =distance(wl,wm)#Ecuación (5.7) 31: else 32: if distance(wl,wm) == MinDistance then K={K,(wl,wm)}end if 33: end if 34: end for 35: if size(K) == 1then (wi,wj) = Kelse (wi,wj) = random.pair(K)end if 36: end if 37: end if 38: return (wi,wj) 39: end 112 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia Cada columna o fila de la matriz Ψrepresenta el agrupamiento de un objeto con el resto de objetos de W. Por lo tanto, comparando dos filas o columnas se puede comparar el agrupamiento de dos objetos con respecto a todos los objetos en W. De esta forma se obtiene un nuevo criterio de similitud para encontrar aquellas parejas de objetos cuyo agrupamiento en Wes más similar. Sobre las parejas (wi,wj)que cumplen el criterio anterior (5.5), elegimos solo aquellas que minimizan la distancia en la matriz Ψentre las filas iyj(véase Procedimiento 2, líneas 16-23): J=(i,j)|(i,j) = argminl,mdistancia(Ψ(l,:),Ψ(m,:),∀(l,m)∈I.(5.6) En el caso de que en Jsolo haya una pareja de objetos se procede al siguiente paso del algoritmo. En caso contrario, se medirá la distancia euclídea entre los objetos de cada pareja restante. Aquella pareja con una menor distancia entre objetos será la elegida. Esto se aplica para todas las parejas de J, seleccionando aquella con menor distancia (véase Procedimiento 2, líneas 27-34): K=(i,j)|(i,j) = argminl,m(distancia(wl,wm)),∀(l,m)∈J.(5.7) Finalmente, en el caso de que en Khaya más de una pareja se escogerá una pareja aleatoriamente (véase Procedimiento 2, línea 35). 5.1.4. Fusión de dos objetos e inclusión del nuevo objeto xz En el paso anterior se seleccionó la pareja (wi,wj)de objetos en Wmás similares según los criterios establecidos. En este paso se fusionan dichos objetos, dejando un espacio libre en Wen el cual introducir el nuevo objeto xz(véase Procedimiento 3). El objeto resultado de la fusión de (wi,wj)se calcula como la media de los vectores de características de los dos objetos. El nuevo objeto resultante de dicha fusión ocupará la posición del antiguo objeto wien W: wi←media(wi,wj),(5.8) donde ←indica que wies reemplazado por la media de wiywj. A continuación, se actualizará Yde tal modo que todos los objetos que estaban representados por los dos objetos pasarán a estar representados por el nuevo objeto. Para ello en Ytodos los objetos con valor jpasan a tener valor i, la posición del objeto fusionado (véase Procedimiento 3, líneas 3-5). 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa (EPN-EAC) 113 Procedimiento 3 Fusión de dos objetos wiywj Input: Lista de objetos W=w1,w2,...,wo,Y, Matrices ΦyΩ, nuevo objeto xz Input: Elementos a fusionar wiywj Output: W,Y,ΦyΩactualizados 1: Function FusiónPareja(W,Y,Φ,Ω,xz,wi,wj) 2: wi←media(wi,wj)#Ecuación (5.8) 3: for l=1to z−1do 4: if yl== jthen yl=iend if 5: end for 6: for l=1to odo 7: Ω(l,i)=Ω(l,i)+Ω(l,j),Φ(l,i)=Φ(l,i)+Φ(l,j)#Ecuación (5.10) 8: Ω(i,l)=Ω(l,i),Φ(i,l)=Φ(l,i) 9: end for 10: wj←xz 11: yz=j 12: for l=1to odo 13: Ω(j,l)=0, Ω(l,j)=0, Φ(j,l)=0, Φ(l,j)=0 14: end for 15: return (W,Y,Φ,Ω) 16: end También precisarán actualización las matrices ΩyΦpara mantener la consistencia con el conjunto Wque representan. La evidencia de los dos antiguos objetos será asignada al objeto fusión (véase Procedimiento 3, líneas 6-9): (Ω(:,i))T=Ω(i,:)←Ω(i,:)+Ω(j,:),(5.9) (Φ(:,i))T=Φ(i,:)←Φ(i,:)+Φ(j,:).(5.10) Posteriormente, se añade el nuevo objeto xzal espacio libre correspondiente al antiguo objeto wj,wj←xz, Actualizando consecuentemente Y, de modo que yz=j. El objeto wjes nuevo y no hay ninguna evidencia sobre él. En consecuencia, las correspondientes columnas y filas de ΩyΦson puestas a cero (Ω(:,j)=0, Ω(j,:)=0, Φ(:,j)=0 y Φ(j,:)=0). 5.1.5. Recolección de evidencia Tras añadir el nuevo objeto xzes preciso recopilar evidencia para el conjunto actualizado W, (para los objetos antiguos tenemos evidencia ya recolectada, pero para el nuevo no hay 114 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia ninguna). Para recolectar la evidencia se crean rparticiones nuevas con los objetos de W. La elección de rdebe realizarse cuidadosamente ya que este número de particiones es creado cada vez que se recibe un nuevo objeto, por lo que influirá notablemente en el rendimiento del algoritmo. Un valor de relevado puede acarrear un coste computacional excesivo mientras que uno demasiado bajo puede impedir recopilar suficiente evidencia para agrupar correctamente el nuevo objeto. El proceso de recolección de evidencia es similar al utilizado en la inicialización. De las r particiones creadas se extrae evidencia que se añade a la ya obtenida mediante la expresión: Ω(i,j)=Ω(i,j)+ai j,(5.11) donde ai j es la cantidad de veces que la pareja de objetos wi,wjha aparecido en el mismo grupo en las rparticiones creadas. Al igual que en la inicialización, se podrían crear opcionalmente al mismo tiempo sparticiones para extraer evidencia negativa: Ω(i,j)=Ω(i,j)−bi j,(5.12) siendo bi j el número de veces en que la pareja de objetos wi,wjno aparece en el mismo grupo a lo largo de las sparticiones. También se debe actualizar la matriz Φpara reflejar las rnuevas particiones de evidencia positiva creadas: Φ(i,j)=Φ(i,j)+r.(5.13) De este modo se han generado con el nuevo objeto de W r particiones de las que se habrá extraído evidencia. 5.1.6. Extracción de la partición final Este último paso se ejecuta únicamente cuando sea necesario conocer un resultado parcial del agrupamiento de los objetos procesados hasta un momento dado (véase Procedimiento 4). Toda la información de evidencia acumulada hasta el momento está contenida en las dos matrices ΩyΦ, que serán utilizadas para obtener la partición final. El primer paso para obtener la partición final es calcular Ψ, según (5.4). Con esta operación se obtiene una matriz de similitud sobre la que aplicar un algoritmo de agrupamiento diseñado para matrices de similitud. En nuestro caso se aplica el algoritmo jerárquico de enlace medio, tal y como hicimos en la acumulación de evidencia estática 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa (EPN-EAC) 115 Procedimiento 4 Extraer la partición final Pz ∗ Input: Lista de objetos W=w1,w2,...,wo,Y, Matrices ΦyΩ Output: Pz ∗partición final de los datos 1: Function ExtraerPartFinal(W,Y,Φ,Ω) 2: Ψ=Ω/Φ(elemento a elemento) #Ecuación (5.4) 3: PW ∗=enlace.medio(Ψ)#Sección 4.1.3 4: Pz ∗=PW ∗ 5: for all (wl)∈PW ∗do 6: L={} 7: for m=1to zdo 8: if ym== lthen L={L,(xm)}end if 9: end for 10: Pz ∗.replace(wl).by(L)#Ecuación (5.1) 11: end for 12: return (Pz ∗) 13: end (véase Procedimiento 4, líneas 2-3). Como resultado obtenemos el agrupamiento final PW ∗de los objetos del conjunto W. Si se quiere extender dicho agrupamiento a todos los objetos de Xzse puede hacer mediante el array Y(véase Procedimiento 4, líneas 5-11). Para ello cada objeto de Wserá reemplazado por los objetos correspondientes de Xz(véase Ecuación (5.1)). Con este último paso finalizaría la iteración del algoritmo y quedaría en espera de un nuevo objeto para procesar. 5.1.7. Ejemplo ilustrativo del funcionamiento del algoritmo Para ilustrar el funcionamiento de la técnica se presenta un ejemplo de su ejecución paso a paso. Dada una serie de datos X={11,14,18,49,3,94,...}, elegimos como tamaño de la lista W,o=4. El número de particiones generadas en la inicialización será p=10 y las particiones que se generan en cada iteración será r=5. Solo se generará evidencia positiva. Inicialización El conjunto Wes inicializado con los primeros oobjetos de X,W={11,14,18,49}. Con dichos objetos generamos p=10 particiones de datos de las cuales extraemos evidencia (véase Ecuación (5.2)). Supongamos que la matriz de evidencia Ωobtenida es: 116 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia Ω=      10 5 4 2 5 10 5 2 4 5 10 2 2 2 2 10      . Inicializamos la matriz Φcon el número de particiones: Φ=      10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10      . También se inicializa el array Ycon los primeros oobjetos de Xvinculados, respectivamente, a los primeros oobjetos de W: Y= (1,2,3,4). Primera iteración Llega un nuevo objeto a procesar, x5=3. El primer paso es encontrar la pareja de objetos más similares en W. Para ello obtenemos la matriz Ψaplicando (5.4): Ψ=      1 0.5 0.4 0.2 0.5 1 0.5 0.2 0.4 0.5 1 0.2 0.2 0.2 0.2 1      . Seleccionamos las celdas con el valor máximo sobre la matriz triangular superior, como se indica en (5.5). De las celdas con valor máximo obtenemos las parejas de índices que presentan similitud máxima: I={(1,2),(2,3)}. Dado que hay más de una pareja en I, medimos la distancia entre las correspondientes columnas o filas de la matriz Ψ, y seleccionamos la pareja (o parejas) con una distancia menor (véase Ecuación (5.6)): 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa (EPN-EAC) 117 distancia(Ψ(1,:),Ψ(2,:)) = q(1−0.5)2+(0.5−1)2+(0.4−0.5)2+(0.2−0.2)2 =√0.51, distancia(Ψ(2,:),Ψ(3,:)) = q(0.5−0.4)2+(1−0.5)2+(0.5−1)2+(0.2−0.2)2 =√0.51, J={(1,2),(2,3)}. Utilizando este criterio ambas parejas obtienen el mismo valor de similitud, lo que no nos permite elegir una de ellas. Por tanto, procederemos a medir la distancia euclídea entre los objetos de Wque representan las parejas de J, quedándonos con la pareja con menor distancia (véase Ecuación (5.7)): distancia(w1,w2) = q(11−14)2=√9=3, distancia(w2,w3) = q(14−18)2=√16 =4, K={(1,2)}. La pareja con menor distancia será la de los objetos w1yw2, que procederemos a fusionar aplicando (5.8): w1←media(w1,w2) = media(11,14) = 12.5, W={12.5,14,18,49}. Actualizaremos Ypara indicar que los objetos x1yx2ahora son representados por w1: Y= (1,1,3,4). Realizaremos las actualizaciones necesarias también en las matrices ΩyΦpara fusionar la evidencia de los dos antiguos objetos, como se indica en (5.10): 118 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia (Ω(:,1))T=Ω(1,:)←Ω(1,:)+Ω(2,:) ={10+5,5+10,4+5,2+2}={15,15,9,4}, (Φ(:,1))T=Φ(1,:)←Φ(1,:)+Φ(2,:) ={10+10,10+10,10+10,10+10}={20,20,20,20}, Ω=      15 15 9 4 15 10 5 2 9 5 10 2 4 2 2 10      ,Φ=      20 20 20 20 20 10 10 10 20 10 10 10 20 10 10 10      . En el espacio libre dejado en Wintroduciremos el nuevo objeto x5=3, w2←x5. Al ser un objeto nuevo actualizamos las matrices de evidencia, poniendo a cero las filas y columnas correspondientes a w2: W={12.5,3,18,49}, Ω=      15 0 9 4 0 0 0 0 9 0 10 2 40210      ,Φ=      20 0 20 20 0 0 0 0 20 0 10 10 20 0 10 10      . Se actualiza Ypara indicar que el objeto x5está representado en Wpor el objeto w2: Y= (1,1,3,4,2). Con el conjunto Wactualizado generaremos r=5 particiones de datos nuevas. Obtendremos la evidencia de dichas particiones y la sumaremos a la ya existente aplicando (5.11) y (5.13): Ω=      15 0 9 4 0 0 0 0 9 0 10 2 40210      +      5141 1521 4251 1115      =      20 1 13 5 1 5 2 1 13 2 15 3 51315      . Φ=      20 0 20 20 0 0 0 0 20 0 10 10 20 0 10 10      +      5555 5555 5555 5555      =      25 5 25 25 5 5 5 5 25 5 15 15 25 5 15 15      . 5.1. Agrupamiento dinámico mediante acumulación de evidencia positiva y negativa (EPN-EAC) 119 Para finalizar la iteración podríamos extraer la partición final de los datos si lo deseáramos. Para ello en primer lugar crearíamos la matriz Ψsobre la que aplicaríamos el algoritmo de agrupamiento correspondiente para obtener PW ∗: Ψ=      0.8 0.2 0.52 0.2 0.2 1 0.4 0.2 0.52 0.4 1 0.2 0.2 0.2 0.2 1      , PW ∗={(w1,w3),(w2),(w4)}. Por medio de Ypodríamos trasladar este agrupamiento al conjunto X(véase (5.1): PW ∗={(w1,w3),(w2),(w4)}. Y= (1,1,3,4,2), wm←{xl}∀xl|yl=m, PX ∗={(x1,x2,x3),(x5),(x4)}. Segunda Iteración Llega un nuevo objeto x6=94. Aplicamos (5.4): Ψ=      0.8 0.2 0.52 0.2 0.2 1 0.4 0.2 0.52 0.4 1 0.2 0.2 0.2 0.2 1      . Seguidamente (5.5): I={(1,3)}. Al tener solo una pareja de objetos en Ipodemos fusionarlos aplicando (5.8): w1←media(w1,w3) = media(12.5,18) = 15.25, W={15.25,3,18,49}. Actualizamos Y: Y= (1,1,1,4,2). 126 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia Tabla 5.2: P-valor del test Wilcoxon de significancia para las distintas estrategias dinámicas a pares. 1aEst. vs 2aEst. 1aEst. vs 3aEst. 2aEst. vs 3aEst. 25 Grupos 0.0002 <0.0001 0.1449 Criterio Tiempo Vida 0.0002 <0.0001 0.0575 Tabla 5.3: P-valor del test Wilcoxon de significancia para cada estrategia entre la versión estática y la versión dinámica. 1aEstrategia 2aEstrategia 3aEstrategia 25 Grupos 0.0408 0.0115 0.0039 Criterio Tiempo Vida 0.9033 <0.0001 0.0124 tiempo de vida y un número fijo de 25 grupos). En este caso también se dividió el número de errores por el número de latidos totales (109966 latidos) para obtener el porcentaje de error. Para los resultados con el número fijo de 25 grupos se obtuvieron unos porcentajes de error de 2.78%, 1.69% y 1.41% para la primera, segunda y tercera estrategias, respectivamente; mientras que en el caso del criterio del tiempo de vida los porcentajes son de 2.95%, 1.21% y 0.87%, respectivamente, teniendo en cuenta que en este caso el número de grupos creado es variable, con una media de 9.22, 53.22 y 76.66 grupos por registro para la primera, segunda y tercera estrategia, respectivamente. Sobre los resultados de la Tabla 5.1, al igual que en la estrategia estática, se aplicaron test estadísticos para verificar la significancia de las diferencias entre las distintas estrategias. Los resultados obtenidos por el test de Wilcoxon son mostrados en la Tabla 5.2. También se aplicaron test estadísticos para verificar la significancia de las diferencias en cada estrategia entre la versión estática y la versión dinámica. Los resultados se muestran en la Tabla 5.3. Además, para la tercera estrategia y un número de grupos fijo de 25 (error del 1.41%), se calculó también la matriz de confusión. Se eligió esta estrategia por ser la que obtuvo un mejor resultado con un número de grupos por registro reducido. Dicha matriz de confusión se muestra con las anotaciones originales de la base de datos MIT-BIH Arrhythmia Database y con las anotaciones recomendadas por la AAMI en las Tablas 5.4 y 5.5, respectivamente. En estas matrices las columnas representan el tipo real de los latidos según las anotaciones y las filas el tipo en que se agrupó. 5.3. Resultados 127 Tabla 5.4: Matriz de confusión para la base de datos MIT-BIH Arrhythmia Database utilizando las anotaciones originales de la base de datos. N L R a V F J A S E j P Q ! e f N74781 0 1 9 70 199 4 175 2 0 104 0 10 0 14 18 L0 8068 0 0 0 0 0 0 0 0 0 0 1 0 0 0 R0 0 7185 0 5 1 29 29 0 2 5 0 0 7 0 0 a0 0 0 135 4 0 0 7 0 0 0 0 0 0 0 0 V33 0 0 3 7018 82 0 0 0 1 0 0 2 6 0 0 F10 0 0 0 22 521 0 0 0 0 1 0 0 0 0 0 J0 0 0 0 0 0 50 0 0 0 0 0 0 0 0 0 A117 0 70 3 1 0 0 2227 0 0 0 0 0 0 2 0 S0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E0 0 0 0 1 0 0 0 0 103 0 0 0 97 0 0 j72 0 0 0 0 0 0 0 0 0 119 0 0 0 0 0 P0 0 0 0 0 0 0 0 0 0 0 6873 0 0 0 39 Q0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 1 !0 4 0 0 8 0 0 106 0 0 0 0 0 362 0 0 e0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 f3 0 0 0 1 0 0 0 0 0 0 151 15 0 0 924 Se 99.69 99.95 99.02 90.00 98.43 64.88 60.24 87.54 0.00 97.17 51.97 97.85 15.15 76.69 0.00 94.09 P+ 99.20 99.99 98.93 92.47 98.22 94.04 100.00 92.02 - 51.24 62.30 99.44 83.33 75.42 - 84.46 128 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia Tabla 5.5: Matriz de confusión para la base de datos MIT-BIH Arrhythmia Database utilizando las anotaciones recomendadas por la AAMI. N S V F Q N90035 371 84 200 29 S259 2543 5 0 0 V37 109 7596 82 2 F10 1 22 521 0 Q3 0 1 0 8008 Se(%) 99.66 84.09 98.55 64.88 99.61 P+(%) 99.25 90.59 97.06 94.04 99.95 5.4. Discusión Es posible apreciar en los resultados con 25 grupos una mejora significativa entre la primera estrategia y la segunda estrategia, pasando del 2.78% al 1.69% (p=0.0002). De la segunda estrategia a la tercera se produce otra mejora del 1.69% al 1.41%, pero no es significativa estadísticamente (p=0.1449). Ocurre algo semejante al utilizar el criterio del tiempo de vida, siendo significativa la diferencia entre la primera y la segunda estrategia del 2.95% al 1.69% (p=0.0002), mientras que la diferencia entre la segunda y tercera estrategia no supera el umbral de 0.05 (p=0.0575). Los resultados son similares a los obtenidos anteriormente para el agrupamiento estático (véase Sección 4.2.2), reforzando las conclusiones ya alcanzadas en dicho capítulo con la técnica PN-EAC. En esta ocasión se aprecia una mayor diferencia entre la primera estrategia y la segunda, al contrario que en la versión estática donde las diferencias más significativas se obtenían entre la segunda y tercera estrategias. Por otra parte, se corrobora, una vez más, el pobre comportamiento del criterio del tiempo de vida en la segunda y tercera estrategias. En este caso se estableció un límite superior de 99 grupos por registro que fue alcanzado en un alto número de registros en la segunda estrategia y por la mayoría en la tercera estrategia. Por ello de nuevo se reafirma la no idoneidad de este criterio para el problema del agrupamiento de latidos, especialmente al combinar evidencia de distintas representaciones e incluir evidencia negativa. Por ello, el resto del análisis lo basaremos sobre los resultados obtenidos utilizando 25 grupos. Comparando los resultados con los obtenidos en el Capítulo 4 es posible apreciar que la adaptación de la técnica (de estática a dinámica) no parece haber empeorado el resultado en 5.4. Discusión 129 términos de error. De hecho, los resultados para la segunda y tercera estrategias han mejorado ligeramente. En la técnica estática (PN-EAC) obteníamos porcentajes de error del 2.25%, 1.81% y 1.44% para la primera, segunda y tercera estrategias, respectivamente, mientras que para la versión dinámica (EPN-EAC) se obtienen porcentajes de error del 2.78%, 1.69% y 1.41%, siendo significativas todas estas diferencias entre las versiones estáticas y dinámicas significativas (véase Tabla 5.3). Esta ligera mejora para la segunda y tercera estrategias se debe principalmente a que en la versión dinámica recopilamos evidencia muchas más veces, incluso si la hacemos sobre un conjunto de latidos más pequeño. Es posible comparar los resultados mostrados en la Tabla 5.4 con los obtenidos por [77] y los obtenidos por [20]. En esta comparación podemos ver cómo EPN-EAC obtiene mejor sensibilidad en 7 (L, a, V, J, E, Q, f) de las 14 clases (clases S y e no fueron consideradas al no estar suficientemente representadas), siendo mejores [20] y [77] en 5 (R, F, A, j, P) y 2 (N, !) de las clases, respectivamente. En cuanto al valor predictivo positivo, EPN-EAC obtuvo el mejor resultado para 7 (N, L, R, V, F, J, P) clases, mientras que [20] y [77] lo obtuvieron para 3 (A, !, f) y 4 (a, E, j, Q) clases, respectivamente. Finalmente, el error total obtenido es del 1.41%, ligeramente inferior a los obtenidos por [20] y [77] de 1.44% y 1.50%, respectivamente. Los peores resultados se obtienen para los latidos de fusión de marcapasos (F), los latidos prematuros supraventriculares (S) y atrioventriculares (J), los latidos de escape auricular (e), y los latidos inclasificables (Q), con los cuales el algoritmo se comporta de forma poco satisfactoria. El tipo F corresponde a la fusión de latidos de marcapasos y latidos normales, una distinción que es inherentemente difícil de hacer, incluso para los especialistas (véase Figura 5.1). Durante la ejecución, más de un 22% de los latidos de este tipo fueron agrupados incorrectamente, principalmente como latidos normales (N). Los latidos de escape auricular (e) y los latidos prematuros supraventriculares (S), auriculares (A) y atrioventriculares (J), además de estar muy poco representados en la base de datos, tienen una morfología muy similar a la de los latidos normales, por lo que resulta muy difícil para un algoritmo distinguirlos (véanse Figuras 4.3 , 5.2 y 5.3). Por otra parte, el tipo Q aglutina a aquellos latidos que en la base de datos no se han podido clasificar en ninguna de las categorías, debido a una morfología no reconocida o al ruido de la señal, por lo que se trata de una clase sin una característica morfológica reconocible. Los mejores resultados se obtienen para los latidos con bloqueo de rama izquierda (L) y para los latidos ventriculares (V), obteniendo en ambos casos mejor sensibilidad y valor predictivo positivo que los otros dos algoritmos. En ambos tipos de latido la morfología es 130 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia N F N F N F N F N MLII Figura 5.1: Fragmento de un ECG en el que se alternan latidos normales con latidos del tipo F. (Fuente: MIT-BIH Arrhythmia Database, registro 213, entre 0:02:07 y 0:02:12) claramente distinta a la del complejo QRS normal y asimétrica, por lo que la representación elegida es especialmente adecuada para distinguir estos tipos de latido. Además, ambos grupos de latidos están entre los más numerosos, después de los latidos normales. En general, los tres trabajos, al estar basados en la morfología del complejo QRS, muestran un rendimiento más pobre al distinguir latidos con una morfología similar entre sí, que para distinguirse requieren de información de la onda P o, en su defecto, información derivada de la distancia entre latidos. El algoritmo aquí presentando podría verse beneficiado de incorporar información de la onda P como forma de evidencia negativa, de forma que ayudara a separar estos casos y al mismo tiempo el resultado del agrupamiento no se viera afectado por la inestabilidad de su detección. En [77] con un número fijo de 25 grupos por registro se obtuvo un error del 1.51% para la base de datos MIT-BIH Arrhythmia Database completa. Los resultados de este capítulo con la primera y segunda estrategias obtienen porcentajes de error mayores, mientras que la tercera estrategia obtiene un error ligeramente menor (1.41%). Los resultados de sensibilidad obtenidos son algo inferiores a los obtenidos por [77], pero se compensa con un mejor resultado en el valor predictivo positivo. Al comparar los resultados con los obtenidos por [20] es posible ver una gran semejanza, variando muy poco el resultado del error total sobre la base de datos MIT-BIH Arrhythmia Database, del 1.44% obtenido por [20] al 1.41% obtenido aquí. Los valores de sensibilidad obtenidos por [20] son ligeramente superiores, mientras que los de valor predictivo positivo son ligeramente inferiores. Es necesario recalcar que en [20] no se utiliza un número fijo de grupos de latidos, sino que dicho número se decide dinámicamente para cada registro. Sin embargo, a diferencia de la técnica aquí presentada, la información derivada de la distancia entre latidos no está integrada en el algoritmo principal de agrupamiento, sino que se aplica 5.4. Discusión 131 NNNNN S N N S V V V V VF MLII Figura 5.2: Fragmento de un ECG en el que se puede apreciar como los tipos de latido S y N son casi idénticos en la morfología del complejo QRS, distinguiéndose únicamente por tener la onda P invertida. (Fuente: MIT-BIH Arrhythmia Database, registro 208, entre 0:17:45 y 0:17:55) N N N N e AN N N N NN N N MLII Figura 5.3: Fragmento de un ECG con latidos del tipo e y A entre latidos normales (N). (Fuente: MIT-BIH Arrhythmia Database, registro 223, entre 0:21:00 y 0:21:10) en una etapa posterior. Esta etapa en muchos casos incrementa considerablemente el número de grupos llegando hasta 96 grupos en el registro 207. Por ello, allí se aplica una última fase offline de fusión de grupos, limitando a un máximo de 25 el número de grupos por registro. También indicar que en [20] se ha parametrizado y adaptado el método utilizando conocimiento experto de cardiólogos y valores extraídos de la bibliografía de cardiología, algo que se ha evitado en este método y que podría haber mejorado el rendimiento. Casi la totalidad de latidos del tipo E, que representa a los latidos ventriculares de escape, aparecen en el registro 207 en el que es posible apreciar que, aun siendo del mismo tipo, los latidos pueden variar considerablemente en morfología e incluso en la distancia entre latidos (véase Figura 5.4). Por ello, este tipo de latidos resultan extremadamente difíciles para un clasificador [123]. Sin embargo, un algoritmo de agrupamiento, como el que aquí se presenta, puede separar fácilmente en distintos grupos las diversas morfologías de este tipo de latidos (véase Tabla 5.4). Nótese que EPN-EAC puede asumir adecuadamente un incremento del número de características que representan al latido. Esto permite incluir muchas otras características como la potencia espectral, la media, la curtosis, etc. Estas características podrían agruparse en los vectores que fueran necesarios y que posteriormente servirían para generar particiones mediante K-means. La complejidad del proceso de extracción de evidencia depende únicamente del parámetro oelegido y del número de particiones generadas, pero no del 132 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia MLII E E E E E E Figura 5.4: Fragmento de un ECG con latidos del tipo E, con distintas morfologías y distancia entre latidos. (Fuente: MIT-BIH Arrhythmia Database, registro 207, entre 0:27:20 y 0:27:30) número de características que representan el latido. Por otra parte, a la hora de buscar la pareja más similar la complejidad del tercer criterio utilizado (5.7), en el que se calcula la distancia euclídea, sí que depende del número de características utilizadas para representar al latido, pero generalmente no será aplicado (en las pruebas realizadas se aplicó en menos de un 1% de las ocasiones). También indicar que en el algoritmo por defecto se extrae evidencia con cada objeto nuevo que se recibe, pero esta acción podría ejecutarse solo cuando un número determinado de nuevos objetos estén disponibles para mejorar la eficiencia computacional. En el caso de utilizar las anotaciones de la AAMI (véase Tabla 5.5), es posible comparar nuestro trabajo con [20], [30], [31], [77] y [85]. También se comparan los resultados obtenidos con los de PN-EAC del capítulo anterior. Podemos ver esta comparación en la Tabla 5.6, donde se aprecia que EPN-EAC y PN-EAC obtienen unos resultados comparables al resto de resultados en la bibliografía y mejores en varios casos. EPN-EAC obtiene mayor sensibilidad para latidos ventriculares (V), mayor valor predictivo positivo para los latidos de fusión (F), mayor F1 para latidos normales (N) y mayor precisión total. PN-EAC por otra parte obtiene mayor sensibilidad para latidos normales, mayor valor predictivo positivo para latidos supraventriculares (S), latidos ventriculares y latidos inclasificables (Q), junto con mayor F1 para latidos ventriculares, de fusión e inclasificables. Por tanto, entre ambos algoritmos obtienen los mejores resultados para casi todos los tipos de latido. El valor F1 nos da una medida general del rendimiento del algoritmo en una determinada clase, combinando sensibilidad y valor predictivo positivo. El algoritmo EPN-EAC obtiene el mejor resultado para los latidos normales, que son el grupo más numeroso en la base de datos. Por ello, la precisión total será mayor para este algoritmo. Además, para los latidos ventriculares e inclasificables solo es superado en este valor por el algoritmo PN-EAC. Estas tres clases de latidos juntas representan el 96.65% de los latidos de la base de datos. Los peores resultados se obtienen para los latidos supraventriculares (S), que requieren de información de la onda P o información de la distancia entre latidos. 5.4. Discusión 133 Por otra parte, el algoritmo PN-EAC obtiene el mejor F1 para los latidos ventriculares, de fusión e inclasificables, pero al ser estos tipos de latidos menos numerosos en la base de datos la precisión total es algo más baja, aunque cercana a la obtenida por el resto de métodos de la bibliografía. 134 Capítulo 5. Agrupamiento dinámico de latidos mediante acumulación de evidencia Tabla 5.6: Comparación con trabajos previos utilizando las anotaciones recomendadas por la AAMI. “Se” representa la sensibilidad, “P+” el valor predictivo positivo y “F1” el Valor-F. † son los resultados obtenido por EPN-EAC en este capítulo y ‡ los resultados obtenidos por PN-EAC en el capítulo anterior. N S V FQ Se P+ F1 Se P+ F1 Se P+ F1 Se P+ F1 Se P+ F1 Acc † 99.66 99.25 99.45 84.09 90.59 87.22 98.55 97.06 97.80 64.88 94.04 76.78 99.61 99.95 99.78 98.89 ‡99.84 98.77 99.31 71.00 95.38 81.40 97.09 98.81 97.95 81.20 86.13 83.59 99.68 99.96 99.82 98.71 Castro2015 99.58 99.25 99.41 87.54 92.01 89.72 96.34 96.49 96.41 75.56 87.57 81.12 99.32 99.7 99.51 98.84 Lagerholm2000 99.05 99.66 99.35 90.89 83.1 86.82 98.01 96.05 97.02 87.32 75.06 80.73 99.95 99.56 99.75 98.77 DeChazal2004 86.86 99.16 92.60 75.94 38.53 51.12 77.74 81.59 79.62 89.43 8.57 15.64 00 - 85.88 DeChazal2006 94.3 99.36 96.76 87.72 46.95 61.16 94.34 94.3 94.32 73.97 29.15 41.82 00 - 93.89 Llamedo2011 77.55 99.47 87.15 72.88 41.34 52.76 91.25 95.94 93.54 94.69 3.67 7.07 -- - 78.3 Conclusiones y trabajo futuro En esta memoria hemos presentado un conjunto de soluciones que tienen por objeto el procesamiento eficiente del electrocardiograma con el fin último de su interpretación, y que abordan desde su representación hasta el agrupamiento morfológico del complejo QRS. En este capítulo se analizarán las principales aportaciones del trabajo realizado y su continuación en el futuro próximo. A nivel de representación se ha optado por utilizar una base de funciones: la de polinomios de Hermite. Esta base de funciones constituye un conjunto ortonormal que proporciona un modelo paramétrico del complejo QRS. Si bien esta propuesta no es nueva en el ámbito del procesamiento electrocardiográfico, su uso adolece en la bibliografía científica de una cierta arbitrariedad en la selección del conjunto de polinomios a utilizar, cuya decisión carece de una argumentación objetiva, y se basa por regla general en términos meramente visuales. En la presente memoria se muestra un estudio de la representación óptima del complejo QRS mediante polinomios de Hermite, a partir de criterios basados en la teoría de la información, y utilizando como referencia medidas como AIC o BIC. Dicho estudio se ha aplicado a la MIT-BIH Arrhythmia Database, base de datos de referencia en el análisis computacional del electrocardiograma, y que presenta la mayor complejidad morfológica de latidos que se puede encontrar en una base de datos de referencia. El estudio aquí propuesto muestra una excesiva simplificación en la representación del complejo QRS que se realiza habitualmente en la bibliografía. En tanto que AIC y BIC son robustos al ruido blanco, esa simplificación excluye necesariamente de la representación a algunos de los fenómenos fisiológicos que subyacen en el trazado electrocardiográfico, quedando fuera de cualquier análisis posterior. Así todo, el coste computacional del cálculo de la representación del complejo QRS basada en funciones de Hermite crece de un modo no lineal con el número de funciones. Por