专栏首页脑机接口Python-EEG工具库MNE中文教程(10)-信号空间投影SSP数学原理

Python-EEG工具库MNE中文教程(10)-信号空间投影SSP数学原理

projector(投影)和投影背景

projector(投影)(简称proj),也称为信号空间投影(SSP),定义了应用于空间上的EEG或MEG数据的线性操作。

可以将该操作看做是一个矩阵乘法,通过将数据投影到较低维度的子空间来降低数据的秩。

这种投影算子可以同时应用于数据和正向运算,用于源定位。注意,可以使用这样的投影算子来完成EEG平均参考。

它存储在info['projs']中的测量信息中。

下面就结合代码来解释投影原理

导入工具库

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

# 定义绘制3d图像
def setup_3d_axes():
    ax = plt.axes(projection='3d')
    ax.view_init(azim=-105, elev=20)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-1, 5)
    ax.set_ylim(-1, 5)
    ax.set_zlim(0, 5)
    return ax

什么是projector(投影)?

在最基本的术语中,投影是将一组点转换为另一组点的操作,在这些点上重复投影操作没有效果。

给一个简单的几何示例,请想象三维空间中的点(3,2,5)。

如果太阳在x,y平面正上方,则该点在x,y平面上的投影看起来很像该点投射的阴影:

ax = setup_3d_axes()

# 绘制向量(3, 2, 5)
origin = np.zeros((3, 1))
point = np.array([[3, 2, 5]]).T
vector = np.hstack([origin, point])
ax.plot(*vector, color='k')
ax.plot(*point, color='k', marker='o')

# 将向量投影到x,y平面上并将其绘制
xy_projection_matrix = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
projected_point = xy_projection_matrix @ point
projected_vector = xy_projection_matrix @ vector
ax.plot(*projected_vector, color='C0')
ax.plot(*projected_point, color='C0', marker='o')

# 添加显示投影的虚线箭头
arrow_coords = np.concatenate([point, projected_point - point]).flatten()
ax.quiver3D(*arrow_coords, length=0.96, arrow_length_ratio=0.1, color='C1',
            linewidth=1, linestyle='dashed')

注意,使用矩阵乘法来计算点(3,2,5)平面的投影:

=

再次将投影应用于结果,只会再次返回结果:

=

从信息的角度来看,此投影采用了点x,y,z,并删除了有关该点在z方向上位置的信息。现在我们所知道的只是它在x,y平面上的位置。此外,将投影矩阵应用于x,y,z空间中的任何点,都会将其缩小为x,y平面上的对应点。术语是子空间:投影矩阵将原始空间中的点投影到比原始空间低维的子空间中。我们的子空间是x,y平面(而不是y,z平面)的原因是投影矩阵中特定值的直接结果。

示例:投影作为降噪

描述这种“信息丢失”或“投影到子空间”的另一种方式是,投影将测量的秩(或“自由度”)从三维降低到二维。如果您知道z方向上的测量分量只是由于您的测量方法而产生的噪声,而您所关心的只是x和y分量,那么将三维测量投影到x,y平面中就可以看作是一种降噪的形式。

当然,如果所有测量噪声都集中在z方向上,那确实是非常幸运的。您可以直接舍弃z分量,而不必费心构造投影矩阵或进行矩阵乘法。假设为了进行该测量,您必须在测量设备上触动触发器(pull a trigger on a measurement device),并且触动触发器的动作会使设备移动一点。如果您测量触发器的触动是如何影响测量设备的位置,则可以"校正"实际测量值,以“投射”触动触发器的效果。在这里,我们假设触发器的平均效果是将测量设备移动(3,−1,1):

trigger_effect = np.array([[3, -1, 1]]).T

计算正交平面

知道了这一点,我们可以计算一个与触发器作用正交的平面(利用一个穿过原点的平面在给定法向量(A,B,C)的情况下具有方程Ax+By+Cz=0),并将实际测量投影到该平面上。

# 计算与trigger_effect正交的平面
x, y = np.meshgrid(np.linspace(-1, 5, 61), np.linspace(-1, 5, 61))
A, B, C = trigger_effect
z = (-A * x - B * y) / C
# 切断z=0以下的平面(只是为了使绘图更精细)
mask = np.where(z >= 0)
x = x[mask]
y = y[mask]
z = z[mask]

使用SVD计算投影矩阵

使用奇异值分解(SVD)从trigger_effect向量计算投影矩阵;

放置好投影矩阵后,我们可以投影原始向量(3,2,5)以消除触发器的影响,然后对其进行绘制:

# 计算投影矩阵
U, S, V = svd(trigger_effect, full_matrices=False)
trigger_projection_matrix = np.eye(3) - U @ U.T

# 将向量投影到正交平面上
projected_point = trigger_projection_matrix @ point
projected_vector = trigger_projection_matrix @ vector

# 绘制trigger_effect及其正交平面
ax = setup_3d_axes()
ax.plot_trisurf(x, y, z, color='C2', shade=False, alpha=0.25)
ax.quiver3D(*np.concatenate([origin, trigger_effect]).flatten(),
            arrow_length_ratio=0.1, color='C2', alpha=0.5)

# 绘制原始向量
ax.plot(*vector, color='k')
ax.plot(*point, color='k', marker='o')
offset = np.full((3, 1), 0.1)
ax.text(*(point + offset).flat, '({}, {}, {})'.format(*point.flat), color='k')

# 绘制投影向量
ax.plot(*projected_vector, color='C0')
ax.plot(*projected_point, color='C0', marker='o')
offset = np.full((3, 1), -0.2)
ax.text(*(projected_point + offset).flat,
        '({}, {}, {})'.format(*np.round(projected_point.flat, 2)),
        color='C0', horizontalalignment='right')

# 添加投影的虚线箭头
arrow_coords = np.concatenate([point, projected_point - point]).flatten()
ax.quiver3D(*arrow_coords, length=0.96, arrow_length_ratio=0.1,
            color='C1', linewidth=1, linestyle='dashed')

与以前一样,投影矩阵会将x,y,z空间中的任何点映射到该平面上,并且一旦某个点投影到该平面上,再次应用投影将不会产生任何效果。因此,很明显,尽管投影点在所有三个x、y和z方向上都有所不同,但是投影点集只有两个有效尺寸(即,它们被约束在一个平面上)。

EEG或MEG信号的投影几乎以相同的方式工作:点x,y,z对应于单个时间点上每个传感器的值,并且投影矩阵根据信号的哪些方面(比如,什么样的噪声)。唯一真正的区别是, 在实际案例中,要处理一系列N维“点”的时间序列(每次采样一个点),而不是单个三维点(x,y,z),其中N通常是几十个或几百个(取决于实验中EEG/MEG系统有多少个传感器)。

由于投影是矩阵运算,因此即使在具有数百个维度和数万个时间点的信号上也可以非常快速地完成投影。

本文分享自微信公众号 - 脑机接口社区(Brain_Computer),作者:Rose

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-12-17

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • EEG数据、伪影的查看与清洗

    在开始脑电(EEG)数据收集和分析之前,一定要确保你的数据尽可能的干净,这意味着收集的数据只是反映了大脑的活动。理论上听起来很简单,但实际上要注意“但是”。由于...

    脑机接口社区
  • eeglab中文教程系列(14)-绘制独立成分ERP贡献

    要完成该操作,必须保证已加载数据和电极位置数据,同时还要对数据进行提取epoch,并对数据进行ICA处理,操作如下:

    脑机接口社区
  • eeglab中文教程系列(4)-预处理工具

    本教程为脑机学习者Rose发表于公众号:脑机接口社区(微信号:Brain_Computer)

    脑机接口社区
  • 协议森林02 小喇叭开始广播 (以太网与WiFi协议)

    “小喇叭开始广播啦”,如果你知道这个,你一定是老一辈的人。“小喇叭”是五十年代到八十年代的儿童广播节目。在节目一开始,都会有一段这样的播音:“小朋友,小喇叭开始...

    Vamei
  • Python气象绘图教程(十五)—Cartopy_5

    本节提要:仿制中央气象台气象服务图片、关于cartopy里的投影与转换、cartopy中extent与boundary。

    气象学家
  • 使用HttpClient的优解

    新工作入职不久,目前仍然还在适应环境当中,笔者不得不说看别人的源码实在是令人痛苦。所幸前些日子终于将工作流畅地看了一遍,接下来就是熟悉框架技术的阶段了。 也正是...

    潘成涛
  • Openshift 4.4 静态 IP 离线安装系列(一):准备离线资源

    本系列文章描述了离线环境下以 UPI (User Provisioned Infrastructure) 模式安装 Openshift Container Pl...

    米开朗基杨
  • 如何在 Windows 10 中移除 Internet Explorer 浏览器 如何通过控制面板删除 Internet Explorer 浏览器通过 PowerShell 删

    现在 Internet Explorer (IE)已经过时了,可以通过控制面板移除这个古老但是依然是一个伟大的浏览器

    林德熙
  • 轻量级Qt键盘-中文输入

    中文输入候选栏ChineseWidget使用QListWidget和样式表实现:

    Qt君
  • 详解 MNIST 数据集

    MNIST 数据集已经是一个被"嚼烂"了的数据集, 很多教程都会对它"下手", 几乎成为一个 "典范". 不过有些人可能对它还不是很了解, 下面来介绍一下.

    用户1558438

扫码关注云+社区

领取腾讯云代金券