转自：https://blog.csdn.net/MissXy_/article/details/81264953
建议点击原文进行阅读，拷贝过来的时候很多公式出现问题，博主是为了复习方便，故而转载过来。。。。
共空间模式丨Common Spatial Pattern (CSP)
特征提取算法之一，在运动想象应用较多。
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 ~= 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;- 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
%
%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
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
