低コストで心拍間隔変動を測ってみる - 続き

はじめに

前回はRaspberryPiで測定した心拍をPCに移してMatlabで処理したのだが、やっていること自体は簡単だということで、RaspberryPi上で、全てまかなうべく、あれこれぶち当たりながら、試してみた。

パワースペクトル密度（PSD；power spectral density）の計算法　

ストレスと自律神経の科学では、自己回帰モデルが推奨されているわけだが、モデルを作るのに最もあり得そうな次数の決定が必要になるらしい。赤池情報量基準（AIC）がよいらしいというのも分かった。が、実際のデータを当てはめると、この次数決定は一筋縄では行かなさそう。データごとに倍、半分どころか、それ以上の次数の差が出てくる。今回のデータをプロットさせてみても、ガタガタで単にAICが最低の次数を選んでいいのか、と疑問に思う。というわけで、あちこち探していたら、Lomb-Scargle 法というのに行き当たった。matlab だとここ。何がいいって、一様にサンプリングされていない信号をそのまま使うだけでよい。というわけで、こちらに飛びついた。心拍間隔変動の計算にも、ちゃんと使われているようだ。

RaspberryPi上に必要moduleをダウンロード

トライアンドエラーの結果、必要なのはnumpy,matplotlib,astropy の３つと判明。とにかくRaspberry Piにpython3 + numpy + matplotlibをインストールするを参照。ブランドほぼニューのRaspberryPiで効いたのは、結局以下のコマンド。
pip3 install bumpy
sudo apt-get install python3-matplotlib
pip3 install astropy

Scipyの方が有名な感じなのだけど、このデータ、この関数に限って言えば、（そして、ひょっとするとユーザーが不慣れな私だからか？）scipy.signal.lombscargleの結果は（ここには出していない他の方法を使った結果と比べてみて）ちょっと？と思えたので、astropy.stats.LombScargleを使うことにした。

あと、非常に有用だったのがQiitaのこの記事
Sudo apt-get install python3-pyqt5
backend : qt5agg
を真似したら、書けなかったグラフが描かれるようになった。ありがとうございます。

実際にグラフを書かせる

比較的安定にデータは取れるのだが、最初いくつかは荒れるので、それは整形しておくのが無難。手動で。
で、それを使って、python3 plot_hrv.py input.npy output.jpg で、図画セーブされる。

こんな感じ。データが混乱してしまったので、前回のデータとは違ったデータ。

LF、HFがパワースペクトルのLF成分の領域（0.05Hzから0.15Hzまで）とHF成分の領域（0.15Hzから0.40Hzまで）の強度。単位はsec^2/Hz。

またまた、超初心者のコードを乗せておく。**間違い等のご指摘お願いします。**って、その前に、配列の操作が直感的でないというのか、分からない。結果オーライ的。そのうち、勉強する。

plot_hrv.py

import sys 
import numpy as np
import matplotlib.pyplot as plt
from astropy.stats import LombScargle 
print('Usage - python3 ana.py input.npy output.jpg')
def main(argv):
    args = sys.argv
	
    ibi = np.load(args[1])
    ibi = ibi/1000

    a = np.shape(ibi)[0]
    time = [];
    for i in range(a):
        b = ibi[:i]
        time.append(b.sum())
    time = np.asarray(time) #list to array

    hr = 60/np.mean(ibi)

    frequency, pgram = LombScargle(time,ibi).autopower()
    print(max(frequency))
    print(max(pgram))
    
    n1 = min((frequency>0.05).nonzero())[0]
    n3 = min((frequency>0.15).nonzero())[0]
    n4 = min((frequency>0.4).nonzero())[0]
    n2 = n3-1
    n4 = n4-1
    
    lf_pow = 0
    hf_pow = 0
    for i in (n1,n2):	
        lf_pow = lf_pow + pgram[i]
    for i in (n3,n4):
        hf_pow = hf_pow + pgram[i]
    lf_pow = (lf_pow - pgram[n1]/2 - pgram[n2]/2) / ((frequency[n2] - frequency[n1])/(n2 - n1))
    hf_pow = (hf_pow - pgram[n3]/2 - pgram[n4]/2) / ((frequency[n4] - frequency[n3])/(n4 - n3))
    lf = round(lf_pow,2)
    hf = round(hf_pow,2)
     
    print("LF_POW = :", lf_pow)
    print("HF_POW = :", hf_pow)
    print("HR = :",hr)
    plt.subplot(121)
    plt.plot(time,ibi,lw=0.5)
    plt.title("Time-RRI")
    plt.subplot(122)
    plt.plot(frequency,pgram,lw=0.5)
    plt.xlim(0,0.5)
    plt.ylim(0,0.2)
    plt.title("Freq-Pow: LF = " + str(lf) + " HF = " + str(hf))
    plt.savefig(args[2])
    plt.show()

if __name__ == "__main__":
    main(sys.argv[1:])

考察というよりなんだかごちゃごちゃと

1. なんとか、RaspberryPiのみで心拍計測、HRVの計算をすることが可能になった。
2. ここでは載せていないが、PSDを求めるべきいくつかの方法をトライしても、ある個人に限って言えば、ストレス指標であるLF,HF,LF/HFの状態による大小関係は変わらなかった。つまり、個人Aの状態X,YでHFがX>Y,LFがXYなら、どの方法でもこの大小関係は一緒だった。（って、n = 3 だけね。）
3. 一方、個人間の大小は、方法によって簡単に変わるみたいだ。（n = 1 で変わってしまった。）
4. ので、一人の人間の状態を見るのにはある程度は使えるかもしれない。
5. 私の携帯の充電器でRaspberryPiに電力供給が可能だった。あとは電力消費の少ないスクリーンがあれば、持ち運べるはず。探す。