#?分析mne中CSP的計(jì)算，為自己實(shí)現(xiàn)CSP算法做準(zhǔn)備
#?協(xié)方差矩陣計(jì)算，這里有幾個(gè)點(diǎn)需要注意下，這里使用的是?Optimizing?Spatial?Filters?for?Robust?EEG?Single-Trial?Analysis?文章中的實(shí)現(xiàn)流程。
import?pickle
import?numpy?as?np
import?scipy.io?as?io
from?mne.decoding?import?CSP
import?matplotlib.pyplot?as?plt

import?mne

mne.set_log_level（False）

with?open（r‘F:\BaiduNetdiskDownload\BCICompetition\BCICIV_2b_gdf\Train\CueLeftRight\trainData.pkl‘‘rb‘）?as?f:
????epochsTrain?=?pickle.load（f）
trainData?=?epochsTrain.get_data（[‘EEG:C3‘?‘EEG:C4‘]）
trainLabels?=?epochsTrain.events[:?-1]



csp?=?CSP（n_components=2?reg=None?cov_est=‘epoch‘?norm_trace=False）
dataCSP?=?csp.fit_transform（trainData?trainLabels）
print（dataCSP.shape）
plt.figure（0）
plt.plot（dataCSP）



EEG_Channels?=?trainData.shape[1]
EEG_Trials?=?trainData.shape[0]
classLabels?=?np.unique（trainLabels）
EEG_Classes?=?len（classLabels）
covMatrix?=?list（range（EEG_Classes））
trialCov?=?np.zeros（[EEG_Channels?EEG_Channels?EEG_Trials]）
for?i?in?range（EEG_Trials）:
????E?=?trainData[i]
????EE?=?np.dot（E?E.T）
????trialCov[:?:?i]?=?EE/np.trace（EE）

for?i?in?range（EEG_Classes）:
????covMatrix[i]?=?np.mean（trialCov[:?:?trainLabels?==?classLabels[i]]?2）

print（‘a(chǎn)‘?covMatrix[0]）
print（covMatrix[1]）
covTotal?=?covMatrix[0]?+?covMatrix[1]

#########特征值降序
eigenvalues?featurevectors?=?np.linalg.eig（covTotal）
egIndex?=?np.argsort（eigenvalues）[::-1]
eigenvalues?=?np.sort（eigenvalues）[::-1]
Ut?=?featurevectors[:?egIndex]
#?矩陣白化
P?=?np.dot（np.diag（np.sqrt（1/eigenvalues））?Ut.T）
#?矩陣P作用求公共特征向量transformedCov1
transformedCov1?=?P.dot（covMatrix[0]）.dot（P.T）
#?計(jì)算公共特征向量transformedCov1的特征向量和特征矩陣
eigenvalues?featurevectors?=?np.linalg.eig（transformedCov1）
egIndex?=?np.argsort（eigenvalues）[::-1]
eigenvalues?=?np.sort（eigenvalues）[::-1]
U1?=?featurevectors[:?egIndex]
#?計(jì)算投影矩陣W
CSPMatrix?=?U1.T.dot（P）

#?計(jì)算特征矩陣
FilterPairs?=?2
feature_train?=?np.zeros（[EEG_Trials?2*FilterPairs+1]）
feature_test?=?np.zeros（[EEG_Trials?2*FilterPairs+1]）
Filter?=?np.r_[CSPMatrix[:FilterPairs]?CSPMatrix[-FilterPairs:]]

#從每一次實(shí)驗(yàn)中提取CSP特征
for?t?in?range（EEG_Trials）:
????#?將數(shù)據(jù)投影到CSP空間
????projectedTrial_train?=?Filter.dot（trainData[t]）
????projectedTrial_test?=?Filter.dot（trainData[t]）
????#?投影后信號(hào)的方差作為特征
????variances_train?=?np.var（projectedTrial_train?1）
????variances_test?=?np.var（projectedTrial_test?1）
????for?f?in?range（len（variances_train））:
????????feature_train[t?f]?=?np.log（variances_train[f]）

????for?f?in?range（len（variances_test））:
????????feature_test[t?f]?=?np.log（variances_test[f]）



CSP_Train_faetures?=?feature_train[:?0:4]
CSP_Test_features?=?feature_test[:?0:4]



plt.figure（1）
plt.plot（CSP_Train_faetures[:?:2]）

plt.show（）