更多技术干货第一时间送达
共空间模式（Common Spatial Pattern, CSP）是一种对两分类任务下的空域滤波特征提取算法，能够从多通道的脑机接口数据里面提取出每一类的空间分布成分。公共空间模式算法的基本原理是利用矩阵的对角化，找到一组最优空间滤波器进行投影，使得两类信号的方差值差异最大化，从而得到具有较高区分度的特征向量。
假设和分别为两分类想象运动任务下的多通道诱发响应时-空信号矩阵，他们的维数均为,为脑电通道数，为每个通道所采集的样本数。为了计算其协方差矩阵，现在假设$N
上式中，和分别代表两种类型任务。假设两种信号源是相互线性独立的；代表两种类型任务下所共同拥有的源信号，假设是由个源所构成的，是由个源所构成，则和便是由和相关的和个共同空间模式组成的。由于每个空间模式都是一个维的向量，现在用该向量来表示单个的源信号所引起的信号在个导联上的分布权重。表示的是与相应的共有的空间模式。CSP算法的目标就是要设计空间滤波器和得到空间因子。
时空信号矩阵和归一化的协方差矩阵和:
上式(2)中，表示矩阵的转置，表示对矩阵对角线上元素求和。之后求解混合空间的协方差矩阵:
上式中，分别为任务1,2的平均协方差矩阵。
由于混合空间协方差矩阵是正定矩阵，由奇异值分解定理进行特征分解：
上式中，是特征向量矩阵，为对应的特征值的对角阵，按特征值按降序排列，白化转换可得：
将矩阵作用于和可得：
、具有公共特征向量，且存在两个对角矩阵、和相同的特征向量矩阵, 对、进行主分量分解，可得：
且两个特征值的对角阵和之和为单位矩阵：
由上式可知，若中的特征值按照降序排列，则中对应的特征值按升序排列。由于、为、的对角矩阵，所以对于特征向量矩阵，当有最大的特征值时，具有最小的特征值。因此可以利用矩阵实现两类问题的分类，由此得到投影矩阵:
投影矩阵就是对应的空间滤波器。
将训练集的运动想象矩阵和经过滤波器滤波可得特征和:
对于测试数据，其特征向量提取方式如下，
将与和进行比较以确定第次想象为想象左还是想象右。根据CSP算法在多电极采集脑电信号特征提取的定义，其中和的定义如下：
clc;clear;EEGSignals = load('graz_data/CSP_train.mat'); % 加载带通滤波后的脑电数据%check and initializationsEEG_Channels = size(EEGSignals.x_train,2);EEG_Trials = size(EEGSignals.x_train,3);classLabels = unique(EEGSignals.y_train);% Return non-repeating valuesEEG_Classes = length(classLabels);covMatrix = cell(EEG_Classes,1); % 协方差矩阵% Computing the normalized covariance matrices for each trialtrialCov = 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); % 计算协方差矩阵endclear 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);% 计算投影矩阵WCSPMatrix = 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 trialfor 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)); % 修改后对应公式 endendCSP_Train_feature = features_train(:,1:4);CSP_Test_feature = features_test(:,1:4);save('CSP_feature.mat','CSP_Train_feature','CSP_Test_feature');联系客服
