REPRESENTACIONES DE LA SEÑAL DE VARIABILIDAD

El estudio de la HRV se basa en el análisis de la serie de ocurrencias temporales, que denotaremos por tk producidos por los complejos QRS, es decir, se basa en el análisis de las variaciones del ritmo cardíaco. Un estudio más adecuado implicaría analizar los instantes de ocurrencia de la onda P, dado que el objetivo principal es el análisis de la influencia de la actividad del SNA sobre el NSA, cuyo reflejo eléctrico en el ECG es la onda P. El principal inconveniente de este enfoque es la baja amplitud de dicha onda, lo que no permite una detección precisa. Por este motivo se utiliza como punto de ocurrencia del latido la onda R del complejo QRS. Esta elección es ampliamente aceptada debido a que el intervalo PR se puede considerar fijo.

Un objetivo importante en el estudio de la HRV es encontrar una representación del ritmo cardíaco que refleje con precisión sus variaciones. Hay que tener también en cuenta el hecho de que la secuencia de ocurrencias de los latidos tk constituye un proceso no uniformemente muestreado; esta es la razón por la que existe la variabilidad del ritmo cardíaco. En ocasiones será necesario remuestrar uniformemente la representación del ritmo cardíaco para que esta sea compatible con las técnicas de análisis, que requieren señales uniformemente muestreadas, como es el caso de algunas técnicas espectrales.

Debe quedar clara la distinción entre el muestreado de la representación del ritmo cardíaco y el muestreado de digitalización de la señal de ECG. El muestreado del ECG tiene que ver con la resolución de la propia señal de ECG, y por lo tanto con la resolución de los instantes temporales de los complejos de QRS. Las frecuencias de muestreo habituales para la señal de ECG se encuentran en el rango 500-1000 Hz, de forma que permitan una detección de complejos QRS exacta. El muestreo del ritmo cardíaco es mucho menor, sólo unos pocos HZ; teniendo en cuenta un ritmo cardíaco de 60 latidos por minuto, eso implica un latido cada segundo y por tanto una frecuencia de muestreo de 1 Hz.

El requerimiento de precisión de la representación elegida implica algún mecanismo para evaluar el comportamiento de las distintas representaciones. Un modelo matemático adecuado para comparar las distintas representaciones es el modelo Integral Pulse Frequency Modulation, IPFM, que reproduce convenientemente la influencia del SNA sobre el ritmo cardíaco, es decir, sobre el NSA.

Figura 1: Modelo IPFM.

Modelo IPFM

El modelo IPFM se utiliza para producir una serie de eventos, como puede ser la serie de ocurrencias de latidos cardíacos. Se asume que el modelo tiene una señal de entrada de tiempo continuo, que posee un significado fisiológico determinado.
El esquema de bloques del modelo se puede observar en la Figura 1. La señal de entrada del modelo se integra hasta que alcanza un umbral, R=1, en ese instante de tiempo, tk, se produce un evento y el integrador es reiniciado. El proceso se repite continuadamente. R representa la longitud de intervalo media entre sucesivos eventos. La señal de entrada del modelo posee dos componentes, m0=HR, que es una constante y representa el nivel de componente continua, y m(t), que es la señal moduladora, y que no posee componente continua. Un requisito de la señal de entrada del modelo es que sea siempre positiva, por lo que la señal moduladora debe cumplir que jm(t)j << m0. En la Figura 2 se puede observar el comportamiento del modelo para una señal de entrada sintética.

El objetivo del análisis de HRV es recuperar la información que transporta la señal m(t) a partir de la serie de eventos que es la única accesible, y que denotaremos (t), donde:

Desde el punto de vista fisiológico m0 representa el ritmo cardíaco medio, m(t) representa la variabilidad del ritmo cardíaco, provocada por la influencia del SNA sobre el NSA. Normalmente m(t) es de banda limitada, de forma que las componentes superiores a 0.5 Hz pueden ser despreciadas. El requisito jm(t)j << m0 representa el hecho de que la HRV es pequeña comparada con el ritmo cardíaco medio. El comportamiento del modelo IPFM se puede expresar como:

donde tk-1 y tk son los instantes de ocurrencia de dos latidos consecutivos. Para entender mejor el sentido de la señal m(t) como modulación del ritmo, es decir, como determinante de la variación entre los intervalos entre eventos, supongamos que m(t) _ 0, de esta forma:

es decir el intervalo temporal entre eventos, es decir, la ocurrencia de los latidos y por tanto el ritmo cardíaco, sería fijo. Dicho de otra forma no habría variación del ritmo (m(t) _ 0). La frecuencia de repetición de los intervalos sería:

Usualmente m0 se escoge con valor 1, de forma que R = TI representa el valor medio de los intervalos RR. De esta forma si se escoge como umbral R = 1, el ritmo cardíaco que se conseguiría con el modelo IPFM sería de 60 latidos por minuto. Si se asume que el primero latido ocurre en el instante de tiempo t0, se puede escribir la ecuación del modelo IPFM como:

Figura 2: Ejemplo de funcionamiento del modelo IPFM. La señal superior muestra la entrada al modelo, en el panel medio se muestra la integral acumulada así como los instantes en los que se alcanza el valor de umbral R representado como crculos. El panel inferior muestra la serie de eventos que representan los instantes de ocurrencia de los latidos.

Señales de variabilidad

Las señales de variabilidad persiguen el objetivo de reejar de forma precisa las variaciones en la frecuencia cardíaca. Las distintas señales de variabilidad deben basarse en los instantes de ocurrencia de los complejos QRS, en concreto del pico de la onda R si se toma esta como referencia. Denotamos los distintos instantes de latidos como t[n].

Señal de Período Cardíaco

Se denomina Señal de Período Cardíaco, denotada por PC[n], a la serie temporal discreta en la que, cada muestra, es el intervalo de tiempo transcurrido entre dos latidos consecutivos, es decir:

Esta señal se conoce también por el nombre de tacograma de intervalos o interval tachogram, también es conocida como Heart Period, HP ver Figura 1

Señal de Ritmo Cardíaco

Se denomina Señal de Ritmo Cardíaco, denotada por RC[n], a la serie temporal discreta en la que, cada muestra, es el inverso del intervalo entre dos latidos consecutivos, es decir

En la literatura también es conocida por Heart Rate, HR.

Señal de Temporización Cardíaca

Se denomina Señal de Temporización Cardíaca, denotada por TC[n], a la serie temporal discreta que se construye como:

donde T es el período cardíaco medio. En la literatura también es conocida por Heart Timming, HT.

MÉTODOS TEMPORALES ESTADÍSTICOS

Las medidas estadísticas consisten, en general, en medidas de dispersión (desviaciones típicas) de los intervalos NN o de sus diferencias.

Estas medidas se pueden obtener sobre registros de larga duración (habitualmente 24 horas), o sobre segmentos de menor longitud (generalmente periodos 5 minutos).
Una consideración importante es si la medida refleja la variabilidad del ritmo cardíaco a largo plazo o a corto plazo, de forma que proporciona principalmente información referente a la actividad parasimpática o a la actividad simpática.

Las medidas basadas en la desviación típida de los intervalos NN en registros de 24 horas reflejan las variaciones de largo plazo del ritmo cardíaco. Además, se han propuesto otras medidas de dispersión basadas en las diferencias de intervalos NN sucesivos, que reflejan principalmente las variaciones de corto plazo del ritmo cardíaco, debido a que el efecto de la operación de diferenciación es acentuar el contenido de alta frecuencia.
Los índices más comunes en el análisis con métodos estadísticos temporales de la señal de HRV se describen brevemente a continuación.

AVNN

Descripción: valor promedio de todos los intervalos NN. Puede obtenerse en registros cortos de 5 minutos, o largos de 24 horas.
Metodología:

donde NNi representa el valor de cada uno de los intervalos NN, y K es el número total de intervalos.

SDNN

Descripción: desviación estándar de los intervalos NN del registro. La varianza es matemáticamente igual a la potencia total del espectro, por lo tanto, esta medida refleja todas las componentes cíclicas responsables de la variabilidad del periodo de registro.
Metodología:

donde NNi representa el valor de cada uno de los intervalos NN, AV NN es la media de todos ellos y K es el número total de intervalos NN.

Propiedades: esta medida calculada sobre periodos de 24 horas sólo proporciona una caracterización aproximada de la señal de HRV, ya que la media del ritmo cardíaco varia considerablemente de las partes activas del día, a la noche, de esta forma refleja principalmente la componente de baja frecuencia. No se trata de una medida robusta y por lo tanto, para poder realizar comparaciones los registros de datos deben tener la misma longitud.

SDANN

Descripción: desviación estándar de las medias de los intervalos NN calculadas sobre periodos de 5 min.
Metodología:

donde AV NN5i representa cada una de las medias de los intervalos de 5 minutos, AV NN es el promedio de todos los NN, y k es el número de intervalos de 5 minutos de los que se dispone.

Propiedades: esta medida refleja principalmente las variaciones lentas del ritmo cardíaco debido a las medias realizadas sobre los periodos de 5 minutos.

SDNNIDX

Descripción: media de la desviación estándar de los intervalos NN de cada segmento de 5 minutos del registro completo.

Metodología:

donde SDNN5i es cada una de las desviaciones estándar de los intervalos de 5 minutos y k es el número de intervalos de 5 minutos.

RMSSD

Descripción: raíz cuadrada de la media de las diferencias al cuadrado entre sucesivos intervalos NN

donde NNi representa a cada uno de los intervalos NN, siendo NNi_1 el intervalo NN anterior y K el número total de intervalos.

Propiedades: al ser una medida basada en la diferencia de intervalos sucesivos, la operación de diferencia acentúa las altas frecuencias contenidas en las series NN.

NN50

Descripción: número de pares de intervalos NN adyacentes que difieren más de 50 ms en el registro total.

Metodología:

donde NNi representa el intervalo NN actual, NNi_1 el intervalo NN anterior, y K es el número total de intervalos.

Propiedades: esta medida no proporciona información tan detallada de las componentes de alta frecuencia como las dos anteriores, pero es menos vulnerable frente a la presencia de artefactos.

pNN50

Descripción: porcentaje de intervalos NN adyacentes que distan más de 50 ms en el registro total.

Metodología:

donde K es el número total de intervalos NN.

Propiedades: la medida pNN50 pertenece a la de estadísticos pNNx, donde x indica umbrales distintos a 50ms que en algunos proporcionan mejores resultados como discriminantes de edad, periodos de sueño o consciencia y salud o enfermedad.

Índice de Cambios de Aceleraciones

Definición: el Indice de Cambios de Aceleraciones (Acceleration Change Index, ACI) es una medida que cuántica el signo de las diferencias de las series de intervalos RR.
Metodología: para obtener el ACI se siguen los siguientes pasos:

• Se obtienen las diferencias de la serie temporal RR:

donde RR(n) es el intervalo entre el latido n y el latido n + 1, y N es el número total deintervalos NN.

• Se define SDRR como el signo de la serie DRR:

Desde n = 1 hasta n = N _1, se genera la serie de cambio de signo (SC) de la siguiente forma:

donde M es el n_umero de cambios de signo. SC(j) = n implica que RR(n) es el j-ésimo
máximo o mínimo local del tacograma.

• Se define DSC como la diferencia entre dos valores de SC consecutivos, para obtener la distancia (en latidos):

• Finalmente el ACI se calcula como la siguiente proporción:

ACI = k/M
donde k representa el número de veces que la serie temporal DSC es igual a 1, es decir, cuenta el número de veces que dos cambios de signo consecutivos son iguales, y M es el número total de muestras total de la serie DSC.

Pruebas

En los siguientes apartados se presentan diferentes pruebas sintéticas para cada una de las medidas estadísticas presentadas en este capítulo, con señales simples cuyo resultado es conocido o fácil de calcular, con el objetivo de verificar la implementación de dichas medidas. Las señales sintéticas empleadas en estas pruebas se detallan a continuación:

Señal ejemplo 1: señal determinista sinusoidal cuyos parametros media y desviacion estandar pueden ser modicados

donde A representa el valor medio y B
raiz de 2=2 es el valor de la desviacion estandar de la señal.

Señal ejemplo 2: señal aleatoria con funcion de densidad de probabilidad gaussiana

donde se representan la media de la señal y la desviacion estandar, estos seran los parametros libres para las pruebas

Señal ejemplo 3: señal suma de las señales ejemplo 1 y 2.

Señal ejemplo 4: señal aleatoria con función de densidad de probabilidad uniforme,

cuyo valor medio y desviación estándar serán los parámetros libres que utilizaremos en las pruebas.

Señal ejemplo 5: mapa logístico, señal caótica determinista que puede ser expresada como una ecuación recursiva

donde r es un parámetro que varia entre 0 y 4 y el valor inicial de x está restringido entre 0 y 1. Para valores de r entre 0 y 3, la ecuación no linear puede tener una solución que se aproxima a un determinado estado y permanece en él (punto fijo en el espacio de estados). Para r entre 3 y 3.5, la solución puede presentar ciclos y oscilaciones entre diferentes valores. A medida que r aumenta aparecen varias bifurcaciones y cuando r alcanza el valor crítico (r = 3:569944:::) la serie temporal resultante tiene un comportamiento caótico.

Señal ejemplo 6: el movimiento Browniano fraccionario (fractional Brownian motion, fBm) es un proceso aleatorio gaussiano no estacionario de media cero. Una de sus principales propiedades es la autosemejanza. El fBm está caracterizado por el parámetro H, exponente de Hurst. Mandelbrot demostró en que la diferenciación del fBm da como resultado un proceso estacionario gaussiano autosemejante conocido por ruido fraccionario Gaussiano (fractional gaussian noise, fGn). A este tipo de procesos autsemejantes se les puede asociar un espectro en ley de potencia, 1=f_, donde _, en el caso del fBm se relaciona con el exponente de Hurst como B = 2H + 1 y en el caso del fGm, como B = 2H _ 1

Señal ejemplo 7: señal real de intervalos NN obtenida de un registro Holter de 24 horas.

Pruebas AVNN

Prueba señal ejemplo 1: en esta prueba se varía el valor medio de la señal de prueba y se comprobamos en cada caso que el índice AV NN calculado coincide con el valor medio de la señal. En la Figura Fig 3 (a), se presenta un ejemplo de una señal con media 4:5, en la que se presenta también el valor obtenido del índice AV NN. A continuación, se han escogido 10 valores para la media y para cada uno los valores se ha obtenido una realización de la señal ejemplo 1.

En la Figura 3 (b) se representa el índice AV NN obtenido para cada uno de los 10 valores, como era de esperar, los valores del índice coinciden con los valores medios de las señales.

Prueba señal ejemplo 2: siguiendo la estructura de la prueba anterior, en la Fig 3 (c) se representa una señal de media m = 4:5 y el valor AV NN obtenido para la misma. A continuación, se escogen 10 valores para la media y para cada valor de la media se obtienen 10 realizaciones de la señal ejemplo 2, cada una de ellas con 1000 muestras. Para cada valor de la media se representa la media y la desviación estándar del índice AV NN, superpuesta a la recta de valores teóricos (Figura 3 (d)). Se puede observar que los valores medios del índice AV NN se ajustan al valor esperado con un error muy pequeño.

Prueba señal ejemplo 3: el valor medio de esta señal prueba es la suma de las medias de las señales ejemplo 1 y ejemplo 2. Por lo tanto, el índice AV NN de la señal suma será igual a la suma de los índices AV NN de estas dos señales. En la Figura 3 (e), se presenta el resultado de la suma de una señal ejemplo 1 de media 2 y una señal ejemplo 2 de media 3, por lo tanto, se obtiene un índice AV NN = 5.

Prueba señal ejemplo 7: Finalmente en la Figura 3 (f) se presenta el índice AV NN obtenido en una señal real de HRV.

Figura 3: Pruebas con AVNN: (a) calculo del AVNN en una señal sinusoidal, la linea roja representa el valor medio de la señal. (b) Calculo del AVNN para distintos valores de la media de una señal sinusoidal.(c) Calculo del AVNN en una señal aleatoria gaussiana. (d) Media y desviacion estandar del AVNN para distintos valores de la media de una señal aleatoria gaussiana, la linea roja representa el valor teórico esperado. (e) Calculo del AVNN en una se~nal suma de señales. (f ) Ejemplo del calculo del AVNN en una señal real.

Pruebas SDNN

Prueba señal ejemplo 1: en esta prueba, el parámetro libre es la desviación estándar de la señal, por lo tanto, se varía el valor de este parámetro y se comprueba, en cada caso, que coincide con el del índice SDNN. En la Figura 4 (a), se presenta un ejemplo de una señal con desviación estándar 1, así como el valor obtenido para el índice SDNN, que coincide con la desviación estándar de la señal ejemplo.
A continuación, se escogen 10 valores de la desviación estándar y para cada uno de estos valores se obtiene una realización de la señal ejemplo 1. En la Figura 4 (b) se representa el índice SDNN obtenido para cada uno de los 10 valores, y como era de esperar, los valores del índice obtenidos coinciden con las desviaciones estándar de las señales.

Prueba señal ejemplo 2: en la Figura 4 (c) se representa una señal ejemplo 1 con desviación estándar 1, también se representa el valor del índice SDNN obtenido para la misma. A continuación, se toman 10 valores distintos para la desviación estándar, y para cada uno de ellos, se obtienen 10 realizaciones de la señal ejemplo 2 con 1000 muestras cada una.

En la Figura 4 (d) se representa la media y la desviación estándar del índice SDNN obtenido para cada valor de la desviación estándar, superpuesto a la recta teórica. Los valores medios del índice SDNN se ajustan al valor esperado, además, se puede observar que la desviación estándar de la estimación del índice aumenta con la desviación estándar de la señal.

Prueba señal ejemplo 7: en la Figura 4 (e) se muestral el valor del índice SDNN obtenido en una señal real de HRV.

Figura 4: Pruebas con SDNN: (a) calculo del SDNN en una señal sinusoidal, las lineas rojas representan la desviacion estandar de la señal. (b) Calculo del SDNN para distintos valores de desviacion estandar de una señal sinusoidal. (c) Calculo del SDNN en una señal aleatoria gaussiana. (d) Media y desviacion estandar del SDNN para distintos valores de la desviacion estandar de una señal aleatoria gaussiana, la
linea roja representa el valor teorico esperado. (e) Ejemplo del calculo del SDNN en una señal real.

Pruebas SDANN

Prueba señal ejemplo 1: en estas pruebas, los parámetros libres son tanto la media como la desviación estándar de la señal ejemplo. El índice SDANN debe ser igual a cero para cualquier valor estos parámetros, ya que en una señal sinusoidal el valor medio es constante para cualquier fragmento de la señal, por lo tanto, la desviación estándar de los valores medios de los intervalo de 5 minutos del registro total será cero para cualquier valor de los parámetros de la señal.
En la Figura 5 (a) se presenta un ejemplo de una señal con media 0:5 y desviación estándar 1 en la que se comprueba que el índice SDANN es cero.

Prueba señal ejemplo 2: para la señal ejemplo 2 el índice SDANN debe ser también cero, ya que, para una señal aleatoria con distribución gaussiana el valor medio es también idéntico para cualquier fragmento de la señal(Fig. 5 (b)).

Prueba señal ejemplo 3: el índice SDANN debe ser cero tambien para la señal ejemplo3 para cualquier combinacion de los parametros libres de la señal, ya que, el valor medio de la señal resultante es la suma de las medias de las dos señales individuales, el cual es un valor constante para todos los fragmentos de la señal suma, es decir, de la señal ejemplo 3 (Fig 5 (c)).

Prueba señal ejemplo 7: por ultimo, en la Figura 5 (d) se presenta el índice SDANN para una señal real de HRV. En este caso, el valor medio no es constante para todos los intervalos de 5 minutos del registro total, por lo que obtenemos un ndice SDANN diferente de cero.

Figura 5: Pruebas con SDANN: (a) calculo del SDANN en una señal sinusoidal, las lineas rojas representan la media en el centro y la desviacion estandar de la señal en los extremos. (b) Calculo del SDANN en una señal aleatoria gaussiana. (c) Indice SDANN de una señal suma. (d) Ejemplo del calculo del SDANN en una señal real.

Pruebas SDNNIDX

Prueba señal ejemplo 1: para la señal ejemplo 1 el valor del SDNNIDX debe ser igual al valor de la desviación estándar de la señal, ya que en una señal sinusoidal el valor de la desviación estándar es el mismo para todos los fragmentos de 5 minutos del registro total de la señal, por lo tanto la media de estos valores será igual al valor de la desviación estándar de la señal. En la Figura 6 (a) se presenta un ejemplo de una señal con desviación estándar 5 y por lo tanto, con índice SDNNIDX = 5.

Prueba señal ejemplo 2: el valor del índice SDNNIDX obtenido es igual al valor de la desviación estándar de la señal ejemplo 2, ya que en una señal aleatoria con distribución gaussiana, el valor de la desviación estándar es también idéntico para cualquier fragmento de la señal (Figura 6 (b)).

Prueba señal ejemplo 7: En la Figura 6 (c) se presenta el índice SDNNIDX para una señal real HRV. El índice SDNNIDX es diferente del valor de la desviación estándar de la señal, ya que los valores de la desviación estándar en los fragmento de 5 minutos del registro Holter total (aproximadamente 24 horas) no son iguales.

Figura 6: Pruebas con SDNNIDX: (a) Calculo del SDNNIDX en una señal sinusoidal, las lineas rojas representan la media en el centro y la desviacion estandar de la señal en los extremos. (b) Calculo del SDNNIDX en una señal aleatoria gaussiana. (c) Ejemplo del calculo del SDNNIDX en una señal real.

Pruebas RMSSD

Prueba señal ejemplo 1: para esta prueba se obtiene diez realizaciones de la señal ejemplo 1, cada una con un valor diferente de la desviación estándar, y se calcula el índice RMSSD para cada una de las realizaciones. Se cumple que el valor del índice RMSSD varía de forma directamente proporcional con la desviación estandar de la señal (Figura 7 (a)).

Prueba señal ejemplo 2: en esta prueba se generan 10 realizaciones de la señal ejemplo 2 para cada valor de la desviación estándar y se calculan la media y la desviación estándar de losíndices RMSSD obtenidos. En la Figura 7 (b) se observa que lo valores medios de los índices se aproximan a la recta que representa una variación lineal con la desviación estándar, se observa también que el error aumenta a medida que aumenta la desviación estándar de la señal.

Prueba señal ejemplo 3: para esta pruebas, se generan 10 realizaciones de la señal ejemplo 3 para 10 valores diferentes de la desviación estándar. En la Figura 7 (c) se presentan la media y la desviación estándar de los índices RMSSD obtenidos para cada realización, se observa que el índice RMSSD aumenta cuando la desviación estándar de la señal aumenta.

Prueba señal ejemplo 7: por último, calculamos el índice RMSSD para 2 señales reales diferentes; la primera perteneciente a un paciente sano y la segunda a un paciente con insuficiencia cardíaca. Se observa (Figura 7 (d) y (e)) que el índice RMSSD es menor para el caso patológico que para el caso normal.

Figura 7: Pruebas con RMSSD: (a) calculo del RMSSD en señales sinusoidales con diferentes valores de la desviaciones estandar. (b) Calculo del RMSSD en señales aleatorias gaussianas con diferentes valores de la desviaciones estandar, para cada valor de la desviacion estandar se han obtenido 10 realizaciones de la señal, se presenta la media y desviacion estandar de los índices RMSSD. (c) Calculo del RMSSD en
señales suma con diferentes valores de la desviaciones estandar, para cada valor de la desviacion estandar se han obtenido 10 realizaciones de la señal, se presenta la media y desviacion estandar de los índices RMSSD. (d) y (e) Ejemplo del calculo del RMSSD en señales reales, normal y patologica respectivamente.

Pruebas NN50

Prueba (a) señal ejemplo 2: para esta primera prueba del índice NN50, se toman 10 valores diferentes del parámetro desviación estándar, y para cada uno de estos valores, se han obtenido 10 realizaciones de la señal ejemplo 2, cada una de ellas con 1000 muestras. En la Figura 8 (a), se representan la media y la desviación estándar del índice NN50 obtenidos para cada uno de los 10 valores. Se puede observar que cuanto mayor es el valor de la desviación estándar de la señal mayor es el número de pares de intervalos NN adyacentes que distan más de 50 ms, y por lo tanto mayor es el valor del índice NN50.

Prueba señal ejemplo 7: en esta prueba, se calcula el valor del índice NN50 de una señal ejemplo 7, es decir, de una señal real de HRV. Se observa (Figura 8 (b)) que el índice NN50 es mucho mayor en este caso que los valores obtenidos en la prueba anterior, esto se debe a que la desviación estándar de la señal real es mayor que las utilizadas para las señales aleatorias gaussianas y sobre todo, a que la longitud de esta señal es mucho mayor, dado que, cuanto mayor sea la longitud de la señal, existen más oportunidades de tener pares de intervalos NN adyacentes que disten en más de 50 ms.

Prueba (b) señal ejemplo 2 : de la prueba anterior, se puede sospechar el índice NN50 es dependiente de la longitud de los datos, para comprobarlo se escogen 10 valores de la longitud de los datos diferentes para un valor de la desviación estándar fijo, y se obtienen una vez más 10 realizaciones de la señal ejemplo 2 para cada uno de estos valores.
En la Figura 8 (c) se representan la media y desviación estándar de los índices NN50 obtenidos. Se puede observar que el valor del índice crece con la longitud de los datos. Por lo tanto, cuando se realicen comparaciones con este índice deberían hacerse sobre registros de igual tamaño.

Figura 8: Pruebas con NN50: (a) Calculo del NN50 en señales aleatorias gaussianas para diferentes valores de la desviacion estandar. (b) Calculo del NN50 en una se~nal real, las lineas rojas representan, la media en el centro y la desviacion estandar en los extremos. (c) Calculo del NN50 en señales aleatorias gaussianas para diferentes valores de la longitud de los datos.

Pruebas pNN50

Prueba (a) señal ejemplo 2: en esta primera prueba del índice pNN50, se escogen 10 valores de la desviación estándar, y para cada uno de estos valores, se obtienen 10 realizaciones de la señal ejemplo 1, cada una de ellas con 1000 muestras. En la Figura 9 (a) se representan la media y la desviación estándar del índice pNN50 obtenidos para cada uno de los 10 valores. Se observa un comportamiento similar al obtenido para el índice NN50, pero en este caso el resultado es ofrecido en forma de porcentaje.

Prueba señal ejemplo 7: calculamos ahora el índice pNN50 para dos señales ejemplo 7 diferentes, una de ellas corresponde a un paciente con insuficiencia cardíaca y la otra corresponde aun paciente normal. El registro del paciente patológico posee mayor número de muestras y menor desviación estándar que el correspondiente al paciente normal (Figura 9 (c) y (d)). A pesar de que la señal correspondiente al registro patológico tiene un mayor número de muestras, su índice pNN50 es mucho menor que el del registro normal debido a que la desviación estándar notablemente menor.

Prueba (b) señal ejemplo 2: en esta prueba, se estudia el comportamiento de este índice respecto a la longitud de los datos, para ello se repite la prueba (b) señal. La Figura 9 (b) muestra que en este caso que el índice pNN50 es aproximadamente independiente la longitud de los datos, como era de esperar, debido a la normalización respecto a la longitud de los datos que se lleva a cabo en este índice.

Figura 9: Pruebas con pNN50: (a) calculo del pNN50 en señales aleatorias gaussianas para diferentes valores de la desviacion estandar. (b) Calculo del pNN50 en señales aleatorias gaussianas para diferentes valores de la longitud de los datos. (c) y (d) Calculo del pNN50 para dos señales reales normal y patologica respectivamente, con diferente desviacion estandar y diferente numero de muestras, las lineas rojas
representan, la media en el centro y la desviacion estandar en los extremos.

Pruebas ACI

Prueba señal ejemplo 2 y señal ejemplo 4: para entender la información proporcionada por el ACI, se comienza por estudiar su comportamiento en variables aleatorias incorrelacionadas. Se puede calcular el valor ACI teórico de una serie temporal mediante los coeficientes de su función de autocorrelación

donde N representa el numero de muestras, (k) y d(k) representan las funciones de autocorrelacion de las series temporales y de las series temporales diferenciadas respectivamente, para el retardo k. Esta ecuacion permite estudiar el comportamiento del ACI para cualquier serie temporal cuyos tres primeros coecientes de su funcion de autocorrelacion sean conocidos.
En el caso de variables aleatorias incorrelacionadas, la ecuacion queda como sigue,

Para verificar la ecuación mediante simulación y observar la varianza del ACI debida al uso de un número limitado de muestras, se generan una señal ejemplo 2 y una señal ejemplo 4. La ecuación predice un valor para el índice de ACI1 = 0:625, en ambos casos. Se calculan los valores del ACI para N = [256; 1024; 16384] números de muestras distintos. Para cada una de las señales ejemplos y cada N diferente, se generan 1000 realizaciones de datos aleatorios incorrelacionados de varianza unidad y se calcula el ACI de cada realización. Por último, se calcula la media y desviación estándar del ACI para cada señal (Figura 10 (a)). Se observa que a pesar de utilizar un número finito de muestras los valores medios se aproximan al valor teórico del ACI, y cuanto mayor es N, menor es el error que se obtiene en esta aproximación.

Prueba señal ejemplo 5: la siguiente prueba se realiza para mostrar la forma en la que el ACI caracteriza un sistema dinámico conocido. Para esta prueba se generaba 4000 valores equiespciados del valor r entre 0.001 y 4, para 50 valores iniciales de x entre 0.01 y 0.99. Para cada valor de r y de x, se obtiene una serie temporal de 1000 muestras a partir de la ecuación. Sólo utilizan la últimas 500 muestras de cada serie temporal en la estimación del ACI, para eliminar así los posibles efectos transitorios. Para cada r, se calcula y representa la media y la desviación estándar del ACI.

La Figura 10 (b) muestra el comportamiento del ACI respecto al parámetro r del mapa logístico. Para r < 3, el ACI es cero, ya que se trata de un índice que caracteriza el signo de las diferencias y para valores de r < 3 el mapa logístico tiende a un punto fijo, por lo que las diferencias son siempre cero y no existe ningún cambio de signo. Entre 3 < r < 3:68 aproximadamente, el ACI es uno, ya que para estos valores de r el mapa logístico oscila entre dos valores y el signo de las diferencias varía a cada muestra (periodicidad 2), obteniéndose así el máximo valor posible para el ACI. Para r > 3:68, dependiendo del valor de r, podemos obtener un comportamiento caótico o comportamientos periódicos (de periodo mayor que dos) del mapa logístico como el que obtenemos para r = 3:84, donde ACI = 0:5.

Prueba señal ejemplo 6: el fBm y fGn poseen una función de autocorrelación definida por,

En esta prueba, primero se calcula el ACI1 para el caso del fBm, de forma que se obtiene el valor teórico del índice sustituyendo la ecuación 1 en la ecuación 2.
Para estudiar el comportamiento del estimador del ACI, se utilizan 100 valores de equiespaciados entre 1:01 y 2:00. A continuación, se han obtenido 100 realizaciones para el RBF con el simulador wfbm del wavelet toolbox de Matlabr 7.4.0 con _ = [1:1; 1:3; 1:5; 1:7; 1:9], y para 256, 1024 y 16384 muestras. Finalmente, se obtenien las medias y la desviaciones estándar del ACI.
La Figura 10 (c) muestra que a medida que el valor de aumenta el valor del ACI se aproxima más al valor teórico calculado. Por otro lado, el error es mayor cuando se utiliza un número reducido de muestras y disminuye a medida que el número de muestras aumenta.

Prueba señal ejemplo 7: en esta última prueba se obtenien el ACI para dos señales ejemplo 7, una de ellas perteneciente a un paciente sano y la otra a un paciente con insuficiencia cardíaca (Figura 10 (d) y (e)). Se obtiene un ACI mayor para el caso patológico que para el caso sano, lo cual concuerda con los resultados obtenidos en [44]. Estos valores podrían deberse al hecho de que algunas patologías están asociadas a uctuaciones erráticas del ritmo cardíaco con propiedades estadísticas que se
asemejan a las del ruido incorrelacionado, , y en el caso de señales realis, más ruidosas, la señal diferencia podría contener más cambios de signo y por tanto un valor debería ser ACI mayor.

Figura 10: Pruebas con ACI: (a) calculo del ACI en variables aleatorias incorrelacionadas para diferentes longitudes de los datos, se representa la media y desviacion estandar para cada distribucion y el valor teórico esperado viene dado por la linea negra. (b) Calculo del ACI para el mapa logstico, se representa para cada valor de r la media de los valores del ACI para los 50 valores iniciales escogidos. (c) Calculo del ACI para el ruido browniano fraccionario, para cada longitud de los datos se representa la media y la desviacion estandar del ACI, la linea azul representa el valor teorico esperado. (d) y (e) Calculo del ACI en una señal real de un paciente sano y de un paciente con insuciencia cardaca respectivamente.

MÉTODOS TEMPORALES GEOMÉTRICOS

Introducción

El análisis de la HRV mediante métodos estadísticos presenta una limitación importante: dependencia de la calidad de los datos de las series de intervalos NN. La calidad de los datos empeora debido a la presencia de outliers, es decir datos fuera de margen cuya magnitud es muy grande respecto a la media del resto de los datos; debido a artefactos, es decir ruido generado por los instrumentos de medida; e incluso por la capacidad de colaboración del paciente. Todo ello obliga a buscar otros métodos de análisis que se vean menos afectados por la calidad de los datos. Los métodos geométricos surgen con la intención de mejorar la robustez de las medidas frente a los problemas de calidad que presentan las observaciones.

Los métodos geométricos se derivan de las propiedades geométricas de la distribución de los intervalos NN. Generalmente, a partir del histograma de una secuencia de intervalos NN dada se construye una forma geométrica, sobre la que se obtiene una medidas que sirve para caracterizar la HRV. Los métodos para la construcción de la forma geométrica están basados en la mayoría de los casos en:

1. La función densidad de probabilidad de la duración de los intervalos NN (histogramas).

2. La función densidad de probabilidad de la diferencia entre intervalos NN adyacentes (histogramas).

3. El diagrama de Lorenz, también llamado mapa de Poincaré, donde se representa la duración de cada intervalo NN frente a la duración del intervalo NN inmediatamente anterior. En la Figura 11 se presentan dos diagramas de Lorenz con diferente dispersión.

A partir de la forma geométrica resultante de cualquiera de las representaciones anteriores, se puede obtener una medida de la HRV. Existen principalmente tres formas de actuación:

1. Se utiliza alguna medida de la figura geométrica como estimación de la HRV, por ejemplo, la anchura de la base del histograma, para el intervalo de confianza del 95 %. De esta forma se consigue eliminar las colas del histograma, donde se encuentran los outliers.

2. El modelo geometrico se interpola mediante una figura matemática, como puede ser la aproximación de la distribución del histograma por un triángulo, o la aproximación del histograma de la diferencia entre intervalos NN mediante una curva exponencial; la medida de la HRV se deduce entonces de los parámetros de estas figuras geométricas.

Figura 11: Diagramas de Lorenz. (a) Diagrama de Lorenz con muy poca dispersion lo que implica una variabilidad baja. (b) Diagrama de Lorenz con una dispersion mayor, lo que implica que la variabilidad es mayor.

En la Figura 12 se puede apreciar el histograma de una secuencia de intervalos NN aproximado por una gura geometrica, en este caso por un triangulo. A partir de las propiedades del triangulo se pueden establecer distintas medidas para la HRV.

La aplicación de estos métodos debe realizarse preferiblemente sobre registros de 24 horas o mas, con el objetivo de conseguir suficientes datos para construir un histograma fiable. Por otro lado, debido a que los registros de 24 horas contienen periodos de actividad durante el día y periodos de descanso durante la noche, el histograma puede mostrar un aspecto bimodal o incluso multimodal. Por lo tanto, el uso de métodos triangulares dejara de ser adecuado ya que tenderan a sobreestimar la variabilidad del ritmo cardíaco.

Las principales índices geometricos se describen brevemente en las siguientes secciones.

Indice Triangular

Descripción: se trata de la integral de la densidad de distribución (es decir, el número total de los intervalos NN) dividido por el máximo de la densidad de distribución (altura del histograma de los intervalos NN).

Metodología: se emplea una escala discreta en el eje horizontal para la medida de intervalos NN (generalmente de 1/128 segundos), esta medida se puede aproximar por la siguiente expresión:

Propiedades: utilizando esta aproximación, la medida depende del tamaño de las clases, o lo que es lo mismo, de la precisión de la escala discreta de la medida. Con este índice se obtiene una medida global de HRV sobre un registro de 24 horas, que esta mas influenciado por bajas que por altas frecuencias.

Figura 12: (a) Presenta el histograma de una serie NN. (b) Aproximación del histograma por un triángulo. Con este tipo de aproximaciones se eliminan las colas del histograma.

TINN

Descripción: ancho de la línea de base de la distribución de intervalos N, medida como la base de un triángulo que aproxima el histograma de intervalos NN. La figura geométrica se ajusta mediante el método de los mínimos cuadrados.

Metodología: de acuerdo con la Figura 13, se establecen los puntos M y N en el eje temporal, y se construye una función multilineal tal que q(t)=0 para t_N y t_M, q(X)=Y, y tal que la integral

sea mínima para los distintos valores de M y N. Entonces, la medida TINN, expresada en milisegundos, viene dada por la fórmula,
TINN = M – N

Propiedades: igual que con el método anterior, se obtiene una medida global de HRV sobre un registro de 24 horas, que está más inuenciado por bajas que por altas frecuencias.

Índice Diferencial

Descripción: diferencia entre los anchos del histograma de las diferencias entre intervalos NN consecutivos, seleccionadas para dos alturadas diferentes.

Metodología: tomando las diferencias entre intervalos NN adyacentes, se construye un histograma, y se comparan sus distintas anchuras, por ejemplo, en niveles tales como 1000 y 10000 muestras. Se obtiene el resultado en milisegundos.

Figura 13: Representación la distribución de densidad muestral. Para la obtención del parámetro TINN, en primer lugar se debe obtener el máximo de la distribución. A continuación se buscan los valores N y M tal que el triángulo de lados YN, NM, MY, sea el que mejor se ajuste a la distribución de los intervalos NN en el sentido de mínimos cuadrados.

Índice Logarítmico

Descripción: coeficiente de la exponencial negativa que mejor aproxima el histograma de las diferencias entre intervalos NN consecutivos, en sentido de mínimos cuadrados.

Metodología: se construye el histograma del valor absoluto de las diferencias de intervalos NN adyacentes y se ajusta una exponencial negativa, según la ecuación,

El algoritmo realiza un ajuste de una curva exponencial al histograma obtenido, aunque en la literatura también se ha propuesto el ajuste mediante interpolación lineal del histograma en escala semilogarítimica, ver Figura 14.

Dispersión del Diagrama de Lorenz

Definición: método gráfico en el que se representa la duración de cada intervalo NN frente a la duración del intervalo NN inmediatamente anterior.

Metodología: la variabilidad se cuantifica mediante las dispersión de los diagramas de Lorenz. Para caracterizar estos diagramas de forma matemática se puede ajustar una elipse a la nube de puntos de la gráfica como se muestra en la Figura 15. La desviación estándar de los puntos en torno al eje x1 se denomina SD1, este valor mide la anchura de la nube de puntos e indica el nivel de las variaciones de corto plazo de la HRV.
La dispersión de la nube a lo largo de la linea identidad indica el nivel de largo plazo de la HRV y se mide mediante el parámetro SD2, el cual es la desviación estándar alrededor del eje x2. Esta medidas están relacionadas con medidas clásicas de la HRV de
la siguiente manera:

Figura 14: En la figura se muestra el histograma del valor absoluto de las diferencias entre intervalos NN sucesivos, as como la curva exponencial que mejor se ajusta a dicho histograma. En la figura también se muestra el valor del parámetro del exponente que se utiliza como medida de la HRV.

SD1 esta, por lo tanto, relacionado con la funcion de autocovarianza mediante,
De la misma forma, SD2 esta tambien relacionado con la funcion de autocovarianza del siguiente modo,

Sumando las dos ecuaciones anteriores, se obtiene:

y por lo tanto,

Es posible hacer una analisis más geometrico del diagrama de Lorenz, para ello, se definen los siguientes vectores:

donde NNi, representa el i-ésimo intervalo NN y n es el numero de puntos en el mapa. El diagrama de Lorenz de una secuencia de intervalos NN queda as denido por los pares de puntos (xi; yi). El centroide o momento de primer orden (xc; yc) de una distribución de n puntos (x; y) se puede definir como,

que en este caso resulta ser el centroide del diagrama de Lorenz, ya que:

La distancia del i-esimo punto del diagrama a los ejes x1 e x2 (Figura 13) se puede calcular como sigue,

Finalmente, se obtienen las siguientes realciones para los ndices SD1 y SD2, que son equivalentes a las anteriores:

La linea identidad posee una relevancia especial en el diagrama de Lorenz, figura 14. Se define el momento de segundo orden en torno a esta linea identidad como,

donde Di es la distancia de cada punto del mapa a la linea identidad.
Si se tienen en cuenta los puntos que se encuentran, respecto a la linea identidad en el diagrama de Lorenz, por encima (que corresponde a un incremento de los intervalos NN) y por debajo (que corresponden a un acortamiento de los intervalos NN), se puede escribir SD1I como,

donde, DiU es el i-esimo punto sobre la linea identidad y lo mismo para los puntos que se encuentran por debajo, DiD , o sobre, DiO, la linea identidad. De la misma forma, nU representa el numero de puntos sobre la linea identidad, nD es el numero de puntos por debajo de la linea identidad y nO sobre la linea identidad.

Se cumple que:

Figura 13: Ajuste de la elipse al diagrama de Lorenz.

Se denen las contribuciones de los puntos sobre y debajo de la linea identidad a SD12I:

donde,

Dado que la variabilidad de la frecuencia cardaca diere dependiendo de los sujetos, podra ser util algun tipo de normalización sobre los parámetros anteriores. Se define la contribución relativa de las aceleraciones, CUP , y deceleraciones, CDOWN, del ritmo cardíaco al parámetro SD12I como:

Estos parámetros pueden ser comparados entre sujetos para un analisis exploratorio de los datos:

Figura 14: Linea identidad del diagrama de Lorenz.

Pruebas

En los siguientes apartados se muestran diferentes pruebas realizadas para cada uno de los indices presentados en este capítulo. Para la realización de dichas pruebas se utilizan tanto señales sintéticas como señales reales, en concreto se emplean las señales ejemplo 1, ejemplo 2, ejemplo 3 y ejemplo 4 descritas en el capítulo anterior, por lo que no volverán a ser descritas aquí, aunque se introduce aquí una nueva señal sintética que se llama señal ejemplo 8.

Señal ejemplo 8 : esta es una señal sintética creada específicamente para probar el índice logarítmico.
La señal se define de forma que el histograma del valor absoluto de su derivdad posea forma exponencial con un exponente regulable. El procedimiento de creación de la señal es el siguiente:

1. Se crea un histograma con forma exponencial

Donde el valor de la constante K y el dominio de la variable r se ha escogido de acuerdo con los valores habituales que se observan en señales reales de variabilidad. Los valores escogidos son: K = 18e3 y r 2 [0; 500].

2. Se crea una señal, con nombre dy(l), de forma que, el histograma de la se~nal dy(l) sea h(r), por lo tanto, esta señal posee tantos puntos con valor ri como indique el histograma h(ri).

3. Finalmente la señal sintética y(n) se obtiene integrando (sumando) la señal dy(l). Para obtener un resultado más realista, se aleatorizan las posiciones de los valores de la señal dy(l), y se asigna signo positivo o negativo a cada valor de forma aleatoria, antes de la integración,

donde M es el número de muestras de la secuencia dy(l)

Figura 15: Prueba con el índice triangular para una señal sinusoidal mas señal de ruido con valor de la desviación estandar variable.

Pruebas Indice Triangular

Prueba señal ejemplo 3: en esta prueba, se genera una señal sinusoidal con valor medio igual a 600 y amplitud igual a 30, para obtener valores del histograma del orden de los presentes en las señales reales. A la señal sinusoidal se le añade ruido Gaussiano de media cero y desviación estándar variable, que ser_a el factor controlable en la prueba. Se emplean cinco desviaciones estándar distintas, desv = [1; 5; 10; 50; 100; 500], para cada desviación se obtienen diez realizaciones distintas del ruido que se sumam a la señal sinusoidal, y para cada una de estas señales resultantes se computa el
índice triangular.

En la Figura 15 se representan el valor medio y la desviación estándar del índice triangular para cada valor de la desviación estándar del ruido. Como era de esperar, cuanto mayor es la desviación estándar del ruido mayor resulta el histograma de la señal resultante, por lo que es coherente el resultado de un aumento del índice triangular.

Prueba señal ejemplo 7: en esta prueba, se analia el comportamiento cuantitativo delíndice triangular en un caso de señal real de HRV. El resultado se puede observar en la Figura 16. Cabe destacar la bimodalidad presente en el histograma de las señales reales de HRV en registros Holter de 24 horas, probablemente debida a los ciclos circadianos.

Pruebas TINN

Prueba señal ejemplo 3: en esta prueba, se utiliza la señal ejemplo 3 con los mismos valores ya expuestos en la prueba señal de Pruebas Indice Triangular. Por lo tanto, se genera una señal sinusoidal y diez realizaciones de ruido Gaussiano para cinco desviaciones estándar diferentes. Se calcula el valor del TINN para la señal resultante de la suma de la sinusoidal y el ruido Gaussiano, de forma que se obtienen diez valores diferentes del TINN para cada desviación estándar de la señal aleatoria. En la Figura 17 se representan el valor medio y la desviación estándar del TINN para cada valor de desviaci_on est_andar del ruido. Como se puede observar el valor medio del TINN es mayor cuanto mayor es la desviaci_on del ruido, aunque el error de la estimación también es mayor debido a la mayor potencia del ruido.

Figura 16: Prueba con el ndice triangular para una señal real de HRV.

Prueba señal ejemplo 7: en esta prueba, se muestra gráficamente el cálculo del TINN en una señal real de HRV. La Figura 18 presenta la interpolación triangular del histograma, en trazo rojo grueso superpuesto al histograma, a partir de la cual se obtiene el valor resultante del TINN.

Pruebas Indice Logarítmico

Prueba señal ejemplo 8: en esta prueba, se obtienen cinco realizaciones distintas de la señal ejemplo 8. El parámetro se puede modificar de forma que se pueden obtener señales sintéticas con distintas caídas exponenciales del histograma, por lo tanto, con distinto índice logarítmico. Para esta prueba, se generan cinco señales con cinco valores de diferentes, = [_0:015;_0:02;_0:03;_0:04;_0:05], para poder estudiar la precisión del método de estimación del índice logarítmico. En la Figura
se pueden observar los resultados del cómputo del índice logarítmico para cada una de las cinco señales sintéticas.
En cada gráfica se representa el histograma del valor absoluto de las diferencias, y superpuesta la curva exponencial que mejor ajusta (en sentido de mínimos cuadrados), en trazo rojo grueso.
Prueba señal ejemplo 7: en esta prueba, se muestra un ejemplo del cálculo del índice logarítmico, en una señal real de HRV.

Pruebas Diagrama de Lorenz

Prueba señal ejemplo 6: en la primera prueba, para los índices basados en el diagram de Lorenz, se computan los índices SD1; SD2; cup y cdown, a partir del diagrama de Lorenz de dos realizaciones del fBm, es decir, señal ejemplo 6, con dos valores diferentes del exponente de Hurst, a saber, H = 0:3 y H = 0:7.

Para la realización con el exponente de menor valor, se obtiene una señal cuyo comportamiento es muy irregular, por lo que es esperable una variabilidad a corto plazo mayor que para la señal con exponente mayor. Por lo tanto, los valores de los índice cup y cdown deben ser parecidos debido a la simetría entre las aceleraciones y deceleraciones (correlación negativa entre incrementos).
En cuanto a la señal de exponente H = 0:7; es de esperar una variabilidad de largo plazo mayor que para la señal con exponente menor, debido a las tendencias presentes en las señales con H grande. En cuanto a los índices de simetría, deberían resultar valores muy distintos, es decir, que valores que revelen una asimetría en cuanto a aceleraciones y deceleraciones, debido a la persistencia que presentan los fBm con H > 0:5, es decir, debido a la correlación positiva entre incrementos (ver Fig. 19)

Prueba señal real: en esta prueba, se utilizan dos señales reales, una correspondiente a un paciente sano y la otra correspondiente a un paciente con insuficiencia cardíaca. Los resultados se presentan en la Figura 21. Es posible apreciar visualmente las diferencias en los diagramas de Lorenz pertenecientes al caso sano y al caso patológico. Así mismo, los índices muestran la disminución de variabilidad en el caso del paciente con insuficiencia.

Figura 17: Prueba con el TINN, para una seÑal sinusoidal y distintas realizaciones de ruido Gaussiano.

Figura 18: Prueba con el TINN para una señal real de registro Holter de 24 hh.

Figura 18: Prueba con el índice logartmico para cinco señales sintéticas con diferente

Figura 19: Prueba con el índice logarítmico para una señal real.

Figura 20. Prueba con señal fBm para el diagrama de Lorenz.

Figura 21: Pruebas con el diagrama de Lorenz para dos señales reales. (a) Diagrama de Lorenz para una señal de un paciente sano. (b) Diagrama de Lorenz correspondiente a un paciente con insuciencia cardíaca. (c) Valores de los parametros obtenidos a traves del diagrama de Lorenz para la señal (a). (d) Valores de los parametros obtenidos a través del diagrama de Lorenz para la señal (b).

MÉTODOS ESPECTRALES

El estudio de la HRV mediante métodos espectrales se basa en el hecho de que las inuencias de las dos ramas del SNA sobre el corazón, así como algunos otros sistemas como puede ser el respiratorio, tiene comportamientos oscilatorios bien definidos y distintos. De esta forma se puede establecer la siguiente división en bandas espectrales según qué sistema es el que tiene una mayor contribución en dicha banda.
DC-0.03 Hz, se conoce como banda de muy baja frecuencia (Very Low Frequency, VLF), en señales de 24 horas se hace una división en ultra baja frecuencia (Ultra Low Frequency, ULF) y VLF. Se asocian a mecanismos de regulación de larga duración (long-term), probablemente relacionados con la termorregulación, el sistema renina-angiotensina o a otros factores humorales, aunque no se ha podido establecer una relación directa y bien definida.

0.03-0.15 Hz, conocida como banda de baja frecuencia (Low Frequency, LF), generalmente centrada alrededor de la frecuencia 0.1 Hz. La interpretación fisiológica aceptada es la de contribuciones tanto de la rama simpática como de la rama parasimpática, aunque se ha probado que un aumento en la potencia de dicha banda se atribuye a aumento de la actividad simpática.

0.15-0.4 Hz, conocida como banda de alta frecuencia (High Frequency, HF), asociada al ritmo respiratorio, y por tanto relacionada a cambios de presión y movimientos mecánicos en el pecho debidos a la respiración, y se debe a la actuación del nervio vago sobre el corazón, por lo que se esta banda está asociada a una contribución de la rama parasimpática.

En la Figura 22 se puede observar el espectro de una representación de la señal de HRV para una persona sana, donde se identifican las tres bandas de frecuencia. Para obtener valores cuantitativos cuando se efectúa análisis de la HRV mediante métodos espectrales, se han establecido una serie de medidas que se deben obtener a partir de la densidad espectral de potencia, y que son:

1. Potencia Total: Se define como la integral bajo la curva de la densidad espectral de potencia (DEP) sobre la gama de frecuencias que va desde DC hasta 0.4 Hz. Se corresponde con la varianza de todos los intervalos RR. Las unidades son ms2.

2. Potencia en la banda ULF: Integral de la DEP en la banda ULF, definida en el rango de frecuencias [DC; 0:003] Hz. Las unidades son ms2.

3. Potencia en la banda VLF: Integral de la DEP en la banda VLF, es decir, en el rango de frecuencias [0:003; 0:04] Hz. Las unidades son ms2.

4. Potencia en la banda LF: Integral de la DEP en la banda LF, definida en la gama de frecuencias [0:04; 0:15] Hz. Las unidades son ms2.

5. Potencia en la banda HF: Integral de la DEP en la banda HF, por tanto en el rango de frecuencias [0:15; 0:4] Hz. Las unidades son ms2.

6. LF=HF: Razón de las potencias pertenecientes a las bandas LF y HF respectivamente.

Figura 22: Densidad espectral de potencia para una persona sana. donde se aprecian las tres bandas espectrales a las que se asocia la influencia de un sistema distinto sobre la actividad cardíaca.

Un hecho importante a tener en cuenta es que la mayoría de métodos espectrales tienen como condición de aplicación que la señal sea uniformemente muestreada, es decir, que los intervalos entre sucesivas muestras sean constantes, pero esto no ocurre en las señales de HRV, ya que estas, por su naturaleza poseen un espaciado entre muestras irregular. Por lo tanto, cuando se aplican las técnicas espectrales se debe realizar una elección entre los dos enfoques siguientes:

• Ignorar los requerimientos, y por lo tanto, suponer que las observaciones ocurren de hecho a intervalos uniformes. Cuando se utiliza este enfoque hay que tener en cuenta que realmente no nos estamos moviendo al dominio de las frecuencias, si no a un dominio de análisis del número de latido, que en la literatura aparece en ocasiones por el nombre de dominio beatquency.

• El otro enfoque es remuestrear la serie temporal escogida para representar la señal de variabilidad, de forma, que se obtenga una serie temporal con muestras a intervalos uniformes, uno de los métodos más empleados es el remuestreo mediante modelos de Integrated Pulse-Frequency Modulation (IPFM).

A continuación se presentan los métodos de estimación de la densidad espectral más utilizados en la literatura para el análisis de señales de HRV [12, 14, 16]. Dentro de los métodos de estimación podemos establecer una primera división entre métodos paramétricos y métodos no paramétricos. En general, los métodos no paramétricos poseen las ventajas de ser simples (basados en la FFT) [62] y poseer una alta
velocidad de procesamiento, mientras que las principales ventajas de los métodos paramétricos son que permiten una mayor precisión en la estimación de la PSD incluso con un número reducido de muestras y presentan las componentes espectrales suavizadas, la cuales pueden ser fácilmente distinguibles, lo cual facilita la etapa de postprocesado, por contra lo métodos paramétricos son más complejos que los no
paramétricos ya que se debe verificar la adecuación del modelo escogido.

Periodograma

Definición: el periodograma es una estimación de la densidad espectral de potencia, asumiendo que la señal es una realización de un proceso estocástico estacionario.

Metodología: dada una serie temporal x[n], de longitud, N, se puede estimar la función de autocorrelación como,

Aplicando la transformada de Fourier a ambos miembros, y realizando algunas manipulaciones algebraicas, se obtiene el periodograma como,

Es decir, podemos obtener el periodograma a partir de los datos, tomando el modulo al cuadrado de la transformada de Fourier de los datos,

Hay que tener en cuenta que estamos utilizando una secuencia de datos finita, y por lo tanto es equivalente a multiplicar la serie temporal por una ventana cuadrada, este enventanado produce unos efectos indeseados sobre la estimación de la densidad espectral, un modo de mejorar dicha estimación es la utilización de ventanas no rectangulares, de forma que el corte que se produce en los bordes de la secuencia no sea tan abrupto. A los periodogramas obtenidos utilizando ventanas no rectangulares se les conoce con el nombre de Periodogramas modificados.

Propiedades: las principales propiedades del periodograma son:

1. Es un estimador sesgado, es decir,

donde Sxx(f) es la verdadera densidad espectral de potencia de la serie x[n], mientras que W(f) es la transformada de Fourier de la ventana utilizada, por lo tanto el sesgo se debe al efecto de la ventana.

2. Es un estimador inconsistente, y se puede mostrar que,

Otro de los parámetros a tener en cuenta a la hora de realizar la estimación espectral a partir del periodograma es el número de puntos que se utilizarán en la FFT.

La Figura muestra el periodrograma de una señal de HRV de una persona sana.

Periodograma de Welch

Definición: el periodograma de Welch se obtiene como la media de periodogramas calculados sobre segmentos de la secuencia de datos. Este método fue inicialmente propuesto por Barlett y posteriormente modificado por Welch, que introdujo el promediado de periodogramas modificados y el solape de segmentos sobre los que se computa cada uno de los periodogramas.

Metodología: para obtener el periodograma de Welch, primeramente se divide la secuencia de datos en K segmentos de L muestras de longitud y se aplica a cada segmento una ventana [50] de la misma longitud que la de los segmentos, es decir, L, lo que resulta en una serie:

donde R representa el parámetro que controla el solape, este se produce cuando R < L, mientras que no existe solape si R = L. El periodograma para el segmento r-esimo es:

donde U es un factor de normalizacion para la ventana. El periodograma de Welch resulta de tomar la media de los periodogramas sobre las K estimaciones de los distintos segmentos:

Propiedades: con esta tecnica de estimacion se consigue reducir el sesgo de la estimación, mientras que respecto a la varianza, obtenemos un estimador asintoticamente consistente,

Es posible demostrar que con un solapamiento del 50% se pude conseguir reducir la varianza en un factor de 2, y que un solapamiento mayor no contribuye a reducir más la varianza.
Como parámetros libres se pueden destacar los mismos que para el periodograma, además de los impuestos por el promediado, como son: el solapamiento y el número de datos en cada segmento.

En la Figura se puede observar la estimación de la densidad espectral de potencia, por el método del periodograma de welch, para le mismo paciente que se muestra en la Figura 23.

Figura 23. Periodograma de Welch de una señal de HRV correspondiente una persona sana.

Método Autoregresivo

Definición: es un método paramétrico, también conocido como método de máxima entropía (Maximum Entropy Method, MEM), o método \todo{polos” (\all-poles” method). El método autoregresivo (AR) de estimación espectral, se basa en crear un modelo para la generación de la señal con un conjunto de parámetros que se pueden estimar a partir de los datos, es decir, la serie temporal, de forma que se puede estimar la densidad espectral de potencia a partir del modelo creado y a partir de los parámetros obtenidos.

Metodología: en el método AR, la señal x[n], que es la señal de HVR, se interpreta como la salida de un sistema lineal, cuya entrada xe[n] es ruido blanco, ver Figura 24.

La secuencia de salida se escribe como una ecuación en diferencias:

La funcion de transferencia resultante es,

La estimación de la densidad espectral de potencia se puede computar a partir de los parámetros obtenidos mediante,

donde p se conoce como el orden del modelo. Las cuestiones que quedan por resolver son el cálculo de los parámetros del modelo, y la elección del orden del modelo. Para la determinación de los parámetros existen distintas técnicas que se describen brevemente a continuación:

• Método de la autocorrelación o Yule-Walker: en este método se busca minimizar el error de predicción, que se define como, ep[n] = x[n]  – ^x[n]

Para este propósito utiliza la matriz de autocorrelación de los datos, que se relacionan con los parámetros del modelo mediante las ecuaciones normales, conocidas por el nombre de Yule- Walker. Una pequeña variación de este método es el conocido por el método de la covarianza, en el que se utiliza la matriz de datos en lugar de la matriz de autocorrelación.

• Método de la covarianza modificada: se trata de un método similar al de la covarianza, pero en este caso se trata de minimizar la media de los errores cuadráticos de predicción lineal hacia delante y hacia atrás,

donde el error de predicción lineal hacia delante ef p [n] se dene como el error de la ecuación, mientras que el error de predicción lineal hacia atras se define como,

Los superndices f y b vienen del ingles forward y backward respectivamente.

• Metodo de Burg: se basa en la minimizacion conjunta de los errores de prediccion hacia delante y hacia atrás, al igual que en el metodo de covarianza modicada, pero en el método de Burg se hace uso explcito de la recursión de Levinson-Durbin para el proceso de minimización.

Figura 24: Estimacion de la DEP mediante un modelo AR de covarianza modicada de orden 12.

Propiedades: La principal ventaja que supone la utilizacion de metodos AR es que el espectro estimado es visualmente mucho mas limpio que el caso de los metodos no paramétricos, lo que permite un postprocesado del espectro mas sencillo. El principal problema que se plantea es la elección del orden del modelo, y por tanto las asunciones en cuanto a complejidad que se llevan a cabo, ya que el espectro puede no sólo reflejar la estructura espectral de la se~nal de entrada, sino ademas
reflejar las asunciones establecidas. En la Figura 24 se puede observar la estimación de la densidad espectral de potencia mediante el metodo AR de covarianza modicada.
Como principal parámetro libre se puede establecer la elección del orden del modelo, para su elección existe una serie de tecnicas estadísticas, como por ejemplo el criterio de información de Akaike (Akaike’s Information Criterion, AIC).

Método autorregresivo, Descomposición en Sub-bandas.

Definición: método cuyo objetivo es la caracterización del espectro de potencia de la señal en función de sus componentes espectrales intrínsecas, lo que se consigue por medio de la expansión en fracciones parciales del espectro de la función de transferencia del modelo AR, H(z).

Metodología: en el modelo AR, el espectro de potencia de la señal de salida (el objetivo del problema en este caso), se puede escribir como,

donde A(z) = 1H(z) . Dado que A(z) es una funcion racional de z, el espectro de potencia de x se puede expresar en funcion de los polos de A(z),

El objetivo es caracterizar Sx(z) en funcion de sus componentes espectrales intrínsecas, para ello se utiliza el metodo de expansión de fracciones parciales, mediante el cual, la función de transferencia H(z) se puede expresar como una suma de términos de primer orden,

Donde los coecientes de cada termino simple se puede obtener como,

Considerando que el modelo AR esta denido por pares de polos complejos conjugados, de forma se que el espectro de potencia Sx(z) sea una funcion par, la expansión en fracciones parciales se puede expresar de forma alternativa como una suma de funciones de transferencia de segundo orden,

Cada función de transferencia individual Hi(z) dene un modelo ARMA de segundo orden,

donde A2i-1 = A2i ya que los pares de polos conjugados como c2i-1 = c2i, provocan en la expansión que los correspondientes coeficientes sean a su vez pares de coeficientes complejos conjugados.

Mediante algunas operaciones algebraicas se llega a

Para el resto del desarrollo es necesario establecer la hipotesis de que las funciones de transferencia Hi(z) no se solapan, es decir,

Esta asunción no es muy restrictiva ya que el espectro de la señal de HRV posee componentes espectrales bien separada unas de otras, de hecho esa es la motivación este método. Empleando la hipotesis, el espectro de potencia se puede expresar como,

En la Figura 25 se puede observar un ejemplo de descomposición en dos componentes.
La densidad espectral de potencia se puede expresar como,

Sx(ej!) = Sx1(ej!) + Sx2(ej!) = H1(ej!)2e + H2(ej!)2

Figura 25: Descomposición del espectro de potencia de un modelo AR mediante la expansión de fracciones parciales.

Propiedades: la potencia para cada una de las componentes espectrales obtenidas mediante la expansion de fracciones simples se puede computar mediante la transformada integral inversa que relación donde el contorno de integracion C debe estar dentro de la region de convergencia, en este caso dentro del crculo unidad. Por lo tanto, la potencia total se obtiene.
Este metodo presenta los mismos parametros libres que el metodo AR convencional.

Periodograma de Lomb

Definición: el Periodograma de Lomb es un método de estimación espectral aplicable a secuencias irregularmente muestreadas, por lo que parece un método adecuado para aplicar a la señal de HRV.

La idea que subyace al método es que es posible ajustar una onda sinusoidal de cualquier frecuencia a un conjunto de observaciones utilizando mínimos cuadrados, una vez realizado el ajuste se obtiene un punto del espectro, posteriormente se substrae dicha componente y se ajusta otra onda sinusoidal, el procedimiento se repite hasta que se obtiene el espectro con la resolución frecuencial deseada. En la Figura 6.8 se muestra la estimación del periodograma de Lomb sobre la señal que se ha utilizado en los ejemplos anteriores.

Propiedades: la principal ventaja que presenta el metodo de estimacion mediante el periodograma de Lomb es que se puede llevar a cabo aunque los datos no este uniformemente muestreados o incluso si existen \huecos” en los datos, sin necesidad de realizar estimaciones de las muestras que faltan o interpolaciones. Se pueden presentar problemas si existen demasiados \huecos” en los datos.

Figura 26: Descomposición en sub-bandas mediante el método de la expansión en fraccinoes parciales de un modelo AR

Figura 27: Estimación de la densidad espectral de potencia medianta el periodograma de Lomb.

MÉTODOS DE TIEMPO-FRECUENCIA

La mayor limitación de los métodos basados en análisis espectral es que no proporcionan información temporal sobre la ocurrencia de las diferentes frecuencias de la señal. La transformada de Fourier únicamente refleja las frecuencias que existen durante todo el periodo de observación de la señal. Este tipo de análisis es adecuado para señales estacionarias, pero no es adecuado para señales no estacionarias con contenido espectral dependiente del tiempo. Por esta razón surgen los métodos que analizan las señales considerando tanto el tiempo como la frecuencia, en los que se representa la frecuencia de la señal en cada instante temporal [133, 6, 122, 27, 43, 23, 115].

Short-Time Fourier Transform

Definición: el método Short-Time Fourier Transform (STFT) consiste en dividir la señal en segmentos multiplicándola por una ventana deslizante en el tiempo y calcular después la transformada de Fourier de cada segmento.

Metodología: en esta aproximación la definición de la transformada de Fourier es modificada de forma que una ventana deslizante w(t) define cada segmento de tiempo a analizar, el resultado es una función X(t; ), que viene dada por,

Análogamente al cálculo del periodograma, el espectrograma de la señal x(t) se obtiene de la siguiente forma,

El espectrograma es una distribución de valores reales y no negativos que proporciona una representación tiempo-frecuencia de la señal.

Propiedades: la longitud de la ventana determina la resolución espacial y temporal; una ventana corta proporcionará una buena resolución temporal pero una resolución en frecuencia pobre y con una ventana larga ocurrirá lo contrario. Se trata de un método computacionalmente eficiente y de un estimador robusto, pero es preciso asumir estacionariedad dentro de la ventana temporal y establecer un compromiso entre las resoluciones frecuencial y espacial.

Distribucion de Wigner-Ville

Definición: la distribución de Wigner-Ville (Wigner-Ville Distribution, WVD) es una distribución cuadratica de la cual se derivan caractersticas temporales y frecuenciales de la señal, pertenece a la clase de Cohen con Kernel igual a uno.

Metodología: se define La WVD en terminos de la señal, x(t), o de su espectro x(!) como sigue,

Propiedades: la propiedad más importante de esta distribución es su excelente resolución combinada en frecuencia y tiempo, pero posee tambien otras propiedades interesantes, Los valores de la WVD son reales, ya que,

W(t; !) = W(t; !)
Si la señal x(t)=0 para jtj > t0, entonces W(t; !) = 0 para jtj > t0, lo que quiere decir que si la señal es cero fuera un cierto intervalo, la distribución tambien lo será.
De la misma forma, en el dominio frecuencial, si la señal X(!) = 0 para j!j > !0, entonces W(t; !) = 0 para j!j > !0. Esta propiedad en tiempo y frecuencia no se cumple en el espectrograma.
Si integramos la WVD respecto al tiempo, obtenemos el espectro total,

Esta propiedad se conoce como condición marginal en frecuencia.
La operación correspondiente respecto a la frecuencia, proporciona la energía instantánea de la señal,

Esto se conoce a su vez como condicion marginal en el tiempo.
Sin embargo, la WVD posee tambien algunos inconvenientes que se deben tener en cuenta:

Al contrario de lo que sucede en el espectrograma, la WVD no contiene unicamente valores positivos, incluso cuando lo que refleja es la energía de la señal en el dominio tiempo-frecuencia.
En la practica este problema no resulta especialmente engorroso ya que las regiones con valores positivos de la WVD resultan corresponder bien con la estructura tiempo-frecuencia esperada.

AR Modo Bloque

Definición: se trata de un modelo paramétrico para secuencias de datos pequeñas pertenecientes a una serie temporal mayor. Siguiendo el mismo principio de la STFT, PSD en cualquier instante de tiempo n, puede ser también obtenida aplicando algoritmos AR clásicos a secuencias de datos cortas centradas en n. El método forward-backward least square (FBLS) es uno de los mas efectivos, este se basa en la minimización del error cuadrático medio de predicción linear hacia adelante y hacia atrás.

Metodología: Si designamos como n a la muestra central de la secuencia x(1); :::x(N) y suponemos que el periodo de muestreo es T = 1, entonces, el espectro asociado al instante n se puede obtener de la siguiente ecuación,

 

 

Propiedades: la resolución temporal viene dada por la duración de la secuencia de datos, sin embargo, al contrario de lo que ocurre con la STFT, la duración de la secuencia, en principio no afecta a la resolución espectral. Si la señal es no estacionaria, la elección del orden del modelo puede ser crítica.

AR Modo Secuencial

Definición: una alternativa a los métodos en modo bloque es, la estimación de los parámetros del modelo, actualizándolos para cada muestra disponible, lo que permite una monitorización en instantánea. Uno de los métodos más conocidos es el Recursive Least Square (RLS), este algoritmo minimiza el error cuadrático medio ponderado exponencialmente de todos los datos disponibles hasta el instante n.

Metodología: el espectro asociado al instante temporal n, puede ser obtenido de la siguiente manera,

 

 

Propiedades: estos métodos son válidos para la estimación de PSD en señales con características que no varíen considerablemente en el tiempo, la señal debe ser casi estacionaria dentro de la ventana exponencial. Poseen una buena resolución frecuencial y mejoran la resolución temporal. Sin embargo, esta resolución temporal depende del tamaño de la ventana y es preciso buscar un compromiso entre dicha resolución y la estabilidad del sistema. El punto crítico es la elección del factor que controla el tamaño de la ventana que pesa los errores y por tanto, la cantidad de memoria del sistema, un valor de grande, podrá dar problemas en la convergencia de los parámetros.

Transformada Wavelet

Definición: las wavelets son funciones matemáticas que descomponen las señales en sus diferentes componentes frecuenciales, para después estudiar cada componente con la resolución adecuada a su escala. Se trata por tanto, de una descomposición de la señal en un conjunto de funciones base obtenidas mediante desplazamientos y escalador de una función madre [80, 122]
Metodología: la transformada Wavelet continua (Continuous Wavelet Transform CWT) de una señal continua x(t) se define como la correlación entre x(t) y una versión escalada y desplazada de (t), donde (t) es la función wavelet madre,

 

 

 

 

 

 

 

 

 

 

Propiedades: la transformada wavelet es adecuada para el análisis de señales de naturaleza no estacionaria, como las señales biomédicas, ya que no existen requisitos en ese aspecto. El análisis mediante wavelets consigue una buena resolución tanto en tiempo como en frecuencia, con dos nuevos grados de libertad de las funciones base, el escalado y la traslación, que permiten el análisis de la presencia conjunta de las formas de onda globales (de escala grande) así como de las estructuras finas (de escala pequeña) en las señales.

MÉTODOS NO LINEALES BASADOS EN LA TEORIA DEL CAOS

Un posible enfoque para el análisis de la HRV es pensar que el origen de las complejas fluctuaciones del ritmo cardíaco se deben en parte al caos determinista, y que las distintas patologías provocan una disminución de dicha variabilidad no lineal.

Partiendo de esta idea es natural intentar caracterizar la dinámica de la HRV mediante los métodos habituales utilizados en el análisis de sistemas caóticos. Los métodos más usuales y que además han sido propuestos en la literatura para aplicarse al estudio de la HRV son, el Exponente de Lyapunov y la Dimensión de Correlación.

Ambos índices están basados en medidas sobre el atractor del espacio de fases reconstruido a partir de la señal de HRV.

Dimensión de Correlación

Definición: la dimensión de correlación es una estimación de la dimensión fractal del atractor que surge en el espacio de fases reconstruido, dado que dicho atractor, si está originado por un sistema dinámico caótico, presenta estructura fractal. La dimensión fractal está relacionada con el número mínimo de variables del sistema necesarias para modelar el atractor, y por tanto relacionado con la complejidad del proprio sistema.

Metodología: la dimensión de correlación puede ser escrita como,

 

 

 

 

Propiedades: la principal característica es que permite una estimación de la dimensión fractal del atractor del sistema a partir de una serie temporal o señal originada por dicho sistema, por tanto es de esperar que el resultado mejore el conocimiento que se tiene sobre el sistema subyacente. 

Cabe decir que es una medida que hay que aplicar con un cierto cuidado dado que es una medida que resulta muy perturbada cuando se aplica a señales no estacionarias o datos con componentes estocásticas, como por ejemplo procesos con espectro 1=f. Teóricamente este tipo de procesos deberán tener una dimensión de correlación infinita, pero la aplicación los algoritmos de estimación pueden ofrecer resultados finitos.

En conclusión, la dimensión de correlación es una herramienta válida para cuantificar la auto-similaridad del atractor cuando se sabe de antemano que está presente, pero debe ser utilizada con mucho más cuidado cuando se quiere establecer dicha característica fractal, es decir, cuando no se se tiene la certeza de que los datos son de baja dimensionalidad y deterministicos.

Desde el punto de vista de la utilización para el análisis de la señal de HRV, se asocia una disminución de la dimensión de correlación con distintos estados patológicos, aunque la aplicabilidad tanto teórica como práctica (clínica) no ha sido aún establecida, debido principalmente a que el número de datos necesario es muy grande y es difícil cumplir las premisas de estacionariedad y determinismo.

Máximo Exponente de Lyapunov

Definición: el máximo exponente de Lyapunov es una medida que evalúa cómo de caótico es un sistema, en realidad, mide la velocidad a la que se separan dos trayectorias que se encuentran cercanas en el espacio de fases en un momento dado. Una característica de los sistemas caóticos es que esta divergencia de trayectorias se produce a una velocidad exponencial, provocada por la propiedad de la dependencia sensible a las condiciones iniciales, y resulta en valores de los exponentes de Lyapunov positivos, por lo que una características de los sistemas caóticos es valores del máximo exponente de Lyapunov positivos, o dicho de otra forma, un valor positivo del máximo exponente de Lyapunov es un fuerte indicador de presencia caos.

Metodología: el algoritmo propuesto para el cálculo del máximo exponente de Lyapunov a partir de la serie temporal, analiza como es la separación exponencial de trayectorias cercanas en el espacio de fases. El algoritmo se apoya, como es habitual, en el espacio de fases reconstruido a partir de la serie temporal.

Por lo tanto se debe elegir en primer lugar la dimensión embebida y a continuación el retardo temporal adecuados para poder realizar la reconstrucción del espacio de fases. Seguidamente se toma un punto del espacio de fases yn0 , y a partir de este punto se elige un vecindario de radio y se toman los puntos que se encuentran dentro del vecindario, es decir, aquellos puntos que se encuentran a una distancia menor que de yn0 . Una vez hecho esto, se calcula un promedio de las distancias de todas las coordenadas para distintos instantes de tiempo n y, se calcula su logaritmo neperiano. Finalmente se promedia el valor resultante sobre un número N de puntos,

 

 

Los puntos de referencia yn0 son los vectores embebidos, U _yn0_ es el vecindario con radio del punto yn0 . El tamaño del vecindario debe ser tan pequeño como sea posible, pero suficientemente grande para que, en promedio, cada punto de referencia tenga un número suficiente de vecinos.

Propiedades: una de las principales ventajas del cómputo del máximo exponente de Lyapunov es la sencillez de su estimación a partir de los datos, además es invariante frente a transformaciones de los datos, bajo la hipótesis de que estas transformaciones sean suaves. De manera análoga a como ocurrá con la dimensión de correlación, se asocia una disminución de los valores del máximo exponente de Lyapunov para casos patológicos. La aplicabilidad de este método adolece de los mismos problemas que presentaba la dimensión de correlación, a saber: la gran demanda de datos y la necesidad de por un lado de determinismo y por otro de estacionariedad de los datos.

MÉTODOS FRACTALES

Otro enfoque posible para al análisis de la señal de HRV es el análisis de las propiedades fractales de la propia señal. Estos métodos se basan en que hay ciertas características de la señal de HRV que sugieren la presencia de una estructura fractal, a saber: correlaciones de larga duración, estructura espectral 1=f y propiedades de auto-semejanza. Es esperable que en estados patológicos exista una ruptura (perdida) de la estructura fractal como reflejo de la pérdida de complejidad del sistema.

Exponente del Espectro 1/f

Definición: se define como el exponente presente en el espectro de los procesos con características de auto-semejanza (fractales), se conocen con el nombre de ley de potencia espectral, o espectro 1=f.

El exponente B caracteriza la estructura fractal de la señal, y por tanto se puede utilizar como una medida fractal de dicha señal. El rango de valores para el parámetro B es: [_1; 3].

Metodología: la estimación del exponente B se realiza mediante la pendiente de la recta de regresión que mejor ajusta, en el sentido de mínimos cuadros, a la densidad espectral en escala logarítmica. Esta estimación se basa en que si tomamos logaritmos a un espectro del tipo 1=f encontramos que, 

log S(f) = B log f + K

de donde es fácil observar que la pendiente de una recta ajustada a la densidad espectral permite realizar una estimación del parámetro B.

Propiedades: desde el punto de vista práctico, en cuanto a su aplicación en el análisis de la señal de HRV, cabe destacar que la estimación del parámetro B se lleva a cabo utilizando apenas el rango de las VLF y ULF, dado que en este rango de frecuencias es más acusada la tendencia espectral 1=f, como se puede apreciar en la Figura.

Se ha encontrado, en distintos estudios, que dicho parámetro sufre una determinada alteración en pacientes con algún tipo de patología cardiovascular. En pacientes normales el parámetro B resulta ser igual a 1, mientras que los pacientes con infarto de miocardio poseen pendientes más acusadas y por lo tanto valores de B mayores.

Parámetros libres del algoritmo: dado que la estimación del parámetro B se realiza a través de una estimación espectral de la señal de HRV, es claro, que las cuestiones algorítmicas más importantes a tener en cuenta son las relacionadas con la estimación, como son: el método de estimación espectral, los parámetros del método (ventana, muestras, . . . )

Se puede establecer como parámetro libre la forma en la que se realiza el ajustes de la recta de regresión al espectro, aunque habitualmente se utiliza la técnica de mínimos cuadrados, pero podrá considerarse otro tipo de ajuste, por ejemplo un ajuste pesado, en el que se tuviese en cuenta la distribución de las muestras (puntos) en la región del espectro que se emplea en la regresión.

Figura: Espectro 1=f en una señal real de HRV.

Exponente de Hurst

Definición: el exponente de Hurst, H es una medida cuantitativa de las propiedades de auto-semejanza (fractales) de la señal. La auto-semejanza para un proceso XH(t), con exponente de Hurst H se puede describir como
XH(t) d= c_HXH(ct)

donde d= debe interpretarse en el sentido de igualdad de distribuciones de probabilidad. La propriedad de auto-semejanza implica que si observamos un proceso XH(t) en una determinada escala de tiempo vamos a observar las mismas propriedades estadísticas que el proceso original, pero normalizadas por un factor de escala.

Metodología: en los procesos auto-semejantes las propiedades fractales tienen un reflejo en el dominio de las frecuencias, en concreto, poseen un espectro 1=f, de forma que, el exponente de Hurst tiene una relación analítica con el exponente del espectro 1=f

Una forma de estimar el exponente de Hurst es a través de la estimación espectral. Se han propuestos métodos de estimación espectral más adecuados para la estimación del parámetro H, aprovechando las propiedades de la transformada Wavelet, que se adecuan mejor a las características de los procesos auto-semejantes.

Los distintos métodos de estimación propuestos son:

• Estimación espectral clásica.
• Reconstrucción del espectro mediante la Transformada de Wavelet Contínua (CWT).
• Reconstrucción del espectro mediante la Transformada de Wavlet Discreta (DWT).
• Método basado en la varianza de los detalles de la DWT.

De forma que, tomando logaritmos se puede ajustar una recta de regresión cuya pendiente se utiliza como estimación del parámetro H donde R es la diferencia entre al valor máximo y el valor mínimo que toma la serie en un período de tiempo T y S es la desviación típica de la serie en ese mismo período. Este método posee una serie de dificultades prácticas que desaconsejan su utilización, como por ejemplo, un importante sesgo o falta de robustez con señales reales.

Propiedades: el exponente de Hurst posee unas propiedades similares a las expuestas para el parámetro del espectro 1=f. La gama de valores posibles para H es [0; 1]. Debido a la relación que existe entre el exponente H y la correlación entre distintos incrementos del proceso, es posible categorizar los procesos de la siguiente forma:

1. 0 < H < 1=2, los sucesivos incrementos poseen una correlación negativa, lo que implica que las gráficas oscilan rápidamente, por tanto muy son rugosas

2. H = 1=2, los incrementos están incorrelacionados, es el caso del movimiento Browniano ordinario.

3. 1=2 < H < 1, los incrementos poseen una correlación positiva, es decir, exista una persistencia, lo que resulta en una gráfica suave.

Parámetros libres del algoritmo: los métodos de estimación basados en la transformada wavelet, poseen además de las peculiaridades de la estimación espectral, algún grado de libertad, como puede se la elección de la wavelet, dado que esta se debe ajustar lo mejor posible a cada señal en particular, o el número de niveles de descomposición.

Figura: Tendecia de cada segmento de la señal integrada.

Análisis de Fluctuaciones Sin Tendencia

Definición: el análisis de fluctuaciones sin tendencia (Dentrended Fluctuations Analysis, DFA) permite cuantificar las propiedades de correlación a largo plazo asociadas a procesos auto-semejantes, y en concreto en series temporales fisiológicas no estacionarias. Intenta superar las dificultades que a veces se presentan a la hora de estimar el exponente de Hurst en procesos que contienen tendencias o bien que no presentan una estructura auto{similar clara.

Metodología: el algoritmo para realizar el computo del DFA, se puede resumir en los siguientes pasos:

1. El primer paso es integrar la serie que estamos tratando, de forma que se pasa a tener una nueva serie que corresponde a un proceso auto-similar.

2. El siguiente paso es dividir la nueva serie temporal y(k) en segmentos de igual longitud, n. Para cada uno de estos segmentos se ajusta una recta por mínimos cuadrados. Esta recta representa la tendencia (trend) de cada segmento (ver Figura 9.2).

3. Seguidamente se elimina la tendencia (detrended) de la serie integrada y(k), substrayendo la tendencia local de cada segmento, que viene dada por yn(k).

4. La fluctuación, i.e., la raíz cuadrada de la serie integrada sin tendencia

5. Este cómputo se repite para todo tamaño de escalas, i.e., longitud de los segmentos, de manera, que se puede ofrecer una relación entre el promedio de las fluctuaciones como función de la longitud de los segmentos n, en el caso de la señal de HRV se establecerá relación con el número de latidos que caen dentro del segmento utilizado.

Propiedades: habitualmente el valor de F(n) crece con n de manera que es posible establecer una relación lineal, en escala logarítmica, lo que indica la presencia de un comportamiento de escalado o auto-semejante. El exponente de dicho escalado (la pendiente de la recta que mejor ajusta la relación lineal en escalas logarítmicas),permite caracterizar la estructura fractal. En la aplicación de este tipo de análisis a las señales de HRV se aprecia la existencia de al menos dos rectas con pendientes diferentes, lo que resulta ser un indicio de la de necesidad de utilizar dos exponentes para la caracterizaci_on fractal de la señal, uno que corresponde con longitudes de segmento pequeñas, es decir pocos latidos y por lo tanto caracteriza el comportamiento a corto plazo y otro que se asocia con longitudes de segmento mayores y por lo tanto caracteriza el comportamiento a largo plazo (ver Figura).

Figura: Exponentes resultantes de aplicar el DFA a una señal real.

Espectro Multifractal

Definición: las señales fisiológicas presentan estructura temporal fractal, siempre bajo condiciones de salud. Recientemente se ha comprobado que las uctuaciones del ritmo cardíaco también poseen características de lo que se conoce como multifractalidad, es decir, estas señales poseen estructura temporal multifractal. Grosso modo se puede decir que los conjuntos multifractales no utilizan un único patrón, que se repite a distintas escalas, sino varios. El motivo de analizar las señales cardíacas desde el punto de vista de la multifractalidad es que se ha comprobado que existe una pérdida de ésta cuando existen problemas cardíacos, esto es, la señal fisiológica que se analiza tras el problema cardíaco resulta ser menos multifractal, y por lo tanto, más monofractal.

Como se viene observando a lo largo de todo el estudio, lo que ocurre es que se produce una pérdida de complejidad en el sistema que controla el ritmo cardíaco, que se refleja, desde este nuevo punto de vista, en una disminución de la multifractalidad.

Metodología: se conocen como señales monofractales a aquellas que presentan las mismas propiedades a distintas escalas a través de toda la señal, por lo que las señales monofractales se pueden caracterizar por un único exponente global de Hurst, es decir, el exponente que calculemos para la señal describe el comportamiento de ésta a todas las escalas. Las señales multifractales, por contra, se pueden descomponer en distintos subconjuntos que estarán caracterizados, cada uno de ellos, por diferentes exponentes locales de Hurst. Estos exponentes cuantifican localmente el comportamiento singular de la señal y, por lo tanto, se relacionan con la escala local de la serie temporal.

De esta forma, para caracterizar completamente las señales multifractales se necesitarán varios exponentes de Hurst. Se pueden cuantificar las propiedades de los distintos subconjuntos de la señal, que se caracterizan por su exponente particular, mediante lo que se conoce como función D(h) que es la dimensión fractal, donde D(h0) es la dimensión fractal del subconjunto caracterizado por el exponente local de Hurst h0.

Propiedades: de los distintos resultados presentados en la literatura consultada se puede extraer la conclusión de que los pacientes con algún tipo de patología cardiovascular presentan un espectro multifractal reducido, más pequeño y con menor apertura entre sus brazos, se puede decir, que se aprecia una perdida de multifractalidad, es decir, una perdida en la complejidad del sistema que controla el sistema cardíaco.

MÉTODOS DE TEORÍA DE LA INFORMACIÓN

La presencia de dinámicas no lineares en las señales fisiológicas sugiere la aplicación de métodos adecuados a este dominio, entre los índices de este tipo propuestos, los métodos basados en el cálculo de la entropia de las señales están siendo ampliamente estudiados en los últimos años.

Los métodos basados en entropia cuantifican la irregularidad presente en las series temporales. El uso de estos métodos para cuantificar la irregularidad en señales cardíacas está motivado por las significativas diferencias encontradas en el nivel de irregularidad de estas señales dependiendo de los diversos estados de salud, lo cual puede reflejar importante información fisiológica.

A continuación se presentan los métodos más utilizados para el cálculo de la entropia en señales fisiológicas.

Entropía Aproximada

Definición: la entropía aproximada (Approximate Entropy, ApEn), es un estadístico que cuantifica el nivel irregularidad presente en una señal, fue propuesto por Pincus en 1991. La ApEn cuantifica la verosimilitud de que patrones similares de observaciones no sean seguidos por observaciones similares adicionales. Por lo tanto, series con patrones repetitivos tendrán un valor ApEn pequeño mientras que series menos predecibles tendrán un valor de ApEn mayor. Clínicamente valores relativamente bajos de ApEn pueden estar relacionados con alguna patología.

Metodología: para el cálculo de la ApEn es preciso definir dos parámetros a priori: la longitud de los vectores de datos a ser comparados, m, y el criterio de semejanza para la comparación entre dichos vectores, r. Dados N puntos u(i) de una señal, se define la ApEn(m; r) como sigue:

1. Se forman los vectores x(1); :::; x(N _m+1), definidos por: x(i) = [u(i); :::; u(i+m_1)] para i = 1; :::;N _ m + 1.

2. Se define la distancia entre los vectores x(i) y x(j), d[x(i); x(j)], como la diferencia máxima en valor absoluto, entre los correspondientes elementos de los vectores x(i) e x(j), esto es, d[x(i); x(j)] = máx k=1;:::;m (ju(i + k – 1) – u(j + k – 1)j)

3. Basada en esta distancia se define la siguiente medida de correlación,

 

 

 

 

 

Parámetros libres: la comparación de medidas de ApEn debe hacerse siempre con los parámetros m y r fijos. Para la aplicación en señales de HRV, en la bibliografía se recomienda utilizar m como 1, 2 o 3 [116, 90].

En cuanto al parámetro r Pincus recomienda utilizar valores entre el 10% y el 20% de la desviación estandar de cada señal [116, 117], obteniendo de esta forma una medida invariante en escala. Otros estudios revelan sin embargo, mejores resultados escogiendo un valor de r fijo [90, 83], independiente de la desviación estandar de cada señal. Este hecho puede deverse a que fijar r como un porcentaje de la desviación estandar puede tornar la medida sensible frente a datos atípicos.
Por otro lado, para filtrar el máximo ruido posible presente en los datos, r debe tener un valor superior al de la mayor parte del ruido.Típicamente el valor usado para el número de muestras N, varia entre 100 y 5000, no obstante, para obtener estimativas razonables es preferible no usar menos de 30m muestras.

Propiedades: las siguientes características de la ApEn hacen que su uso en procesamiento de señales biomédicas sea interesante: Obtención de una estimativa robusta incluso con un número reducido de muestras. Se trata de una medida robusta frente a interferencias transitorias bruscas y de corta duración (p.e., outliers). Se ve poco afectada por ruido de magnitud inferior a r. Puede ser aplicada en señales tanto deterministas como en señales estocásticas o en combinaciones de ambas.
De estas propiedades, las tres primeras la hacen apta para el análisis de series de datos cortas, contaminadas con ruido y/o obtenidas experimentalmente. La última propiedad es adecuada para el estudio de señales biológicas, ya que las salidas de los sistemas biológicos habitualmente poseen tanto componentes deterministas como aleatorias.

A pesar de las buenas características que presenta la ApEn para la caracterización de datos biológicos, esta posee también inconvenientes: Se trata de un estadístico sesgado, esto se debe a que en el algoritmo de cálculo se compara cada vector consigo mismo para evitar la ocurrencia del ln (0). Este hecho provoca que la ApEn sea fuertemente dependiente de la longitud de la secuencia de datos. Carece de consistencia relativa, es decir, que si la ApEn de un conjunto de datos es mayor que la de otro conjunto de datos, deberá permanecer mayor para cualquier condición (indeptendientemente de valor de m y r ), pero esta condición no se siempre se cumple.

Entropía Muestral

Definición: la entropía muestral (Sample Entropy, SampEn) se define como el negativo del logaritmo natural de la probabilidad condicionada de que dos secuencias similares para m puntos continúen siendo similares en el siguiente punto, donde las auto-comparaciones no están incluidas en el cálculo de las probabilidades.
Para estar definida la SampEn únicamente requiere que dos vectores parecidos para m observaciones continúen siendo parecidos para m + 1 observaciones. Este método surge para intentar mejorar el algoritmo frete a los problemas que presenta la ApEn, J.S Richman e J.R Moorman desarrollaron la familia de estadísticos SampEn, cuyas principales diferencias respecto a la ApEn son:

1. No realiza autocomparaciones de vectores.

2. Son considerados únicamente los primeros N _ m vectores de longitud m.

3. La medida de probabilidad es calculada directamente como el logaritmo de la probabilidad condicionada, en vez de como la razón de las sumas logarítmicas.

Metodología: el algoritmo de cálculo de la SampEn es el siguiente:

 

 

 

 

 

 

 

 

 

Propiedades: la SampEn consigue reducir el sesgo de los estadísticos y torna la medida más independiente de la longitud de los datos [124]. Además, presenta consistencia relativa en casos en los que la ApEn no lo hace, no obstante, no se puede afirmar que esta consistencia relativa se mantenga en todos los casos ya que la SampEn es en esencia, un estadístico que cuenta eventos, y si estos eventos son dispersos, el estadístico será inestable, lo que puede significar falta de consistencia relativa.

Entropía Multiescala

Definición: la entropía multiescala (Multiscale Entropy, MSE) propone un análisis de las series temporales fisiológicas en sus diferentes escalas espacio-temporales [30, 32].

Los algoritmos tradicionales basados en entropías cuantifican la regularidad de las series temporales, la entropía aumenta con el grado de irregularidad y es máxima para los sistemas aleatorios. Sin embargo, un incremento de la entropía no siempre está asociado a un incremento de la complejidad dinámica. Muchas patologías asociadas a un comportamiento más regular, proporcionan valores de entropía bajos en comparación con las dinámicas de estados de salud. No obstante, algunas patologías como la fibrilación auricular, están asociadas a fluctuaciones erráticas muy grandes con propiedades que
se asemejan a las del ruido icorrelacionado. Los algoritmos tradicionales estimaran, en estos casos, valores de entropía mayores para estas señales patológicas ruidosas que para las dinámicas de estados de salud.

Esta inconsistencia puede deberse a que estos algoritmos analizan las señales a una única escala y no tienen en cuenta las múltiples fluctuaciones temporales inherentes a los sistemas fisiológicos.

Metodología: dada una serie temporal discreta x1; :::; xi; :::; xN, se construye una serie temporal coarsegrained, siguiendo los siguientes pasos:

 

 

 

 

Propiedades: el método MSE requiere una longitud de datos adecuada para obtener estadísticas fiables de las medidas de entropía para cada escala. Típicamente se han utilizado series de 2 x 10 (4) para análisis en 20 escalas, de forma que la serie temporal coarse-grained mas pequeña tenga 1 x 10 (3) puntos.