Tutoriales PIC: Monitor del ritmo cardíaco

Vamos a volver a usar el sensor de luz TSL235 en una aplicación más interesante: monitorizar nuestro ritmo cardíaco. ¿Cómo hacerlo con un sensor de luz?

La idea viene explicada en una "application note" dedicada a estos sensores:

http://smartdata.usbid.com/datasheets/usbid/2000/2000-q2/230evm2taos.pdf

La configuración propuesta se muestra en la imagen adjunta. Un LED (u otra fuente de luz) ilumina un dedo, detrás del cuál se encuentra nuestro sensor. A cada pulsación del corazón, un nuevo flujo de sangre llega al dedo, cambiando su "opacidad" a la luz. El objetivo es detectar dichas oscilaciones, determinar su frecuencia fundamental y a partir de ahí el ritmo cardíaco en pulsaciones por minuto (ppm).

En la aplicación citada antes, el microcontrolador solo recoge los datos del sensor y los manda al PC para su proceso. La principal mejora de mi proyecto es que el PIC no sólo captura la señal, sino que la procesa y obtiene el ritmo cardíaco. De esta forma tenemos un sistema autónomo.

Al igual que en otras entradas el código del PIC se complementa con un programa MATLAB que nos permite ver gráficamente los datos sobre los que trabajamos y los resultados obtenidos. Recalcar que el programa MATLAB sólo visualiza los datos: todos los cálculos se llevan a cabo en el PIC.

En este corto video podemos visualizar los datos obtenidos tras un pre-proceso de los datos del sensor. Es sobre esta señal sobre la que el PIC trabaja para estima su periodo:

Archivos de código asociados a esta entrada:

PIC: pulse_monitor2.c

MATLAB: heart_monitor.m

Hardware (TSL235R):

Vamos a usar el mismo sensor de luz (TSL235R) que presentamos en la entrada anterior. Entonces usábamos el modo CAPTURE para medir la duración de un periodo (frecuencia instantánea). En este proyecto, para evitar los efectos de aliasing que vimos entonces vamos a usar el otro enfoque. En vez de medir la duración de un periodo, vamos a contar el número de pulsos en recibidos en un intervalo más grande. Obtendremos así la frecuencia media y por lo tanto la luminosidad media durante el intervalo. Para contar los pulsos usaremos la interrupción INT0 (igual que con los codificadores en cuadratura). Para poder usar INTO conectaremos la salida del sensor al pin RB0 (ver figura).

Respecto a la configuración LUZ / DEDO / SENSOR, he decidido no usar un LED sino la fuente de luz (natural o artificial) de la habitación (ver figura adjunta). Obviamente esto requiere que haya suficiente luz ambiente para que llegue algo de luz al sensor.

La opción del LED es preferible si se desea construir un dispositivo (por ejemplo para fijarlo en un dedo) y que no dependa del nivel de luz ambiente. La ventaja de mi enfoque es simplemente que es más sencillo de usar en las pruebas: basta tapar el sensor (colocado plano sobre una breadboard) con el dedo.

De todas formas, de cara al software, ambas configuraciones pueden usarse de forma similar.

Software:

Como hemos descrito el hardware empleado se reduce al sensor, por lo que el grueso del trabajo recaerá en el programa.

En primer lugar debemos contar con una ISR para la interrupción INT0 que se incremente cada vez que se reciba un pulso:

uint16 pulsos=0; // Pulse counter

// High priority interruption

#pragma interrupt high_ISR

void high_ISR (void)

{

if (INT0_flag)

{

pulsos++;

INT0_flag=0;

}

// Code @ 0x0008 -> Jump to ISR for high priority ints

#pragma code high_vector = 0x0008

void high_interrupt (void){_asm goto high_ISR _endasm}

#pragma code

La variable pulsos (uint16) acumula los pulsos recibidos. Cada cierto tiempo (usando timer TMR3) consulto el número de pulsos recibidos. Programaré el contador del timer TMR3 para saltar cada 50 milisegundos. La elección de dicho espaciado se basa en varios criterios:

1. Como la frecuencia máxima del sensor es de 1 MHz podríamos llegar a tener un máximo de 50000 pulsos durante esos 50 msec. Eso asegura que el contador de pulsos no se desborda.

2. Un periodo de 50 mseg es un tiempo suficientemente largo para no ver efectos de aliasing si usamos luz artificial. Estamos midiendo la "luz media" durante 5 ciclos de encendido-apagado.

3. Un espaciado de 50 msec corresponde a una frecuencia de muestreo de (1000/50) 20 muestras por segundo. Valores "normales" del ritmo cardíaco oscilan entre 30 y 150 pulsaciones por minuto. Dividiendo por 60 (pulsaciones por segundo) corresponden a un rango entre 0.5 y 2.5 Hz. Una frecuencia de muestreo de 20 Hz es más que suficiente para captarlas adecuadamente.

En la interrupción del TMR3 consultamos el número de pulsos recibidos y reseteamos el contador. El número de pulsos recibidos constituye nuestra medida del sensor. A partir de ahí pasamos a procesar los datos.

Procesado de señal:

Nuestro objetivo en el proyecto es detectar y medir una oscilación regular con una frecuencia de entre 0.5 y 3 Hz (si tu pulso está por debajo de 30 ppm o por encima de 180 ppm deberías hacertelo mirar).

La señal medida (luz transmitida a través del dedo) además de esa oscilación fundamental tendrá otras componentes:

frecuencias bajas, oscilaciones lentas del nivel de luz, tal vez causadas por un pequeño movimiento del dedo)
frecuencias más altas (ruido en la medida del sensor).

Para mejorar los resultados posteriores vamos a hacer un preproceso inicial de los datos, consistente en un filtrado paso-banda, ajustado para tratar de eliminar frecuencias fuera de la zona de interés (0.5-3 Hz = 30-180 ppm).

Por simplicidad se ha implementado con un filtro FIR (Finite Impulse Response) con coeficientes enteros. Esencialmente un filtro FIR consiste en que la salida y[n] es combinación de unas cuantas entradas anteriores. En este caso he usado las 8 entradas anteriores:

y[n] = sum_k {b(k) x[n-k]} k=0,1,...,7

Los coeficientes usados son b(k) = {1 2 2 1 -1 -2 -2 -1}. Por eficiencia (para poder usar aritmética entera) hemos usado coeficientes enteros. De todas formas, la función de transferencia de este filtro se corresponde bastante con lo que quiero:

Se observa que en el intervalo [1,3] Hz = [60,180] ppm la atenuación es pequeña. Por el contrario, frecuencias fuera de esa banda son notablemente atenuadas.

Para implementar el filtrado anterior en el PIC uso un buffer circular last_x[], con las últimas 8 medidas. En cada paso combino los resultados del buffer para obtener la salida y[n]

// Band-pass filtering b(k) = [1 2 2 1 -1 -2 -2 -1]

next_x++; next_x&=7; last_x[next_x]=res;

p=next_x+1; p&=7; res-=last_x[p];

p=next_x-1; p&=7; res+=(last_x[p]<<1);

p=next_x+2; p&=7; res-=(last_x[p]<<1);

p=next_x-2; p&=7; res+=(last_x[p]<<1);

p=next_x+3; p&=7; res-=(last_x[p]<<1);

p=next_x-3; p&=7; res+=last_x[p];

p=next_x+4; p&=7; res-=last_x[p];

Al final en la variable res está el valor de la salida y[n]. En la figura adjunta se muestran ambas señales (x[n] e y[n]) tomadas durante un periodo de 8 segundos:

En la parte de arriba tenemos los datos obtenidos directamente del sensor (x[n]). Se observa la oscilación principal (unos 8 ciclos en esos 8 segundos), pero también la contaminación con oscilaciones lentas de luz (caída) más un ruido de alta frecuencia (pequeños picos).

En la parte de abajo tenemos la salida del filtro paso-banda, y[n]. La oscilación en la que estamos interesados queda claramente resaltada. Sobre esta señal será mucho más sencillo medir la duración de un periodo.

Detección del ritmo cardíaco

Si mandamos los datos (filtrados) al PC sería fácil en MATLAB (analizando cierto número de muestras) detectar la frecuencia predominante. Por ejemplo podríamos usar análisis de Fourier, a través de la conocida FFT (Fast Fourier Transform).

Sin embargo estoy interesado en que sea el PIC quien haga los cálculos. Aunque sería posible programar una FFT en el PIC, veremos que no es necesario, ya que puede hacerse de una forma más sencilla de programar y computacionalmente menos costosa.

En este proyecto voy a usar la auto-correlación. La auto-correlación de una señal (o correlación consigo misma) se define como:

C[n] = sum_k {y[k] y[k-n]}

C[n] mide la similitud de las muestras de una señal separadas por n muestras. La idea es que si y[k] e y[k-n] son similares su producto será siempre positivo (+ x + o – x - = +) y la suma de dichos valores tenderá a ser alta. Si y[k] e y[k-n] no están relacionadas, su producto será unas veces positivo y otras negativo y su suma tenderá a cancelarse.

En la imagen adjunta tenemos la señal recogida por el sensor y su autocorrelación (abajo) desde n=0 hasta n=100.

El máximo de la autocorrelación lo encontramos para n=0, ya que nada se parece más a una señal que ella misma. En el caso de señales cuasi-periódicas como las que nos ocupan, también encontraremos un máximo (un poco menor) para n = T (periodo de la señal). La idea es que pasadas T muestras la señal se repite.

En este caso se puede observar que T es del orden de 18-19 muestras.

Hay que tener en cuenta que la correlación también muestra máximos en 2T, 3T, etc. Tendremos que tener cuidado con esto.

Para calcular la autocorrelación necesito guardar unas cuantas muestras del pasado y[k-n], para "comparar" con la muestra más reciente y[n]. Iré guardando en un buffer las muestras obtenidas. En cada paso descartaré la más antigua y guardaré la más reciente, teniendo siempre disponibles los N últimos valores:

El tamaño de este buffer debe ser el suficiente para capturar la correlación del T más grande en el que pudiera estar interesado. He usado N_past=32 muestras que permiten detectar una correlación con T = 32 muestras = 32/20 = 1.6 sec. Esto corresponde a f = 1/1.6 = 0.625 Hz = 37 ppm.

Esta vez, en vez de un buffer circular uso un buffer de desplazamiento. A cada paso corro las N-1 muestras y hago hueco para guradar la última recibida en la posición [0]:

for(p=N_past;p>=1;p--) past[p]=past[p-1]; // shift buffer

past[0]=res; // last sample

Esto es poco eficaz, pero al ritmo de muestras recibidas (20 por segundo) no tengo problemas. La ventaja es que, al no ser circular, será más fácil (más legible en el código) direccionar las muestras pasadas:

past[0] = y[k] (última) past[1]=y[k-1] , . . . , past[n]=y[k-n]

Busqueda inicial:

En un primer momento (fase de búsqueda) no tenemos ni idea de nuestro ritmo cardíaco, por lo que debemos calcular todas las correlaciones en un cierto margen. Los resultados se guardan en la variable corr[].

Debo decidir que rango de correlaciones (n) debo calcular. No necesito calcular correlaciones muy próximas porque no me interesa el máximo que se que existe para n=0. He decidido buscar entre n=5 y n=30 muestras, correspondientes a ritmos entre 240 ppm (T=5/20) y 40 ppm (T=30/20).

En el programa un par de #defines determinan esta elección:

// Mem for complete correlation lag=delta,delta+1, ..., delta+N_lag-1

#define N_lag 25

#define delta 5

int16 corr[N_lag]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};

N_lag indica el número de correlaciones a estimar (25) y delta el punto de inicio (5)

Si tu pulso está por debajo de 40 ppm deberías extender el tamaño de Na (p.e a 30). Si está por encima de 240 deberías acudir a un hospital (y solo después, disminuir delta).

Durante un cierto número de muestras (N_corr=30) se van acumulando en corr[] los valores de la auto-correlación. Al cumplirse el tiempo, una bandera (corr_ready) se pone a 1 para indicar al programa principal que las correlaciones ya están calculadas.

// Accumulate correlations

for(i=0,p=delta;i<N_lag;i++,p++) corr[i]+= (res*past[p]);

if (n_sample==NT_corr)

{

T = center; T*=FRAC_SCALE; T+=frac;

#ifdef MATLAB_MODE

send_msg(T,254);

#endif

n_sample=0; corr_ready=1;

}

En este código:

1. Se calculan las correlaciones desde p=delta a delta+N_lag sumando y[n]*y[n-p] a los correladores.

2. Cuando llevamos acumuladas NT_corr (por defecto 30) muestras, ponemos corr_ready a 1, para indicar que la correlación está lista y podemos pasar a calcular su máximo, etc. (en programa principal). En su caso, también enviamos el periodo hallado al PC a través del puerto serie.

En el programa principal, si las correlaciones están listas pasamos a buscar su máximo y a tomar algunas decisiones:

if (corr_ready==1) //

{

M=find_max(); cont++;

if ( (M>0) && (abs8(M-last_max)<=1)) aciertos++;

last_max=M;

if(M>0) center=M+delta; else center=0; frac=0;

if (cont==3)

{

if (aciertos==3) tracking=1;

cont=0; aciertos=0;

}

corr_ready=0; for(k=0;k<N_lag;k++) corr[k]=0;

}

El código anterior realiza las siguientes funciones:

1. Halla el máximo más cercano a cero de las correlaciones y (caso de encontrarlo) lo guarda en la variable center. Si no lo encuentra devuelve 0. Un máximo se define como un punto que supera un cierto umbral y que es mayor que sus 4 vecinos (dos a cada lado)

2. Se verifica si el máximo encontrado no se ha movido mucho respecto al anterior (aciertos++). Si eso pasa varias veces seguidas, confiamos en que hemos encontrado un un pico de correlación fiable. En ese caso se declara tracking=1 y entramos en la fase de tracking.

3. Se pone corr_ready a 0 y se resetean las correlaciones para volver a calcularlas.

Tracking:

Si sistemáticamente encontramos el máximo de correlaciones en la misma muestra pasamos a la siguiente fase: el seguimiento fino o tracking. La idea es que partimos de una buena idea inicial para el periodo T (el valor estimado en la fase de búsqueda). La precisión de esta aproximación es de aproximadamente 1 muestra.

Imaginad una situación con un pulso de 80 ppm. Para esa frecuencia el periodo es de 0.75 segundos, o lo que es lo mismo, 15 muestras. En ese caso el error causado por 1 muestra correspondería a

T = 14 -> f = 85.7 ppm

T = 16 -> f = 75 ppm

Tenemos saltos del orden de 5 ppm en la estimación del ritmo cardíaco, lo que claramente es mucho.

El objetivo en esta fase de tracking es mejorar la estimación del ritmo cardíaco y permitir un rápido seguimiento de posibles cambios.

El sistema que vamos a usar y que describo a continuación es similar al usado en un receptor GPS cuando está midiendo la distancia a los satélites.

En esta fase se calcularan únicamente tres correlaciones (durante NT_track= unas 10 muestras):

· corr[0], correspondiente al periodo ideal T

· corr[1] y corr[2], correspondiente a T+1 y T-1 muestras respectivamente.

La idea es que si la correlación máxima se da en T los valores de correlación en T-1 y T+1 deberían ser iguales.

En la gráfica adjunta se ve que el máximo está en T=18. Se aprecia que los valores de correlación en T-1=17 y T+1=19 son muy parecidos entre sí, lo que indica que T=18 es una buena estimación del máximo.

Sin embargo, si nos fijamos con cuidado vemos que corr[17] un poco mayor que corr[19], indicando que el máximo de correlación deberíamos moverlo un poco a la izquierda. La idea es que si las correlaciones en T-1 y T+1 no son iguales, moveremos ligeramente el periodo ideal T hacia el lado que muestra una correlación mayor que el otro.

En principio podría parecer que mover T supondría añadir o restar una muestra. Sin embargo, como queremos tener más precisión que la que nos da 1 muestra, vamos a permitir que T tome valores fraccionarios. El valor de T se codificará con dos variables de tipo uint8. La primera (center) contendrá la parte entera y la segunda (frac) la parte fraccionaria (escalada por un factor FRAC_SCALE=32).

Por ejemplo, si center=17 y frac = 7 el periodo de correlación es:

T = center + (frac/32) = 17 + 7/32 = 17.22 muestras

Con este enfoque la precisión con la que podemos especificar el periodo es de 1/FRAC_SCALE = 1/32 muestras, lo que (para un T del orden de 15) no da una precisión de unos 0.2 ppm, más que suficiente en nuestra aplicación.

Al usar valores de T no enteros nos surge otro problema. Para calcular la correlación tendría que comparar la muestra recibida y[k] con p.e. y[k-17.22]. Obviamente yo no puedo preguntar por past[17.22]. La solución es usar una simple interpolación lineal y estimar past[17.22] en función de past[17] y past[18]. De forma general:

y[n-T]=past[T]= past[center+frac/32]

= past[center]*(1-frac/32) + past[center+1](frac/32)

= ( past[center]*(32-frac) + past[center+1]*frac )/32

De esta forma podemos hacer los cálculos en aritmética entera.

La función interp_past interpola el array past[] para el retraso T determinado por (center,frac):

// Interpolates data @ c+f/FRAC_SCALE using past[c] y past[c+1]

int16 interp_past(uint8 c,uint8 f)

{

int16 temp;

temp=past[c]*(FRAC_SCALE-f)+past[c+1]*f;

temp/=FRAC_SCALE;

return temp;

}

El código, dentro de la interrupción del timer, para el caso de que tracking=1 es:

corr[0]+= res*interp_past(center,frac);

corr[1]+= res*interp_past(center+1,frac);

corr[2]+= res*interp_past(center-1,frac);

// Steers tracking based on acums

if(n_sample==NT_track)

{

tracking=steer_tracking();

n_sample=0;

}

Se calculan las tres correlaciones y cada cierto tiempo (NT_track muestras) se usan para modificar la estimación del máximo de correlación T (que está codificado en center + frac). La función steer_tracking se encarga de ello:

uint8 steer_tracking(void)

{

float c;

uint8 res=1;

int16 T;

c = (abs16(corr[1])+ abs16(corr[2]));

if (c==0) res=0;

else

{

c = (corr[1]-corr[2])/c; c*=FACTOR_STEERING;

if (c>0) c+=0.5; if(c<0) c-=0.5; // Rounding

}

frac+=(int8)c;

T=(int16)c; send_msg(T,250);

send_msg(corr[0],0);

send_msg(corr[1],1);

send_msg(corr[2],2);

while(frac>=FRAC_SCALE) { frac-=FRAC_SCALE; center++; }

while(frac<0) { frac+=FRAC_SCALE; center--; }

if ( (center>=N_lag) || (center<8) ) res=0; // Lost tracking

// checks if tracking is able to follow changes

if ( (corr[0]>corr[1]) && (corr[0]>corr[2]) ) bad=0; else bad++;

PORTD=bad; //if (bad==6) res=0;

corr[0]=0;corr[1]=0;corr[2]=0; // Reset acums

return res;

}

La corrección usada es: (corr(T+1) – corr(T-1))/(|corr(T+1)|+ |corr(T-1)|).

Este valor está normalizado entre [-1 1] y se multiplica por FACTOR_STEERING para obtener el cambio que aplicaremos a la parte fraccionaria del periodo (frac). FACTOR_STEERING es un #define que marca el máximo cambio que es capaz de hacer el algoritmo sobre frac (estamos usando un valor 8).

Dentro de la función también comprobamos si hay problemas para seguir la señal. La función devuelve 0 si detecta que se ha perdido el tracking. En ese caso se vuelve a la fase anterior de búsqueda, basada en encontrar máximo de la función de correlación.

Resultados:

Como en proyectos anteriores los datos se comunican al PC a través del puerto serie.

Se puede elegir entre una interfaz de texto (#define TEXT_MODE) o a través de un programa de MATLAB (#define MATLAB_MODE)para visualizar datos.

En el modo texto simplemente se vuelca el status (NO DATA, SEARCH o TRACK) y en caso de tener datos la estimación del ritmo cardíaco. El modo NO DATA se define cuando los valores (número de pulsos) recibidos por el sensor son demasiado altos, indicando que no hay un dedo encima del sensor:

if (res>1000) res=last_x[next_x]; // Way too much light: no finger placed on sensor

Si intentáis hacer este proyecto, dependiendo del nivel de luz con el que esteis trabajando, tendréis que modificar ese valor.

El video final muestra los resultados finales del proyecto usando la interfaz de MATLAB (heart_monitor.m) para visualizar los datos:

De nuevo, resaltar que el programa MATLAB solo se ocupa de presentar de una forma más vistosa los datos, status y resultados del proceso. Todos los cálculos se hacen en el PIC y podríamos configurarlo como un sistema autónomo volcando los resultados a través de un LCD.