脑机接口SSVEP经典算法 TRCA任务相关成分分析 matlab实战
文章目录
- 前言
- 一、TRCA是什么？
- 二、TRCA方法
- 三、matlab实战
- 1.数据集
- 2. test_trca.m
- 3. trca.m
- 4.运行结果
- 总结
- 参考文献
前言
传统的典型相关分析（CCA）方法在SSVEP信号识别中取得了一定成效，但容易受到自发脑电活动的干扰，且未充分利用相位信息。为此，研究者提出了TRCA方法，通过最大化任务期间神经影像数据的可重复性，提取任务相关成分，增强信号的信噪比（SNR），抑制背景脑电噪声。
一、TRCA是什么？
在SSVEP-BCI系统中，任务相关成分分析（Task-Related Component Analysis，TRCA）被用于提高SSVEP信号的检测性能。通过增强任务相关成分，TRCA能够提高信号的信噪比，抑制背景脑电活动的干扰，从而提升系统的准确性和稳定性。
二、TRCA方法
TRCA（任务相关成分分析）是一种通过最大化任务期间的可复现性来高效提取任务相关成分的方法。假设观测到的多通道 EEG 信号X(t)由以下线性生成模型表示：
x
j
(
t
)
=
a
1
,
j
s
(
t
)
+
a
2
,
j
n
(
t
)
,
j
=
1
,
2
,
…
,
N
c
x_j(t) = a_{1,j} s(t) + a_{2,j} n(t), \quad j = 1,2,\dots,N_c
xj(t)=a1,js(t)+a2,jn(t),j=1,2,…,Nc其中，
j
j
j 是通道索引，
a
1
,
j
a_{1,j}
a1,j 和
a
2
,
j
a_{2,j}
a2,j 是混合系数，用于将源信号投影到 EEG 观测信号中。
问题在于如何从观测信号 x(t) 的线性组合中恢复任务相关成分 s(t)，即：
y
(
t
)
=
∑
j
=
1
N
c
w
j
x
j
(
t
)
=
∑
j
=
1
N
c
(
w
j
a
1
,
j
s
(
t
)
+
w
j
a
2
,
j
n
(
t
)
)
y(t) = \sum_{j=1}^{N_c} w_j x_j(t) = \sum_{j=1}^{N_c} \left( w_j a_{1,j} s(t) + w_j a_{2,j} n(t) \right)
y(t)=j=1∑Ncwjxj(t)=j=1∑Nc(wja1,js(t)+wja2,jn(t))
理想情况下，该问题的解满足：
∑
j
=
1
N
c
w
j
a
1
,
j
=
1
\sum_{j=1}^{N_c} w_j a_{1,j} = 1
∑j=1Ncwja1,j=1 以及
∑
j
=
1
N
c
w
j
a
2
,
j
=
0
\sum_{j=1}^{N_c} w_j a_{2,j} = 0
∑j=1Ncwja2,j=0
从而得到最终解：
y
(
t
)
=
s
(
t
)
y(t) = s(t)
y(t)=s(t)
该问题可以通过跨试次协方差最大化（inter-trial covariance maximization）来求解。
C
h
1
h
2
=
Cov
(
y
(
h
1
)
(
t
)
,
y
(
h
2
)
(
t
)
)
C_{h_1 h_2} = \text{Cov} \left( y^{(h_1)}(t), y^{(h_2)}(t) \right)
Ch1h2=Cov(y(h1)(t),y(h2)(t))
第
h
ℎ
h 次试验的 EEG 信号和估计的任务相关成分分别表示为
x
(
h
)
(
t
)
和
y
(
h
)
(
t
)
,
其中
h
=
1
,
2.....
N
t
x^{(h)}(t)和y^{(h)}(t), 其中 h =1,2.....N_t
x(h)(t)和y(h)(t),其中h=1,2.....Nt
这里的
y
(
h
)
(
t
)
y^{(h)}(t)
y(h)(t)的时间区间固定为
t
∈
[
t
h
,
t
h
+
T
]
t∈[t_h,t_h+T]
t∈[th,th+T]，此处 𝑇 是每个任务试次的持续时间。
y
(
t
)
y(t)
y(t) 在试次
h
1
h_1
h1和
h
2
h_2
h2 之间的协方差表示：
C
h
1
h
2
=
Cov
(
y
(
h
1
)
(
t
)
,
y
(
h
2
)
(
t
)
)
=
∑
j
1
,
j
2
=
1
N
c
w
j
1
w
j
2
Cov
(
x
j
1
(
h
1
)
(
t
)
,
x
j
2
(
h
2
)
(
t
)
)
C_{h_1 h_2} = \text{Cov} \left( y^{(h_1)}(t), y^{(h_2)}(t) \right) = \sum_{j_1,j_2=1}^{N_c} w_{j_1} w_{j_2} \text{Cov} \left( x^{(h_1)}_{j_1}(t), x^{(h_2)}_{j_2}(t) \right)
Ch1h2=Cov(y(h1)(t),y(h2)(t))=j1,j2=1∑Ncwj1wj2Cov(xj1(h1)(t),xj2(h2)(t))
所有可能的试次组合被总结如下：
∑
h
1
,
h
2
=
1
h
1
≠
h
2
N
t
C
h
1
h
2
=
∑
h
1
,
h
2
=
1
h
1
≠
h
2
N
t
∑
j
1
,
j
2
=
1
N
c
w
j
1
w
j
2
Cov
(
x
j
1
(
h
1
)
(
t
)
,
x
j
2
(
h
2
)
(
t
)
)
=
w
T
S
w
\sum_{\substack{h_1, h_2 = 1 \\ h_1 \neq h_2}}^{N_t} C_{h_1 h_2} = \sum_{\substack{h_1, h_2 = 1 \\ h_1 \neq h_2}}^{N_t} \sum_{j_1,j_2=1}^{N_c} w_{j_1} w_{j_2} \text{Cov} \left( x^{(h_1)}_{j_1}(t), x^{(h_2)}_{j_2}(t) \right) = w^T S w
h1,h2=1h1=h2∑NtCh1h2=h1,h2=1h1=h2∑Ntj1,j2=1∑Ncwj1wj2Cov(xj1(h1)(t),xj2(h2)(t))=wTSw
这里，矩阵 S = ( S j 1 j 2 ) 1 ≤ j 1 , j 2 ≤ N c S = (S_{j_1 j_2})_{1 \leq j_1, j_2 \leq N_c} S=(Sj1j2)1≤j1,j2≤Nc会被定义为 S j 1 j 2 = ∑ h 1 , h 2 = 1 h 1 ≠ h 2 N t Cov ( x j 1 ( h 1 ) ( t ) , x j 2 ( h 2 ) ( t ) ) S_{j_1 j_2} = \sum_{\substack{h_1, h_2 = 1 \\ h_1 \neq h_2}}^{N_t} \text{Cov} \left( x^{(h_1)}_{j_1}(t), x^{(h_2)}_{j_2}(t) \right) Sj1j2=∑h1,h2=1h1=h2NtCov(xj1(h1)(t),xj2(h2)(t))
为了获得有限解， y ( t ) y(t) y(t)的方差被约束为：
Var ( y ( t ) ) = ∑ j 1 , j 2 = 1 N c w j 1 w j 2 Cov ( x j 1 ( t ) , x j 2 ( t ) ) = w T Q w = 1 \text{Var}(y(t)) = \sum_{j_1,j_2=1}^{N_c} w_{j_1} w_{j_2} \text{Cov} (x_{j_1}(t), x_{j_2}(t)) = w^T Q w = 1 Var(y(t))=j1,j2=1∑Ncwj1wj2Cov(xj1(t),xj2(t))=wTQw=1
该约束优化问题可以求解为：
w ^ = arg max w w T S w w T Q w \hat{w} = \arg \max_w \frac{w^T S w}{w^T Q w} w^=argwmaxwTQwwTSw
最优的系数向量可以通过矩阵 Q − 1 S Q^{-1} S Q−1S的特征向量获得。
一般来说，从一个 N c × N c N_c \times N_c Nc×Nc矩阵中可以得到 N c N_c Nc个特征值和特征向量。矩阵 Q − 1 S Q^{-1} S Q−1S的特征值 λ \lambda λ 可以按降序排列，它表示相应特征向量 w ^ \hat{w} w^ 的成本函数值。因此，这些特征值反映了多个试次之间的任务一致性。
图 (a) 显示了基于 TRCA 方法的流程图。在基于 SSVEP 的脑机接口（BCI）中，TRCA 可用于设计空间滤波器，以去除头皮记录中的背景 EEG 活动。
针对第
n
n
n 个刺激的空间滤波器，
w
n
(
m
)
∈
R
N
c
w_n^{(m)} \in \mathbb{R}^{N_c}
wn(m)∈RNc 通过 TRCA 从个体校准数据
χ
n
(
m
)
\chi_n^{(m)}
χn(m)中获取，并在应用滤波器组分析后计算得到。在这里，公式
Var
(
y
(
t
)
)
=
w
T
Q
w
=
1
\text{Var}(y(t)) = w^T Q w = 1
Var(y(t))=wTQw=1中的
Q
Q
Q 矩阵是使用所有训练试次的拼接矩阵
χ
n
(
m
)
\chi_n^{(m)}
χn(m)计算得到的。
单试次测试数据
X
(
m
)
∈
R
N
c
×
N
s
X^{(m)} \in \mathbb{R}^{N_c \times N_s}
X(m)∈RNc×Ns与跨试次平均的第
n
n
n个视觉刺激训练数据
χ
ˉ
n
(
m
)
∈
R
N
c
×
N
s
\bar{\chi}_n^{(m)} \in \mathbb{R}^{N_c \times N_s}
χˉn(m)∈RNc×Ns之间的相关系数计算如下：
T
r
n
(
m
)
=
ρ
(
(
X
(
m
)
)
T
w
n
(
m
)
,
χ
ˉ
n
(
m
)
w
n
(
m
)
)
T r_n^{(m)} = \rho \left( \left( X^{(m)} \right)^T w_n^{(m)}, \bar{\chi}_n^{(m)} w_n^{(m)} \right)
Trn(m)=ρ((X(m))Twn(m),χˉn(m)wn(m))
其中，
p
(
a
,
b
)
p(a,b)
p(a,b) 表示信号
a
a
a 和
b
b
b 之间的皮尔逊相关分析。
三、matlab实战
1.数据集
使用清华大学脑机接口团队的公开数据集，benchmark dataset
下载地址： https://bci.med.tsinghua.edu.cn/
关于该数据集的介绍和使用，可以看这篇文章：
清华大学脑机接口研究组Benchmark dataset数据集介绍
2. test_trca.m
trca的测试代码如下，主要是读取数据文件，然后对数据滤波处理后，送入trca算法。
% 加载所需的数据文件
load ('Freq_Phase.mat'); % 载入频率和相位数据
data = load('S1.mat'); % 载入EEG数据
eeg_data = data.data; % 获取EEG数据矩阵
% 处理数据
% 原始数据，维度是 [64, 1500, 40, 6]（64 个电极，1500 个时间点, 40目标刺激, 6个试次）
% 提取顶叶和枕叶区域（Pz、PO5、PO3、POz、PO4、PO6、O1、Oz和O2）的SSVEP信号
stim_data = eeg_data([48, 54, 55, 56, 57, 58, 61, 62, 63], 251:1250, :, :);
% 对stim_data进行5-90Hz带通滤波
% stim_data的维度是 [9, 1500, 40, 6]
% 设置滤波器参数
low_freq = 5; % 下边界频率
high_freq = 90; % 上边界频率
fs = 250; % 采样频率为250Hz
% 获取数据的维度
[N_channels, N_timepoints, N_targets, N_blocks] = size(stim_data);
% 对每个数据块和每个目标进行滤波
for block_idx = 1:N_blocks
for target_idx = 1:N_targets
for channel_idx = 1:N_channels
% 提取每个通道、目标和数据块的信号
signal = squeeze(stim_data(channel_idx, :, target_idx, block_idx));
% 使用bandpass滤波器对信号进行5-90Hz滤波
stim_data(channel_idx, :, target_idx, block_idx) = bandpass(signal, [low_freq, high_freq], fs);
end
end
end
eeg = stim_data;
% 获取EEG数据的维度：通道数、时间点数、目标刺激数、数据块数
[N_channel,~, N_target, N_block] = size(eeg);
%% ------------分类部分-------------
tic
% 留一交叉验证（LOO Cross-Validation）
for loocv_i = 1:N_block
% 提取测试数据（留一交叉验证中，当前块作为测试数据）
Testdata = squeeze(eeg(:, :, :, loocv_i));
% 提取训练数据（去除当前块数据作为训练集）
Traindata = eeg;
Traindata(:, :, :, loocv_i) = [];
% 对每个目标频率，计算训练数据的平均值
for targ_i = 1:N_target
aver_Traindata(:, :, targ_i) = squeeze(mean(squeeze(Traindata(:,:,targ_i,:)),3)); % 对每个目标频率进行平均
end % end targ_i
% 标签赋值，真实标签是频率
truelabels = freqs;
% 获取测试数据的试次数量
N_testTrial = size(Testdata, 3);
% 对每个测试试次进行处理
for trial_i = 1:N_testTrial
coefficience = zeros(1, length(truelabels)); % 初始化系数矩阵
% 对每个频率进行处理
for targ_j = 1:length(freqs)
% 使用训练数据计算空间滤波器 wn
wn = TRCA(squeeze(Traindata(:,:,targ_j,:))); % 调用TRCA函数计算空间滤波器
% 计算加权训练数据与测试数据的相关性
weighted_train = wn' * aver_Traindata(:,:,targ_j); % 加权训练数据
weighted_test = wn' * Testdata(:,:,trial_i); % 加权测试数据
% 计算测试数据和训练数据之间的相关系数
coefficienceMatrix = corrcoef(weighted_test, weighted_train);
coefficience(targ_j) = coefficienceMatrix(1, 2); % 获取相关系数
end % end targ_i
% 目标检测：选择具有最大相关系数的频率作为输出标签
[~, index] = max(coefficience); % 获取最大相关系数对应的频率
outputlabels(trial_i) = freqs(index); % 输出预测的标签
end % end trial_i
% 计算该次测试的正确分类数量
trueNum = sum((outputlabels - truelabels) == 0); % 统计预测标签与真实标签相同的个数
acc(loocv_i) = trueNum / length(truelabels); % 计算该次交叉验证的准确率
% 打印当前交叉验证的准确率
fprintf('第%d试次的交叉验证准确率是: %.4f, 样本数: %d/%d\n', loocv_i, acc(loocv_i), trueNum, N_testTrial)
end % end looCv_i
% 计算总的运行时间
t = toc;
% 打印总的运行时间和平均准确率
fprintf('\n-----------------------------------------\n')
disp(['总运行时间: ', num2str(t), ' s']);
fprintf('6试次数据交叉验证的平均准确率是: %.4f\n', mean(acc));
3. trca.m
trca的计算代码
% -------------------------------------------------------------------
% Task-related component analysis[1].extract task-related components
% by maximizing the reproducibility during the task period.
%
% Input:
% eeg : input eeg training data
% (# channels, # points, # trials)
% Output:
% W: w means eigenvectors corresponding to max eigenvalues
%
% Reference:
% [1] H. Tanaka, T. Katura, H. Sato,
% "Task-related component analysis for functional neuroimaging and
% application to near-infrared spectroscopy data",
% NeuroImage, vol. 64, pp. 308-327, 2013.
% --------------------------------------------------------------------
function [w] = TRCA(eeg)
% 输入：eeg：一个三维矩阵，包含 EEG 数据（通道数 × 样本数 × 试次数）
% 输出：w：任务相关成分的空间滤波器（特征向量）
% 获取 EEG 数据的维度
[num_chans, num_smpls, num_trials] = size(eeg);
% 初始化协方差矩阵 S 为零矩阵
S = zeros(num_chans);
% 计算协方差矩阵 S，双重循环遍历每一对试次
for trial_i = 1:1:num_trials-1 % 遍历第一个试次
x1 = squeeze(eeg(:,:,trial_i)); % 提取第 trial_i 试次的数据（去掉第三维）
x1 = bsxfun(@minus, x1, mean(x1, 2)); % 对每个通道去均值，避免偏差
for trial_j = trial_i+1:1:num_trials % 遍历第二个试次（只考虑 trial_i 和 trial_j 的组合）
x2 = squeeze(eeg(:,:,trial_j)); % 提取第 trial_j 试次的数据
x2 = bsxfun(@minus, x2, mean(x2, 2)); % 对每个通道去均值
% 计算协方差矩阵的增量，并更新 S
S = S + x1*x2' + x2*x1';
end % end trial_j
end % end trial_i
% 将 EEG 数据重塑为一个二维矩阵，每一列表示一个试次的所有样本
UX = reshape(eeg, num_chans, num_smpls*num_trials);
% 对重塑后的数据进行去均值处理
UX = bsxfun(@minus, UX, mean(UX, 2));
% 计算 UX 的自相关矩阵 Q
Q = UX*UX';
% 通过广义特征值问题求解空间滤波器 w
[V, ~] = eigs(S, Q); % 求解 S 和 Q 的广义特征值问题，V 为特征向量
% 选择最大特征值对应的特征向量作为空间滤波器
w = V(:, 1); % 返回第一个特征向量作为空间滤波器（任务相关成分）
end %end function
4.运行结果
总结
TRCA作为一种有效的任务相关信号提取方法，在SSVEP信号检测中发挥了重要作用。通过最大化任务相关成分的可重复性，TRCA提高了信号的信噪比，增强了SSVEP信号的检测性能。与其他方法相比，TRCA在处理背景脑电噪声方面表现出色，为SSVEP-BCI系统的优化提供了新的思路。
参考文献
M. Nakanishi, Y. Wang, X. Chen, Y. -T. Wang, X. Gao and T. -P. Jung, “Enhancing Detection of SSVEPs for a High-Speed Brain Speller Using Task-Related Component Analysis,” in IEEE Transactions on Biomedical Engineering, vol. 65, no. 1, pp. 104-112, Jan. 2018, doi: 10.1109/TBME.2017.2694818.