基于分类任务的信号(EEG)处理--代码分步解析

最新推荐文章于 2024-03-07 14:50:00 发布

脑机接口社区

最新推荐文章于 2024-03-07 14:50:00 发布

阅读量3.3k

点赞数 8

分类专栏：脑机接口社区 # 脑电信号处理文章标签：脑机接口脑电数据 BCI EEG 脑电波

原文链接：https://blog.csdn.net/Kay_Xiaohe_He/

版权

脑机接口社区同时被 2 个专栏收录

402 篇文章 623 订阅

订阅专栏

脑电信号处理

143 篇文章 409 订阅

订阅专栏

本文由网友Jon_Snow_Stark授权分享

本篇文章是对作者对另一篇文章《基于分类任务的信号(EEG)处理》的扩展，之前文章是从宏观方面介绍利用EEG信号进行分类任务，本文利用详细的代码进行解析，大家可以两篇结合着来看，希望可以帮助到大家。

读取脑电信号

在读取设备采集的脑电信号上EEGLAB是一个非常强大的工具包，我在本文中就是使用这一工具包。首先在MATLAB的命令行输入eeglab（前提是你已经在MATLAB中添加了EEGLAB工具包），则会弹出EEGLAB的GUI界面，大家可以通过GUI界面上的按钮和调用相关函数进行操作，调用函数大家可以通过help按钮进行查阅。下面我们通过按钮读取相关数据。

依次点击File -> Import data -> Using EEGLAB funtions and pludins -> From biosemi BDF file(最后一个是选择符合你保存的数据的格式的读入方式)。

接下来会弹出两个框，直接选择OK即可。导入完成后导入的数据的信息会显示在GUI界面上，而导入的数据则会保存在工作区的EEG结构中。我们也可以打开EEG查看我们导入的数据，脑电数据就保存在data中，后续对脑电信号的处理就是对EEG.data进行处理。至此，我们读取数据的过程就全部完成了，就得到了可以用于计算的数据了。

预处理（分频带提取特征，功率谱特征PSD）

首先将得到的脑电信号拿出来，成为一个矩阵：

data = [];
data = [data;EEG.data'];
data = double(data);

然后获得脑电数据矩阵的通道数和样本数，从上边图片中EEG.data变量可以看到是按照一个通道一行进行排列的，但是在取出EEG.data时我进行了转置（该步可以不转，后续处理按行向量处理即可），那么我们读到的矩阵大小行数即为采样点数，列数即为通道数：

[samples,channals] = size(data);

然后初始化采样率（从GUI界面可以看到采样率为1000Hz），由于我的标签是每30s一个，因此我设定时间窗口为30s：

frequency = 1000;
time_win = 30;

通过读取标签文件获得标签的数量，定义为label_num。
接下来就是按通道分别对每个通道进行分频带和提取特征，用一个for循环进行：

for i = 1:channals
  .......
end

for循环里边就是预处理过程。首先对信号进行分频带，此处使用的是巴特沃斯带通滤波器：

%提取四个频带
delta = butter_bandpass_filter(data(:,i), 1, 4, frequency,3);
theta = butter_bandpass_filter(data(:,i), 4, 8, frequency,3);
alpha = butter_bandpass_filter(data(:,i), 8, 14, frequency,3);
beta  = butter_bandpass_filter(data(:,i), 14, 30, frequency,3);

此时将信号分别滤波到四个频带上，然后分别提取四个频带的功率谱密度（PSD）特征。首先我们要明白我提取特征要得到一个什么样的结果：我们提取特征是要用这个特征进行分类，那么提取之后就是一个分类标签对应一个特征，每个分类标签都对应自己的一个特征，然后分类器学习相同特征之间的相似性，区分不同特征之间的不同，得到一个分类模型，后续你再向这个模型输入一个分类未知的数据的时候它就可以根据学习到的内容对这个数据进行预测了。那么我的标签有label_num个，每一个标签对应30s的数据，那么我就要以30s为切片，得到一个PSD特征，最终每个频带得到label_num个特征。那么我就循环label_num次，30s的动态窗口依次滑动label_num节，每一次取30s*frequency个点（即30s的数据），计算得到一个特征。

%提取四个频带的PSD特征
psd_delta = [];
psd_theta = [];
psd_alpha = [];
psd_beta  = [];
for j = 1:label_num
    psd_delta = [psd_delta;compute_psd(delta(((j-1)*frequency*time_win+1):j*frequency*time_win, 1))]; %提取delta频带PSD
    psd_theta = [psd_theta;compute_psd(theta(((j-1)*frequency*time_win+1):j*frequency*time_win, 1))]; %提取theta频带PSD
    psd_alpha = [psd_alpha;compute_psd(alpha(((j-1)*frequency*time_win+1):j*frequency*time_win, 1))]; %提取alpha频带PSD
    psd_beta  = [psd_beta;compute_psd(beta(((j-1)*frequency*time_win+1):j*frequency*time_win, 1))];   %提取beta 频带PSD
end

经过计算，分别得到了四个频带label_num个特征，此时还是四个分开的特征向量，我们按照delta、theta、alpha、beta的顺序将该通道的数据组合成一个label_num*4的矩阵。

psd_temp = [];
psd_temp = [psd_temp; psd_delta,psd_theta,psd_alpha,psd_beta];

那么每个通道的四个频带的特征就保存在psd_temp这个矩阵中了，矩阵的大小为label_num*4，各列分别对应四个频带，每一行对应对应标签的四个频带的特征。这样计算我们得到的是各个通道分开的特征矩阵，可是我们想把每个人的所有通道的所有特征都保存在一个矩阵里，这该怎么操作呢？

首先我们在按通道进行for循环的前边创建一个空矩阵psd_decomposed = [];然后将psd_temp并入到该空矩阵中。这样循环完后所有特征就都保存在一个矩阵里了，最后要记得检查以下保存的对不对，reshape就行了，reshape到正确的格式，如果不正确程序会终止。

psd_decomposed = [];
for i = 1:channals
  ......
  psd_decomposed = [psd_decomposed,psd_temp];
end
psd_decomposed = reshape(psd_decomposed,label_num,channals*4);

再之后就是把预处理得到的特征矩阵进行保存啦（及时保存数据，避免重复计算，耗费大量之间），

PSD = [];
PSD = [PSD,psd_decomposed];
save('F:\EEG_Thermal_Sleep\Subject3_2\EEGsleep6.mat','PSD');

分类器分类

我使用的LIBSVM分类器，相关信息大家可以参考网上其他人写的。在利用SVM进行分类时一定要对数据进行归一化，不然你可能会发现无论怎么搞都得不到理想分类精度，作为一个曾经怀疑人生的小白的忠告！

% 对数据进行归一化，有minmax和zscore两种
normalized_data = mapminmax(signal);
% normalized_data1 = zscore(signal);

简单的可以对数据进行4/6分或者2/8分进行训练和测试看一下分类精度，其实最好的还是进行10折交叉验证。
这里给出我利用KFold进行交叉验证的代码，关于KFold的用法和原理大家可以自行在网上检索。

% 利用SVM建模,用KFOLD方式进行10折交叉验证
K = 10; % 折叠次数，即要分包的个数
indices = crossvalind('Kfold',label,K); % 将数据分为5个包，并为每个包分配自己对应的索引
acc =[];
predicted = [];
for k = 1:K
    test = (indices == k); % 第k个包作为测试集
    train = ~test;         % 除了第k个包作为训练集
    trainlabel = label(train,1);          % 从标签中获得训练标签
    testlabel = label(test,1);            % 从标签中获得测试标签
    traindata = normalized_data(train,:); % 从数据集中获得训练数据集
    testdata  = normalized_data(test,:);  % 从数据集中获得测试数据集
    model = svmtrain(trainlabel, traindata,'-b 1');
    [predicted_label, acc_one, decision_values] = svmpredict(testlabel, testdata, model,'-b 1');
    acc = [acc,acc_one];
    predicted = [predicted;predicted_label];
end
acc_mean = mean(acc')';
accuracy(:,n) = acc_mean;

参考：

https://blog.csdn.net/Kay_Xiaohe_He/