MAX3010x Pulse Oximeter Modules (Part 4)

The MAX3010x sensors can do more than just heart rate measurements. Let's do some pulse oximetry.

Pulse Oximetry

In the last parts of this MAX3010x series we focused on heart rate measurements. However, for this task a much simpler sensor would be sufficient. The specialty of the MAX3010x sensors is that they can be used for pulse oximetry.

Pulse oximetry is a method to determine the peripheral oxygen saturation also known as SpO2 value. The peripheral oxygen saturation is normally similar to the arterial blood oxygen saturation within an accuracy of about 2%. The latter is what really matters. However, the huge advantage of using pulse oximetry is that it is a non-invasive method.

How does it work? All the MAX3010x sensors have at least two LEDs with different wavelengths. A red LED (660 nm) and an IR LED (880 nm). The red blood cells use the protein hemoglobin for transporting oxygen. Hemoglobin that is carrying oxygen (oxygenated hemoglobin) has a different color than deoxygenated hemoglobin. In consequence, they absorb light of specific wavelengths differently. By using two LEDs of different wavelengths and measuring the absorbance, we can in theory determine the percentage of oxygenated hemoglobin using Lambert-Beer-Law. This percentage is what is the SpO2 value. Normal values at sea level are between 95 % and 99 % at rest. Below 90 % you suffer from hypoxemia. Values below 80 % are critical for humans.

WARNING: In this tutorial a non-calibrated sensor is used. The measured values can deviate from the real ones. The measurements should by no means be used to indicate any medical conditions.

SpO2 Calculation

How can we calculate the SpO2 value? Well, there is no real formula for that. There are to many unknown factors. We don't know how thick the finger is, what skin tone it has, heck, we don't even known the intensity of the LEDs. Luckily, Takuo Aoyagi discovered the so called R-value that correlates with the SpO2 value. He is known as the inventor of the pulse oximeter.

The R-value is defined as following: \(R = {{AC_{red} \over DC_{red}} \over {AC_{ir} \over DC_{ir}}}\)

\(AC_{red}\) and \(AC_{ir}\) are the pulsatile components of the measured signals while \(DC_{red}\) and \(DC_{ir}\) describe the DC signal components. The AC component is caused by the changes of the blood volume during the cardiac cycle. The DC component is determined by the LED intensity, the surrounding tissue and many more factors. By creating a ratio between them, our unknown factors get eliminated to a certain degree. If we calculate the ratio between the ratios of both LEDs, we get a value that correlates to the oxygen saturation of the blood. Because deoxygenated hemoglobin absorbs red light stronger than IR light, the R-value increases the lower the oxygen saturation is.

The R-value is dimensionless and meaningless on its own. For getting the SpO2 value, the sensor needs to be calibrated. The calibration requires clinical tests with multiple probands. However, Maxim Integrated, the manufacturer of the MAX3010x sensors, has published example calibration values for the MAX30101 sensor. They might not be accurate to all the sensors from the MAX3010x range, but they are good enough for a test with the MAX30105. We want to try out the technique, we don't need accurate results.

Maxim Integrated used a quadratic regression for calibration. This leaves us with the following formula for calculating the SpO2 value:
\(SpO_2 = a \cdot R^2 + b \cdot R + c\)

The calibration factors for the MAX30101 are:
\(a = 1.5958422\)
\(b = -34.6596622\)
\(c = 112.6898759\)

The Code

Now, that we know how to calculate the SpO2 value we need to adjust the code to gather the required information. We need to add the code for determining the AC and DC components, before the high pass filter, as this filter removes the DC component. The DC value is simply the average value, while the AC value is the amplitude of the AC component in our signal. To calculate the amplitude we need to determine the minimum and maximum value in each cardiac cycle. We can the use these three values - average, minimum and maximum - to calculate the R-value each time we detect a heart beat. From the R-value we can then compute the SpO2 value.

The code is again quite complex, and I don't want to get into every detail. The main change is the new MaxMinAvgStatistics class, which is used to gather the minimum, maximum and average value. Averaging can be enabled by setting kEnableAveraging to true. In this tutorial it is set to false, so that we can always see the most recent value. With that being said, here is the code:

Main program
filters.h

#include <MAX3010x.h>
#include "filters.h"

// Sensor (adjust to your sensor type)
MAX30105 sensor;
const auto kSamplingRate = sensor.SAMPLING_RATE_400SPS;
const float kSamplingFrequency = 400.0;

// Finger Detection Threshold and Cooldown
const unsigned long kFingerThreshold = 10000;
const unsigned int kFingerCooldownMs = 500;

// Edge Detection Threshold (decrease for MAX30100)
const float kEdgeThreshold = -2000.0;

// Filters
const float kLowPassCutoff = 5.0;
const float kHighPassCutoff = 0.5;

// Averaging
const bool kEnableAveraging = false;
const int kAveragingSamples = 5;
const int kSampleThreshold = 5;

void setup() {
  Serial.begin(9600);

  if(sensor.begin() && sensor.setSamplingRate(kSamplingRate)) { 
    Serial.println("Sensor initialized");
  }
  else {
    Serial.println("Sensor not found");  
    while(1);
  }
}

// Filter Instances
LowPassFilter low_pass_filter_red(kLowPassCutoff, kSamplingFrequency);
LowPassFilter low_pass_filter_ir(kLowPassCutoff, kSamplingFrequency);
HighPassFilter high_pass_filter(kHighPassCutoff, kSamplingFrequency);
Differentiator differentiator(kSamplingFrequency);
MovingAverageFilter<kAveragingSamples> averager_bpm;
MovingAverageFilter<kAveragingSamples> averager_r;
MovingAverageFilter<kAveragingSamples> averager_spo2;

// Statistic for pulse oximetry
MinMaxAvgStatistic stat_red;
MinMaxAvgStatistic stat_ir;

// R value to SpO2 calibration factors
// See https://www.maximintegrated.com/en/design/technical-documents/app-notes/6/6845.html
float kSpO2_A = 1.5958422;
float kSpO2_B = -34.6596622;
float kSpO2_C = 112.6898759;

// Timestamp of the last heartbeat
long last_heartbeat = 0;

// Timestamp for finger detection
long finger_timestamp = 0;
bool finger_detected = false;

// Last diff to detect zero crossing
float last_diff = NAN;
bool crossed = false;
long crossed_time = 0;

void loop() {
  auto sample = sensor.readSample(1000);
  float current_value_red = sample.red;
  float current_value_ir = sample.ir;
  
  // Detect Finger using raw sensor value
  if(sample.red > kFingerThreshold) {
    if(millis() - finger_timestamp > kFingerCooldownMs) {
      finger_detected = true;
    }
  }
  else {
    // Reset values if the finger is removed
    differentiator.reset();
    averager_bpm.reset();
    averager_r.reset();
    averager_spo2.reset();
    low_pass_filter_red.reset();
    low_pass_filter_ir.reset();
    high_pass_filter.reset();
    stat_red.reset();
    stat_ir.reset();
    
    finger_detected = false;
    finger_timestamp = millis();
  }

  if(finger_detected) {
    current_value_red = low_pass_filter_red.process(current_value_red);
    current_value_ir = low_pass_filter_ir.process(current_value_ir);

    // Statistics for pulse oximetry
    stat_red.process(current_value_red);
    stat_ir.process(current_value_ir);

    // Heart beat detection using value for red LED
    float current_value = high_pass_filter.process(current_value_red);
    float current_diff = differentiator.process(current_value);

    // Valid values?
    if(!isnan(current_diff) && !isnan(last_diff)) {
      
      // Detect Heartbeat - Zero-Crossing
      if(last_diff > 0 && current_diff < 0) {
        crossed = true;
        crossed_time = millis();
      }
      
      if(current_diff > 0) {
        crossed = false;
      }
  
      // Detect Heartbeat - Falling Edge Threshold
      if(crossed && current_diff < kEdgeThreshold) {
        if(last_heartbeat != 0 && crossed_time - last_heartbeat > 300) {
          // Show Results
          int bpm = 60000/(crossed_time - last_heartbeat);
          float rred = (stat_red.maximum()-stat_red.minimum())/stat_red.average();
          float rir = (stat_ir.maximum()-stat_ir.minimum())/stat_ir.average();
          float r = rred/rir;
          float spo2 = kSpO2_A * r * r + kSpO2_B * r + kSpO2_C;
          
          if(bpm > 50 && bpm < 250) {
            // Average?
            if(kEnableAveraging) {
              int average_bpm = averager_bpm.process(bpm);
              int average_r = averager_r.process(r);
              int average_spo2 = averager_spo2.process(spo2);
  
              // Show if enough samples have been collected
              if(averager_bpm.count() >= kSampleThreshold) {
                Serial.print("Time (ms): ");
                Serial.println(millis()); 
                Serial.print("Heart Rate (avg, bpm): ");
                Serial.println(average_bpm);
                Serial.print("R-Value (avg): ");
                Serial.println(average_r);  
                Serial.print("SpO2 (avg, %): ");
                Serial.println(average_spo2);  
              }
            }
            else {
              Serial.print("Time (ms): ");
              Serial.println(millis()); 
              Serial.print("Heart Rate (current, bpm): ");
              Serial.println(bpm);  
              Serial.print("R-Value (current): ");
              Serial.println(r);
              Serial.print("SpO2 (current, %): ");
              Serial.println(spo2);   
            }
          }

          // Reset statistic
          stat_red.reset();
          stat_ir.reset();
        }
  
        crossed = false;
        last_heartbeat = crossed_time;
      }
    }

    last_diff = current_diff;
  }
}

#ifndef FILTERS_H
#define FILTERS_H

/**
 * @brief Statistic block for min/nax/avg
 */
class MinMaxAvgStatistic {
  float min_;
  float max_;
  float sum_;
  int count_;
public:
  /**
   * @brief Initialize the Statistic block
   */
  MinMaxAvgStatistic() :
    min_(NAN),
    max_(NAN),
    sum_(0),
    count_(0){}

  /**
   * @brief Add value to the statistic
   */
  void process(float value) {  
    min_ = isnan(min_) ? value : min(min_, value);
    max_ = isnan(max_) ? value : max(max_, value);
    sum_ += value;
    count_++;
  }

  /**
   * @brief Resets the stored values
   */
  void reset() {
    min_ = NAN;
    max_ = NAN;
    sum_ = 0;
    count_ = 0;
  }

  /**
   * @brief Get Minimum
   * @return Minimum Value
   */
  float minimum() const {
    return min_;
  }

  /**
   * @brief Get Maximum
   * @return Maximum Value
   */
  float maximum() const {
    return max_;
  }

  /**
   * @brief Get Average
   * @return Average Value
   */
  float average() const {
    return sum_/count_;
  }
};

/**
 * @brief High Pass Filter 
 */
class HighPassFilter {
  const float kX;
  const float kA0;
  const float kA1;
  const float kB1;
  float last_filter_value_;
  float last_raw_value_;
public:
  /**
   * @brief Initialize the High Pass Filter
   * @param samples Number of samples until decay to 36.8 %
   * @remark Sample number is an RC time-constant equivalent
   */
  HighPassFilter(float samples) :
    kX(exp(-1/samples)),
    kA0((1+kX)/2),
    kA1(-kA0),
    kB1(kX),
    last_filter_value_(NAN),
    last_raw_value_(NAN){}

  /**
   * @brief Initialize the High Pass Filter
   * @param cutoff Cutoff frequency
   * @pram sampling_frequency Sampling frequency
   */
  HighPassFilter(float cutoff, float sampling_frequency) :
    HighPassFilter(sampling_frequency/(cutoff*2*PI)){}

  /**
   * @brief Applies the high pass filter
   */
  float process(float value) { 
    if(isnan(last_filter_value_) || isnan(last_raw_value_)) {
      last_filter_value_ = 0.0;
    }
    else {
      last_filter_value_ = 
        kA0 * value 
        + kA1 * last_raw_value_ 
        + kB1 * last_filter_value_;
    }
    
    last_raw_value_ = value;
    return last_filter_value_;
  }

  /**
   * @brief Resets the stored values
   */
  void reset() {
    last_raw_value_ = NAN;
    last_filter_value_ = NAN;
  }
};

/**
 * @brief Low Pass Filter 
 */
class LowPassFilter {
  const float kX;
  const float kA0;
  const float kB1;
  float last_value_;
public:
  /**
   * @brief Initialize the Low Pass Filter
   * @param samples Number of samples until decay to 36.8 %
   * @remark Sample number is an RC time-constant equivalent
   */
  LowPassFilter(float samples) :
    kX(exp(-1/samples)),
    kA0(1-kX),
    kB1(kX),
    last_value_(NAN){}

  /**
   * @brief Initialize the Low Pass Filter
   * @param cutoff Cutoff frequency
   * @pram sampling_frequency Sampling frequency
   */
  LowPassFilter(float cutoff, float sampling_frequency) :
    LowPassFilter(sampling_frequency/(cutoff*2*PI)){}

  /**
   * @brief Applies the low pass filter
   */
  float process(float value) {  
    if(isnan(last_value_)) {
      last_value_ = value;
    }
    else {  
      last_value_ = kA0 * value + kB1 * last_value_;
    }
    return last_value_;
  }

  /**
   * @brief Resets the stored values
   */
  void reset() {
    last_value_ = NAN;
  }
};

/**
 * @brief Differentiator
 */
class Differentiator {
  const float kSamplingFrequency;
  float last_value_;
public:
  /**
   * @brief Initializes the differentiator
   */
  Differentiator(float sampling_frequency) :
    kSamplingFrequency(sampling_frequency),
    last_value_(NAN){}

  /**
   * @brief Applies the differentiator
   */
  float process(float value) {  
      float diff = (value-last_value_)*kSamplingFrequency;
      last_value_ = value;
      return diff;
  }

  /**
   * @brief Resets the stored values
   */
  void reset() {
    last_value_ = NAN;
  }
};

/**
 * @brief MovingAverageFilter
 * @tparam buffer_size Number of samples to average over
 */
template<int kBufferSize> class MovingAverageFilter {
  int index_;
  int count_;
  float values_[kBufferSize];
public:
  /**
   * @brief Initalize moving average filter
   */
  MovingAverageFilter() :
    index_(0),
    count_(0){}

  /**
   * @brief Applies the moving average filter
   */
  float process(float value) {  
      // Add value
      values_[index_] = value;

      // Increase index and count
      index_ = (index_ + 1) % kBufferSize;
      if(count_ < kBufferSize) {
        count_++;  
      }

      // Calculate sum
      float sum = 0.0;
      for(int i = 0; i < count_; i++) {
          sum += values_[i];
      }

      // Calculate average
      return sum/count_;
  }

  /**
   * @brief Resets the stored values
   */
  void reset() {
    index_ = 0;
    count_ = 0;
  }

  /**
   * @brief Get number of samples
   * @return Number of stored samples
   */
  int count() const {
    return count_;
  }
};

#endif // FILTERS_H

The program consist of two files: the main program and the filters.h file. To create a second file click on the arrow on the right side of the Arduino IDE and click New Tab or just press Ctrl+Shift+N. Name the auxiliary file filters.h and paste the corresponding code from above. Opening a new tab in the Arduino IDE

The Result

The result is quite impressive. With my sensor I get an SpO2 reading of 96 % at rest. This result is probably a bit low, but still within the expected value range. I made the measurements at an altitude of 900 m and not at sea level. The values to expect decrease with rising altitude. This means the value could be okay. However, in the end, it that is impossible to say without a properly calibrated pulse oximeter for comparison.

Let's do a small experiment. What happens if I hold my breath? Does the oxygen saturation decrease? Here is the result:

Decreasing SpO2 values while breath-holding

Indeed, the oxygen saturation does decrease when holding your breath for a long time. However, this starts at the point, at which one already feel like running out of air. As soon as you start breathing again, the value rapidly increase back to the normal level.

WARNING: Don't try to achieve the lowest possible values. I don't want you to pass out and injure yourself during this experiment.

Previous Post Next Post