稳态视觉诱发电位(SSVEP)识别| Task-Related Component Analysis, TRCA

Task-Related Component Analysis(TRCA)是目前SSVEP识别最流行的方法之一，在此方法提出之前多用CCA或是MSI进行频率识别，这些方法最大的好处就是不需要受试者的训练数据。随着联合频相刺激的SSVEP范式提出，如何有效利用SSVEP的相位信息变得很重要，但标准的CCA和MSI是没有利用到相位信息的（因为人工构造的正余弦参考信号没有包含相位项）。如果能够有效的利用相位信息，那么SSVEP的识别必然会提升较大的幅度，这也是为什么subject independent template的引入能够显著提高识别能力的原因。但是这又导致了另外的一个问题，那就是需要训练数据集了！毕竟要在ACC和ITR上妥协嘛~
Masaki Nakanishi 等人率先将TRCA方法引入到SSVEP识别上来1，而该方法最早是于2013年Hirokazu Tanaka 等人提出并应用于NIRS的识别上2 3。 该方法的能够通过最大化每个task中神经影像数据的复现性(reproducibility)，从而提取任务相关成分。对于SSVEP这种time-locked的信号，该方法非常合适。在SSVEP中即最大化多个trial之间的复现性，从而提高信噪比，抑制自发脑电活动。下面的介绍我均以SSVEP中的说法为例，给出的示意图来自这三篇参考文献。

首先，我们认为一个多通道的EEG信号​是由任务相关(task-related)信号和任务无关信号组合而成的，即

其中，分别表示混合系数。以单通道的EEG信号为例，在连续的三个trial（蓝色部分）中可以图示为

可以看到，任务相关部分在trial之间保持着不变性，而任务无关部分在trial之间是变化的。它们两者存在着这样一个性质，即之间的协方差为一正常数，而与之间的协方差为0。
我们希望从EEG信号中仅得到任务相关的成分s(t)，这样即能增加信号的信噪比。首先我们对多通道的EEG信号进行加权求和，表示为一个线性一维模型

为了达到我们所期待的目标，使得估计的任务相关成分，那么就必须使0。如果我们得到了这样的一个解，那么多通道的EEG信号经过线性加权求和后，便能够得到在多个trial之间具有高度相关性的y(t)，表示为下图

蓝色部分表示不同的trial，可以看到原本trial之间变化较大的N通道EEG信号，经过处理后得到的y ( t ) y(t)y(t)在trial之间相关性较高。这里的w i w_iwi​即空间滤波器，这个过程可以表示为，接下来的问题就是如何对W WW求解！可以通过两个类似的求解方法实现：协方差最大化(CovMax)或相关性最大化(CorrMax)。

定义来自k-th trial的EEG信号和其估计的任务相关成分为，同理来自l-th trial的信号分别为和，它们之间的相关性可以通过皮尔逊相关系数公式求解

在实际应用过程中，我们需要找到的是所有trial可能的组合中使相关系数之和最大的解，所以我们的目标函数为

上式即为Correlation maximization(CorrMax)。为了对系数进行约束，我们使y(t)的方差归一化到1（注：这里的N表示channel）

其中，。
采用CorrMax计算并没有什么问题，但是它的解并不是封闭解析的，最重要的是它只能得到一个任务相关成分！所以在SSVEP识别的过程中，并没有采用该方法，而是使用下面所述的最大化trial之间的协方差。

这里和CorrMax不同的仅是将trial之间相关系数的求解变为计算协方差矩阵。类似地，计算两个trial之间的协方差

多个trial之间的协方差之和为我们的目标函数，该公式为Covariance maximization(CovMax)

上式中对称矩阵S SS可以表示为

和CorrMax中一样，我们需要引入对y(t)的归一化约束，我们的约束优化问题转变为Rayleigh–Ritz特征值问题

根据Rayleigh–Ritz理论，最佳的系数向量可以通过求解的特征向量得到。的特征向量为一N维的向量矩阵，每一维对应一特征值，并且按特征值大小降序排列。特征值的大小表示了按其对应的特征向量对信号进行空间滤波后，y(t)在trial之间的task consistency。在SSVEP识别过程中，仅使用了最大特征值对应的特征向量。

以上，我们已经知道了TRCA的原理，也知道了具体的求解方法——CovMax。
下面进行实际代码操作的时候还是先梳理一下整个SSVEP识别过程：
首先我们需要明白，CovMax的计算要用到已有的trial数据，表明了TRCA需要受试者训练数据。我们利用受试者某一刺激n的训练数据，通过TRCA得到对应的空间滤波器w n w_nwn​；利用w n w_nwn​分别对test trial X和训练集平均信号χ ˉ n \bar\chi_nχˉ​n​进行空间滤波；然后计算两者的皮尔逊相关系数；最后找到最大相关系数对应的刺激频率即可。多说一句，目前类别数较多时，SSVEP的识别主要是根据相关系数来的，各方法的改进之处在于如何找到一个更好的空间滤波器w n w_nwn​。下面是TRCA和ensemble TRCA的流程图。

<span style="color:#000000"><code>function [w] = TRCA(eeg)
% -------------------------------------------------------------------
% 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.
% --------------------------------------------------------------------
[num_chans, num_smpls, num_trials]  = size(eeg);
S = zeros(num_chans);
for trial_i = 1:1:num_trials-1x1 = squeeze(eeg(:,:,trial_i));x1 = bsxfun(@minus, x1, mean(x1,2));for trial_j = trial_i+1:1:num_trialsx2 = squeeze(eeg(:,:,trial_j));x2 = bsxfun(@minus, x2, mean(x2,2));S = S + x1*x2' + x2*x1';end % end trial_j
end % end trial_i
UX = reshape(eeg, num_chans, num_smpls*num_trials);
UX = bsxfun(@minus, UX, mean(UX,2));
Q = UX*UX';
[V,~] = eigs(S, Q);
w = V(:,1);
end %end function
</code></span>

对上述程序做简要解释

• 输入参数eeg为三维的信号矩阵，表示受试者某一刺激频率的训练数据；
• 输出参数为得到的空间滤波器w；
• CovMax的关键在于Q − 1 S Q^{-1}SQ−1S的求解；
• Q矩阵来自于对所有trial进行横向联结(# channel, # points*trials)；
• w仅取了最大特征值对应的特征向量。

<span style="color:#000000"><code>% -------------------------------------------------------------------------
% Main for task-related component analysis
%
% Dataset (Sx.mat):
%   A 40-target SSVEP dataset recorded from a single subject. The stimuli
%   were generated by the joint frequency-phase modulation (JFPM) 
%     - Stimulus frequencies    : 8.0 - 15.8 Hz with an interval of 0.2 Hz
%     - Stimulus phases         : 0pi, 0.5pi, 1.0pi, and 1.5pi
%     - # of channels           :Oz
%     - # of recording blocks   : 6
%     - Data length of epochs   : 1.5 [seconds]
%     - Sampling rate           : 250 [Hz]
%     - Data format             : # channels, # points, # targets, # blocks
% 
% See also:
%   TRCA.m
% -------------------------------------------------------------------------clear all
close allload ('Freq_Phase.mat')
load('subject1.mat')
eeg = subject1;
[N_channel,~, N_target, N_block] = size(eeg);%% ------------classification-------------
tic
% LOO cross-validation
for loocv_i = 1:N_blockTestdata = squeeze(eeg(:, :, :, loocv_i));Traindata = eeg;Traindata(:, :, :, loocv_i) = [];for targ_i = 1:N_targetaver_Traindata(:, :, targ_i) = squeeze(mean(squeeze(Traindata(:,:,targ_i,:)),3));end % end targ_i% labels assignment according to testdatatruelabels=freqs;N_testTrial=size(Testdata, 3);for trial_i=1:N_testTrialcoefficience = zeros(1,length(truelabels));for targ_j=1:length(freqs)             % compute spatial filter wn using training datawn = TRCA(squeeze(Traindata(:,:,targ_j,:)));% compute correlation between test and averaged training dataweighted_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% target detection[~, index] = max(coefficience);outputlabels(trial_i) = freqs(index);end % end trial_itrueNum = sum((outputlabels-truelabels)==0);acc(loocv_i) = trueNum/length(truelabels);fprintf('The %d-th CV accuracy is: %.4f, samples: %d/%d\n',loocv_i,...acc(loocv_i),trueNum, N_testTrial)
end % end looCv_i
t=toc;
% data visualization
fprintf('\n-----------------------------------------\n')
disp(['total time: ',num2str(t),' s']);
fprintf('6-fold CV average accuracy is: %.4f\n',mean(acc))
</code></span>

对上述程序做简要解释

• 采用的数据为清华benchmark dataset数据，导入的是对受试者1进行1.5s数据截取和滤波后的数据；
• 该数据包含6 blocks，故进行了6次LOO；
• 该数据为40类，8~15.8Hz；
• 完整可执行程序见我的TRCA仓库。

运行结果为

<span style="color:#000000"><code>The 1-th CV accuracy is: 0.9250, samples: 37/40
The 2-th CV accuracy is: 0.9750, samples: 39/40
The 3-th CV accuracy is: 1.0000, samples: 40/40
The 4-th CV accuracy is: 0.9250, samples: 37/40
The 5-th CV accuracy is: 0.9750, samples: 39/40
The 6-th CV accuracy is: 0.9500, samples: 38/40-----------------------------------------
total time: 53.7772 s
6-fold CV average accuracy is: 0.9583
</code></span>

以上，是我关于TRCA的总结。看到这里如果有疑问欢迎下方留言，如果有什么见解麻烦不吝赐教，一起交流学习~

1. Nakanishi M, Wang YJ, Chen XG, Wang YT, Gao XR, Jung TP (2018) Enhancing detection of ssveps for a high-speed brain speller using task-related component analysis. Ieee Transactions on Biomedical Engineering 65:104-112.
2. Tanaka H, Katura T, Sato H (2013) Task-related component analysis for functional neuroimaging and application to near-infrared spectroscopy data. NeuroImage 64:308-327.
3. Tanaka H, Katura T, Sato H (2014) Task-related oxygenation and cerebral blood volume changes estimated from NIRS signals in motor and cognitive tasks. NeuroImage 94:107-119.