●Spectral Analysis 2
セボフルラン全身麻酔中の脳波のスペクトル解析
・実際にセボフルラン全身麻酔中にEEG Analyzerを用いて、BISモニタから取得したデジタル脳波データファイルeeg_bis.tsvを用いて、Pythonによるスペクトル解析を実践する。窓関数の影響、マルチテーパー法の適用、そしてスペクトログラムへ発展させる。
データファイル：eeg_bis_10_20_2.tsv
Jupyter Notebook: spectral_analysis_1.ipynb
深い麻酔の状態の脳波分析
|
1 2 3 4 5 6 7 8 |
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.fftpack import fft from scipy import signal as sig #SciPyパッケージからsignal関数インポート from IPython.display import set_matplotlib_formats set_matplotlib_formats('svg') |
・EEG Analyzerを用いて、BISモニタから取得したデジタル脳波データファイルeeg_bis_10_20_2.tsvから脳波データを、全部PandasのDataFrameであるdf_tsvに取り込ませる。
|
1 |
df_tsv = pd.read_table('eeg_bis_10_20_2.tsv') |
・取り込んだデータの先頭５件を表示
|
1 |
df_tsv.head(5) |
Ch Time ch[0] ch[1] ch[2] ch[3] ch[4] ch[5] ch[6] ch[7] ch[8] ch[9] ch[10] ch[11] ch[12] ch[13] ch[14] ch[15]
0 ch1: 09:18:52 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.0
1 ch1: 09:18:52 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.0
2 ch1: 09:18:52 27.35 26.90 27.00 23.05 16.15 9.7 6.35 8.30 13.30 24.65 37.40 41.55 33.15 15.40 -9.80 -21.0
3 ch1: 09:18:52 -5.30 22.50 40.80 42.45 30.45 18.0 11.60 5.65 -0.95 -3.85 0.25 12.20 28.35 40.05 43.20 31.7
4 ch1: 09:18:52 12.60 5.15 9.65 13.10 14.80 17.7 21.80 27.55 29.25 23.85 15.65 9.30 3.80 5.25 17.85 31.4
・取り込んだデータのラスト５件を表示
|
1 |
df_tsv.tail(5) |
Ch Time ch[0] ch[1] ch[2] ch[3] ch[4] ch[5] ch[6] ch[7] ch[8] ch[9] ch[10] ch[11] ch[12] ch[13] ch[14] ch[15]
109607 ch1: 13:13:31 5.80 5.90 5.75 5.65 5.80 5.90 6.15 6.20 5.85 5.55 5.60 5.70 5.50 5.55 5.65 5.70
109608 ch1: 13:13:31 5.80 5.85 5.85 6.00 5.75 5.65 5.80 5.65 5.70 6.15 6.10 5.85 5.75 5.75 5.90 6.10
109609 ch1: 13:13:31 5.70 5.60 5.65 5.65 5.60 5.60 5.65 5.85 5.65 5.80 6.00 5.70 5.50 5.65 5.55 5.75
109610 ch1: 13:13:31 5.90 5.85 5.90 5.95 6.00 6.00 5.90 6.00 5.95 5.70 5.85 6.10 5.75 5.70 5.95 6.10
109611 ch1: 13:13:31 5.85 5.95 6.15 6.05 6.05 6.20 6.00 5.85 5.85 5.80 5.80 6.05 5.80 5.75 5.85 5.85
・PandasのDataFrameであるdf_eegを新たに作成し、コラム名”Time”と”eeg”を設定する。
|
1 2 |
cols = ['Time', 'eeg'] df_eeg = pd.DataFrame(data=None, index=[], columns=cols) |
・df_tsv（109611行）のうち、深いセボフルラン麻酔状態にある部分の電圧?Vデータを1024行をdf_eegに写す。
|
1 2 3 4 5 6 7 8 9 |
#for row in range(len(df_tsv)): #かなりの時間がかかります。 for row in range(80000, 81024): #1024行のデータだけ取り込む。8行/秒なので、128秒＝2分8秒 data数＝16,384 #1時間分のデータは 8x60x60＝28800行、データ数：460800 #ここのrangeを弄ることで、どこの脳波を取り出せるかが変わる。ただし1時間分くらいでないと時間がかかりすぎる。 t_tmp = df_tsv.iat[row,1] for column in range(2, 18): eeg_tmp = df_tsv.iat[row, column] data_tmp = pd.Series( [t_tmp, eeg_tmp], index=df_eeg.columns ) df_eeg = df_eeg.append(data_tmp, ignore_index=True) |
・df_eegに取り込んだ1024行x16=16383件の?Vデータを確認表示する。128秒分。
|
1 |
df_eeg |
Time eeg
0 12:10:08 21.05
1 12:10:08 17.45
2 12:10:08 14.20
3 12:10:08 7.45
4 12:10:08 -0.25
5 12:10:08 -3.45
6 12:10:08 2.15
7 12:10:08 14.35
8 12:10:08 23.45
9 12:10:08 25.25
10 12:10:08 21.95
11 12:10:08 14.55
12 12:10:08 6.30
13 12:10:08 2.30
14 12:10:08 4.35
15 12:10:08 13.65
16 12:10:08 23.70
17 12:10:08 23.90
18 12:10:08 16.65
19 12:10:08 11.50
20 12:10:08 9.45
21 12:10:08 10.75
22 12:10:08 17.40
23 12:10:08 19.45
24 12:10:08 9.45
25 12:10:08 0.45
26 12:10:08 2.90
27 12:10:08 12.00
28 12:10:08 19.20
29 12:10:08 21.95
… … …
16354 12:12:21 11.25
16355 12:12:21 8.45
16356 12:12:21 0.85
16357 12:12:21 -9.25
16358 12:12:21 -24.45
16359 12:12:21 -41.90
16360 12:12:21 -59.45
16361 12:12:21 -75.10
16362 12:12:21 -78.40
16363 12:12:21 -67.50
16364 12:12:21 -50.40
16365 12:12:21 -34.15
16366 12:12:21 -23.95
16367 12:12:21 -18.80
16368 12:12:21 -14.55
16369 12:12:21 -11.65
16370 12:12:21 -13.90
16371 12:12:21 -30.20
16372 12:12:21 -46.35
16373 12:12:21 -34.35
16374 12:12:21 -2.85
16375 12:12:21 21.90
16376 12:12:21 36.65
16377 12:12:21 39.30
16378 12:12:21 31.70
16379 12:12:21 21.85
16380 12:12:21 11.85
16381 12:12:21 9.00
16382 12:12:21 14.05
16383 12:12:21 10.60
16384 rows × 2 columns
・df_eegに取り込んだデータをcsvファイルとして保存しておく。
|
1 |
df_eeg.to_csv("eeg_linealized_3−R3.csv") |
・df_eegの時間と?Vデータをそれぞれ新たなPython配列Timeとeegに写す。
|
1 2 3 |
#データフレームの脳波データを配置 time = df_eeg['Time'].values eeg = df_eeg['eeg'].values |
・配列eegに取り込んだデータをplotで表示させる。
|
1 2 3 4 |
#脳波データを時系列グラフ表示 fig, ax = plt.subplots(figsize=(18,6)) plt.plot(eeg) plt.savefig('spectral_analysis_2_fig_1.svg') |
データ前処理
・eegの平均値と標準偏差を算定。
|
1 2 |
m=np.mean(eeg) sd=np.std(eeg) |
|
1 |
・null(NA)データ数をカウントする。 |
|
1 |
df_eeg['eeg'].isnull().sum() |
Out: 0
・外れ値の除去：平均値から標準偏差の３倍をはみ出るデータをNoneに変更する。
|
1 2 3 4 5 6 7 |
for i in range(len(eeg)): if df_eeg['eeg'][i] > m+3*sd: df_eeg['eeg'][i] = None for i in range(len(eeg)): if df_eeg['eeg'][i] < m-3*sd: df_eeg['eeg'][i] = None |
・外れ値を除去したあと、再度null(NA)データ数をカウントする。
|
1 |
df_eeg['eeg'].isnull().sum() |
Out: 152
・df_eegのコピーを作成、コピーに対してnull(NA値)を前後の値でスプライン補完する。
|
1 2 |
df_eeg_copy = df_eeg.copy() df_eeg_copy.interpolate(inplace=True) |
・df_eegのコピーのnull(NA)データ数をカウントする。
|
1 |
df_eeg_copy['eeg'].isnull().sum() |
Out: 0
・外れ値を前後のデータでスプライン補完修正したファイルをcsv保存
|
1 |
df_eeg.to_csv("eeg_linealized_3-R3_ip.csv") |
・eegデータを配列eeg_ipへコピーする。
|
1 |
eeg_ip = df_eeg_copy['eeg'].values |
・eegのデータをplotして表示。
|
1 2 3 |
#欠損値を補完補正した脳波データを時系列グラフ表示 plt.subplots(figsize=(18,6)) plt.plot(eeg_ip) |
・eegのデータから、512件（４秒, eeg_512）と、8192件(64秒, eeg_8192)を取り出す。
|
1 2 |
eeg_512 = eeg_ip[1800:2312] eeg_8192 = eeg_ip[1800:9992] |
eeg_512をグラフ表示
|
1 2 3 4 5 |
#欠損値を補完補正した脳波データを時系列グラフ表示 plt.subplots(figsize=(18,6)) plt.ylim(-100, 100) plt.plot(eeg_512, linewidth =0.5) plt.savefig('spectral_analysis_2_fig_3.svg') |
・eeg_512の１秒分128件を表示
|
1 2 3 4 5 |
#欠損値を補完補正した脳波データを時系列グラフ表示 plt.subplots(figsize=(18,6)) plt.ylim(-100, 100) plt.plot(eeg_512[0:128], linewidth =0.5) plt.savefig('spectral_analysis_2_fig_4.svg') |
・eeg_8192のデータをグラフ表示
|
1 2 3 4 |
#欠損値を補完補正した脳波データを時系列グラフ表示 plt.subplots(figsize=(18,6)) plt.plot(eeg_8192, linewidth =0.5) plt.savefig('spectral_analysis_2_fig_5.svg') |
・FFT関数のパラメータをセット。ここではサンプリング周波数128とする。
|
1 2 3 4 5 6 |
#FFT関数のパラメータセット N=len(eeg_512) #総データポイント数 N2=len(eeg_8192) #総データポイント数 f_s = 128 #サンプリング周波数 T = 1.0/f_s #データ幅(秒) s_rate = 1/f_s #サンプリングレート |
・デシベル変換式を関数定義する
|
1 2 3 4 5 6 7 |
#dB変換 def db(x, dBref): y = 20 * np.log10(x / dBref) #変換式 return y #dB値を返す #FFT関数の中身の演算を確認 1.0/(2.0*T), N/2 |
Out: (64.0, 256.0)
・Scipyのfft()関数でeeg_512(4秒分）データをFFT演算する。
|
1 2 |
xfft = np.linspace(0, int(1.0/(2.0*T)), int(N/2)) yfft=fft(eeg_512) |
・単純なFFTでのパワースペクトルを求める。
|
1 2 3 4 5 6 |
#それぞれのスプリット脳波のパワースペクトル表示 plt.plot(xfft, 2.0/N*abs(yfft[0:int(N/2)]), linewidth =0.5, label="periodogram") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.grid() plt.yscale("log") plt.savefig('spectral_analysis_2_fig_6.svg') |
・デシベル変換する
|
1 2 |
dBref = 1 yfft_db = db(2.0/N*abs(yfft[0:int(N/2)]), dBref) |
・単純なFFT（ピリオドグラム）によるSPDをデジベル変換表示する。
|
1 2 3 4 5 6 |
#それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, yfft_db[0:int(N/2)], linewidth =0.5, label="Periodogram") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) #plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_7.svg') |
・ScipyのSignalモジュールを利用してHann窓を作成する。
|
1 2 3 4 |
han = sig.hann(N) #ハニング窓の作成 acf_han = 1/(sum(han)/N) #Amplitude Correction Factorを計算 print(acf_han) |
Out: 2.0039138943248536
・eeg_512にHann窓を適応して、FFT演算を行う。
|
1 2 |
#窓関数補正 fft_data_acf_han = acf_han*np.abs(fft(eeg_512*han)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 |
・Hann窓関数を適応したピリオドグラムによるPSDを表示する。
|
1 2 3 4 5 6 |
#それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_han[0:int(N/2)], linewidth =0.5, label="Hann") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_8.svg') |
・デシベル変換
|
1 2 |
dBref = 1 yfft_acf_db_han = db(fft_data_acf_han[0:int(N/2)], dBref) |
・Hann窓関数を適応したピリオドグラムによるデシベルPSDを表示する。
|
1 2 3 4 5 6 7 |
#それぞれのスプリット脳波のパワースペクトル表示 plt.plot(xfft, yfft_acf_db_han[0:int(N/2)], linewidth =0.5, label="Hann") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) #plt.yscale("log") plt.ylim(-50,25) plt.grid() plt.savefig('spectral_analysis_2_fig_9.svg') |
・単純なピリオドグラムとHann窓の適応とのPSDを比較する。
|
1 2 3 4 5 6 7 8 |
#それぞれのスプリット脳波のパワースペクトル表示 plt.plot(xfft, yfft_db, linewidth =0.5, label="periodogram", alpha=0.8) plt.plot(xfft, yfft_acf_db_han, linewidth =0.5, label="hann", alpha=0.8) plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.grid() #plt.yscale("log") plt.savefig('spectral_analysis_2_fig_10.svg') |
・Hann窓、Hamming窓、Blackman窓を作成する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
#複数n窓関数を用意 今回は、numpyで。 window_n = np.hanning(N) window_m = np.hamming(N) window_b = np.blackman(N) x = np.arange(0,N,1) #横軸 #窓関数の補正値 acf_bm=1/(sum(window_b)/N) F_abs_amp_bm = acf_bm*np.abs(fft(eeg_512*window_b)/(N/2)) #baclman窓：FFT後の数値に掛ければOK dBref = 1 yfft_acf_db_bm = db(F_abs_amp_bm[0:int(N/2)], dBref) |
・Blackman窓関数をアプライしたPSDをHann窓と比較するグラフ表示する。
|
1 2 3 4 5 6 7 8 9 10 |
#それぞれのスプリット脳波のパワースペクトル表示 plt.plot(xfft, yfft_acf_db_han, linewidth =0.5, label="hann", alpha=0.7) plt.plot(xfft, yfft_acf_db_bm, linewidth =0.5, label="blackman", alpha=0.7) plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.grid() plt.ylim(-50,25) #plt.yscale("log") plt.savefig('spectral_analysis_2_fig_11.svg') |
脳波複数解析セグメントのパワースペクトルの実践
・窓関数の選択ができるようにする。Hann窓を選択する。解析データウィンドウは512で、ウインドウ間は50%オーバーラップとする。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
dt = 1/128 #サンプリングレート split = 512/N2 #分析枠データの全体データ比率 overlap = 0.5 #オーバーラップ率 window = "hanning" #窓関数選択: hanning, hamming, blackman #窓の用意 if window == "hanning": window_select = "hanning" # ハニング窓 print("Hanningを選択") elif window == "hamming": window_select = "hamming" # ハミング窓 print("Hammingを選択") elif window == "blackman": window_select = "blackman" # ブラックマン窓 print("Blackmanを選択") else: window_select = "hanning" # ハニング窓 print("Hanningを選択") #fq_ave, F_abs_amp_ave, time, freq, power = FFT_main(time, eeg, dt, split, overlap, window_select) |
Hanningを選択
・脳波記録を解析セグメントごとにスプリットして解析できるようにする関数data_split()を定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
#データをFFT解析幅にスプリットする関数で、スプリットされたデータを配列として返す。 def data_split(_t, _x, _split, _overlap): split_data = [] one_frame_cnt = int(len(_t)*_split) #1フレームサンプル数 overlap_cnt = int(one_frame_cnt*_overlap) #オーバーラップサンプル数 print("フレームサンプル数: ", one_frame_cnt) print("オーバーラップサンプル数: ", overlap_cnt) start_S = 0 end_S = start_S + one_frame_cnt while True: t_cnt = _t[start_S:end_S] x_cnt = _x[start_S:end_S] split_data.append([t_cnt, x_cnt]) start_S = start_S + (one_frame_cnt - overlap_cnt) end_S = start_S + one_frame_cnt if end_S > len(_t): break return np.array(split_data) |
・FFT関数を定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
#FFT関数 def FFT(_data_input, _dt, _window_select): data_cnt_adj = len(_data_input[0]) #窓関数処理 if _window_select == "hanning": window_select = np.hanning(data_cnt_adj) #Hanning elif _window_select == "hamming": window_select = np.hamming(data_cnt_adj) #Hamming elif _window_select == "blackman": window_select = np.blackman(data_cnt_adj) #Blackman else: window_select = np.hanning(data_cnt_adj) #Hanning #窓関数後信号 x_windowed = _data_input[1]*window_select #FFT計算 F = np.fft.fft(x_windowed) F_abs = np.abs(F) F_abs_amp = F_abs / data_cnt_adj * 2 fq = np.linspace(0, 1.0/dt, data_cnt_adj) #窓関数補正 acf=1/(sum(window_select)/data_cnt_adj) F_abs_amp = acf*F_abs_amp #ナイキスト定数まで抽出 fq_out = fq[:int(data_cnt_adj/2)+1] F_abs_amp_out = F_abs_amp[:int(data_cnt_adj/2)+1] return [fq_out, F_abs_amp_out] |
・FFT結果をプロットする関数plot_TTFを定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
#FFTの結果をグラフ表示（平均値の表示に使用） def plot_FFT(_t, _x, _fq, _F_abs_amp): fig = plt.figure(figsize=(16, 6)) ax2 = fig.add_subplot(121) plt.title("time") plt.xlabel("time [s]") plt.ylabel("amplitude"+"[V]") plt.plot(_t, _x, linewidth =0.5) ax2 = fig.add_subplot(122) plt.title("freq") plt.xlabel('freqency(Hz)') plt.ylabel("amplitude"+"[V/rtHz]") #plt.xscale("log") plt.yscale("log") plt.plot(_fq, _F_abs_amp, linewidth =0.5) |
・スプリットされた脳波を表示させる関数plot_split_eeg()を定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
#スプリットされた脳波を描写 def plot_split_eeg(_split_data, _FFT_result_list, _frame_cnt, _FFT_result_cnt): #スプリットされた脳波データ fig3 = plt.figure(figsize=(24, 8)) plt.title("time") plt.xlabel("time [s]") plt.ylabel("amplitude"+"[V]") loc = 0 for _split_data_cnt, _FFT_result_cnt in zip(_split_data, _FFT_result_list): plt.xscale("linear") plt.yscale("linear") ax3 = fig3.add_subplot(_frame_cnt, 1, 1+loc) plt.plot(_split_data_cnt[0], _split_data_cnt[1]) loc +=1 plt.savefig('spectral_analysis_2_fig_12.svg') |
・スプリットされた脳波を全部まとめて表示させる関数plot_FFT_all()を定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
#個々のスプリット脳波データ def plot_FFT_all(_split_data, _FFT_result_list, _frame_cnt, _FFT_result_cnt, _overlap_cnt, _dBref): #個々のスプリットデータのスペクトラム fig4 = plt.figure(figsize=(12, 4)) ax4 = fig4.add_subplot(122) plt.title("freq") plt.xlabel('freqency(Hz)') plt.ylabel("amplitude"+"[V/rtHz]") #plt.xscale("log") plt.yscale("log") noc=1 for _split_data_cnt, _FFT_result_cnt in zip(_split_data, _FFT_result_list): plt.plot(_FFT_result_cnt[0], _FFT_result_cnt[1], label=noc, linewidth =0.5) plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=8) noc +=1 |
・データを三次元配置する関数set_3D_dataを定義する。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 |
#三次元データ配置 def set_3D_data(_split_data, _FFT_result_list, _frame_cnt, _FFT_result_cnt, _overlap_cnt): #三次元データ： T＝時間、Y＝周波数、z＝パワーを、リスト形式の配列に収める。 #時間軸配列へのデータ配置 T_pre= np.arange(0, 6, 6/_frame_cnt) T=[] for i in T_pre: T.append([i] *(1+_overlap_cnt)) print("時間配列2次元データ数：", len(T)) print("時間配列1次元データ数：", len(T[0])) #print(T) #周波数軸配列へのデータ配置 Y_pre= _FFT_result_cnt[0].tolist() Y=[] for i in T_pre: Y.append(Y_pre) #print(len(Y_pre)) print("周波数配列2次元データ数：", len(Y)) print("周波数配列1次元データ数：", len(Y[0])) #print(Y) #周波数パワー軸配列へのデータ配置 Z=[] for split_data_cnt, FFT_result_cnt in zip(_split_data, _FFT_result_list): Z_pre=10*np.log10(FFT_result_cnt[1]) Z.append(Z_pre.tolist()) print("周波数パワー配列2次元データ数：", len(Z)) print("周波数パワー配列2次元データ数：", len(Z[0])) #print(Z) #３次元データをnumpy配列へ変換 T=np.array(T) Y=np.array(Y) Z=np.array(Z) return T, Y, Z |
・4秒間512データの時間軸配列を設定する。
|
1 |
Time_adj= np.linspace(0, 4, N2) |
・eeg_8192（64秒分)のデータについて、SPD算定の処理を行う。ここでは31セグメントに分割される。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
#EEGデータのオーバーラップ分割。 split_data = data_split(Time_adj, eeg_8192, split, overlap) frame_cnt =len(split_data) print("総フレーム数: ", frame_cnt) #FFT演算を行う。 FFT_result_list = [] for split_data_cnt in split_data: FFT_result_cnt = FFT(split_data_cnt, dt, window) FFT_result_list.append(FFT_result_cnt) one_frame_cnt = int(len(time)*split) overlap_cnt = int(one_frame_cnt*overlap) #オーバーラップサンプル数 #print("FFT_result_list[0]: ", FFT_result_list[0]) #print("FFT_result_list[1]: ", FFT_result_list[1]) #パワー平均化 fq_ave = FFT_result_list[0][0] F_abs_amp_ave = np.zeros(len(fq_ave)) for i in range(len(FFT_result_list)): F_abs_amp_ave = F_abs_amp_ave + FFT_result_list[i][1] F_abs_amp_ave = F_abs_amp_ave/(i+1) |
フレームサンプル数: 512
オーバーラップサンプル数: 256
総フレーム数: 31
・スプリットされた脳波データを表示する。
|
1 2 |
#スプリット脳波の描画 plot_split_eeg(split_data, FFT_result_list, frame_cnt, FFT_result_cnt) |
・全体で平均化されたアンプリチュードとスペクトルを表示する。
|
1 2 3 |
#平均化されたスペクトラム plot_FFT(Time_adj, eeg_8192, fq_ave, F_abs_amp_ave) plt.savefig('spectral_analysis_2_fig_13.svg') |
・セグメントに分けて解析されたSPDを同じグラフ上にまとめて表示させる。
|
1 2 3 4 |
#個々のスプリットデータのスペクトラム dBref=1 plot_FFT_all(split_data, FFT_result_list, frame_cnt, FFT_result_cnt, overlap_cnt, dBref) plt.savefig('spectral_analysis_2_fig_14.svg') |
スペクトログラムの構築
|
1 2 |
#3Dデータのセット: 時間、周波数、周波数パワーの２次元配列：各次元のデータ数は一致する必要あり。 time2, freq, power = set_3D_data(split_data, FFT_result_list, frame_cnt, FFT_result_cnt, overlap_cnt) |
時間配列2次元データ数： 31
時間配列1次元データ数： 513
周波数配列2次元データ数： 31
周波数配列1次元データ数： 257
周波数パワー配列2次元データ数： 31
周波数パワー配列2次元データ数： 257
・単純なピリオドグラムからのSPDで描くスペクトログラム構築
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
nFFT = 512 # the length of the windowing segments f_s =128 # the sampling frequency noverlap=nFFT*0.5 # Pxx is the segments x freqs array of instantaneous power, freqs is # the frequency vector, bins are the centers of the time bins in which # the power is computed, and im is the matplotlib.image.AxesImage # instance window = np.blackman(nFFT) plt.figure(figsize=(4, 4)) Pxx, freqs, bins, im = plt.specgram(eeg_8192, NFFT=nFFT, Fs=f_s, scale='dB', window=window, noverlap=noverlap, cmap='jet') plt.colorbar() plt.savefig('spectral_analysis_2_fig_15.svg') |
・セグメントに分けて解析されたSPDを同じグラフ上にまとめて表示させる。
|
1 2 3 4 5 |
import matplotlib.pyplot as plt #plt.plot(freqs, 10*np.log10(Pxx)) Bref=1 plt.plot(freqs, db(Pxx, Bref), linewidth =0.5) plt.savefig('spectral_analysis_2_fig_16.svg') |
・より高度なスペクトログラム構築のため、libtfrライブラリ、librosaライブラリ等を導入する
|
1 2 3 4 |
import libtfr from librosa import display from scipy.signal import chirp import librosa |
・libtfrを用いたスペクトログラム算定のためのパラメータ定義を行う。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
nfft = 512 f_s = 128 shift = nfft/2 K = 6 tm = 6.0 flock = 1 tlock = 5 #nfft = 1024 #f_s = 128 #shift = nfft/16 #K = 6 #tm = 6.0 #flock = 1 #tlock = 5 # load signal of dimension (npoints,) S = libtfr.tfr_spec(eeg_8192, nfft, shift, f_s, K, tm, flock, tlock) print(np.shape(S), np.max(S), np.min(S)) S2 = librosa.amplitude_to_db(S, ref=np.max, top_db =100) print(np.shape(S2), np.max(S2), np.min(S2)) #S3 = librosa.power_to_db(S, ref=np.max, top_db =100) #print(np.shape(S3), np.max(S3), np.min(S3)) |
(257, 32) 231388.4331020175 0.0
(257, 32) 0.0 -100.0
(257, 32) 0.0 -100.0
・libtfrのtfr_spec()関数を用いた周波数時間再割当て法によるSPDによるスペクトログラム構築
|
1 2 3 4 |
fig, ax = plt.subplots(figsize=(4,4)) display.specshow(S2, x_axis='time', y_axis='hz', cmap='jet', sr=f_s, hop_length=nfft) plt.colorbar() plt.savefig('spectral_analysis_2_fig_17.svg') |
|
1 2 3 4 5 |
list=[] for i in range(0, 256): list.append(np.mean(S2[i])) Fq= np.linspace(0, 64, 256) |
・はじめのデータセグメントのPSDを表示させる。
|
1 2 |
plt.plot(Fq, list, linewidth =0.8) plt.savefig('spectral_analysis_2_fig_18.svg') |
・libtfrのマルチテーパ法mfft_dpss()を用いて、スペクトログラムを作成する。
|
1 2 3 4 5 6 7 8 9 |
# generate a transform object with size 512 samples and 5 tapers for short-time analysis D = libtfr.mfft_dpss(512, 12, 5) # complex STFT with frame shift of 10 samples Z = D.mtstft(eeg_8192, 64) # spectrogram with frame shift of 10 samples P = D.mtspec(eeg_8192, 64) print(np.shape(Z), np.max(Z), np.min(Z)) print(np.shape(P), np.max(P), np.min(P)) |
(257, 121, 5) (290.40539223527196+0j) (-284.31933845774296+21.55352222591088j)
(257, 121) 50412.716942419465 0.002685783319892681
|
1 2 3 4 5 6 |
P = librosa.amplitude_to_db(P, ref=np.max, top_db =100) fig, ax = plt.subplots(figsize=(4,4)) display.specshow(P, x_axis='time', y_axis='hz', cmap='jet', sr=f_s, hop_length = nfft) plt.colorbar() plt.savefig('spectral_analysis_2_fig_21.svg') |
・マルチテーパ法による第一セグメントのSPDを表示させる。
|
1 2 3 4 5 6 |
list3=[] for i in range(0, 256): list3.append(np.mean(P[i])) plt.plot(Fq, list3, linewidth =0.8) plt.savefig('spectral_analysis_2_fig_22.svg') |
・第一セグメントのPSDについて、libtfrの２つの方法（時間再割当て法、マルチテーパ法）の２つを比較する。
|
1 2 3 4 5 |
fig = plt.figure() plt.plot(Fq, list, linewidth =0.8) #plt.plot(Fq, list2) plt.plot(Fq, list3, linewidth =0.8) plt.savefig('spectral_analysis_2_fig_23.svg') |
マルチテーパ法の分解
・マルチテーパ法ライブラリmtspecを導入
|
1 2 3 |
from mtspec import dpss tapers, lamb, theta = dpss(N, 2.5, 5) |
・Slepianシークエンス5つをを描画。
|
1 2 3 |
for i in range(5): plt.plot(tapers[:, i]) plt.savefig('spectral_analysis_2_fig_24.svg') |
・第1Slepianテーパ関数処理
|
1 2 |
plt.plot(tapers[:, 0]) plt.savefig('spectral_analysis_2_fig_25.svg') |
|
1 2 |
plt.plot(tapers[:, 0]*eeg_512) plt.savefig('spectral_analysis_2_fig_26.svg') |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
acf_0 = 1/(sum(tapers[:, 0])/N) #Amplitude Correction Factorを計算 print(acf_0) fft_data_acf_0 = np.abs(fft(tapers[:, 0]*eeg_512)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 #dBref = 1 #yfft_acf_db_0 = db(fft_data_acf_0[0:int(N/2)], dBref) #それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_0[0:int(N/2)], linewidth =0.5, label="taper #1") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_27.svg') |
28.843575280051258
・第2Slepianテーパ関数処理
|
1 2 |
plt.plot(tapers[:, 1]) plt.savefig('spectral_analysis_2_fig_28.svg') |
|
1 2 |
plt.plot(tapers[:, 1]*eeg_512) plt.savefig('spectral_analysis_2_fig_29.svg') |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
acf_1 = 1/(sum(tapers[:, 1])/N) #Amplitude Correction Factorを計算 print(acf_1) fft_data_acf_1 = np.abs(fft(tapers[:, 1]*eeg_512)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 #dBref = 1 #yfft_acf_db_1 = db(fft_data_acf_1[0:int(N/2)], dBref) #それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_1[0:int(N/2)], linewidth =0.5, label="taper #2") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_30.svg') |
・第3Slepianテーパ関数
|
1 2 |
plt.plot(tapers[:, 2]) plt.savefig('spectral_analysis_2_fig_31.svg') |
|
1 2 |
plt.plot(tapers[:, 2]*eeg_512) plt.savefig('spectral_analysis_2_fig_32.svg') |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
acf_2 = 1/(sum(tapers[:, 2])/N) #Amplitude Correction Factorを計算 print(acf_2) fft_data_acf_2 = np.abs(fft(tapers[:, 2]*eeg_512)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 #dBref = 1 #yfft_acf_db_2 = db(fft_data_acf_2[0:int(N/2)], dBref) #それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_2[0:int(N/2)], linewidth =0.5, label="taper #3") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_33.svg') |
・第4Slepianテーパ関数処理
|
1 2 |
plt.plot(tapers[:, 3]) plt.savefig('spectral_analysis_2_fig_34.svg') |
|
1 2 |
plt.plot(tapers[:, 3]*eeg_512) plt.savefig('spectral_analysis_2_fig_35.svg') |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
acf_3 = 1/(sum(tapers[:, 3])/N) #Amplitude Correction Factorを計算 print(acf_3) fft_data_acf_3 = np.abs(fft(tapers[:, 3]*eeg_512)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 #dBref = 1 #yfft_acf_db_3 = db(fft_data_acf_3[0:int(N/2)], dBref) #それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_3[0:int(N/2)], linewidth =0.5, label="taper #4") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_36.svg') |
・第5Slepianテーパ関数処理
|
1 2 |
plt.plot(tapers[:, 4]) plt.savefig('spectral_analysis_2_fig_37.svg') |
|
1 2 |
plt.plot(tapers[:, 4]*eeg_512) plt.savefig('spectral_analysis_2_fig_38.svg') |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
acf_4 = 1/(sum(tapers[:, 4])/N) #Amplitude Correction Factorを計算 print(acf_4) fft_data_acf_4 = np.abs(fft(tapers[:, 4]*eeg_512)/(N/2)) #dataは時間波形, FFT,正規化,振幅補正を実施 dBref = 1 yfft_acf_db_4 = db(fft_data_acf_4[0:int(N/2)], dBref) #それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, fft_data_acf_4[0:int(N/2)], linewidth =0.5, label="taper #5") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) plt.yscale("log") plt.grid() plt.savefig('spectral_analysis_2_fig_39.svg') |
・５つのテーパ処理からFFTされた結果を平均化して、PSD表示
|
1 2 3 4 |
yfft_mt=(fft_data_acf_0+fft_data_acf_1+fft_data_acf_2+fft_data_acf_3+fft_data_acf_4) dBref = 1/5 yfft_acf_db_mt = db(yfft_mt, dBref) |
|
1 2 3 4 5 6 7 |
#それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, yfft_acf_db_mt[0:int(N/2)], linewidth =0.5, label="multitaper") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) #plt.yscale("log") plt.grid() plt.ylim(-50,25) plt.savefig('spectral_analysis_2_fig_40.svg') |
・マルチテーパ法によるPSDと、ピリオドグラムによるPSDを比較表示させる。
|
1 2 3 4 5 6 7 8 9 |
#それぞれのスプリット脳波のパワースペクトル表示（ここでは、影響ない） plt.plot(xfft, yfft_db, linewidth =0.5, label="periodogram", alpha=0.8) plt.plot(xfft, yfft_acf_db_mt[0:int(N/2)], linewidth =0.5, label="multitaper") plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=12) #plt.yscale("log") plt.grid() plt.ylim(-50,25) plt.savefig('spectral_analysis_2_fig_41.svg') |
