心拍データの収集と活用その1
AppleWatchで心電図データを収集する
AppleWatchで心電図 App で心電図 (ECG) をとることができるとの情報をもとに
測定方法などに興味があり調べてみることにしました。
AppleWatchで収集したデータをpythonの心電図ライブラリを使って解析した
報告が下記のURLで見つかりました。
pythonの心電図解析ライブラリ(BioSPPy)の利用例
光学式心拍センサの概要
光学式心拍センサは、日本で開発され現在特許が切れています。心拍センサの構造は単純です。
指先などに光源光を当て帰ってくる光を強度信号を取り出します。光源光は周囲の環境影響
が少ない緑色が使われます。血液中のヘモグロビンは光を吸収します。多くのヘモグロビンが
あると光の吸収が大きくなります。この吸収信号が心拍信号にほぼ同じということを利用しています。正式な心電図センサは電極を数か所体に張り付ける必要がありますが、光学式センサは不要です。一方で、安定した測定が難しい場合があります。
心拍センサ(APDS9008)
ArduinoNanoとの接続
AppleWatchより周波数分解性能を上げるため、
1KHzでセンサデータをサンプリングします。(iWatch=512Hzで約2倍)
ArduinoNanoのプログラムコード
心拍データの収集と保存プログラム
1msecで10秒サンプリングしセンサ信号を表示保存するプログラムを作りました。
この時点では、pythonの心電図解析ライブラリ(BioSPPy)は使用していません。
プログラム概要と処理の流れ
データ保存ファイル名を入力し⏎すると測定を開始
10秒間x1000=10000データをcsv形式で保存
測定波形(10秒分)をグラフ表示
DFT処理を行い心拍数に応じたピーク検出を行い心拍数を推定
心拍波形のグラフ表示
DFTによる心拍数推定
心拍データの収集と活用その2
心拍データサンプリングツールとプログラムの改善
プロト機での問題点は以下2点です。
・測定時間を延ばさないと周波数分解能が確保できないこと
・センサ部と親指の接触が不安定で、測定途中で信号振幅が大きく変化すること
今回は、センサ部の改善とプログラムの修正を行いました、
測定装置の改善(左→右)と測定時間の変更(10→30秒)
測定プログラム(ECGpre.ipynb)
データ解析プログラム(ECGpost.ipynb)
Python心電図ライブラリ(biosppy)を使用するとcsvデータからフィルター処理した波形グラフ、心拍数の時間変化グラフ、PQRSTU波形のTemplate波形グラフを一度にサマリーとして描きだすことができますが、その後のデータ処理や特徴値抽出にはひと手間が必要です。
Scipyを使って信号処理をひとつひとつ納得しながら処理する法を選びました。
Python心電図ライブラリ(biosppy)を使用した一括出力
心拍データの収集と活用その3
洞調律(一貫したパターンの拍動)の測定
心拍データからの有効な特徴値として、洞調律(RRI)や心拍波形を加重平均した加重平均波形
があります。今回は心臓の健常性制を示洞調律の抽出を行いました。
ScipyによるFilter処理について
心房細動があると高周波(120以上)の拍動が観察されます。
まず、R波の位置が測定時の不安定要因によるノイズを除く目的でバンドパスフィルタ
処理を行いました。 Scipyでは、引数filtered[]でRawデータから抽出されます。
また、filtered[]波形のピーク抽出はChristov法により検出されてINDEXに格納しました。
Rawデータ
BPF処理後のデータ
BPFにより高周波ノイズが低減され波形がそろってきているように思います。
洞調律(心拍数の時間変化)グラフ
BPF処理後、filtered[]波形のピーク抽出はChristov法により検出されてINDEX[]に
格納しています。INDEX[]の中身をグラフ表示したものを示します。
約0.8秒間隔で脈を打っていることが判ります。
DFTでは周波数軸でのピークを検出して平均心拍数を求めることができるが、1秒単位の時間変化のグラフは作成できないためピーク検出法を使用する合理性を確認することができました。更にピーク検出の制度を上げる目的でScipyのfind_peak()関数を検討しました。
Christov法によるピーク検出はピークの高さが整ってない場合は、ピークをカウントしません。 ピーク高が半分以下になっても検出できるfind_peak()関数を試しました。
SciPyのfind_peaks()関数を使用した検討
Scipy find_peak()関数によるピーク検出
ピーク高が半分以下になっても検出が可能出ることを確認しました。
RRIグラフ
HeartRate= 71 (/min) stdev(%)= 2
Python心電図ライブラリ(biosppy)でRRIグラフが出なかった原因
次のステップとして、交感神経の活性度指標を得るためにRRI波形のAC解析を行います。
そこでの問題は、RRI(R-R Interval)データの時間軸が等間隔でないことです。
FFT処理を行うため等間隔にする必要があり等時間間隔に補完しなければなりません。
心拍データの収集と活用その4
RRI波形からストレスと自律神経指標の抽出
自律神経のバランス指標として、心拍データから抽出した交感神経の活性化指標が使われます。
RRI波形のAC解析から、HF成分:副交感神経の活性指標とLF/HF比:交感神経の活動指標を
抽出します。LF値/HF値は以下のように定義されています。
[LF値:0.05-0.15Hzのパワースペクトル密度 HF値:0.15-0.45Hzのパワースペクトル密度}
RRI波形の等時間間隔補完(アッパーサンプリング)プログラム
200秒間の心拍波形からRRI波形抽出
①心拍波形の生データにピーク検出処理を実施
サンプリング:1KHz 測定時間:200秒(全データ数:200000)
②RRI波形の抽出 平均心拍数:67/分 バラツキ(標準偏差):8%
RRI波形の補間処理とFFT処理
③RRI波形のスプライン補完区分的3次エルミート補完を使用
FFT処理のため200秒間のRRI波形を0.1秒間隔で再サンプリングを行いました。
④再サンプルRRI波形のFFT処理
200→2000に再サンプリングしたRRI波形にFFT処理を実施しました。
周波数軸で左右対称ですが5Hzより右側は鏡像です。
光学式心拍センサの波形から交感神経の活動指標を取り出すことができました。
ここまでくれば、AppleWatchの機能を超えたものになったかなと思います。
もうひとつの応用は波形テンプレートです。長さも振幅も異なる心拍波形を
どのようにテンプレート化するのか、これができれば脈波形の機械学習による
分類ができるように思います。 波形データの重心平均について調べてみました。
心拍データの収集と活用その5
心拍波形の加重平均処理
心拍波形をいくつかのtemplate波形として抽出するライブラリがあります。
測定した人の心拍モードを自動で抽出する機能のようです。
波形を切りだし時間軸を統一して、数種類の波形テンプレートを出力する機能
のようです。
Scipyのecg.ecg()関数を使用して波形Template を抽出する
測定を代表する心拍波形抽出について
統一時間軸で括られた複数のtemplate波形の平均を求めました。
DTW(Dynamic Time Warping)による非同期波形の重心平均化検討
「tslearn.barycenters」ライブラリを使用するることで、複数の波形の重心平均を求める
ことができます。以下、チュートリアルにあるプログラムコードと重心平均結果を示します。
「tslearn.barycenters」ライブラリのsoftdtw_barycenter()関数
pickleを使って(バイナリファイルを経由して)2つの環境間で処理を行う
複数心拍波形の重心平均処理プログラム
30秒間測定した心電図波形“nanoc30.csv”のtemplate波形の重心平均
処理パラメータとしてγがあり微調整が可能ですが、波形ピークに角が発生します。
パラメータγとtolによる角消去効果の確認
心拍データの収集と活用その6
プロトタイプ機とGUIプログラムの製作
光学式心拍センサを一般的に使っていただけるように、2022年7月に
センサ部構造とマイコンを変更しました。また、python環境のないWindows-PC
でも使用できるように、GUIプログラムをpysimpleGUIライブラリで製作し直し、
Pythonのプログラムコードをpyinstallerを使ってコンパイルし実行形式(.exe)で
配布することができました。プロトタイプ測定装置とプログラム操作の概要をまとめ
ました。
Interigentsensor1との接続
InteligentSensor1のコマンドとGUIプログラム一覧表
Inteligentsensor1には信号処理用のコマンドプログラムを9つ書き込んでいます。
また、それぞれのコマンドごとに接続するセンサに合わせたGUIプログラムを準備しました。
今回はコマンドモード「e」のGUIプログラムを紹介します。
心拍測定プログラム mode_eGUI.exeの操作方法
mode_eGUI.exeのソースコード(1/2)
mode_eGUI.exeのソースコード(2/2)