转载自公众号：路同学
作者：路同学
Hello,
这里是行上行下，我是喵君姐姐~
在之前所发布的“MNE-Python的简易中文教程简单入门”中，我们介绍了使用MNE来对EEG/MEG进行预处理；
另外，在发布的“MNE进阶教程|实例详解EEG/MEG数据读取”教程中，主要介绍了以一个公开数据为案例，详解基于MNE-Python的从“原数据”到“可处理数据”的读取流程。
随后，发布的“MNE-Python简易教程|解析EEG/MEG数据中的事件信息”教程中，详解如何从EEG/MEG数据中解析出事件信息，即介绍MNE-Python中的事件信息的读取。
以及，在发布的“MNE-Python简易教程|EEG/MEG数据的伪影识别、修复坏道以及去除坏段”中，详细的介绍了数据预处理的核心阶段：伪影的处理。
接下来，我们详解利用MNE-Python进行数据滤波与降采样。
数据滤波
首先还是导入example数据方便后续说明：
import numpy as npimport mnefrom mne.datasets import sampledata_path = sample.data_path()raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'# 使用0-20ms的数据tmin, tmax = 0, 20raw = mne.io.read_raw_fif(raw_fname)raw.crop(tmin, tmax).load_data()# 计算时排除两个坏道raw.info['bads'] = ['MEG 2443', 'EEG 053']# 设置频率范围：2-300Hzfmin, fmax = 2, 300# 选择右侧颞叶的导联selection = mne.read_selection('Left-temporal')# 挑选导联：仅MEG 且 去除坏道 且 右侧颞叶导联picks = mne.pick_types(raw.info, meg='mag', eeg=False, eog=False, stim=False, exclude='bads', selection=selection)
可以画出功率谱图如下：
raw.plot_psd(area_mode='range', picks=picks, average=False)
使用陷波滤波器去噪：
之前介绍过电力线噪音，数据在60Hz、120Hz、180Hz和240Hz存在窄频率峰值，我们使用notch_filter()进行陷波滤波器对数据进行滤波：
raw.notch_filter(np.arange(60, 241, 60), picks=picks)raw.plot_psd(area_mode='range', picks=picks, average=False)
使用低通滤波器去噪：
使用滤波方法filter()进行50Hz的低通滤波：
raw.filter(None, 50.)raw.plot_psd(area_mode='range', picks=picks, average=False)
使用高通滤波器去噪：
使用滤波方法filter()进行1Hz的高通滤波：(EEG通常使用0.1Hz高通滤波)
raw.filter(1., None)raw.plot_psd(area_mode='range', picks=picks, average=False)
使用带通滤波器去噪：
直接使用一个1-50Hz的带通滤波实现上述两步的同样过程：
raw.filter(1., 50.)
数据降采样
使用resample()进行降采样：(将数据采样率降为100Hz)
# 这里传入的参数是新的采样频率：100Hzraw.resample(100, npad="auto")raw.plot_psd(area_mode='range', picks=picks, average=False)
resample()不仅仅在Raw对象里有，Epochs对象和Evoked对象下均有resample()方法。
接下来，我们详解利用MNE-Python进行独立成分分析(ICA)。
mne.preprocessing.ICA类
为了方便对于MNE-Python实现ICA部分的理解，首先要强调一下利用MNE-Python实现ICA的思路。
MNE-Python实现独立成分分析(ICA)的功能主要通过ICA这个类来实现，在进行ICA的时候思路类似使用scikit-learn一样，首先初始化一个ICA对象，然后用这个ICA对象去fit数据，得到一个“训练后“的ICA对象，之后可以调用ICA对象的各种方法进行数据重建、画图等等功能。
此外，在MNE-Python中不仅仅对Raw对象可以进行ICA，对Epochs和Evoked对象都可以进行ICA。
ICA的实现
首先还是导入example数据：
import osimport mnefrom mne.preprocessing import ICA# 导入数据sample_data_folder = mne.datasets.sample.data_path()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)# 1Hz的高通滤波filt_raw = raw.copy()filt_raw.load_data()filt_raw.filter(1., None)
ICA：
ica = ICA(n_components=15)ica.fit(filt_raw)
现在数据已经分成了15个独立成分了，我们看看这些成分：
各成分的波动图：
ica.plot_sources(filt_raw)
各成分的地形图：
ica.plot_components()
我们可以看看去除第一个成分前后信号的差异：
ica.plot_overlay(filt_raw, exclude=[0])
也可以查看单个成分的情况：
ica.plot_properties(filt_raw, picks=[0])
当然MNE-Python的ICA类中提供了三种计算ICA的算法，可以通过
method构造参数进行修改选择，具体可以参考：
https://mne.tools/stable/generated/mne.preprocessing.ICA.html#mne.preprocessing.ICA.plot_scores
去除ICA的成分
MNE-Python中通过ica.exclude属性来设置要去除的ICA成分，再用ica.apply()完成去除这个步骤。我们以去除前两个成分为例：
ica.exclude = [0, 1]reconst_raw = filt_raw.copy()ica.apply(reconst_raw)
这样就已经对数据去除了前两个ICA的成分，可以画出数据信号图对比一下取出前后的信号差异
filt_raw.plot()reconst_raw.plot()
MNE-Python也提供了通过EOG、ECG或者参考导联来识别与其高相关的成分方便去除，三个方法分别是：
find_bads_ecg()
find_bads_eog()
find_bads_ref()
这里以用过眼动EOG导联进行ICA成分选择为示例：
ica.exclude = []ecg_indices, ecg_scores = ica.find_bads_ecg(raw, method='correlation')ica.exclude = ecg_indices
可以看到调用find_bads_eog()方法后返回的是成分序号和成分与眼动信号的相关系数，将前者传给exclude即可以，后续直接apply()就可以实现剔除。
使用模板匹配选择ICA成分
当我们处理多个被试的数据时，MNE-Python提供了一种基于模板匹配的方案，其核心思想就是：
跨被试的伪迹是存在一些模式上的相似的，我们可以通过一个被试中坏的ICA成分来作为模板和其他被试的ICA成分进行快速匹配，这样能更高效地选择出各个被试数据的不好的ICA成分。
为了更好地用代码说明，这里先导入一个包含多个被试的example数据并对四个被试的数据跑ICA：
import mnefrom mne.preprocessing import ICA, corrmapmapping = { 'Fc5.': 'FC5', 'Fc3.': 'FC3', 'Fc1.': 'FC1', 'Fcz.': 'FCz', 'Fc2.': 'FC2', 'Fc4.': 'FC4', 'Fc6.': 'FC6', 'C5..': 'C5', 'C3..': 'C3', 'C1..': 'C1', 'Cz..': 'Cz', 'C2..': 'C2', 'C4..': 'C4', 'C6..': 'C6', 'Cp5.': 'CP5', 'Cp3.': 'CP3', 'Cp1.': 'CP1', 'Cpz.': 'CPz', 'Cp2.': 'CP2', 'Cp4.': 'CP4', 'Cp6.': 'CP6', 'Fp1.': 'Fp1', 'Fpz.': 'Fpz', 'Fp2.': 'Fp2', 'Af7.': 'AF7', 'Af3.': 'AF3', 'Afz.': 'AFz', 'Af4.': 'AF4', 'Af8.': 'AF8', 'F7..': 'F7', 'F5..': 'F5', 'F3..': 'F3', 'F1..': 'F1', 'Fz..': 'Fz', 'F2..': 'F2', 'F4..': 'F4', 'F6..': 'F6', 'F8..': 'F8', 'Ft7.': 'FT7', 'Ft8.': 'FT8', 'T7..': 'T7', 'T8..': 'T8', 'T9..': 'T9', 'T10.': 'T10', 'Tp7.': 'TP7', 'Tp8.': 'TP8', 'P7..': 'P7', 'P5..': 'P5', 'P3..': 'P3', 'P1..': 'P1', 'Pz..': 'Pz', 'P2..': 'P2', 'P4..': 'P4', 'P6..': 'P6', 'P8..': 'P8', 'Po7.': 'PO7', 'Po3.': 'PO3', 'Poz.': 'POz', 'Po4.': 'PO4', 'Po8.': 'PO8', 'O1..': 'O1', 'Oz..': 'Oz', 'O2..': 'O2', 'Iz..': 'Iz'}raws = list()icas = list()for subj in range(4): fname = mne.datasets.eegbci.load_data(subj + 1, runs=[3])[0] raw = mne.io.read_raw_edf(fname) # 传入mapping信息 raw.rename_channels(mapping) raw.set_montage('standard_1005') # ICA ica = ICA(n_components=30, random_state=97) ica.fit(raw) raws.append(raw) icas.append(ica)
MNE-Python中进行模板匹配的函数是mne.preprocessing.corrmap()，这里我们选择第一个被试(subject 0)作为我们选择ICA坏成分的模板：
raw = raws[0]ica = icas[0]eog_inds, eog_scores = ica.find_bads_eog(raw, ch_name='Fpz')corrmap(icas, template=(0, eog_inds[0]))
可以看到对subject 0利用find_bads_eog()方法找到了不好的成分ICA000，然后用这个成分与其他被试的ICA成分匹配找到了subject 2和subject 3各有一个类似的坏的ICA成分。
我们可以看看四个被试的ICA成分信号图：
for index, (ica, raw) in enumerate(zip(icas, raws)): fig = ica.plot_sources(raw) fig.suptitle('Subject {}'.format(index))
可以看到确实subject 1没有和subject 0的ICA000相似的成分，不过subject 3的ICA000和ICA002两个成分似乎都与subject 0的ICA000很相似，但我们刚才使用模板只识别出来了一个。
这里我们设置一下匹配阈值来放低匹配标准：
corrmap(icas, template=(0, eog_inds[0]), threshold=0.9)
我们可以把通过模板找到的成分设置一个标签为“blink“，然后我们看看匹配出来的信息：
corrmap(icas, template=(0, eog_inds[0]), threshold=0.9, label='blink', plot=False)print([ica.labels_ for ica in icas])
‘eog/0/Fpz’和‘eog’信息是由find_bads_eog()得到的，这里我们只看‘blinks’的信息，可以清楚地看到通过模板匹配的方式找到的四个被试的ICA成分是哪几个。
再看看subject 3找到的ICA坏的成分的情况：
icas[3].plot_components(picks=icas[3].labels_['blink'])
icas[3].exclude = icas[3].labels_['blink']icas[3].plot_sources(raws[3])
一点额外需要知道的：
corrmap()函数中的template不仅可以传入类似上面这种(0, eog_inds[0])被试序号和成分号的元祖，也可以直接传入一个ICA成分拓扑矩阵(即一个[n_channels, n_components]的numpy array矩阵。
一个实例：
template_eog_component = icas[0].get_components()[:, eog_inds[0]]corrmap(icas, template=template_eog_component, threshold=0.9)
可以看到与template=(0, eog_inds[0])是一样的结果。
PS：后台回复关键词“MNE数据滤波”，既可获得所述的文字版介绍及相关资料啦！
转载自公众号：路同学
作者：路同学
排版：shirly
MNE-Python的简易中文教程简单入门
MNE 进阶教程|实例详解EEG/MEG数据读取
MNE-Python教程|解析EEG/MEG数据的事件信息
MNE-Python教程|EEG/MEG数据伪影识别、修复
转播
