脑机接口SSVEP经典算法 ITCCA个体模板典型相关分析 matlab实战

前言

        itCCA最开始用于C-VEP(code modulated VEP)信号的解码，这种信号的特征难以用正余弦波去描述，但确实存在稳定的类间差异性和类内统一性。因此基于正余弦信号的标准CCA算法，不适用于该情况。

一、ItCCA的进步

        需要把CCA的正余弦模板 Yk，替换成训练数据的叠加平均 X̅k

$${\hat{u}}_k,{\hat{v}}_k=\underset{u_k,v_k}{\mathrm{argmax}}\frac{\mathrm{Cov}\left(u_k\mathcal{X},v_k{\bar{X}}_k\right)}{\sqrt{\mathrm{Var}\left(u_k\mathcal{X}\right)}\sqrt{\mathrm{Var}\left(v_k{\bar{X}}_k\right)}}=\underset{u_k,v_k}{\mathrm{argmax}}\frac{u_k{C}_{\mathcal{X}{\bar{X}}_k}{v_k}^T}{\sqrt{u_k{C}_{\mathcal{X}\mathcal{X}}{u_k}^T}\sqrt{v_k{C}_{{\bar{X}}_k{\bar{X}}_k}{v_k}^T}}$$
总的来说，itCCA是对CCA的改进，但只是对参考信号进行了额外处理，对结果带来了改进。

二、在C-VEP的应用

1.C-VEP介绍

       C-VEP（Code Modulated Visual Evoked Potential，编码调制视觉诱发电位）信号是一种通过编码调制技术生成的视觉诱发电位（VEP）。VEP是神经生理学中用于评估视觉系统功能的信号。C-VEP在传统的VEP信号基础上，结合了编码技术，使得信号更加稳定，并且对噪声更具抗干扰性。它广泛应用于临床神经学和视觉科学中，特别是在视觉系统疾病的诊断中。

2.应用方法

        第一篇文章是应用在C-VEP的做法，是在(2)通过在K个周期上取平均，来获得参考模板。

        (1)在训练阶段，用户需要专注于参考目标。在我们的实验中，参考目标是目标T20。N个刺激周期内的EEG数据被收集为 $$x_n(t),n=1,2,\dots ,N.(1)$$

        (2)参考模板Mr(t)是通过在k个周期上求平均值而获得的。在我们的实验中，参考模板是T20的模板，即M20(t):
$$M_{20}(t)=\frac{1}{N}\sum _{n=1}^N{x}_n(t)(2)$$

        (3)通过移动参考模板获得所有目标的模板：
$$M_k(t)=M_{20}(t-({\tau }_k-{\tau }_{20}))k=0,1,2,\dots ,31(3)$$

        其中τk-τ20表示目标k和参考目标T20之间的时间滞后。

        (4)对于一段EEG数据x(t)，x(t)和模板Mk(t)之间的相关系数ρk计算为
$${\rho }_k=\frac{⟨M_k(t),x(t)⟩}{\sqrt{⟨M_k(t),M_k(t)⟩⟨x(t),x(t)⟩}}(4)$$

其中，<x, y> 表示x和y的乘积。
(5)通过选择使相关系数最大化的目标来识别注视目标。

三、在SSVEP的应用

       SSVEP（Steady-State Visual Evoked Potential，稳态视觉诱发电位）是一种由视觉刺激引起的大脑电活动反应，通常用于脑电图（EEG）研究和脑机接口（BCI）系统中。SSVEP 是大脑在看到周期性或频率固定的视觉刺激时所产生的电位反应，其特征是当一个人的眼睛注视着一个频率调制的视觉刺激源时，大脑的电活动会以与刺激频率相同的频率反应。

1.标准CCA模板的局限

        从受试者的大脑提取的基于SSVEP的EEG信号通常不是纯的，即包含某种形式的伪影。信号通常由三个分量组成：纯SSVEP、自发EEG信号和由其他非目标刺激的干扰引起的噪声。成功地将标准的CCA方法用于脑电信号的识别。

        如果用正余弦信号做模板，不包含自发EEG信号和噪声信号的特征，即不包含EEG信号所有分量的特征，只有纯SSVEP的特征。从这点上看，也就是说由于缺乏训练，CCA中使用的参考信号并不包含基于SSVEP的EEG信号的所有分量的特征。它只包含纯SSVEP的特征，而不包含自发EEG信号和噪声信号的特征。因此，标准CCA方法可能不会给出最佳结果。ITCCA方法是一种旨在解决标准CCA问题的方法。所以得到的结果不是最佳的。需要使用ITCCA，取EEG信号的特征，并将其添加到标准参考信号中，来优化人工参考信号。

2. 实验结果：ITCCA的优势

        实验结果如下图所示，在各方面，ITCCA都优于CCA，且总体准确率更高。原因是从受试者的大脑提取的基于SSVEP的EEG信号通常不是纯的，即包含某种形式的伪影。信号通常由三个分量组成：纯SSVEP、自发EEG信号和由其他非目标刺激的干扰引起的噪声。成功地将标准的CCA方法用于脑电信号的识别。

四、matlab实现

       由上可知，itCCA其实没有对CCA改变什么，只是对CCA输入的样本进行了处理。即以下代码的前8行。因为一个文件代表一个受试者。所以读取S002到S010共9个文件，9个人的实验数据，进行求和平均。将9个人的数据变成一个人的，这个就是 X̅k ，代码中是data_re，这是训练数据。 然后只读取第一个人的数据即x，在代码中是data_2，这是测试数据。要重复这个过程把所有频率的训练数据都得出来。最后所有得训练数据（已知频率），依次和这个训练数据做CCA，相关系数P大的就是当前训练数据对应频率。

数据集使用清华大学脑机接口团队的公开数据集，这里使用S1~S10
下载地址： https://bci.med.tsinghua.edu.cn/download.html

matlab代码

%********   itCCA.m  ******************
data_sum = zeros([8, 710]); % 初始化一个存储总和的矩阵
for i = 2:10    %取S002、S003....一直到S010，一共10个实验者数据
    filename = sprintf('S%03d.mat', i);
    mat_contents = load(filename); % 加载.mat文件
    data_1 = mat_contents.data(:, :, 1, 1, 1); % 提取同一块
    data_sum = data_sum + data_1; % 累加数据
end
data_re = data_sum / 9; % 求平均值

% 加载.mat文件
load('S001.mat');
% 获取data变量部分数据
data_2 = data(:,:,1,1,1);

x = data_re'; % 710x8
y = data_2';  % 710x6

[A, B, C1] = CCA(x, y);
P = C1;  % 相关系数p等于λ，这里特征值是λ的平方，后面开根号
fprintf('最大相关系数 p = %.5f\n', sqrt(P)); % %.5f 表示保留5位小数
A
B
if P > 0.9
    fprintf('当前注视的图像频率为 9.25 Hz\n');
end

%**********   CCA.m    ************
function [A,B,C1]=CCA(X,Y)
%CCA canonical correlation analysis
% [A,B]=cca(X,Y)
dimx=size(X,2); dimy=size(Y,2);  %列数(变量数)
Sall = cov([X Y]);  %cov求协方差矩阵。原阵X大小为M*N，则cov(X)大小为N*N的矩阵
dim=min(dimx,dimy);
Sxx=Sall(1:dim,1:dim);
Syy=Sall(dim+1:end,dim+1:end);
Sxy=Sall(1:dim,dim+1:end);
Syx=Sall(dim+1:end,1:dim);

iSxx=inv(Sxx);  %矩阵求逆
iSyy=inv(Syy);

M1 = iSxx*Sxy*iSyy*Syx;
M2 = iSyy*Syx*iSxx*Sxy;

[vec1,val1]=eig(M1);  %val对角阵，元素是特征值。 vec各列为对应特征值的特征向量
[vec2,val2]=eig(M2);

% A......找到每列的最大值和对应的行索引
[maxVals1, rowIndices1] = max(val1);
% 找到最大值的列标
[~, colIndex1] = max(maxVals1);
A =vec1(:, colIndex1);
C1 = max(maxVals1);
A = A * diag(diag((eye(size(A, 2)) ./ sqrt(A' * Sxx * A))))* -1;


% B......找到每列的最大值和对应的行索引
[maxVals2, rowIndices2] = max(val2);
% 找到最大值的列标
[~, colIndex2] = max(maxVals2);
B = vec2(:, colIndex2);
B = B * diag(diag((eye(size(B, 2)) ./ sqrt(B' * Syy * B))))* -1;

运行结果：

五、参考文献
Bin G, Gao X, Wang Y, Li Y, Hong B, Gao S. A high-speed BCI based on code modulation VEP. J Neural Eng. 2011 Apr;8(2):025015. doi: 10.1088/1741-2560/8/2/025015. Epub 2011 Mar 24. PMID: 21436527.
Nwachukwu, Sandra Ebele et al. “An SSVEP Recognition Method by Combining Individual Template with CCA.” Proceedings of the 2019 3rd International Conference on Innovation in Artificial Intelligence (2019): n. pag.