MAX3010x Pulseoximetermodule (Teil 4)

Mit den MAX3010x-Sensoren sind nicht bloß Pulsmessungen möglich. Heute geht es um das Thema Pulsoximetrie.

Pulsoximetrie

In den letzten Teilen dieser MAX3010x-Serie haben wir uns auf Pulsmessungen beschränkt. Allerdings reicht für diese Aufgabe ein wesentlich einfacherer Sensor aus. Das Besondere an den MAX3010x-Sensoren ist, dass sie auch die Pulsoximetrie ermöglichen.

Pulsoximetrie ist eine Methode zur Bestimmung der peripheren Sauerstoffsättigung, auch bekannt als SpO2-Wert. Die periphere Sauerstoffsättigung stimmt normalerweise mit der arteriellen Blutsauerstoffsättigung mit einer Genauigkeit von etwa 2% überein. Letztere ist das, was wirklich zählt. Der große Vorteil der Pulsoximetrie ist jedoch, dass es sich um eine nicht-invasive Methode handelt.

Wie funktioniert Pulsoximetrie? Alle MAX3010x-Sensoren verfügen über mindestens zwei LEDs mit unterschiedlichen Wellenlängen. Eine rote LED (660 nm) und eine IR-LED (880 nm). Die roten Blutkörperchen verwenden das Protein Hämoglobin für den Transport von Sauerstoff. Sauerstoffhaltiges Hämoglobin, hat eine andere Farbe als sauerstoffarmes Hämoglobin. Folglich absorbieren sie Licht bestimmter Wellenlängen unterschiedlich. Wenn man zwei LEDs mit unterschiedlichen Wellenlängen verwendet und die Absorption misst, kann man theoretisch den Prozentsatz des sauerstoffhaltigen Hämoglobins mithilfe des Lambert-Beer-Gesetzes bestimmen. Dieser Prozentsatz ist der SpO2-Wert. Normale Werte liegen auf Meereshöhe zwischen 95 % und 99 % in Ruhe. Unter 90 % leidet man an Sauerstoffarmut (Hypoxämie). Werte unter 80 % sind für den Menschen kritisch.

WARNUNG: In diesem Tutorial wird ein nicht kalibrierter Sensor verwendet. Die gemessenen Werte können von den realen Werten abweichen. Die Messungen sollten auf keinen Fall als Indikation für den eigenen Gesundheitszustand verwendet werden.

Berechnung des SpO2-Wertes

Wie können wir den SpO2-Wert berechnen? Nun, dafür gibt es keine genaue Formel. Es gibt zu viele unbekannte Faktoren. Wir wissen nicht, wie dick der Finger ist, welchen Hautton er hat, ja, wir kennen nicht einmal die Helligkeit der LEDs. Glücklicherweise hat Takuo Aoyagi den sogenannten R-Wert entdeckt, der mit dem SpO2-Wert korreliert. Er ist als Erfinder des Pulsoximeters bekannt.

Der R-Wert ist wie folgt definiert: \(R = {{AC_{rot} \over DC_{rot}} \over {AC_{ir} \over DC_{ir}}}\)

\(AC_{rot}\) und \(AC_{ir}\) sind die pulsierenden Komponenten der gemessenen Signale, während \(DC_{rot}\) und \(DC_{ir}\) die DC-Signalkomponenten beschreiben. Die AC-Komponente wird durch die Änderungen des Blutvolumens während des Herzzyklus verursacht. Die DC-Komponente wird durch die LED-Intensität, das umgebende Gewebe und viele weitere Faktoren bestimmt. Durch die Bildung eines Verhältnisses zwischen ihnen werden die unbekannten Faktoren zu einem gewissen Grad eliminiert. Wenn wir den Quotienten zwischen den Verhältnissen der beiden LEDs berechnen, erhalten wir einen Wert, der mit der Sauerstoffsättigung des Blutes korreliert. Da sauerstoffarmes Hämoglobin rotes Licht stärker absorbiert als IR-Licht, steigt der R-Wert, je niedriger die Sauerstoffsättigung ist.

Der R-Wert ist dimensionslos und für sich genommen ohne Aussagekraft. Um den SpO2-Wert zu ermitteln, muss der Sensor kalibriert werden. Für die Kalibrierung sind klinische Tests mit mehreren Probanden erforderlich. Maxim Integrated, der Hersteller der MAX3010x-Sensoren, hat jedoch beispielhaft Kalibrierungswerte für den MAX30101 veröffentlicht, die wir für dieses Tutorial verwenden werden. Diese Werte sind möglicherweise nicht für alle Sensoren aus dem MAX3010x-Sensorsortiment passend, aber für einen Test sind sie allemal genug. Deshalb werde ich sie im folgenden auch für den hier genutzten MAX30105 verwenden. Wir wollen die Technik lediglich ausprobieren, wir brauchen keine exakten Ergebnisse.

Maxim Integrated benutzt eine quadratische Regression für die Kalibrierung. Daraus ergibt sich die folgende Formel zur Berechnung des SpO2-Wertes:
\(SpO_2 = a \cdot R^2 + b \cdot R + c\)

Die Kalibrierungsfaktoren für den MAX30101 sind:
\(a = 1,5958422\)
\(b = -34,6596622\)
\(c = 112,6898759\)

Der Code

Da wir nun wissen, wie der SpO2-Wert zu berechnen ist, müssen wir den Code anpassen, um die erforderlichen Werte zu ermitteln. Der Code zur Bestimmung der AC- und DC-Komponenten muss vor dem Hochpassfilter eingefügt werden, da dieser Filter die DC-Komponente entfernt. Der DC-Wert ist einfach der Durchschnittswert, während der AC-Wert die Amplitude der AC-Komponente in unserem Signal ist. Um die Amplitude zu berechnen, müssen wir den Minimal- und Maximalwert in jedem Herzzyklus bestimmen. Anhand dieser drei Werte - Durchschnitt, Minimum und Maximum - können wir bei jedem Herzschlag den R-Wert berechnen. Anhand des R-Werts können wir dann den SpO2-Wert bestimmen.

Der Code ist wieder recht komplex, und ich möchte hier nicht auf jedes Detail eingehen. Die wichtigste Änderung ist die neue Klasse MaxMinAvgStatistics, die verwendet wird, um den Minimal-, Maximal- und Durchschnittswert zu ermitteln. Es ist möglich die Endergebnisse zu mittlen, indem man das Flag kEnableAveraging auf true setzt. Für dieses Tutorial ist es auf auf false gesetzt, sodass wir immer den aktuellsten Wert sehen können. Das komplette Programm sieht wie folgt aus:

Hauptprogramm
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

Das Programm besteht aus zwei Dateien: dem Hauptprogramm und der Datei filters.h. Um eine zweite Datei zu erstellen, klicke auf den Pfeil rechts in der Arduino IDE und wähle Neues Tab oder drücke einfach Strg+Shift+N. Nenne die Hilfsdatei filters.h und füge den Code von oben ein. Öffnen eines neuen Tabs in der Arduino IDE

Das Ergebnis

Das Ergebnis ist ziemlich beeindruckend. Mit meinem Sensor erhalte ich einen SpO2-Wert von 96 % in Ruhe. Dieses Ergebnis ist wahrscheinlich etwas zu niedrig, aber immer noch im erwarteten Wertebereich. Ich habe die Messungen in einer Höhe von 900 m und nicht auf Meereshöhe durchgeführt. Die zu erwartenden Werte nehmen mit zunehmender Höhe ab. Es ist also möglich, dass dieser Wert stimmt. Letztlich kann man dies ohne ein korrekt geeichtes Pulsoximeter zum Vergleichen nicht sicher sagen.

Machen wir ein kleines Experiment. Was passiert, wenn man die Luft anhält? Nimmt die Sauerstoffsättigung dann ab? Hier ist das Ergebnis:

Abnehmende SpO2-Werte beim Luftanhalten

In der Tat sinkt die Sauerstoffsättigung, wenn man den Atem lange Zeit anhält. Dies beginnt aber erst in dem Moment, in dem man bereits das Gefühl hat, dass einem die Luft ausgeht. Sobald man wieder zu atmen beginnt, steigt der Wert rasch wieder auf das normale Niveau an.

WARNUNG: Probier bitte nicht, einen möglichst niedrigen Wert zu erreichen. Ich möchte nicht, dass du ohnmächtig wirst und dich während dieses Experiments verletzt.

Vorheriger Beitrag Nächster Beitrag