💓

Apple Watchの心電図データをPythonで分析して遊ぶ

2021/02/07に公開約10,600字11件のコメント

概要

日本国内でも、Apple Watchの心電図測定機能が解放されました。測定したデータをCSVで書き出し、Pythonの生体情報ライブラリであるBioSPPyによって分析してみました。
心電図の波形を見たり、ピーク検出を使ってRRI・RRVといった心電指標を計算したりできるものの、LF/HFなどの周波数領域資料の算出は難しいようです。

サンプルコードをGoogle Colabで公開しています。ご自身がアップロードしたCSVデータで心電図の解析をお試しいただけます!

Apple Watchによる心電図(ECG)測定

背景

2020年1月27日、日本国内で販売されたApple Watchを対象に、心電図アプリによる心電図(electrocardiogram; ECG)測定を可能にするアップデートが行われました。iOS 14.4とwatchOS 7.3にアップデートすれば、すでに購入済みのApple Watch Series 4、5、6でも自分のECGを見られます。

https://support.apple.com/ja-jp/HT208955

そもそもApple WatchによるECG測定は、2018年発売のSeries 4から搭載され、米国などでは利用できました。しかし、日本での展開に必要な医療機器認証をすぐには取得できず、機能が長期間封印されていました。わざと海外でApple Watchを購入するなどの手段を取れた人を除けば、多くの人にとって、待ちに待ったECG測定解禁となります。

ECGデータを解析して遊びたい

Apple Watchで測定したECGは、心臓に異常を感じている方にとっては間違いなく重要なデータでしょう。しかし、特に異常を感じておらず、心電図を読めるだけの知識もない多数のユーザからすると「測定してなんとなく波形を眺めて終了」になりかねません。せっかく日常的にECGを測定できるセンサを初めて手にしているのですから、取得したデータであれこれ遊びたいものです。

ECG機能の国内リリース直後に、三重大学の奥村先生がECGデータをPythonに読み込んでプロットした様子をツイートされていました。約512 Hzで30秒間の測定データを、数値データとして得られるようです。本記事では、この数値データを使って、PythonライブラリBioSPPyによる簡単なECGデータ分析を行ってみます。

ECGの測定とデータの取り出し

ECG測定が初めてなら、Apple公式の案内を読みながら「心電図App」を設定しましょう。アプリをWatchで開き、Digital Crownと呼ばれるボタンに触れると測定が始まります。測定自体は30秒で簡単に終わるものの、きれいな波形にするためには工夫が必要です。

上の画像は腕をしっかり固定して測定した場合、下の画像は測定中に腕を動かした場合です。腕を動かしてしまうと、波形が上下にブレたり、細かいノイズが乗ってしまったりすることが分かります。

Apple WatchでのECG測定では、Watchの背面から体内を通ってDigital Crownに触れた指までを電気回路として扱うことで電圧を測定します。測定する電圧はとても小さいため、少しでも腕などが動くと波形が大きくブレてしまうのです。机の上に腕を置くなど、30秒間しっかりと上半身を固定できる状況で測定しましょう。

測定を終えたら、iPhoneの「ヘルスケア」アプリに移ってデータを取り出す作業をします。奥村先生のWebサイトでも解説されているとおり、「すべてのヘルスケアデータを書き出す」機能でZipファイルを生成します。
この機能はECGに限らずあらゆるヘルスケアデータを書き出すので、過去に蓄積したデータが多いほど時間がかかります。気長に待ちましょう。

Zipファイルができたら、iOSの共有メニューが開きます。Macをお使いならAirDrop、そうでない場合はiCkoud Drive・Google Drive・Dropboxなどのクラウドストレージに保存するのが楽です。

ECGデータは、electrocardiogramsフォルダ内にCSVファイルとして格納されています。これをお好みのフォルダに展開すれば、分析の準備は完了です。

プロットしてみる

※本節以降で使用したコードをGoogle Colabで公開しています

まずは、奥村先生のコードを参考にmatplotlibで値をプロットしてみます。CSVの先頭に記載されているように、サンプリング周波数は512 Hz、単位はµVです。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
plt.rcParams['font.family'] = 'IPAexGothic'

ecg_df = pd.read_csv('ecg_2021-02-06.csv', skiprows=13, header=None)

fig = plt.figure(figsize=(10,2))

plt.plot(ecg_df[0])
plt.xlabel('sample')
plt.xlim(0, 15360)
plt.xticks(np.arange(0, 15360, 512 * 4))
plt.ylabel('ECG [μV]')
plt.ylim(-500, 1000)
plt.grid()

plt.show()

しっかり30秒分取れていることが分かります。

# (省略)
plt.xlim(512 * 5, 512 * 10)
plt.xticks(np.arange(512 * 5, 512 * 10, 512))
# (省略)

plt.show()

少し拡大してみると、よく見る心電図っぽい波形が見えてきました。

ECGの波形には、P波~U波までの名称が付けられています。もっとも高く鋭いピークはR波と呼ばれ、心臓が血液を送り出すタイミングを示します。
R波から次のR波までの間隔を RR間隔(R-R Interval; RRI) といい、心拍数などの計算に用いられます。

BioSPPyを使ってみる

Pythonで生体情報を扱えるライブラリは複数ありますが、本記事ではBioSPPy[1]をご紹介します。信号処理やパターン認識などの技術が搭載されており、光電容積脈波(PPG)・心電(ECG)・皮膚電気活動(EDA)・脳波(EEG)・筋電(EMG)・呼吸をサポートしています。
ポルトガルの研究機関Instituto de Telecomunicaçõesに所属するPattern and Image Analysis (PIA) Groupによってメンテナンスされており、研究者にとっても安心して使える品質のライブラリです。

https://biosppy.readthedocs.io/en/stable/

ECGデータの読み取り

ライブラリをpipからインストールします。

$ pip install biosppy

ecgモジュールを使ってECG信号を処理します。showオプションをTrueにすると、データの概要がまとめてプロットされます。ecg.ecg()では、データの周波数フィルタリング・R波ピークの検出・心拍数の算出が行われています。

from biosppy.signals import ecg

ecg_data = ecg.ecg(signal=ecg_df[0], sampling_rate=512., show=True)

概要プロットは、左列が上から生の波形・フィルタリングされた波形と検出したR波・心拍数の変動、右列がECGの1波形を示しています。

周波数フィルタ

BioSPPyのフィルタリングによる効果を確かめるために、生波形とフィルタリングされた波形を見比べてみましょう。

(ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate) = ecg_data

fig = plt.figure(figsize=(10, 4))
plt.subplots_adjust(wspace=0.4, hspace=0.7)

for (i, data, title) in [(1, ecg_df[0], '生データ'), (2, filtered, 'フィルタリング')]:
    ax = fig.add_subplot(2, 1, i)
    ax.plot(ts, data)
    ax.set_title(title, fontsize=10)
    ax.set_xlim(0, 30)
    ax.set_ylim(-500, 1000)
    ax.set_xlabel('time [s]')
    ax.set_ylabel('ECG [μV]')
    ax.grid()

plt.show()

画像はなるべく腕を固定して測定したECG波形です。生波形では0付近の波形がブレているのに対して、フィルタリング後は波以外の箇所がしっかり0に揃っているのが分かります。ソースコードを見ると、3 Hzから45 Hzまでのバンドパスフィルタを掛けているようです。

周波数フィルタリングは、体動によって波形がブレてしまったECGデータに対して威力を発揮します。画像はわざと腕を動かして測定したデータで、生波形は大きく動揺しています。フィルタ処理後も細かいノイズは残っているものの、R波を見やすい波形に変換できました。

R波のピーク検出

BioSPPyのECGモジュールには、5種類のR波検出アルゴリズムが実装されています。それぞれの論文は、一部を除いて公式ドキュメントに引用されています
標準のecg.ecg()で呼び出されているのはhamilton_segmenter()です。

  • christov_segmenter() : (Christov, 2004)
  • engzee_segmenter() : (Engelse and Zeelenberg, 1979; Lourenco et al., 2012)
  • gamboa_segmenter() : (Gamboa, 2008)
  • hamilton_segmenter() : (Hamilton, 2002)
  • ssf_segmenter() : (Zong et al., 2003)

あえてピーク検出を難しくするため、わざと腕を動かした時のECGデータに対して、それぞれのアルゴリズムを適用してみます。

# 体動入りECGデータ
(ts, filtered, _, _, _, _, _,) = ecg.ecg(signal=ecg_df_2[0], sampling_rate=512., show=False)

rpeaks_christov, = ecg.christov_segmenter(filtered, 512.)
rpeaks_engzee, = ecg.engzee_segmenter(filtered, 512.)
rpeaks_gamboa, = ecg.gamboa_segmenter(filtered, 512.)
rpeaks_hamilton, = ecg.hamilton_segmenter(filtered, 512.)
rpeaks_ssf, = ecg.ssf_segmenter(filtered, 512.)
rpeaks_ssf_2000, = ecg.ssf_segmenter(filtered, 512., threshold=2000.)

検出されたピークを散布図にプロットしてみました。Christov、Engzee、Hamiltonのアルゴリズムはすべて同じピークを検出できています。

一方で、GamboaのアルゴリズムとSSFアルゴリズムは、R波以外のピークやノイズまで拾ってしまったようです。SSFの方は、オプションのthreshold値を100倍にすることでかなり改善しました(それでも、29秒付近でピークを2重に検出してしまっています)。

BioSPPyを開発しているPIA Groupによるレビュー論文[2]によれば、公開データベースのECGについてChristovのアルゴリズムとSSFアルゴリズムがもっとも高い性能を発揮したようです。閾値を決める必要がないので、Christovアルゴリズムを使うのが楽そうです。

RRIと心電指標の計算

BioSPPyで可能な分析は、R波の検出までです。ECGから得られる指標は他にもある[3]ので、代表的なものを算出してみましょう。

R-R間隔(RRI)

R波から次のR波までの間隔です。単位はミリ秒(ms)を使うことが多いようです。

検出したR波ピーク同士の差分を取ることで簡単に求まります。

ts_peaks = ts[rpeaks_christov]
rri = np.diff(ts_peaks) * 1000
# [929.62697347,  912.04999288,  925.72097778,  ... ]

ただ差分を計算しただけでは、RRIの値同士が等間隔になっていないので、周波数領域での処理ができません。そこで、スプライン補間などによって等間隔(例えば1秒間隔)のデータに再サンプリングする処理をよく行います[3:1]

import scipy

spline_func = scipy.interpolate.interp1d(ts_peaks[:-1], rri, kind='cubic')
ts_1sec = np.arange(ts_peaks[0], ts_peaks[-2], 1)
rri_1sec = spline_func(ts_1sec).round(6)

スプライン補間した曲線からは、RRIが数秒の周期で変化していることが見て取れます。

心拍変動(R-R interval variability; RRV)

RRIの分散で、心拍がどの程度変動しているかを示します。RRVの指標として、SDNN・RMSSD・NN50などが提案されています

ここでは、RRIの標準偏差であるSDNNを計算してみます。

np.std(rri_1sec)
# 26.51 ms

Apple Watchには、定期的に測定した心拍データから心拍変動を推定する機能があります(RRIから求めたものではないため、厳密にはRRVではない)。心電図データを取得したのとなるべく近い時刻の心拍変動は33 msと記録されており、今回RRIから求めた26.51 msと大きくは違わない値でした。

周波数領域指標

RRIをフーリエ変換によってパワースペクトル密度(power spectral density; PSD)に変換すると、周波数領域での指標を計算できます。例えば交感神経活動と副交感神経活動のバランスを表すLF/HFなどがあります。

ただ、Apple Watchのデータでは測定時間が短すぎたのか、LF/HFを正しく算出できませんでした(低周波であるLFが全く出てきませんでした)。
心電測定のガイドラインでは、LF/HFの算出には最低限2分間が必要で、他の研究との比較のために5分間のデータを使うのが望ましいとされています[4]
Apple WatchのECG測定時間は30秒なので、この基準に遠く及びません。10回ほど連続で測定すれば時間は確保可能ですが、その間Apple Watchを何度も操作しながら、なるべく動かずに測定を続けるのは困難でしょう。

おわりに

Apple Watchで測定したECGデータをPythonに読み込み、心電波形の確認や心電指標の算出をしてみました。一度自分のECGデータをあれこれ分析してみると、iPhoneでなんとなく波形を眺めるだけよりもさらに測定できている実感が湧くのではないでしょうか。

脚注
  1. Carreiras C, Alves AP, Lourenço A, Canento F, Silva H, Fred A, et al. "BioSPPy - Biosignal Processing in Python", 2015-, https://github.com/PIA-Group/BioSPPy/ [Online].
    ↩︎

  2. Canento, Filipe, et al. "Review and comparison of real time electrocardiogram segmentation algorithms for biometric applications." Proc. 6th Int. Conf. Health Inform. (HEALTHINF). 2013. ↩︎

  3. 藤原 幸一, ヘルスモニタリングのための心拍変動解析, システム/制御/情報, 2017, 61 巻, 9 号, p. 381-386. https://doi.org/10.11509/isciesci.61.9_381
    ↩︎ ↩︎

  4. Task Force of the European Society Electrophysiology, "Heart Rate Variability: Standards of Measurement, Physiological Interpretation, and Clinical Use," Circulation, vol. 93, no. 5, pp. 1043–1065, Mar. 1996. https://www.ahajournals.org/doi/10.1161/01.CIR.93.5.1043
    ↩︎

GitHubで編集を提案

Discussion

初めまして。
丁寧な説明をしていただき、とても参考になりました。
RRIからLF/HFを算出したいのですが、以下のコードで大丈夫かと思ったのですが、エラーが返ってきます。
もし、ご存じでしたら、LF/HFの算出の仕方のコードを教えていただけますか?

$ pip install hrv
from hrv.classical import frequency_domain
results = frequency_domain(
rri=input_data,
fs=4.0,
method='welch',
interp_method='cubic',
detrend='linear'
)

よろしくお願いいたします。

そのライブラリ(hrv)は使ったことがないので、ドキュメントを読んでの一般的な回答をします。
本記事の最後でも触れたように、Apple WatchのECGデータでは低周波指標を算出するのに十分な測定時間がないため、以下の方法で正しいLF/HFが計算できているかどうかは分かりません。もしお使いになる際はライブラリの説明やコードをよくお調べください。


hrvパッケージは、RRIデータを通常の配列ではなく"RRi object"として管理するそうです。frequency_domain関数にも、RRiオブジェクトを渡さなければエラーになります。

https://hrv.readthedocs.io/en/latest/statistics.html

hrv.rriからRRiをインポートしたうえで、RRi([ {RRIデータの配列} ])のように呼び出すとRRiオブジェクトが得られます。

従って、以下のようなコードであれば動作すると思います。

from hrv.classical import frequency_domain
from hrv.rri import RRi

results = frequency_domain(
  rri=RRi(rri_1sec), # この記事の方法で算出した1秒周期のRRI
  fs=4.0,
  method='welch',
  interp_method='cubic',
  detrend='linear'
)
print(results)

(出力例)

/usr/local/lib/python3.7/dist-packages/scipy/signal/spectral.py:1966: UserWarning: nperseg = 256 is greater than input length  = 119, using nperseg = 119
  .format(nperseg, input_length))
{'hf': 144.00894487898282,
 'hfnu': 37.89223391610113,
 'lf': 236.03976166558095,
 'lf_hf': 1.6390631975252354,
 'lfnu': 62.10776608389889,
 'total_power': 535.0450720926806,
 'vlf': 154.9963655481168}

大変遅くなりましたが、無事実行できました。
とても参考になりました、ありがとうございます。

素晴らしい記事ありがとうございます!

RRIをフーリエ変換によってパワースペクトル密度(power spectral density; PSD)に変換する

この部分なのですが、私も同じことやったらやっぱり微妙なPSDになりました。で、LFHFの算出はストレスと自律神経の科学 によるとFFTではなくARでPSDを求めるっぽいです(中身10%も理解できませんでしたが、、)。なかなか手を出しづらいですが、数少ない有益な指標なのでチャレンジしてみたいと思っています。

興味深い内容で、大変参考になります。
下記のコードなのですが、figをもっと画像が良くなるようにしたいと考えています。
<<<<<<<<<<<<
plt.rcParams['figure.figsize'] = [12, 7]
ecg_data=ecg.ecg(signal=ecg_df[0], sampling_rate=512., show=True)
<<<<<<<<<<<<

上記コードに下記のコードを付け加えたのですが、うまくsaveできません。
どうすれば良いか教えていただけないでしょうか?
plt.savefig("file_name.jpeg",dpi=500,bbox_inches="tight")
files.download("file_name.jpeg")

biosppy.ecg.ecg()は内部でbiosppy.plotting.plot_ecg()というプロット用のメソッドを呼び出しているのですが、これの最後でplt.show()を実行しているのが画像が真っ白になる原因と思われます。plt.show()が終わると現在のグラフ画像を閉じてしまうので、その後plt.savefig()を実行しても同じグラフを出力できないようです。

biosppy.plotting.plot_ecg()のドキュメントを見ると、pathというオプションにファイル名を指定すると画像ファイルを保存してくれるそうです。

というわけで、こんな感じでいかがでしょうか

from biosppy.plotting import plot_ecg

plt.rcParams['figure.figsize'] = [12, 7]
ecg_data=ecg.ecg(signal=ecg_df[0], sampling_rate=512., show=False)
(ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate) = ecg_data
plot_ecg(ts=ts, raw=ecg_df[0], filtered=filtered, rpeaks=rpeaks, templates_ts=templates_ts, templates=templates, heart_rate_ts=heart_rate_ts, heart_rate=heart_rate, path='file_name.jpg')

ご連絡ありがとうございました。

下記のコードを実行してみたのですが、やはり白く抜けてしまいました。
粗くなりますが、画像を直接コピーするしかないでしょうか?
なかなか悩ましいです。

from biosppy.plotting import plot_ecg

plt.rcParams['figure.figsize'] = [12, 7]
ecg_data=ecg.ecg(signal=ecg_df[0], sampling_rate=512., show=False)
(ts, filtered, rpeaks, templates_ts, templates, heart_rate_ts, heart_rate) = ecg_data
plot_ecg(ts=ts, raw=ecg_df[0], filtered=filtered, rpeaks=rpeaks, templates_ts=templates_ts, templates=templates, heart_rate_ts=heart_rate_ts, heart_rate=heart_rate, path='file_name.jpg')

plt.savefig("file_name.jpg")
files.download("file_name.jpg")

あ、plot_ecgを実行した時点でファイルが保存されている(内部でplt.savefigを実行している)はずなので下から2行目のplt.savefig("file_name.jpg")は消してみてください

返信ありがとうございます。
plt.savefig("file_name.jpg")を消してみたのですが、今度は下記のエラーが出てしまいました。

FileNotFoundError: Cannot find file: "file_name.jpg"

file_name.jpgのところをpathにしてみたりしたのですが、同様にfileを見つけられないようです。
plot_ecgの内部で何をやっているのかよくわからず、悩みます。

plot_ecg のソースコードはこちらで読めます↓

https://github.com/PIA-Group/BioSPPy/blob/b391119e656bb5d3f6ddbd345493256867e46571/biosppy/plotting.py?_pjax=%23js-repo-pjax-container%2C div[itemtype%3D"http%3A%2F%2Fschema.org%2FSoftwareSourceCode"] main%2C [data-pjax-container]#L1160-L1267

一番最後のあたりがファイル保存関連なのですが、私も検証してみたところ拡張子をjpgに指定しても勝手にpngとして保存されてしまうコードになってしまってます。
従って"file_name.jpg"ではなく"file_name.png"が生成されているので、files.download("file_name.png")とすればファイルを取得できました。

ちなみに、Google Colabで実行されている場合は、左のファイルタブでどんなファイルが生成されているか確認できます。

ありがとうございました。
大変ご迷惑おかけしました。
無事、綺麗な画像を保存できました。

ログインするとコメントできます