はじめに
修士論文の研究で、ウェアラブルセンサからストレス指標を取得する必要がありました。
そこで、スマートシャツHexoskinを購入し、Hexoskin APIを利用してストレス指標の一つであるLF/HFを計算するスクリプトを書きました。
今回はそのとき作ったスクリプトと、その周辺知識について、解説していきます。
(Hexoskin APIの概要と事始めについては、こちらの記事をご参照ください)
ストレス指標LF/HFとは
LF/HF(Low-Frequency/High-Frequency)とは、交感神経の活性度を表す指標で、多くの研究や計測機器でストレス指標として用いられています。
詳しくは、以下のサイトが参考になります。
- ストレスと自律神経の科学 - http://hclab.sakura.ne.jp/
LF/HFの算出方法
Hexoskinを使って、ストレス指標LF/HFを算出する方法は、以下の通りです。
- Hexoskin APIから心拍間隔(RRI;R-R Interval)の時系列データを取得する
- RRI時系列データからパワースペクトル密度(PSD;Power Spectral Density)を推定する
- 推定したPSDの低周波成分(LF成分, 0.04~0.15Hz)と高周波成分(HF成分, 0.15~0.45Hz)を算出する
- LF成分とHF成分の比であるLF/HFを算出する
また今回は、2.のPSDの推定にWelch法(平均修正ピリオドグラム法)を用いました。
メモ:Welch法によるPSDの推定
Welch法によるPSDの推定モデルは、おそらくこんな感じです。
\tilde{P_{xx}}(f) = \frac{1}{S} \sum_{t=0}^{S-1} \tilde{P_{xx}^{(t)}}(f) \\
\tilde{P_{xx}^{(t)}}(f) = \frac{1}{NU} \begin{vmatrix}\sum_{m=0}^{N-1}x(m)\,e^{-i \frac{2\pi km}{N}}\end{vmatrix}^2 \\
U = \frac{1}{N} \sum_{m=0}^{N-1} \begin{vmatrix}w(m)\end{vmatrix}^2 \\
\begin{pmatrix}
S = セグメント数, \;\; \\
N = サンプル数, \; \\
\hspace{10pt} m = 時間のインデックス \, (0 \leq m \leq N-1), \hspace{10pt} \;\;\;\;\;\;\;\; \\
\hspace{10pt} k = 周波数のインデックス \, (0 \leq k \leq N-1), \hspace{10pt} \;\;\;\;\;\;\;\; \\
w(m)=窓関数 \, (今回はハニング関数を利用) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{pmatrix}
Welch法は、時系列データから複数のセグメントを分割し、各セグメントについて修正ピリオドグラム法でPSDを推定し、最後にそれらの平均を計算する手法です(上式)。
修正ピリオドグラム法は、通常のピリオドグラムにおける計算の前処理として、時系列データに対して窓関数を用いた正規化(下式)をする手法で、得られるPSDを平滑化することができます(中式)。
- 参考資料
- スペクトル密度 - wikipedia
- Welch's method - wikipedia
- Percival D. B. and A. T. Walden, 1993: Spectral analysis for physical applications: Multaper and
conventioinal univariate techniques. Cambridge University Press, pp. 583.
以上の計算の性質上、時系列データの異常値や欠損値に(少し)強い特性があるため、今回はWelch法を採用しました。
ちなみに、Scipyのsignal.welch関数を使えば、簡単にWelch法を利用することができます。
- scipy.signal.welch — SciPy v0.14.0 Reference Guide
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html
LF/HFを算出するスクリプト
今回書いたスクリプトは、Hexoskin公式のPythonクライアントを利用したものです。
以下のURLからPythonクライアントをcloneまたはdownloadしておいてください。
- Pythonクライアント - https://github.com/lguerdan/HexoBio
実行環境(動作確認済み)
- OS: macOS High Sierra
- python: 2.7.14
事前準備①: 外部ライブラリのインストール
必要なライブラリをインストールします。
$ pip install requests urllib numpy scipy
事前準備②: Hexoskinの認証情報をまとめる
今回は、認証情報などを外部ファイルにまとめて管理します。
各認証情報をまだ用意してない場合は、こちらの記事を参考に取得してください。
{
"HexAPIkey": "",
"HexAPIsecret": "",
"HexUserName": "",
"HexPassword": ""
}
ファイル構造
HexoBio
├── credentials.json
├── get_lf_hf.py
├── HxApi2_0.py
└── ...
ソースコード
今回はget_lf_hf.py
というスクリプトを用意しました。
get_lf_hf.py
は、Hexoskin APIからRRI時系列データを取得するmain関数と、RRI時系列データを引数にしてストレス指標LF/HFの値を返すコールバック関数とで構成したPythonファイルです。
このスクリプトを使えば、直近のActivityのLF/HFを算出することができます。
以下、ソースコード。
# -*- coding: utf-8 -*-
import HxApi2_0
import urllib, json
import numpy as np
from numpy import arange, cumsum
from scipy.signal import welch
from scipy.interpolate import splrep, splev
def main():
# JSONから認証情報を読み込む
with open("credentials.json", "r") as f:
credentials = json.load(f)
username = credentials["HexUserName"]
password = credentials["HexPassword"]
publicKey = credentials["HexAPIkey"]
privateKey = credentials["HexAPIsecret"]
# userIDを取得
userURL = "https://" + username + ":" + password + "@api.hexoskin.com/api/user/"
userInfo = urllib.urlopen(userURL).read()
userInfoJSON = json.loads(userInfo)
userID = userInfoJSON["objects"][0]["id"]
# OAuth2.0認証を行なう
auth = HxApi2_0.SessionInfo(publicKey=publicKey,privateKey=privateKey,username=username,password=password)
# 記録したActivityのIDを取得
records = HxApi2_0.getRecordList(auth, limit = "2000")
all_recordID = [record['id'] for record in records]
last_recordID = all_recordID[len(all_recordID)-1] # 直近のActivityのID
# 直近のActivityのデータを取得
dataset = HxApi2_0.getRecordData(auth,recordID=last_recordID)
# RRIを含んだ時系列データをRRI_2_LFHF関数に渡してLF/HFを計算させる
lf_hf = round(RRI_2_LFHF(dataset), 2)
print "LF/HF = " + lf_hf
def RRI_2_LFHF(dataset):
# データからRRI時系列データを抽出
rri = [data[1] for data in dataset["rrinterval"]]
ibi = np.array(rri)
ibi = ibi * 1000
# Welch法によるPSD推定
Fxx, Pxx = welch(ibi, fs=4.0, window='hanning', nperseg=256, noverlap=128)
# パラメータを設定
vlf = 0.04
lf = 0.15
hf = 0.4
Fs = 250
# LFとHFの周波数の範囲に対応するインデックスを設定
lf_freq_band = (Fxx >= vlf) & (Fxx <= lf)
hf_freq_band = (Fxx >= lf) & (Fxx <= hf)
# LFとHFの曲線の下の面積を計算
dy = 1.0 / Fs
LF = np.trapz(y=abs(Pxx[lf_freq_band]), x=None, dx=dy)
HF = np.trapz(y=abs(Pxx[hf_freq_band]), x=None, dx=dy)
LF_HF = float(LF) / HF
return LF_HF
if __name__ == "__main__":
main()
実行結果の例がこちら。
$ python get_lf_hf.py
https://api.hexoskin.com
Downloading hr_quality
Downloading rrinterval
LF/HF = 1.11
この例では、「LF/HF = 1.11」という結果が得られました。
- 参考リンク
- Python: SciPy のパワースペクトル密度推定の関数
https://org-technology.com/posts/power-spectral-density.html - PSD estimation of a RRi serie with Python
https://rhenanbartels.wordpress.com/2014/04/06/first-post/
算出したLF/HFの解釈
ストレス指標LF/HFはあくまでも指標なので、この指標をどう解釈するかは、さらなる検討が必要です。
先行研究でも、デスクワークや運転など、与えた課題によってLF/HFの判断基準は異なっていました。
- Healey, J. A., Picard, R. W.:Detecting stress during real-world driving tasks using physiological sensors. IEEE Transactions on intelligent transportation systems, Vol. 6, No. 2, pp. 156-166, 2005.
- Hjortskov, N., Rissén, D., Blangsted, A. K., Fallentin, N., Lundberg, U., Søgaard, K.:The effect of mental stress on heart rate variability and blood pressure during computer work. European journal of applied physiology, Vol. 92, No. 1-2, pp. 84-89, 2004.
これら以外にもLF/HFをストレス指標とした研究は数多くあるので、まずは目的に沿った巨人の肩を探してみると良いでしょう。
ちなみに下記の特許では、LH/HFの値を基準に、疲労・ストレス度を以下の表のように定量的に判定しています。
LF/HF | 判定 |
---|---|
0.0〜2.0 | 良好 |
2.0〜5.0 | 注意 |
5.0以上 | 要注意 |
- 疲労度の判定処理システム
https://patents.google.com/patent/JP5491749B2/ja
この特許は総合的な意味での疲労・ストレス度を判定していると謳っているので、巨人の肩が見つからない時は、便宜的に上記の基準を参考にしても良いかもしれません。
終わりに
この記事では、Hexoskinでストレス指標LF/HFを算出する手順と数式、そしてソースコードを示しました。
ちなみにストレス指標はLF/HF以外にも、アミラーゼ値やフリッカー値などがストレス研究史上よく使われてきました。
それらの他にも、色々と提唱はされているみたいです。
正直どの指標が有意な結果をもたらすかは、実際に測って比べてみないと分からないので、色々計測しておくのが吉でしょう。