Python专栏 | ICA应用：如何识别伪影信号？（三）

最新推荐文章于 2021-11-09 09:33:48 发布

脑机接口研习社

最新推荐文章于 2021-11-09 09:33:48 发布

阅读量963

点赞数

分类专栏： Python与脑机接口文章标签： python 人工智能深度学习

原文链接：https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.html#fitting-and-plotting-the-ica-solution

版权

Python与脑机接口专栏收录该内容

9 篇文章 23 订阅

订阅专栏

在这里插入图片描述

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

系列文章目录

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

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

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

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

Python专栏 | ICA应用：如何识别伪影信号？（一）

Python专栏 | ICA应用：如何识别伪影信号？（二）

持续更新中……

文章目录

系列文章目录
Selecting ICA components using template matching
总结

Selecting ICA components using template matching

当我们需要处理多个被试数据时，也可以通过手动选择的方式选择要排除的IC，具体方法是：首先手动选取某一个被试数据中要排除的IC，然后使用mne.preprocessing.corrmap将该IC作为“模板（template）”，接着，根据此“模板”从其他被试的数据中选择要排除的IC。

corrmap背后的假设是，被试之间的伪影模式（artifacts pattern）足够相似，以致于符合条件的IC可以通过与公认的“模板”相对应，以及通过选择具有最高相关程度IC的方式来进行识别。

corrmap包含ICA solution列表和template参数，该参数将指定哪个ICA对象和其中的哪个components用作模板。

由于我们的样本数据集仅包含一个被试的数据，因此我们将使用包含多个被试的另一个数据集：EEGBCI数据集。

详细见：https://doi.org/10.1109/TBME.2004.827072；https://doi.org/10.1161/01.CIR.101.23.e215

该数据集包含109个被试，我们只需从前4个被试中下载一个任务（a left/right hand movement task，左/右手移动任务）。

代码示例：

#通过模板匹配template matching来选择ICs
mapping = {
    '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):
    #EEGBCI subjects are 1-indexed; run 3 is a left/right hand movement task
    #对4个被试的数据依次进行ica.fit，将结果依次保存在raws和icas的list里面
    fname = mne.datasets.eegbci.load_data(subj + 1, runs=[3])[0]
    raw = mne.io.read_raw_edf(fname)
    # remove trailing `.` from channel names so we can set montage
    raw.rename_channels(mapping)
    raw.set_montage('standard_1005')
    # fit ICA
    ica = ICA(n_components=30, random_state=97)
    ica.fit(raw)
    raws.append(raw)
    icas.append(ica)

Note：
上面的字典mapping里的值FC5, FC3, FC1……，指的是实验时EEG/MEG设备在被试脑袋上插的电极位置标记。

10-20系统电极放置法是国际脑电图学会规定的标准电极放置法。额极中点至鼻根的距离和枕点至枕外粗隆的距离各占此连线全长的10％，其余各点均以此连线全长的20％相隔，如下图所示。因此命名为10-20系统。

在这里插入图片描述
▲图片来源：https://www.diytdcs.com/2012/07/1020-system-electrode-distances/

详情请参考：https://my.oschina.net/u/989000/blog/4665357

输出结果：

Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R03.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Fitting ICA to data using 64 channels (please be patient, this may take a while)
Selecting by number: 30 components
Fitting ICA took 7.1s.
Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S002/S002R03.edf...
EDF file detected
Setting channel info structure...
Creating Ending factory farming structure...
Fitting ICA to data using 64 channels (please be patient, this may take a while)
Selecting by number: 30 components
Fitting ICA took 4.3s.
Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S003/S003R03.edf...
EDF file detected
Setting channel info structure...
Creating Ending factory farming structure...
Fitting ICA to data using 64 channels (please be patient, this may take a while)
Selecting by number: 30 components
Fitting ICA took 10.6s.
Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S004/S004R03.edf...
EDF file detected
Setting channel info structure...
Creating Ending factory farming structure...
Fitting ICA to data using 64 channels (please be patient, this may take a while)
Selecting by number: 30 components
Fitting ICA took 2.8s.

接着，让我们来运行corrmap:

代码示例：

# use the first subject as template; use Fpz as proxy for EOG
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]))

输出结果：
在这里插入图片描述

在这里插入图片描述

Using channel Fpz as EOG channel
... filtering ICA sources
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: 1600 samples (10.000 sec)

... filtering target
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: 1600 samples (10.000 sec)

Median correlation with constructed map: 0.983
Displaying selected ICs per subject.
No maps selected for subject [1], consider a more liberal threshold.

上面输出结果里的第一张图显示了模板图，而第二张图里显示了被认为与模板“匹配”的所有图（包括模板本身）。

由第二张图可以看到，四个被试只有三个match，没有subject[1]的原因可以在输出结果中看到：“No maps selected for subject [1], consider a more liberal threshold.（未选择被试1，请考虑更宽松的阈值）”

默认情况下，阈值是通过尝试多个值自动设置的。在这里，系统可能选择的阈值太高了。

接下来，让我们看一下每个被试的ICA资料：

代码示例：

for index, (ica, raw) in enumerate(zip(icas, raws)):
    fig = ica.plot_sources(raw, show_scrollbars=False)
    fig.subplots_adjust(top=0.9)  # make space for title
    fig.suptitle('Subject {}'.format(index))

输出结果：
在这里插入图片描述

Creating RawArray with float64 data, n_channels=30, n_times=20000
    Range : 0 ... 19999 =      0.000 ...   124.994 secs
Ready.
Creating RawArray with float64 data, n_channels=30, n_times=19680
    Range : 0 ... 19679 =      0.000 ...   122.994 secs
Ready.
Creating RawArray with float64 data, n_channels=30, n_times=20000
    Range : 0 ... 19999 =      0.000 ...   124.994 secs
Ready.
Creating RawArray with float64 data, n_channels=30, n_times=19680
    Range : 0 ... 19679 =      0.000 ...   122.994 secs
Ready.

由上图可看出，被试1确实具有看起来像反映眨眼伪影的IC（组件ICA000）。同时还要注意，被试3似乎具有反映眼部伪影的两个ICs（ICA000和ICA002），但是却只有一个被corrmap捕获。

所以，接下来让我们尝试手动设置阈值（threshold）：

代码示例：

corrmap(icas, template=(0, eog_inds[0]), threshold=0.9)

输出结果：
在这里插入图片描述

Median correlation with constructed map: 0.951
Displaying selected ICs per subject.
At least 1 IC detected for each subject.

现在，每个被试至少检测到了1个IC。

然后，我们将使用corrmap里的参数label =‘blink’，plot = False重新运行corrmap，以标记每个被试里捕获眨眼伪影的每个IC。（plot=False表示无需再次绘制图）。

代码示例：

corrmap(icas, template=(0, eog_inds[0]), threshold=0.9, label='blink',
        plot=False)
print([ica.labels_ for ica in icas])

输出结果：

Median correlation with constructed map: 0.951
At least 1 IC detected for each subject.
[{'eog/0/Fpz': [0], 'eog': [0], 'blink': [0]}, {'blink': [0]}, {'blink': [0]}, {'blink': [0, 2]}]

请注意上面的输出结果，第一个被试具有3个不同的IC标签：“eog / 0 / Fpz”，“eog”和“blink”。前两个由find_bads_eog添加， “blink”标签是在最后一次调用corrmap时添加的。

还要注意，被试3有两个ICs被标记为“blink”：0和2，这也与上面画出的图一致。

也可以手动编辑ICA对象的labels_属性，以使用自定义标签注释IC。

在绘图时自定义labels_很方便：

代码示例：

icas[3].plot_components(picks=icas[3].labels_['blink'])
icas[3].exclude = icas[3].labels_['blink']
icas[3].plot_sources(raws[3], show_scrollbars=False)

输出结果：

在这里插入图片描述

Creating RawArray with float64 data, n_channels=30, n_times=19680
    Range : 0 ... 19679 =      0.000 ...   122.994 secs
Ready.

最后，可以使用ICA对象的get_components方法以数字方式提取IC。这将返回一个可以传递给corrmap的NumPy数组，而不是我们之前传递的（subject_index，component_index）的元组，并且将产生相同的结果：

代码示例：

template_eog_component = icas[0].get_components()[:, eog_inds[0]]
corrmap(icas, template=template_eog_component, threshold=0.9)
print(template_eog_component)

输出结果：
在这里插入图片描述

Median correlation with constructed map: 0.951
Displaying selected ICs per subject.
At least 1 IC detected for each subject.
[-0.44495318 -0.40983842 -0.38930046 -0.37648581 -0.4081585  -0.43522276
 -0.47475364 -0.29240548 -0.28406613 -0.26882023 -0.25526935 -0.28871567
 -0.30326613 -0.30316262 -0.21016141 -0.22111437 -0.21807985 -0.207049
 -0.22375209 -0.24014921 -0.25321907 -1.68027772 -1.4824676  -1.62847163
 -1.40547681 -1.36613268 -0.85900872 -1.07936549 -1.52916767 -0.69233214
 -0.77286864 -0.65947637 -0.61444207 -0.59479909 -0.64739204 -0.65451653
 -0.76721442 -0.78455491 -0.38975992 -0.43540168 -0.24198818 -0.29395946
 -0.12909375 -0.16455719 -0.17778243 -0.18440583 -0.13619307 -0.16361049
 -0.17566436 -0.16126002 -0.16756067 -0.16994137 -0.18353601 -0.17522907
 -0.14012023 -0.09664031 -0.12390408 -0.12631576 -0.1284134  -0.09195589
 -0.09006861 -0.09628503 -0.06248746 -0.0248017 ]

总结

今天用到的代码总结：

#通过模板匹配template matching来选择ICs
mapping = {
    '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):
    #EEGBCI subjects are 1-indexed; run 3 is a left/right hand movement task
    #对4个被试的数据依次进行ica.fit，将结果依次保存在raws和icas的list里面
    fname = mne.datasets.eegbci.load_data(subj + 1, runs=[3])[0]
    raw = mne.io.read_raw_edf(fname)
    # remove trailing `.` from channel names so we can set montage
    raw.rename_channels(mapping)
    raw.set_montage('standard_1005')
    # fit ICA
    ica = ICA(n_components=30, random_state=97)
ica.fit(raw)
    raws.append(raw)
    icas.append(ica)
#接着运行corrmap
#use the first subject as template; use Fpz as proxy for EOG
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]))

#Let’s take a look at the ICA sources for each subject:
for index, (ica, raw) in enumerate(zip(icas, raws)):
    fig = ica.plot_sources(raw, show_scrollbars=False)
    fig.subplots_adjust(top=5)  # make space for title
    fig.suptitle('Subject {}'.format(index)) #设置每张图的标题
#Notice that subject 1 does seem to have an IC that looks like it reflects blink artifacts (component ICA000). 
#Notice also that subject 3 appears to have two components that are reflecting ocular artifacts (ICA000 and ICA002), but only one was caught by corrmap. 
#Let’s try setting the threshold manually:
corrmap(icas, template=(0, eog_inds[0]), threshold=0.9)

#At this point we’ll re-run corrmap with parameters label='blink', plot=False to label the ICs from each subject that capture the blink artifacts (without plotting them again).
corrmap(icas, template=(0, eog_inds[0]), threshold=0.9, label='blink',
        plot=False)
print([ica.labels_ for ica in icas])

#The labels_ attribute of ICA objects can also be manually edited to annotate the ICs with custom labels. They also come in handy when plotting:
icas[3].plot_components(picks=icas[3].labels_['blink'])
icas[3].exclude = icas[3].labels_['blink']
icas[3].plot_sources(raws[3], show_scrollbars=False)

#使用ICA对象的get_components方法以数字方式提取IC。
template_eog_component = icas[0].get_components()[:, eog_inds[0]]
corrmap(icas, template=template_eog_component, threshold=0.9)
print(template_eog_component)

参考链接
https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.html#fitting-and-plotting-the-ica-solution
https://my.oschina.net/u/989000/blog/4665357
图源/谷歌图片

在这里插入图片描述

脑机接口研习社

关注

0
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
Python专栏 | ICA应用：如何识别伪影信号？（三）

系列文章目录Python专栏 | 脑电图和脑磁图(EEG/MEG)的数据分析方法之载入数据Python专栏 | MNE脑电数据(EEG/MEG)可视化Python专栏 | MNE数据预处理方法——独立成分分析Python专栏 | 独立成分分析(ICA)的实例应用：消除伪影信号Python专栏 | ICA应用：如何识别伪影信号？（一）Python专栏 | ICA应用：如何识别伪影信号？（二）持续更新中……文章目录系列文章目录Selecting ICA components using te.
复制链接

扫一扫