LoginSignup
16
17

More than 3 years have passed since last update.

@Hirosaji

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

Last updated at Posted at 2018-03-01

はじめに

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

（Hexoskin APIの概要と事始めについては、こちらの記事をご参照ください）

ストレス指標LF/HFとは

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

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の推定モデルは、おそらくこんな感じです。

\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を平滑化することができます（中式）。

以上の計算の性質上、時系列データの異常値や欠損値に(少し)強い特性があるため、今回はWelch法を採用しました。

ちなみに、Scipyのsignal.welch関数を使えば、簡単にWelch法を利用することができます。

LF/HFを算出するスクリプト

今回書いたスクリプトは、Hexoskin公式のPythonクライアントを利用したものです。
以下のURLからPythonクライアントをcloneまたはdownloadしておいてください。

実行環境（動作確認済み）

  • OS: macOS High Sierra
  • python: 2.7.14

事前準備①: 外部ライブラリのインストール

必要なライブラリをインストールします。

Terminal
$ pip install requests urllib numpy scipy

事前準備②: Hexoskinの認証情報をまとめる

今回は、認証情報などを外部ファイルにまとめて管理します。
各認証情報をまだ用意してない場合は、こちらの記事を参考に取得してください。

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を算出することができます。

以下、ソースコード。

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
$ python get_lf_hf.py
https://api.hexoskin.com
Downloading hr_quality
Downloading rrinterval
LF/HF = 1.11

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

算出した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以上 要注意

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

終わりに

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

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

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

16
17
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
16
17