完全自作 arduinoでパルスオキシメーター 〜プログラム解説編〜

前回の記事では”完全自作パルスオキシメーター”のハードウェアについて説明してきました。

動画は↓

今回は、ハードウェアで計測した赤色LEDと赤外LEDの透過光から得られた動脈波形から心拍(HeatRate)と酸素飽和度(SPO2)を実際に算出する、パルスオキシメーターのキモであるプログラムを解説しようと思います。

ちなみに私はプロのプログラマーではなく、地方で働く臨床工学技士という職種の医療従事者です。

ですので、プログラムは可読性が悪いかもしれませんし、上手い方法ではなかもしれません。こうすればもっと良いよなどあれば教えてください。

自作パルスオキシメーターのおさらい

前回の記事ではハードウェアについて説明していますが、少しだけ制御について説明していますが、今回も少しだけプログラムの動き方について、つまり制御について触れておきます。

制御フロー

もう制御に関してはフローチャートのまんまです。

赤色光や赤外光LEDの光を指に照射しフォトダイオードで透過光を受光するとフォトダイオードで検知できる光強度をオーバーしたり、脈を検知するには低い透過度の場合があります。

指など簡単にうごく対象の場合、室内の光などの外来光の影響をフォトダイオードはバリクソ影響を受けます。

ですので、フォトダイオードにいい状態の透過光を常に安定的に照射する制御にしなければ安定した酸素飽和度は算出できません。

常時同じ光強度で照射するだけでは指を少し動かしたり、センサ位置を変えただけで簡単に光強度が強くなったり弱くなったりしてうまくいきません。

ですので、センサが読み取る光強度をMAX強度を100%として60%から90%の範囲で受光できるように常にLED駆動電圧を調節する制御をしています。

さらにですが、LEDの光強度を調節しても透過光が強かったり、弱かったりする場合はA/D変換コントローラーのゲイン=プログラムゲインアンプを変更させてより細かなパルスを読み取るようにしています。

ざっくりした制御の説明はこれで終わりです。

まぁフローチャートの方がより詳しく説明しているかな。

プログラムについて

プログラムの流れ

プログラムの全体像についてですが、まずはADS1115のA/DコントローラーやMCP4725のD/Aコントローラーを使う為のライブラリ、波形のDC成分を算出する移動平均ライブラリを使用する為の設定を行なう初期設定、いわゆるsetupをおこないます。

そして、arduino(ESP8266)で使用するIOpinの設定やI2C通信で使うPIN設定、赤色・赤外LEDを駆動するPIN設定などをおこないます。これもsetupですね。

上記が終われば測定モードに移行し、指が挟まれていないと自動的にスリープモードになり、指が挟まれるまで待ちます。

指が挟まればすぐに測定モードに移行しフローチャートの如く心拍と酸素飽和度を算出するという流れです。

ちなみにですが、脈波で心拍を一泊ごとに検知するのですが、この検知ができないと酸素飽和度が計算できません。

ですので脈波の波高値が小さいと結果酸素飽和度は算出できません。

このあたりがプログラム上での数値の加減で調節しましたが、どうでしょう。精度に影響ができますので、今後ブラッシュアップが必要かもしれまん。

実際のプログラム

なにわともあれ、実際にプログラムをみてください。

少しでもarduinoIDEやExcelのif分岐などをいじった事があればすぐに理解できると思います。

とりあえずarduinoのinoファイルをUPしておきます。ダウンロードしたいかたはどうぞ。

HandmadePulseoximeter.inoのダウンロード

ESP8266をお持ちならダウンロードファイルを書き込んで前回記事のハードウェアを作れば自作パルスオキシメーターの出来上がりです。

ここからはプログラムソースを分解して説明していきまうす。

まずは設定にかかわるソースを説明しまーす。というかソース中に説明書いてるのでその通りです。

//ADS1115でDRを設定できるライブラリ
//ESP8266を使用

//ライブラリの読み込み 事前にarduinoIDEでライブラリのインストールが必要です。
#include <Wire.h>//I2Cを使う
#include <ADS1X15.h>//16bitADC1115のライブラリ
#include <Adafruit_MCP4725.h>//LEDに出力する電圧を変化させるDAC
#include<movingAvg.h>//移動平均を取る為のライブラリ

//ローパスフィルタの定数 今回は使ってません
#define a 0.9


ADS1115 ADS(0x48);//ADS1115のインスタンスを作成()内はi2cのアドレス
Adafruit_MCP4725 dac; // MCP4725のインスタンスを作成 電圧を設定

//ESP8266のLED出力pinを決定 他のarduinoにする場合番号にかえる
int R_LED = 14;//ESP8266でいうとD5
int IR_LED =12;//ESP8266でいうとD6

//fotoダイオードで読み取る変数
float R_LEDSence;
float ROOMSence;
float IR_LEDSence;

//fotoダイオードで読み取る変数
float red_v;
float ir_v;

//MCP4725で電圧決める変数
//LEDに供給する電圧のデジタル値0-4095段階、LEDの光強度より2200-4095の範囲を使用
unsigned short int voltR=2200;
//unsigned short int voltIR=1680;//2300-700=1600 こっちは使ってない

//sleep変数  スリープ処理する時に使う
bool taiki=false;

//ads1115のゲインを調節する時に使用
int adsgain=0;



//時間計測の為の変数 プログラムする時のデバックに使った変数 使ってません。
unsigned long time1;
unsigned long time2;
unsigned long time3;
unsigned long time4;
unsigned long time5;

//ローパスフィルタの変数 使ってない
float rc =0;
float rc2 =0;

//------ 定数定義
const uint32_t TH_Terget = 30000;   //指が置いてあるかどうかの閾値30000
const int32_t  TH_AMOUNT = 100; //パルスの起点となる検出変化量

//測定対象がある事を感知する為の変数0ー100%表示40-90の間で測定するようにする
float TH_RED_Tergetrate;
float TH_IR_Tergetrate;


//pulse sop2 
double hr;
double spo2;


const int32_t MIN_INIT = 9999999; //最小値の初期値
const int32_t MAX_INIT = 0;       //最大値の初期値
 
//表示する心拍数とSpO2の範囲(用途により適宜変更)
const uint32_t DISP_MIN_HR   = 30;
const uint32_t DISP_MAX_HR   = 500;
const uint32_t DISP_MIN_SPO2 = 70;
const uint32_t DISP_MAX_SPO2 = 100;
 
//パルス検出保持変数
long l_time = millis();    //最後のパルス検出時間の保持
int32_t before_ir_v = 0;   //最後のir_v値保持
int32_t b_diff = 0;        // 最後の差の保持
long pulse_interval = -1;  // パルス間隔時間
 
//1パルス中の生データの最大・最小値の保持
int32_t min_ir_v = MIN_INIT ,max_ir_v  = MAX_INIT;
int32_t min_red_v= MIN_INIT ,max_red_v = MAX_INIT;
 
//移動平均値(IR_DC、RED_DC) 30(20~50)サンプル
movingAvg avgIr_v(30);
movingAvg avgRed_v(30);
 
//移動平均値(心拍数は直近3ビート、SPO2は5ビートの平均)
movingAvg avgHR(3);
movingAvg avgSPO2(5);

次のソースが実際のセンサデータ取得やメイン処理に関わるプログラムになります。

スポンサードリンク

もっと細分化すればよかったですが、ちょっとめんどくさくなったので一気に出します。

例の如くソース内で説明してます。

//近赤外光を照射した時のセンサ読み取り
void irread(){
  time3 = micros();
  digitalWrite(IR_LED,HIGH);
  //A/Dは2回の加算平均してます。意味があるかどうかわかりませんが。
  for(int i=0;i<2;i++){
    IR_LEDSence += ADS.readADC_Differential_0_1();
    }
  IR_LEDSence = IR_LEDSence/2;
  time4= micros()-time3;
  time5=micros()-time1;
  digitalWrite(IR_LED,LOW);
  }


//なにも光を当てていない時のセンサ読み取り //使ってません、この処理で複雑な指検知をしようと思ってました。
 void roomread(){
  for(int i=0;i<2;i++){
    ROOMSence += ADS.readADC_Differential_0_1();
    }
  ROOMSence = ROOMSence/2;
  }
  
//セットアップです。ここでシリアル通信の設定とi2cで使うpin設定をします。
void setup() 
{
 Serial.begin(115200);
 Wire.begin(4,5);//SDA,SCL ESP8266ではD1が5、D2が4です。

 //LED setup 出力に設定ね
 pinMode(IR_LED,OUTPUT);
 pinMode(R_LED,OUTPUT);
 
//初期は0vに。
 digitalWrite(IR_LED,LOW);//IR_LED
 digitalWrite(R_LED,LOW);//R_LED


//ここからADS1115の設定です
/* ----- Setting of ADS1115 AD converter ----- */
  ADS.begin();//スタートね
  ADS.setMode(1);//mode 0:continuous 1:single シングルモードを使用。
  ADS.setDataRate(4);//datarate860sps ads1115 time 1/860=1.16279msec  A/Dのスピード4は128で変換時間はだいたい8msec
  ADS.setGain(4);//PGA value +-1.024V  ゲイン設定 初期は4倍で他に8倍、16倍にできる。詳細はライブラリの説明を読んで

//MCP4725のi2cアドレスを設定
/* ----- Setting of MCP4725 DA converter ----- */
  dac.begin(0x60);


  //初期値として一回読んでおく(ir_vをパルス計測に用いる)
  irread();
  before_ir_v =  IR_LEDSence;
  
  //パルス最終取得時間の初期化
  l_time = millis();
 
  //移動平均ライブラリの初期化
  avgIr_v.begin();
  avgRed_v.begin();
  avgHR.begin();
  avgSPO2.begin();
}

//ローパスフィルタ 波高が下がるし波形が変形するので使ってない
void lowpassfilter()
  {   
      //ソフトウェアローパスフィルター
      rc = a * rc + (1-a) *red_v;
      //ソフトウェアローパスフィルター2
      rc2 = a * rc2 + (1-a) *ir_v;
  }


//このTransparentcalcで透過光の閾値をいい感じにして測定できる状態になったら測定へ
void Transparentcalc(){
  //センサ情報を取得 光強度は最低から2300*0.805664064=1.853v 1600*0.805664064=1.289v
  dac.setVoltage(voltR,false);//3.3V VDD max4095 
  redread();
  red_v = R_LEDSence;
 

  dac.setVoltage(voltR-700,false);//2.728V VDD max3410
  irread();
  ir_v = IR_LEDSence;

  //lowpassfilter();
  //red_v=rc;
  //ir_v=rc2;
  //taisyouhannteiha赤色LEDde判断
  TH_RED_Tergetrate=red_v/32768*100;
  TH_IR_Tergetrate=ir_v/32768*100;

  
  if(50 <= TH_RED_Tergetrate && TH_RED_Tergetrate <= 90){
    
  //透過50-90の間の場合
   //ここで心拍と酸素飽和度を計算にいく
    pulsespo2calc();

    }else if(TH_RED_Tergetrate<50){
      //透過50%以下 redledの電圧を決めるvoltRのMAXは4095 MINは2300
      if(voltR>=4095){
        //光強度MAXで透過が50以下なので測定不能
        //測定対象がごつい可能性
        //Serial.print("Cannot measure");
        voltR=4095;
        switch(adsgain)
          {
            //from gain4
            case 0:
            ADS.setGain(8);
            adsgain=1;
            break;
            
            //from gain8
            case 1:
            ADS.setGain(16);
            adsgain=2;
            break;
            
            //from gain16
            case 2:
            //Serial.print("Cannot measure1");
            ADS.setGain(4);
            adsgain=0;
            break;  
          };
        }else{
          //voltRをプラスして電圧を上昇させる
          //voltR=voltR+5;
            //voltR=4095;
            //voltR=(voltR*90)/TH_RED_Tergetrate;
            voltR=voltR*2*90/100;
          }
      
      }else if(voltR<=2100){
        
        //透過が90%以上の場合 redledの電圧を決めるvoltRのMAXは4095 MINは2150
        //光強度がMINなのに透過が90%以上の状態
        //測定対象が薄っぺらい若しくは測定対象が挟まっていない状態
        voltR=2100;
        switch(adsgain)
          {
            //from gain16
            case 2:
            ADS.setGain(8);
            adsgain=1;
            break;
            
            //from gain8
            case 1:
            ADS.setGain(4);
            adsgain=0;
            break;
            
            //from gain4
            case 0:
            //Serial.print("Cannot measure2");
            ADS.setGain(4);
            adsgain=0;
            taiki=true;
            sleep();
            break;  
          };
        
        }else{
          //voltR=voltR-5;
          //voltR=(voltR*90)/TH_RED_Tergetrate;
          voltR=voltR*90/100;
          }
  }


//心拍と酸素飽和度を算出する関数
void pulsespo2calc(){
/*この処理はテストで動かしてただけなのでコメントアウトしてます。使ってません。
 //----------Comment out except during testing-------
  //RLEDSence
  dac.setVoltage(voltR,false);//3.3V VDD 4095
  redread();
  red_v = R_LEDSence;
  TH_RED_Tergetrate=red_v/32768*100;
  
  //IRLEDSence
  dac.setVoltage(voltR-907,false);//2.728V VDD 3410
  irread();
  ir_v = IR_LEDSence;
  TH_IR_Tergetrate=ir_v/32768*100;
//--------------------------------------------------
*/
  

  //移動平均値(IR_DC RED_DCの算出) 移動平均がたったこの2行で終わるなんて。。
 //計算で使うDC成分ね。
  double ir_v_dc = avgIr_v.reading(ir_v);
  double red_v_dc = avgRed_v.reading(red_v);
  
  //IRとREDのACを求める為に最大値・最初値の更新 こっちが計算で使うAC成分になりまーす
  if(ir_v<min_ir_v) min_ir_v = ir_v; if(ir_v>max_ir_v) max_ir_v = ir_v;
  if(red_v<min_red_v) min_red_v = red_v; if(red_v>max_red_v) max_red_v = red_v;
 
  //パルス検出にはir_vを利用。前回との差分(変化値)を求める
  int32_t diff =  before_ir_v - ir_v;

//シリアルプロットで波形を数値を見るために記述 別に心拍と酸素飽和度だけでもOK
      Serial.print(TH_IR_Tergetrate);
      Serial.print(",");
      Serial.print(TH_RED_Tergetrate);
      Serial.print(",");
      Serial.print("adsgain=");
      Serial.print(adsgain);
      Serial.print(",");
      Serial.print("voltR=");
      Serial.print(voltR);
      Serial.print(",");
      Serial.print("HR=");
      Serial.print(hr);
      Serial.print(",");
      Serial.print("SPO2=");
      Serial.println(spo2);
      
//ここが心拍を検出するプログラム、ここで心拍が検出されないと酸素飽和度も計算されない
//もうすこし心拍の検出する精度をあげれたらなぁと思ってます。
  //TH_AMOUNT以上の変化がある場合、パルス起点とする
  if(b_diff < TH_AMOUNT && diff > TH_AMOUNT){
    
    //1パルスの時間を計算 心拍数の計算に使う
    pulse_interval = millis() - l_time;
    l_time = millis();


    //1パルスの時間より心拍数の計算。avgHRは整数のみなので1000倍して1000で割り少数保持
     hr =  (double)avgHR.reading(60000*1000/pulse_interval) / 1000.0;
 
    //SPO2の計算
    //IR・REDのACを求める(振幅最大-振幅最小)
    int32_t ir_v_ac = max_ir_v-min_ir_v;
    int32_t red_v_ac = max_red_v-min_red_v;
 
    // R = (AC_RED / DC_RED) / (AC_IR / DC_IR)の計算式より
    double red_div = double(red_v_ac)/red_v_dc;
    double ir_div = double(ir_v_ac)/ir_v_dc;
    double R = red_div / ir_div;
    
    // SPO2 = -45.060*R^2 + 30.354*R + 94.845 これはspo2_algorithm.cppにあったのを準用
    // 乗除1000は少数保持のため
     spo2 = (double)avgSPO2.reading((-45.060*R*R + 30.354*R + 94.845)*1000.0) / 1000.0;
 
    //最大値・最小値の初期化
    min_ir_v = MIN_INIT;
    max_ir_v = MAX_INIT;
    min_red_v = MIN_INIT;
    max_red_v = MAX_INIT;

    //心拍数とSPO2の表示(範囲は定義で変更可)
    //if(hr <= DISP_MAX_HR && hr >= DISP_MIN_HR && spo2 <= DISP_MAX_SPO2 && spo2 >= DISP_MIN_SPO2){    
      
      //Serial.print(TH_RED_Tergetrate);
      //Serial.print(",");
      //Serial.print(TH_IR_Tergetrate);
      //Serial.print(",");
      //Serial.print("HR=");
      //Serial.print(hr);
      //Serial.print(",");
      //Serial.print("SPO2=");
      //Serial.println(spo2);
    //}

  }


  //パルス検出用値の保持
  before_ir_v = ir_v;
  b_diff = diff;
  
  }

//ここからスリープの関数
void sleep(){
    while(taiki){
  //RLEDSence
  dac.setVoltage(2100,false);//
  redread();
  red_v = R_LEDSence;

  delay(10);
  
  //IRLEDSence
  dac.setVoltage(1600,false);//
  irread();
  ir_v = IR_LEDSence;

  delay(1000);
  //Serial.print(ir_v);
  //Serial.print(",");
  //Serial.println(red_v);
  
  //指を置いているか15bit32767であれば30000は21%程度
  if(red_v< TH_Terget || ir_v < TH_Terget){
    taiki=false;
    R_LEDSence=0;
    IR_LEDSence=0;
    ROOMSence=0;
    ir_v=0;
    red_v=0;
    
    }
 R_LEDSence=0;
 IR_LEDSence=0;
 ROOMSence=0;
 ir_v=0;
 red_v=0;
      
      
      } 
  }


void loop() 
{
  
 Transparentcalc();

 //test
 //pulsespo2calc();

 R_LEDSence=0;
 IR_LEDSence=0;
 ROOMSence=0;
 ir_v=0;
 red_v=0;
  
}



// -- END OF FILE --

プログラム的には順を追って読んでいけば難しくないです。

それではまた!!

スポンサードリンク
カテゴリー: 制作事例 パーマリンク

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください