Hexoskinでストレス指標LF/HFを算出する

# はじめに
修士論文の研究で、ウェアラブルセンサから**ストレス指標**を取得する必要がありました。
そこで、スマートシャツHexoskinを購入し、Hexoskin APIを利用して**ストレス指標**の一つである**LF/HF**を計算するスクリプトを書きました。
今回はそのとき作ったスクリプトと、その周辺知識について、解説していきます。

（Hexoskin APIの概要と事始めについては、[こちら](https://qiita.com/Hirosaji/items/edb8ce323603371110db)の記事をご参照ください）

# ストレス指標LF/HFとは
LF/HF（Low-Frequency/High-Frequency）とは、交感神経の活性度を表す指標で、多くの研究や計測機器でストレス指標として用いられています。
詳しくは、以下のサイトが参考になります。

- ストレスと自律神経の科学 - http://hclab.sakura.ne.jp/

# LF/HFの算出方法
Hexoskinを使って、ストレス指標LF/HFを算出する方法は、以下の通りです。

1. Hexoskin APIから心拍間隔（RRI；R-R Interval）の時系列データを取得する
2. RRI時系列データからパワースペクトル密度（PSD；Power Spectral Density）を推定する
3. 推定したPSDの低周波成分（LF成分, 0.04～0.15Hz）と高周波成分（HF成分, 0.15～0.45Hz）を算出する
4. LF成分とHF成分の比であるLF/HFを算出する

また今回は、2.のPSDの推定にWelch法（平均修正ピリオドグラム法）を用いました。

## メモ：Welch法によるPSDの推定
Welch法によるPSDの推定モデルは、おそらくこんな感じです。

```math
\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](https://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%9A%E3%82%AF%E3%83%88%E3%83%AB%E5%AF%86%E5%BA%A6)
 - [Welch's method - wikipedia](https://en.wikipedia.org/wiki/Welch%27s_method)
 - 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<br>
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

## 事前準備①: 外部ライブラリのインストール
必要なライブラリをインストールします。

```terminal:Terminal
$ pip install requests urllib numpy scipy
```

## 事前準備②: Hexoskinの認証情報をまとめる
今回は、認証情報などを外部ファイルにまとめて管理します。
各認証情報をまだ用意してない場合は、[こちらの記事]()を参考に取得してください。

```json:credentials.json
{
	"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を算出することができます。

以下、ソースコード。

```Python:get_lf_hf.py
# -*- 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()
```

実行結果の例がこちら。

```terminal:Terminal
$ python get_lf_hf.py
https://api.hexoskin.com
Downloading hr_quality
Downloading rrinterval
LF/HF = 1.11
```

この例では、「LF/HF = 1.11」という結果が得られました。

- 参考リンク
 - Python: SciPy のパワースペクトル密度推定の関数<br>
https://org-technology.com/posts/power-spectral-density.html
 - PSD estimation of a RRi serie with Python<br>
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以上 | 要注意 |

- 疲労度の判定処理システム<br>
 https://patents.google.com/patent/JP5491749B2/ja

この特許は総合的な意味での疲労・ストレス度を判定していると謳っているので、巨人の肩が見つからない時は、便宜的に上記の基準を参考にしても良いかもしれません。

# 終わりに
この記事では、**Hexoskin**で**ストレス指標LF/HF**を算出する手順と数式、そしてソースコードを示しました。

ちなみに**ストレス指標**は**LF/HF**以外にも、アミラーゼ値やフリッカー値などがストレス研究史上よく使われてきました。
それらの他にも、色々と提唱はされているみたいです。

- [えっ、こんな部位でストレス測定？-息、唾液、音](http://healthcare-biz.jp/2017/07/%E3%81%88%E3%81%A3%E3%80%81%E3%81%93%E3%82%93%E3%81%AA%E9%83%A8%E4%BD%8D%E3%81%A7%E3%82%B9%E3%83%88%E3%83%AC%E3%82%B9%E6%B8%AC%E5%AE%9A%EF%BC%9F-%E6%81%AF%E3%80%81%E5%94%BE%E6%B6%B2%E3%80%81%E9%9F%B3/)

正直どの指標が有意な結果をもたらすかは、実際に測って比べてみないと分からないので、色々計測しておくのが吉でしょう。