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

Posted 2021-03-22 rosevorchid

目录

• projector(投影)和投影背景
• 案例解释投影原理
  • 导入工具库
  • 什么是projector(投影)？
  • 计算正交平面
  • 使用SVD计算投影矩阵

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

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系统有多少个传感器）。

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

参考：
信号空间投影SSP数学原理

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