FFT - yana_lab

FFT(Fast Fourier Transform)

FFT(日:高速フーリエ変換)とは波形に含まれる周波数成分を数値化する周波数解析。
下の例は、左の0.15Hzと0.3Hzの合成波のデータに対してFFTの処理を行い右図のように各周波数帯を数値化をした。0.15と0.3Hzのみの合成波を解析したため他の周波数帯に成分が存在しない結果が表れた。このサンプルプログラムは後の「ストレス指標LF/HFの算出」の項目に載せており、そこから同様の結果が得られる。

詳細な計算式については以下のサイトを参照
https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79

ストレスの評価

　心拍間隔(RRI)データを経過時間ごとにまとめた心拍数変動時系列データに含まれる周波数成分は交感神経や副交感神経の影響を受ける。そのためデータに対して周波数解析を行い、後述する周波数帯の成分を解析することでストレス状態を数値化する。

FFTの前処理

　FFTを心拍変動時系列データに掛ける前にはリサンプリングとデトレンドの処理を行う必要がある。
　リサンプリングとはサンプリングしたデータを別の標本点系列でサンプリングを行う処理。周波数解析をデータに行うためには解析されるデータのサンプリングは一定の間隔である必要がある。拍動の間隔は不規則であるため、心拍データが一定の間隔でサンプリングされないため行う。
　デトレンドとは波形に含まれるトレンド成分を除去する処理。トレンド成分とは波形に含まれる体系的な増加傾向や減少傾向のこと。このトレンド成分が残ったまま解析を行ってしまうと余計なエネルギー量が含まれてしまうため行う。
　下の図は実際に計測した20代の男性のデータであり、左から順番に心拍変動時系列データ、線形補間データ、デトレンドデータである。心拍変動時系列データに対して線形補間を掛けている様子やデトレンドによりトレンド成分が除去されたことによる微細な変化が確認できる。

プログラム

　下記のプログラムはExcelのRRIデータに対してリサンプリングとデトレンドの処理を行い、データを保存するプログラムである。
　リサンプリングはサンプリングしたデータ間の欠損した部分を線形補間を行うことで補った後、一定の間隔でデータを抽出している。そのため、補完関数fを用いて線形補間を行った後リサンプリング波形を生成する。fftを行う場合は解析されるデータ数が2のn乗個でなければならないため、numの値を128としている。サンプリング周期を1秒にするとサンプル数は128個あるため、128秒ごとに解析を行うプログラムとなる。また、デトレンドはscipyのデトレンド関数を使用して処理を行った。

注意
　実際にプログラムを動かす際は参照するファイル名が正しいことやExcelがcsvファイルで保存されているか確認する。また、リサンプリング波形を生成する際のtsとteのプログラムはデータによって修正する必要がある。tsは解析するデータの最初の計測開始時間、teは解析するデータの最後の計測時間を指定しており、心拍を測定しながら連続で解析を行う場合はtsとteの値が解析ごとに再指定されるようにプログラムする必要がある。

プログラム

#RRIデータに対してリサンプリングとデトレンドの処理を行うプログラム

from scipy import interpolate
from matplotlib import pyplot as plt
import numpy as np
from scipy import signal
 
#Excelファイルの参照 
data = np.loadtxt("RRI_resultA_1113am.csv",delimiter=',')
t = data[:,0]#Excelの1列目全ての要素
y = data[:,1]#Excelの2列目全ての要素

# 補間関数fを作成
#線形補間によりデータ間の欠損した部分を補う
f = interpolate.interp1d(t, y, kind='linear')
 
# 補間した結果からリサンプリング波形を生成
ts = 1    
te = 129
num = 128
t_resample = np.linspace(ts, te, num)
#「ts」から「te」まで「num」個の要素数で等間隔の配列を生成。このnumはサンプリング間隔が0.1sになるように調整している。
y_resample = f(t_resample)#補間関数fから直前に生成した等間隔の配列でデータを抽出することでリサンプリングを行う
print(y_resample)
 
yd = signal.detrend(y_resample)#デトレンド

savedata = np.vstack([t_resample, yd])
#print(savedata)
np.savetxt("re_test3.csv",savedata.T,delimiter=',')#データの保存
#print(savedata.T)
plt.subplot(2,1,1)
plt.plot(t,y,'b')
plt.subplot(2,1,2)
plt.plot(t_resample,yd,'r')
plt.show()
 

ストレス指標LF/HFの算出

　心拍変動時系列データに含まれる低周波成分LF(0.05~0.15Hz)は交感神経と副交感神経の機能を反映し、高周波成分HF(0.15~0.4Hz)は副交感神経の機能を反映する。人間はリラックス時には副交感神経が活発化しLFが増大する。反対にストレスを感じると交感神経が活発化しLFが増大、HFが減少する。そのため、LF/HFをストレスの指標とした。ストレスが高まるとリラックス時に比べてLF/HFは値が高くなる。しかし、LF/HFの値は個人差があり明確にストレス状態と判断する閾値がないため、判断を行うためには実際に計測を行い個人の傾向を調べる必要がある。
　下の図は左が心拍変動時系列データにリサンプリングとデトレンドを行ったデータ、右がそれにFFTの処理を行ったデータである。また、左の図は先程も挙げた20代の男性のデータと同様のものである。この実験データからLF/HFは約1.49と算出された。実験者はリラックス状態およそで0.5~0.8、平常時で0.8~1.1付近の値を算出していたため、若干のストレスや疲れがみられることが分かった。

プログラム

　下のプログラムは上記のリサンプリングとデトレンドの処理を行ったExcelファイルに対して行うFFTプログラムである。このプログラム単体でもコメントアウトを消去することで確認用のサンプリングに対してFFTを掛けることができる。サンプル波形のFFTの様子はこのページのトップ部分の図から確認できる。

プログラム

#周波数解析FFTを行うプログラム

####注意##########################################################
#このプログラムはリサンプリングとデトレンドを行ったデータに行うこと
###################################################################

import numpy as np
import matplotlib.pyplot as plt
 
data = np.loadtxt("re_test3.csv",delimiter=',')
t = data[:,0]
x = data[:,1]
 
N = 128            # サンプル数
dt = 1          # サンプリング周期 [s]

##確認用#########################################################
#f1, f2 = 0.15, 0.3    # 周波数 [Hz]
#t = np.arange(0, N*dt, dt) # 時間 [s]
#x = 1.5*np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + 3 # 信号
#################################################################

F = np.fft.fft(x) # 変換結果
freq = np.fft.fftfreq(N, d=dt) # 周波数
Amp = np.abs(F/(N/2)) # 振幅

#####LF_HF算出用########################
#print(freq)
LF=HF=0
z=0
while freq[z] < 0.05:
    z=z+1
while freq[z] < 0.15:
    LF=LF+Amp[z]
    z=z+1
while freq[z] < 0.4:
    HF=HF+Amp[z]
    z=z+1
LF_HF=LF/HF
print(LF_HF)
############################

plt.plot(freq[1:int(N/2)], Amp[1:int(N/2)])
plt.xlabel("Freqency [Hz]")
plt.ylabel("Amplitude")
plt.xlim([0,0.5])##横軸の設定

plt.show()