提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

        在脑机接口(BCI)领域，有SSVEP方向，中文叫做稳态视觉诱发电位，当人观看闪烁的视觉刺激（比如闪烁的灯光或图像）时，大脑会在与刺激频率一致的频率上，产生电位波动。这段脑电波是一种稳定且容易检测的脑电信号。通过确定ssvep波的频率，可以准确地判断用户正专注于哪个视觉刺激目标，帮助BCI系统识别用户意图，从而控制相应的设备或应用。例如，用户可以通过注视不同频率的方块或图标来控制计算机光标、游戏角色、轮椅等设备。还可以通过确定频率成分，对脑部健康状况进行初步筛查或监控，进而识别一些可能的脑部异常活动。

---

一、CCA的用途

        过去广泛使用的频率检测方法是PSDA，基于功率谱密度的分析，例如傅里叶变换，小波变换等。CCA（典型相关分析，Canonical Correlation Analysis） 是一种多变量统计分析方法，用于研究两组变量之间的相关性。在SSVEP-BCI（稳态视觉诱发电位-脑机接口）系统中，通常利用CCA来提取信号的频率特征并解码用户的意图。CCA的核心目标是通过求解典型相关系数来寻找两个信号集之间的关联性。

The most widely used frequency detection method in SSVEP-based BCIs is power spectral density based analysis。(PSDA)在基于SSVEP的脑机接口中，最广泛使用的频率检测方法是基于功率谱密度的分析。

---

二、频率求解思路

        使用CCA，需要输入两个变量，$X$ 和 $Y$。$X$ 是采集到的数据，$Y$ 是用已知的标准正余弦信号做为模板进行预设的数据。CCA旨在找到一对变换矩阵 $A$ 和 $B$ ，使得变换后的两个信号集之间的相关性最大化。相关性就体现在相关系数 $p$ 的求解上：
$$\text{max}\rho =\frac{\text{cov}(AX,BY)}{\sqrt{\text{cov}(AX,AX)\cdot \text{cov}(BY,BY)}}$$
        在得到的所有系数中，找到最大的系数。最大系数对应的频率，与SSVEP波频率相关性最大，就代表了SSVEP波的频率就是这个值。

---

三、输入数据结构

        使用数据结构说明如下描述，看图辅助文字理解。输入 $x$ 的行数为通道数8， $x$ 是采集到的数据。输入 $y$ 行数为使用的谐波数6， $y$ 是标准正余弦信号。 $x$ 和 $y$ 两者的列数为时间t。
        输入的 $x$ 和 $y$ 数据结构如下图

         $y$ 是标准正余弦信号，是下图这样的形式

CCA works on two sets of variables. In our method, variables in one set are the signals, x(t),recorded from several channels within a local region and the second set is from stimulus signals.
CCA处理两组变量。在我们的方法中，一组变量x(t)是从局部区域内的几个通道记录的信号，第二组变量来自刺激信号。

---

四、判断方法

        相关系数最大时的频率，就是我们要找的频率。而相关系数 $p={\lambda }^2$ ，这个关系由CCA理论推导得到。篇幅所限这里只记结论，后续会发另一篇文章，零基础理论学习CCA典型相关分析，内有详细推导，包教包会。这里想求出 $p$ ，只要找到对应特征值，开根号就好了。

Suppose there are K stimulus frequencies f1,f2,…fk and that the analyzed signal has been acquired from N channels within an Ls window. Our strategy for recognition is as follows. where p(f) is the CCA coefficient of xLN and y.
假设有K个刺激频率 f1，f2，…fk，并且所分析的信号是在一个长度为Ls的窗口内从N个通道中获取的。我们的识别策略如下。刺激频率fs满足 p(f) 是 xLN 和 y 的 CCA 系数。

In this approach, EEG signals from multiple channels are used to calculate the CCA coefficients with all stimulus frequencies in the system. The frequency with the largest coefficient is the one of SSVEP.
在该方法中，使用EEG信号来自多个通道，来计算系统中所有刺激频率的CCA系数。具有最大系数的频率是SSVEP的频率。

---

五、matlab实践

1.数据集获取及处理

        使用清华大学脑机接口研究组的公开数据集，链接 https://bci.med.tsinghua.edu.cn/
        此处下载使用Wearable SSVEP BCI Dataset的S1~S10。

        该数据集刺激信息如下图，一共12个频率。

        根据An Open Dataset for Wearable SSVEP-BasedBrain-ComputerInterfaces，这篇文章对该数据集的描述。

        该数据集包括102个MATLAB MAT文件，对应于102个受试者的脑电图数据。所有文件都是根据参与者的索引命名的（即S001.mat、S102.mat）。每个文件都包含一个名为“data”的5-D矩阵，维度为[8，710，2，10，12]，并存储为双精度浮点值。五个维度分别表示“通道标号”、“时间点”、“电极标号”、”块标号“和”目标标号“。每个矩阵对应240个历元（12个目标×10个块×2个电极），每个历元由未经任何处理的8个通道的原始脑电图数据组成，长度为2.84s（2.84s×250=710个时间点）
        根据连续EEG数据的事件通道中记录的刺激起始点，可以提取数据历元。每个数据历元的长度为2.84s，包括刺激开始前0.5s、视觉反应延迟0.14s、刺激后2s和刺激后0.2s。为了降低存储和计算成本，所有数据都被下采样到250 Hz。

        这里因为采样频率为250hz，所以一秒可以采集250个点，所以2.84s采集 $2.84s\times 250=710$ 点。 或者这样计算，采样点数 $N$，采样时间 $T$，采样间隔 $Ts$ ，采样频率 $fs$ ,
$$N=T/Ts$$
$$Ts=1/fs$$
所以
$$N=2.84/(1/250)=2.84x250=710点$$
        我们只使用采集的数据验证下CCA，不需要像文章里，对不同的通道/电极做对比评测。所以只需要拿第二个维度上，710个采集数据点使用即可。在这里用的 $x$ 是data_1只取8个通道刺激后2秒的数据，即2~2.84秒，只取最后0.84s的数据，对应210个点。所以 $y$ 在公式中的T要用到210个值。 $T/S=210/250=0.84$ 这也刚好说明是对应时间长度的。

2.matlab代码

        测试代码test_CCA.m如下，根据文中信息，数据 $x$ 要使用8个通道的数据。对 $y$ 的构建，只使用6个谐波。测试代码配合CCA算法，两个文件配合使用。这里仅读取S001.mat。
         test_CCA.m :

%*******  test_CCA.m  **********%
% 加载.mat文件
load('S001.mat');% 获取data变量部分数据
data_1 = data(:,501:710,1,1,1);  % 8x210% t  6x84
t = 1:210;  % 定义要用来循环的fre的值
fre_values = [9.25, 11.25, 13.25, 9.75, 11.75, 13.75, 10.25, 12.25, 14.25, 10.75, 12.75, 14.75]; % 可以根据需要添加更多值%% cca
x = data_1'; % 210x8
A_list = zeros([6, 12]); % 初始化存储 A 的多维数组
B_list = zeros([8, 12]); % 初始化存储 B 的多维数组
C1_list = zeros(1, length(fre_values)); % 初始化存储 C1 的单元数组% 循环通过fre_1的值并生成相应的data_2矩阵
for i = 1:length(fre_values)fre = fre_values(i);data_2 = [sin(2*pi*fre*t/250);cos(2*pi*fre*t/250);sin(4*pi*fre*t/250);cos(4*pi*fre*t/250);sin(6*pi*fre*t/250);cos(6*pi*fre*t/250)];y = data_2'; % 210x6[A, B, C1] = CCA(x, y);C1_list(i) = C1;A_list(:, i) = A;B_list(:, i) = B;
end
[max_value, index] = max(C1_list);
P = sqrt(max_value);  % 相关系数p等于λ，这里特征值是λ的平方，后面开根号
fprintf('最大相关系数 p = %.5f\n', P); % %.5f 表示保留5位小数
A = A_list(:, index)
B = B_list(:, index)
fprintf('当前注视的图像频率为 %.2f Hz\n', fre_values(index)); % %.2f 表示保留两位小数

---

        CCA算法, CCA.m :

%********  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;% 使用 SVD 分解
[U1, S1, V1] = svd(M1);  % M1 的奇异值分解
[U2, S2, V2] = svd(M2);  % M2 的奇异值分解% 手动单位化 U1 和 U2：使得每列（奇异向量）的模长为 1
for i = 1:size(U1, 2)U1(:, i) = U1(:, i) / sqrt(U1(:, i)' * U1(:, i));  % 每列单位化
endfor i = 1:size(U2, 2)U2(:, i) = U2(:, i) / sqrt(U2(:, i)' * U2(:, i));  % 每列单位化
end% 手动归一化 S1 和 S2，使得最大奇异值为 1
S1 = S1 / max(diag(S1));  % 将 S1 的最大奇异值归一化为 1
S2 = S2 / max(diag(S2));  % 将 S2 的最大奇异值归一化为 1% A ...... 选择最大奇异值对应的奇异向量
[maxVals1, colIndex1] = max(diag(S1));  % 找到最大奇异值所在的列
A = U1(:, colIndex1);  % 选择对应的奇异向量
C1 = maxVals1;  % 最大奇异值
A = A * diag(1 ./ sqrt(A' * Sxx * A));  % 归一化，注意矩阵的大小匹配% B ...... 选择最大奇异值对应的奇异向量
[maxVals2, colIndex2] = max(diag(S2));  % 找到最大奇异值所在的列
B = U2(:, colIndex2);  % 选择对应的奇异向量
B = B * diag(1 ./ sqrt(B' * Syy * B));  % 归一化，注意矩阵的大小匹配

---

3.运行及结果

        将测试代码test_CCA.m和CCA.m文件与数据S001.mat放到同一目录下，方便读取。

        运行结果为9.25hz，相关系数和 $A,B$ 具体值，matlab版本不同会有轻微差异。

六、参考文献

CCA学习：
Z. Lin, C. Zhang, W. Wu and X. Gao, “Frequency Recognition Based on Canonical Correlation Analysis for SSVEP-Based BCIs,” in IEEE Transactions on Biomedical Engineering, vol. 53, no. 12, pp. 2610-2614, Dec. 2006, doi: 10.1109/TBME.2006.886577.

数据集来源：
Zhu F, Jiang L, Dong G, Gao X, Wang Y. An Open Dataset for Wearable SSVEP-Based Brain-Computer Interfaces. Sensors (Basel). 2021 Feb 10;21(4):1256. doi: 10.3390/s21041256.