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

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

1 理论介绍

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

  ​ 假设 X1 X 1 <script type="math/tex" id="MathJax-Element-1">X_1</script>和 X2 X 2 <script type="math/tex" id="MathJax-Element-2">X_2</script>分别为两分类想象运动任务下的多通道诱发响应时-空信号矩阵，他们的维数均为 N∗T N ∗ T <script type="math/tex" id="MathJax-Element-3">N*T</script>， N N <script type="math/tex" id="MathJax-Element-4">N</script>为脑电通道数， T <script type="math/tex" id="MathJax-Element-5">T</script>为每个通道所采集的样本数。为了计算其协方差矩阵，现在假设 N＜T N ＜ T <script type="math/tex" id="MathJax-Element-6">N＜T</script>。在两种脑电想象任务情况下，一般采用复合源的数学模型来描述 EEG E E G <script type="math/tex" id="MathJax-Element-7">EEG</script>信号，为了方便计算，一般忽略噪声所产生的影响。 X1 X 1 <script type="math/tex" id="MathJax-Element-8">X_1</script>和 X2 X 2 <script type="math/tex" id="MathJax-Element-9">X_2</script>可以分别写成：

X1=[C1CM][S1SM],X2=[C2CM][S2SM](1) (1) X 1 = [ C 1 C M ] [ S 1 S M ] , X 2 = [ C 2 C M ] [ S 2 S M ]
<script type="math/tex; mode=display" id="MathJax-Element-10"> 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} </script>

​  （1）式中： S1 S 1 <script type="math/tex" id="MathJax-Element-11">S_1</script>和 S2 S 2 <script type="math/tex" id="MathJax-Element-12">S_2</script>分别代表两种类型任务。不妨假设这两个源信号是相互线性独立的； SM S M <script type="math/tex" id="MathJax-Element-13">S_M</script>代表两种类型任务下所共同拥有的源信号，假设 S1 S 1 <script type="math/tex" id="MathJax-Element-14">S_1</script>是由 m1 m 1 <script type="math/tex" id="MathJax-Element-15">m_1</script>个源所构成的， S2 S 2 <script type="math/tex" id="MathJax-Element-16">S_2</script>是由 m2 m 2 <script type="math/tex" id="MathJax-Element-17">m_2</script>个源所构成。则 C1 C 1 <script type="math/tex" id="MathJax-Element-18">C_1</script>和 C2 C 2 <script type="math/tex" id="MathJax-Element-19">C_2</script>便是由 S1 S 1 <script type="math/tex" id="MathJax-Element-20">S_1</script>和 S2 S 2 <script type="math/tex" id="MathJax-Element-21">S_2</script>相关的 m1 m 1 <script type="math/tex" id="MathJax-Element-22">m_1</script>和 m2 m 2 <script type="math/tex" id="MathJax-Element-23">m_2</script>个共同空间模式组成的，由于每个空间模式都是一个 N∗1 N ∗ 1 <script type="math/tex" id="MathJax-Element-24">N*1</script>维的向量，现在用这个向量来表示单个的源信号所引起的信号在 N N <script type="math/tex" id="MathJax-Element-25">N</script>个导联上的分布权重。 C M <script type="math/tex" id="MathJax-Element-26">C_M</script>表示的是与 SM S M <script type="math/tex" id="MathJax-Element-27">S_M</script>相应的共有的空间模式。 CSP C S P <script type="math/tex" id="MathJax-Element-28">CSP</script>算法的目标就是要设计空间滤波器 F1 F 1 <script type="math/tex" id="MathJax-Element-29">F_1</script>和 F2 F 2 <script type="math/tex" id="MathJax-Element-30">F_2</script>得到空间因子 W W <script type="math/tex" id="MathJax-Element-31">W</script>。

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

   X 1 <script type="math/tex" id="MathJax-Element-32">X_1</script>和 X2 X 2 <script type="math/tex" id="MathJax-Element-33">X_2</script>归一化后的协方差矩阵 R1 R 1 <script type="math/tex" id="MathJax-Element-34">R_1</script>和 R2 R 2 <script type="math/tex" id="MathJax-Element-35">R_2</script>分别为：

R1=X1XT1trace(X1XT1),R2=X2XT2trace(X2XT2)(2) (2) R 1 = X 1 X 1 T t r a c e ( X 1 X 1 T ) , R 2 = X 2 X 2 T t r a c e ( X 2 X 2 T )
<script type="math/tex; mode=display" id="MathJax-Element-36"> 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} </script>

  ​ （2）式中： XT X T <script type="math/tex" id="MathJax-Element-37">X^T</script>表示 X X <script type="math/tex" id="MathJax-Element-38">X</script>矩阵的转置， t r a c e ( X ) <script type="math/tex" id="MathJax-Element-39">trace(X)</script>表示矩阵对角线上元素的和。然后求混合空间协方差矩阵 R R <script type="math/tex" id="MathJax-Element-40">R</script>:

(3) R = R 1 ¯ + R 2 ¯
<script type="math/tex; mode=display" id="MathJax-Element-41"> R=\bar{R_1}+\bar{R_2}\tag{3} </script>

​   （3）式中： Ri¯(i=1,2) R i ¯ ( i = 1 , 2 ) <script type="math/tex" id="MathJax-Element-42">\bar{R_i}(i=1,2)</script>分别为任务1,2实验的平均协方差矩阵。

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

  对混合空间协方差矩阵 R R <script type="math/tex" id="MathJax-Element-43">R</script>按式进行特征值分解:

(4) R = U λ U T
<script type="math/tex; mode=display" id="MathJax-Element-44"> R=U\lambda U^T\tag{4} </script>

​   （4）式中： U U <script type="math/tex" id="MathJax-Element-45">U</script>是矩阵 λ <script type="math/tex" id="MathJax-Element-46">\lambda</script>的特征向量矩阵， λ λ <script type="math/tex" id="MathJax-Element-47">λ</script>是对应的特征值构成的对角阵。将特征值进行降序排列，白化值矩阵为:

(5) P = λ − 1 U T
<script type="math/tex; mode=display" id="MathJax-Element-48"> P=\sqrt{\lambda^{-1}}U^T\tag{5} </script>

1.3 构造空间滤波器

  ​ 对 R1 R 1 <script type="math/tex" id="MathJax-Element-49">R_1</script>和 R2 R 2 <script type="math/tex" id="MathJax-Element-50">R_2</script>进行如下变换：

S1=PR1PT,S2=PR2PT(6) (6) S 1 = P R 1 P T , S 2 = P R 2 P T
<script type="math/tex; mode=display" id="MathJax-Element-51"> S_1=PR_1P^T,S_2=PR_2P^T\tag{6} </script>

  ​ 然后对 S1 S 1 <script type="math/tex" id="MathJax-Element-52">S_1</script>和 S2 S 2 <script type="math/tex" id="MathJax-Element-53">S_2</script>做主分量分解，得到：

S1=B1λ1BT1,S2=B2λ2BT2(7) (7) S 1 = B 1 λ 1 B 1 T , S 2 = B 2 λ 2 B 2 T
<script type="math/tex; mode=display" id="MathJax-Element-54"> S_1=B_1\lambda_{1} B_1^T,S_2=B_2\lambda_2 B_2^T\tag{7} </script>
​   通过上面的式子可以证明矩阵 S1 S 1 <script type="math/tex" id="MathJax-Element-55">S_1</script>的特征向量和矩阵 S2 S 2 <script type="math/tex" id="MathJax-Element-56">S_2</script>的特征向量矩阵是相等的，即：
B1=B2=V B 1 = B 2 = V
<script type="math/tex; mode=display" id="MathJax-Element-57"> B_1=B_2=V </script>
​   与此同时，两个特征值的对角阵 λ1 λ 1 <script type="math/tex" id="MathJax-Element-58">\lambda_1</script>和 λ2 λ 2 <script type="math/tex" id="MathJax-Element-59">\lambda_2</script>之和为单位矩阵，即：
λ1+λ1=I λ 1 + λ 1 = I
<script type="math/tex; mode=display" id="MathJax-Element-60"> \lambda_1+\lambda_1=I </script>
  ​ 由于两类矩阵的特征值相加总是为1，则 S1 S 1 <script type="math/tex" id="MathJax-Element-61">S_1</script>的最大特征值所对应的特征向量使 S2 S 2 <script type="math/tex" id="MathJax-Element-62">S_2</script> 有最小的特征值，反之亦然。

  ​ 把 λ1 λ 1 <script type="math/tex" id="MathJax-Element-63">\lambda_1</script>中的特征值按照降序排列，则 λ2 λ 2 <script type="math/tex" id="MathJax-Element-64">\lambda_2</script>中对应的特征值按升序排列，根据这点可以推断出 λ1 λ 1 <script type="math/tex" id="MathJax-Element-65">\lambda_1</script>和 λ2 λ 2 <script type="math/tex" id="MathJax-Element-66">\lambda_2</script>具有下面的形式：

λ1=diag(I1 σM 0),λ2=diag(0 σM I2) λ 1 = d i a g ( I 1   σ M   0 ) , λ 2 = d i a g ( 0   σ M   I 2 )
<script type="math/tex; mode=display" id="MathJax-Element-67"> \lambda_1=diag(I_1 \space \sigma_M \space 0), \lambda_2=diag( 0 \space \sigma_M \space I_2) </script>

​   白化EEG到与 λ1 λ 1 <script type="math/tex" id="MathJax-Element-68">\lambda_1</script>和 λ2 λ 2 <script type="math/tex" id="MathJax-Element-69">\lambda_2</script>中的最大特征值对应的特征向量的变换对于分离两个信号矩阵中的方差是最佳的。投影矩阵 W W <script type="math/tex" id="MathJax-Element-70">W</script>是所对应的空间滤波器为：

W = B T P
<script type="math/tex; mode=display" id="MathJax-Element-71"> W=B^TP </script>

1.4 特征提取

  ​ 将训练集的运动想象矩阵 XL、XR X L 、 X R <script type="math/tex" id="MathJax-Element-72">X_L 、X_R</script> 经过构造的相应滤波器 W W <script type="math/tex" id="MathJax-Element-73">W</script> 滤波可得特征 Z L 、 Z R <script type="math/tex" id="MathJax-Element-74">Z_L 、Z_R</script> 为:

ZL=W×XLZR=W×XR Z L = W × X L Z R = W × X R
<script type="math/tex; mode=display" id="MathJax-Element-75"> Z_L=W\times X_L\\ Z_R=W\times X_R </script>
  ​ 根据 CSP 算法在多电极采集脑电信号特征提取的定义，本研究选取 fL f L <script type="math/tex" id="MathJax-Element-76">f_L</script> 和 fR f R <script type="math/tex" id="MathJax-Element-77">f_R</script> 为想象左和想象右的特征向量，定义如下:
fL=VAR(ZL)sum(VAR(ZL))fR=VAR(ZR)sum(VAR(ZR)) f L = V A R ( Z L ) s u m ( V A R ( Z L ) ) f R = V A R ( Z R ) s u m ( V A R ( Z R ) )
<script type="math/tex; mode=display" id="MathJax-Element-78"> f_L=\frac{VAR(Z_L)}{sum(VAR(Z_L))}\\ f_R=\frac{VAR(Z_R)}{sum(VAR(Z_R))} </script>
​   对于测试数据 Xi X i <script type="math/tex" id="MathJax-Element-79">X_i</script> 来说，其特征向量 fi f i <script type="math/tex" id="MathJax-Element-80">f_i</script> 提取方式如下，并与 fL f L <script type="math/tex" id="MathJax-Element-81">f_L</script> 和 fR f R <script type="math/tex" id="MathJax-Element-82">f_R</script> 进行比较以确定第 i i <script type="math/tex" id="MathJax-Element-83"> i</script> 次想象为想左或想右。
{ Z i = W × X i f i = V A R ( Z i ) s u m ( V A R ( Z i ) )
<script type="math/tex; mode=display" id="MathJax-Element-84"> \left\{\begin{array}{ll} Z_i=W\times X_i\\ f_i=\frac{VAR(Z_i)}{sum(VAR(Z_i))} \end{array} \right. </script>
  ​ 求出测试的 fi f i <script type="math/tex" id="MathJax-Element-85">f_i</script> 特征，进行下面的分类。

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
Logo
2048 AI社区

有“AI”的1024 = 2048，欢迎大家加入2048 AI社区

更多推荐

cover

Electron实战——从HTML到桌面应用(1)

avatar 2048 AI社区
cover

9款免费AI论文工具实测：维普查重一把过，零AIGC痕迹超省心！

avatar 2048 AI社区

2025年技术总结 | 在Unity游戏开发路上的持续探索与沉淀

摘要： 2025年，向宇系统构建了Unity技术专栏体系，涵盖网络开发、插件推荐、特效实现等核心领域，形成完整学习路径。技术认知实现三大转变：从功能实现转向架构设计、从单一技术拓展至生态融合、从个人开发升级为知识分享。年度突破包括网络同步优化、Compute Shader特效性能提升及自定义渲染管线实践。通过主题式学习与时间管理策略克服挑战，2026年计划深入DOTS技术、AI集成等方向，并升级专

avatar 2048 AI社区