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

  • 2019 年 10 月 6 日
  • 笔记

共空间模式CSP

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

共空间模式理论

假设X_1X_2分别为两分类想象运动任务下的多通道诱发响应时-空信号矩阵,他们的维数均为N*T,N为脑电通道数,T为每个通道所采集的样本数。为了计算其协方差矩阵,现在假设N<T.在两种脑电想象任务情况下,一般采用复合源的数学模型描述EEG信号。为了简化计算,常忽略噪声所带来的影响。X_1X_2分别为:

上式(1)中,S_1S_2分别代表两种类型任务。假设两种信号源是相互线性独立的;S_M代表两种类型任务下所共同拥有的源信号,假设S_1是由m_1个源所构成的,S_1是由m_2个源所构成,则C_1C_2便是由S_1S_2相关的m_1m_2个共同空间模式组成的。由于每个空间模式都是一个N∗1维的向量,现在用该向量来表示单个的源信号所引起的信号在N个导联上的分布权重。C_M表示的是与S_M相应的共有的空间模式。CSP算法的目标就是要设计空间滤波器F_1F_2得到空间因子W

1.求解协方差矩阵

时空信号矩阵X_1X_2归一化的协方差矩阵R_1R_2:

R_i=frac{X_{i}X_{i}^{T}}{traceleft(X_{i}X_{i}^{T}right)} (i=1,2) tag{2}

上式(2)中,X^{T}表示矩阵X的转置,trace(X)表示对矩阵对角线上元素求和。之后求解混合空间的协方差矩阵R:

R=overline{R_1}+overline{R_2} tag{3}

上式中,overline{R_i}(i=1,2)分别为任务1,2的平均协方差矩阵。

2.构造空间滤波器

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

由于混合空间协方差矩阵R是正定矩阵,由奇异值分解定理进行特征分解:

R=Ulambda U^{T} tag{4}

上式中,U是特征向量矩阵,lambda为对应的特征值的对角阵,按特征值按降序排列,白化转换U可得:

P=frac{1}{sqrt{lambda}}U^{T} tag{5}

2.2 构建空间滤波器

将矩阵P作用于C_1C_2可得:

S_1=PR_{1}P^{T},S_2=PR_{2}P^{T} tag{6}

S_1S_2具有公共特征向量,且存在两个对角矩阵lambda_{1}lambda_{2}和相同的特征向量矩阵B, 对S_1S_2进行主分量分解,可得:

S_1=B lambda_{1}B^{T}

S_2=B lambda_{2}B^{T} tag{7}

且两个特征值的对角阵lambda_{1}lambda_{2}之和为单位矩阵:

lambda_{1}+lambda_{2}=I tag{8}

由上式(8)可知,若lambda_{1}中的特征值按照降序排列,则lambda_{2}中对应的特征值按升序排列。由于lambda_{1}lambda_{2}S_1S_2的对角矩阵,所以对于特征向量矩阵B,当S_1有最大的特征值时,S_2具有最小的特征值。因此可以利用矩阵B实现两类问题的分类,由此得到投影矩阵W:

W=B^{T}P tag{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} tag{10}

对于测试数据X_i,其特征向量f_i提取方式如下,

begin{cases} Z_i=W times X_i \ f_i=frac{var(Z_i)}{sum(var(Z_i))}end{cases} tag{11}

f_if_{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))} tag{12}

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

代码来源于网络

更多分享,请关注公众号