Python (Scipy) による神経活動データ解析入門・その１ | .mat ファイルの読み込みからデータのプロットまで
前提として、Python, iPython, NumPy, SciPy, Matplotlib はインストール済み、Neurotycho のデータはダウンロード済みであるとする。
SciPyの起動
まずはシェルプロンプトから
% ipython -pylab
で iPython(Python の対話的シェル)を起動。"-pylab" を付けることで、MATLAB 風のデータプロットのコマンドをサポートするパッケージ(Matplotlib) を組み込んだモードになる。
立ち上がった iPython のプロンプトから
In [*]: from numpy import * In [*]: from scipy import *
で NumPy, SciPy のモジュールをすべてインポート。これでひと通りの機能が使えるようになる。
データの読み込み
MATLAB形式のデータ(.mat)を読み込むには scipy.io.loadmat() を使う。まずこれをインポートする。
In [*]: from scipy.io import loadmat
Neurotycho > download > sleep task では ECoG データを４つのファイルに分割して公開しているが、まず最初のデータファイル(ECoG_100604-1_32.mat) を読み込んでみる。
In [*]: dataset = loadmat( 'ECoG_1006040-1_32.mat' )
とし、データファイルの内容を辞書オブジェクト(＝連想配列) dataset にコピーする。
データの確認
まずおもむろに
In [*]: dataset
とすると、辞書オブジェクト dataset の概要が表示され、このデータが以下の4つの要素を持つことが分かる:
- ラベル 'X': int16型の配列オブジェクト
- ラベル '__globals__': 空のリスト
- ラベル '__header__': ファイル情報を記述した文字列
- ラベル '__version__': バージョン数を記述した文字列
これらのうち、'X' とラベルされた最初の要素がECoGのデータ本体であると思われる。この要素のみを取り出してみる。
In [*]: ecog_data = dataset['X']
以下、この配列 dataset の情報を得るためのメソッドをいくつか紹介していく。まず、配列の次元を知るには dim メソッドを使う。
In [*]: ecog_data.dim Out[*]: 2
となり、2次元配列であることが分かる。各次元の大きさを知るには shape メソッドを使う。
In [*]: ecog_data.shape Out[*]: (32, 2062570)
となり、データファイルの名前や、データに付属している ReadMe.txt の内容から推測すると、第1次元が ECoG のチャンネル、第2次元がデータ点に対応しているらしいと分かる。すなわち、ecog_data[n] を参照すればn番目のチャンネルのデータが取り出せることになる(ただし 0 <= n <= 31 に注意)。
データのプロット
データを入手してまず最初にすべきことは、それを目で見て確認すること。そのためにデータをプロットしてグラフにする。最も単純な方法は
In [*]: plot( ecog_data[0] )
とすれば、配列 ecog_data[0] の要素をそのままプロットできる。だがこれだと、横軸が配列の要素番号となってしまい、時間との対応がつかない。また縦軸の単位もわからない。
Neurotycho > download の記述に、ECoG データの単位はマイクロボルト(uV)、サンプリング周波数は 1kHz だとあるので、まず、縦軸をミリボルト単位でプロットするため、データの単位を変更して新しい配列 ecog_mV に移す。
In [*]: ecog_mV = ecog_data[0] / 1000.0
横軸に関しては、秒単位でプロットするため、以下のようにして横軸のグリッドを表す配列 n_sec を生成する。
In [*]: n_sec = arange( 0, ecog_data[0].shape[0] ) In [*]: n_sec = n_sec / 1000.0
やっていることは、まず 0 から ecog_data[0] の要素数まで、1づつ増える整数の配列を生成し、それをサンプリング周波数で割っている。この配列 n_sec と先程の単位変更したデータ ecog_mV を使うと、
In [*]: plot( n_sec, ecog_mV )
として、横軸が時間(sec)、縦軸が ECoG 信号(mV)のグラフが描ける。軸にラベルを付けるには
In [*]: xlabel( 'Time (sec)' ) In [*]: ylabel( 'ECoG signal (mV)' )
でそれぞれ横軸と縦軸にラベルが付けられる。結果が以下のグラフ。
ReadMe.txt に、測定開始後631秒の時点で麻酔剤を注入したとあるが、その100秒後あたりから ECoG 信号の振幅が増大していることが分かる。また、それ以前・以後の時間帯では信号の強度は概ね定常であることも分かる。
