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

系列文章目录

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

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

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

持续更新中……

文章目录

系列文章目录
举个例子
- 1. 载入数据（Loading data）
- 2. Visualizing the artifacts
- 3. Filtering to remove slow drifts
- 4. Fitting and plotting the ICA solution
总结

举个例子

独立成分分析（ICA）的一个应用例子是利用ICA消除伪影（artifacts）。

伪影是医学影像领域中的专业术语。伪影可以定义为图像中在视场（FOV）里歪曲物体的任何特征。所有MRI图像都包含一些伪影，它们的识别对于避免误诊和获取诊断图像很重要。被试的运动（例如呼吸，肠蠕动、心脏跳动、眼皮跳动等正常的生理活动）、设备、软件等因素均会引起伪影。

下面讲解如何利用ICA消除伪影。

1. 载入数据（Loading data）

关于在MNE里载入数据的代码解释可点击查看下面这篇文章：

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

代码示例：

import os
import mne
from mne.preprocessing
import (ICA, create_eog_epochs, create_ecg_epochs, corrmap)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)

输出结果：

Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
Read a total of 3 projection items:
PCA-v1 (1 x 102)  idle       PCA-v2 (1 x 102)  idle
PCA-v3 (1 x 102)  idle    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.

2. Visualizing the artifacts

让我们首先从可视化我们要修复的伪影开始。在已有数据集中，伪影是足够大的，可以很容易地在原始数据中看到：

代码示例：

#pick some channels that clearly show heartbeats and blinks
regexp = r'(MEG [12][45][123]1|EEG 00.)'artifact_picks = mne.pick_channels_regexp(raw.ch_names, regexp=regexp)raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False)

输出结果：

如上图所示，波段突起的部分（图上红色箭头部分）就是由于眼动或心脏跳动引起的运动性伪影（红色箭头为小编所加）。

下面利用create_eog_epochs显示每个通道中的眼动伪影（ocular artifact），即将上图中的突起波段进行突出显示。

代码示例：

eog_evoked = create_eog_epochs(raw).average()eog_evoked.apply_baseline(baseline=(None, -0.2))eog_evoked.plot_joint()

Electrocardiogram (ECG）心电图
Electrooculogram (EOG) 眼电图

详细代码请参考：
https://mne.tools/stable/auto_tutorials/preprocessing/plot_10_preprocessing_overview.html#tut-artifact-overview

输出结果：

EOG channel index for this subject is: [375]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Setting up band-pass filter from 1 - 10 Hz
FIR 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)
Now detecting blinks and generating corresponding events
Found 46 significant peaks
Number of EOG events detected : 46
Not setting metadata
Not setting metadata
46 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
Loading data for 46 events and 601 original time points ...
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)
0 bad epochs dropped
Applying baseline correction (mode: mean)
Created an SSP operator (subspace dimension = 3)
3 projection items activated
SSP projectors applied...
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
<Figure size 576x302.4 with 7 Axes>,<Figure size 576x302.4 with 7 Axes>,<Figure size 576x302.4 with 7 Axes>

同样的道理，我们利用create_ecg_epochs来可视化心脏跳动伪影。

代码示例：

ecg_evoked = create_ecg_epochs(raw).average()
ecg_evoked.apply_baseline(baseline=(None, -0.2))
ecg_evoked.plot_joint()

输出结果：

Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 8 - 16 Hz
FIR 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: 8.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 7.75 Hz)
- Upper passband edge: 16.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 16.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)
Number of ECG events detected : 283 (average pulse 61 / min.)
Not setting metadata
Not setting metadata
283 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
Loading data for 283 events and 601 original time points ...
0 bad epochs dropped
Applying baseline correction (mode: mean)
Created an SSP operator (subspace dimension = 3)
3 projection items activated
SSP projectors applied...
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>[<Figure size 576x302.4 with 7 Axes>,<Figure size 576x302.4 with 7 Axes>,<Figure size 576x302.4 with 7 Axes>]

3. Filtering to remove slow drifts

在运行ICA之前，比较重要的一个步骤是对数据进行滤波以消除低频漂移（low-frequency drifts），因为低频漂移可能会对ICA的拟合质量产生负面影响。

慢的漂移（slow drifts）是比较棘手的，因为它们会降低独立来源的独立性（independency）。

例如在上一篇提到的“鸡尾酒宴会问题”。详见：

MNE数据预处理方法——独立成分分析

我们想要还原的信号是下图 s，从图中可以看出信号波动明显，信号源非常具有独立性。

但是我们实际得到的信号是下图 x，中间由于存在许多slow drifts，所以信号源丧失了其独立性（红色箭头为小编所加）。

低频漂移的存在使得算法很难找到准确的ICA解决方案（ICA solution）。

这里建议使用截止频率为1Hz的低通滤波器。

由于滤波是线性操作，因此可以将已过滤信号中的ICA解决方案，应用于未过滤的信号。

更多信息请参见：https://doi.org/10.1109/EMBC.2015.7319296

所以我们将保留未过滤的Raw对象的副本，以便我们可以稍后对其应用ICA解决方案。

代码示例：

#过滤低频漂移
filt_raw = raw.copy() #创建副本
filt_raw.load_data().filter(l_freq=1., h_freq=None)

输出结果：

Reading 0 ... 166799  =      0.000 ...   277.714 secs...
Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 1 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1983 samples (3.302 sec)

4. Fitting and plotting the ICA solution

现在，我们准备设置并安装ICA。

▲点击图片可查看大图

既然我们通过观察原始数据，知道EOG和ECG的伪影相当显著，所以我们希望这些伪影在PCA分解的前几个维度中捕获。

我们可能不需要大量的分量（components）就可以很好地完成分解伪影的工作。通常情况下，最好n_components设置的大一些，会分解的较为精准。

n_components是需要不断尝试和测试，来最终确定到底设置成多少。

第一次尝试，我们将以n_components = 15来运行ICA。考虑到我们的Raw data里有超过300个以上的通道，故将n_components设置成15实际上偏小了，但是这样可以运行的比较快。

ICA拟合（ICA.fit）不是确定性的，例如，components在不同的运行过程中可能会出现符号翻转（sign flip），或者不一定总是以相同的顺序返回。因此我们还将指定一个随机种子（random seed），以便每次都能获得相同的结果。

代码示例：

ica = ICA(n_components=15, random_state=97)
ica.fit(filt_raw)

输出结果：

Fitting ICA to data using 364 channels (please be patient, this may take a while)
Selecting by number: 15 components
Fitting ICA took 8.9s.
<ICA | raw data decomposition, fit (fastica): 166800 samples, 15 components, channels used: "mag"; "grad"; "eeg">

接着，我们可以通过查看独立成分（Independent Components, ICs），来看它们捕获了什么。

plot_sources将显示ICs的时间序列。

代码示例：

#ICA结果的可视化呈现
filt_raw.load_data()
ica.plot_sources(filt_raw, show_scrollbars=False)

输出结果：

Creating RawArray with float64 data, n_channels=16, n_times=166800Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.

Note：
在对plot_sources的调用中，我们可以使用未经过滤的原始Raw对象。

代码示例：

raw.load_data()
ica.plot_sources(raw, show_scrollbars=False)

输出结果：

Creating RawArray with float64 data, n_channels=16, n_times=166800Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.

由上图可知，第一个分量（ICA000）很好地捕获了EOG 061信号，第二个组件（ICA001）看起来很像心跳。

P.S.有关在视觉上识别独立成分IC的更多信息，EEGLAB教程是一个很好的参考资料。https://labeling.ucsd.edu/tutorial/labels

我们还可以使用plot_components可视化每个组件的头皮区域（scalp field distribution）分布。

代码示例：

ica.plot_components()

输出结果：

在上面的图中，很明显地可以看出，哪些IC正在捕获被试的EOG和ECG伪影。

同时，也有其他的可视化方案。我们还可以使用plot_overlay绘制原始信号与排除了伪影的重构信号的对比图。

代码示例：

# blinks
ica.plot_overlay(raw, exclude=[0], picks='eeg')
#选择EEG，是因为EEG更容易监测到眨眼活动。exclude=[0]中的0是指上面的ICA000。# heartbeats
ica.plot_overlay(raw, exclude=[1], picks='mag')
#选择MAG，是因为MAG更容易监测到心脏跳动。exclude=[1]中的1是指上面的ICA001。

输出结果：

Applying ICA to Raw instanceTransforming to ICA space (15 components)Zeroing out 1 ICA componentProjecting back using 364 PCA componentsApplying ICA to Raw instanceTransforming to ICA space (15 components)Zeroing out 1 ICA componentProjecting back using 364 PCA components

我们还可以使用plot_properties绘制每个IC的一些诊断信息：

代码示例：

ica.plot_properties(raw, picks=[0, 1])
#选择ICA000和ICA001

输出结果：

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
Not setting metadata
Not setting metadata
138 matching events found
No baseline correction applied
0 projection items activated
0 bad epochs dropped

总结

今天用到的代码总结：

#载入数据
import os
import numpy as np
import mne
from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs,corrmap)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)#raw.crop(tmax=60.) #取60s的数据#可视化伪影
# pick some channels that clearly show heartbeats and blinks
regexp = r'(MEG [12][45][123]1|EEG 00.)'
artifact_picks = mne.pick_channels_regexp(raw.ch_names, regexp=regexp)
raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False)#利用create_eog_epochs显示每个通道中的眼动伪影（ocular artifact），
#详细代码参考：https://mne.tools/stable/auto_tutorials/preprocessing/plot_10_preprocessing_overview.html#tut-artifact-overview
eog_evoked = create_eog_epochs(raw).average()
eog_evoked.apply_baseline(baseline=(None, -0.2))
eog_evoked.plot_joint()
#Electrocardiogram (ECG）心电图；Electrooculogram (EOG) 眼电图#同样的道理，我们利用create_ecg_epochs来可视化心脏跳动伪影
ecg_evoked = create_ecg_epochs(raw).average()
ecg_evoked.apply_baseline(baseline=(None, -0.2))
ecg_evoked.plot_joint()#过滤低频漂移
filt_raw = raw.copy() #创建副本
filt_raw.load_data().filter(l_freq=1., h_freq=None)#ICA拟合
ica = ICA(n_components=15, random_state=97)
ica.fit(filt_raw)#ICA结果的可视化呈现
filt_raw.load_data()
ica.plot_sources(filt_raw, show_scrollbars=False)#Note：在对plot_sources的调用中，我们也可以将ica应用在未经过滤的原始Raw对象上。
raw.load_data()
ica.plot_sources(raw, show_scrollbars=False)#使用plot_components可视化每个组件的头皮区域（scalp field distribution）分布。
ica.plot_components()#可以使用plot_overlay绘制原始信号与排除了伪影的重构信号的对比图。
# blinks 眨眼伪影
ica.plot_overlay(raw, exclude=[0], picks='eeg')
#选择EEG，是因为EEG更容易监测到眨眼活动。exclude=[0]中的0是指上面的ICA000。# heartbeats 心动伪影
ica.plot_overlay(raw, exclude=[1], picks='mag')
#选择MAG，是因为MAG更容易监测到心脏跳动。exclude=[1]中的0是指上面的ICA001。#我们还可以使用plot_properties绘制每个IC的一些诊断信息：
ica.plot_properties(raw, picks=[0, 1])#选择ICA000和ICA001

未完待续……

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

参考链接

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

https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.html#visualizing-the-artifacts

https://zhuanlan.zhihu.com/p/104479669

图源/谷歌图片

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

独立成分分析ICA/FastICA
最近在学习CSP,然后又注意到了ICA,这个算法之前就用过,但是没有系统的整理一下,所以就在这里梳理一下相关内容,方便以后查阅. 独立成分分析ICA/FastICA 1 盲源分离(BlindSourc ...
学习笔记 | 独立成分分析(ICA, FastICA)及应用
学习笔记 | 独立成分分析[ICA, FastICA]及应用 1 背景说明 2 算法原理 2.1 ICA简介 2.2 形式化表达 3 算法步骤与代码 4 算法改进:FastICA 5 ICA实例与应用 ...
R语言使用caret包的preProcess函数进行数据预处理：对所有的数据列进行独立成分分析ICA（Independent components analysis）、设置method参数为ica
R语言使用caret包的preProcess函数进行数据预处理:对所有的数据列进行独立成分分析ICA(Independent components analysis).设置method参数为ica 目 ...
独立成分分析ICA、因子分析、LDA降维、NMF非负矩阵分解
独立成分分析ICA.因子分析.LDA降维.NMF非负矩阵分解目录独立成分分析ICA.因子分析.LDA降维.NMF非负矩阵分解独立成分分析ICA
独立成分分析ICA系列1：意义
1.前言独立成分分析思想和方法最早源于上世纪八十年代几个法国学者的研究工作,尽管当时他们并没有命名其为ICP;在1986年举行的神经网络计算会议上,法国学者Heraull和Jutten提出了一个基于 ...
独立成分分析 ( ICA ) 与主成分分析 ( PCA ) 的区别
1.前言参考资料:https://www.zhihu.com/question/28845451 书上写的是: 1. 主成分分析假设源信号间彼此非相关,独立成分分析假设源信号间彼此独立. 2. 主成 ...
lda 吗样本中心化需要_机器学习 —— 基础整理（四）：特征提取之线性方法——主成分分析PCA、独立成分分析ICA、线性判别分析LDA...
本文简单整理了以下内容: (一)维数灾难 (二)特征提取--线性方法 1. 主成分分析PCA 2. 独立成分分析ICA 3. 线性判别分析LDA (一)维数灾难(Curse of dimensiona ...
sklearn自学指南(part50)--独立成分分析(ICA)
学习笔记,仅供参考,有错必纠文章目录分解信号的分量(矩阵分解问题) 独立成分分析(ICA) 分解信号的分量(矩阵分解问题) 独立成分分析(ICA) 独立成分分析将一个多变量信号分离成最大程度独立的 ...
独立成分分析ICA在MATLAB中的实现
独立成分分析ICA在MATLAB中实现展示 by:YANG LIU 测试: 输入信号输出结果: 特别注意:代码28,29,30行".* "和" * "的区别. ...

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

系列文章目录

文章目录

举个例子

1. 载入数据（Loading data）

2. Visualizing the artifacts

3. Filtering to remove slow drifts

4. Fitting and plotting the ICA solution

总结

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

最新文章

热门文章