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

文章目录

前言
一、ItCCA的进步
二、在C-VEP的应用
- 1.C-VEP介绍
- 2.应用方法
三、在SSVEP的应用
- 1.标准CCA模板的局限
- 2. 实验结果：ITCCA的优势
四、matlab实现

前言

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

一、ItCCA的进步

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

$\hat{\pmb{u}}_k, \hat{\pmb{v}}_k = \underset{\pmb{u}_k, \pmb{v}_k} \argmax {\dfrac{{\rm Cov} \left(\pmb{u}_k \pmb{\mathcal{X}}, \pmb{v}_k \bar{\pmb{X}}_k \right)} {\sqrt{{\rm Var} \left(\pmb{u}_k \pmb{\mathcal{X}} \right)} \sqrt{{\rm Var}\left(\pmb{v}_k \bar{\pmb{X}}_k \right)}}} = \underset{\pmb{u}_k, \pmb{v}_k} \argmax {\dfrac{\pmb{u}_k \pmb{C}_{\pmb{\mathcal{X}} \bar{\pmb{X}}_k} {\pmb{v}_k}^T} {\sqrt{\pmb{u}_k \pmb{C}_{\pmb{\mathcal{X}} \pmb{\mathcal{X}}} {\pmb{u}_k}^T} \sqrt{\pmb{v}_k \pmb{C}_{\bar{\pmb{X}}_k \bar{\pmb{X}}_k} {\pmb{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), \quad n = 1, 2, \dots, N.\quad\quad(1)$

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

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

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

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

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

SSVEP_50">三、在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_64">四、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.9fprintf('当前注视的图像频率为 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.