参考脑机接口社区公众号文章《脑电分析系列[MNE-Python-11]| 信号空间投影SSP 应用》(以下简称《文章》)学习记录
首先感谢脑机接口社区公众号分享的《文章》,该博客仅作为本人在实际学习和理解该文章知识的记录和期间的疑惑、调试过程的记录。
《文章》使用的源码如下:
"""
信号空间投影(SSP)
在前面一篇分享(脑电分析系列[MNE-Python-10]|
信号空间投影SSP数学原理)中提到,
投影矩阵将根据您试图投射出的噪声种类而变化。
信号空间投影(SSP)
是一种通过比较有无感兴趣信号的测量值来估算投影矩阵应该是什么的方法。
例如,您可以进行其他“空房间”测量,以记录没有对象存在时传感器上的活动。
通过查看空房间测量中各MEG传感器的活动空间模式,可以创建一个或多个N维向量,
以给出传感器空间中环境噪声的“方向”(类似于上面示例中“触发器的影响”的向量)。
SSP通常也用于消除心跳和眼睛运动伪影,在用于消除心跳和眼睛运动伪影的案例中,
就不是通过空房间录制,而是通过检测伪影,
提取伪影周围的时间段(epochs)并求平均值来估计噪声的方向。
有关示例,请参见使用SSP修复工件。
一旦知道了噪声向量,就可以创建一个与其正交的超平面,
并构造一个投影矩阵,将实验记录投影到该超平面上。这样,
测量中与环境噪声相关的部分就可以被移除。同样,应该清楚的是,
投影降低了数据的维数-你仍然会有相同数量的传感器信号,
但它们不会都是线性独立的-但通常有数十或数百个传感器,
而你要消除的噪声子空间只有3-5维,因此自由度的损失通常是没有问题的。
"""
"""
在示例数据中,已经使用空房间记录执行了SSP,
但是投影与原始数据一起存储,并且尚未应用(或者说,投影尚未激活)。
在这里,我将加载示例数据并将其裁剪为60秒;
可以在以下read_raw_fif()的输出中看到投影
"""
import os
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import svd
import mne
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()
print(raw)
"""
在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'])
"""
raw.info['projs']是投影对象的普通Python列表,
可以通过索引来访问各个投影。投影对象(Projection object)
本身类似于Python dict,
因此可以使用其.keys()方法查看它包含哪些字段
(通常不需要直接访问其属性,但如有必要,可以这样做):
"""
first_projector = raw.info['projs'][0]
print(first_projector)
print(first_projector.keys())
"""
Raw,Epoch和Evoked对象都有一个布尔类型的 proj属性,
该属性指示对象中是否存储有任何未应用/不活动的投影。
换句话说,如果至少有一个投影并且所有投影都处于活动状态,
则proj属性为True。此外,每个投影还具有一个布尔活动字段:
"""
print(raw.proj)
print(first_projector['active']