全身麻酔中の脳波スペクトル解析１：pythonによる基礎編

Published on by anesth

spectral_analysis_1
データファイル：eeg_bis_10_20_2.tsv

JupyterNotebook：spectral_analysis_1.ipynb

窓関数とマルチテーパ法
パワースペクトル（またはパワースペクトル密度の場合はPSD）として知られる、周波数関数としての信号内のパワー分布は、離散フーリエ変換（DFT）を使用して推定する。高速フーリエ変換アルゴリズム（FFT）を使用したパワースペクトルの単純な推定は、ピリオドグラムと呼ばれる。
ナイーブなピリドグラムによる推定では、解析に用いる時間制限信号内のすべてのサンプルがそのまま取得され（暗黙的に1倍）、この時間制限信号外のすべてのサンプルは無視される（暗黙的に0が乗算）。これは、信号が矩形窓（ボックスカー）ウィンドウでサンプルごとに乗算された場合に発生することと同様である。時間領域で信号をボックスカーウィンドウで乗算することは、周波数ドメインでボックスカーウィンドウのスペクトルで信号を解析することによる「スペクトル漏れ」を分析しながら、窓関数やマルチテーパ法の周波数特性について確認する。
・まずは、解析に必要となる以下のモジュール/関数をインポートする。
stats modelsは統計解析ライブラリ
nitimeは脳科学イメージングライブラリ
scipyは科学技術計算ライブラリ
matplotlibはグラフィックライブラリ
numpyは数学演算ライブラリー

出力図をSVGフォーマットでファイル保存するために、以下のコードを実行する。

・時系列データ数128のうち、zero_fracで指定した比率で両端が0、それ以外は１となる矩形窓（ボックスカー窓）関数boxcar_windowを作成する。
周波数解析における窓関数の役割を理解するために、まず窓を明示的に用いない（＝矩形窓）が
矩形窓（＝窓なし）関数を作成し、矩形波を描く


・Scipyライブラリのfftpackのfft()関数と、ゼロ周波数を中心にシフトさせて相対的な周波数に対するパワースペクトルを求めるfftshift()関数を使って矩形窓（ボックスカー窓）自体のPSDを算定する。


・Scipyライブラリのsignalモジュールを利用して、Hann， Hamming， Blackman窓を取得する。Hann, Blackmanはゼロで始まり、ゼロで終わる。Hammingはゼロに近づくがゼロではない特徴を持つ。


・窓関数自体の周波数解析を行いPSDを求めて、矩形窓（=窓なし ボックスカー）と比較する。中心の周波数の両翼の周波数が抑制されることが、広いダイナミックレンジの構築に役立つことが理解できる。

・中心周波数部分を拡大してみることで、周波数解像度に影響を与える中心周波数周辺のスペクトルの幅（漏れ）を理解する。Blackman窓では、就寝周波数の幅はHann窓やHamming窓と比べて広く、その両翼の周波数は低い。つまり周波数の解像度を犠牲にして、ダイナミックレンジを広げていることが解る。

マルチテーパ法による周波数解析

・マルチテーパ法による周波数解析について理解を進める。マルチテーパ法ライブラリmtspecのdpss()を導入する。 128データ数で、5つのSlepianシークエンスを求める。


・上記Slepian関数を矩形窓（ボックスカー）と同じスケールでグラフ化してみる。


・各Slepian関数の周波数解析を行い、矩形波のそれとPSDを比較する。


・上記の中心部分を拡大すると、各Slepianの特性とマルチテーパという複数のテーパ関数の組み合わせにより、周辺エッジ部のスペクトルも補完されることが理解できる。

自己回帰シークエンス由来の波形解析によるパワースペクトル解析
・自己回帰シークエンスにより真のパワースペクトル密度と、様々な解析方法による影響とを比較することで、利点、欠点を知ることができる。 ・dB変換する関数dB()を自作導入する。（10ではなくて20を乗じているのは、FFTで無視される波形のマイナス部分のパワーを算定するために２倍している）

・自己回帰移動平均モデル y_t= 1+ -2.7607 * y_1 + 3.8106 * y_2 + -2.6535* y_3 + 0.9238** y_4　を導入する。

・自己回帰移動平均モデルで生成された人工の波形を描画してみる。


・この人工波形の真のスペクトルをnitimeライブラリのfreq_response()関数で求める。結果は、複素数フーリエ変換式に基づいて、複素数で返される。そのためパワーを求めるために絶対値変換を行い、実数としてのパワーを求める。

・自己回帰シークエンスの波形自体をnitimeライブラリのナイーブピリオドグラム関数periodogram()でピリオドグラムPSDを求めて、真のスペクトルと比較する。40Hz以上の領域でスペクトルのズレが確認できた。いわゆる矩形窓のエッジの不自然な処理からこの漏れが発生していると考えられる。


・上記方法に、Hann窓関数を適応させてみる。前述の40Hz以上のズレが補正されてなくなることが解る。

・Nitimeのmulti_taper_psd（）関数を利用して、マルチテーパ方によるPSDを求める。得られたPSDは、真のPSDに近づくことが理解できる。
In [35]:


・χ二乗モデルによる95％信頼区間2 * Kmaxの自由度を計算します。


・反復法（[Thomson2007]）を使用することで、特定の信号の実際のスペクトル濃度に応じて、さまざまなテーパーの重みを適応的に設定できる。


上で述べたように、スペクトル自体の推定に加えて、スペクトルの信頼区間の推定は、ジャックナイフ法を使用して生成できます。 関数内でオプション指定 jackknife=Trueとする。
ジャックナイフ法（リーブワンアウト測定のことで、分布内のすべてのサンプル（すべてのテーパースペクトル）から平均化されたパラメーター推定値に対して、 1つの測定値（1つのテーパースペクトル）を除外することから得られる：
simple sample estimate θ? =1/n?iYi （分布内のすべてのサンプル（すべてのテーパースペクトル）から平均化されたパラメーター推定値）
リーブワンアウト測定 θ?−i=1/(n−1)?k≠i Yk （これは推定値のグループを定義し、各サンプルは1つの測定値（1つのテーパースペクトル）を除外） pseudovalues θ?_i=nθ? −(n−1)θ?−i
ジャックナイフ推定量は次のように計算。 θ? =1/n?iθ?_i = nθ? −(n−1)/n ?iθ?−i
この推定量に関して、真のパラメータθについてほぼスチューデントのt分布として分布し、標準誤差は次のように定義される。 s^2=(n−1)/n?i(θ?_i−θ?)^2
そして、使用されるテーパーの数に依存する自由度：Kmax-1

・さらに、「適応」フラグがTrueに設定されている場合、スペクトルのバイアスを修正するために反復適応法を使用する。そして、重みの適応推定とジャックナイフ現象を組み合わせる。


サイン波、コサイン波を用いた周波数解析の特徴理解
・5Hzの周期を持つ256データポイントのサイン波、コサイン波を作成する。

・サイン波を描画：ゼロで始まりゼロで終わる。

・コサイン波を描画：１で始まり、１で終わる。

・上記サイン・コサイン波の周波数解析でSPDを求める。両者は重なってほとんど差が無いように一見みえる。


・中心周波数5Hz領域を拡大してみる。

・デシベル変換して、表示させてみると、両者の中心周波数外の差が理解できる。


・次に前述の以下の解析窓を導入してみる。まずは窓関数をグラフ表示させる。

・窓関数とマルチテーパ法をサイン波形にかけ合わせて、FFTを行いPSDを求める。窓関数により、中心周波数（５Hz)以外の領域の周波数パワーが抑制されることが理解できる。

・窓関数とマルチテーパ法をコサイン波形にかけ合わせて、FFTを行いPSDを求める。窓関数により、中心周波数（５Hz)以外の領域の周波数パワーが抑制されることが理解できる。また上記サイン波とは異なることも解る。

・サイン波とコサイン波のPSDの違いを把握するため、単純なピリオドグラムによるPSDだけを比較してみる。

・サイン波PSDの５Hzの中心周波数部分を拡大表示してみる。


・コサイン波PSDの５Hzの中心周波数部分を拡大表示してみる。


・サイン波PSDの５Hzの中心周波数部分を拡大表示をデジベル変換で表示してみる。マルチテーパ法が中心周波数の解像度を犠牲にしていることが理解できる。

・コサイン波PSDの５Hzの中心周波数部分を拡大表示をデジベル変換で表示してみる。マルチテーパ法が中心周波数の解像度を犠牲にしていることが理解できる。

ノイズの載ったサイン波への周波数解析の影響
・サイン波に対して、ノイズデータを載せて、それに対して周波数解析を行うことでノイズへの影響を調べる。ノイズはnumpyのrandom.normal正規分布（平均値ゼロ、標準偏差0.1）で生成して、サイン波に付加。

・ノイズの載ったサイン波を生成した。


・サイン波＋ノイズを単純に周波数解析して、PSDを求める。


・サイン波＋ノイズを単純に周波数解析してPSDを求めたものの5Hz周辺の拡大図


・上記PSDをデジベル表示する。


・次にサイン波＋ノイズに窓関数を当てはめて、PSDを求める。

・サイン波＋ノイズに対して、矩形窓、窓関数を当てはめた周波数解析に基づくPSDを比較する

・サイン波＋ノイズに対して、矩形窓、窓関数、さらにマルチテーパ法を当てはめた周波数解析に基づくPSDを比較する。マルチテーパ法では、他の窓関数よりも広いダイナミックレンジを確保していることが理解できる。

・サイン波＋ノイズに対して、矩形窓、窓関数、さらにマルチテーパ法を当てはめた周波数解析に基づくPSDを比較する。５Hz周辺の拡大表示。


・サイン波＋ノイズに対して、矩形窓、窓関数、さらにマルチテーパ法を当てはめた周波数解析に基づくPSDを比較する。５Hz周辺の拡大デシベル表示。