关注微信公众号：脑机接口研习社
了解脑机接口最近进展

系列文章目录

Python专栏 | 脑电图和脑磁图(EEG/MEG)的数据分析方法之载入数据

Python专栏 | MNE脑电数据(EEG/MEG)可视化

Python专栏 | MNE数据预处理方法——独立成分分析

Python专栏 | 独立成分分析(ICA)的实例应用：消除伪影信号

持续更新中……

---

文章目录

• 系列文章目录
• 识别伪影信号的方法
• • 1. Selecting ICA components manually （手动识别）
  • 2. Using an EOG channel to select ICA components
• 总结

---

识别伪影信号的方法

1. Selecting ICA components manually （手动识别）

确定要排除的Components后，可以通过设置ica.exclude属性手动指定。

仅设置ica.exclude不会立即执行任何操作，只会将排除的Independent Components (IC)添加到列表中，以便之后在需要时使用。

一旦设置了排除项，即使未传递exclude参数，比如plot_overlay的ICA方法，也会排除指定的Components。

使用mne.preprocessing.ICA.save和mne.preprocessing.read_ica时可以保留被排除的Components列表。

代码示例：

ica.exclude = [0, 1] #indices chosen based on various plots above

现在我们已经设置好了排除项（ica.exclude），接着我们可以使用ica.apply的方法重塑已经去除伪影的传感器信号。

即现在进入了下图中的ica.apply环节。

▲点击图片可查看大图

我们将把原始数据与重塑的数据一起绘制，来表明心跳和眨眼伪影已得到修复。

代码示例：

#ica.apply() changes the Raw object in-place, so let's make a copy first: raw.load_data() reconst_raw = raw.copy() ica.apply(reconst_raw)#原始数据 raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑数据，和原始数据进行对比，来表明心跳和眨眼伪影已得到修复。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #删除变量reconst_raw

输出结果：

Applying ICA to Raw instanceTransforming to ICA space (15 components)Zeroing out 2 ICA componentsProjecting back using 364 PCA components

2. Using an EOG channel to select ICA components

通过上面手动选择要排除的IC的方法似乎比较简单，但是当需要处理成百上千个被试的数据时，手动排除的方式就略显繁琐。

在这种情况下，比较合理的做法是：首先使用专用的EOG或ECG传感器进行实验，把由此测得的信号作为一种“模型（pattern）”，然后将“模型”与每一种IC相对应，利用算法自动标记出和“EOG / ECG模型信号”匹配的所有ICs。

在mne-python包里，我们将使用find_bads_eog自动查找与EOG信号最匹配的ICs，然后再使用plot_scores以及其他一些绘图功能来查看选择了哪些IC。

由于在手动选择步骤里，我们使用了ica.exclude = [0, 1]，手动选择了ICA001和ICA000，因此在使用自动选择功能函数前，我们首先需要将ica.exclude重置为空列表。

代码示例：

#自动选择与ECG和EOG信号对应的ICs ica.exclude = []#首先将ica.exclude变成空列表#find which ICs match the EOG pattern eog_indices, eog_scores = ica.find_bads_eog(raw) ica.exclude = eog_indices# barplot of ICA component "EOG match" scores ica.plot_scores(eog_scores) #eog_scores是计算每个ICs（本例子里共有15个ICs）与EOG信号的match scores。绝对值越接近1代表匹配度越高，越接近0代表匹配度越低。 #由画出的条形图可以看出，ICA000的值是0.94763376，匹配度最高。# plot diagnostics ica.plot_properties(raw, picks=eog_indices)# plot ICs applied to raw data, with EOG matches highlighted ica.plot_sources(raw, show_scrollbars=False)# plot ICs applied to the averaged EOG epochs, with EOG matches highlighted ica.plot_sources(eog_evoked)

输出结果：

filtering target Setting up band-pass filter from 1 - 10 HzFIR filter parameters --------------------- Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter: - Windowed frequency-domain design (firwin2) method - Hann window - Lower passband edge: 1.00 - Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz) - Upper passband edge: 10.00 Hz - Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 10.25 Hz) - Filter length: 6007 samples (10.001 sec)C:\ProgramData\Anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py:331: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.output = mkl_fft.rfft_numpy(a, n=n, axis=axis)Using multitaper spectrum estimation with 7 DPSS windows Not setting metadata Not setting metadata 138 matching events found No baseline correction applied 0 projection items activated 0 bad epochs dropped C:\ProgramData\Anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py:331: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.output = mkl_fft.rfft_numpy(a, n=n, axis=axis)Creating RawArray with float64 data, n_channels=16, n_times=166800Range : 25800 ... 192599 = 42.956 ... 320.670 secs Ready.

以上我们在原始Raw实例以及提取了EOG伪影的Evoked实例上都使用了plot_sources。这可以是另一种确认find_bads_eog已标识正确component的方法。

总结

今天用到的代码总结：

#手动选择要排除的ICs ica.exclude = [0, 1] # indices chosen based on various plots above#现在我们已经设置好了排除项（ica.exclude） #接着我们可以使用ica.apply的方法重塑已经去除伪影的传感器信号。 #ica.apply() changes the Raw object in-place, so let's make a copy first: raw.load_data() reconst_raw = raw.copy() ica.apply(reconst_raw)#原始数据 raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑数据，和原始数据进行对比，来表明心跳和眨眼伪影已得到修复。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #删除变量reconst_raw#自动选择与ECG和EOG信号对应的ICs ica.exclude = []#首先将ica.exclude变成空列表#find which ICs match the EOG pattern eog_indices, eog_scores = ica.find_bads_eog(raw) ica.exclude = eog_indices# barplot of ICA component "EOG match" scores ica.plot_scores(eog_scores) #eog_scores是计算每个ICs（本例子里共有15个ICs）与EOG信号的match scores。绝对值越接近1代表匹配度越高，越接近0代表匹配度越低。 #由画出的条形图可以看出，ICA000的值是0.94763376，匹配度最高。# plot diagnostics ica.plot_properties(raw, picks=eog_indices)# plot ICs applied to raw data, with EOG matches highlighted ica.plot_sources(raw, show_scrollbars=False)# plot ICs applied to the averaged EOG epochs, with EOG matches highlighted ica.plot_sources(eog_evoked)

参考链接

https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.html#fitting-and-plotting-the-ica-solution

未完待续……

欲知后事如何，请关注我们的公众号：脑机接口研习社

总结

以上是生活随笔为你收集整理的Python专栏 | ICA应用：如何识别伪影信号？（一）的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。