概要

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

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

• applewatch_ecg_python - Colaboratory

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, interactive=False)

概要プロットは、左列が上から生の波形・フィルタリングされた波形と検出した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でなんとなく波形を眺めるだけよりもさらに測定できている実感が湧くのではないでしょうか。