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

数据来源（bbci竞赛数据graz两分类）

链接:https://pan.baidu.com/s/1GVPp2UBDiIUng6PlVU23uw
提取码:7jal

代码

def csp_yx(EEGSignals_x_train,EEGSignals_y_train):
    
    EEG_Channels = 3
    EEG_Trials = 280
    EEG_Classes = 2
    class_labels = [1, 2]
    
    
    
    # 计算协方差矩阵 cov1 cov2
    trialCov = np.zeros((EEG_Channels,EEG_Channels,EEG_Trials))
    for i in range(280):
        E = EEGSignals_x_train[:,:,i].T
    #     print(E.shape)
        EE = np.dot(E,E.T)
        trialCov[:,:,i] = EE/np.trace(EE);    

    
    
    cov1 = np.zeros((EEG_Channels,EEG_Channels,140))
    cov2 = np.zeros((EEG_Channels,EEG_Channels,140))
    
    MatrixCov = np.zeros((3,3,2))
    k=0
    s=0
    for i in range(EEG_Channels):
        if EEGSignals_y_train[i] == 1:
            cov1[:,:,k] = trialCov[:,:,i]
            k=k+1
        if EEGSignals_y_train[i] == 2:
            cov2[:,:,s] = trialCov[:,:,i]
            s=s+1
            
    np.set_printoptions(precision=4)
    
    # 计算平均协方差矩阵之和 covTotal
    cov1 = np.mean(cov1,2)
    cov2 = np.mean(cov2,2)
    covTotal = cov1 + cov2
    print(covTotal)
    
    # 计算公共特征向量矩阵和特征值
    [Dt,Uc] = np.linalg.eigh(covTotal)
    # 降序排序
    Uc = Uc[:,Dt.argsort()[::-1]]
    Dt = sorted(Dt,reverse=True)
    
    # 矩阵白化
    Dt = np.diag(Dt)
    Dt = np.sqrt(1./Dt)
  
    # 去 inf值
    c = np.isinf(Dt)
    Dt[c] = 0
    print(Dt)
    P = np.dot(Dt,Uc.T)
    
#     print(P)
   
    # 将P作用于 cov1 cov2
    transformedCov1 = np.dot(np.dot(P,cov1),P.T)
    transformedCov2 = np.dot(np.dot(P,cov2),P.T)  # S1 = P*R1*P.T
   

    # 将 transformedCov1 按主分量分解得到公共特征向量矩阵 B
    
    [D1,U1] = np.linalg.eig(transformedCov1)
    # 降序
    U1 = U1[:,D1.argsort()]
    D1 = sorted(D1,reverse=True)

#     print(U1)
    
    [D2,U2] = np.linalg.eig(transformedCov2)
    # 升序
    U2 = U2[:,D2.argsort()]
    D2 = np.sort(D2)
#     print(U2)

    
    # 计算投影矩阵
    CSPMatrix = np.dot(U1.T,P)
    
    return CSPMatrix

import scipy.io as sio
import numpy as np
dataSetfile = "../L-R hand/dataset_BCIcomp1.mat"
testLabelfile = "../test/labels_data_set_iii.mat"
EEGSignals = sio.loadmat(dataSetfile)
LabelSet = sio.loadmat(testLabelfile)
EEGSignals_x_train = np.concatenate([EEGSignals['x_train'],EEGSignals['x_test']], axis=2)
EEGSignals_y_train = np.concatenate([EEGSignals['y_train'],LabelSet['y_test']])

W = csp_yx.csp_yx(EEGSignals_x_train,EEGSignals_y_train)

FilterPairs = 1       # CSP特征选择参数m
features_train = np.zeros((280, 2*FilterPairs))
# features_test = np.zeros((80, 2*FilterPairs))

m1 = W[0:FilterPairs]
m2 = W[3-FilterPairs:3]

Filter = np.concatenate([m1,m2])
Filter

data = np.concatenate([EEGSignals['x_train'],EEGSignals['x_test']],axis=2)

projectedTrial = np.zeros((2*FilterPairs,1152))
EEG_Trials = 280
for i in range(EEG_Trials):
    # 将源信号投影到 CSP空间上
    projectedTrial = np.dot(Filter,data[:,:,i].T)
        
    # 将投影信号的方差得对数作为特征
    variances_train = np.var(projectedTrial,1)
    for f in range(len(variances_train)):
        features_train[i,f] = np.log(variances_train[f])

CSP = features_train[0,:2]