特征提取方法——共空间模式（Common Spatial Pattern,CSP）算法原理与实现

一、原理

公共空间模式（CSP）算法采用监督的方法创建一个最优的公共空间滤波器，在最大化一类方差的同时最小化另一类方差，采用同时对角化两类任务协方差矩阵的方式，得到可区分程度最大的特征向量。适用于二分类任务的特征提取。

公共空间模式（CSP）算法的过程如下：
（1）将原始数据按照类别进行分段。两类样本数据E可分类为 E 1 E_1 E1​和 E 2 E_2 E2​。 E 1 E_1 E1​为第一类样本数据， E 2 E_2 E2​为第二类样本数据。
（2）计算分段后的原始数据的协方差矩阵，协方差矩阵的计算公式为： C i = E i ⋅ E i T t r a c e ( E i ⋅ E i T ) （ i = 1 ， 2 ） C_i = \frac{E_i\cdot E_i ^ \mathrm{ T }}{trace(E_i\cdot E_i ^ \mathrm{ T })} （i =1，2） Ci​=trace(Ei​⋅EiT​)Ei​⋅EiT​​（i=1，2） t r a c e ( E ) trace(E) trace(E)表示求矩阵E的迹。
分别计算分段后原始数据的协方差矩阵， C 1 C_1 C1​为第一类样本数据的空间协方差矩阵的期望， C 2 C_2 C2​表示第二类样本的空间协方差矩阵的期望， C c C_c Cc​表示两类数据的空间协方差矩阵之和，则有 C c = C 1 + C 2 C_c = C_1+C_2 Cc​=C1​+C2​
(3）正交白化变换并且同时对角化（需理解矩阵白化的定义）
C c C_c Cc​是正定矩阵，由奇值分解定理则：
C c = U c Λ c U c T C_c= U_c Λ_c U_c^ \mathrm{ T } Cc​=Uc​Λc​UcT​
U c U_c Uc​为特征向量矩阵， Λ c Λ_c Λc​表示特征值的对角阵，且特征值是降序排列的。
白化转换 U c U_c Uc​可得到：
P = 1 Λ c ⋅ U c T P= \frac{1}{\sqrt{Λ_c}}\cdot U_c^ \mathrm{ T } P=Λc​ ​1​⋅UcT​
将矩阵 P P P作用于 C 1 、 C 2 C_1 、C_2 C1​、C2​,可得到：
S 1 = P C 1 P T ， S 2 = P C 2 P T S_1=PC_1 P^T，S_2=PC_2 P^T S1​=PC1​PT，S2​=PC2​PT
S 1 、 S 2 S_1 、S_2 S1​、S2​具有公共特征向量，而且存在两个对角矩阵 Λ 1 、 Λ 2 Λ_1、Λ_2 Λ1​、Λ2​和特征向量矩阵 B B B，满足如下条件：
S 1 = B Λ 1 B T S_1=BΛ_1 B^T S1​=BΛ1​BT
S 2 = B Λ 2 B T S_2=BΛ_2 B^T S2​=BΛ2​BT
Λ 1 + Λ 2 = I Λ_1+ Λ_2=I Λ1​+Λ2​=I
其中 I I I为单位阵。由此可以看出 S 1 S_1 S1​和 S 2 S_2 S2​的特征值之和等于1。
（4）计算投影矩阵
对于特征向量矩阵 Q Q Q,当一个类别 S 1 S_1 S1​有最大的特征值时，此时另一个类别 S 2 S_2 S2​有最小的特征值，因此可以利用矩阵Q实现两类问题的分类，可以得到投影矩阵： W = ( Q T P ) T W= (Q^T P)^T W=(QTP)T
（5）经过投影得到特征矩阵
将原始脑电数据E_(M×N)通过投影矩阵W进行投影得到特征矩阵： Z M × N = W M × M E M × N Z_{M\times N}= W_{M\times M} E_{M\times N} ZM×N​=WM×M​EM×N​
可选择 Z M × N Z_{M\times N} ZM×N​的前m行和后m行（2m<M）作为原始输入数据的特征
（6）特征归一化
y i = log ⁡ ( v a r ( Z i ) ∑ n = 1 2 m v a r ( Z n ) ) y_i = \log(\frac{var(Z_i)}{\sum_{n=1}^{2m}{var(Z_n)}}) yi​=log(∑n=12m​var(Zn​)var(Zi​)​)
其中 y i y_i yi​ 为第i个样本归一化后的特征矩阵。
算法生成的CSP特征矩阵，其信息不是等效的。特征信息主要集中在特征矩阵的头部和尾部，而中间的特征信息不明显可以忽略，所以选取m行和后m行数据作为CSP特征提取的特征矩阵。

二 、matlab算法实现

使用CSP算法处理二分类任务的运动想象脑电数据。下面是matlab代码，实现BCI2003运动想象数据集的脑电特征提取，提取的特征数为4个。

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');