共空间模式 Common Spatial Pattern(CSP)原理和实战

共空间模式CSP

共空间模式（Common Spatial Pattern, CSP）是一种对两分类任务下的空域滤波特征提取算法，能够从多通道的脑机接口数据里面提取出每一类的空间分布成分。公共空间模式算法的基本原理是利用矩阵的对角化，找到一组最优空间滤波器进行投影，使得两类信号的方差值差异最大化，从而得到具有较高区分度的特征向量。

共空间模式理论

假设$X_1$和$X_2$分别为两分类想象运动任务下的多通道诱发响应时-空信号矩阵，他们的维数均为$N\ast T$,$N$为脑电通道数，$T$为每个通道所采集的样本数。为了计算其协方差矩阵，现在假设$N<T$.在两种脑电想象任务情况下，一般采用复合源的数学模型描述EEG信号。为了简化计算，常忽略噪声所带来的影响。$X_1$和$X_2$分别为：
$$X_i=\left[C_i{C}_M\right]\left[\begin{array}{c}{S}_i\\ S_M\end{array}\right],(i=1,2)\text{\hspace{1em}(1)}$$
上式(1)中，$S_1$和$S_2$分别代表两种类型任务。假设两种信号源是相互线性独立的；$S_M$代表两种类型任务下所共同拥有的源信号，假设$S_1$是由$m_1$个源所构成的，$S_1$是由$m_2$个源所构成，则$C_1$和$C_2$便是由$S_1$和$S_2$相关的$m_1$和$m_2$个共同空间模式组成的。由于每个空间模式都是一个$N\ast 1$维的向量，现在用该向量来表示单个的源信号所引起的信号在$N$个导联上的分布权重。$C_M$表示的是与$S_M$相应的共有的空间模式。CSP算法的目标就是要设计空间滤波器$F_1$和$F_2$得到空间因子$W$。

1.求解协方差矩阵

时空信号矩阵$X_1$和$X_2$归一化的协方差矩阵$R_1$和$R_2$:
$$R_i=\frac{X_i{X}_i^T}{trace\left(X_i{X}_i^T\right)}(i=1,2)\text{\hspace{1em}(2)}$$
上式(2)中，$X^T$表示矩阵$X$的转置，$trace(X)$表示对矩阵对角线上元素求和。之后求解混合空间的协方差矩阵$R$:
$$R=\overline{R_1}+\overline{R_2}\text{\hspace{1em}(3)}$$
上式中，$\overline{R_i}(i=1,2)$分别为任务1,2的平均协方差矩阵。

2.构造空间滤波器

2.1 正交白化变换求白化特征矩阵P

由于混合空间协方差矩阵$R$是正定矩阵，由奇异值分解定理进行特征分解：
$$R=U\lambda U^T\text{\hspace{1em}(4)}$$
上式中，$U$是特征向量矩阵，$\lambda$为对应的特征值的对角阵，按特征值按降序排列，白化转换$U$可得：
$$P=\frac{1}{\sqrt{\lambda }}{U}^T\text{\hspace{1em}(5)}$$

2.2 构建空间滤波器

将矩阵$P$作用于$C_1$和$C_2$可得：
$$S_1=P{R}_1{P}^T,S_2=P{R}_2{P}^T\text{\hspace{1em}(6)}$$
$S_1$、$S_2$具有公共特征向量，且存在两个对角矩阵${\lambda }_1$、${\lambda }_2$和相同的特征向量矩阵$B$, 对$S_1$、$S_2$进行主分量分解，可得：
$$S_1=B{\lambda }_1{B}^T$$
$$S_2=B{\lambda }_2{B}^T\text{\hspace{1em}(7)}$$
且两个特征值的对角阵${\lambda }_1$和${\lambda }_2$之和为单位矩阵：
$${\lambda }_1+{\lambda }_2=I\text{\hspace{1em}(8)}$$
由上式(8)可知，若${\lambda }_1$中的特征值按照降序排列，则${\lambda }_2$中对应的特征值按升序排列。由于${\lambda }_1$、${\lambda }_2$为$S_1$、$S_2$的对角矩阵，所以对于特征向量矩阵$B$，当$S_1$有最大的特征值时，$S_2$具有最小的特征值。因此可以利用矩阵$B$实现两类问题的分类，由此得到投影矩阵$W$:
$$W=B^{T}P\text{\hspace{1em}(9)}$$
投影矩阵$W$就是对应的空间滤波器。

2.3 特征提取

将训练集的运动想象矩阵$X_L$和$X_R$经过滤波器$W$滤波可得特征$Z_L$和$Z_R$:
$$Z_L=W\times X_L$$
$$Z_R=W\times X_R\text{\hspace{1em}(10)}$$
对于测试数据$X_i$，其特征向量$f_i$提取方式如下，
$$\{\begin{array}{l}{Z}_i=W\times X_i\\ f_i=\frac{var(Z_i)}{\sum (var(Z_i))}\end{array}$$
将$f_i$与$f_L$和$f_R$进行比较以确定第$i$次想象为想象左还是想象右。根据CSP算法在多电极采集脑电信号特征提取的定义，其中$f_L$和$f_R$的定义如下：
$$f_L=\frac{var(Z_L)}{\sum (var(Z_L))}$$
$$f_R=\frac{var(Z_R)}{\sum (var(Z_R))}$$

Matlab实战

clc;
clear;
EEGSignals = load('graz_data/CSP_train.mat');   % 加载带通滤波后的脑电数据
%check and initializations
EEG_Channels = size(EEGSignals.x_train,2);
EEG_Trials = size(EEGSignals.x_train,3);
classLabels = unique(EEGSignals.y_train);% Return non-repeating values
EEG_Classes = length(classLabels);
covMatrix = cell(EEG_Classes,1); % 协方差矩阵
% Computing the normalized covariance matrices for each trial
trialCov = zeros(EEG_Channels,EEG_Channels,EEG_Trials);
for i = 1:EEG_Trials
    E = EEGSignals.x_train(:,:,i)';
    EE = E*E';
    trialCov(:,:,i) = EE./trace(EE);  % 计算协方差矩阵
end
clear E;
clear EE;
% 计算每一类样本数据的空间协方差之和
for i = 1:EEG_Classes
    covMatrix{i} = mean(trialCov(:,:,EEGSignals.y_train == classLabels(i)),3);
end
% 计算两类数据的空间协方差之和
covTotal = covMatrix{1} + covMatrix{2};
% 计算特征向量和特征矩阵
[Uc,Dt] = eig(covTotal);
% 特征值要降序排列
eigenvalues = diag(Dt);
[eigenvalues,egIndex] = sort(eigenvalues, 'descend');% 降序
Ut = Uc(:,egIndex);
% 矩阵白化
P = diag(sqrt(1./eigenvalues))*Ut';
% 矩阵P作用求公共特征向量transformedCov1 
transformedCov1 = P*covMatrix{1}*P';
%计算公共特征向量transformedCov1的特征向量和特征矩阵
[U1,D1] = eig(transformedCov1);
eigenvalues = diag(D1);
[eigenvalues,egIndex] = sort(eigenvalues, 'descend');% 降序排列
U1 = U1(:, egIndex);
% 计算投影矩阵W
CSPMatrix = U1' * P;
% 计算特征矩阵
FilterPairs = 2;       % CSP特征选择参数m    CSP特征为2*m个
features_train = zeros(EEG_Trials, 2*FilterPairs+1);
features_test = zeros(EEG_Trials, 2*FilterPairs+1);
Filter = CSPMatrix([1:FilterPairs (end-FilterPairs+1):end],:);
%extracting the CSP features from each trial
for t=1:EEG_Trials    
    %projecting the data onto the CSP filters    
    projectedTrial_train = Filter * EEGSignals.x_train(:,:,t)';    
    projectedTrial_test = Filter * EEGSignals.x_test(:,:,t)';
    %generating the features as the log variance of the projected signals
    variances_train = var(projectedTrial_train,0,2);  
    variances_test = var(projectedTrial_test,0,2);
    for f=1:length(variances_train)
        features_train(t,f) = log(variances_train(f));
        % features_train(t,f) = log(variances_train(f)/sum(variances_train));   %修改后对应公式
    end
    for f=1:length(variances_test)
        features_test(t,f) = log(variances_test(f));
        %features_test(t,f) = log(variances_test(f)/sum(variances_test));  % 修改后对应公式
    end
end
CSP_Train_feature = features_train(:,1:4);
CSP_Test_feature = features_test(:,1:4);
save('CSP_feature.mat','CSP_Train_feature','CSP_Test_feature');

代码来源于网络

脑际爱好者笔QQ交流群：903290195
更多分享，请关注公众号：