この記事は，脳波からヒトの状態推定を可能にするモデルをどうにか作りたくて獅子奮闘している研究アーカイブ（自分用）です．
前回の記事はこちら．
（参考：「最小二乗法によるパラメータ同定奮闘記1」）
最小二乗法は演算処理にかかるコストがめちゃくちゃ低いので，パラメータ同定にかかる時間がかなり早いので使い勝手が良いです．
プログラムも組みやすいし，かなりメジャーなアルゴリズムなので説明コストも低いのが嬉しい．
具体的な内容をさらっと書くと，
脳波の持つ非線形性を，どうにか数値化できるようにするために非線形振動子によるモデリングを行なっています．
扱う非線形振動子は，１自由度系なので直感的に理解し易いです．
さてさて，前回の記事を読んでいただければ話は早いのですが，あまりうまくいかなかったので，パラメタ同定に工夫を入れます．
以下のように修正します．
修正の内容としては，
- Running window内でFFT
- ピークを１つ決める
- 周波数帯域を上図に示すように分けておいて，ピークが属する帯域でIFFT
- IFFTで得た1秒間の時系列波形を用いてパラメタ同定
- モデル右辺の振動数はピークを与える
的な内容です．
ここ周波数帯域の説明を再度付け加えておくと，以下のような10つの帯域に分けてます．
- δ1 (0.5 – 2 Hz)
- δ1 (2 – 4 Hz)
- θ1 (4 – 6 Hz)
- θ2 (6 – 8 Hz)
- α1 (8 – 10 Hz)
- α2 (10 – 13 Hz)
- β1 (13 – 16 Hz)
- β2 (16 – 20 Hz)
- β3 (20 -25 Hz)
- β4 (25 – 30 Hz)
としました．
なので，IFFTを回した結果は，1秒 (512Hz)ごとにこうなります．
上の図は４秒間の一例ですが，本当は被験者ごとに２４０秒間あります．
1秒ごとの支配的な周波数帯域によるIFFT結果です．
今回は，このIFFT結果を用いてパラメタ同定を行ってみました．
被験者aだけの結果を一例として示します．
各パラメタの集中度の特徴が際立ったような気もしてます．
（赤が集中状態，緑が安静状態）
また，被験者ごとの各パラメタ値の相関に関しては，以下のような結果になりました．
このような結果になりました．
相変わらず線形性強度を示すパラメタA,Bは高い相関を示しているのが悔しいです．
まだまだ対応策を考えていきます．
周波数分解能をあげるか，モデルに与えるピークの数を増やすことで再度修正する方向を考えています．
以下はPythonプログラム
# -*- coding: utf-8 -*-
import csv
import scipy.signal as sps
from scipy.optimize import leastsq
import numpy as np
import matplotlib.pyplot as plt
N = 512 # サンプル数
dt = 0.002 # サンプリング間隔
p_max = 30.000000
p_min = 0.500000
p_band = 1.00000
delta_bound, theta_min, theta_bound, alpha_min, alpha_bound, beta_min, beta_bound1, beta_bound2, beta_bound3 = 2.0, 4.0, 6.0, 8.0, 10.0, 13.0, 16.0, 20.0, 25.0
inputData = 'exp.txt'
outputData = 'peak_band/parameters_d_100mass.txt'
ifftData = 'peak_band/ifft_d_100mass.txt'
eeg = list()
i_eeg = list()
t = list()
xd = list()
xdd = list()
w = list()
p = list()
ieeg = list()
def fft_peak(t, eeg):
peak_idx = list()
# データのパラメータ
#t = np.arange(0, N*dt, dt) # 時間軸
freq = np.linspace(0, 1.0/dt, N) # 周波数軸
# 高速フーリエ変換
F = np.fft.fft(eeg)
Power = F.copy()
# 振幅スペクトルを計算
Power = np.abs(Power)
# 振幅を元の信号のスケールに揃える
Power = Power / (N/2) # 交流成分
Power[0] = Power[0] / 2 # 直流成分
F2 = F.copy()
F2 = F2 / (N/2) # 交流成分
F2[0] = F2[0] / 2 # 直流成分
# FFTデータからピークを自動検出
peak_idx = sps.argrelmax(Power, order=1)[0] # ピーク（極大値）のインデックス取得
# ピーク検出感度調整もどき、後半側（ナイキスト超）と閾値より小さい振幅ピークを除外
peak_cut = 0.1 # ピーク閾値
peak_idx = peak_idx[(Power[peak_idx] > peak_cut) & (peak_idx >= p_min) & (peak_idx <= p_max)]
selected_peakPower = dict(zip(freq[peak_idx], Power[peak_idx]))
peak = dict(sorted(selected_peakPower.items(), key=lambda x:x[1], reverse=True))
# p = list(peak.values())
# p = list(map(float, p))
f = list(peak.keys())
f = list(map(float, f))
del f[1:]
print(f[0])
#p_min <= f[0] < delta_bound
if((f[0] >= p_min) & (f[0] < delta_bound)):
p_min2 = p_min
p_max2 = delta_bound
#delta_bound <= f[0] < theta_min
elif((f[0] >= delta_bound) & (f[0] < theta_min)):
p_min2 = delta_bound
p_max2 = theta_min
#theta_min <= f[0] < theta_bound
elif((f[0] >= theta_min) & (f[0] < theta_bound)):
p_min2 = theta_min
p_max2 = theta_bound
#theta_bound <= f[0] < alpha_min
elif((f[0] >= theta_bound)&(f[0] < alpha_min)):
p_min2 = theta_bound
p_max2 = alpha_min
#alpha_min <= f[0] < alpha_bound
elif((f[0] >= alpha_min)&(f[0] < alpha_bound)):
p_min2 = alpha_min
p_max2 = alpha_bound
#alpha_bound <= f[0] < beta_min
elif((f[0] >= alpha_bound)&(f[0] < beta_min)):
p_min2 = alpha_bound
p_max2 = beta_min
#beta_min <= f[0] < beta_bound1
elif((f[0] >= beta_min)&(f[0] < beta_bound1)):
p_min2 = beta_min
p_max2 = beta_bound1
#beta_bound1 <= f[0] < beta_bound2
elif((f[0] >= beta_bound1)&(f[0] < beta_bound2)):
p_min2 = beta_bound1
p_max2 = beta_bound2
#beta_bound2 <= f[0] < beta_bound3
elif((f[0] >= beta_bound2)&(f[0] < beta_bound3)):
p_min2 = beta_bound2
p_max2 = beta_bound3
#beta_bound3 <= f[0] <= p_max
elif((f[0] >= beta_bound3) & (f[0] <= p_max)):
p_min2 = beta_bound3
p_max2 = p_max
# p_max2 = f[0] + p_band
# p_min2 = f[0] - p_band
# del p[2:]
# w = 2.0 * np.pi * f
w = f
w = [i*2.0*np.pi for i in w]
#ローパスフィルタ
F2[(freq > p_max2)] = 0.0
#ハイパスフィルタ
F2[(freq < p_min2)] = 0.0
#逆高速フーリエ変換
i_eeg = np.fft.ifft(F2)
i_eeg = np.real(i_eeg*N)
# グラフ表示
plt.figure()
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 17
plt.subplot(311)
plt.plot(t[:], eeg[:], label='f(n)')
plt.xlabel("Time", fontsize=20)
plt.ylabel("raw_Signal", fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)
plt.subplot(312)
plt.xlim([0,30])
plt.plot(freq[:250], F2[:250], label='|F(k)|')
plt.xlabel('Frequency', fontsize=20)
plt.ylabel('Amplitude(afterBPF)', fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)
plt.subplot(313)
plt.plot(t[:], i_eeg[:], label='|F(k)|')
plt.xlabel("Time", fontsize=20)
plt.ylabel("ifft_Signal", fontsize=20)
plt.grid()
leg = plt.legend(loc=1, fontsize=25)
leg.get_frame().set_alpha(1)
plt.show()
return(w, p, i_eeg)
return(w, i_eeg)
#def p_cal(p):
# p_sum = p[0]+p[1]+p[2]
# p[0] = p[0]/p_sum#*p_sum
# p[1] = p[1]/p_sum#*p_sum
# p[2] = p[2]/p_sum#*p_sum
# return p
def cal(t, x):
j = 0
x = list(x)
xd.append(float((x[j+1] - x[j])/dt)) #xd[0]
j+=1
for j in range(N-1):
xd.append(float((x[j+1] - x[j])/dt))#xd[1]
xdd.append(float((xd[j] - xd[j-1])/dt))#xdd[0]
t.pop()
x.pop()
xd.pop()
return (t,x)
def func3(param, t, x, xd, xdd):
residual = xdd + (param[0]*xd + param[1]*x + param[2]*x*x*x - param[3]*np.cos(w[0]*t)- param[4]*np.sin(w[0]*t))
return residual
def lsm(t, x, xd, xdd):
para = [0., 0., 0., 0., 0.]
#para = [0, 0, 0, 0]
t = np.array(t)
x = np.array(x)
xd = np.array(xd)
xdd = np.array(xdd)
params = leastsq(func3, para, args=(t, x, xd, xdd))
return params
with open(outputData, mode='w') as OutputFile:
with open(inputData, 'r') as InputFile:
reader = csv.reader(InputFile, delimiter='\t')
i = 0
z = 0
for row in reader:
t.append(float(row[0]))
eeg.append(float(row[4]))
# print("{:d}".format(i))
# print("{:f}\t{:f}".format(float(row[0]), float(row[1])))
if(((i+1)%N==0) & (i>1)):
print("fft_start")
w, p, i_eeg = fft_peak(t, eeg)
ieeg.extend(i_eeg)
print(w)
print(p)
#p = p_cal(p)
t, i_eeg = cal(t, i_eeg)
Result = lsm(t, i_eeg, xd, xdd)
print(Result)
OutputFile.write("{:.3f}\t{:.6f}\t{:.6f}\t{:.6f}\t{:.6f}\t{:.6f}\t{:.6f}\n".format(0.510+z*dt*N,Result[0][0],Result[0][1],Result[0][2],Result[0][3],Result[0][4],w[0]))
print("xdd + {:.6f}xd + {:.6f}x + {:.6f}x*x*x - {:.6f}cos{:.6f}t - {:.6f}sin{:.6f}t\n".format(Result[0][0],Result[0][1],Result[0][2],Result[0][3],w[0],Result[0][4],w[0]))
z+=1
eeg = list()
i_eeg = list()
t = list()
xd = list()
xdd = list()
# if(i > 1538):
# break
i+=1
ieeg = map(str,ieeg)
textdata = '\n'.join(ieeg)
with open(ifftData, 'wt') as f:
f.write(textdata)
# next(reader)
