1.算法背景

脑机接口（Brain-Computer Interface, BCI）因其在神经工程与神经科学中的广泛应用价值而备受研究者们的关注。BCI系统可以在人类或动物被试与外部设备之间提供一个在线交流通道, 而无需依赖于常规的外围神经与肌肉输出通路。因其较高的时间分辨率、相对低的成本与高便携性, 从头皮采集的脑电图(electroencephalograhy, EEG)是被用于构建BCI系统的最广泛使用的数据模态。诸如皮层慢电位(slow cortical potentials)、感觉运动节律(sensorimotor rhythms)、P300诱发电位(P300 evoked potentials)、稳态视觉诱发电位(steady-state visual evoked potentials, SSVEPs)等脑电信号均可用于构建BCI系统。其中, SSVEPs由于其较高的信息传输速率(information transfer rate, ITR)、较高的信噪比(signal-to-noise, SNR)与仅依赖于较少的训练时间/无需训练等优点而备受研究者的青睐, 并已成为BCI系统重要的研究分支之一。

SSVEPs是由受试者注视某个快速重复的闪光或反转刺激后在受试者大脑枕叶区与顶叶区诱发的响应信号。SSVEPs实际上是近乎正弦的信号, 其信号成分中包含一系列离散的频率成分, 包括诱发刺激的基频与谐波成分。在基于SSVEPs的BCI系统中, 目标由单个刺激频率或不同刺激频率之间的组合来进行编码, 不同的目标编码则可以代表不同的控制命令。BCI系统则可以通过识别不同的目标编码以帮助BCI用户控制外围设备, 常见的外围设备有仿生机械腿、无人机、拨号面板、与虚拟键盘。

然而, 不同的用户的SSVEPs在振幅、相位以及数据分布上具有很大的差异。这个问题直接导致许多SSVEP频率识别算法依赖于很多参数调优步骤(如通道挑选或选择合适的数据长度), 以使得算法对于每个被试可发挥出尽可能高的效用。但是, 这些算法的优化需求将限制SSVEP-BCI系统的实际使用。因此, 一个无需依赖于调优步骤或仅依赖于少量调优过程的SSVEP频率识别算法亟待提出, 以使得BCI系统的使用变得便捷高效。

2.算法原理

从头皮采集的EEG信号可以被视作为大脑电生理活动反映在皮质层水平上的各电极的线性混合信号。多元同步索引(Multivariate Synchronization Index, MSI)算法旨在估计实际混合信号和参考信号之间的同步指数，以作为识别刺激频率的潜在指标。 具体而言, MSI算法使用S-estimator作为估计指数。

S-estimator基于多元信号相关矩阵的归一化特征值熵值。因为特征谱的方差与数据维度相关, 故数据维度是一个非常好的同步估计量。测量的同步量反比于数据嵌入维度(embedding dimensionality)。若两个时间序列数据完全无关, 数据嵌入维度将会最大化; 反之, 若两者之间高度相关(即完美同步), 则数据嵌入维度将会最小化。

该机理背后的逻辑在于: 数据嵌入维度越小, 多元信号相关矩阵的原始特征谱越稀疏, 则归一化后的特征值的熵值越大。反之, 数据嵌入维度越大, 多元信号相关矩阵的原始特征谱越稠密, 则归一化后的特征值熵值越小。

对于SSVEP频率识别任务而言, MSI算法中的多元信号对应于待测的SSVEP信号X∈RN×MX \\in \\mathbb{R}^{N \\times M}X∈RN×M与人工生成的正余弦参考信号Y∈R2Nh×MY \\in \\mathbb{R}^{2N_h \\times M}Y∈R2Nh​×M。其中, $N$ 表示通道数量, $M$ 代表时间点个数, 而 NhN_hNh​ 为谐波数量, 正余弦参考信号$Y$生成方式与CCA算法相同。MSI算法则通过计算 $X$ 与 $Y$ 之间的S-estimator指数以估计多元同步量, 并以此识别目标刺激的SSVEP频率。

为了计算$X$ 与 $Y$ 之间的S-estimator指数, MSI算法首先要求$X$ 与 $Y$ 归一化为均值0, 方差1的标准化矩阵。之后, 相关矩阵可被计算为:

C=[C11C12C21C22](1)C=\\begin{bmatrix} C11 & C12 \\\\ C21 & C22 \\\\ \\end{bmatrix} \\tag{1} C=[C11C21​C12C22​](1)
其中,

C11=1MXXT(2)C11=\\frac{1}{M}XX^T \\tag{2}C11=M1​XXT(2)

C22=1MYYT(3)C22=\\frac{1}{M}YY^T \\tag{3}C22=M1​YYT(3)

C12=1MXYT(4)C12=\\frac{1}{M}XY^T \\tag{4}C12=M1​XYT(4)

C21=(C12)T=1MYXT(4)C21=(C12)^T=\\frac{1}{M}YX^T \\tag{4}C21=(C12)T=M1​YXT(4)

在公式 (1) 中, 相关矩阵 $C$ 包含 $X$ 的自相关, $Y$ 的自相关以及 $X$ 与 $Y$ 的跨相关信息。由于自相关会影响同步的测量过程, 故如下的线性转换将在算法的实际实施过程中被采纳:

U=[C11(−1/2)00C22(−1/2)](5)U=\\begin{bmatrix} C11^{(-1/2)} & 0 \\\\ 0 & C22^{(-1/2)} \\\\ \\end{bmatrix} \\tag{5} U=[C11(−1/2)0​0C22(−1/2)​](5)

由此, 经线性转换后的相关矩阵 $R$ 将被表示为:

R=UCUT=[IN×NC11(−1/2)C12C22(−1/2)C22(−1/2)C21C11(−1/2)I2Nh×2Nh](6)R=UCU^T=\\begin{bmatrix} I_{N \\times N }& C11^{(-1/2)}C12C22^{(-1/2)} \\\\ C22^{(-1/2)}C21C11^{(-1/2)} & I_{2N_h \\times 2N_h}\\\\ \\end{bmatrix} \\tag{6} R=UCUT=[IN×N​C22(−1/2)C21C11(−1/2)​C11(−1/2)C12C22(−1/2)I2Nh​×2Nh​​​](6)

其中, IN×NI_{N \\times N}IN×N​ 与 I2Nh×2NhI_{2N_h \\times 2N_h}I2Nh​×2Nh​​分别表示维度为$N$ 与 2Nh2N_h2Nh​的单位矩阵。在经过公式 (5)中的线性转换后, 公式 (6) 中的自相关信息将被抵消。

若λ1,λ2,...,λP\\lambda_1, \\lambda_2, ..., \\lambda_Pλ1​,λ2​,...,λP​ 为相关矩阵 $R$ 的 $P$ 个特征值, 那么归一化后的特征值将被计算为如下结果:

λi′=λi∑i=1Pλi=λitr(R)(7)\\lambda_i^{'}=\\frac{\\lambda_i}{\\sum_{i=1}^{P}\\lambda_i}=\\frac{\\lambda_i}{tr(R)} \\tag{7}λi′​=∑i=1P​λi​λi​​=tr(R)λi​​(7)

其中, P=N+2NhP=N+2N_hP=N+2Nh​。最终, $X$ 与 $Y$ 之间的同步索引由以下公式计算得到:

S=1+∑i=1Pλi′log2(λi′)log2(P)(8)S=1+\\frac{\\sum_{i=1}^{P}\\lambda_i^{'}log_2(\\lambda_i^{'})}{log_2(P)} \\tag{8}S=1+log2​(P)∑i=1P​λi′​log2​(λi′​)​(8)

若 $X$ 与 $Y$ 完全无关, 则跨相关$C12=C21=0$, 相关矩阵 $R$ 将呈对角化, λi′=1/P,S=0\\lambda_i^{'}=1 / P, S=0λi′​=1/P,S=0。相反, 若$X$ 与 $Y$ 完美同步, 则仅有一个归一化的特征值为1, 其余特征值为全0, 此时$S=1$。对于其它情况, $S$ 将介于0与1之间。

在实际的SSVEP频率识别任务中, 若刺激目标数量为$K$, 即f1,f2,...,fKf_1, f_2, ..., f_Kf1​,f2​,...,fK​, 那么正余弦模板信号为Y={Y1,Y2,...,YK}∈RK×N×MY=\\{Y_1, Y_2, ..., Y_K\\}\\in\\mathbb{R}^{K \\times N \\times M}Y={Y1​,Y2​,...,YK​}∈RK×N×M。 由MSI算法识别的刺激目标频率为:

ftarget=argmax⁡fiSi(X,Yi),i=1,2,...,K(9)f_{target}=\\underset {f_i} {\\operatorname {argmax}} S_i(X, Y_i), \\ \\ i=1, 2, ..., K\\tag{9}ftarget​=fi​argmax​Si​(X,Yi​),  i=1,2,...,K(9)

3.Python代码实现

 import numpy as npclass MSI_Base():def __init__(self, opt):super(MSI_Base, self).__init__()# number of trialsself.Nh = opt.Nh# sample frequencyself.Fs = opt.Fs# number of stimulus targetsself.Nf = opt.Nf# window sizeself.ws = opt.ws# number of channelsself.Nc = opt.Nc# number of time pointsself.T = int(self.Fs * self.ws)def get_Reference_Signal(self, num_harmonics, targets):reference_signals = []t = np.arange(0, (self.T / self.Fs), step=1.0 / self.Fs)for f in targets:reference_f = []for h in range(1, num_harmonics + 1):reference_f.append(np.sin(2 * np.pi * h * f * t)[0:self.T])reference_f.append(np.cos(2 * np.pi * h * f * t)[0:self.T])reference_signals.append(reference_f)reference_signals = np.asarray(reference_signals)return reference_signalsdef find_Synchronization_Index(self, X, Y):num_freq = Y.shape[0]num_harm = Y.shape[1]result = np.zeros(num_freq)for freq_idx in range(0, num_freq):y = Y[freq_idx]'''normalized X and Y to have a zero mean and unitary variance'''X = X[:] - np.mean(X).repeat(self.T * self.Nc).reshape(self.Nc, self.T)X = X[:] / np.std(X).repeat(self.T * self.Nc).reshape(self.Nc, self.T)y = y[:] - np.mean(y).repeat(self.T * num_harm).reshape(num_harm, self.T)y = y[:] / np.std(y).repeat(self.T * num_harm).reshape(num_harm, self.T)'''equation (2)-(4), auto-correlation and cross-correlation'''c11 = (1 / self.T) * (np.dot(X, X.T))c22 = (1 / self.T) * (np.dot(y, y.T))c12 = (1 / self.T) * (np.dot(X, y.T))c21 = c12.T'''equation (1), correlation matrix C'''C_up = np.column_stack([c11, c12])C_down = np.column_stack([c21, c22])C = np.row_stack([C_up, C_down])# print("c11:", c11)# print("c22:", c22)v1, Q1 = np.linalg.eig(c11)v2, Q2 = np.linalg.eig(c22)V1 = np.diag(v1  (-0.5))V2 = np.diag(v2  (-0.5))'''equation (5), linear transformation'''C11 = np.dot(np.dot(Q1, V1.T), np.linalg.inv(Q1))C22 = np.dot(np.dot(Q2, V2.T), np.linalg.inv(Q2))# print("Q1 * Q1^(-1):", np.dot(Q1, np.linalg.inv(Q1)))# print("Q2 * Q2^(-1):", np.dot(Q2, np.linalg.inv(Q2)))U_up = np.column_stack([C11, np.zeros((self.Nc, num_harm))])U_down = np.column_stack([np.zeros((y.shape[0], self.Nc)), C22])U = np.row_stack([U_up, U_down])'''equation (6), correlation matrix R'''R = np.dot(np.dot(U, C), U.T)'''equation (7), normalized eigevalues λ'''eig_val, _ = np.linalg.eig(R)# print("eig_val:", eig_val, eig_val.shape)E = eig_val / np.sum(eig_val)'''equation (8), synchronization index S'''S = 1 + np.sum(E * np.log(E)) / np.log(self.Nc + num_harm)result[freq_idx] = Sreturn resultdef msi_classify(self, targets, test_data, num_harmonics=3):reference_signals = self.get_Reference_Signal(num_harmonics, targets)print("segmented_data.shape:", test_data.shape)print("reference_signals.shape:", reference_signals.shape)predicted_class = []labels = []num_segments = test_data.shape[0]num_perCls = num_segments // reference_signals.shape[0]for segment in range(0, num_segments):labels.append(segment // num_perCls)result = self.find_Synchronization_Index(test_data[segment], reference_signals)predicted_class.append(np.argmax(result) + 1)labels = np.array(labels) + 1predicted_class = np.array(predicted_class)return labels, predicted_class

完整实现代码: https://github.com/YuDongPan/Canonical_Classifier

参考资料与文献:

1. 应用于SSVEP脑电信号识别的CCA算法——CSDN博客
2. Zhang Y, Xu P, Cheng K, et al. Multivariate synchronization index for frequency recognition of SSVEP-based brain–computer interface[J]. Journal of neuroscience methods, 2014, 221: 32-40