ハートビート(heart beat)は、心臓の鼓動ですよね？タイトルだけだと何の事やら予想不能だと思いますが、今回作るのは脈拍センサ(と、脈拍計)です。
脈拍数(heart rate)は、健康な人の安静時では60〜90[bpm](beats per minute、1分間当たりの脈拍数)程度らしいです。キツイ運動をすると170bpm程度まで上がる人もいるみたいです。低い方で50bpmを下回ると、「徐脈性不整脈」と呼ばれる症状に該当する可能性があります。
この辺りを根拠に、今回作成する脈拍計では48〜180[bpm]の測定ができるように設計したいと思います。周波数遮断フィルタ(LPFとHPF)を使うので、ターゲット範囲外の脈拍数は計測できなくなります。周波数換算だと0.8〜3[Hz]ですね。
医療機器では、心臓から直接拍動を拾う物を「心拍計」、指先や耳たぶなど心臓から離れたところで計測する物を「脈拍計」と呼ぶみたいです。この定義に従うと、本ページの装置は「脈拍計」に相当します。
箇条書きでまとめてみますね。
主に血液で反射して戻ってくる光を測定するので、赤色光を使うのが基本です。が、理由は定かではありませんが、緑色光(約550nm)を使う機器もあるらしいです……。
電子工作の立場で考えてみると、発光器(エミッタ) → LED、反射光の受光器(ディテクタ) → フォトトランジスタ、なのですが、あいにく手元に可視光のフォトトランジスタが有りません。そこで代用品として、反射型赤外線センサ(TCRT5000L、下図)を使うことにします。
※TCRT5000Lが発光する赤外線の波長は約940nmです。
測定原理を以下図にまとめておきます。
実は、自然光にも赤外線は含まれています。自然光に含まれる赤外線が指を透過してディテクタに入射してしまうと、外乱要素(ノイズ)となってしまいます。よって可能であれば、赤外線を遮蔽できるような材質で影を作っておくのが理想的です。
今回は可視光フォトトランジスタが無かったので止む無く赤外線を使うのですが、別な理由で近赤外線を使う医療計器もある様です。
※残念ながら、可視光フォトトランジスタもないし、校正テーブルも持ち合わせていないので、本ページの装置では%SpO2は測定できません、あしからず。ただ、%SpO2の測定が可能と言うことは、赤外線でも脈拍の計測が可能だという事です。安心して先へ進みます。
本節では、測定システムを設計して行きます。
検出される信号のモデルを解説するために、フライングでIRエミッタ(IR LED)とIRディテクタ(IR photo-transistor)部の回路図をご覧に入れます(下図)。
IRディテクタQ1は、通常、赤外線を検知していなければOFFしているので、Vs点は抵抗RcでプルアップされてVddを指しています。Q1が赤外線を検知するとONして、コレクタ-エミッタ間に電流が流れます。その電流値(すなわち、反射赤外線の強度)に応じて、Vs点はGNDに引っ張られることになります。
Vs点の電圧値を以下イメージ図で説明します。
TCRT5000Lに指を当てるまでは、反射してくる赤外線が無いので、基本的に電源Vddにプルアップされています。上図左端の時刻に指を当てるとIRエミッタから発光された赤外線が指の中で反射して、IRディテクタに戻ってくるので、Vsのレベル(電位)が低下します。
この時、心臓の拍動による血流量の増減によって反射赤外線強度が変動します。これが、赤い矢印で示した信号成分であり、今回計測したい信号となります。
上図では、説明し易い様にかなり大きい変動として図示してありますが、実際には数10mV程度の非常に微小な信号だと思われます。普通にArduino UNOのアナログ入力(10ビットAD変換器)で読み込んでも、変動(すなわち脈拍)を検知することはできないでしょう。
つまり、信号を増幅して読み込むことが不可欠となります。
TCRT5000Lは赤外線リモコンの送信用IR LEDとして使うために購入したのですが、赤外線リモコンで使う場合には、如何にして遠くまで赤外線を飛ばすか？が、使い勝手を左右するため、サージ電流スペックを使ってガツンと電流を流して使いました(「#8-8 反射型赤外線センサ」参照)。
サージ電流スペックを使えたのは、赤外線リモコンがIRの発光に38kHz(デューティ1/3)のキャリアを使っていたためです。今回は脈拍の測定中ずっと赤外線を発光し続けるため、無茶はできません。定格電流スペック内で光らせます。
※If＝60mA(ABS.MAX)が絶対最大定格なのですが、ギリギリまで電流を流すIf＝50mAで実験してみたところ、逆にうまく計測できませんでした。後ほど解説(失敗事例)します。
今回のアプリケーションでは非常に微弱な信号を増幅して読み込む必要があるため、IRディテクタの後段にオペアンプによる増幅回路を接続することが必須です。オペアンプには、正常に動作可能な入力電圧範囲(Input Common Mode Voltage Range)が規定されています。
今回使用するLM358(以下図)では0〜Vcc−1.5Vです。
前節の回路図で説明した様に、赤外線を検知していない間、IRディテクタの出力Vsは電源電圧Vddにプルアップされているので、Vdd＝Vccで使ってしまうと正常動作範囲外の電圧が印加されてしまいます。
※これでLM358が「壊れる」様なことはありませんが、信号の増幅ができなくなるなど、正常に機能しなくなってしまう恐れがあります。
そこで、IRディテクタにはVdd＝3.3VをUNOから供給して使います。同じくVcc＝5VをUNOから供給するので、LM358の入力電圧範囲0V＜Vdd(＝3.3V)＜Vcc−1.5V(＝3.5V)に収まって、良い感じです。
オペアンプを使った増幅回路と言うと、「#6-2 オペアンプにspice」などで解説していますので、復習が必要な様でしたら、ご参照ください(以下、反転増幅回路例を「#6-2」から再掲)。
単電源で反転増幅回路を使う時、入出力電圧の関係式は以下の通りです。
Vo＝Vb−Rf/Rs×(Vi−Vb)
バイアス電圧Vbを0点、それより高い電圧をプラス、低い電圧をマイナスの値と見て、増幅(−Rf/Rs倍)します。増幅率のマイナス符号は、バイアス電位Vbに対して符号が反転して出力されることを示しています。ここで重要なことは、微弱な信号を増幅したい時、その信号の変動範囲の中央付近にうまくバイアスを設定できないと、微弱な信号だけを増幅することができない、と言う事です。
一方、TCRT5000LのIRディテクタで検出される反射赤外線強度は大きくばらつきます。以下図(再掲)で説明します。
上図には電圧波形が3本(水色、橙色、灰色)描いてありますが、例えば、人体の個体差(皮膚が厚いとか、指先の血流量が少ないとか)だけでなく、IRエミッタやIRディテクタの電気的特性に起因する個体差などにも影響して、3本の電圧波形の様にばらつく、と言う事です。
上のイメージ図の場合、水色→橙色→灰色の順に反射赤外線強度が強い想定で描画してあります。必要なのは赤い矢印で示したような、脈拍による微弱な信号変化成分なのですが、条件によってレベルがばらつくため、バイアス点Vbを決め打ちにできないことが最大の課題です。
アイディアとして、観測した信号電圧の低周波数成分をバイアスとして設定する方法を思いつきました。前節の図で、3本(水色、橙色、灰色)の電圧波形の低周波数成分を破線で作図しておきました。これがちょうど、取り出したい信号成分の真ん中あたりを貫いています。つまりこの破線を反転増幅回路のバイアス電圧として設定すれば課題をクリアできる、と、考えました。
低周波数成分を残す → LPF(Low Pass Filter、低域通過フィルタ)ですよね。最初に脈拍の計測対象周波数を0.8〜3[Hz]にすることを決めたので、0.8Hz以下の周波数だけを残して、それより高い周波数の信号は遮断してしまいます。言い換えると0.8Hz以下の低周波数成分をバイアス電圧Vbとするのです。
遮断周波数fc＝1/(2π×R×C)＝0.8Hzですから、この数式を満たすことができる抵抗RとコンデンサCを使ってパッシブLPFを形成します。回路としては以下図の様になります。
※いつもの様に、手持ちにバリエーションの少ないコンデンサC＝10μFから決めて、それに合わせてR＝20kΩを使います。この組み合わせでfc≒0.8Hzです。
入力信号成分が微弱なことと、脈拍による山と谷がわかれば良いことから、増幅率は×1,000倍に設定してあります。これでも、1mVの信号しか出力されなかった場合には1mV×1,000倍＝1Vにしかなりません……。実験してみて増幅率が足りなかったら、抵抗Rf，Rsの比で調整するか、あるいは単なる「比較回路」にしてしまうか？必要に応じて対策しようと思います。
LTspiceを使ったシミュレーションで動作を確認してみました(下図)。
※LTspiceについては、「逆引き索引」から過去記事をあたってみてください。
Vs(緑線)がIRディテクタからの出力信号です。脈拍をエミュレートして約1s周期で脈動させ、30mVほどの振幅に10mV程度のノイズ成分を加えておきました。VsのLPF通過後の信号がVb(青線)で、これがバイアス電圧となります。赤い線Voが、Vsの信号成分だけを×1,000倍に増幅した出力です。
想定では、Vbがもっとフラットになるイメージ(前節のイメージ図破線)を持っていたのですが、fc≒0.8Hzで遮断したところ思いの他脈動しています(上図青線)。ただ、出力VoにはVsの脈動に対応した矩形波らしき信号が出てきているので、まあ、想定の範囲内です。脈拍のセンシングは可能でしょう。
この回路の動作をもう少し考察しておきます。
と、言う事で、出力Voには高周波数のノイズ成分まで増幅されて出力されているため、立ち上がりと立ち下がりにノイズが発生していて、綺麗な矩形波になっていないのです。今度こそ本物のLPFを通して高周波数成分を除去しないと、脈拍数のカウントどころではありません。
オペアンプのフィードバック抵抗Rfに並列に、適切なコンデンサCfを接続すると一次アクティブLPFとなるので、これを使うのも一手です。この場合、遮断周波数fc＝1/(2π×Rf×Cf)であり、増幅率A＝−Rf/Rsと独立して設定できない(Rfが共通)ので、注意が必要です。
今回は、ハードウェアのアナログフィルタではなくて、デジタルフィルタ(ソフトウェア)でLPFを構成してみたいと思います。
ここで改めて、計測対象の脈拍の特徴についておさらいです。最初に脈拍の測定可能範囲を48〜180[bpm](＝0.8〜3[Hz])と決めました。前節の非反転増幅回路はfc≒0.8HzのHPFとなっていて、測定可能範囲の下限である0.8Hz以下の低周波数成分を遮断したことになります。
そう。ここで検討すべきLPFは、測定可能範囲の上限である3Hz以上の高周波数成分を遮断できれば良いことになります。
今回はソフトウェアでデジタルLPFを作ります。デジタルフィルタってちょっと苦手で、大学の授業も良く理解できなかったのだけど、比較的構成が簡単なFIR(Finite Impulse Response)フィルタを使おうと思います(下図)。
「Z-1」の箱は、遅延回路です。横につらつらとつながっているのは「シフトレジスタ」になっているということです。
タップ数(n＋1)に従ってa0〜anまで(n＋1)個の係数が必要となります。この系数にどんな値を設定するか？で、このフィルタの特性が変わってきます。(n＋1)の取り方や各係数anによって、LPF、HPF、BPF、BEFなど様々なデジタルフィルタを設計することができます。
次の課題は、これらのパラメータをどうやって決めるか？です。
デジタルフィルタ(離散値のフィルタ)なので、まずはサンプリング周波数から決定します。ナイキスト周波数(＝サンプリング周波数/2で定義される周波数)以下の信号がフィルタの対象となるので、今回の上限周波数3Hzがナイキスト周波数以下となるようにサンプリング周波数を3Hz×2＝6Hz以上に設定する必要があります。
てか、「たったの6Hzで良いの？」って、感じですね。少し余裕を持たせて、サンプリング周波数を10Hzに決定しました。サンプリング周期では100ms(＝hbt_msec_SMPL)ですので、非力なArduino UNOでも余裕をもって処理することができるでしょう。
タップ数(＝FIRフィルタの次数＋1)を決めないと先に進めないのですが、これが良くわからない……。
など、表面的な性質はわかったのですが、具体的なタップ数の決定プロセスについては解りませんでした。タップ数をNとすると、フィルタ結果を計算するために毎回N回の積和演算が必要となります。Arduino UNOでの処理なので、そんなにタップ数を増やすわけにも行かず、N＝15とすることにしました。← よってここは、ちと、自信無し……＜（１ｘ＿ｘ）＞。
最後に、各タップの係数15個を決定する必要があります。↑程度の理解力で、係数の決定なんてできるでしょうか？
できるんです！便利な世の中になったものですな〜。Mathematicaを使います。
Mathematica ： LeastSquaresFilterKernel[{"Lowpass", <ωc>}, <N>]
※Raspberry PiでMathematicaを動かした記事については「逆引き索引」から過去記事をあたってください。
これで、角周波数ωcを遮断周波数とするLPFの、タップ数N個の係数を求めることかできます。Mathematicaすごいな、ってか、頼り過ぎ（１＾＾；）っ？
なお、ここで指定するωcは「正規化遮断角周波数」、です。以下、ωcの考え方と計算方法を示します。
Mathematicaで係数を決定した時の様子を下図に示します。
大学時代「Z変換」とか意味が解らず、まじめに授業聞いていなかったけど、作り方を知っていてMathematicaみたいなツールが使えたお陰で、「僕にも弾けた(出典：懐かしすぎるCM、ヤマハ)」。
※ちゃんと授業聞いていれば「タップ数の決め方とか、ちゃんと解ったんじゃない？」ってな、突っ込みもありますな……（１＾＾；）。「これから」の方は、是非ちゃんと授業聞いて理解してくださいね。
全体の回路図は以下図の様になりました。
TCRT5000LのIRエミッタ(IR LED)はNPNトランジスタ(手持ちのS8050)でドライブしましたが、電流制限抵抗が2kΩなので流れる電流は2mAを下回ります。よって、直接ArduinoのI/Oポートでドライブしても問題ないレベルだと思われます。
※部品点数を削減したり、製作コストを下げる必要のある時に有効。
前の節で示したFIRフィルタをソフトウェアで実装します。見ての通りシフトレジスタと積和演算だけなので、実装は容易です。
const float fir_coef[fir_MAX]= // fir: coefficient
{
// Mathematica, M_PI*3Hz/(10Hz/2), 15
0.0267283,
-0.0504551,
2.3389e-17,
0.0756827,
-0.062366,
-0.0935489,
0.302731,
0.6,
0.302731,
-0.0935489,
-0.062366,
0.0756827,
2.3389e-17,
-0.0504551,
0.0267283
}
さっきMathematicaで求めた係数を配列に定義しておきます。サンプリング周期が100msもあって余裕なので、思い切ってfloat型で処理しちゃいました。
// fir shift register
fir_sreg[fir_ptr]=anlg;
fir_ptr++;
if (fir_ptr>=fir_MAX)
fir_ptr=0;
シフトレジスタはリングバッファで構成しました。
// fir filter
hbt_fir=0;
for (idx=0; idx<fir_MAX; idx++)
hbt_fir += int (fir_sreg[(idx + fir_ptr)%fir_MAX]*
fir_coef[(fir_MAX - 1) - idx] + 0.5);
毎回N＝15回の積和演算をしないと結果を得られないのが玉に瑕ですね。より良い周波数特性を得ようとタップ数を増加してしまうと、Arduino UNOでは処理しきれなくなってしまうでしょう。シフトレジスタ(リングバッファ)や係数配列もタップ数に合わせて増加するのでメモリも圧迫されます。
脈動する測定データから、山(peak)と谷(trough)を判定して、山の数をカウントする必要があります。山の数と、それに要した時間がわかれば、脈拍数を求めることができるからです。
IRディテクタの出力が微小なため、ピークが高かったり低かったりして一定の閾値で判定すると、カウントを失敗する恐れがあると思いました。そこで、以下の様なアルゴリズムで判定します。
上のアルゴリズムによって山を検出した時間を使って脈拍数を計算できますが、さすがに、1回、2回の測定結果で計算すると誤差が大きすぎるため、何回(hbt_MAX＝16)か分の測定値を平均して脈拍数とします。
具体的なインプリ方法に関しては、完成プログラムをご覧ください。
完成プログラム：heart_rate.ino
コンパイルしてUNOに転送すると、シリアルモニタから以下のコマンドを受け付けるようになっています。
ま、コマンド1個だけですけど……。他にやること(インプリすべき機能)が無かった。
【heart】コマンドを実行している時間だけ、IR LEDが光ります(見えませんけど)。指を当てずにIRディテクタの出力を監視すると、ノイズを拾ってでたらめな数値を出力してきますので注意です。脈拍を測定する時は、１）TCRT5000Lに指を当てる、２）シリアルモニタから【heart】コマンドを実行する、３）測定終了まで待つ、４）指を離す＆結果を確認する、と言う手順を守ってください。
※測定対象は「脈拍」なので、短時間で大きく変動することは、まず、有りません。一方、ノイズを検出している時は検出の周期が一定ではないはずなので、この特性を使うと、指を当てた時(すなわち測定したい時)と、そうではない時を切り分けることは可能だと思われます(未実施)。
以下、実験中の様子です。
実験サンプル：【 】内はシリアルモニタからの送信データ
【heart】
OK: heart 30[s]
HART: 0 21 784 0
HART: 100 -18 784 0
HART: 200 -18 784 0
HART: 300 41 784 0
HART: 400 -7 784 0
HART: 500 -79 784 0
HART: 600 158 784 0
HART: 701 628 784 0
HART: 801 865 784 0
HART: 901 793 784 0
HART: 1001 745 784 0
：
： (中略)
：
HART: 29003 809 784 63
HART: 29103 634 785 63
HART: 29203 219 784 63
HART: 29303 -106 784 63
HART: 29403 -91 3 64
HART: 29503 244 3 64
HART: 29603 641 3 64
HART: 29703 804 3 64
HART: 29803 838 692 64
HART: 29903 815 784 64
HART: 30003 629 784 64
30s間脈拍数を測定して、結果をシリアルモニタへ表示しました。先にTCRT5000Lへ指先を当てておいてからコマンドを実行してみてください。TCRT5000L全体を指先で覆い隠すように接触させてください(後ほどご覧いただく動画で確認できます)。検出した脈拍に応じてLED_BUILTINが点滅するようになっていますので、空いている方の手を胸に当ててみたりして実際の鼓動とLED_BUILTINの点滅を見比べてみてください。
※細かいことを言うと、FIRフィルタの出力遅延(位相差)があるので、タイミングは一致してないかもしれません。急激に脈拍が変化しない限り、鼓動と点滅の間隔は一致するはずです。
シリアルモニタに表示される4つの数値は、「計測した時間[ms]、FIRフィルタの出力[LSB]、瞬間のアナログ読み込み値[LSB]、脈拍数[bpm]」です。この数値をExcelにコピー＆ペーストして作成したグラフが以下図です。
上半分は、指を当ててから【heart】コマンドを実行した時の0 → 30[s]間の測定データです。その内、3〜12[秒]の間を拡大したのが下半分。横軸に時間の経過[ms]、縦軸にアナログデータ[LSB]を取っています。
オレンジ線が「瞬間のアナログ読み込み値」で、その値にデジタルLPF(FIRフィルタ)を適用したものが水色線です。綺麗に脈拍が計測できていますね。計測開始から5s辺りまでは計測データが安定しておらず波形が乱れていますが、それを過ぎると安定して計測できています。
※オレンジ線は0／1のデジタル信号(矩形波)に近いですが、LPF通過後の水色線の方は高周波数成分が遮断されて、なだらかな曲線になっています。オーバシュートとアンダシュート、リンギングが起きるのは、FIRフィルタの設計が今一つなせいかもしれません……。
上側グラフの赤線が、計測した脈拍数[bpm]で、スケールは右側の第二軸です。この時の脈拍数は64bpm(≒1.07Hz)でした。脈拍の「山」がhbt_MAX(＝16)回検出されるまで、脈拍数(第4数値)は0の出力となっています。その後、暫く計測を続けて20〜30[秒]辺りになると、実際の脈拍数に近い値が表示されるはずです。
2〜4[桁]の7SEG LEDやLCD1602等を使って表示するも良し、後は皆様にお任せします。
サンプルデータが良くないのか……、上の測定データだとオレンジ線には高周波数のノイズもあまり乗っていないし、「デジタルLPF要らないんじゃね？」って、感じです。しかも、この位はっきり検出されていれば適切な固定閾値で判定してしまっても計測できたかも？ですね。ちょっと、老婆心になって考え過ぎたかも？
多少、測定波形が崩れても正しく脈拍数を計測できるロバスト性を有している、と言う事ですから、有って困ることは無いでしょう（１＾−＾）ｖ。
測定中の様子を動画にしてみました。
ぱっと見では、ただのLチカと変わりませんが、ちゃんとセンシングした脈拍に合わせてLEDが点滅していて、シリアルモニタの方には測定した脈拍数[bpm]が表示されています。
動画では、「ただのLチカ？」という誤解を与えぬ様、オプションのLCDサブシステムに脈拍値を表示してみました。プログラム(heart_rate.ino)先頭付近の「//#define USE_LCD_SUBSYS」というコメントを外してコンパイルすると、LCDサブシステムに測定結果を表示します。「#10-3 LCDサブシステムを何とかせねば」等の記事を参考に、UNOとLCDサブシステムを接続することが必要です。
回路とプログラムが完成して最初に実験した時、うまく脈拍をセンスできませんでした。初めの10〜20[秒]くらいまではセンスできているのに時間が経つと全く検出できなくなってしまったのです。
原因は「赤外線発光強度が強すぎた」ためでした。
人体の中で反射してきた赤外線強度を測るので、「すごく明るく光らせた方が良いに違いない」と、スペックぎりぎりまで電流値Ifを上げておいたのだけど、これがダメでした。どうやら発光が明るすぎたことで反射光強度も強すぎたため、ハレーションみたいな状態になってしまったものと推察しています。
電流制限抵抗を増やして、今はIf＜2mA程度に絞ってあります。これで安定して計測できるようになりました。赤外線の発光強度は明るすぎても暗すぎてもダメなので、もし、うまく測定できないようならIfを振ってみてください。
装置を測定状態にして(【heart】コマンドを実行)、IRディテクタ(TCRT5000L)にTVの赤外線リモコン信号を送信したところ、脈拍検知を表示するLED_BUILTINがピカピカ光りました（１＾＾；）。そりゃそうだよ〜、赤外線を検知しちゃうんだから。脈拍を測定する時には、このような外乱が紛れ込まないように注意して下さい。
