转自:https://blog.csdn.net/MissXy_/article/details/81264953

建议点击原文进行阅读,拷贝过来的时候很多公式出现问题,博主是为了复习方便,故而转载过来。。。。

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

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

    • 1 理论介绍

      • 1.1 求2类数据的混合空间协方差矩阵
      • 1.2 应用主成分分析法,求出白化特征值矩阵P
      • 1.3 构造空间滤波器
      • 1.4 特征提取
    • 2 MATLAB 实现

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

1 理论介绍

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

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

X1=[C1CM][S1SM],X2=[C2CM][S2SM](1)(1)X1=[C1CM][S1SM],X2=[C2CM][S2SM]

​  (1)式中:S1S1和S2S2分别代表两种类型任务。不妨假设这两个源信号是相互线性独立的;SMSM代表两种类型任务下所共同拥有的源信号,假设S1S1是由m1m1个源所构成的,S2S2是由m2m2个源所构成。则C1C1和C2C2便是由S1S1和S2S2相关的m1m1和m2m2个共同空间模式组成的,由于每个空间模式都是一个N∗1N∗1维的向量,现在用这个向量来表示单个的源信号所引起的信号在NN个导联上的分布权重。CMCM表示的是与SMSM相应的共有的空间模式。CSPCSP算法的目标就是要设计空间滤波器F1F1和F2F2得到空间因子WW。

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

  X1X1和X2X2归一化后的协方差矩阵R1R1和R2R2分别为:

R1=X1XT1trace(X1XT1),R2=X2XT2trace(X2XT2)(2)(2)R1=X1X1Ttrace(X1X1T),R2=X2X2Ttrace(X2X2T)

  ​ (2)式中:XTXT表示XX矩阵的转置,trace(X)trace(X)表示矩阵对角线上元素的和。然后求混合空间协方差矩阵RR:

R=R1¯+R2¯(3)(3)R=R1¯+R2¯

​   (3)式中:Ri¯(i=1,2)Ri¯(i=1,2)分别为任务1,2实验的平均协方差矩阵

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

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

R=UλUT(4)(4)R=UλUT

​   (4)式中:UU是矩阵λλ的特征向量矩阵,λλ是对应的特征值构成的对角阵。将特征值进行降序排列,白化值矩阵为:

P=λ−1−−−√UT(5)(5)P=λ−1UT

1.3 构造空间滤波器

  ​ 对R1R1和R2R2进行如下变换:

S1=PR1PT,S2=PR2PT(6)(6)S1=PR1PT,S2=PR2PT

  ​ 然后对S1S1和S2S2做主分量分解,得到:

S1=B1λ1BT1,S2=B2λ2BT2(7)(7)S1=B1λ1B1T,S2=B2λ2B2T

​   通过上面的式子可以证明矩阵S1S1的特征向量和矩阵S2S2的特征向量矩阵是相等的,即:

B1=B2=VB1=B2=V

​   与此同时,两个特征值的对角阵λ1λ1和λ2λ2之和为单位矩阵,即:

λ1+λ1=Iλ1+λ1=I

  ​ 由于两类矩阵的特征值相加总是为1,则S1S1的最大特征值所对应的特征向量使S2S2 有最小的特征值,反之亦然。

  ​ 把λ1λ1中的特征值按照降序排列,则λ2λ2中对应的特征值按升序排列,根据这点可以推断出λ1λ1和λ2λ2具有下面的形式:

λ1=diag(I1 σM 0),λ2=diag(0 σM I2)λ1=diag(I1 σM 0),λ2=diag(0 σM I2)

​   白化EEG到与λ1λ1和λ2λ2中的最大特征值对应的特征向量的变换对于分离两个信号矩阵中的方差是最佳的。投影矩阵WW是所对应的空间滤波器为:

W=BTPW=BTP

1.4 特征提取

  ​ 将训练集的运动想象矩阵XL、XRXL、XR 经过构造的相应滤波器WW 滤波可得特征ZL、ZRZL、ZR 为:

ZL=W×XLZR=W×XRZL=W×XLZR=W×XR

  ​ 根据 CSP 算法在多电极采集脑电信号特征提取的定义,本研究选取 fLfL 和 fRfR 为想象左和想象右的特征向量,定义如下:

fL=VAR(ZL)sum(VAR(ZL))fR=VAR(ZR)sum(VAR(ZR))fL=VAR(ZL)sum(VAR(ZL))fR=VAR(ZR)sum(VAR(ZR))

​   对于测试数据 XiXi 来说,其特征向量 fifi 提取方式如下,并与 fLfL 和 fRfR 进行比较以确定第ii 次想象为想左或想右。

{Zi=W×Xifi=VAR(Zi)sum(VAR(Zi)){Zi=W×Xifi=VAR(Zi)sum(VAR(Zi))

  ​ 求出测试的fifi 特征,进行下面的分类。

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 ~= 2disp('ERROR! CSP can only be used for two classes');return;
endcovMatrices = 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:nbTrialsE = EEGSignals.x(:,:,t)';                       %note the transposeEE = 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;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67

特征提取: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
%
%initializationsnbTrials = 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 projectedsignalsvariances = var(projectedTrial,0,2);    for f=1:length(variances)features(t,f) = log(1+variances(f));end
end

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

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

    共空间模式CSP 共空间模式CSP 共空间模式理论 1.求解协方差矩阵 2.构造空间滤波器 2.1 正交白化变换求白化特征矩阵P 2.2 构建空间滤波器 2.3 特征提取 Matlab实战 共空间模式 ...

  2. Common Spatial Pattern(CSP)共空间模式(包含Matlab代码)

    文章目录 前言 1.CSP算法 2.代码实现 2.1-learnCSPLagrangian.m 2.2-learnCSP.m 2.3-extractCSPFeatures.m 2.4-Test.m 2 ...

  3. CSP(Common spatial patterns)共空间模式算法简介

    CSP(Common spatial patterns)共空间模式算法简介 CSP算法基础知识 在了解CSP算法之前的相关知识基础 空间滤波器(spatial filters) 滤波是指通过操作接受或 ...

  4. 运动想象| EEG信号、共空间模式算法(CSP)

    摘要 作为一种特殊的人机交互模式,脑-机接口(brain-computer interface, BCI)技术成为了当前信息交互的研究热点.其中脑电信号(electroencephalography, ...

  5. 滤波器组共空间模式(Filter Bank Common Spatial Pattern,FBCSP)

    FBCSP matlab Matlab FBCSP实现 FBCSP原理 测试数据 文件结构结构如下 数据结构如下 运行结果 代码 主文件Test.m FBCSP.m函数 FBCSPOnline.m函数 ...

  6. 脑电信号特征提取常用算法(共空间模式CSP、小波变换DWT、功率谱密度PSD、AR模型)

    1 共空间模式CSP 原理:共空间模式(CSP)是一种对两分类任务下的空域滤波特征提取算法,能够从多通道的脑机接口数据里面提取出每一类的空间分布成分.公共空间模式算法的基本原理是利用矩阵的对角化,找到 ...

  7. CSP(共空间模式)的python实现

    数据来源(bbci竞赛数据graz两分类) 链接:https://pan.baidu.com/s/1GVPp2UBDiIUng6PlVU23uw 提取码:7jal 代码 def csp_yx(EEGS ...

  8. python实现共空间模式CSP

    直接调用库函数 mne.decoding.CSP(n_components=4, reg=None, log=None, cov_est='concat', transform_into='avera ...

  9. matlab csp共空间代码,运动想象中共空间模式算法(CSP)的实现

    最近在研究运动想象算法,其中CSP来提取特征用的比较多,尤为是在二分类的问题中,以前写过一篇如何在MNE库中实现CSP算法的博客,用的是MNE库中已经写好的算法,如今想本身实现该算法,研究了几天发现坑 ...

  10. 第19章 解释器模式(Interpreter Pattern)

    原文 第19章 解释器模式(Interpreter Pattern) 解释器模式 导读:解释器模式,平常用的比较的少,所以在写这个模式之前在博客园搜索了一番,看完之后那叫一个头大.篇幅很长,我鼓足了劲 ...

最新文章

  1. 编译可在Android上运行的qemu user mode
  2. Linux驱动之内核加载模块过程分析
  3. fluent的udf需要c语言环境吗,[转载]FLUENT UDF 使用指导
  4. 洛谷——P1042 乒乓球
  5. Linux中修改登录提示文件
  6. gcc是java的什么意思_为什么gcc支持Java而不是C#
  7. SQL SERVER data tier application 的作用及如何使用SSDT进行SQL数据库的自动化部署到生产环境和版本控制
  8. bzoj1085骑士精神(搜索)
  9. potplayer 多个进程_进程组、会话、控制终端概念,如何创建守护进程?
  10. python高阶_Python高阶学习
  11. symmetry methods for differential equations,exercise 1.4
  12. 有关i++问题,和一些另外的易错点
  13. JavaScript继承理解:ES5继承方式+ES6Class继承对比
  14. 从零开始的全栈工程师——html篇1.4
  15. 数学软件Maple使用教程
  16. 程序员修炼之道——读序
  17. 2021半年度博客总结
  18. java获取微信token_Java微信公众平台开发(六)--微信开发中的token获取
  19. 磁条卡,接触式IC卡,非接触式IC卡的优缺点
  20. Java面试必背八股文[11]:计算机网络

热门文章

  1. 传奇服务器端显示时间问题,架设传奇出现is not a valid date and time时间错误
  2. 费马小定理的两个证明
  3. 58赶集基于 Docker 的自动化部署实践
  4. mysql5.6 relay.info_Relay log 导致复制启动失败
  5. 计算机的应用数据处理,计算机的应用领域:数据处理(或信息处理)
  6. 保镖机器人作文_【保镖的作文】_玛雅作文网
  7. 六级考研单词之路-五十四
  8. [4G5G专题-8]:RRU 峰均比降低技术CFR(波峰系数削减)
  9. 大数据可视化技术应用学习目标与复习小结
  10. 大数据可视化(一)数据可视化概述