MNE-pythonのversionチェック¶
2017年4月の時点でMNE-pythonはpython 2.7で動きます。python 3でも動きますが、三次元画像処理はできません。
In [1]:
import sys
sys.version_info
Out[1]:
MNEの準備とraw-fifファイルと読み書き¶
In [2]:
import mne # mneを組み込む
from mne.datasets import sample
mne.set_log_level('INFO');
data_path=sample.data_path()
raw_fname=data_path+'\\MEG\\sample\\sample_audvis_filt-0-40_raw.fif'
print(raw_fname)
In [3]:
raw=mne.io.read_raw_fif(raw_fname) # raw_fifファイルの読み込み
print(raw)
print(raw.info)
In [4]:
print(raw.ch_names)
In [5]:
start,stop=raw.time_as_index([100,115]) # 100～115秒で切り出し
data,times=raw[:,start:stop]
print(data.shape) # MATLABだとsize(data)
print(times.shape) # MATLABだとsize(times)
data,times=raw[2:20:3,start:stop] # MATLABだと2:3:20
raw.plot();# 初期設定はduration=10秒 starrt=0.0秒 N_channels=20 chずつ・・・などです
jupyterでは絵になってしまうのでできませんが、カーソルバーで表示されるチャンネルや時間幅を変更できます。ProjボタンでSSPのon/offが選択できます。
In [6]:
picks=mne.pick_types(raw.info,meg=True,eeg=False,stim=True,exclude='bads')
raw.save('sample_audvis_meg_raw.fif',tmin=0,tmax=150,picks=picks,overwrite=True) # 0～150秒のデータを保存
raw-fifから加算¶
In [7]:
events=mne.find_events(raw,stim_channel='STI 014'); # トリガー信号のチャンネルを指定します。
print(events[:5])
In [8]:
mne.set_config('MNE_STIM_CHANNEL','STI101',set_env=True) # トリガーチャンネルの初期設定
In [9]:
event_id=dict(aud_l=1,aud_r=2) # トリガー信号 聴覚左が1、聴覚右が2
tmin,tmax=-0.2,0.5 # -0.2～0.5秒
In [10]:
raw.info['bads']+=['MEG 2443','EEG 053'];
In [11]:
picks=mne.pick_types(raw.info,meg=True,eeg=True,eog=True,stim=False,exclude='bads');
In [12]:
baseline=(None,0) # 基線の補正
In [13]:
reject=dict(grad=4000e-13,mag=4e-12,eog=150e-6) # 高振幅は除外
In [14]:
epochs=mne.Epochs(raw,events,event_id,tmin,tmax,proj=True,picks=picks,baseline=baseline,preload=False,reject=reject
)
print(epochs);
In [15]:
epochs_data=epochs['aud_l'].get_data()
print(epochs_data.shape)
In [16]:
from scipy import io
print(type(epochs_data))
print(epochs_data.shape)
io.savemat('epochs_data.mat',dict(epochs_data=epochs_data,oned_as='row')) # MATLABのデータとして保存
MATLABで開と55×365×106の三次元配列epochs_rowという文字列のdataとoned_asが読み込まれます。
In [17]:
epochs.save('sample-epo.fif') # FIFファイル形式で保存できます。
sample-epo.fifはデータ形式が三次元配列のFIFFファイルです。 Neuromagのアプリケーションでは読めません。 Source Modellingで開くとIncorrect number of channels in data matrix!と表示され読めません。 Graphで開くとUNKNOWN PACKING FORMAT 0 in file /home/neurosurgery/Deksotp/sample-epo.fif Reported by Graph routine df_open_continuousと表示されます。
In [18]:
del epochs; # 変数を消去します。
save_epochs=mne.read_epochs('sample-epo.fif') # 読み込みます。
In [19]:
evoked=save_epochs['aud_l'].average() # 加算平均します。
print(evoked)
evoked.plot();
試行毎の最大値を求めます。
In [20]:
print(len(save_epochs['aud_l'])) # aud_lの回数
max_in_each_epoch=[e.max() for e in save_epochs['aud_l']] # eはsave_epochs['aud_l][0～55].get_data()
print(max_in_each_epoch[:4]) # 最初の4つだけ表示
x=save_epochs['aud_l'][0].get_data();
print(x.max())
脳波のスケールが大きいんで、脳波ばかりが選ばれてます。
加算波形の読み書き・波形の差分¶
In [21]:
evoked_fname=data_path+'\\MEG\\sample/sample_audvis-ave.fif'
evoked1=mne.read_evokeds(evoked_fname,condition='Left Auditory',baseline=(None,0),proj=True) # projはSSP
In [22]:
evoked2=mne.read_evokeds(evoked_fname,condition='Right Auditory',baseline=(None,0),proj=True)
In [23]:
print([evoked1.nave,evoked2.nave]) # 左聴覚刺激と右聴覚刺激の回数
contrast=mne.combine_evoked([evoked1,evoked2],weights=[0.5,-0.5])
contrast=mne.combine_evoked([evoked1,-evoked2],weights='equal') # この書き方も可能です
print(contrast)
In [24]:
average=mne.combine_evoked([evoked1,evoked2],weights='nave')
print(contrast)
In [25]:
epochs_eq=save_epochs.copy(); # 参照渡しなのでコピーにします。
epochs_eq=epochs_eq.equalize_event_counts(['aud_l','aud_r'])#左右の刺激回数を同じにする
print(epochs_eq) # ２つの返り値あり
print(epochs_eq[0]) # １つめの帰り値
print(epochs_eq[1]) # ２つめの返り値 余った試行
evoked1,evoked2=epochs_eq[0]['aud_l'].average(),epochs_eq[0]['aud_r'].average()
print(evoked1)
print(evoked2)
contrast=mne.combine_evoked([evoked1,evoked2],weights=[0.5,-0.5])
print(contrast)
In [26]:
epochs_eq=save_epochs.copy(); # 参照渡しなのでコピーにします。
epochs_eq=epochs_eq.equalize_event_counts(['aud_l','aud_r'])
時間周波数解析 induced power and inter trial coherence¶
In [27]:
import numpy as np
n_cycles=2 # Morlet waveletの周波数
freqs=np.arange(7,30,3) # MATLABだと7:3:30
from mne.time_frequency import tfr_morlet
power,itc=tfr_morlet(save_epochs,freqs=freqs,n_cycles=n_cycles,return_itc=True,decim=3,n_jobs=2) # n_jobsはCPUのスレッド数
print(power)
print(itc) # inter trial coherence 試行間でのcoherence
power.plot([power.ch_names.index('MEG 1332')]);
itc.plot([itc.ch_names.index('MEG 1332')]);
逆問題 MNEとdSPM¶
In [28]:
from mne.minimum_norm import apply_inverse,read_inverse_operator
fname_inv=data_path+'\\MEG\\sample\\sample_audvis-meg-oct-6-meg-inv.fif'
inverse_operator=read_inverse_operator(fname_inv)
In [29]:
snr=3.0
lambda2=1.0/snr**2
method="dSPM"
stc=apply_inverse(evoked,inverse_operator,lambda2,method) # まずは加算波形の電流源推定
In [30]:
print(stc) # Source Time Course 電流源の経時変化
print(stc.subject)
print(type(stc.data),stc.data.shape,stc.shape); # 頂点×時間
print(type(stc.times),stc.times.shape); # 時間
print(type(stc.vertices),len(stc.vertices));# １つの頂点に自由度2の電流双極子
stc.save('mne_dSPM_inverse')
In [31]:
fname_label=data_path+'\\MEG\\sample\\labels\\Aud-lh.label'
label=mne.read_label(fname_label)
print(label)
print(label.pos.shape) # 座標
print(label.values.shape) # 頂点の値
自発波形データからメッシュの一部だけの電流源推定を行います。
In [32]:
from mne.minimum_norm import apply_inverse_raw
start,stop=raw.time_as_index([0,15]) # 0～15秒の時間幅を取得
stc=apply_inverse_raw(raw,inverse_operator,lambda2,method,label,start,stop)
In [33]:
print(stc)
stc.save('mne_dSPM_raw_inverse_Aud')
In [34]:
import matplotlib.pyplot as plt
%matplotlib inline
fig=plt.figure();
plt.plot(stc.times,stc.data.T,'b');
plt.xlim([stc.times[0],stc.times[-1]]);
