内容
MNE-Pythonのtutorialのこの部分を試してみました。TutorialのPreprocessingの復習のような内容です。
エポックの切り出しと加算平均
データの読み込みをします。
import
os.path as op
import
numpy as np
import
mne
data_path=mne.datasets.sample.data_path()
fname=op.join(data_path,'MEG','sample','sample_audvis_raw.fif')
raw=mne.io.read_raw_fif(fname,add_eeg_ref=False)
raw.set_eeg_reference()
# set EEG average reference
エポックを経時的に表示します。
order=np.arange(raw.info['nchan'])
order[9]=312
#
order[312]=9
# トリガーチャンネルを10番目のチャンネルとして表示
raw.plot(n_channels=10,order=order)
エラーが出ますが無視します。
C:\Users\akira\Anaconda3\lib\site-packages\mne\viz\raw.py:389:
FutureWarning: elementwise comparison failed; returning scalar
instead, but in the future will perform elementwise comparison
if order in
['selection', 'position']:
STI
014がトリガーチャンネルです。
エポックの切り出しをします。
events=mne.find_events(raw)
print(events)
[[ 27977 0 2]
[ 28345 0 3]
[ 28771 0 1]
...,
[154080 0 1]
[154486 0 4]
[168672 0 32]]
トリガー別にエポックを経時的に色別表示させます。
event_id={'Auditory/Left':1,'Auditory/Right':2,'Visual/Left':3,'Visual/Right':4,'smiley':5,'button':32}
color={1:'green',2:'yellow',3:'red',4:'c',5:'black',32:'blue'}
mne.viz.plot_events(events,raw.info['sfreq'],raw.first_samp,color=color,event_id=event_id)
raw.plotにeventsという引数があれば、トリガー時が水色となります。
raw.plot(events=events,n_channels=10,order=order)
エポックの期間を-0.2〜0.3秒とし、-0.2〜0.0秒を基線とし、magnetometer＞4×10-12、眼電図＞200×10-6を除外した結果を表示します。
tmin,tmax=-0.2,0.5
event_id={'Auditory/Left':1,'Auditory/Right':2,'Visual/Left':3,'Visual/Right':4}
picks=mne.pick_types(raw.info,meg=True,eeg=False,eog=True)
baseline=(None,0.0)
reject={'mag':4e-12,'eog':200e-6}
epochs=mne.Epochs(raw,events=events,event_id=event_id,tmin=tmin,tmax=tmax,reject=reject,picks=picks,add_eeg_ref=False)
epochs.plot()
加算されなかったエポックの統計を表示します。
epochs.plot_drop_log()
加算平均し、波形を表示します。
picks=mne.pick_types(epochs.info,meg=True,eog=False)
evoked_left=epochs['Auditory/Left'].average(picks=picks)
evoked_right=epochs['Auditory/Right'].average(picks=picks)
evoked_left.plot()
evoked_right.plot()
ウインドウを閉じます。
import
matplotlib.pyplot as plt
plt.close('all')
事象関連脱電位とトポグラフィ
いろいろ読み込みます。
import
mne
from
mne.datasets import sample
data_path=sample.data_path()
raw_fname=data_path+'\\MEG\\sample\\sample_audvis_filt-0-40_raw.fif'
event_fname=data_path+'\\MEG\\sample\\sample_audvis_filt-0-40_raw-eve.fif'
raw=mne.io.read_raw_fif(raw_fname,add_eeg_ref=False,preload=True)
raw.set_eeg_reference()
# set EEG average reference
脳波の情報を見てみます。
raw.pick_types(meg=False,eeg=True,eog=True)
print(raw.info)
<Raw |
sample_audvis_filt-0-40_raw.fif, n_channels x n_times : 60 x 41700
(277.7 sec), ~22.1 MB, data loaded>
OSがWindows以外だともう少し情報が出るようです。
眼電図を脳波ということにして情報を見てみます。
raw.set_channel_types(mapping={'EOG
061': 'eeg'})
#
print(raw.info) #Windowsでは不要みたい
<Info
| 20 non-empty fields
bads : 'list | 0 items
buffer_size_sec :
'numpy.float64 | 13.3196808772
ch_names : 'list | EEG 001,
EEG 002, EEG 003, EEG 004, EEG 005, EEG 006, ...
chs : 'list | 60 items (EEG:
60)
comps : 'list | 0 items
custom_ref_applied : 'bool |
False
dev_head_t :
'mne.transforms.Transform | 3 items
dig : 'list | 146 items
events : 'list | 0 items
file_id : 'dict | 4 items
filename : 'str |
C:\Users\a.../sample_audvis_filt-0-40_raw.fif
highpass : 'float |
0.10000000149011612 Hz
hpi_meas : 'list | 1 items
hpi_results : 'list | 1
items
lowpass : 'float | 40.0 Hz
meas_date : 'numpy.ndarray |
2002-12-04 04:01:10
meas_id : 'dict | 4 items
nchan : 'int | 60
projs : 'list | PCA-v1: off,
PCA-v2: off, PCA-v3: off, ...
sfreq : 'float |
150.15374755859375 Hz
acq_pars : 'NoneType
acq_stim : 'NoneType
ctf_head_t : 'NoneType
description : 'NoneType
dev_ctf_t : 'NoneType
experimenter : 'NoneType
hpi_subsystem : 'NoneType
kit_system_id : 'NoneType
line_freq : 'NoneType
proj_id : 'NoneType
proj_name : 'NoneType
subject_info : 'NoneType
xplotter_layout : 'NoneType
眼電図のチャンネル名をEOG 61という名前に変更します。
raw.rename_channels(mapping={'EOG
061':'EOG'})
元に戻します。
raw.set_channel_types(mapping={'EOG':'eog'})
脳波のチャンネル位置を見てみます。チャンネル情報で0番チャンネルのlocで指定します。
print(raw.info['chs'][0]['loc'])
[-0.03737009 0.10568011 0.07333875 ..., 0.
0.
1. ]
12個の数字が出てきます。FIFFなんでたぶん、座標、緯度方向の単位ベクトル、経度方向の単位ベクトル、法線方向の単位ベクトルです。
センサ配列の表示
センサ配列を鳥瞰図表示します。
raw.plot_sensors()
三次元表示します。
raw.plot_sensors('3d')
モンタージュの設定
モンタージュを表示させます。国際10-20法を表示させます。
montage=mne.channels.read_montage('standard_1020')
print(montage)
<Montage
| standard_1020 - 97 channels: LPA, RPA, Nz ...>
電位補正の取り消し
raw_no_ref,_=mne.io.set_eeg_reference(raw,[])
#電位補正eeg_referenceを無しにする
事象関連電位の計算
電位補正なしの状態で、脳波＞180×10-6、眼電図＞150×10-6のエポックを除外し、-0.2〜0.3秒の加算平均を行い、波形を表示させます。
reject=dict(eeg=180e-6,eog=150e-6)
event_id,tmin,tmax={'left/auditory':1},-0.2,0.5
events=mne.read_events(event_fname)
epochs_params=dict(events=events,event_id=event_id,tmin=tmin,tmax=tmax,reject=reject,add_eeg_ref=False)
evoked_no_ref=mne.Epochs(raw_no_ref,**epochs_params).average()
del
raw_no_ref # メモリ節約
title='EEG
Original reference'
evoked_no_ref.plot(titles=dict(eeg=title))
警告が出ますが無視します。
<stdin>:1: RuntimeWarning: An average reference projection was
already added. The data has been left untouched.
0.1秒のトポグラフィを表示させます。
evoked_no_ref.plot_topomap(times=[0.1],size=3.,title=title)
電位補正ありの状態で加算平均し、波形・トポグラフィを表示させます。初期状態ではこちらが適応されます。
EEG
001とEEG
002で電位補正してみます。
raw_custom,_=mne.io.set_eeg_reference(raw,['EEG
001','EEG 002'])
evoked_custom=mne.Epochs(raw_custom,**epochs_params).average()
del
raw_custom # メモリ節約
title='EEG Custom reference'
evoked_custom.plot(titles=dict(eeg=title))
evoked_custom.plot_topomap(times=[0.1],size=3.,title=title)
加算平均同士の加減
最初に４つの課題別の加算平均を行います。
event_id={'left/auditory':1,'right/auditory':2,'left/visual':3,'right/visual':4}
epochs_params=dict(events=events,event_id=event_id,tmin=tmin,tmax=tmax,reject=reject)
epochs=mne.Epochs(raw,**epochs_params)
print(epochs)
警告が出ますが無視します。
<stdin>:1: DeprecationWarning: add_eeg_ref defaults to True in
0.13, will default to False in 0.14, and will be removed in 0.15. We recommend
to use add_eeg_ref=False and set_eeg_reference() instead.
>>>
<Epochs |
n_events : 288 (good & bad), tmin : -0.199795213158 (s), tmax :
0.499488032896 (s), baseline : (None, 0), ~3.1 MB, data not loaded,
'left/auditory': 72, 'left/visual': 73,
'right/auditory': 73, 'right/visual': 70>
重みづけをつけて合計し、トポグラフィ表示させます。左聴覚誘発＋左視覚誘発―右聴覚誘発―右視覚誘発をみてみます。
left,right=epochs['left'].average(),epochs['right'].average()
mne.combine_evoked([left,-right],weights='equal').plot_joint()
４つの課題別加算平均は以下のようにもできますが・・・
aud_l=epochs['auditory','left'].average()
aud_r=epochs['auditory','right'].average()
vis_l=epochs['visual','left'].average()
vis_r=epochs['visual','right'].average()
all_evokeds=[aud_l,aud_r,vis_l,vis_r]
print(all_evokeds)
[<Evoked |
comment : 'auditory+left',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 185, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'auditory+right',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 174, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'visual+left',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 179, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'visual+right',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 185, n_channels x n_times :
59 x 106, ~3.1 MB>]
Pythonのlist機能を使えば以下の文で１行で表せます。
all_evokeds=[epochs[cond].average()
for cond in sorted(event_id.keys())]
print(all_evokeds)
[<Evoked |
comment : 'auditory+left',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 185, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'auditory+right',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 174, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'visual+left',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 179, n_channels x n_times :
59 x 106, ~3.1 MB>,
<Evoked |
comment : 'visual+right',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 185, n_channels x n_times :
59 x 106, ~3.1 MB>]
左聴覚1/4、右聴覚-1/4、左視覚1/4、右視覚-1/4の重みづけで合計した波形とトポグラフィを示します。
mne.combine_evoked(all_evokeds,weights=(0.25,-0.25,0.25,-0.25)).plot_joint()
加算平均データの保存
ファイルとして保存します。
grand_average=mne.grand_average(all_evokeds)
mne.write_evokeds('c:\\Users\\akira\\tmp-ave.fif',all_evokeds)
#Windows用に変更
辞書オブジェクトにします
all_evokeds=dict((cond,epochs[cond].average())
for cond in event_id)
print(all_evokeds['left/auditory'])
<Evoked |
comment : 'left/auditory',
kind
: average, time : [-0.199795, 0.499488], n_epochs : 56, n_channels x n_times :
59 x 106, ~3.1 MB>
この辞書オブジェクトを使うと条件別４課題の波形とトポグラフィの表示が簡単になります。
for
cond in all_evokeds:
all_evokeds[cond].plot_joint(title=cond)
ウインドウを全部閉じます。
import
matplotlib.pyplot as plt
plt.close('all')