声明：转载自https://blog.csdn.net/weixin_42765703/article/details/105604481，感谢博主详尽的分析，最近正在看，学习一下
Note： 没事搬一下砖，重温一下TRCA的计算过程。在知网上检索过SSVEP的硕博论文，发现到2020了更多的仍然是关于CCA及其改进方法的介绍（可能是因为标准CCA毕竟不需要训练，而且计算速度较快），所以我还是只有看期刊上的推导公式。这里我简单谈谈TRCA在SSVEP处理上的实现过程，有不对的地方恳请各位指正。
已有基础仅仅是来找代码的可以直接转Masaki Nakanishi的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)。
Correlation maximization(CorrMax)
定义来自k-th trial的EEG信号和其估计的任务相关成分为
，同理来自l-th trial的信号分别为
和
，它们之间的相关性可以通过皮尔逊相关系数公式求解
在实际应用过程中，我们需要找到的是所有trial可能的组合中使相关系数之和最大的解，所以我们的目标函数为
上式即为Correlation maximization(CorrMax)。为了对系数进行约束，我们使y(t)的方差归一化到1（注：这里的N表示channel）
其中，
。
采用CorrMax计算并没有什么问题，但是它的解并不是封闭解析的，最重要的是它只能得到一个任务相关成分！所以在SSVEP识别的过程中，并没有采用该方法，而是使用下面所述的最大化trial之间的协方差。
Covariance maximization(CovMax)
这里和CorrMax不同的仅是将trial之间相关系数的求解变为计算协方差矩阵。类似地，计算两个trial之间的协方差
多个trial之间的协方差之和为我们的目标函数，该公式为Covariance maximization(CovMax)
上式中对称矩阵S SS可以表示为
和CorrMax中一样，我们需要引入对y(t)的归一化约束，我们的约束优化问题转变为Rayleigh–Ritz特征值问题
根据Rayleigh–Ritz理论，最佳的系数向量
可以通过求解
的特征向量得到。
的特征向量为一N维的向量矩阵，每一维对应一特征值，并且按特征值大小降序排列。特征值的大小表示了按其对应的特征向量对信号进行空间滤波后，y(t)在trial之间的task consistency。在SSVEP识别过程中，仅使用了最大特征值对应的特征向量。
MATLAB代码
以上，我们已经知道了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的流程图。
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-1
x1 = squeeze(eeg(:,:,trial_i));
x1 = bsxfun(@minus, x1, mean(x1,2));
for trial_j = trial_i+1:1:num_trials
x2 = 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仅取了最大特征值对应的特征向量。
Classification
<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 all
load ('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_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
% labels assignment according to testdata
truelabels=freqs;
N_testTrial=size(Testdata, 3);
for trial_i=1:N_testTrial
coefficience = zeros(1,length(truelabels));
for targ_j=1:length(freqs)
% compute spatial filter wn using training data
wn = TRCA(squeeze(Traindata(:,:,targ_j,:)));
% compute correlation between test and averaged training data
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
% target detection
[~, index] = max(coefficience);
outputlabels(trial_i) = freqs(index);
end % end trial_i
trueNum = 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的总结。看到这里如果有疑问欢迎下方留言，如果有什么见解麻烦不吝赐教，一起交流学习~
References
- 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.
- 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.
- 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.
1359
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
