Python (Scipy) による神経活動データ解析入門・その２ | プロットあれこれ

前回定義した配列 ecog_mV を、いろいろプロットを変えて詳しく調べていく。

プロット範囲の指定

前回の例では2,000秒間ほどのデータがすべてプロットされていたが、より細かい時間スケールの変動を見るため、プロット範囲を短く設定する。

IN [*]: xlim( 0, 1 )

これでデータの最初の1秒間がプロットされる。同様に、

IN [*]: xlim( 1000, 1001 )

とすると、1,000秒から1,001秒までの1秒間がプロットされる。
(y軸の範囲の指定は ylim( y_min, y_max ) で、x軸y軸の同時指定は axis( [x_min, x_max, y_min, y_max] ) で行う。)

0−1秒の時間帯は麻酔剤注入前の覚醒状態(awake period)、1,000−1,001秒は麻酔が十分効いている状態(sleep period)と考えられるが、これらの間の ECoG 信号を比較してみると、前回確認した、変動の振幅の大きさだけでなく、変動の時間スケールも全く異なっていることが分かる。すなわち、0−1秒では1秒間に20回ほどの規則的な振動が見られるが、1,000−1,001秒ではそれが見られず、大振幅のゆっくりとした振動に変化している。

対数プロット

信号の変動の時間スケールを定量的に知りたい場合、特に、特定の周波数の振動の有無を知りたい場合、まず調べるのは信号のパワースペクトル密度である。SciPy(+Matplotlib) では、コマンド一つでスペクトル密度の計算とプロットがいっぺんにできてしまう。
今回は覚醒状態・麻酔状態で別々にスペクトル密度を見たいので、元のデータのうち最初の500秒間を覚醒状態のサンプル、1,000秒から1,500秒までを麻酔状態のサンプルとして別の配列にコピーし、それぞれについてスペクトル密度を計算する。

IN [*]: ecog_awake = ecog_mV[0:500*1000]
IN [*]: ecog_sleep = ecog_mV[1000*1000:1500*1000]
IN [*]: psd( ecog_awake, Fs=1000, label='awake' )
IN [*]: psd( ecog_sleep, Fs=1000, label='sleep' )

psd に与えているデータ以外の引数は、Fs がサンプリング周波数、label は後で凡例を作るときに使うラベル。
この方法では、縦軸の単位は自動的に dB/Hz でプロットされる。デシベル(dB)は信号強度を対数変換したものなので、つまりこれは実質的に片対数プロットになっている。横軸も対数スケールにするには、

IN [*]: xscale( 'log' )

とする。その結果が以下のグラフ。（ただし上記の設定以外に、grid()でグリッドを、legend()で凡例を表示させている。）

覚醒状態のグラフ(awake)は、右下がりの直線上に 15-20 Hz の周波数帯でピークが乗った形になっている。両対数グラフで直線になるということは、元の信号がいわゆる1/fゆらぎを持っていることを示している。15-20 Hz のピークは、最初に示したグラフで見えた細かい振動に対応しており、ヒトEEGで閉眼安静時に顕著に見られるアルファ波に相当するものと思われる。麻酔状態(sleep)ではこのピークは消失し、代わりに 40-80 Hzあたりに肩のような盛り上がりが見える。