MATLABで心電図を分析しよう Ep. 1
Introduction
私は心電図が好きです．心電図から得られるただの値が，自分の身体の活動を表現していると思うとわくわくします．ただそれだけです．この気持ちをこの記事を読んでいただいている方にも知ってほしいという思いでこの記事を書いています．私自身もこの道を極めているわけではないので，至らない部分があると思いますがその時はご指摘いただければありがたいです．
Goal
- よくみる心電図波形をプロットする
- 心拍の検知
- 時間領域の解析で得られるパラメータの導出
- To Be Determined
大まかな目標としてこれらを挙げますが，4以降は未定です．今回は1を取り上げます．2以降のトピックは次回の記事に取り上げる予定です．
MATLAB
MATLABは数値解析などに用いられるソフトウェアで，私もよく使います．雑な説明にはなってしまいますがMATLABを使えばなんでもできます．MATLABを使うと部屋が綺麗になったり，人間関係も良好に．万事(よろず)のイベントに効果があります(個人の意見です)．今回のような信号処理も例外ではありません．ただし配列などのインデックスが0からでなく1から始まることには注意ですね．
Pythonを使おうとも考えましたが，別の機会にやります．
心電図データ
心電図のデータをなかなか自分で測定，取得することはないと思いますが，今回は私の心電図データを用います．私のGithubにECG_dataというリポジトリがあるので，そこからテキストデータとして入手できます(リポジトリ)．2列目が心電図のデータでそれ以外の列は不要なデータです．Kaggleのデータセットにも心電図データのようなものがありました(ECG Heartbeat Categorization Dataset)．測定機器としてOpenBCI Ganglionというバイオセンシングボードを使用しました(Ganglion)．測定方法は一般的な方法である第Ⅱ誘導，サンプルレートは200.0Hzです(つまり1秒間に200個のデータが測定できます)．第Ⅱ誘導という方法はQRS波が大きく現れるのでよく使われるようです(心電図測定方法)．
Figure 1. wikipediaより[1]
じゃあ，やってみましょう
前置きが長くてすみません，おしゃべり好きなんですよ，いつも独り言ばかりで．
それはそうと，炊飯器で簡単に焼き芋というか蒸し芋を作れるということを最近知って，作ってみたらおいしかったです．焼き芋買わなくてもいいんじゃないのかなと思うくらい．
1. テキストデータの読み込み
MATLABにはファイルからデータを読み込むための関数が色々あります．今回はテキストデータを読み込むためにreadtable関数を使いましたが，読み込めればなんでもいいです．変数file_nameは任意のファイル名に変更してください．
file_name = "ecg.txt";
T = readtable(file_name);
%% 全てのデータではなく，2列目のデータが欲しい
raw_ecg = T.Var2;
ここで出てくるT.Var2が今回欲しい心電図のデータです．というのも，MATLABのtable型というのはカラム名がないとVar1, Var2...のようにカラム名をつけてくれるので，それに従った形になります．
これでとりあえず生データは読み込めました．小さい範囲でプロットしてみましょう．(配列インデックスを指定する時，[ ]ではなく( )でくくるというのもMATLABの特徴ですね[2])
plot(raw_ecg(50000:50500));
Figure 2.
まあこんな感じです．見栄えのために綺麗な部分を選んで50000個目から50500個目のデータをプロットしています．こうしてみるとすでにみたことあるような心電図ですね．y軸は振幅(mV)で，x軸は何番目の値かを表しています．人に見せる時はFigureにちゃんとタイトルとか軸にラベルをつけないとですね．上記のFigure 2は反面教師としてはいい例です．ちなみに最初の20000個らへんまでのデータが安定して測定できていないので捨てちゃってもOKです．
raw_ecg = raw_ecg(20001:end);
こんな感じ．私も捨てちゃいました．
ところで説明のために，"xx個目のデータ"と書いていますが，もう少し詳細に表現すると200Hzのサンプルレートなので200個のデータで1秒間ということになります．例えば，不安定なデータとして1から20000個目までを使わないようにしましたが，つまりこれは0~100秒までのデータを使わないようにしたということです．
2. データの前処理
心電図の生データには多くのノイズが含まれています．例えば筋肉の活動によるノイズ(筋電)や，商用電源などから発生する50Hzや60Hzのノイズ(ハムノイズ)などです．綺麗な心電図を得るにはこれらを取り除いたり小さくしたりする必要があります．それを達成するには様々な方法があると思いますが，今回はお手軽に使えるbandpass関数[3]を使います．
簡単に説明すると，指定した周波数の範囲以外の信号を減衰させる(小さくする)関数です．ここで気になることが，欲しい心電図の情報もその処理によって小さくなってしまうのではないか，ということです．ではこれを明確にします．
Figure 3. Research Gateより[4]
これは生体信号のスペクトル密度を表現しています．砕いて言うと，信号がどんな周波数で作られていて，その周波数がどのくらいの割合で信号を構成しているのか，を表しています[5]．
それを踏まえた上でFigure 3を見てみましょう．ECGが心電図の成分です．3~5Hzあたりから30~35Hzあたりが主な周波数帯なのが分かりますね．つまりその周波数の範囲以外を小さくする処理をすればいいということになります．
passed_ecg = bandpass(raw_ecg, [3 35], 200);
引数として左から順に，生データ，残しておきたい周波数帯，サンプルレートです．
結果をプロットしてみます．
plot(passed_ecg(29500:30000));
% 綺麗なところを選びました
Figure. 6[6]
これはもうよくみる心電図になったんじゃないでしょうか．なりました．なりましたね．QRS波もはっきりと表れています．これでひとまず目標としてあげていたうちの"1. よくみる心電図波形をプロットする"は達成できました．今回はここまでとします．コードはGithubのECG_dataとは別のリポジトリにあります(リポジトリ)．
Conclusion
コードの量は少ないですが，いろいろと処理ができるのがMATLABです．便利ですね．心電図の生データにちょっと手を加えるだけで比較的綺麗な波形を得ることができます．病院にあるような心電図計測器(心電計)はもっとたくさん処理をしていると思います，たぶん比にならないくらい．それはそうと心電図を測ってそれを読む医療系のお仕事をされている方々はすごいですね，信号を読むって個人的にかっこいいと思っています．だんだん話が逸れそうなので一旦お開きにします．次回は目標の2以降を取り上げます．それでは良い1日を．
-
ちなみに[ ]はどう呼ぶべきなんですかね．わからずにいつも四角かっことかスクエアブラケットとか呼んでいます ↩︎
-
https://www.researchgate.net/figure/Relative-power-spectra-of-QRS-complex-P-and-T-waves-and-muscle-noise-and-motion_fig2_275260452 ↩︎
-
適切な表現ではないかもしれないので，より深掘りしたい方は"パワースペクトル密度"でググってみてください． ↩︎
-
人に見せる時はFigureにちゃんとタイトルとか軸にラベル(略) ↩︎
Discussion