点击上面"脑机接口社区"关注我们
更多技术干货第一时间送达
事件相关去同步化与同步化(ERD/S)和运动相关皮质电位(MRCP)在下肢康复的脑机接口(BCI)中，特别是在站立和坐姿中，起着重要的作用。然而，人们对站立和坐着的大脑皮层活动的差异知之甚少，尤其是大脑的意图是如何调节运动前的感觉运动节奏的。在本研究中，研究人员旨在研究在站立和坐着的动作观察(AO)、运动想象(MI)和运动执行(ME) 期间连续性EEG节奏的解码。研究人员开发了一项行为任务，在该任务中，参与者被指示对坐立和站坐的动作执行AO和MI/ME。实验结果表明，在AO期间ERD比较显著，而在MI期间ERS在感觉运动区域的alpha带较为典型。结合常用空间模式(FBCSP)和支持向量机(SVM)进行离线和分类器测试分析。离线分析表明，AO和MI的分类在站-坐转换时的平均准确率最高，为82.73±2.54%。通过分类器测试分析，研究人员证明了MI范式比ME范式具有更高的解码神经意图的性能。
为了研究在连续脑电图记录下的运动执行过程中解码MI信号(包括ERD/S)和MRCPs的可行性，整个实验过程由MI和ME两个阶段组成。每一阶段包括3次运行过程(每次5次试验)，共包含30次试验。实验以坐姿开始，然后重复5次坐立和站坐交替试验。图1显示了每次试验中四个状态的序列:R、AO、idle和任务执行状态(MI或ME)。在R状态期间，显示器上显示黑色屏幕长达6秒，参与者被要求保持放松和静止。为了避免指令的模糊性，我们提供了持续4 ~ 5秒的坐-站或站-坐视频任务的视频刺激来指导参与者处于AO状态。参与者被要求在听到音频提示后立即完成两个阶段的任务。在ME阶段，参与者被要求在音频提示后完成一系列自主节奏的动作。在MI阶段，参与者在听到音频提示后开始想象指定的动作。在MI期间，可以通过音频或视觉提示获得运动启动，而在ME期间，可以从EMG信号获得运动启动。
Fig 1.每个实验的时间安排
每个试验的的时间安排表，显示的四种状态包括休眠(0-4秒)、AO(4-8秒)、idle(8-9秒)和任务执行(MI或ME(9-13秒)。
数据采集
搭建传感系统，在整个实验过程中同时记录EEG、EOG、EMG信号，如图2所示。
Fig 2
下图为国际10-20系统通道配置(11个EEG和2个EOG记录电极)。左面板上每个电极的对应位置；右边的面板表示索引。
Fig 3
EEG和EOG信号
使用g.USBampRESEARCH来重新记录EEG和EOG信号，如上图所示。
采样率设置为1200 Hz。
EEG：将11个电极放置在FCz，C3，Cz，C4，CP3，CPz，CP4，P3，Pz，P4和POz上
EOG：将2个电极放在右眼下方(VEOG)和(HEOG)上
在整个实验过程中，EEG和EOG信号的阻抗均保持在10kΩ以下
在MI和ME会话中均提供了EEG和EOG信号
MI / ME会话在每次从坐-站/从站-坐的转换时，原始EEG和EOG的尺寸为participants×runs×trials×channels×timepoints(8×3×5×6×16800)。
Fig 4.6个EMG记录电极的通道配置，显示每个电极的索引对应位置
EMG信号
使用OpenBCI记录EMG信号。
采样率设置为250 Hz。
将6个电极置于两下肢的股直肌(RF)、胫骨前肌(TA)和腓肠肌外侧(GL)上，如上图所示。
在实验会话上只有记录肌电图
每次坐-站/站-坐转换的原始肌电图形成在
participants×runs×trials×channels×timepoints (8×3×5×6×3500)。
数据预处理
离线信号处理使用MNE-python软件包(0.20.0版)执行。在MI和ME期间，预处理过程分为两个主要步骤：基于EEG的MI和基于EEG的MRCP。下图说明了EEG，EOG和EMG数据处理的过程。
Fig 5.EEG、EOG和肌电图数据预处理概述。
EEG、EOG和肌电图数据预处理概述:
展示了MI在EEG和EOG数据上的预处理过程;
展示了从ME中提取MRCPs的预处理步骤。
运动想象：陷波滤波器设置为50Hz，以减少电噪声。使用2阶非因果Butterworth滤波器对记录的EEG和EOG信号进行1–40 Hz之间的带通滤波。两个信号都被下采样到250 Hz。
将基于独立成分分析(ICA)的眼动伪影校正算法应用于EOG信号，对脑电数据中剔除的伪影分量进行识别。将脑电图信号按每类(R、AO、MI)开始的epoch分段4 s，如下图所示，然后采用移0.2 s的2 s滑动窗口进行预处理。
来自每个参与者的每个类别的处理数据包含试验×窗口×通道×时间点的集合(15×11×11×500)。
使用scikit-learn对15个试验使用left-one(trial) out交叉验证(LOOCV)实现受试者独立分类任务，如图6所示。在训练过程中，首先对训练集进行信号预处理，如图5所示。利用滤波器组公共空间模式(FBCSP)从下采样训练集中提取空间特征，生成用于分类任务的特征向量。重要的是，FBCSP在MI分类任务中通常表现良好。引入FBCSP作为公共空间模式(CSP)的扩展，从多个滤波器组中自主选择识别的脑电图特征。在本研究中，9个滤波器组带通滤波器的带宽为4 Hz从4到40 Hz (4 - 8 Hz, 8-12Hz，…，36–40 Hz)。构建了6个滤波器组带通滤波器，其中第一个频带的带宽为0.4 Hz，其他频带带宽为0.5 Hz (0.1-0.5 Hz,0.5 - 1 Hz，…，2.5-3Hz)。
Fig 6.
图6所示为基于网格搜索算法的二分类模型的leave-one(trial)-outcross-validation (LOOCV)架构。受试者逐个独立执行LOOCV。
Fig 7.
图7所示。在MI会话中使用的伪在线分析流程图。当连续检测5次(AO >= 5)产生动作观察时，应用网格搜索算法辅助确定动作观察与运动想象(AOvs.MI)分类模型。
图8.对MI任务的神经反应。在整个试验中，针对事件的频谱相关摄动(ERSP)在4–40 Hz之间进行分组，在(a)sit-to-stand和(b) stand-to-sit任务相比，R基线状态(1 - 0 s)。从0-4 s的时间间隔对应于R状态，从4-8 s的间隔对应于AO状态，从8-9 s的间隔对应于空闲状态，从9 s开始的时间间隔对应于执行状态。为了可视化，采样率设置为600 Hz。与基线相比，所有目前的ERSP值均具有统计学意义(p = 0:05)。
Fig 8. 对MI任务的神经反应
图9所示。从-1.5-1 s到实际运动开始，8个参与者sit-to-stand (a) andstand-to-sit (b)转换的Grandaverage MRCP波形(11个通道)。
Fig 9.
下图头皮地形图显示了MRCP振幅随时间变化的空间表征。
研究人员在这项研究中开发的任务中，参与者被指示对坐立和站坐的动作执行AO和MI/ME。实验结果表明，在AO期间ERD比较显著，而在MI期间ERS在感觉运动区域的alpha带较为典型。结合常用空间模式(FBCSP)和支持向量机(SVM)进行离线和分类器测试分析。离线分析表明，AO和MI的分类在站-坐转换时的平均准确率最高，为82.73±2.54%。通过分类器测试分析，研究人员证明了MI范式比ME范式具有更高的解码神经意图的性能。
这些观察让我们看到了使用基于AO和MI集成的开发任务构建未来基于外骨骼的康复系统的前景。
项目案例代码
fbcsp处理运动想象数据源码
run_fbcsp_mi
import numpy as npimport osimport sysimport csvfrom sklearn.model_selection import KFoldfrom pysitstand.model import fbcspfrom pysitstand.utils import sliding_window, sliding_window2from pysitstand.eeg_preprocessing import apply_eeg_preprocessing"""Binary classification model.We apply FBCSP-SVM (9 subbands from 4-40 Hz) on the subject-dependent scheme (leave a single trial for testing) for EEG-based MI classification.x sec window size with y% step (0.1 means overlap 90%)# How to run>> python run_fbcsp_mi.py >> python run_fbcsp_mi.py 2 0.1 4 stand R_vs_AO rASR>> python run_fbcsp_mi.py 2 0.1 4 sit R_vs_AO rASR>> python run_fbcsp_mi.py 2 0.1 4 stand AO_vs_MI rASR>> python run_fbcsp_mi.py 2 0.1 4 sit AO_vs_MI rASR>> python run_fbcsp_mi.py 2 0.1 4 stand AO_vs_MI ICA && python run_fbcsp_mi.py 2 0.1 4 stand R_vs_AO ICA>> python run_fbcsp_mi.py 2 0.1 4 sit AO_vs_MI ICA && python run_fbcsp_mi.py 2 0.1 4 sit R_vs_AO ICA"""def load_data(subject, task, prediction_model, artifact_remover, filter_order, window_size, step, sfreq): #load data the preprocessing # filter params notch = {'f0': 50} bandpass = {'lowcut': 1, 'highcut': 40, 'order': filter_order} ica = {'new_sfreq': sfreq, 'save_name': None, 'threshold': 2} rASR = {'new_sfreq': sfreq} # it will perform preprocessing from this order if artifact_remover == 'ICA': filter_medthod = {'notch_filter': notch, 'butter_bandpass_filter': bandpass, 'ica': ica} elif artifact_remover == 'rASR': filter_medthod = {'notch_filter': notch, 'butter_bandpass_filter': bandpass, 'rASR': rASR} # apply filter and ICA data = apply_eeg_preprocessing(subject_name=subject, session='mi', task=task, filter_medthod=filter_medthod) # data : 15 sec # define data R_class = data[:,:,int(2*sfreq):int(6*sfreq)] # rest 2 to 6 s AO_class = data[:,:,int(6*sfreq):int(10*sfreq)] MI_class = data[:,:,int(11*sfreq):int(15*sfreq)] # beep at 11s len_data_point = R_class.shape[-1] num_windows = int(((len_data_point-win_len_point)/(win_len_point*step))+1) # define class if prediction_model == 'R_vs_AO': # sliding windows R_class_slided = np.zeros([15, num_windows, 11, window_size*sfreq]) AO_class_slided = np.zeros([15, num_windows, 11, window_size*sfreq]) for i, (R,AO) in enumerate(zip(R_class, AO_class)): R_class_slided[i,:,:,:] = np.copy(sliding_window2(np.array([R]), win_sec_len=window_size, step=step, sfreq=sfreq)) AO_class_slided[i,:,:,:] = np.copy(sliding_window2(np.array([AO]), win_sec_len=window_size, step=step, sfreq=sfreq)) X0 = np.copy(R_class_slided) X1 = np.copy(AO_class_slided) elif prediction_model == 'AO_vs_MI': # sliding windows AO_class_slided = np.zeros([15, num_windows, 11, window_size*sfreq]) MI_class_slided = np.zeros([15, num_windows, 11, window_size*sfreq]) for i, (AO,MI) in enumerate(zip(AO_class, MI_class)): AO_class_slided[i,:,:,:] = np.copy(sliding_window2(np.array([AO]), win_sec_len=window_size, step=step, sfreq=sfreq)) MI_class_slided[i,:,:,:] = np.copy(sliding_window2(np.array([MI]), win_sec_len=window_size, step=step, sfreq=sfreq)) X0 = np.copy(AO_class_slided) X1 = np.copy(MI_class_slided) del data, R_class, AO_class, MI_class y0 = np.zeros([X0.shape[0], X0.shape[1]]) y1 = np.ones([X1.shape[0], X1.shape[1]]) assert len(X0) == len(y0) assert len(X1) == len(y1) return X0, y0, X1, y1if __name__ == "__main__": window_size = int(sys.argv[1]) # 1,2,3 sec. step = float(sys.argv[2]) # 0.1 --> overlap(90%) filter_order = int(sys.argv[3]) # 2 order of all fillter task = sys.argv[4] # stand, sit prediction_model = sys.argv[5] # R_vs_AO, AO_vs_MI artifact_remover = sys.argv[6] # ICA, rASR sfreq = 250 # new sampling rate [max = 1200 Hz] win_len_point = int(window_size*sfreq) for x in sys.argv: print("Argument: ", x) subjects = [ 'S01', 'S02', 'S03', 'S04', 'S05', 'S06', 'S07', 'S08'] if task == 'stand': save_name = 'sit_to_stand_mi' elif task == 'sit': save_name = 'stand_to_sit_mi' if prediction_model == 'R_vs_AO': save_path = 'MI-v2-'+artifact_remover+'-FBCSP-cv-'+str(window_size)+'s_'+task+'_'+prediction_model+'_filter_order_'+str(filter_order) elif prediction_model == 'AO_vs_MI': save_path = 'MI-v2-'+artifact_remover+'-FBCSP-cv-'+str(window_size)+'s_'+task+'_'+prediction_model+'_filter_order_'+str(filter_order) header = [ 'fold', 'accuracy', '0.0 f1-score', '1.0 f1-score', 'average f1-score', '0.0 recall', '1.0 recall', 'average recall', '0.0 precision', '1.0 precision', 'average precision', 'sensitivity', 'specificity' ] sum_value_all_subjects = [] for subject in subjects: from joblib import dump, load print('===================='+subject+'==========================') for directory in [save_path, save_path+'/model', save_path+'/y_slice_wise']: if not os.path.exists(directory): os.makedirs(directory) #load data the preprocessing X0, y0, X1, y1 = load_data(subject=subject, task=task, prediction_model=prediction_model, artifact_remover=artifact_remover, filter_order=filter_order, window_size=window_size, step=step, sfreq=sfreq) with open(save_path+'/'+save_path+'_result.csv', 'a') as csvFile: writer = csv.writer(csvFile) writer.writerow([str(subject)]) writer.writerow(header) kf = KFold(n_splits=15, shuffle=False) # Define the split - into 15 folds print(kf) accuracy_sum, precision_0_sum, recall_0_sum, f1_0_sum, precision_1_sum, recall_1_sum, f1_1_sum, precision_sum, recall_sum, f1_sum = [], [], [], [], [], [], [], [], [], [] sen_sum, spec_sum = [], [] predict_result = [] X_csp_com = [] for index_fold, (train_idx, test_idx) in enumerate(kf.split(X0)): print("=============fold {:02d}==============".format(index_fold)) print('fold: {}, train_index: {}, test_index: {}'.format(index_fold, train_idx, test_idx)) X0_train, X1_train = X0[train_idx], X1[train_idx] y0_train, y1_train = y0[train_idx], y1[train_idx] X0_test, X1_test = X0[test_idx], X1[test_idx] y0_test, y1_test = y0[test_idx], y1[test_idx] X_train = np.concatenate((X0_train.reshape(-1, X0_train.shape[-2], X0_train.shape[-1]), X1[train_idx].reshape(-1, X1_train.shape[-2], X1_train.shape[-1])), axis=0) y_train = np.concatenate((y0_train.reshape(-1), y1_train.reshape(-1)), axis=0) X_test = np.concatenate((X0_test.reshape(-1, X0_test.shape[-2], X0_test.shape[-1]), X1[test_idx].reshape(-1, X1_test.shape[-2], X1_test.shape[-1])), axis=0) y_test = np.concatenate((y0_test.reshape(-1), y1_test.reshape(-1)), axis=0) print("Dimesion of training set is: {} and label is: {}".format (X_train.shape, y_train.shape)) print("Dimesion of testing set is: {} and label is: {}".format( X_test.shape, y_test.shape)) # classification accuracy, report, sen, spec, X_test_csp, y_true, y_pred, classifier = fbcsp(X_train=X_train, y_train=y_train, X_test=X_test, y_test=y_test, filter_order=filter_order, session='mi') dump(classifier, save_path+'/model/'+subject+save_name+'_'+str(index_fold+1).zfill(2)+'.gz') # saving precision_0 = report['0.0']['precision'] recall_0 = report['0.0']['recall'] f1_0 = report['0.0']['f1-score'] precision_1 = report['1.0']['precision'] recall_1 = report['1.0']['recall'] f1_1 = report['1.0']['f1-score'] precision = report['weighted avg']['precision'] recall = report['weighted avg']['recall'] f1 = report['weighted avg']['f1-score'] accuracy_sum.append(accuracy) precision_0_sum.append(precision_0) recall_0_sum.append(recall_0) f1_0_sum.append(f1_0) precision_1_sum.append(precision_1) recall_1_sum.append(recall_1) f1_1_sum.append(f1_1) precision_sum.append(precision) recall_sum.append(recall) f1_sum.append(f1) sen_sum.append(sen) spec_sum.append(spec) row = [index_fold+1, accuracy, f1_0, f1_1, f1, recall_0, recall_1, recall, precision_0, precision_1, precision, sen, spec] predict_result.append([y_true, y_pred]) X_csp_com.append(X_test_csp) with open(save_path+'/'+save_path+'_result.csv', 'a') as csvFile: writer = csv.writer(csvFile) writer.writerow(row) print(subject, 'save DONE!!!!') print('***************************************') print('***************************************') print('***************************************') print('***************************************') mean_value = [np.mean(accuracy_sum), np.mean(f1_0_sum), np.mean(f1_1_sum), np.mean(f1_sum), np.mean(recall_0_sum), np.mean(recall_1_sum), np.mean(recall_sum), np.mean(precision_0_sum), np.mean(precision_1_sum), np.mean(precision_sum), np.mean(sen_sum), np.mean(spec_sum)] sum_value_all_subjects.append(mean_value) np.savez(save_path+'/y_slice_wise/'+subject+save_name+'.npz', x = np.array(X_csp_com), y = np.array(predict_result)) with open(save_path+'/'+save_path+'_result.csv', 'a') as csvFile: writer = csv.writer(csvFile) writer.writerow(['mean', mean_value[0], mean_value[1], mean_value[2], mean_value[3], mean_value[4], mean_value[5], mean_value[6], mean_value[7], mean_value[8], mean_value[9], mean_value[10], mean_value[11]]) writer.writerow([]) mean_all = np.mean(sum_value_all_subjects, axis=0) print(mean_all) with open(save_path+'/'+save_path+'_result.csv', 'a') as csvFile: writer = csv.writer(csvFile) writer.writerow(['accuracy', '0.0 f1-score', '1.0 f1-score', 'average f1-score', '0.0 recall', '1.0 recall', 'average recall', '0.0 precision', '1.0 precision', 'average precision', 'sensitivity', 'specificity' ]) writer.writerows(sum_value_all_subjects) writer.writerow(mean_all)后台发送"MIME"即可获取相关代码资源
论文信息
Decoding EEG Rhythms During Action Observation, Motor Imagery, and Execution for Standing and Sitting
文章仅限学习使用，不用于商业行为，若有侵权及疑问，请后台留言！
文章转载请联系后台管理人员！
更多阅读
脑机接口让脊髓损伤患者重新获得手部触觉
一种基于EEG和sEMG的假手控制策略
信号处理之频谱原理与python实现
EMD算法之Hilbert-Huang Transform原理详解和案例分析
EEGLAB 使用教程 3 -参考电极和重采样
第2期 | 国内脑机接口领域专家教授汇总
精彩长文 | 脑机接口技术的现状与未来！
ICA处理脑电资料汇总
收藏 | 脑电EEG基础与处理汇总
脑机接口BCI学习交流QQ群：903290195
微信群请扫码添加，Rose拉你进群
(请务必填写备注，eg. 姓名+单位+专业/领域/行业)
长按加群
欢迎点个在看鼓励一下
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
