往期文章

希望了解更多的道友点这里

0. 分享【脑机接口 + 人工智能】的学习之路

1.1 . 脑电EEG代码开源分享 【1.前置准备-静息态篇】

1.2 . 脑电EEG代码开源分享 【1.前置准备-任务态篇】

2.1 . 脑电EEG代码开源分享 【2.预处理-静息态篇】

2.2 . 脑电EEG代码开源分享 【2.预处理-任务态篇】

3.1 . 脑电EEG代码开源分享 【3.可视化分析-静息态篇】

3.2 . 脑电EEG代码开源分享 【3.可视化分析-任务态篇】

4.1 . 脑电EEG代码开源分享 【4.特征提取-时域篇】

4.2 . 脑电EEG代码开源分享 【4.特征提取-频域篇】

4.3 . 脑电EEG代码开源分享 【4.特征提取-时频域篇】

4.4 . 脑电EEG代码开源分享 【4.特征提取-空域篇】

5 . 脑电EEG代码开源分享 【5.特征选择】

6.1 . 脑电EEG代码开源分享 【6.分类模型-机器学习篇】

6.2 . 脑电EEG代码开源分享 【6.分类模型-深度学习篇】

汇总. 专栏：脑电EEG代码开源分享【文档+代码+经验】

0 . 【深度学习】常用网络总结

一          、前言

本文档旨在归纳BCI-EEG-matlab的数据处理代码          ，作为EEG数据处理的总结 
  
 
  
 
  
 
  
 
  
，方便快速搭建处理框架的Baseline









，实现自动化 
  
 
  
 
  
 
  
 
  
、模块插拔化









、快速化

 
  
  
  
  
 。本文以任务态（锁时刺激 
  
  
  
  
 ，如快速序列视觉呈现）为例 
  
  
  
  
 ，分享脑电EEG的分析处理方法
  
  
  
  
  。

脑电数据分析系列
 


 


 


 


 

。分为以下6个模块： 前置准备 数据预处理 数据可视化 特征提取（特征候选集） 特征选择（量化特征择优） 分类模型

本文内容：【4. 特征提取-频域篇】

提示：以下为各功能代码详细介绍

 


 


 


 


 
，若节约阅读时间 
 
 
 
 
，请下滑至文末的整合代码

二 
  
  
  
  
 、特征提取 框架介绍

特征提取作为承上启下的重要阶段 
   
   
   
   
  ，是本系列中篇幅最长的部分 
 
 
 
 
。承上
 

 

 

 

 
，紧接预处理结果和可视化分析
 

 

 

 

 
，对庞大的原始数据进行凝练          ，用少量维度指标表征整体数据特点；启下 

 
 

 
 

 
 

 
 

 
，这些代表性 
  
  
  
  
 、凝练性的特征指标量化了数据性能









，为后续的认知解码

 


 


 


 


 
、状态监测 
 
 
 
 
、神经调控等现实需求提供参考
   
   
   
   
   。

特征提取的常用特征域为时域 
   
   
   
   
  、频域
 

 

 

 

 
、时频域
 

 

 

 

 
、空域等
 

 

 

 

 
。本文特征主要为手动设置的经验特征               ，大多源于脑科学及认知科学的机制结论 
 
 
 
 
 
 
 
 
 
，提取的特征具有可解释的解剖          、认知 

 
 

 
 

 
 

 
 

 
、物理含义；也有部分是工程人员的实践发现














，在模型性能提升中效果显著 

 

 

 

 

。

特征提取的代码框图









、流程如下所示：

时频域-特征提取的主要功能          ，分为以下2部分：

小波变换 本证模态分解

小波变换：小波变换是时频变换最经典的方法之一 
  
 
  
 
  
 
  
 
  
，通过滑动的小波基拆解信号时频信息          。时频变换最常用的可视化方式为 雨谱图 or 瀑布图（不同领域名称不同）









，但总体思路就是将信号展开为 时域x频域

的 二维图像 
  
  
  
  
 ，可以清晰观察到随时间变化相应的频域变化

 
 

 
 

 
 

 
 

 
 。

我在这里画了4种视觉姓名刺激诱发的Cz电极的脑响应雨谱图 
  
  
  
  
 ，可以看出视觉刺激诱发了在500ms的3-5Hz响应

 


 


 


 


 
，脑激活强度与姓名内容有关 
 
 
 
 
，强度大小为自身>熟人>陌生人>空白











。

本文未应用小波变换代码 
   
   
   
   
  ，主要因为小波变换对小波基选择十分敏感
 

 

 

 

 
，对脑电信号和小波基选择间缺少固定经验
 

 

 

 

 
，因此鼓励大家自己尝试          ，探索适合自己脑电任务的小波基               。 本证模态分解：(Empirical Mode Decomposition,EMD) 本征模态分解也称经验模态分解 

 
 

 
 

 
 

 
 

 
，目前广泛应用于脑电处理领域
 
 
 
 
 
 
 
 
 
 。不同于时频分解方法固定的小波基               、短时傅里叶变换等









，本征模态分解优势在于其 自适应的分解基底               ，可以让一个领域内毫无经验的小白快速上手














。matlab 2017之后的版本已自带官方

EMD算法函数          。

EMD分解大致思路为 
 
 
 
 
 
 
 
 
 
，1.分别连接信号的极大值点 
 
 
 
 
 
 
 
 
 
、极小值点














，形成上包络














、下包络
  
 
  
 
  
 
  
 
  
 。2.三次样条差值分别拟合上包络和下包络曲线          ，这一步主要为了光滑









。3.上包络和下包络曲线均值 
  
  
  
  
 。4.原始信号减去包络均值 
  
 
  
 
  
 
  
 
  
，结果为1阶的本征模态函数

（IMF）









，完成1次分解
  
  
  
  
  。重复1-4对信号继续分解
 


 


 


 


 

。分解公式如下图：

我画了各阶IMF的频域图如下 
  
  
  
  
 ，可以看出EMD一遍遍滤掉了信号中的相对高频 
  
  
  
  
 ，IMF逐渐向低频聚拢 
 
 
 
 
。我第一次看结果的时候鸡皮疙瘩都起来了

 


 


 


 


 
，极值包络想法实在太巧妙的过滤掉了信号中的相对高频


   
   
   
   
   。

按个人经验及文献阅读结论 
 
 
 
 
，脑电信号一般选择2阶IMF进行下一步的特征处理 
   
   
   
   
  ，即下图红色虚线对应的结果
 

 

 

 

 
。一方面
 

 

 

 

 
，2阶IMF频率基本集中在0-60Hz
 

 

 

 

 
，是正常脑电能量集中的频带；另一方面          ，一阶IMF初步滤掉了高频噪声



 
 

 
 

 
 

 
 

 
，数据质量有改善 

 

 

 

 

。

下图 EMD基础理论引用自： EMD基础理论









，也可看出各阶IMF时域内振荡周期（高频成分）逐渐降低          。

个人对于时频域的薄见供大家参考：时域和频域都内含着丰富的潜在特征               ，并且两个特征域是互联互通          、互相转换的

 
 

 
 

 
 

 
 

 
 。但纠其一面则失全貌 
 
 
 
 
 
 
 
 
 
，就像好看的姑娘正面和侧脸都有不一样的美感









。特征提取需要我们将原数据掰开了 
  
 
  
 
  
 
  
 
  
、嚼碎了把最简单凝练的特征设计出来














，因此时域多一点









、还是频域多一点          ，就成了时频域常面临的平衡问题               。这种你中有我 
  
  
  
  
 、我中有你的 时 
  
  
  
  
 、频共生关系很有趣 
  
 
  
 
  
 
  
 
  
，就像我的老师第一次总结小波变换









，她说时域和频域必然损失一个 
  
  
  
  
 ，学会取舍
 
 
 
 
 
 
 
 
 
 。

三

 


 


 


 


 
、代码格式说明

本文非锁时任务态（下文以静息态代替）范例为：ADHD患者 
 
 
 
 
、正常人群在静息状态下的脑模式分类

代码名称：代码命名为Festure_ candidate_xxx (time / freq/ imf/ space) 参数设置：预处理结果\采样率\时域是否非线性熵特征（耗时）\频域均分分辨度\imf阶数\space对比通道数及频带范围














。 输入格式：输入格式承接规范预处理最后一项输出 
  
  
  
  
 ，Proprocess_xxx(预处理最终步骤)_target/nontarget          。 输出及保存格式：输出格式为试次数*特征个数

 


 


 


 


 
，按照除空域特征外 
 
 
 
 
，按照通道的特征拼接 
   
   
   
   
  ，先为1通道内的所有特征
 

 

 

 

 
，接着2通道的所有特征
  
 
  
 
  
 
  
 
  
 。保存文件名称为Festure_candidate_xxx(特征域名称)_target/nontarget









。

三 
   
   
   
   
  、脑电特征提取 代码

提示：代码环境为 matlab 2018

3.0 参数设置

可视化内容可以选择
 

 

 

 

 
，把希望可视化特征域写在Featute_content 中

一次进行10人次的批处理          ，subject_num = [1;10] 特征提取内容： Featute_content = [‘time’,‘freq’,‘time_freq’,‘space’]; 时域
 

 

 

 

 
、频域
 

 

 

 

 
、时频域          、空域均分析 时频特征选择 本征模态算法 + 微分熵 ： Featute_time_freq_content = [‘emd’,‘DE’]; 2阶的本征模态函数提取特征：imf_level = 2; %% 0.特征候选集-参数设置 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data_path = C:\Users\EEG\Desktop\basetest_flod\save_fold\; svae_path = C:\Users\EEG\Desktop\basetest_flod\save_fold\; subject_num = [1;10]; freq_resolut = 2; freq_scale = [1;60]; Featute_content = [time\,freq\,time_freq\,space]; Featute_time_freq_content = [emd\,DE\]; imf_level = 2; disp([||特征候选集-参数设置||]); disp([特征域内容: , Featute_content]); disp([时域-候选集: , Featute_time_content]);

3.1 标准输入赋值

导入上一步预处理阶段处理后的数据：

%% 1.标准输入赋值 Proprocess_target_file = load([data_path ,Proprocess_target_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))]); Proprocess_nontarget_file = load([data_path ,Proprocess_nontarget_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))]); stuct_target_name = Proprocess_target; stuct_nontarget_name = Proprocess_nontarget; Proprocess_target_data = Proprocess_target_file.(stuct_target_name).data; Proprocess_nontarget_data = Proprocess_nontarget_file.(stuct_nontarget_name).data; subject_num = Proprocess_target_file.(stuct_target_name).subject_num; fs_down = Proprocess_target_file.(stuct_target_name).fs_down; remain_trial_target = Proprocess_target_file.(stuct_target_name).remain_trial; remain_trial_nontarget = Proprocess_nontarget_file.(stuct_nontarget_name).remain_trial; disp([目标试次剩余: , num2str(remain_trial_target),||平均： , num2str(mean(remain_trial_target))]); disp([非目标试次剩余: , num2str(remain_trial_nontarget),||平均： , num2str(mean(remain_trial_nontarget))]);

3.2 时频域-特征提取

主函数中 调用时频域提取函数

主函数 调用 时频域 特征提取函数Festure_candidate_time_freq

%% 4.时频域特征候选集 if contains(Featute_content,time_freq) disp([时频域特征计算中...]); tic; [Festure_time_freq_target,Festure_time_freq_candidate_num_target]= Festure_candidate_time_freq(Proprocess_target_data,Featute_time_freq_content,remain_trial_target,fs_down,imf_level); [Festure_time_freq_nontarget,Festure_time_freq_candidate_num_nontarget]= Festure_candidate_time_freq(Proprocess_nontarget_data,Featute_time_freq_content,remain_trial_nontarget,fs_down,imf_level); if contains(Featute_time_freq_content,DE) Festure_time_freq_target = log(Festure_time_freq_target); Festure_time_freq_nontarget = log(Festure_time_freq_nontarget); end t_time_freq_candidate_cost = toc; disp([时频域特征计算完毕 

 
 

 
 

 
 

 
 

 
，耗时(秒)： ,num2str(t_time_freq_candidate_cost)]); Festure_candidate_time_freq_target = []; Festure_candidate_time_freq_target.data = Festure_time_freq_target; Festure_candidate_time_freq_target.Featute_time_freq_content = Featute_time_freq_content; Festure_candidate_time_freq_target.remain_trial_target = remain_trial_target; Festure_candidate_time_freq_target.Festure_time_freq_candidate_num_target = Festure_time_freq_candidate_num_target; Festure_candidate_time_freq_target.fs_down = fs_down; Festure_candidate_time_freq_nontarget = []; Festure_candidate_time_freq_nontarget.data = Festure_time_freq_nontarget; Festure_candidate_time_freq_nontarget.Featute_time_freq_content = Featute_time_freq_content; Festure_candidate_time_freq_nontarget.remain_trial_nontarget = remain_trial_nontarget; Festure_candidate_time_freq_nontarget.Festure_time_freq_candidate_num_nontarget = Festure_time_freq_candidate_num_nontarget; Festure_candidate_time_freq_nontarget.fs_down = fs_down; disp([时频域特征保存中...]); save([ svae_path , Festure_candidate_time_freq_target_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))],Festure_candidate_time_freq_target); save([ svae_path , Festure_candidate_time_freq_nontarget_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))],Festure_candidate_time_freq_nontarget); disp([时频域特征保存完毕]); end

3.2.1 时频域特征提取函数

时频域 特征提取函数Festure_candidate_time_freq









，调用的EMD函数为matlab2018自带

function [Festure_time_freq,Festure_time_freq_candidate_num]= Festure_candidate_time_freq(Standard_input_data,Featute_time_freq_content,remain_trial,fs_down,imf_level) Festure_time_freq = []; %% 1.emd只进行传统5频带 five_band if contains(Featute_time_freq_content,emd) fest_five_band = []; count_trial = 1; for sub_loop = 1:size(remain_trial,2) for trial_loop = 1:remain_trial(1,sub_loop) five_band_temp = []; for channel_loop = 1:size(Standard_input_data{1,1},1) fft_temp = []; imf_temp = []; imf_temp = emd(Standard_input_data{trial_loop, sub_loop}(channel_loop,:),Interpolation,pchip,MaxNumIMF,imf_level ,Display,0); fft_temp = abs(fft(imf_temp(:,imf_level),fs_down)); five_band_temp(channel_loop,:) =sum_five_band(fft_temp,fs_down); end fest_five_band(count_trial,:) = reshape(five_band_temp,1,size(five_band_temp,1)*size(five_band_temp,2)); count_trial = count_trial+1; end end end %% 2.汇总视频特征 Festure_time_freq = [ fest_five_band]; Festure_time_freq_candidate_num = size(Festure_time_freq,2); end

3.2.3 传统5频带 方法

function five_band_sum = sum_five_band(fft_temp,fs_down) %% 这只是一行的5频带求和               ，请在外面加Channe_loop循环 delta =[1;4]; %δ theta =[4;8]; %θ alpha =[8;12]; %α? beta = [12;30]; %β ? gamma =[30;60]; %γ ? five_band = [delta theta alpha beta gamma]; fft_resolut = fs_down/size(fft_temp,2); epoch_num = size(five_band,2); epoch_length = five_band.*fft_resolut; five_band_sum = []; for cut_loop = 1:epoch_num fft_sum_temp = sum(fft_temp(:, epoch_length(1,cut_loop): epoch_length(2,cut_loop))); five_band_sum = [five_band_sum fft_sum_temp]; end end

总结

时频域特征融合了各自的长处 
 
 
 
 
 
 
 
 
 
，交叉

了时域频域的信息














，

方便研究人员更全面的了解信号特点 
  
  
  
  
 。

目前时频特征还是在长时任务中应用较多          ，归因于时频分解还是注重频带

的信息 
  
 
  
 
  
 
  
 
  
，

长时任务有较宽的频带能量分布









，而任务态脑电的频域集中在低频

本文着重介绍的EMD算法 
  
  
  
  
 ，突出了自适应

的基底优势 
  
  
  
  
 ，

建议新入门的朋友可以尝试使用
  
  
  
  
  。

也推荐大家尝试更多的分解模式如：变分模态分解

 


 


 


 


 
，动态模态分解

 
 
 
 
 
，

据其他领域的师兄称效果也出奇的好
 


 


 


 


 

。

同时 
   
   
   
   
  ，对经典特征的融合 

 
 

 
 

 
 

 
 

 
、组合也是发掘更优

混合特征的常用方式 
 
 
 
 
。

大家可以探索和发掘是用自己研究的优质特征策略
   
   
   
   
   。

目前多样性的特征还在不断发展









、丰富
 

 

 

 

 
，新的特征提取方法逐渐多元化
 

 

 

 

 
。

进阶特征如脑网络               、拓扑图等
 

 

 

 

 
，基于人工智能的端到端特征提取方法          ，会在新的专栏中介绍 

 

 

 

 

。

囿于能力 

 
 

 
 

 
 

 
 

 
，挂一漏万









，如有笔误请大家指正~

感谢您耐心的观看               ，本系列更新了约30000字 
 
 
 
 
 
 
 
 
 
，约3000行开源代码














，体量相当于一篇硕士工作          。

往期内容放在了文章开头          ，麻烦帮忙点点赞 
  
 
  
 
  
 
  
 
  
，分享给有需要的朋友~

坚定初心









，本博客永远：

免费拿走 
  
  
  
  
 ，全部开源 
  
  
  
  
 ，全部无偿分享~

To：新想法 
 
 
 
 
 
 
 
 
 
、鬼点子的道友：

自己：脑机接口+人工智领域

 


 


 


 


 
，主攻大脑模式解码














、身份认证          、仿脑模型…

在读博士第3年 
 
 
 
 
，在最后1年 
   
   
   
   
  ，希望将代码 
  
 
  
 
  
 
  
 
  
、文档









、经验 
  
  
  
  
 、掉坑的经历分享给大家~

做的不好请大佬们多批评 
  
  
  
  
 、多指导~ 虚心向大伙请教！

想一起做些事情 or 奇奇怪怪点子 or 单纯批评我的
 

 

 

 

 
，请至Rongkaizhang_bci@163.com