Python-EEG工具库MNE中文教程(11)-信号空间投影SSP 应用

@[toc]

本分享为脑机学习者Rose整理发表于公众号：脑机接口社区(微信号：Brain_Computer).QQ交流群：903290195

信号空间投影（SSP）

在前面一篇分享信号空间投影SSP数学原理中提到，投影矩阵将根据您试图投射出的噪声种类而变化。信号空间投影（SSP）是一种经过比较有无感兴趣信号的测量值来估算投影矩阵应该是什么的方法。例如，您能够进行其余“空房间”测量，以记录没有对象存在时传感器上的活动。经过查看空房间测量中各MEG传感器的活动空间模式，能够建立一个或多个N维向量，以给出传感器空间中环境噪声的“方向”（相似于上面示例中“触发器的影响”的向量）。SSP一般也用于消除心跳和眼睛运动伪影，在用于消除心跳和眼睛运动伪影的案例中，就不是经过空房间录制，而是经过检测伪影，提取伪影周围的时间段(epochs)并求平均值来估计噪声的方向。有关示例，请参见使用SSP修复工件。

一旦知道了噪声向量，就能够建立一个与其正交的超平面，并构造一个投影矩阵，将实验记录投影到该超平面上。这样，测量中与环境噪声相关的部分就能够被移除。一样，应该清楚的是，投影下降了数据的维数-你仍然会有相同数量的传感器信号，但它们不会都是线性独立的-但一般有数十或数百个传感器，而你要消除的噪声子空间只有3-5维，所以自由度的损失一般是没有问题的。

MNE Python中的投影(projector)

在示例数据中，已经使用空房间记录执行了SSP，可是投影与原始数据一块儿存储，而且还没有应用（或者说，投影还没有激活）。
在这里，我将加载示例数据并将其裁剪为60秒；

能够在如下read_raw_fif（）的输出中看到投影：

1.导入工具库

import os
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa
from scipy.linalg import svd
import mne

2.加载数据

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).load_data()

在MNE-Python中，环境噪声矢量是经过主成分分析（一般缩写为PCA）来计算的，这就是为何SSP投影仪一般有“PCA-v1”这样的名称。（顺便说一句，因为执行主成分分析的过程在幕后使用了奇异值分解，所以在已发表的论文中也常常看到相似“投影仪是使用SVD计算的”这样的短语。）投影仪存储在raw.info的projs字段中：

在MNE-Python中，使用主成分分析(一般缩写为"PCA")来计算环境噪声向量，这就是SSP投影一般使用"PCA-v1"之类的名称的缘由。 (顺便说一下，因为执行PCA的过程在后台使用了奇异值分解，因此在已发表的论文中常常会看到相似"使用SVD计算投影"之类的短语。） 投影(projector)存储在raw.info的projs字段中：

print(raw.info['projs'])

[<Projection | PCA-v1, active : False, n_channels : 102>, <Projection | PCA-v2, active : False, n_channels : 102>, <Projection | PCA-v3, active : False, n_channels : 102>]

raw.info['projs']是投影对象的普通Python列表，能够经过索引来访问各个投影。投影对象(Projection object)自己相似于Python dict，所以可使用其.keys()方法查看它包含哪些字段(一般不须要直接访问其属性，但若有必要，能够这样作)：

first_projector = raw.info['projs'][0]
print(first_projector)
print(first_projector.keys())

<Projection | PCA-v1, active : False, n_channels : 102> dict_keys(['kind', 'active', 'desc', 'data', 'explained_var'])

Raw，Epoch和Evoked对象都有一个布尔类型的 proj属性，该属性指示对象中是否存储有任何未应用/不活动的投影。换句话说，若是至少有一个投影而且全部投影都处于活动状态，则proj属性为True。此外，每一个投影还具备一个布尔活动字段：

print(raw.proj)
print(first_projector['active'])

False False

3.计算投影

在MNE Python中，SSP向量可使用如下通用函数计算： mne.compute_proj_raw； mne.compute_proj_epochs； mne.compute_proj_evoked。 这些函数所作的通常假设是，传递的数据包含要经过投影修复的工件的原始数据、时间段或平均值。 在实践中，这一般涉及空房间记录或平均ECG或EOG伪影的连续原始数据。

"""
经过比较使用和不使用投影的曲线图，能够看到投影仪对测量信号的影响。
默认状况下，`raw.plot（）`将在绘图前在后台应用投影仪(不修改：class：`mne.io.Raw`对象);
能够经过以下所示的布尔``proj``参数来控制它，
也能够经过绘图窗口右下角的：kbd：`Proj`按钮访问投影界面，
以交互方式打开和关闭它们.
这里咱们只看磁力计，还有一个文件开头的2秒样本。
"""
mags = raw.copy().crop(tmax=2).pick_types(meg='mag')
for proj in (False, True):
    fig = mags.plot(butterfly=True, proj=proj)
    fig.subplots_adjust(top=0.9)
    fig.suptitle('proj={}'.format(proj), size='xx-large', weight='bold')

如上图，未进行投影的数据proj=Flase的效果展现，点击"proj"红框弹出SSP projectors vectors能够发现都未选中激活；进行了投影的数据图为proj=True点击"proj"红框能够发现投影的信息。

4.加载和保存投影

SSP除了能够减小环境噪声外，还能够用于其余类型的信号清洗。能够发如今上一个图中的磁力计信号中有两个较大的偏移，这些偏移没有被空房间的投影消除，这是受试者心跳的伪影。SSP也能够用于移除这些工件。样本数据包括用于下降心跳噪声的投影，这些投影与原始数据保存在单独的文件中，可使用mne.read_proj()函数加载该文件：

ecg_proj_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                             'sample_audvis_ecg-proj.fif')
ecg_projs = mne.read_proj(ecg_proj_file)
print(ecg_projs)

"""
利用mne.write_proj()函数，可用于将投影数据以.fif格式保存到磁盘：

MNE Python推荐使用以-proj.fif(或-proj.fif.gz)来保存投影数据
"""
mne.write_proj('heartbeat-proj.fif', ecg_projs)

5.添加和移除投影

上面，当咱们打印从文件加载的ecg_projs列表时，它显示了两台用于梯度计的投影（前两台，标为"planar"），两台用于磁力计的投影（中两台，标为"axial"），两台用于EEG 传感器（最后两个，标记为"eeg"）。咱们可使用add_proj()方法将它们添加到Raw对象：

raw.add_proj(ecg_projs)

要删除投影，能够利用del_proj()方法，它是根据raw.info['projs']列表中的索引删除投影。 若是想要用新的投影替换现有投影，可使用raw.add_proj(ecg_projs,remove_existing = True)来实现。

想要了解心电图(ECG)投影仪如何影响测量的信号，咱们能够再次对使用投影和不使用投影的数据进行绘图（注：plot()方法只能临时应用投影进行可视化，而不会永久更改基础数据）。咱们将上面建立的mags变量（只有空房间SSP投影）与空房间和ECG投影仪的数据进行比较：

mags_ecg = raw.copy().crop(tmax=2).pick_types(meg='mag')
for data, title in zip([mags, mags_ecg], ['Without', 'With']):
    fig = data.plot(butterfly=True, proj=True)
    fig.subplots_adjust(top=0.9)
    fig.suptitle('{} ECG projector'.format(title), size='xx-large',
                 weight='bold')

在without ECG projector中，meg数据中的ECG部分没有进行projector，而在with ECG projector中，meg数据中ECG部分进行了projector，结果要平滑一些。

参考 Uusitalo MA and Ilmoniemi RJ. (1997). Signal-space projection method for separating MEG or EEG into components. Med Biol Eng Comput 35(2), 135–140. doi:10.1007/BF02534144

Python-EEG工具库MNE中文教程(11)-信号空间投影SSP 应用

脑机学习者Rose笔记分享,QQ交流群：903290195 更多分享，请关注公众号