MNE学习笔记（六）：Epoched data的可视化
参考文章：https://mne.tools/stable/auto_tutorials/epochs/20_visualize_epochs.html
准备工作
包括
- 导包
- 加载
这个部分之前已经有详细地说明了，如果有疑惑可以去看MNE学习笔记（三）：Epoch数据结构，这里就直接给出代码：
import os
import mne
from matplotlib import pyplot as plt
# 如果没有数据则用这个自动下载
# sample_data_folder = mne.datasets.sample.data_path()
# 已有数据，则直接加载即可
sample_data_folder = "D:\Data\MNE-sample-data"
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False).crop(tmax=120)
events = mne.find_events(raw, stim_channel='STI 014')
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
'visual/right': 4, 'face': 5, 'button': 32}
epochs = mne.Epochs(raw, events, tmin=-0.2, tmax=0.5, event_id=event_dict,
preload=True)
del raw
from matplotlib import pyplot as plt官方文档中是没有这个的，建议还是加上，不然画出的图像会闪退
输出为类似这样的即可：
将Epochs可视化为时间序列
对于每一个channel都用一个时间序列来表示，下面我们来绘制磁力仪通道中的心跳伪影（heartbeat artifacts in the magnetometer channels）
代码：
# 可视化为时间序列
catch_trials_and_buttonpresses = mne.pick_events(events, include=[5, 32])
epochs['face'].plot(events=catch_trials_and_buttonpresses, event_id=event_dict,
event_color=dict(button='red', face='blue'))
plt.show()
结果：
如果要一次查看所有sensors，则可以使用butterfly mode和group by selection
代码：
# 一次显示所有sensors
epochs['face'].plot(events=catch_trials_and_buttonpresses, event_id=event_dict,
event_color=dict(button='red', face='blue'),
group_by='selection', butterfly=True)
plt.show()
结果：(在上面的图的基础上还会出现下图)
通过epochs对象绘制projectors
从磁盘上装入心电图projectors，并将其应用到数据中（load ECG projectors from disk and apply them to the data）
代码：
ecg_proj_file = os.path.join(sample_data_folder, 'MEG', 'sample',
'sample_audvis_ecg-proj.fif')
ecg_projs = mne.read_proj(ecg_proj_file)
epochs.add_proj(ecg_projs)
epochs.apply_proj()
结果：
由于原来的三个empty-room magnetometer projectors继承自Raw文件，我们为每个sensor类型增加了两个ECG projectors，我们应该看到9个projector地形图:
代码：
# 绘制地形图
epochs.plot_projs_topomap(vlim='joint')
print(all(proj['active'] for proj in epochs.info['projs']))
plt.show()
结果：（print结果为true）
绘制sensor位置
代码：
# 绘制sensor位置
epochs.plot_sensors(kind='3d', ch_type='all')
epochs.plot_sensors(kind='topomap', ch_type='all')
plt.show()
结果：
绘制Epochs的功率谱
代码：
# 绘制Epochs的功率谱
epochs['auditory'].plot_psd(picks='eeg')
plt.show()
结果：
还可以使用plot_psd_topomap，将sensors之间的光谱估计值作为头皮地形图（scalp topography）绘制出来。
代码：
# 绘制scalp topography
epochs['visual/right'].plot_psd_topomap()
plt.show()
结果：
与plot_projs_topomap一样，plot_psd_topomap有一个vlim='joint'选项，用于在所有子图上共同固定颜色栏限制，以便更好地了解每个波段的相对大小。可以通过ch_type参数更改 channel 类型，如果想查看与默认值不同的频带，bands参数有一个元组的列表，每个元组包含一个 single frequency 和 subplot title ，或 lower / upper 频率限制和一个subplot title:
代码：
bands = [(10, '10 Hz'), (15, '15 Hz'), (20, '20 Hz'), (10, 20, '10-20 Hz')]
epochs['visual/right'].plot_psd_topomap(bands=bands,
vlim='joint', ch_type='grad')
plt.show()
结果：
如果喜欢untransformed power estimates，也可以设置dB=False；通过normalize=True可以通过除以所有频率的总功率，使功率估计归一化。
图像映射（image map）
这部分我就直接翻译MNE官方教程文档了
同时可视化多个Epochs的一种方便的方法是将它们绘制成图像映射，图像中的每一行像素代表一个历元，横轴代表时间，每个像素的颜色代表该历元在该时间样本中的信号值。
当然，这需要为每个通道提供一个单独的图像映射，或者以某种方式组合跨通道的信息。 后者可以使用plot_image方法; 前者可以通过plot_image方法(一次一个通道)或plot_topo_image方法(一次所有传感器)来实现。
默认情况下，由plot_image生成的图像映射将伴随着一个指示颜色映射范围的标尺，一个时间序列，显示Epochs间的平均信号，以及一个围绕平均值的引导95%置信区间。
Plot_image是一个具有许多参数的高度可定制的方法，包括辅助颜色条和平均时间序列子图的定制。 参见plot_image和mne.viz的文档字符串。
Plot_compare_evokeds(用于绘制平均时间序列)以获得完整的细节。 在这里，我们将展示在有听觉刺激的情况下所有时期的磁强计平均值:
代码：
# 有听觉刺激的情况下所有时期的磁强计平均值
epochs['auditory'].plot_image(picks='mag', combine='mean')
plt.show()
结果：
若要绘制单个sensors或一小组sensors的image map，请使用picks参数。
传入combine=None(默认值)将为picks中的每个sensor生成单独的plots; 通过combine='gfp'将绘制出全局场功率(用于组合极性相反的sensors)。
代码：
# combine = 'gfp'
epochs['auditory'].plot_image(picks=['MEG 0242', 'MEG 0243'])
epochs['auditory'].plot_image(picks=['MEG 0242', 'MEG 0243'], combine='gfp')
plt.show()
结果：
为了给所有sonsers绘制image map，可以使用对同时绘制大量image maps进行优化过且允许你对每一个小image map进行点击从而生成一个单独的图与全尺寸图像图的plot_topo_image。
这部分中显示的小范围中，很难在这些图中看到很多有用的细节; 在交互式绘图时，最好将拓扑图最大化为全屏。默认情况下是一个带有黑色背景的图形，所以这里我们指定一个白色背景和黑色前景文本。 默认情况下，plot_topo_image将在同一plot上显示磁力计和梯度计(magnetometers and gradiometers)(不显示颜色条，因为传感器的比例不同)，所以我们还将传递一个Layout，将每个plot限制为一种channel类型。 然而，我们得先删除任何具有异常高信号水平的epochs，因为它们可能导致颜色图限制太极端，因此掩盖了我们感兴趣的较小信号波动。
代码：
# 交互图
reject_criteria = dict(mag=3000e-15, # 3000 fT
grad=3000e-13, # 3000 fT/cm
eeg=150e-6) # 150 µV
epochs.drop_bad(reject=reject_criteria)
for ch_type, title in dict(mag='Magnetometers', grad='Gradiometers').items():
layout = mne.channels.find_layout(epochs.info, ch_type=ch_type)
epochs['auditory/left'].plot_topo_image(layout=layout, fig_facecolor='w',
font_color='k', title=title)
plt.show()
结果：
点击其中一个后：
点击其中一个后：
这里只截取了部分
要绘制所有EEG sensors的image map，需要将EEG layout作为plot_topo_image的layout参数。 还要注意这里对每个image map在竖直维度进行了平滑(across epochs)的sigma参数的使用，它可以让我们更容易看见small image maps的patterns(通过smear嘈杂的epochs到他们的邻居,同时加强相邻epochs相似的部分图像)。 然而，sigma也可以掩盖有持续的极端值的时代，也许应该被排除，所以应该谨慎使用。
代码：
结果：
单击其中一种后：
完整代码
import os
import mne
from matplotlib import pyplot as plt
# 如果没有数据则用这个自动下载
# sample_data_folder = mne.datasets.sample.data_path()
# 已有数据，则直接加载即可
sample_data_folder = "D:\Data\MNE-sample-data"
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False).crop(tmax=120)
events = mne.find_events(raw, stim_channel='STI 014')
event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
'visual/right': 4, 'face': 5, 'button': 32}
epochs = mne.Epochs(raw, events, tmin=-0.2, tmax=0.5, event_id=event_dict,
preload=True)
del raw
# 可视化为时间序列
catch_trials_and_buttonpresses = mne.pick_events(events, include=[5, 32])
epochs['face'].plot(events=catch_trials_and_buttonpresses, event_id=event_dict,
event_color=dict(button='red', face='blue'))
plt.show()
# 一次显示所有sensors
epochs['face'].plot(events=catch_trials_and_buttonpresses, event_id=event_dict,
event_color=dict(button='red', face='blue'),
group_by='selection', butterfly=True)
plt.show()
# projectors
ecg_proj_file = os.path.join(sample_data_folder, 'MEG', 'sample',
'sample_audvis_ecg-proj.fif')
ecg_projs = mne.read_proj(ecg_proj_file)
epochs.add_proj(ecg_projs)
epochs.apply_proj()
# 绘制地形图
epochs.plot_projs_topomap(vlim='joint')
print(all(proj['active'] for proj in epochs.info['projs']))
plt.show()
# 绘制sensor位置
epochs.plot_sensors(kind='3d', ch_type='all')
epochs.plot_sensors(kind='topomap', ch_type='all')
plt.show()
# 绘制Epochs的功率谱
epochs['auditory'].plot_psd(picks='eeg')
plt.show()
# 绘制scalp topography
epochs['visual/right'].plot_psd_topomap()
plt.show()
bands = [(10, '10 Hz'), (15, '15 Hz'), (20, '20 Hz'), (10, 20, '10-20 Hz')]
epochs['visual/right'].plot_psd_topomap(bands=bands,
vlim='joint', ch_type='grad')
plt.show()
# 有听觉刺激的情况下所有时期的磁强计平均值
epochs['auditory'].plot_image(picks='mag', combine='mean')
plt.show()
# combine = 'gfp'
epochs['auditory'].plot_image(picks=['MEG 0242', 'MEG 0243'])
epochs['auditory'].plot_image(picks=['MEG 0242', 'MEG 0243'], combine='gfp')
plt.show()
# 交互图
reject_criteria = dict(mag=3000e-15, # 3000 fT
grad=3000e-13, # 3000 fT/cm
eeg=150e-6) # 150 µV
epochs.drop_bad(reject=reject_criteria)
for ch_type, title in dict(mag='Magnetometers', grad='Gradiometers').items():
layout = mne.channels.find_layout(epochs.info, ch_type=ch_type)
epochs['auditory/left'].plot_topo_image(layout=layout, fig_facecolor='w',
font_color='k', title=title)
plt.show()
layout = mne.channels.find_layout(epochs.info, ch_type='eeg')
epochs['auditory/left'].plot_topo_image(layout=layout, fig_facecolor='w',
font_color='k', sigma=1)
plt.show()
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
