脑电波特征提取（脑电EEG代码开源分享 【4.特征提取-频域篇】）
往期文章
希望了解更多的道友点这里
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. 特征提取-频域篇】
提示：以下为各功能代码详细介绍 ，若节约阅读时间 ，请下滑至文末的整合代码
二 、特征提取 框架介绍
特征提取作为承上启下的重要阶段 ，是本系列中篇幅最长的部分 。承上 ，紧接预处理结果和可视化分析 ，对庞大的原始数据进行凝练 ，用少量维度指标表征整体数据特点；启下 ，这些代表性 、凝练性的特征指标量化了数据性能 ，为后续的认知解码 、状态监测 、神经调控等现实需求提供参考 。
特征提取的常用特征域为时域 、频域 、时频域 、空域等 。本文特征主要为手动设置的经验特征 ，大多源于脑科学及认知科学的机制结论 ，提取的特征具有可解释的解剖 、认知 、物理含义；也有部分是工程人员的实践发现 ，在模型性能提升中效果显著 。
特征提取的代码框图 、流程如下所示：
频域-特征提取的主要功能 ，分为以下3部分：
传统5频带分段法 均分频带法 微分熵传统5频带分段法：5个公认脑电频段范围如下表所示 。由于5频带分段法具有一定的神经科学理论依据 ，目前5频带分段的频域能量提取方法应用最广泛 。主要流程为分别计算各频段内的能量叠加 。需要注意的是 ，各频段的边界频率仍存在争论 ，下表中4\8\14\30Hz的边界是普遍认可的 。此外 ，建议脑电数据时间窗截取长一些 ，可以更准确刻画各频点的能量值 ，由于包括0.5Hz的频率分辨率 ，建议时间窗在2秒
以上 。
均分频带法：本人在实际操作中也常用5频带法 ，但是不免发现一些工程应用上的限制 。例如 ，当导联电极较少时 ，例如1导联 ，那么频域只能提取5个频域特征 ，难以满足分类器的输入需求 。此外 ，5频带的划分相对宽泛 、粗糙 ，大脑复杂丰富的功能仅划分5个频率区间 ，如果有宽度1Hz的细节特征会被淹没在宽频带中 。
因此 ，我们提出一种偏工程类的均分频带方法 ，将脑电能量集中在0-60Hz的频带均分为30段 ，每段频带宽度为2Hz 。均分方法可能缺少神经科学依据 ，但是在工程实践中大大提升了特征维度 ，用频域细节弥补空间信息 ，测试效果确实提升分类性能 。均分频带法还能定位到精准频点 ，例如通过特征选择方法筛选出16-17Hz特征效果较好 ，有助于反向推动该16-17Hz频段与实验任务的机制研究。
当然 ，均分频带法对于总体频带范围是0-60Hz还是10-80Hz ，以及均分为30段还是40频带 ，都是可调整的 ，可根据采集的频率分辨率 、导联数量 、任务主要频带进行设置 。
如下图蓝色虚线是5频带法 ，红色虚线是均分频带方法：
微分熵：微分熵并不像上文单指一种频带分段方法 ，而是对特征附加的非线性变换 ，尤其在频域能量特征中效果显著 。起初微分熵的应用在情绪分类等方向 ，目前应用范围逐渐扩大至运动想象 、静息态检测等领域 。计算的原理也较为简单 ，仅对能量特征求log 。个人经验建议大家尝试一下 ，有的任务中确实有效 ，但有的数据中和简单能量特征的性能接近 。
三 、代码格式说明
本文非锁时任务态（下文以静息态代替）范例为：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’]; 时域 、频域 、时频域 、空域均分析 时域特征内容：均分频带法 、5频带法 、微分熵 。Featute_freq_content = [‘average_band’,‘five_band’,‘DE’]; 均分频带的窗口长度2Hz：freq_resolut = 2; 均分频带的总体区间 1-60 Hz：freq_scale = [1;60]; %% 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_freq_content = [average_band\,five_band\,DE\]; 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_freq
%% 3.频域特征候选集 if contains(Featute_content,freq) disp([频域特征计算中...]); tic; [Festure_freq_target,Festure_freq_candidate_num_target]= Festure_candidate_freq(Proprocess_target_data,Featute_freq_content,remain_trial_target,freq_resolut,fs_down,freq_scale); [Festure_freq_nontarget,Festure_freq_candidate_num_nontarget]= Festure_candidate_freq(Proprocess_nontarget_data,Featute_freq_content,remain_trial_nontarget,freq_resolut,fs_down,freq_scale); if contains(Featute_freq_content,DE) Festure_freq_target = log(Festure_freq_target); Festure_freq_nontarget = log(Festure_freq_nontarget); end t_freq_candidate_cost = toc; disp([频域特征计算完毕 ，耗时(秒)： ,num2str(t_freq_candidate_cost)]); Festure_candidate_freq_target = []; Festure_candidate_freq_target.data = Festure_freq_target; Festure_candidate_freq_target.Featute_freq_content = Featute_freq_content; Festure_candidate_freq_target.remain_trial_target = remain_trial_target; Festure_candidate_freq_target.Festure_freq_candidate_num_target = Festure_freq_candidate_num_target; Festure_candidate_freq_target.fs_down = fs_down; Festure_candidate_freq_nontarget = []; Festure_candidate_freq_nontarget.data = Festure_freq_nontarget; Festure_candidate_freq_nontarget.Featute_freq_content = Featute_freq_content; Festure_candidate_freq_nontarget.remain_trial_nontarget = remain_trial_nontarget; Festure_candidate_freq_nontarget.Festure_freq_candidate_num_nontarget = Festure_freq_candidate_num_nontarget; Festure_candidate_freq_nontarget.fs_down = fs_down; disp([频域特征保存中...]); save([ svae_path , Festure_candidate_freq_target_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))],Festure_candidate_freq_target); save([ svae_path , Festure_candidate_freq_nontarget_,num2str(subject_num(1,1)),_,num2str(subject_num(2,1))],Festure_candidate_freq_nontarget); disp([频域特征保存完毕]); end3.2.1 频域特征提取函数
频域 特征提取函数Festure_candidate_freq
function [Festure_freq,Festure_freq_candidate_num]= Festure_candidate_freq(Standard_input_data,Featute_freq_content,remain_trial,freq_resolut,fs_down,freq_scale) Festure_freq = []; %% 1.均分频带average_band 或 传统5频带 five_band fest_average_band = []; fest_five_band = []; count_trial = 1; for sub_loop = 1:size(remain_trial,2) for trial_loop = 1:remain_trial(1,sub_loop) average_band_temp = []; five_band_temp = []; for channel_loop = 1:size(Standard_input_data{1,1},1) fft_temp = []; fft_temp = abs(fft(Standard_input_data{trial_loop,sub_loop}(channel_loop,:),fs_down)); if contains(Featute_freq_content,average_band) average_band_temp(channel_loop,:) =sum_average_band(fft_temp,fs_down,freq_resolut,freq_scale); end if contains(Featute_freq_content,five_band) five_band_temp(channel_loop,:) =sum_five_band(fft_temp,fs_down); end end if contains(Featute_freq_content,average_band) fest_average_band(count_trial,:) = reshape(average_band_temp,1,size(average_band_temp,1)*size(average_band_temp,2)); end if contains(Featute_freq_content,five_band) fest_five_band(count_trial,:) = reshape(five_band_temp,1,size(five_band_temp,1)*size(five_band_temp,2)); end count_trial = count_trial+1; end end %% 2. Festure_freq = [fest_average_band fest_five_band]; Festure_freq_candidate_num = size(Festure_freq,2); end3.2.2 均分频带 方法
function average_band_sum = sum_average_band(fft_temp,fs_down,freq_resolut,freq_scale) fft_resolut = fs_down/size(fft_temp,2); epoch_num = (freq_scale(2,1) - freq_scale(1,1) + 1)/freq_resolut; epoch_length = freq_resolut/fft_resolut; average_band_sum = []; for cut_loop = 1:epoch_num fft_sum_temp = sum(fft_temp(:, (cut_loop-1)*epoch_length +1: cut_loop*epoch_length)); average_band_sum = [average_band_sum fft_sum_temp]; end end3.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总结
大脑神经元放电产生的振荡 、节律信息 ，对节律的探索仍是起步和模糊
状态 ，
个人认为脑电EEG具有特征处理风格 ，锁时任务特征偏时域 ，长时任务特征偏频域，
本系列将脑电任务分为锁时 、长时的原因也在此脑电作为一种随机性+节律性
的神经信号 ，目前对节律性的频谱分析较多
但对随机性的探索还不足 ，推荐大家结合随机信号分析推出新的见解同时 ，对经典特征的融合 、组合也是发掘更优
混合特征的常用方式 。
大家可以探索和发掘是用自己研究的优质特征策略 。目前多样性的特征还在不断发展 、丰富 ，新的特征提取方法逐渐多元化 。
进阶特征如脑网络 、拓扑图等 ，基于人工智能的端到端特征提取方法 ，会在新的专栏中介绍 。囿于能力 ，挂一漏万 ，如有笔误请大家指正~
感谢您耐心的观看 ，本系列更新了约30000字 ，约3000行开源代码 ，体量相当于一篇硕士工作 。
往期内容放在了文章开头 ，麻烦帮忙点点赞 ，分享给有需要的朋友~
坚定初心 ，本博客永远：
免费拿走 ，全部开源 ，全部无偿分享~To：新想法 、鬼点子的道友：
自己：脑机接口+人工智领域 ，主攻大脑模式解码 、身份认证 、仿脑模型…
在读博士第3年 ，在最后1年 ，希望将代码 、文档 、经验 、掉坑的经历分享给大家~
做的不好请大佬们多批评 、多指导~ 虚心向大伙请教！
想一起做些事情 or 奇奇怪怪点子 or 单纯批评我的 ，请至Rongkaizhang_bci@163.com创心域SEO版权声明：以上内容作者已申请原创保护，未经允许不得转载，侵权必究！授权事宜、对本内容有异议或投诉，敬请联系网站管理员，我们将尽快回复您，谢谢合作！