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

一、原理

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

公共空间模式（CSP）算法的过程如下：
（1）将原始数据按照类别进行分段。两类样本数据E可分类为$E_1$和$E_2$。$E_1$为第一类样本数据，$E_2$为第二类样本数据。
（2）计算分段后的原始数据的协方差矩阵，协方差矩阵的计算公式为：$$C_i=\frac{E_i\cdot E_i^{\mathrm{T}}}{trace(E_i\cdot E_i^{\mathrm{T}})}（i=1，2）$$ $trace(E)$表示求矩阵E的迹。
分别计算分段后原始数据的协方差矩阵， $C_1$为第一类样本数据的空间协方差矩阵的期望， $C_2$表示第二类样本的空间协方差矩阵的期望， $C_c$表示两类数据的空间协方差矩阵之和，则有$$C_c=C_1+C_2$$
(3）正交白化变换并且同时对角化（需理解矩阵白化的定义）
$C_c$是正定矩阵，由奇值分解定理则：
$$C_c=U_c{\Lambda }_c{U}_c^{\mathrm{T}}$$
$U_c$为特征向量矩阵，${\Lambda }_c$表示特征值的对角阵，且特征值是降序排列的。
白化转换$U_c$可得到：
$$P=\frac{1}{\sqrt{{\Lambda }_c}}\cdot U_c^{\mathrm{T}}$$
将矩阵$P$作用于$C_1、C_2$,可得到：
$$S_1=P{C}_1{P}^T，S_2=P{C}_2{P}^T$$
$S_1、S_2$具有公共特征向量，而且存在两个对角矩阵${\Lambda }_1、{\Lambda }_2$和特征向量矩阵$B$，满足如下条件：
$$S_1=B{\Lambda }_1{B}^T$$
$$S_2=B{\Lambda }_2{B}^T$$
$${\Lambda }_1+{\Lambda }_2=I$$
其中$I$为单位阵。由此可以看出$S_1$和$S_2$的特征值之和等于1。
（4）计算投影矩阵
对于特征向量矩阵$Q$,当一个类别$S_1$有最大的特征值时，此时另一个类别$S_2$有最小的特征值，因此可以利用矩阵Q实现两类问题的分类，可以得到投影矩阵：$$W=(Q^{T}P{)}^T$$
（5）经过投影得到特征矩阵
将原始脑电数据E_(M×N)通过投影矩阵W进行投影得到特征矩阵：$$Z_{M\times N}=W_{M\times M}{E}_{M\times N}$$
可选择$Z_{M\times N}$的前m行和后m行（2m<M）作为原始输入数据的特征
（6）特征归一化
$$y_i=\log (\frac{var(Z_i)}{{\sum }_{n=1}^{2m}var(Z_n)})$$
其中$y_i$ 为第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');