所要時間
おそらく長くても1時間ほど
おそらく心拍ゆらぎを用いたストレス計算のプログラムの中では最も簡単でかなり高速
精度は多少劣るかもしれないが遊ぶには十分
前提条件
心拍数とストレスに相関があると言われていることを知っている。
自律神経評価について多少は調べている。
LF/HFが何かわかっている
たぶんLF/HFについて調べてからの方がいいと思います。
概要
心拍ゆらぎにおける低周波と高周波の割合で自律神経評価をする。
おおまかに、低周波が多く、高周波が少ない時に大きなストレスを感じていることが知られている
低周波のゆらぎは交換神経の作用が働いている可能性があることを利用する。
データについて
ArduinoやRaspberryPiと心拍センサを用いた計測
CSVファイル形式で保存することにする。(※今回は、"heart_rate_data.csv")
基本的には1行目に変数名、2行目からはデータとする。
データは1列目が時間t[s]、2列目がRRI(R波とR波の感覚)[s]
RRIは、拍動1回にかかる時間と考えれば良い。
つまり、心拍数≒60÷RRI
データの読み込み
data <- read.csv("~/heart_rate_data.csv",header = 1)
time_math <- x$time
as.vector(time_math)
rri_math <- x$rri
as.vector(rri_math)
簡単に言うとCSVファイルの読み込み→ベクトル化という流れ。
ベクトル化をしないと後の計算が不可能となるので要注意
1行目の文字データを数字として読み込ませない工夫
read.csv("ファイル名",header=1) で、1行目を変数名として設定する
※ファイル名が日本語だと稀におかしな挙動をするので注意
※同様にディレクトリ名も日本語を使用するのはオススメしない
ベクトル化
as.vector(変数名)
を用いる。Rは変数の宣言がないのでしっかり型を指定してあげる。
たぶんこれをやらないとdaraframeかlistあたりになってしまう。
とりあえずプロットする
plot(time_math,rri_math)
lines(time_math,rri_math)
plotで点を描画できる。linesでは線を描画できる。
2次元データの表示は非常に簡単。
データの補間
自分で作ったセンサとかの場合は誤差がひどいけどデータが少ないよりはマシ
ということでかなりのデータのかさましを行う。
今回は2分程度の計測なので100個ほどのデータだがこれを2000個ほどに拡張する
(このあと、周波数領域の評価を行うので2のn乗個のデータにしておく)
実際、5分ほど計測した後処理をするのであればこの部分は省略できる。
sp <- smooth.spline(time_math, rri_math)
i <- seq(0,125000, length=2048)
pred <- predict(sp,i)
一応表示してみる
lines(pred)
たぶんひどいデータになる。もし結果が怪しいのであれば補間の部分は省略した方が良いかもしれない。
Rでは補間も非常に簡単なのでデータ点数をかさ増しする場合には非常に便利。
ここから重要
自己回帰モデルを使ってみる
自己回帰モデルとは
心拍間隔のデータは時系列なので自己回帰モデルを使う
スペクトル分析ができれば何でもいいけど・・・
FFTでは、周期関数ではないデータのスペクトル分解が出来ない
また、例え出来たとしてもフーリエ変換は比較的時間がかかる(Rでは)
これらを踏まえると
自己再帰モデルが良い(もちろんほかのやり方もある)
spec <- spec.ar(pred$y, plot = FALSE)
lf_data <- 0
hf_data <- 0
spec_lf <- spec$spec[50:150] #0.05Hz to 0.15Hz ... LF
spec_hf <- spec$spec[150:350] #0.15Hz to 0.35Hz ...HF
今回は0.05-0.15HzをLF、0.15-0.35HzをHFとする。
今回のプログラムでは、spec[i]に、0.01*i[Hz]のスペクトル密度が格納される。
Rは[j:k]で、j番からk番までのデータを抜き出せる。
一応、パワースペクトル密度の公式と合わせるために2乗する。
ぶっちゃけ2乗しなくても"そこそこ"の評価はできる
spec_lf <- spec_lf * spec_lf
spec_hf <- spec_hf * spec_hf
ベクトル計算
今回はRのいいところ、ベクトル計算を使っている
当たり前のようにspec_lf*spec_lfと書いているが、これをC言語風に書くと
for(int i = 0;i < N; i++){
spec_lf[i] = spec_lf[i]*spec_lf[i];
}
こんな感じ。C言語風に解説すると、同じ番号の配列をそのままいっきに掛け算できる。
このベクトル計算が簡単にできるのがRのいいところ
ベクトル計算ができなかったら正直恐ろしいほど遅いプログラムになる
さらに素晴らしいことに、C言語風に言うと配列の要素を全部一気に足せる
とても簡単
lf_data <- sum(spec_lf) #sum spectrum
hf_data <- sum(spec_hf) #sum spectrum
ストレス指標LHの計算
ストレス指標はLHとする。
LHは、全体における低周波(LF)と高周波(HF)の割合なので
lh <- lf_data / hf_data #Calculate LH Ratio
で計算できる。このLHが大きいとストレスが大きいと言われている。
一応いらない変数は毎回消しておく
rm(hf_data,lf_data,spec_hf,spec_lf,spec,data)
最後にLHの値を表示する。
lh
おわりに
非常に簡単にストレスチェックができた
AR法を用いるとかなり高速化できる
ほかにもたくさんの自律神経評価のプログラムがあるがかなり高速化できているはず。