ヘルスケアロボットの開発 - yana_lab

ヘルスケアロボットの開発

緒言

　厚生労働省の調査結果より, 仕事や職業生活でストレスを感じる人の割合は50% を超える. ストレスの解消を行わず無理をして仕事を続けてしまうと, うつ病などのメンタル不調に陥る. そのため,ストレスが溜まったことを早めに認識し対処する必要がある.このような問題に対してストレスを評価する研究が盛んに行われている.しかしこれらの研究では, ユーザ自身のストレス状態を知らせるのに時間がかかるという問題がある. そこで短時間にフィードバックしユーザにストレスの状態を意識させることが必要である.本研究では, 脈波を用いてストレス状態をフィードバックするヘルスケアロボットのためのストレス評価方法について検討を行う.

脈波とストレスの関係

人が精神的や肉体的にストレスを感じた時、大脳が刺激され交感神経に影響する．結果として、心拍数が上昇し、動悸を感じたり血圧の上昇が起こる．つまり、脈波を計測し、心拍数を測定することでストレスを感じているかどうかを判断できると考えられる．

ローレンツプロット

ローレンツプロットは、n秒のデータをx軸，n+1秒をy軸に信号をプロットする．プロットの面積を求めることで，信号のばらつき（カオス性）を求めることができる．
脈波には、心臓が血液を送り出す際に生じるR波という波形が存在する。n回目のR波とn+1回目のR波の間隔のことをRRI(R-R Interval)という．RRIの間隔はリラックス時はばらつきが大きいが、ストレス負荷時は心拍数が増加し、ばらつきが小さくなる。
このことから、RRIを計測し、ローレンツプロットでプロット、プロットの楕円面積を求めることで、ストレスを感じているかどうかを判断することができる．リラックス時はRRIのばらつきが大きいため楕円面積が大きくなり、ストレス負荷時はばらつきが小さく、楕円面積が小さくなる．

唾液アミラーゼ

ストレスを評価する別の方法として唾液のアミラーゼ量を測定する方法がある．唾液内のアミラーゼ量が少ないとリラックスしており、多いとストレスを感じていることが知られている．

ロボットの製作

ストレスが溜まった時にユーザに伝えるロボットの製作を行った。
製作過程は図を用いて説明をしているため,以下の手順でダウンロードをする。
NasNavi＞LS210D87Eをクリック＞綱島研究室＞栁澤研ロボット班2018＞ロボットの製作

ここでは,ロボットに動作を持たせるためのサーボモータのプログラミングを以下に示す。今回はマイクロコンピュータのESP32を使用したプログラミングとなっている。

サーボモータのプログラミング [ESP32]

#include "esp_system.h"   //esp32を使用することを宣言

//初期設定
//パルス幅を変えることでサーボモータの角度を変えることが出来る
 
int min = 26;   //パルス幅の計算 (26/1024)*20ms ≒ 0.5 ms  (角度0°)
int max = 123;  //パルス幅の計算(123/1024)*20ms ≒ 2.4 ms  (角度180°)
int b = 74;     //パルス幅の計算 (74/1024)*20ms ≒ 1.45 ms (角度90°)
int n = min;    //minをnに置き換える
 
 
void setup() {
  Serial.begin(9600); //9600bpsで通信
 
  // ピンのセットアップ
  ledcSetup(0, 50, 10); // 0ch:高速チャンネル 50Hz 10bit:0～1023
  // ピンのチャンネルをセット
  ledcAttachPin(15, 0); // 15pin:15番ピン, 0ch:高速チャンネル
  // PWM出力を行う.analogWrite()の代用ができる.
  ledcWrite(0, b);      // 0ch:0チャンネル 最初にサーボホーンを90°回転させる
  delay(10000);         //10秒待つ
 }
 
void loop() {
  delay(9000);       //9秒待つ
  ledcWrite(0, max); //0度にする
  delay(1000);       //1秒待つ
  ledcWrite(0, n);   //180°にする
  delay(1000);       //1秒待つ
  ledcWrite(0, b);   //90°にする
  delay(1000);       //1秒待つ
  ledcWrite(0, max); //0度にする
  delay(1000);       //1秒待つ
  ledcWrite(0, n);   //180°にする
  delay(1000);       //1秒待つ
  ledcWrite(0, b);   //90°にする
  delay(9000);       //9秒待つ
 }
 

心拍センサの仕組み

心拍センサには，安価な反射型脈拍センサ(https://pulsesensor.com/)を使用する．反射型脈拍センサは緑色波長(565[nm])の光を照射し，生体内を反射した光を計測する．動脈の血液内の酸化ヘモグロビンには，光を吸収する特性があり，心臓の脈動によって動脈の血液量が変化するため，生体内を反射する光の変化量から脈拍の測定を行うことができる．

プログラム

心拍数の測定方法

心拍数を計測するには、R波を検出することができれば求めることができる。R波の検出ができれば、RRIを求めることができるため、RRIを移動平均し、60[s]÷RRI[s]＝心拍数[回]から心拍数を求めることができる．

心拍数の測定 [Arduino]

#define Pulse_pin 25

volatile int rate[60];
volatile unsigned long sampleCounter = 0;
volatile unsigned long lastBeatTime = 0;
volatile unsigned long StressCheckTime = 0;
volatile unsigned long TotalTime = 0;   q
volatile int P = 1600;
volatile int T = 1600;
volatile int thresh = 1800;
volatile int amp = 0;
volatile boolean firstBeat = true;
volatile boolean secondBeat = false;
volatile boolean Pulse = false;
volatile int BPM = 0;
volatile int Signal;
volatile int RRI = 800;
volatile int On  = 0;
volatile int Stress  = 0;
 
volatile int interruptCounter;
hw_timer_t * timer = NULL;
portMUX_TYPE timerMux = portMUX_INITIALIZER_UNLOCKED;
 
void IRAM_ATTR onTimer() {/* this function must be placed in IRAM */
  /* we will do it inside a critical section ,
     these marcos use portMUX_TYPE variable */
 
  portENTER_CRITICAL_ISR(&timerMux);
  interruptCounter++;
  portEXIT_CRITICAL_ISR(&timerMux);
 
  /* other handling code if needed */
 

void setup() {

 // initialize serial communication at 9600 bits per second:
 Serial.begin(9600);
 pinMode(Pulse_pin,INPUT);
 
 /* set prescaler 80 -> tick every 1us (80MHz / 80) */
 timer = timerBegin(0, 80, true);
 /* set the function name to call when interruption occured */
 timerAttachInterrupt(timer, &onTimer, true);
 /*  set timer-interruption-interval to 1000000 counts -> every 1s */
 timerAlarmWrite(timer, 50000, true);
 /* enable timer */
 timerAlarmEnable(timer);
 

}

void loop() {

 if (interruptCounter > 0) {
       
   Signal = analogRead(Pulse_pin);
 
   sampleCounter += 50;
   TotalTime     += 50;
   StressCheckTime += 50;
   int N = sampleCounter - lastBeatTime;

   if(Signal < thresh && Signal < T){
     T = Signal
   }

   if(Signal > thresh && Signal > P){
     P = Signal;
   }

   if (N > 250)
     if (Signal > thresh && Pulse == false){
       Pulse = true;
       On = 1000;
       RRI = sampleCounter - lastBeatTime;
       lastBeatTime = sampleCounter;

       if(secondBeat){
         secondBeat = false;
         for(int i=0; i<=59; i++){
           rate[i]= RRI;
         }
       }

       if(firstBeat){
         firstBeat = false;
         secondBeat = true;
         return;
       }

       word runningTotal = 0;

       for(int i=0; i<=58; i++){
         rate[i] = rate[i+1];
         runningTotal += rate[i];
       }

       rate[59] = RRI;
       runningTotal += rate[59];
       runningTotal /= 60;
       BPM = 60000/runningTotal;
     }

     if (Signal < thresh && Pulse == true){
       Pulse = false;
       amp = P - T;
       thresh = (amp*2)/3 + T;
       P = thresh;
       T = thresh;
     }

     if (N > 30000){
       thresh = 1800;
       P = 1600;
       T = 1600;
       lastBeatTime = sampleCounter;
       firstBeat = true;
       secondBeat = false;
     }

   if(On == 1000){
     On = 0;
   }

   delay(50);
 }

}

}

心拍数の測定 [Python]

今回センサはArduinoにつなぎ、それをパソコン上でシリアル通信を用いてpython上で処理を行っている。

import numpy as np
import serial
import time
import csv
from matplotlib import pyplot
 
######################################  変数設定  ######################################

ser=serial.Serial('COM3',9600)                    ## Serial通信先の設定

a = 100                                           ## センサ値の移動平均をとる個数
b = 10                                            ## 差の移動平均をとる個数

############ ↓Arduinoの計測プログラムのサンプリング時間を変更したら変える↓ ###########

########  注意1：t1 / t3が上の移動平均の個数a,bより小さくならないようにすること  #######
########  注意2：t1,t2,t4,60 > t3 になるようにすること  　　　　　　　　　　　　 #######

t1 = 2                                            ## 閾値設定のサンプリング時間 [秒]
t2 = 60                                           ## 心拍測定のサンプリング時間 [秒]
t3 = 0.01                                         ## arduinoの心拍センサのサンプリング時間 [秒]
t4 = 30                                           ## センサの値が0の時に停止するまでの時間設定 [秒]

######################################  関数一覧  ######################################

 
def Measurement(ser,t3,t4):                               ## Arduinoから心拍の計測情報を受け取る関数

    t = int(t4 / t3)                                      ## forのrange範囲 = 操作なしで停止する時間 / センサのサンプリング時間 小数点以下いらないのでintで整数化

    for i in range(0,10000):                              ## while文は無限ループが怖いのでfor文で対応 10000は適当 数値が大きければ何でもok

        try:                                              ## エラーが出た場合，大抵もう一度読み込めば読めるのでcontinueする
            line = ser.readline()                         ## シリアル通信先の出力を読み込む
            V = int(line.strip().decode('utf-8'))         ## 出力をutf-8に変換.末端行削除.strなのでintにして計算できるようにした.
        except:
            continue
 
        if 1 <= V <= 1000 or V >= 3000:                   ## 値が異常のためやり直し
            pass
 
        elif V == 0:                                      ## V = 0は手が離れているため，for文で手が置かれるまでループ
            for n in range(0,t):
                try:
                    line = ser.readline()                 ## シリアル通信先の出力を読み込む
                    V = int(line.strip().decode('utf-8')) ## 出力をutf-8に変換.末端行削除.strなのでintにして計算できるようにした.
                except:
                    continue
 
                if V >= 3000:                             ## 手が置かれた瞬間はライトの光量がセンサに返ってくるまで完全に真っ暗になり，センサ値が3000以上になる(しばらくすると1000～2000くらいになる)
                    break                                 ## 3000以下でbreakする設定すると部屋の光でbreakする可能性あり

                V = 0                                     ## 一定時間計測できなかった場合はV = 0 でセンサ値を返す
                    
            if V ==0:
                break
            else:
                pass
        else:
            break
 
    return V                                              ## センサの計測値を返す(センサ値は電圧)

def threshold_setting(ser,t1,t3,t4,a,b):                  ## 閾値設定関数

    print("閾値設定開始")
 
    V2 = Measurement(ser,t3,t4)                           ## センサ値読み取り

    t = int(t1 / t3)                                      ## for文の繰り返し回数 = 閾値設定時間 / センサのサンプリング時間 intで整数化

    stack_V = [V2]                                        ## センサ値の平均を求めるのに使用   (V:Voltage)
    stack_D = [0]                                         ## センサ値の差．センサ値の閾値を決める計算に使用 (D:Difference)
 
    for i in range(0,t):                                  ## 心拍のサンプリング時間×t回 = 設定時間

        V1 = Measurement(ser,t3,t4)                       ## センサ値読み取り
        D = V2 - V1                                       ## 一個前のセンサ値 - 現在のセンサ値        

        if D <= 0:                                        ## 差がマイナスならプラスに直す
            D = D * -1                            
 
        stack_D = np.hstack([stack_D,D])                  ## 差をリストに格納
        stack_V = np.hstack([stack_V,V1])                 ## 測定値をリストに格納
        
        V2 = V1                                           ## 次のループ用に一個前のセンサ値として現在のセンサ値を登録

    start_V,upper_limit,down_limit = threshold_calculating(stack_V,stack_D,a,b)  ## 閾値計算関数

    V_list           = [V2]                               ## Voltage計測結果のリスト
    upper_limit_list = [upper_limit]                      ## 上限閾値のリスト
    down_limit_list  = [down_limit]                       ## 下限閾値のリスト
    max_V            = max(stack_V)                       ## 閾値設定中の最大電圧

    print("閾値設定終了")
    print("______________________________")
    print("                              ")
    print("現在の心拍センサ電圧値: " + str(V2))
    print("閾値設定中の最大センサ電圧値" + str(max_V))
    print("電圧平均: " + str(start_V))
    print("閾値上限: " + str(upper_limit))
    print("閾値下限: " + str(down_limit))
    print("______________________________")
    print("                              ")
    print("で測定開始します")
    print("                              ")
 
    return V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list,max_V 
 
    ## 返し値:左から一個前のセンサ値，Vの平均，差の平均，上限閾値，下限閾値，センサ値のリスト，上限閾値のリスト，下限閾値のリスト

def threshold_calculating(stack_V,stack_D,a,b):           ## 閾値計算関数

    move_ave_V       = np.convolve(stack_V,np.ones(a)/a,'valid')   ## センサ値の移動平均を算出
    start_V          = move_ave_V[-1]                              ## センサ値の移動平均のリストの一番最新の値を出力 
    move_ave_D       = np.convolve(stack_D,np.ones(b)/b,'valid')   ## 差の移動平均を算出
    ave_D            = move_ave_D[-1]                              ## 差の移動平均のリストの一番最新の値を出力 
    upper_limit      = start_V + (ave_D * (2/3))                   ## 閾値の上限を設定
    down_limit       = start_V + (ave_D * (1/3))                   ## 閾値の下限を設定

    return start_V,upper_limit,down_limit
 
## ↓開始前に1分計測する関数(BPMを求めるのに，1分間の脈動回数が必要なため)

def start_setting(ser,t3,t4,V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list):
 
    n = int(60 / t3)                                      ## for文の繰り返し回数 = 60秒 / センサ値のサンプリング時間 intで整数化
    c = 0                                                 ## 1秒の脈拍回数の保存用
    e = 0                                                 ## 1秒を図るのに使用
    t = 0                                                 ## 計測時間の保存用
    g = 0                                                 ## 閾値の切替に使用
    
    stack_R = [0]                                         ## 1秒の脈拍回数の保存用リスト

    print("1分後に測定開始します。")
 
    for i in range(0,n):                                  
 
        V1 = Measurement(ser,t3,t4)                       ## センサ値を読み取り
        
        if V1 == 0:                                       ## 一定時間計測できなかった場合，V = 0が出力されるので停止
            
            print("一定時間心拍が読み取れなかったため，中止しました")
            break
 
        t = t + t3                                        ## サンプリング時間を足していき、総計測時間を求める (センサ読み取り関数のエラーしたときなどのcontinueした時間を足すとより正確な計測時間になる)
        e = e + t3                                        ## サンプリング時間を足していき、1秒を求める      　(同上)
        
        Extreme_Value1 = V1 - V2                          ## センサ値の傾きを求める

        if g == 0:                                        ## 閾値が下限の時

            stack_V = np.hstack([stack_V,float(V1)])      ## 脈動していない時のセンサ値を保存するリストを更新(閾値を求めるのに使用)

            if V2 <= down_limit:                          ## 下限の閾値を下回ったら上限の閾値に移行
                g = 1
 
        elif g == 1:                                      ## 閾値が上限の時
            
            if V2 >= upper_limit and Extreme_Value2 >= 0 and Extreme_Value1 <= 0: ## 脈拍したとき
                
                c = c + 1                                 ## 脈拍回数を1増やす

                print(t)
 
                D = V2 - start_V                          ## 脈拍してない時としてる時のセンサ値の差を求める(閾値の設定に使用)
    
                if D <= 0:                                ## 差が負の場合，正に直す
                    D = D * -1
 
                stack_D = np.hstack([stack_D,float(D)])   ## センサ値の差を保存するリストを更新(閾値を求めるのに使用)
                
                g = 0                                     ## 閾値を下限に移行

            elif V1 <= down_limit:                        ## センサ値が下限以下の時
                
                stack_V = np.hstack([stack_V,float(V1)])  ## 脈動していない時のセンサ値を保存するリストを更新(閾値を求めるのに使用)

        if e >= 1:                                        ## 1秒経過したら

            stack_R = np.hstack([stack_R,c])              ## 脈拍回数を保存

            c = 0                                         ## 脈拍回数をリセット
            e = 0                                         ## 1秒間の計測をリセット
    
        start_V,upper_limit,down_limit = threshold_calculating(stack_V,stack_D,a,b) ## 閾値を計算
    
        V2 = V1                                           ## 1回前のセンサ値を保存しておく(センサ値の差を求めるのに使用)
        Extreme_Value2 = Extreme_Value1                   ## 1回前のセンサ値の傾きを保存(脈拍を計測するのに使用)

        if 15.01 >= t >= 15.00:                           ## 15秒経過
            print("15秒経過")
        elif 30.01 >= t >= 30.00:                         ## 30秒経過
            print("30秒経過")
        elif 45.01 >= t >= 45:                            ## 45秒経過
            print("45秒経過")
 
    print("計測開始")
    print("　　　　")
 
    return V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list,stack_R
 
def BPM(ser,t1,t2,t3,t4,a,b):                             ## BPMの計測関数

    t  = 0                                                ## 計測時間の保存用
    c  = int(t2 / t3)                                     ## for文の繰り返し回数 = 計測時間[秒] / センサ値のサンプリング時間 intで整数化
    n  = 0                                                ## 閾値の切替に使用
    f = 0                                                 ## 1秒を図るのに使用
    g = 0                                                 ## 1秒の脈拍回数の保存用
    Extreme_Value2 = 0                                    ## 1回前のセンサ値の傾きを保存(脈拍を計測するのに使用)
    
    BPMdata  = ['time','BPM']                             ## BPMの保存用リスト
    original = ['time','Voltage','V_average','upper_limit','down_limit']  ## 原信号保存用リスト

    ## ↓閾値設定関数
    V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list,max_V = threshold_setting(ser,t1,t3,t4,a,b)
 
    ## ↓開始前に1分計測する関数(BPMを求めるのに，1分間の脈動回数が必要なため)
    V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list,stack_R = start_setting(ser,t3,t4,V2,stack_V,stack_D,start_V,upper_limit,down_limit,V_list,upper_limit_list,down_limit_list)
 
    for i in range(0,c):
 
        V1 = Measurement(ser,t3,t4)                       ## センサ値を読み取り

        if V1 == 0:                                       ## 一定時間計測できなかった場合，V = 0が出力されるので停止
            print("一定時間心拍が読み取れなかったため，中止しました")
            break
 
        t = t + t3                                        ## サンプリング時間を足していき、総計測時間を求める (センサ読み取り関数のエラーしたときなどのcontinueした時間を足すとより正確な計測時間になる)
        f = f + t3                                        ## サンプリング時間を足していき、1秒を求める      　(同上)
        
        Extreme_Value1 = V1 - V2                          ## センサ値の傾きを求める

        if n == 0:                                        ## 閾値が下限の時

            stack_V = np.hstack([stack_V,float(V1)])      ## 脈動していない時のセンサ値を保存するリストを更新(閾値を求めるのに使用)

            if V2 <= down_limit:                          ## 下限の閾値を下回ったら上限の閾値に移行
                n = 1
 
        elif n == 1:                                      ## 閾値が上限の時
            
            if V2 >= upper_limit and Extreme_Value2 >= 0 and Extreme_Value1 <= 0: ## 脈拍したとき
                
                g = g  + 1                                ## 脈拍回数を1増やす

                D = V2 - start_V                          ## 脈拍してない時としてる時のセンサ値の差を求める(閾値の設定に使用)
    
                if D <= 0:                                ## 差が負の場合，正に直す
                    D = D * -1
 
                stack_D = np.hstack([stack_D,float(D)])   ## センサ値の差を保存するリストを更新(閾値を求めるのに使用)
                
                n = 0                                     ## 閾値を下限に移行

            elif V1 <= down_limit:                        ## センサ値が下限以下の時
                
                stack_V = np.hstack([stack_V,float(V1)])  ## 脈動していない時のセンサ値を保存するリストを更新(閾値を求めるのに使用)

        if f >= 1:                                        ## 1秒経過したら

            stack_R = np.hstack([stack_R,g])              ## 脈拍回数を保存

            BPM = BPM_calculating(stack_R)                ## BPMを計算
            
            print(BPM,t,V1)
 
            data2 = [t,BPM]                               ## BPMをリスト化
            BPMdata = np.vstack([BPMdata,data2])          ## リスト化したBPMを保存

            f = 0                                         ## 脈拍回数をリセット
            g = 0                                         ## 1秒間の計測をリセット
    
        start_V,upper_limit,down_limit = threshold_calculating(stack_V,stack_D,a,b) ## 閾値を計算

        data1 = [t,V1,start_V,upper_limit,down_limit]     ## 原信号をリスト化
        original = np.vstack([original,data1])            ## リスト化した原信号を保存
    
        V2 = V1                                           ## 1回前のセンサ値を保存しておく(センサ値の差を求めるのに使用)
        Extreme_Value2 = Extreme_Value1                   ## 1回前のセンサ値の傾きを保存(脈拍を計測するのに使用)

    ser.close()                                           ## シリアルポートを閉じる

    return original,BPMdata                               ## 原信号とBPMのリストを出力

def BPM_calculating(stack_R):
 
    b = len(stack_R)                                      ## 脈拍回数のリストの要素数を出力
    c = 0                                                 ## 脈拍回数を出力した際に使用
    d = 0                                                 ## 脈拍回数の総和を求めるのに使用

    ## 脈拍回数のリストに最初0が入っているため，for文の回数はリストの要素数-1しなくていい
    
    if b >= 60:                                           ## 要素が60以上なら
        a = b - 60                                        ## 要素から60引いた数を求める(for文で最新60秒の脈拍回数の総和を求めたいため)               
    else:
        a = 0                                             ## 60以下だと引くとマイナスになってしまうため,0 (基本ないはず)
    
    for i in range(a,b):                                  ## 最新60秒間の脈拍回数からBPMを求める

        c = stack_R[i]                                    ## 1秒間の脈拍回数を出力
        d = d + c                                         ## 60回1秒間の脈拍回数を足せば自動的にBPMになる

    return d                                              ## BPMを出力

def savedata(original,BPMdata):                           ## 結果保存用関数

    with open('Original.csv','a') as f:                   ## csv製作 csvを開く
        writer = csv.writer(f,lineterminator='\n')        ## /nで改行
        writer.writerows(original)                        ## 複数行の配列の出力　原信号を保存

    f.close()                                             ## csvを閉じる

    with open('BPM.csv','a') as f:                        ## csv製作 csvを開く
        writer = csv.writer(f,lineterminator='\n')        ## /nで改行
        writer.writerows(BPMdata)                         ## 複数行の配列の出力　BPMを保存

    f.close()                                             ## csvを閉じる

    print("end")
 
 
def main():                                               ## メインループ関数
    original,BPMdata = BPM(ser,t1,t2,t3,t4,a,b)           ## 計測関数
    savedata(original,BPMdata)                            ## 結果保存用関数

######################################  実行処理  ######################################

 
if __name__ == "__main__":                                ## 他のプログラムにimportしたときに勝手に実行しないためのおまじない
    main()                                                ## メインループ関数

 

心電図で胸部から測定した心拍数と上記のプログラムと心拍センサを用いて手首から測定した心拍数の比較実験を行った．結果、心電図の心拍数71[回]に対し、上記のセンサ及びプログラムは70[回]で誤差率0.7[%]の精度であった．心電図が胸部で測定を行っているのに対し、心拍センサは手首で測定を行っており、心臓の脈動を検知するタイミングがずれていると考えられることから、±1[回]の心拍数の誤差は許容範囲内に収まっていると考えられる．

原信号のCSVファイルから心拍数を求める

原信号➡心拍数[python]

import csv
import numpy as np
 
csv_file = open("Original.csv", "r", encoding="ms932", errors="", newline="" )
#リスト形式
f = csv.reader(csv_file, delimiter=",", doublequote=True, lineterminator="\r\n", quotechar='"', skipinitialspace=True)
 
data = [ v for v in f]
 
t = 0
t1 = 0
Extreme_Value2 = 0
m = 1
n = 1
a = 0
 
BPMdata = ['time','Voltage','BPM']
x =  [0]
x2 = [0]
y2 = [0]
sa_stack = [0]
 
print("閾値設定")
 
c =data[1]
V2 = int(float(c[1]))
start_value=[V2]
 
for i in range(2,201):
    c = data[i]
 
    V1 = int(float(c[1]))
    sa = V2 - V1
    if sa <= 0:
        sa = sa * -1
 
    if V1 == 0:
        pass
    else:
        start_value = np.hstack([start_value,V1])
        sa_stack = np.hstack([sa_stack,sa])
        V2 = V1
 
 
print("start")
 
V_start = np.average(start_value)
V_start = int(V_start)
sa_ave = np.average(sa_stack)
alpha = V_start + (sa_ave * (2/3))
alpha2 = V_start + (sa_ave * (1/3))
y =  [V1]
alpha_line =  [alpha]
alpha_line2 =  [alpha2]
 
print(V_start,sa_ave,alpha,alpha2)
 
for i in range(201,8601):
    c = data[i]
 
    V1 = int(float(c[1]))
    t = t + 0.05
    Extreme_Value1 = V1 - V2
 
    if V2 >= alpha and Extreme_Value2 >= 0 and Extreme_Value1 <= 0 and a==0:
        BPM = 60/(t - t1)
        print(BPM,t)
        data2 = [t,V2,BPM]
        BPMdata = np.vstack([BPMdata,data2])
        t1 = t
 
        sa = V2 - V_start
 
        if sa <= 0:
            sa = sa * -1
 
        sa_stack = np.hstack([sa_stack,float(sa)])
        move_ave_sa = np.convolve(sa_stack,np.ones(5)/5.0,'valid')
        sa_ave = move_ave_sa[-1]
 
        x2 = np.vstack([x2,t])
        y2 = np.vstack([y2,BPM])
        a = 1
 
    if V2 <=alpha2 and a == 1:
        a=0
 
    if a == 0:
        start_value = np.hstack([start_value,float(V1)])
        move_ave = np.convolve(start_value,np.ones(50)/50.0,'valid')
        V_start = move_ave[-1]
 
    alpha = V_start + (sa_ave * (2/3))
    alpha2 = V_start + (sa_ave * (1/3))
    x = np.vstack([x,t])
    y = np.vstack([y,V1])
    alpha_line = np.vstack([alpha_line,alpha])
    alpha_line2 = np.vstack([alpha_line2,alpha2])
 
    V2 = V1
    Extreme_Value2 = Extreme_Value1
 
    n = n + 1
 
 
print("end")
 
with open('BPM.csv','a') as f:
    writer = csv.writer(f,lineterminator='\n')  #/nで改行
    writer.writerows(BPMdata)                  #複数行の配列の出力

f.close()
 
 
 

ローレンツプロットの面積

ローレンツプロット

　今回はArudino esp32を用いてRRIからローレンツプロットの面積を求める。下記に書くプログラムは
ローレンツプロットの面積を求めるプログラムであってRRIを求めてないのでこのプログラムでは動かない。
また行っているのは大まかに以下の5つの手順を行っている。

①RRIデータをスタック
②RRIデータからローレンツプロットをプロット
③ローレンツプロットのプロット点を-45度回転させる
④回転後のプロット点のx・yの標準偏差を求める
⑤標準偏差から楕円の公式を用いて面積を求める

volatile float  angle = -1*(3.141592653589793/4);      //ローレンツプロット点の回転角度
volatile int RRI_data[650];                            //ここでローレンツプロットのプロット点を保存する物を650個以下とするもしそれを超えるとエラーが出る
volatile float Area_ave = 0;                           //面積の平均
volatile float Area_sum = 0;                           //面積の合計
volatile float Area = 0;                               //面積
volatile int zz  = 1;                                  //ローレンツプロットの面積を求める範囲の変更
volatile int x  = 0;                                   //プロット点の個数を確認
volatile float Average_xdata = 0;                      //x座標の平均値
volatile float Average_ydata = 0;                      //y座標の平均値
volatile float RRI_std_xdata = 0;                      //x軸の標準偏差
volatile float RRI_std_ydata = 0;                      //y軸の標準偏差
void setup() {
 }
void loop() {
 
      //略//
 
        //ローレンツプロットの面積を求めるためにRRIをスタック
        if(TotalTime>=60000 && 500<=RRI && RRI<=1000){      //始めの1分間はRRIの閾値設定に使うのでその後にRRIが500～1000の間の物を1時保存
            RRI_data[x] = RRI;                              //RRIをRRI_dataにスタック
            x = x+1;                                        //個数を確認
        }
 
	//略//
 
    //ローレンツプロットの面積を求める
    if(TotalTime>=60000+(300000*zz)){              //最初の1分間(60000)を除いた後5分間(300000ms)おきにArea面積を求める 
      zz = zz+1;                                   //1回計算するごとに上記の条件を変えるよう
      int RRI_plot_xdata[x-1];                     //ローレンツプロットのプロット点のx座標（x-1は最後に1回余計にx+1しているため）
      int RRI_plot_ydata[x-1];                     //ローレンツプロットのプロット点のx座標（x-1は最後に1回余計にx+1しているため）
      for(int z = 1; z <= (x-1); z++){
        RRI_plot_xdata[z] = RRI_data[z-1];         //xのデータは0から始まっておりzは1からなので最初のデータから所得する
        RRI_plot_ydata[z] = RRI_data[z];           //ローレンツプロットはn+1した値をy軸データとするため1点ずらす
        }
 
      float rotation_RRI_xdata[x-1];               //回転した後のx座標
      float rotation_RRI_ydata[x-1];               //回転した後のy座標
      float RRI_sum_xdata = 0;                     //1回1回リセットする
      float RRI_sum_ydata = 0;                     //1回1回リセットする
 
      for(int z = 1; z <= (x-1); z++){             //スタックしてあるものを読み込むために繰り返す
        rotation_RRI_xdata[z] = RRI_plot_xdata[z]*cos(angle)-1*RRI_plot_ydata[z]*sin(angle);        //x座標の値を-45度回転させプロット点を回転
        RRI_sum_xdata = RRI_sum_xdata+rotation_RRI_xdata[z];                                        //Averageを求めるために合計を求める
        rotation_RRI_ydata[z] = RRI_plot_xdata[z]*sin(angle)+1*RRI_plot_ydata[z]*cos(angle);        //y座標の値を-45度回転させプロット点を回転
        RRI_sum_ydata = RRI_sum_ydata+rotation_RRI_ydata[z];                                        //Averageを求めるために合計を求める
        }
      Average_xdata = RRI_sum_xdata/(x-1);        //標準偏差を求めるために回転後のxdataの平均値を求める
      Average_ydata = RRI_sum_ydata/(x-1);        //標準偏差を求めるために回転後のydataの平均値を求める
 
      float s_xdata = 0;                          //1回1回リセットする
      float s_ydata = 0;                          //1回1回リセットする
 
      for(int z = 1; z <= (x-1); z++){
        s_xdata = s_xdata+sq((rotation_RRI_xdata[z]-Average_xdata));          //それぞれの回転後のxdataから平均値を引いた値を2乗
        s_ydata = s_ydata+sq((rotation_RRI_ydata[z]-Average_ydata));          //それぞれの回転後のydataから平均値を引いた値を2乗
        }
      RRI_std_xdata = sqrt(s_xdata/(x-1));                                    //xの標準偏差を求める
      RRI_std_ydata = sqrt(s_ydata/(x-1));                                    //yの標準偏差を求める
 
      Area = (RRI_std_xdata*RRI_std_ydata*3.141592653589793)/4;               //ローレンツプロットの楕円の面積を求める
 
      x = 0;          //個数をリセット
    Serial.println(Area); //面積の値を出力
 

注意・・・//略//の部分でセンサーからの信号を得て経過時間TotalTimeと心拍間隔RRIを算出したりしている。

ローレンツプロットの面積からストレスが評価できるか検証するため、唾液アミラーゼと比較実験を行った．結果として、唾液アミラーゼとローレンツプロットの面積のストレス評価が同じような結果になったことを確認した．

活動量計による心拍数の測定(Fitbit Charge 2)

fitbit APIの心拍数データ取得方法 [Python]

プログラム製作後、更新予定

研究の進捗状況

2019年度

• ロボット製作拠点であるハイテクリサーチセンターの整備
• 使用する心拍センサと心電図の比較実験を行い、精度確認の実施
• ロボット試作1号機の製作
• 心拍数の測定プログラムの製作
• RRI間隔のローレンツプロット面積を用いたストレス評価プログラムの製作
• 唾液アミラーゼとローレンツプロット面積によるストレス値の比較実験

今後の展望

●ローレンツプロットによるストレス評価実験のN増し

●ロボットの動きや色の変化を用いてストレス値を評価する
➡現在はロボットのしっぽの動作でストレスがあるかないかを評価しているが，色の変化を用いる場合，青色ならストレスが低く赤色に近づくほどストレス値が高いといったように段階的にストレス値の評価を行えるようにする

●ロボットの動作や外見による癒し効果やストレス値が高い場合にロボットが介入を行うことでストレスの軽減効果があるかを検証する

参考文献

マスコットロボットを活用したリアルタイムによるテクノストレス解消の検証,荒木亮磨,土肥紳一,情報処理学会第80回全国大会,(2018),pp.217-218

• 研究している内容が近いため，読んでおくといいと思います

柔らかいぬいぐるみロボットの動作生成,石川 達也,長谷川 晶一,研究報告 エンタテインメントコンピューティング（EC）,(2011),pp.1-6

• ロボットの動作機構に利用している糸駆動に関する先行研究です

身体性が人工ペットとのふれあいによるセラピー効果に与える影響,林 里奈,加藤 昇平,日本感性工学会論文誌16巻1号,(2017),pp.75-81

• 実物のロボットとバーチャルロボットの違いについて研究が行われています

心拍揺らぎによる精神的ストレス評価法に関する研究，松本 佳昭，森 信彰,三田尻 涼,江 鐘偉,ライフサポート22巻3号,(2010),pp.105-111

• 心拍と精神的ストレスについて詳しく書かれています

平成29年　厚生労働省労働安全衛生調査（実態調査）結果の概況

• 論文ではないですが現代社会においてストレスを感じている社会人が年々増え続けていることがわかります．厚生労働省のHPから見ることができます．

こころの耳　精神障害の労災補償件数の推移と主なできごと

• こちらも論文ではないですが，こころの耳というサイトの数字と絵でわかる職場のメンタルヘルスというページ（http://kokoro.mhlw.go.jp/infographics/02.html）からメンタルヘルスについて考える必要性があることがわかります