共空間模式 Common Spatial Pattern(CSP)原理和實戰
- 2019 年 10 月 6 日
- 筆記
共空間模式CSP
共空間模式（Common Spatial Pattern, CSP）是一種對兩分類任務下的空域濾波特徵提取算法，能夠從多通道的腦機接口數據裏面提取出每一類的空間分佈成分。公共空間模式算法的基本原理是利用矩陣的對角化，找到一組最優空間濾波器進行投影，使得兩類信號的方差值差異最大化，從而得到具有較高區分度的特徵向量。
共空間模式理論
假設X_1和X_2分別為兩分類想像運動任務下的多通道誘發響應時-空信號矩陣，他們的維數均為N*T,N為腦電通道數，T為每個通道所採集的樣本數。為了計算其協方差矩陣，現在假設N<T.在兩種腦電想像任務情況下，一般採用複合源的數學模型描述EEG信號。為了簡化計算，常忽略噪聲所帶來的影響。X_1和X_2分別為：
上式(1)中，S_1和S_2分別代表兩種類型任務。假設兩種信號源是相互線性獨立的；S_M代表兩種類型任務下所共同擁有的源信號，假設S_1是由m_1個源所構成的，S_1是由m_2個源所構成，則C_1和C_2便是由S_1和S_2相關的m_1和m_2個共同空間模式組成的。由於每個空間模式都是一個N∗1維的向量，現在用該向量來表示單個的源信號所引起的信號在N個導聯上的分佈權重。C_M表示的是與S_M相應的共有的空間模式。CSP算法的目標就是要設計空間濾波器F_1和F_2得到空間因子W。
1.求解協方差矩陣
時空信號矩陣X_1和X_2歸一化的協方差矩陣R_1和R_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_1和C_2可得：
S_1=PR_{1}P^{T},S_2=PR_{2}P^{T} tag{6}
S_1、S_2具有公共特徵向量，且存在兩個對角矩陣lambda_{1}、lambda_{2}和相同的特徵向量矩陣B, 對S_1、S_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_1、S_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_i與f_{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');
代碼來源於網絡
