前回の記事では”完全自作パルスオキシメーター”のハードウェアについて説明してきました。
動画は↓何回もいいか・・（笑）
今回は、ハードウェアで計測した赤色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をお持ちならダウンロードファイルを書き込んで前回記事のハードウェアを作れば自作パルスオキシメーターの出来上がりです。
ここからはプログラムソースを分解して説明していきまうす。
分けてますがすべて繋がっていて一つのソースです。
まずは設定にかかわるソースを説明しまーす。というかソース中に説明書いてるのでその通りです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 |
//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); |
次のソースが実際のセンサデータ取得やメイン処理に関わるプログラムになります。
もっと細分化すればよかったですが、ちょっとめんどくさくなったので一気に出します。
例の如くソース内で説明してます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 |
//近赤外光を照射した時のセンサ読み取り 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(5,4);//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 -- |
プログラム的には順を追って読んでいけば難しくないです。ただ長いだけ・・
わからない事があればなんなりと聞いてみてください。
それではまた！！
スポンサードリンク
