脑机接口实例二:脑电信号CSP处理

文章目录二、实例代码及结果分析 三、问题与反思四、下一阶段计划
前言
本文在实例一的基础之上,选用BCI大赛08年数据集,基于CSP的特征提取,运用算法进行集成分类,以期实现运动想象中左、右手动作脑电信号的分类 。
一、实例说明 1.1 期望目标1.2 数据解读 实验范式
数据集由9名受试者的脑电图数据组成 。BCI范式由四个不同的运动想象任务组成,即左手(类别1)、右手(类别2)、双脚(类别3)和舌头(类别4)的运动想象 。范式记录不同日期的两天数据 。(每天6圈,每圈) 。(共计576个,每个任务)
正式范式之前,记录5分钟估计EOG值,分三块:看屏幕上的十字(2分钟)、闭眼(1分钟)、眼动(1分钟),注:由于技术问题被试A04T的眼动数据较短,仅包含眼动状态 。如下图所示:
被试者坐在电脑屏幕前面的舒适的扶手椅上,在试验开始时(t=0s),黑屏上出现了一个十字 。此外,还提出了一种短声告警音 。两秒后(t=2s),以箭头的形式指向左、右、下或上(对应于四个类中左手、右手、脚或舌头的一个)的提示,出现在屏幕上1.25秒 。这促使受试者完成所需的运动想象任务 。不提供任何反馈 。被试在t=6s时被要求执行运动想象任务,直到屏幕上的十字消失为止 。随后屏幕再次为黑色的短暂中断,流程如下图所示:(有效数据为2-6s,即运动想象任务记录区间,亦为需提取数据的区间)
数据记录

脑机接口实例二:脑电信号CSP处理

文章插图
250Hz采样,并在0.5Hz和100Hz之间进行带通滤波 。除了22个脑电图通道,3个单极EOG通道,以250Hz的频率记录和采样(见下图) 。它们在0.5Hz和100Hz之间进行带通滤波,放大器的灵敏度设置为1mV 。EOG通道是为伪像处理方法的后续应用提供的,不得用于分类 。
数据文件描述
事件类型描述如下 (event type 769/770 对应左右手的动作,以此为基准进行数据提取) :
二、实例代码及结果分析 2.1 代码
主函数(main):clc;clear;%% 读取数据并对数据预处理Subjects = 9;%被试数Fs = 250;%采样率windowLength = 4;%单次采样时间chanSelect = [8,10,12]; %通道选择c3,c4,cztotalFlt = [4 40]; %总的滤波频段选择load rawdata1.mat[Data_train,label_train] = preProccess(Fs,windowLength,EEG,chanSelect,totalFlt); %preProccess执行数据提取、分段、滤波处理load rawdata2.mat[Data_test,label_test] = preProccess(Fs,windowLength,EEG,chanSelect,totalFlt);%% CSP特征提取EEGSignals.x = Data_train;EEGSignals.y = label_train;Y = label_train;classLabels = unique(EEGSignals.y); CSPMatrix = learnCSP(EEGSignals,classLabels);nbFilterPairs = 1;X = extractCSP(EEGSignals, CSPMatrix, nbFilterPairs); EEGSignals_test.x=Data_test;EEGSignals_test.y=label_test;T = extractCSP(EEGSignals_test, CSPMatrix, nbFilterPairs); save dataCSP.mat X Y Tcolor_L = [0 102 255] ./ 255;color_R = [255, 0, 102] ./ 255;pos = find(Y==1);plot(X(pos,1),X(pos,2),'x','Color',color_L,'LineWidth',2);hold onpos = find(Y==2);plot(X(pos,1),X(pos,2),'o','Color',color_R,'LineWidth',2);legend('Left Hand','Right Hand')xlabel('C3','fontweight','bold')ylabel('C4','fontweight','bold')
CSP算法:function CSPMatrix = learnCSP(EEGSignals,classLabels)%%Input:%EEGSignals: the training EEG signals, composed of 2 classes. These signals%are a structure such that:%EEGSignals.x: the EEG signals as a [Ns * Nc * Nt] Matrix where%Ns: number of EEG samples per trial%Nc: number of channels (EEG electrodes)%nT: number of trials%EEGSignals.y: a [1 * Nt] vector containing the class labels for each trial%EEGSignals.s: the sampling frequency (in Hz)%%Output:%CSPMatrix: the learnt CSP filters (a [Nc*Nc] matrix with the filters as rows)%%See also: extractCSPFeatures%check and initializationsnbChannels = size(EEGSignals.x,2);nbTrials = size(EEGSignals.x,3);nbClasses = length(classLabels);if nbClasses ~= 2disp('ERROR! CSP can only be used for two classes');return;endcovMatrices = cell(nbClasses,1); %the covariance matrices for each class%% Computing the normalized covariance matrices for each trialtrialCov = zeros(nbChannels,nbChannels,nbTrials);for t=1:nbTrialsE = EEGSignals.x(:,:,t)';%note the transposeEE = E * E';trialCov(:,:,t) = EE ./ trace(EE);endclear E;clear EE;%computing the covariance matrix for each classfor c=1:nbClassescovMatrices{c} = mean(trialCov(:,:,EEGSignals.y == classLabels(c)),3); %EEGSignals.y==classLabels(c) returns the indeces corresponding to the class labelsend%the total covariance matrixcovTotal = covMatrices{1} + covMatrices{2};%whitening transform of total covariance matrix[Ut Dt] = eig(covTotal); %caution: the eigenvalues are initially in increasing ordereigenvalues = diag(Dt);[eigenvalues egIndex] = sort(eigenvalues, 'descend');Ut = Ut(:,egIndex);P = diag(sqrt(1./eigenvalues)) * Ut';%transforming covariance matrix of first class using PtransformedCov1 =P * covMatrices{1} * P';%EVD of the transformed covariance matrix[U1 D1] = eig(transformedCov1);eigenvalues = diag(D1);[eigenvalues egIndex] = sort(eigenvalues, 'descend');U1 = U1(:, egIndex);CSPMatrix = U1' * P;