Luis David Avendaño Valencia¹, Luis Enrique Avendaño², José María Ferrero³ y Germán Castellanos Domínguez⁴

¹ Ingeniero electrónico. Estudiante maestría, automatización industrial, Universidad Nacional de Colombia, Manizales. Investigador, Grupo de Control y Procesamiento Digital de Señales. ldavendanov@unal.edu.co
² Ingeniero electrónico. DEA bioelectrónica, Universidad Politécnica de Valencia, España. M.Sc., Ingeniería, Universidad de los Andes, Colombia. Profesor titular, Universidad Tecnológica de Pereira, Colombia. leavenda@utp.edu.co
³ Ph.D., en bioingeniería. Investigador, Grupo de Bioelectrónica, Universidad Politécnica de Valencia, España. jmferrer@eln.upv.es
⁴ Ph.D. en ingeniería, Nauchno Isseledovatelskiy Institut, MTUCI, Russia. Profesor titular, Universidad Nacional de Colombia, Manizales. gcastell@telesat.com.co, cgcastellanosd@unal.edu.co

RESUMEN

En este artículo se presenta el desarrollo de un filtro para la reducción de la interferencia de línea de potencia en señales electrocardiográficas (ECG), basado en estimación dual de parámetros y de estado, empleando la filtración Kalman, en el cual se consideran modelos independientes entre la interferencia de línea de potencia y la señal ECG. Ambos modelos son combinados para simular la señal ECG medida sobre la que se realiza la estimación de estado para separar la señal de la interferencia. El algoritmo propuesto es sintonizado y comparado en un conjunto de pruebas realizadas sobre la base de datos QT de electrocardiografía. Inicialmente se hacen pruebas de sintonización del algoritmo para el rastreo de la señal ECG limpia, cuyos resultados son utilizados después para las pruebas de filtrado. Luego se llevan a cabo pruebas exhaustivas sobre la base de datos QT en la filtración de interferencia de línea de potencia, la cual ha sido introducida artificialmente en los registros, para una relación de señal a ruido (SNR) dada, obteniendo así curvas del desempeño del algoritmo, que permiten a su vez comparar con el desempeño de otros algoritmos de filtración, a saber, un filtro notch recursivo de respuesta infinita al impulso (IIR) y un filtro de Kalman, basado en un modelo más simple para la señal ECG. Como resultado, se demuestra que el algoritmo de filtrado obtenido es robusto a los cambios de amplitud de la interferencia; además, conserva sus propiedades para los diferentes tipos de morfologías de señales ECG normales y patológicas.

Palabras clave: filtro de Kalman, reducción de línea de potencia, transformada delta, electrocardiograma.

ABSTRACT

This paper presents a filter for reducing powerline interference in electrocardiographic signals (ECG), based on dual parameter and state estimation using with a Kalman filter. Two models were used to represent power-line interference and ECG signal. Both models were combined to simulate the ECG signal whose state was estimated for separating the ECG signal from the interference. The proposed algorithm was fine-tuned and compared using a set of tests relying on the QT arrhythmia database. Tuning tests were done for tracking clean ECG; these results were used for setting the algorithm’s parameters for later filtering tests. Exhaustive filtering tests were carried out on artificially corrupted database registers for given signal to noise ratios; performance curves were thus obtained, leading to comparing the proposed algorithm with other filtering methods. The proposed algorithm was compared to an recursive infinite impulse response filter (IIR) and a Kalman filter based on a simpler model. A filtering algorithm was thus obtained which is robust for changes in interference amplitude and keeps these properties for different types of ECG morphologies.

Keywords: Kalman filtering, power line reduction, dual state and parameter estimation, delta operator, electrocardiogram.

Recibido: agosto 21 de 2006
Aceptado: octubre 31 de 2007

Introducción

Un problema básico en el análisis automático de electrocardiografía está en la reducción de perturbaciones que afectan sus registros, tales como la interferencia de línea de potencia, de línea base, el ruido electromiográfico, artefactos por desconexión de electrodos, etc. En particular, se analiza la interferencia de línea de potencia causada porque los cables que transportan la señal ECG del paciente al dispositivo de monitoreo son susceptibles a la interacción electromagnética de la línea de alimentación (con frecuencia nominal 60 Hz), haciendo que el análisis y la interpretación del ECG sean más complicados, ya que la delineación de ondas de baja amplitud no es confiable y se introducen ondas falsas (Huhta y Webster, 1973).

La complejidad del filtrado de la interferencia de alimentación está condicionada por la variabilidad en el tiempo de su frecuencia, cuyo rango espectral se traslapa con el rango de frecuencias de la señal ECG. Además, la influencia de los armónicos superiores de la interferencia no puede ser disminuida en la medida en que se tiene una tasa baja de muestreo para las señales electrocardiográficas. En el proceso de reducción de la interferencia de la línea de potencia se utilizan diversas técnicas, por ejemplo, seleccionando el lugar de registro de las señales electrocardiográficas rodeado de pocos dispositivos eléctricos en operación o aislando y conduciendo adecuadamente a tierra el dispositivo de registro; sin embargo, aún sigue siendo necesario el procesamiento digital de la señal, para lo que se han planteado diversos métodos: filtros recursivos de respuesta infinita al impulso (IIR), filtros adaptativos, filtro de Kalman.

Los filtros IIR son la estrategia más simple para esta tarea; no obstante, tienen el inconveniente de la relación entre el ancho de banda del filtro y la duración del período transitorio del mismo. En (Er, 1991; Jang & Chicharo, 1993) se describe un procedimiento para la inicialización del filtro y reducir el transitorio; no obstante, esto sólo es aplicable al comienzo del proceso de filtrado, ya que este problema persiste en los segmentos QRS, que pueden ser vistos por el filtro como impulsos, después de los cuales se genera de nuevo un período transitorio, que se ve reflejado en la aparición de distorsiones al final del mismo, y que pueden eventualmente ser interpretadas como algún tipo de patología (Sörnmo y Laguna 2005).

Al utilizar técnicas de filtrado no lineal, tales como los filtros adaptativos, se puede generar un proceso menos sensible a los transitorios. Los filtros adaptativos se basan en la idea de sustraer una sinusoide, esta se adapta a la interferencia de línea de potencia presente en la señal observada (Ahlstrom y Tompkins, 1985; Hamilton, 1996; Pei y Tseng, 1995). El proceso de adaptación se ve reflejado en el ajuste de los polos de un filtro lineal, dependiente de las condiciones de la señal (Sörnmo y Laguna, 2005); sin embargo, el ajuste del parámetro de adaptación conlleva un compromiso entre la adaptabilidad del filtro, la introducción de distorsión a la señal y el tiempo de adaptación o de convergencia.

Otra estrategia consiste en estimar la amplitud y la fase de la señal interferente en el segmento isoeléctrico, seguida de la substracción de la sinusoide estimada en el latido completo (Levkov et al., 1984; Dotsinsky y Daskalov, 1996). El procedimiento se basa en el cálculo de coeficientes para ajustar las funciones base, seno y coseno en la frecuencia de la interferencia, minimizando el error cuadrático medio. Este método es efectivo, pero requiere la segmentación previa del latido y además no se puede realizar en tiempo real, ya que, como se dijo antes, se necesita al menos un latido completo de la señal ECG para poder segmentarlo y filtrarlo. Además de las dificultades mencionadas, se puede agregar que este método de filtrado pierde efectividad cuando la interferencia sinusoidal presenta variaciones en los parámetros de amplitud, frecuencia y fase. La estrategia planteada por Ziarani y Honrad (2002) consiste en la estimación adaptativa de los parámetros de la interferencia, incluyendo su frecuencia, la cual ha sido previamente aislada de la señal ECG a través de un filtro pasa altas. La sinusoide estimada posteriormente es removida de la señal ECG, con la corrección previa de la fase, agregada por el filtro. Aunque el filtro tiene características que lo hacen muy adecuado para la estimación de sinusoides en ambientes ruidosos, pierde efectividad debido a la baja frecuencia de muestreo de la señal y a la pobre discretización que se aplica para llevar el modelo de forma continua a discreta.

El filtro de Kalman ha sido aplicado a la resolución de este problema, por su proceso robusto. En Mneimneh et al. (2006) se propone un filtro de Kalman adaptativo para la eliminación de la interferencia de línea base, basado en un modelo autorregresivo que modela tanto la señal ECG como la interferencia de línea base. Este filtro funciona adecuadamente para cambios de la interferencia de línea base lentos; sin embargo, en cambios rápidos pierde la sintonía. En Mamen et al. (2006) se utiliza un modelo no lineal, basado en la suma de campanas gaussianas, para modelar la señal ECG. Sobre este modelo se realiza la estimación de la señal electrocardiográfica observada con el filtro extendido de Kalman. Dicho método tiene gran potencia y permite eliminar cualquier tipo de interferencia; solo que su costo computacional es bastante alto, pues requiere la obtención de un modelo diferente en cada derivación del ECG.

Basados en los anteriores trabajos con filtración Kalman, en el presente trabajo, se plantea la reducción de la interferencia de la línea de potencia en dos partes: con relación a la señal de 60 Hz, se emplea un oscilador sinusoidal para su modelado; al modelo del oscilador se aplica el operador delta, con el fin de mejorar su robustez en el tratamiento digital (Middleton y Goodwin, 1986; Middleton y Goodwin, 1990). Con respecto a la señal electrocardiográfica, se usa como modelo un proceso autorregresivo variante en el tiempo, que sirve como predictor lineal, cuyo orden se obtiene a partir del criterio de información de Akaike. Los parámetros de este modelo se ajustan de forma adaptativa mediante un filtro de Kalman, el cual minimiza el error cuadrático medio MSE de predicción (Tarvainen et al., 2005). Luego se aplica el filtro de Kalman dual (Wan y Nelson, 1997), con el objetivo de estimar de forma conjunta los parámetros del proceso autorregresivo que modela la señal electrocardiográfica y el estado, consistente de la señal ECG limpia y del oscilador que modela la interferencia.

Métodos

Transformada delta

Una forma para describir sistemas en tiempo discreto es utilizando el operador de desplazamiento temporal q, definido por

Usando el operador (1), un sistema discreto lineal en espacio de estados se escribe como

El operador de desplazamiento temporal q, (1), es ampliamente usado para describir sistemas en tiempo discreto (Aström y Wittenmark, 1997). No obstante, una desventaja de este operador es que no funciona como su homólogo en tiempo continuo d/dt. De forma intuitiva, se infiere una mejor correspondencia entre tiempo continuo y discreto si el operador de desplazamiento temporal fuera reemplazado por un operador de diferencia, obteniendo una mejor aproximación de la derivada. De esta forma, se define el operador delta a través de la diferencia directa

Así, la forma equivalente de la ecuación (1) es

La relación entre δ y q es una función lineal y, por lo tanto, δ ofrece la misma flexibilidad para modelar sistemas discretos, de la forma como q lo hace. La elección del operador es función de la aplicación particular. Generalmente, q lleva a expresiones más simples y enfatiza la naturaleza secuencial de las señales muestreadas (Middleton y Goodwin, 1986). Por otra parte, dado que δ es una diferencia, lleva a modelos muy similares a aquellos modelados con el operador d/dt, por lo que los resultados obtenidos en tiempo continuo, también se pueden utilizar en aplicaciones de tiempo discreto.

Las ecuaciones de estado se pueden entonces expresar en forma δ como sigue:

En general la señal se encuentra descrita en forma discreta. Para transformarla a la forma δ se aplica la relación

De esta forma, pueden ser establecidas las siguientes equivalencias:

Donde A_q, B_q y C_q son las matrices en el modelo de desplazamiento temporal dadas en (2a) y (2b) e I es la matriz identidad.

Aunque la forma dada en (7) es correcta, no es el mejor método para evaluar el modelo delta, dado que se generan problemas asociados con el proceso de conversión del desplazamiento temporal (Middleton & Goodwin, 1990). Una mejor estrategia consiste en la evaluación directa de la transformada delta a partir de las ecuaciones en espacio de estado en tiempo continuo:

Donde

Este resultado se genera directamente a partir de la consideración de una expansión en series de potencia de A_δ=(e^Ah-I)/h. Para determinar se aplican métodos convencionales, tales como el de valores propios, teorema de Cayley–Hamilton, etc. A partir de (8) y (9) se ve que hay una relación cercana entre los modelos en el dominio delta y los modelos en el dominio del tiempo continuo, dado que →I cuando h→0.

Modelo de interferencia sinusoidal

La interferencia de la línea de potencia se puede modelar como una señal sinusoidal no estacionaria, con cambios aleatorios de amplitud, fase y frecuencia. La aplicación del filtro de Kalman para estimar y luego eliminar esta interferencia requiere la obtención de un modelo discreto que describa su dinámica.

Una forma de obtener un modelo de oscilador sinusoidal se basa en la representación general para un sistema discreto de segundo orden, descrito por las siguientes ecuaciones en espacio de estado:

donde r es el radio de los polos y ω_n=2πf/f_s es la frecuencia de oscilación del sistema. En este caso, si r=1, el sistema no tiene atenuación y se obtendrá un oscilador. Aunque este modelo es simple de obtener, no posee la robustez necesaria para rastrear una señal sinusoidal, como se mostrará más adelante. Un modelo mejorado se puede obtener cuando se parte de un oscilador sinusoidal continuo, descrito por la siguiente ecuación diferencial en espacio de estado

donde ω_n=2πf, corresponde a la frecuencia de oscilación. Aplicando las ecuaciones (8) y (9), se obtiene:

Y

con lo cual se obtiene la matriz de transición de estado en el dominio delta dada por

La matriz de medición C_δ es equivalente a la del modelo continuo, por lo que

Por lo tanto, el modelo planteado en (11b) se transforma al dominio delta como sigue:

Ahora: para efectos de proceso computacional, aplicando la relación (6), el modelo delta en (13b) se lleva a la forma discreta:

El modelo dado en (14b) posee mejores características que el descrito en (10b) debido a que mediante la transformada delta se logra una mejor aproximación de la dinámica continua subyacente en el proceso discreto (Middleton y Goodwin, 1986).

Reducción de interferencia de línea de potencia en señales electrocardiográficas

Puesto que la señal electrocardiográfica puede presentar cambios de ritmo (variabilidad en el ritmo cardíaco, HRV) o de forma, debido a la presencia de patologías o a cambios en la derivación que se esté analizando, no se puede aplicar directamente el filtro de Kalman. Para ello se requiere modelar la señal. Se han probado diferentes modelos haciendo uso de las técnicas de las series temporales.

Una forma simple para modelar la señal cardíaca es utilizar un modelo autorregresivo de primer orden (AR(1)), que tiene características espectrales similares a las del electrocardiograma. Esta suposición hace que el tamaño del espacio de estado conste sólo de tres estados: dos correspondientes al oscilador, y un tercero, modelo AR(1) que emula la dinámica del electrocardiograma, es decir:

Con

y donde v=[v₁ 0 v₂]^T corresponde al ruido del proceso, un ruido blanco gaussiano de media cero y matriz de covarianza R_v, mientras que n(k) corresponde al ruido en la medida, igualmente blanco gaussiano, de varianza r_n. Así, si se aplica el filtro de Kalman sobre la señal electrocardiográfica contaminada con ruido de 60 (Hz), aplicando el modelo (15b) en el estado x₁ se obtendría el estimado de la interferencia de línea de potencia, mientras que en el estado x₃ se seguiría el estimado de la señal electrocardiográfica libre de ruido.

Tomando como punto de partida a la representación anterior se puede llegar a un modelo más completo de la señal electrocardiográfica, utilizando un proceso autorregresivo variante en el tiempo de mayor orden, que teóricamente tendría mayor capacidad de emular la dinámica del electrocardiograma. Se puede utilizar una estrategia de estimación dual de parámetros y de estado, en la cual se realiza la estimación de los parámetros del proceso AR(p) (donde p es el orden del proceso) y, con base en este conjunto de parámetros estimados, efectuar la estimación del estado de la señal contaminada.

De acuerdo con lo anterior, se supone que la señal electrocardiográfica se genera por la siguiente relación de recurrencia:

donde a_j(k) es el valor del parámetro j en el instante k y e(k) corresponde a una secuencia de ruido blanco gaussiano, de media cero y varianza r_e, que representa el error de estimación de la secuencia del ECG. Para hacer la notación más compacta, se hacen las siguientes definiciones:

obteniendo la siguiente representación de un modelo AR(p) variante en el tiempo:

A su vez, se supone que los parámetros del proceso AR(p) están sujetos a variaciones aleatorias de la siguiente forma:

donde r(k) corresponde a un proceso aleatorio gaussiano de media cero y covarianza

haciendo que cada uno de los pesos a_i,i=1,...,p varíe de forma independiente de los demás.

Si sólo se desea modelar la señal electrocardiográfica, las ecuaciones (16) y (17) sirven como representaciones en el espacio de estado, suficientes para modelarla. Por lo tanto, el problema de estimación de parámetros se resuelve al utilizar el filtro de Kalman como se muestra en el algoritmo 1 (Tarvainen et al, 2005; Haykin, 2001).

Por el contrario, si se desea reducir la interferencia de línea de potencia, es necesario combinar las representaciones dadas en (14b) y (16) para obtener el modelo conjunto de evolución de la dinámica del electrocardiograma y de la interferencia dados por:

Donde

y donde n=p-1, 0_mxn es una matriz de ceros de m filas y n columnas; δ_m(i) es la función delta, compuesta por un vector de longitud m de ceros con un uno en la posición i x(_pl) son los estados correspondientes a la interferencia sinusoidal; x(_ecg) son los estados correspondientes al proceso autorregresivo AR(p) que modela la señal sinusoidal; y v₁ y v₂ son secuencias de ruido blanco gaussiano independientes, alimentadas al oscilador y al proceso autorregresivo, respectivamente, a través de la matriz de acople de ruido B.

De acuerdo con lo anterior, la señal electrocardiográfica medida se genera mediante la siguiente ecuación:

donde x₁(_ecg)(k) corresponde a la señal ECG estimada con el modelo autorregresivo, x₁(_pl)(k) corresponde a la interferencia de línea de potencia estimada, C=[1 0 1 0 ... 0] es el vector de medición y n(k) es ruido de medición, el cual se supone es blanco gaussiano de media cero y varianza r_n.

Las ecuaciones (17), (18) y (19) constituyen la representación en espacio de estado de la evolución de la interferencia y del ECG, de los pesos con los que se ajusta al proceso autorregresivo para modelar al ECG y la señal ECG medida, y son utilizados para el esquema de estimación dual que permitirá separar la señal ECG de la interferencia sinusoidal, mediante el conjunto de ecuaciones que se muestra en el algoritmo 2 (Wan & Nelson, 1997).

En el algoritmo anterior se supone que las covarianzas correspondientes al estimador de parámetros, R_r y r_e son constantes, lo cual hace que el filtro en todo instante responda de la misma forma a cambios suaves o bruscos de la dinámica del ECG. Esto además acarrea el problema de obtener los valores adecuados para estas covarianzas, de forma que el filtro pueda responder de forma suave pero adecuada a los cambios en la dinámica de la señal. En este sentido, una mejor forma de manejar estas covarianzas es hacer que se comporten de forma adaptativa a la dinámica de la señal, haciendo que el filtro se acople a la no estacionariedad de la misma. En la siguiente sección se explican dos métodos para adaptar las covarianzas.

Otro aspecto importante en el algoritmo 2 es el cálculo del vector C_a, el cual acopla las correcciones a cada uno de los parámetros del filtro. Se calcula como sigue (Nelson, 2000; Haykin, 2001):

Se debe notar que el cálculo de esta derivada no es directo ya que, dada la estructura recurrente del sistema, es función de , y ambos son función de . Por lo tanto, el vector C_a debe ser calculado usando derivadas recurrentes con una rutina similar al aprendizaje recursivo en tiempo real, usado en redes neuronales. Tomando la derivada de las ecuaciones del filtro de señal, resulta el siguiente sistema de ecuaciones recursivas:

El último término en la ecuación (22) puede ser despreciado si se asume que la ganancia de Kalman K_x(k) es independiente del vector de parámetros a. Aunque esto simplifica bastante al algoritmo, el valor exacto de puede ser calculado; sin embargo, esto sólo se hace necesario cuando la señal es altamente no lineal o está perturbada por un nivel alto de ruido.

Así, se obtienen las siguientes ecuaciones para calcular C_a:

Donde

Adaptación de las covarianzas

Una parte importante para el acople del estimador de parámetros a la no estacionariedad de la señal, ya sea del conjunto de ecuaciones dadas en el algoritmo (1) o en el algoritmo (2), es hacer que las covarianzas R_r y r_e sean adaptativas.

La adaptabilidad de la matriz de covarianza R_r hace que los parámetros del modelo tengan mayor influencia de los datos recientes, mientras que la adaptabilidad de la varianza r_e permite estimar la varianza del ruido que es aplicado al proceso autorregresivo, lo cual tiene directa relación con la dispersión observada de los datos.

En (Tarvainen et al 2005), se propone estimar el ruido de observación r_e, iterativamente en cada paso del filtro de Kalman, como sigue:

donde γ≤1 es un coeficiente de adaptación y ε(k) es el error de predicción, definido como

para el estimador del algoritmo (2), o

para el estimador del algoritmo (1).

Por otra parte, la matriz de covarianza R_r se supone de la forma

donde el parámetro σ_r² se ajusta en cada paso de acuerdo con:

donde es la varianza estimada de la señal ECG en el instante k y α es un parámetro de actualización a través del cual la adaptabilidad del algoritmo puede ser ajustada. Se incluye la varianza de la señal ECG con el fin de remover la influencia de la amplitud de la señal en las estimaciones. De esta forma, no es necesario escoger un factor α independiente para señales ECG con diferentes niveles de variación. Se puede incrementar la adaptabilidad del algoritmo al incrementar α. La varianza de los parámetros estimados es, sin embargo, inversamente proporcional al valor de α y, por lo tanto, α debería ser especificado de tal forma que se alcance el balance deseado entre la adaptabilidad del filtro y la varianza de los parámetros.

Otra forma de realizar la adaptación de la matriz R_r es definiéndola como:

donde λ<1 se define como el factor de olvido del estimador. Como se demuestra en (Nelson, 2000), cuando se elige R_r(k) de esta forma se genera una ventana exponencial en el tiempo que sopesa los datos estimados, dándole mayor ponderación a las estimaciones de estado recientes con respecto a las antiguas, haciendo de esta forma que los datos antiguos tengan menor influencia en los parámetros estimados .

La introducción de este factor de olvido hace al filtro más flexible, voviéndolo más sensible a los datos nuevos, con estimados de señal más precisos. Nótese, sin embargo, que aunque el factor de olvido es implementado dentro del estimador de parámetros, a través de la covarianza del proceso R_r, esto no implica que el sistema en esencia sea variante en el tiempo.

Desarrollo experimental

A continuación se describen los procedimientos realizados junto con la descripción de la base de datos de señales ECG utilizada para probar los algoritmos descritos previamente.

Base de datos

Los experimentos realizados para probar los algoritmos se basan en la base de datos QT, la que incluye señales ECG escogidas como representativas de una amplia variedad de morfologías de los segmentos QRS y ST–T, con el fin de probar los algoritmos de detección de segmentos QT con variabilidad de señales del mundo real. Los registros son escogidos principalmente de las bases ECG existentes, incluyendo la base de datos MIT–BIH de arritmias, la base de datos de la Sociedad Europea de Cardiología ST–T, y otras bases de datos provenientes del Centro Médico Beth Israel Deaconess, de Boston.

Las bases de datos existentes son una excelente fuente de datos variados y bien caracterizados, para los cuales se han adicionado anotaciones de referencia que marcan la localización de los límites de las ondas. Los registros adicionales fueron escogidos para representar extremos en las patologías cardíacas. Además se han reunido datos de registros Holter de pacientes que han experimentado muerte súbita durante las grabaciones, y de pacientes seleccionados en edad y género sin una enfermedad cardíaca diagnosticada.

La base de datos QT contiene un total de 105 registros de 15 (min) de dos canales, seleccionados para que no contengan desviación de línea base u otros artefactos. Todos los registros fueron muestreados a 250 (Hz). Aquellos que no fueron muestreados originalmente con esta tasa fueron convertidos utilizando un software de submuestreo.

Dentro de cada registro fueron etiquetados manualmente, por cardiólogos, entre 30 y 100 latidos representativos, quienes identificaron el inicio, el pico y el fin de la onda P, el inicio y fin del segmento QRS (la marca fiducial del QRS, típicamente el pico de la onda R fue dada por un detector de QRS automático), el pico y fin de la onda T y, en caso de presentarse, el pico y fin de la onda U.

Para comparar las propiedades de adaptación de cada uno de los algoritmos con diferentes conjuntos de parámetros se utilizaron señales ECG tomadas de la base de datos QT (primeros 8(s) de sele0129 y sel230), mostradas en la Figura 1. El registro sele0129, mostrado en la Figura 1(a), contiene una señal ECG normal, mientras que el registro sel230,indicado en la Figura 1(b), contiene una señal que posee un cambio de morfología, correspondiente a una patología. El objetivo de tomar estos registros es comparar la dinámica de los parámetros estimados en los casos de un ECG normal y de un ECG patológico.

Pruebas de modelado adaptativo

Mediante la pruebas llevadas a cabo en esta sección se pretende estudiar cómo el algoritmo de estimación de parámetros para el modelo autorregresivo es afectado por los cambios en los parámetros de adaptación, es decir, ante las variaciones de α, γ y λ, para adaptar las matrices de covarianza. Esto se realiza con el objeto de tener un conjunto de parámetros adecuados para las pruebas consecutivas de filtrado.

Inicialmente, con el fin de conocer el orden óptimo para el modelo autorregresivo, se calcula el criterio de información de Akaike (Akaike information criterion: AIC) para cada una de las señales ECG de la base de datos QT; luego, el orden óptimo que se ajusta mejor a la base de datos. La curva promedio se muestra en la Figura 2. Se observa que la mayor dispersión se concentra al final de la misma y en un pequeño segmento de la curva al inicio. De acuerdo con esto, se puede decir que el orden óptimo para representar el ECG con un modelo AR es p=11.

Una vez obtenido el orden, se efctúan pruebas sobre las señales sele0129 y sel230 con el algoritmo de estimación de parámetros. En las Figuras 3 y 4 se registra la varianza de la medida estimada y el primer parámetro del modelo AR estimado para la señal sele0129 y sel230, respectivamente, con los parámetros (α=1, γ=0.95), (α=10^-4, γ=0.95), (α=1, γ=0.75) y (α=10^-4, γ=0.75).

Al comparar cada una de las curvas se nota que al disminuir el valor de α, la tasa de variación de los parámetros disminuye, reflejándose en la menor velocidad con la que varían los parámetros. De igual forma disminuye la amplitud máxima de la varianza de la medida estimada, conservando el tiempo de caída desde el valor máximo al valor mínimo. En el caso del parámetro γ se observa que, al disminuirlo, los parámetros estimados tienen menor dispersión, pero varían de forma rápida también. A su vez, vemos que la varianza aumenta un poco en magnitud y tiene un menor tiempo de caída desde el valor máximo al mínimo.

En la Figura 4, donde aparecen los resultados para la señal sel230, se observan las mismas tendencias en las señales al modificar los parámetros, con la diferencia que la dispersión en el tamaño de los picos es mucho mayor que en el caso anterior. A diferencia de lo que se podría esperar, no se observan cambios en los parámetros, ni en las varianzas, en los instantes en que la señal cambia de forma, a excepción de un pequeño cambio en el primer parámetro del modelo AR.

En la Figura 5 se muestran las varianzas y los parámetros estimados por el algoritmo modificado, utilizando como parámetros (λ=1, γ=0.95), (λ=0.995, γ=0.95), (λ=0.9, γ=0.75) y (α=0.995, γ=0.75) para cada una de las señales. En este caso se ve que a medida que el factor de olvido se aproxima a 1, los parámetros estimados se hacen más estables y al igual que lo hace la varianza estimada.

Finalmente, en la Figura 7 se muestra la señal ECG junto con la evolución de la varianza estimada r_e, definida en la ecuación (24), durante un latido.

Allí se puede observar la correspondencia entre los puntos fiduciales del complejo QRS y la señal de la evolución de la covarianza r_e. El inicio del complejo QRS se marca por el incremento abrupto de la covarianza, el pico R se encuentra en la mitad de los picos, mientras que el fin del complejo QRS se marca también por la disminución de la covarianza.

Pruebas del algoritmo de filtrado

Las pruebas realizadas con los algoritmos de filtrado siguen el esquema descrito en Figura 8:

que corresponde al esquema que se describe a continuación:

1. Obtención de la señal ECG limpia de la base de datos QT: se toman segmentos de 20(s) de registros de la base de datos QT que no posean interferencia de línea de potencia.

2. Adición de interferencia: las señales son contaminadas con interferencia de 60(Hz) para diferentes relaciones de señal a ruido (SNR), las cuales se definen como

donde x(n) es la señal original, es la señal contaminada y es la media de x(n).

3. Eliminación de perturbación: las señales son filtradas con el filtro de Kalman, basado en la representación en espacio de estado descrito por las ecuaciones (15b), o con el filtro de Kalman dual, basado en la representación en espacio de estado dada en las ecuaciones (17), (18) y (19), utilizando los dos tipos de adaptación de la covarianza.

4. Medición del desempeño: se compara el desempeño de los diferentes algoritmos midiendo el índice de correlación entre la señal original y la señal filtrada, que se define como

donde C_xy es la covarianza cruzada de x e y y C_xx es la covarianza de x.

En la Figura 9 se exhiben los resultados de filtrado en la base de datos QT. Se muestran secuencialmente los resultados para el filtro IIR notch, para el filtro de Kalman simple, y para el filtro de Kalman dual con adaptación basada en factor de olvido. De izquierda a derecha, la relación señal a ruido aumenta, por lo tanto, la potencia de la señal ECG se incrementaen relación con la de la interferencia.

Los resultados obtenidos reflejan un buen comportamiento de todos los filtros, dado que en general se observa que el índice de correlación es superior a . Comparando uno a uno los filtros, se aprecia que el más estable es el filtro IIR notch, dado que los índices de correlación resultantes se agrupan mejor en un punto más cercano a 1 que los demás; sin embargo, como se hará notar más adelante, estos resultados esconden las interferencias que genera este filtro al final del segmento QRS en algunas señales. De aquí se deduce también que el filtro de Kalman dual tiene un desempeño superior al del filtro de Kalman simple.

En la Figura 10 se registran los resultados de los algoritmos de filtrado en las señales que habían sido tomadas como referencia. En la parte superior se muestra la señal original; en la parte media, la señal con ruido de 60 Hz; y en la parte inferior, la señal filtrada con cada uno de los algoritmos, donde KF representa al filtro de Kalman simple, DKF1 al filtro de Kalman dual con adaptación basada en el ajuste con α, y DKF2 al filtro de Kalman dual con adaptación basada en factor de olvido.

En la Figura 11 se compara el latido de una señal ECG limpia (señal de la parte superior) con las señales filtradas mediante un filtro notch IIR (parte media) y el método desarrollado (parte inferior). En este caso se indicael resultado del filtro de Kalman dual con adaptación mediante factor de olvido. Se aprecia que mediante el filtro notch IIR se obtiene una distorsión al final del complejo QRS, mientras que con el filtro de Kalman se obtiene una señal sin distorsiones.

Discusión y Conclusiones

En este artículo se han presentado dos métodos para la remoción de la interferencia de línea de potencia basados en el filtro de Kalman. El primero consiste en la implementación de un estimador de estado para un modelo consistente de un oscilador sinusoidal y un modelo autorregresivo de primer orden. El segundo consiste en la implementación de un estimador de parámetros y de estados, mediante una estrategia de estimación dual, para un modelo construido con un oscilador sinusoidal y un modelo autorregresivo de orden 11. Para este estimador dual se plantean dos estrategias de adaptación, de las cuales se evidenciaque se comportan de forma similar y que no presentan diferencias para el proceso de filtrado.

La estrategia de modelado adaptativo de la señal ECG mediante un proceso autorregresivo de orden 11 para el proceso de remoción de perturbaciones permite obtener resultados mejores que los obtenidos mediante un modelo autorregresivo invariante de primer orden. A lo largo del artículo esto se demuestra de dos formas: primero, observando la curva promedio del criterio de información de Akaike, mostrada en Figura 2, donde es clara la diferencia que hay entre un modelo de orden 1 y un modelo de orden 11; segundo en la prueba realizada sobre la base de datos QT, cuyos resultados se observan en Figura 4, donde los resultados obtenidos con el filtro de Kalman tienen mayor dispersión y menor mediana que los obtenidos con el filtro de Kalman dual.

El filtro de Kalman tiene como gran ventaja que no genera distorsiones en las señales procesadas, lo cual depende enteramente de la elección de los parámetros del filtro y del modelo de adaptación, no obstante los resultados son inestables ante ciertos valores aleatorios de arranque en comparación con el filtro IIR, lo cual se manifiesta en la mayor dispersión y en la presencia de outliers en los boxplots.

El ajuste de los parámetros del filtro de Kalman es uno de los puntos vitales de esta metodología y depende totalmente de las condiciones a las que se enfrentará. En el caso en el que el filtro se enfrente a una interferencia con parámetros poco variables, se utilizan parámetros de adaptación más laxos, mientras que en el caso que la interferencia varíe mucho, el proceso de adaptación debe ser más rígido. Las covarianzas del ruido del proceso y del ruido de medición R_n y R_v determinan la forma como el filtro aproxima la señal ECG y rastrea las variaciones de la interferencia sinusoidal. El uso de un ruido de medición mayor al ruido del proceso hace que el filtro reaccione ante cambios en la forma de la señal cardíaca como otro tipo de perturbación, lo cual provoca distorsión en la señal procesada.

Los parámetros y la varianza r_e calculados por el estimador de parámetros se encuentran relacionados directamente con la evolución de la señal cardíaca. Sobre todo, la evolución de la varianza r_e, muestra una estrecha relación con los eventos del complejo QRS, lo cual puede ser útil en la obtención de los puntos fiduciales, la segmentación y posiblemente en la clasificación de las señales ECG.

Como parte del trabajo futuro se plantea seguir estudiando las posibilidades que ofrece el filtro de Kalman dual en el procesamiento de señales ECG, entre ellos, la remoción de perturbación de línea de potencia y de línea base, estimando de forma la frecuencia de la interferencia, la extensión de este método a procesos de ubicación de puntos fiduciales, mediante la estimación de la varianza r_e y de los parámetros del proceso autorregresivo.

Agradecimientos

Este trabajo se realiza en el marco del proyecto de la DIMA Técnicas de Computación de Alto Rendimiento en la Interpretación Automatizada de Imágenes Médicas y Bioseñales.

Bibliografía

Ahlstrom, M., Tompkins, J., Digital filters for real-time ECG signal processing using microprocessors., IEEE Transactions on Biomedical Engineering, 32, 9, 1985, pp. 708–713.         [ Links ]

Aström, K. J., Wittenmark, B., Computer-Controlled Systems: Theory and Design, third edn., Prentice Hall Information and System Sciences, 1997.         [ Links ]

Dotsinsky, I., Daskalov, I., Accuracy of 50(Hz) interference subtraction from an electrocardiogram., Medicine, Biology, Engineering and Computing, 34, 1996, pp. 489–494.         [ Links ]

Er, M. H., Designing notch filter with controlled null width., IEEE Signal Processing, 24, 1991, pp. 319–329.         [ Links ]

Hamilton, P., A comparison of adaptive and nonadaptive filters for the reduction of powerline interference in the ECG., IEEE Transactions on Biomedical Engineering, 43, 1996, 105–109.         [ Links ]

Haykin, S., Kalman Filtering and Neural Networks., first edn, Wiley Interscience, 2001.         [ Links ]

Huhta, J., Webster, J., 60(Hz) interference in electrocardiography., IEEE Transactions on Biomedical Engineering, 20, 1973, pp. 12–28.         [ Links ]

Jang, Y. K., Chicharo, F., Adaptive IIR comb filter for harmonic signal cancellation., Int. J. Electronics, 75, 1993, pp. 241–250.         [ Links ]

Levkov, C., Michov, G., Ivanov, R., Daskalov, I., Subtraction of 50(HZ) interference form the electrocardiogram., Medicine, Biology, Engineering and Computing, 22, 1984, pp. 371–373.         [ Links ]

Middleton, R. H., Goodwin, G., Improved finite word length characteristics in digital control using delta operators, IEEE AC-31, 1986.         [ Links ]

Middleton, R. H., Goodwin, G., Digital Control and Estimation: A unified Approach, first edn., Prentice-Hall International, Inc., 1990.         [ Links ]

Mneimneh, M., Yaz, E., Johnson, M., Povinelli, R., An adaptive Kalman filter for removing baseline wandering in ECG signals., Computers in Cardiology, 33, 2006, pp. 253–256.         [ Links ]

Nelson, A., Nonlinear estimation of noisy time–series by dual Kalman filtering methods., PhD thesis, Oregon Graduate Institute of Science and Technology, 2000.         [ Links ]

Pei, S., Tseng, C., Elimination of AC interference in electrocardiogram using IIR notch filter with transient suppression., IEEE Transactions on Biomedical Engineering, 42, 1995, pp. 1128–1132.         [ Links ]

Sameni, R., Shamsollahi, M., Jutten, C., Multi-channel electrocardiogram denoising using a bayesian filtering framework., Computers in Cardiology, 33, 2006, pp. 185–188.         [ Links ]

Sörnmo, L., Laguna, P., Bioelectrical Signal Processing in Cardiac and Neurological Applications., Elsevier Academic Press, 2005.         [ Links ]

Tarvainen, M., Georgiadis, S., Karjalainen, P., Time-varying analysis of heart rate variability with Kalman smoother algorithm., 27^th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2005.         [ Links ]

Wan, E., Nelson, A., Dual Kalman filtering methods for nonlinear prediction, estimation and smoothing., Advances in Neural Information Processing Systems, 1997.         [ Links ]

Ziarani, A., Konrad, A., A nonlinear adaptive method of elimination of power line interference in ECG signals., IEEE Transactions on Biomedical Engineering, 49, 6, 2002, pp. 540–547.        [ Links ]