Monitor de frequência cardíaca Arduino usando MAX30102 e oximetria de pulso

Tempo de leitura: 10 minutes

A oximetria de pulso monitora a saturação de oxigênio no sangue medindo a magnitude da luz vermelha e infravermelha refletida [leia mais sobre oximetria de pulso aqui e aqui]. Os oxímetro de pulso também podem aproximar a freqüência cardíaca analisando a resposta da série de tempo da luz vermelha e infravermelha refletida. O oxímetro de pulso MAX30102 é um sensor compatível com Arduino e barato que permite o cálculo da frequência cardíaca usando o método descrito acima. Neste tutorial, o sensor MAX30102 será apresentado junto com várias análises aprofundadas dos dados de reflexão de vermelho e infravermelho que serão usados para calcular parâmetros como frequência cardíaca e saturação de oxigênio no sangue.

 

Lista de peças e Ligação

Os principais componentes necessários para este tutorial são o oxímetro de pulso MAX30102 e um microcontrolador Arduino. Também estarei usando um Raspberry Pi para ler dados seriais de alta velocidade impressos pelo Arduino. Estou usando o Raspberry Pi porque analisarei os dados de pulso com programas Python robustos e bibliotecas que não estão disponíveis na plataforma Arduino. A lista de peças para os experimentos é mostrada abaixo:

  1. Oxímetro de pulso MAX30102
  2. Arduino Uno
  3. Raspberry Pi + cartão SD de 8 GB
  4. Fios de jumpers
  5. Mini Protoboard

O MAX30102 usa comunicação I2C de dois fios para fazer interface com a placa Arduino Uno. Eu uso as portas I2C nas portas A4/A5 na placa Arduino. A ligação é mostrada abaixo:

 

MAX30102 Arduino Library

O Sparkfun possui uma biblioteca que lida com a comunicação entre o Arduino e o MAX30102. Estaremos usando a biblioteca Sparkfun para lidar com a leitura em alta velocidade dos dados de refletância de infravermelho e vermelho.

No IDE do Arduino:

Vá para Sketch -> Incluir Biblioteca -> Gerenciar Bibliotecas

Digite “max30” na barra de pesquisa

Baixe a “Sparkfun MAX3010x Pulse and Proximity Sensor Library”

Uma configuração simples de alta velocidade baseada no Arduino Uno é mostrada abaixo. Ele obtém amostras do MAX30102 a 400 Hz e imprime na porta serial. A 400 Hz e uma taxa de transmissão serial de 115200, o Raspberry Pi é capaz de ler cada ponto de dados sem problemas.

#include <Wire.h>
#include "MAX30105.h"

MAX30105 particleSensor; // inicializa MAX30102 com I2C

void setup()
{
  Serial.begin(115200);
  while(!Serial);// Devemos esperar que Teensy fique online
  delay(100);
  Serial.println("");
  Serial.println("MAX30102");
  Serial.println("");
  delay(100);
  // Initialize sensor
  if (particleSensor.begin(Wire, I2C_SPEED_FAST) == false) // Use a porta I2C padrão, velocidade de 400 kHz
  {
    Serial.println("MAX30105 não foi encontrado. Verifique a ligação/alimentação.");
    while (1);
  }

  byte ledBrightness = 70; //Options: 0=Off to 255=50mA
  byte sampleAverage = 1; //Options: 1, 2, 4, 8, 16, 32
  byte ledMode = 2; //Options: 1 = Red only, 2 = Red + IR, 3 = Red + IR + Green
  int sampleRate = 400; //Options: 50, 100, 200, 400, 800, 1000, 1600, 3200
  int pulseWidth = 69; //Options: 69, 118, 215, 411
  int adcRange = 16384; //Options: 2048, 4096, 8192, 16384

  particleSensor.setup(ledBrightness, sampleAverage, ledMode, sampleRate, pulseWidth, adcRange); //Configure sensor with these settings
}

void loop() {
  particleSensor.check(); //Check the sensor
  while (particleSensor.available()) {
      // ler IR armazenado
      Serial.print(particleSensor.getFIFOIR());
      Serial.print(",");
      // ler o vermelho armazenado
      Serial.println(particleSensor.getFIFORed());
      // ler o próximo conjunto de amostras
      particleSensor.nextSample();      
  }
}

Com um dedo preso ao MAX30102 (seja por elástico, fita ou encapulação), a impressão para o plotter serial Arduino deve ser o seguinte:

Não precisamos nos preocupar com a incapacidade de rastrear a forma de cada gráfico, pois vamos lê-los em Python usando o leitor serial e analisar totalmente os dados de vermelho e infravermelho do sensor MAX30102. Na próxima seção, estaremos lendo os dados de 400 Hz em tempo real em Python e analisando os dados usando vários algoritmos complexos que variam de análise de domínio de frequência e análise de wavelet.

 

Salvando dados MAX30102 com Python

Podemos ler os dados de saída serial do Arduino no Raspberry Pi usando a biblioteca serial do Python. No código do Arduino acima, a única mudança que precisamos fazer é adicionar uma impressão da função ‘micros()’ para anexar um carimbo de data/hora às leituras de dados de valores de refletividade de infravermelho e vermelho. O código do Arduino é mostrado abaixo:

#include "MAX30105.h"

MAX30105 particleSensor; // inicializar MAX30102 com I2C

void setup()
{
  Serial.begin(115200);
  while(!Serial); //Devemos esperar que Teensy fique online
  delay(100);
  Serial.println("");
  Serial.println("MAX30102");
  delay(100);
  // Inicializar sensor
  if (particleSensor.begin(Wire, I2C_SPEED_FAST) == false) //Use a porta I2C padrão, velocidade de 400 kHz
  {
    Serial.println("MAX30105 was not found. Please check wiring/power. ");
    while (1);
  }

  byte ledBrightness = 70; //Options: 0=Off to 255=50mA
  byte sampleAverage = 1; //Options: 1, 2, 4, 8, 16, 32
  byte ledMode = 2; //Options: 1 = Red only, 2 = Red + IR, 3 = Red + IR + Green
  int sampleRate = 400; //Options: 50, 100, 200, 400, 800, 1000, 1600, 3200
  int pulseWidth = 69; //Options: 69, 118, 215, 411
  int adcRange = 16384; //Options: 2048, 4096, 8192, 16384

  particleSensor.setup(ledBrightness, sampleAverage, ledMode, sampleRate, pulseWidth, adcRange); //Configure sensor with these settings
}

void loop() {
  particleSensor.check(); // Verifique o sensor
  while (particleSensor.available()) {
      // ler IR armazenado
      Serial.print(micros());
      Serial.print(",");
      Serial.print(particleSensor.getFIFOIR());
      Serial.print(",");
      // ler o vermelho armazenado
      Serial.println(particleSensor.getFIFORed());
      // ler o próximo conjunto de amostras
      particleSensor.nextSample();      
  }
}

Um algoritmo de leitura serial de alta velocidade para Python é mostrado abaixo para ler os valores de impressão serial do Arduino. O código Python irá adquirir os dados e salvá-los em um arquivo .csv para análise posterior. O motivo pelo qual os salvamos imediatamente é que os dados chegam em velocidades tão altas que queremos minimizar a quantidade de processamento feita entre as aquisições em série. Este código para salvar os dados de impressão do Arduino em Python é mostrado abaixo:

import serial,time,csv,os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

plt.style.use('ggplot')

## inicializar a porta serial (ttyUSB0 ou ttyACM0) a 115200 baud rate
ser = serial.Serial('/dev/ttyUSB0',
                    baudrate=115200)
## set filename and delete it if it already exists
datafile_name = 'test_data.csv'
if os.path.isfile(datafile_name):
    os.remove(datafile_name)

## looping através de impressões seriais e aguarde a reinicialização do Arduino Uno
## com palavra inicial "MAX30102"
all_data = []
start_word = False
while True:
    try:        
        curr_line = ser.readline() # leia linha
        if start_word == False:
            if curr_line[0:-2]==b'MAX30102':
                start_word = True
                print("Início do programa")
                continue
            else:
                continue
        all_data.append(curr_line) # append to data vector
    except KeyboardInterrupt:
        break   
print("Saiu do Loop")

# percorrer o vetor de dados e remover dados inválidos
# então, crie vetores para as variáveis de tempo, vermelho e IV
t_vec,ir_vec,red_vec = [],[],[]
ir_prev,red_prev = 0.0,0.0
for ii in range(3,len(all_data)):
    try:
        curr_data = (all_data[ii][0:-2]).decode("utf-8").split(',')
    except:
        continue
    
    if len(curr_data)==3:

        if abs((float(curr_data[1])-ir_prev)/float(curr_data[1]))>1.01 or\
        abs((float(curr_data[2])-red_prev)/float(curr_data[2]))>1.01:
            continue
        
        t_vec.append(float(curr_data[0])/1000000.0)
        ir_vec.append(float(curr_data[1]))
        red_vec.append(float(curr_data[2]))
        ir_prev = float(curr_data[1])
        red_prev = float(curr_data[2])

print('Taxa de amostragem: {0:2.1f}Hz'.format(1.0/np.mean(np.abs(np.diff(t_vec)))))

## salvando dados
with open(datafile_name,'a') as f:
    writer = csv.writer(f,delimiter=',')
    for t,x,y in zip(t_vec,ir_vec,red_vec):
        writer.writerow([t,x,y])
        
## plotagem de vetores de dados
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(111)
ax1.set_xlabel('Time [s]',fontsize=24)
ax1.set_ylabel('IR Amplitude',fontsize=24,color='#CE445D',labelpad=10)
ax1.tick_params(axis='both',which='major',labelsize=16)
plt1 = ax1.plot(t_vec,ir_vec,label='IR',color='#CE445D',linewidth=4)
ax1_2 = plt.twinx()
ax1_2.grid('off')
ax1_2.set_ylabel('Red Amplitude',fontsize=24,color='#37A490',labelpad=10)
ax1_2.tick_params(axis='y',which='major',labelsize=16)
plt2 = ax1_2.plot(t_vec,red_vec,label='Red',color='#37A490',linewidth=4)
lns = plt1+plt2
labels = [l.get_label() for l in lns]
ax1_2.legend(lns,labels,fontsize=16)
plt.xlim([t_vec[0],t_vec[-1]])
plt.tight_layout(pad=1.2)
plt.savefig('max30102_python_example.png',dpi=300,facecolor=[252/255,252/255,252/255])
plt.show()

O gráfico final produzido pelo código Python acima deve ser o seguinte:

Na próxima seção, explorarei vários métodos para analisar os dados da série temporal mostrados no gráfico acima. Também discutirei algumas análises de domínio de frequência e wavelet para determinar a periodicidade dos pulsos para aproximação da frequência cardíaca e oxigenação do sangue.

 

Transformação rápida de Fourier para frequência cardíaca aproximada

Abaixo estão alguns links para pesquisas científicas realizadas sobre oximetria de pulso e a relação com a oxigenação no sistema circulatório:

Nesta seção, apresentarei a relação básica entre os valores de refletância de vermelho e infravermelho medidos pelo oxímetro de pulso MAX30102 e a frequência cardíaca. Na publicação intitulada “Oximetria de pulso livre de calibração baseada em dois comprimentos de onda no infravermelho – um estudo preliminar”, a figura a seguir é apresentada como um pulso PPG onde a luz transmitida através do tecido diminui durante um evento chamado sístole (o o coração se contrai e bombeia sangue de suas câmaras para as artérias) e aumenta durante a diástole (o coração relaxa e suas câmaras se enchem de sangue).

De: Oximetria de pulso livre de calibração baseada em dois comprimentos de onda no infravermelho – um estudo preliminar: O pulso PPG. A luz transmitida através do tecido diminui durante a sístole e aumenta durante a diástole. ID e IS representam a transmissão máxima e mínima de luz através do tecido.

Se fôssemos aumentar o zoom em um de nossos pulsos MAX30102, veríamos quase exatamente o mesmo perfil nas respostas de vermelho e infravermelho:

Podemos usar esse comportamento cíclico para aproximar o intervalo entre as “batidas” do coração para determinar a frequência cardíaca aproximada de um indivíduo. A maneira mais simples de calcular a freqüência cardíaca é registrar alguns segundos de dados de refletância de infravermelho ou vermelho e calcular o conteúdo da freqüência dominante do sinal. Se usarmos a transformada rápida de Fourier de Python (FFT em Numpy), o pico da FFT se aproxima da frequência do ciclo de contração e relaxamento do coração – o que chamamos de frequência cardíaca.

A maneira mais simples de calcular a freqüência cardíaca é registrar alguns segundos de dados de refletância de infravermelho ou vermelho e calcular o conteúdo da freqüência dominante do sinal.

O código e o gráfico abaixo mostram o método FFT para aproximar a freqüência cardíaca para uma amostra de 9 segundos de dados MAX30102.

import csv
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')

## reading data saved in .csv file
t_vec,ir_vec,red_vec = [],[],[]
with open('test_data.csv',newline='') as csvfile:
    csvreader = csv.reader(csvfile,delimiter=',')
    for row in csvreader:
        t_vec.append(float(row[0]))
        ir_vec.append(float(row[1]))
        red_vec.append(float(row[2]))

s1 = 0 # change this for different range of data
s2 = len(t_vec) # change this for ending range of data
t_vec = np.array(t_vec[s1:s2])
ir_vec = ir_vec[s1:s2]
red_vec = red_vec[s1:s2]

# sample rate and heart rate ranges
samp_rate = 1/np.mean(np.diff(t_vec)) # average sample rate for determining peaks
heart_rate_range = [0,250] # BPM
heart_rate_range_hz = np.divide(heart_rate_range,60.0)
max_time_bw_samps = 1/heart_rate_range_hz[1] # max seconds between beats
max_pts_bw_samps = max_time_bw_samps*samp_rate # max points between beats

## plotting time series data
fig = plt.figure(figsize=(14,8))
ax1 = fig.add_subplot(111)
ax1.set_xlabel('Time [s]',fontsize=24)
ax1.set_ylabel('IR Amplitude',fontsize=24,color='#CE445D',labelpad=10)
ax1.tick_params(axis='both',which='major',labelsize=16)
plt1 = ax1.plot(t_vec,ir_vec,label='IR',color='#CE445D',linewidth=4)
ax1_2 = plt.twinx()
ax1_2.grid('off')
ax1_2.set_ylabel('Red Amplitude',fontsize=24,color='#37A490',labelpad=10)
ax1_2.tick_params(axis='y',which='major',labelsize=16)
plt2 = ax1_2.plot(t_vec,red_vec,label='Red',color='#37A490',linewidth=4)
lns = plt1+plt2
labels = [l.get_label() for l in lns]
ax1_2.legend(lns,labels,fontsize=16,loc='upper center')
plt.xlim([t_vec[0],t_vec[-1]])
plt.tight_layout(pad=1.2)

## FFT and plotting frequency spectrum of data
f_vec = np.arange(0,int(len(t_vec)/2))*(samp_rate/(len(t_vec)))
f_vec = f_vec*60
fft_var = np.fft.fft(red_vec)
fft_var = np.append(np.abs(fft_var[0]),2.0*np.abs(fft_var[1:int(len(fft_var)/2)]),
                    np.abs(fft_var[int(len(fft_var)/2)]))

bpm_max_loc = np.argmin(np.abs(f_vec-heart_rate_range[1]))
f_step = 1
f_max_loc = np.argmax(fft_var[f_step:bpm_max_loc])+f_step
print('BPM: {0:2.1f}'.format(f_vec[f_max_loc]))
fig2 = plt.figure(figsize=(14,8))
ax2 = fig2.add_subplot(111)
ax2.loglog(f_vec,fft_var,color=[50/255,108/255,136/255],linewidth=4)
ax2.set_xlim([0,f_vec[-1]])
ax2.set_ylim([np.min(fft_var)-np.std(fft_var),np.max(fft_var)])
ax2.tick_params(axis='both',which='major',labelsize=16)
ax2.set_xlabel('Frequency [BPM]',fontsize=24)
ax2.set_ylabel('Amplitude',fontsize=24)
ax2.annotate('Heart Rate: {0:2.0f} BPM'.format(f_vec[f_max_loc]),
             xy = (f_vec[f_max_loc],fft_var[f_max_loc]+(np.std(fft_var)/10)),xytext=(-10,70),
             textcoords='offset points',arrowprops=dict(facecolor='k'),
             fontsize=16,horizontalalignment='center')
fig2.savefig('max30102_fft_heart_rate.png',dpi=300,facecolor=[252/255,252/255,252/255])
plt.show()

De minha própria experiência com este método, recomendo pelo menos 8 segundos de dados de refletividade para calcular a frequência cardíaca verdadeira. Abaixo disso, o FFT começará a ver conteúdos de diferentes frequências no sinal. Isso poderia ser evitado por meio de filtragem, mas não incluí isso aqui.

 

Aplicação de uma aproximação de gradiente para encontrar a frequência cardíaca

A dificuldade em usar o FFT para cálculo da freqüência cardíaca é o número necessário de ciclos. Vários ciclos são necessários para uma aproximação de frequência precisa. Portanto, outro método é apresentado aqui, o qual usa uma função gradiente de segunda ordem para aproximar a taxa de mudança do pulso. Como o ponto mais acentuado durante o ciclo circulatório é o ponto sistólico (contração do coração), podemos usar esse fato para desenvolver um algoritmo de localização de pico que procura cada pico do gradiente sistólico.

O código de gradiente em Python é mostrado abaixo. Ele faz o seguinte:

  1. Registre 4 segundos de dados MAX30102
  2. Suavize os dados com uma curta convolução
  3. Calcule o gradiente com a função ‘gradiente()’ do Numpy
  4. Procure picos onde o gradiente sistólico é máximo
  5. Frequência cardíaca aproximada do período entre os picos
import serial,time,os
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')

ser = serial.Serial('/dev/ttyUSB0',
                    baudrate=115200)
    
start_word = False

heart_rate_span = [10,250] # max span of heart rate
pts = 1800 # points used for peak finding (400 Hz, I recommend at least 4s (1600 pts)
smoothing_size = 20 # convolution smoothing size

# setup live plotting
plt.ion()
fig = plt.figure(figsize=(14,8))
ax1 = fig.add_subplot(111)
line1, = ax1.plot(np.arange(0,pts),np.zeros((pts,)),linewidth=4,label='Smoothed Data')
line2, = ax1.plot(0,0,label='Gradient Peaks',marker='o',linestyle='',color='k',markersize=10)
ax1.set_xlabel('Time [s]',fontsize=16)
ax1.set_ylabel('Amplitude',fontsize=16)
ax1.legend(fontsize=16)
ax1.tick_params(axis='both',which='major',labelsize=16)
plt.show()

while True:
    t_vec,y_vals = [],[]
    ser.flushInput()  
    try:
        print('Place Finger on Sensor...')
        while len(y_vals)<pts:
            curr_line = ser.readline()
            if start_word == False:
                if curr_line[0:-2]==b'MAX30102':
                    start_word = True
                    print("Program Start")
                    continue
                else:
                    continue
            curr_data = (curr_line[0:-2]).decode("utf-8").split(',')
            if len(curr_data)==3:
                try:
                    t_vec.append(float(curr_data[0])/1000.0)
                    y_vals.append(float(curr_data[2]))
                except:
                    continue
        ser.flushInput() # flush serial port to avoid overflow
        
        # calculating heart rate
        t1 = time.time()
        samp_rate = 1/np.mean(np.diff(t_vec)) # average sample rate for determining peaks
        min_time_bw_samps = (60.0/heart_rate_span[1])
        
        # convolve, calculate gradient, and remove bad endpoints
        y_vals = np.convolve(y_vals,np.ones((smoothing_size,)),'same')/smoothing_size
        red_grad = np.gradient(y_vals,t_vec)
        red_grad[0:int(smoothing_size/2)+1] = np.zeros((int(smoothing_size/2)+1,))
        red_grad[-int(smoothing_size/2)-1:] = np.zeros((int(smoothing_size/2)+1,))  

        y_vals = np.append(np.repeat(y_vals[int(smoothing_size/2)],int(smoothing_size/2)),y_vals[int(smoothing_size/2):-int(smoothing_size/2)])
        y_vals = np.append(y_vals,np.repeat(y_vals[-int(smoothing_size/2)],int(smoothing_size/2)))

        # update plot with new Time and red/IR data
        line1.set_xdata(t_vec)
        line1.set_ydata(y_vals)
        
        ax1.set_xlim([np.min(t_vec),np.max(t_vec)])
        if line1.axes.get_ylim()[0]<0.95*np.min(y_vals) or\
           np.max(y_vals)>line1.axes.get_ylim()[1] or\
           np.min(y_vals)<line1.axes.get_ylim()[0]:
                        ax1.set_ylim([np.min(y_vals),np.max(y_vals)])
        plt.pause(0.001)

        # peak locator algorithm
        peak_locs = np.where(red_grad<-np.std(red_grad))
        if len(peak_locs[0])==0:
            continue
        
        prev_pk = peak_locs[0][0]
        true_peak_locs,pk_loc_span = [],[]
        for ii in peak_locs[0]:
            y_pk = y_vals[ii]
            if (t_vec[ii]-t_vec[prev_pk])<min_time_bw_samps:
                pk_loc_span.append(ii)
            else:
                if pk_loc_span==[]:
                    true_peak_locs.append(ii)
                else:
                    true_peak_locs.append(int(np.mean(pk_loc_span)))
                    pk_loc_span = []

            prev_pk = int(ii)

        t_peaks = [t_vec[kk] for kk in true_peak_locs]
        if t_peaks==[]:
            continue
        else:
            print('BPM: {0:2.1f}'.format(60.0/np.mean(np.diff(t_peaks))))
            ax1.set_title('{0:2.0f} BPM'.format(60.0/np.mean(np.diff(t_peaks))),fontsize=24)
        # plot gradient peaks on original plot to view how BPM is calculated
        scatter_x,scatter_y = [],[]
        for jj in true_peak_locs:
            scatter_x.append(t_vec[jj])
            scatter_y.append(y_vals[jj])
            
        line2.set_data(scatter_x,scatter_y)

        plt.pause(0.001)

        savefig = input("Save Figure? ")
        if savefig=='y':
            plt.savefig('gradient_plot.png',dpi=300,facecolor=[252/255,252/255,252/255])
    except KeyboardInterrupt:
        break

O código acima representa a seção 4s dos dados de frequência cardíaca e coloca pontos pretos sobre o pico do gradiente sistólico. O período entre cada ponto sucessivo é usado para calcular a freqüência cardíaca usando a taxa de amostragem. Um exemplo de gráfico de saída é mostrado abaixo com a aproximação de BPM.

Conclusão

Neste tutorial, o sensor de oxímetro de pulso MAX30102 foi apresentado junto com sua biblioteca compatível com o Arduino e funcionalidade básica. Os dados do MAX30102 foram então analisados usando Python para aproximar o comportamento cíclico da contração do coração e período de relaxamento para medir a frequência cardíaca. Em seguida, uma análise de gradiente mais precisa e estável dos dados de refletividade foi usada para encontrar picos no comportamento sistólico do coração para calcular a frequência periódica do coração. O objetivo deste tutorial era demonstrar o poder de um oxímetro de pulso barato para monitoramento de saúde, enquanto também explorava métodos avançados de análise de dados para aplicação no mundo real.