特征提取丨共空间模式 Common Spatial Pattern (CSP)

共空间模式丨Common Spatial Pattern (CSP)

特征提取算法之一，在运动想象应用较多。

1 理论介绍

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

  ​ 假设$X_1$和$X_2$分别为两分类想象运动任务下的多通道诱发响应时-空信号矩阵，他们的维数均为$N*T$，$N$为脑电通道数，$T$为每个通道所采集的样本数。为了计算其协方差矩阵，现在假设$N＜T$。在两种脑电想象任务情况下，一般采用复合源的数学模型来描述$EEG$信号，为了方便计算，一般忽略噪声所产生的影响。$X_1$和$X_2$可以分别写成：

$$X_1=\begin{bmatrix}C_1C_M\\ \end{bmatrix}\begin{bmatrix}S_1\\S_M\\ \end{bmatrix},X_2=\begin{bmatrix}C_2C_M\\ \end{bmatrix}\begin{bmatrix}S_2\\S_M\\ \end{bmatrix} \tag{1}$$

​  （1）式中：$S_1$和$S_2$分别代表两种类型任务。不妨假设这两个源信号是相互线性独立的；$S_M$代表两种类型任务下所共同拥有的源信号，假设$S_1$是由$m_1$个源所构成的，$S_2$是由$m_2$个源所构成。则$C_1$和$C_2$便是由$S_1$和$S_2$相关的$m_1$和$m_2$个共同空间模式组成的，由于每个空间模式都是一个$N*1$维的向量，现在用这个向量来表示单个的源信号所引起的信号在$N$个导联上的分布权重。$C_M$表示的是与$S_M$相应的共有的空间模式。$CSP$算法的目标就是要设计空间滤波器$F_1$和$F_2$得到空间因子$W$。

1.1 求2类数据的混合空间协方差矩阵

  $X_1$和$X_2$归一化后的协方差矩阵$R_1$和$R_2$分别为：

$$R_1=\frac{X_1X_1^T}{trace(X_1X_1^T)},R_2=\frac{X_2X_2^T}{trace(X_2X_2^T)}\tag{2}$$

  ​ （2）式中：$X^T$表示$X$矩阵的转置，$trace(X)$表示矩阵对角线上元素的和。然后求混合空间协方差矩阵$R$:

$$R=\bar{R_1}+\bar{R_2}\tag{3}$$

​   （3）式中：$\bar{R_i}(i=1,2)$分别为任务1,2实验的平均协方差矩阵。

1.2 应用主成分分析法，求出白化特征值矩阵P

  对混合空间协方差矩阵$R$按式进行特征值分解:

$$R=U\lambda U^T\tag{4}$$

​   （4）式中：$U$是矩阵$\lambda$的特征向量矩阵，$λ$是对应的特征值构成的对角阵。将特征值进行降序排列，白化值矩阵为:

$$P=\sqrt{\lambda^{-1}}U^T\tag{5}$$

1.3 构造空间滤波器

  ​ 对$R_1$和$R_2$进行如下变换：

$$S_1=PR_1P^T,S_2=PR_2P^T\tag{6}$$

  ​ 然后对$S_1$和$S_2$做主分量分解，得到：

$$S_1=B_1\lambda_{1} B_1^T,S_2=B_2\lambda_2 B_2^T\tag{7}$$
​   通过上面的式子可以证明矩阵 $S_1$的特征向量和矩阵 $S_2$的特征向量矩阵是相等的，即：
$$B_1=B_2=V$$
​   与此同时，两个特征值的对角阵 $\lambda_1$和 $\lambda_2$之和为单位矩阵，即：
$$\lambda_1+\lambda_1=I$$
  ​ 由于两类矩阵的特征值相加总是为1，则 $S_1$的最大特征值所对应的特征向量使 $S_2$ 有最小的特征值，反之亦然。

  ​ 把$\lambda_1$中的特征值按照降序排列，则$\lambda_2$中对应的特征值按升序排列，根据这点可以推断出$\lambda_1$和$\lambda_2$具有下面的形式：

$$\lambda_1=diag(I_1 \space \sigma_M \space 0), \lambda_2=diag( 0 \space \sigma_M \space I_2)$$

​   白化EEG到与$\lambda_1$和$\lambda_2$中的最大特征值对应的特征向量的变换对于分离两个信号矩阵中的方差是最佳的。投影矩阵$W$是所对应的空间滤波器为：

$$W=B^TP$$

1.4 特征提取

  ​ 将训练集的运动想象矩阵$X_L 、X_R$ 经过构造的相应滤波器$W$ 滤波可得特征$Z_L 、Z_R$ 为:

$$Z_L=W\times X_L\\ Z_R=W\times X_R$$
  ​ 根据 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))}$$
​   对于测试数据 $X_i$ 来说，其特征向量 $f_i$ 提取方式如下，并与 $f_L$ 和 $f_R$ 进行比较以确定第 $i$ 次想象为想左或想右。
$$\left\{\begin{array}{ll} Z_i=W\times X_i\\ f_i=\frac{VAR(Z_i)}{sum(VAR(Z_i))} \end{array} \right.$$
  ​ 求出测试的 $f_i$ 特征，进行下面的分类。

2 MATLAB 实现

构造空间滤波器：learnCSP.m

function CSPMatrix = learnCSP(EEGSignals,classLabels)
%
%Input:
%EEGSignals: the training EEG signals, composed of 2 classes. These signals
%are a structure such that:
%   EEGSignals.x: the EEG signals as a [Ns * Nc * Nt] Matrix where
%       Ns: number of EEG samples per trial
%       Nc: number of channels (EEG electrodes)
%       nT: number of trials
%   EEGSignals.y: a [1 * Nt] vector containing the class labels for each trial
%   EEGSignals.s: the sampling frequency (in Hz)
%
%Output:
%CSPMatrix: the learnt CSP filters (a [Nc*Nc] matrix with the filters as rows)
%
%See also: extractCSPFeatures

%check and initializations
nbChannels = size(EEGSignals.x,2);      % 通道
nbTrials = size(EEGSignals.x,3);        % 实验次数
nbClasses = length(classLabels);        % 类别

if nbClasses ~= 2
    disp('ERROR! CSP can only be used for two classes');
    return;
end

covMatrices = cell(nbClasses,1); %the covariance matrices for each class

%% Computing the normalized covariance matrices for each trial
%% 为每个试验计算标准化的协方差矩阵。
trialCov = zeros(nbChannels,nbChannels,nbTrials);
for t=1:nbTrials
    E = EEGSignals.x(:,:,t)';                       %note the transpose
    EE = E * E';
    trialCov(:,:,t) = EE ./ trace(EE);
end
clear E;
clear EE;

%computing the covariance matrix for each class
for c=1:nbClasses      
    %EEGSignals.y==classLabels(c) returns the indeces corresponding to the class labels 
    covMatrices{c} = mean(trialCov(:,:,EEGSignals.y == classLabels(c)),3);  
end

%the total covariance matrix
covTotal = covMatrices{1} + covMatrices{2};

%whitening transform of total covariance matrix
%caution: the eigenvalues are initially in increasing order注意:特征值最初是递增的
[Ut Dt] = eig(covTotal); 
eigenvalues = diag(Dt);
[eigenvalues egIndex] = sort(eigenvalues, 'descend');
Ut = Ut(:,egIndex);
P = diag(sqrt(1./eigenvalues)) * Ut';

%transforming covariance matrix of first class using P 
%用P变换第一类协方差矩阵
transformedCov1 =  P * covMatrices{1} * P';

%EVD of the transformed covariance matrix 变换协方差矩阵的EVD
[U1 D1] = eig(transformedCov1);
eigenvalues = diag(D1);
[eigenvalues egIndex] = sort(eigenvalues, 'descend');
U1 = U1(:, egIndex);
CSPMatrix = U1' * P;

特征提取：extractCSP.m

function features = extractCSP(EEGSignals, CSPMatrix, nbFilterPairs)
%
%extract features from an EEG data set using the Common Spatial Patterns (CSP) algorithm
%
%Input:
%EEGSignals: the EEGSignals from which extracting the CSP features. These signals
%are a structure such that:
%   EEGSignals.x: the EEG signals as a [Ns * Nc * Nt] Matrix where
%       Ns: number of EEG samples per trial
%       Nc: number of channels (EEG electrodes)
%       nT: number of trials
%   EEGSignals.y: a [1 * Nt] vector containing the class labels for each trial
%   EEGSignals.s: the sampling frequency (in Hz)
%CSPMatrix: the CSP projection matrix, learnt previously (see function learnCSP)
%nbFilterPairs: number of pairs of CSP filters to be used. The number of
%   features extracted will be twice the value of this parameter. The
%   filters selected are the one corresponding to the lowest and highest
%   eigenvalues
%
%Output:
%features: the features extracted from this EEG data set 
%   as a [Nt * (nbFilterPairs*2 + 1)] matrix, with the class labels as the
%   last column   
%
%initializations

nbTrials = size(EEGSignals.x,3);
features = zeros(nbTrials, 2*nbFilterPairs+1);
Filter = CSPMatrix([1:nbFilterPairs (end-nbFilterPairs+1):end],:);

%extracting the CSP features from each trial
for t=1:nbTrials    
    %projecting the data onto the CSP filters    
    projectedTrial = Filter * EEGSignals.x(:,:,t)';    

    %generating the features as the log variance of the projectedsignals

    variances = var(projectedTrial,0,2);    
    for f=1:length(variances)
        features(t,f) = log(1+variances(f));
    end
end