往期文章

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

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