简介：这是一份面向EEGLAB用户的脑电频谱参数化工具插件，基于FOOOF算法对神经功率谱进行自动拟合与绘图，适用于需要从脑电数据中分离周期性振荡成分与非周期性背景成分的研究人员。资源包共27个文件，包含18个MATLAB函数或脚本、6张结果示例图、1份Markdown说明文档以及若干Git版本控制文件，压缩包仅287KB，轻量且易于部署。目前已有678人学习下载，属于脑电频谱分析方向较受关注的小型工具资源。插件提供菜单级调用与命令行接口，支持批量执行频谱拟合、参数提取和结果出图；说明文档中列出了依赖安装、环境版本校验以及插件挂载步骤，能帮助研究者在自有数据集上复现参数化流程，减少从原始脑电数据到可发表图谱之间的中间环节。 做脑电的朋友应该都有这种体验：数据预处理熬了大半天，终于画出被试的平均功率谱，老板过来看了一眼，问了一句——“你这个alpha峰的变化，到底是振荡活动本身增强了，还是背景活动变了？”
这个问题，用传统的频带功率分析很难回答。因为经典做法是把某个频段（比如8~12 Hz）的能量直接拎出来算均值，可这个值里其实混着两类完全不同的生理成分：周期性振荡峰，以及宽频的1/f背景活动。俩东西搅在一起，结果既不好解释，也容易被低频漂移带偏。FOOOF（Fitting Oscillations & One Over F）就是为拆开这两部分而生的。它把功率谱参数化成“非周期背景 + 若干个高斯峰”，每个峰对应一个振荡成分，背景对应1/f趋势，这样你就能理直气壮地回答老板的问题。
我这篇就基于一个真实可复现的工作流：用MATLAB做主力环境，EEGLAB完成脑电预处理，MATLAB版的FOOOF对功率谱做参数化拟合，最后用MATLAB把拟合结果画成能直接放进论文的图。整条链路不需要切到Python，数据从原始信号到最终成图，一个软件全包。适合正在做EEG/ERP研究、手头有数据但不知道频谱参数化怎么落地的朋友参考。
1. 项目整体思路与工作流设计
1.1 为什么选MATLAB + EEGLAB + FOOOF这套组合
先说结论：这套组合的每一环，都是当前生态里最稳的选择之一。
EEGLAB在EEG预处理上的积累太深了，从数据导入、电极定位、滤波分段、坏段剔除，到ICA去伪迹，都有成熟的图形界面和命令行函数。对新手来说，GUI点一点能上手；对老手来说，写成脚本跑批量，也不拖后腿。更重要的是EEGLAB处理完的数据结构（EEG struct）可以直接拿来算功率谱，不需要自己再折腾数据格式转换。
FOOOF解决的问题，前面已经提了：把功率谱拆成周期性成分和非周期性成分。传统频谱分析做完，你只能画一条PSD曲线，说“哦，alpha频段功率很高”；FOOOF做完了，你能得到三个可量化的参数——alpha峰的峰值高度（peak height）、中心频率（center frequency）、带宽（bandwidth），以及背景的offset和exponent。这几个参数可以直接进统计模型，做组间比较、条件对比，甚至与行为数据做相关。
MATLAB作为中间桥梁，好处是“一栈式”。EEGLAB是MATLAB工具箱，FOOOF有官方MATLAB移植版，数据流动不用经过CSV或文件转换，内存里直接传递。这一点在批量处理几十个被试时非常省心。
1.2 从原始数据到最终成图的完整流程
我建议把整套流程拆成五个阶段，每个阶段有明确的输入输出，方便排查问题。
原始EEG数据（.set / .raw / .edf） ↓ 阶段一：EEGLAB预处理（滤波、分段、ICA等） 干净的连续/分段EEG数据 ↓ 阶段二：功率谱密度计算（每个电极分别算） freqs（频率向量） + psd（功率谱） ↓ 阶段三：FOOOF参数化拟合 aperiodic_params、peak_params、拟合谱 ↓ 阶段四：可视化绘图 单电极/群体/地形图等多类型图片 ↓ 阶段五：统计分析与结果解释我实际跑的时候，阶段二和阶段三之间最容易出问题。原因在于EEGLAB的spectopo函数默认输出的是dB单位，而FOOOF内部会再做一次log变换，单位搞错，拟合出来的结果就是废的。后面我会专门讲这个坑。
2. 预处理关键环节：EEGLAB中的操作与踩坑
2.1 导入数据与电极定位
这个步骤看起来基础，但真不能跳。FOOOF拟合的是每个电极上的功率谱，如果电极位置信息缺失，后面画地形图时topoplot直接罢工。
%% 1. 导入数据（以.set为例） EEG = pop_loadset('filepath', 'D:\eeg_data\', 'filename', 'sub01.set'); %% 2. 电极定位（如果数据里没有chanlocs） EEG = pop_chanedit(EEG, 'lookup', 'standard-10-5-cap385.elp');这里standard-10-5-cap385.elp是EEGLAB自带的国际10-5系统电极坐标文件，路径一般在EEGLAB安装目录下的plugins/dipfit/里。如果电极帽是64导的，通常pop_chanedit会自动匹配标准坐标；如果用的是自定义电极帽，就要自己检查坐标文件对不对。
我的建议是：建模时用标准坐标，但特定情况要手动核对电极名。比如有些厂商把Cz命名为Cz1，EEGLAB可能不会自动对齐，最后拟合出来的功率谱错位到隔壁电极上，这种情况我用一次性脚本核对通道顺序，确保channel名称与EEG.data的行号对应。
2.2 滤波、分段与剔除伪迹的关键参数
预处理参数这事，不同实验室习惯差异大。我不打算替所有人下结论，只说我实测下来比较稳的组合，以及参数背后的理由。
%% 1. 带通滤波：0.1~45 Hz EEG = pop_eegfiltnew(EEG, 0.1, 45); %% 2. 分段：以事件为锚点，取[-1, 3]秒 EEG = pop_epoch(EEG, {'Stimulus'}, [-1 3]); %% 3. 基线校正 EEG = pop_rmbase(EEG, [-1000 0]); %% 4. 坏段剔除（联合概率法，阈值±5个标准差） EEG = pop_jointprob(EEG, 1, [1 1], [5 5]); %% 5. ICA去伪迹 EEG = pop_runica(EEG, 'icatype', 'runica'); EEG = pop_iclabel(EEG, 'default'); % 然后根据ICLabel的分类结果，人工确认后剔除眨眼和肌肉成分滤波带宽这块要特别小心。FOOOF拟合的是1/f背景，如果高通截止频率设置得太高（比如1 Hz以上），会把低频段真实的1/f斜率削掉一部分，导致拟合出的exponent虚高。我通常设0.1 Hz，虽然会引入一些缓慢漂移，但ICA和基线校正能处理一部分，更重要的是保留了低频段的形状信息。
分段长度也影响PSD频率分辨率。频率分辨率=1/窗口长度，如果分段是4秒，分辨率是0.25 Hz，FOOOF对窄峰（比如~10 Hz的alpha）的估计精度还可以；如果分段只有1秒，分辨率是1 Hz，两个相邻的峰（比如theta 6 Hz和alpha 9 Hz）就容易被糊成一个宽峰，拟合时peak数量会被低估。我的原则是：做FOOOF的epoch至少2秒，如果实验设计允许，4秒更好。
2.3 计算功率谱密度的正确姿势
这是整个工作流里最容易翻车的地方，我单独拿出来讲。
EEGLAB里计算每个电极平均功率谱最方便的函数是spectopo：
%% 计算每个电极的平均功率谱 % 输入是EEG.data（通道×时间×试次），返回spectra（通道×频率）和freqs（频率向量） [spectra, freqs] = spectopo(EEG.data, 0, EEG.srate, ... 'chanlocs', EEG.chanlocs, 'freqrange', [1 45], 'plot', 'off');注意，spectopo返回的spectra单位是dB（即10×log10(power)），但FOOOF的MATLAB版在内部会再次对输入PSD做log变换。所以如果你直接把dB值丢给FOOOF，就相当于做了两次log，拟合结果会完全错乱。
正确的做法是：把dB先转回线性功率，再交给FOOOF。
%% 转回线性单位（μV²/Hz） psd_lin = 10.^(spectra / 10);或者，如果你的FOOOF版本支持log_psd参数，也可以直接传入dB值并设置'log_psd', true。但我更推荐前者，因为直接在外部处理单位，逻辑透明，换版本也不容易踩雷。
另外，spectopo默认用的是Welch平均法，通过多段重叠窗平均得到平滑的PSD，窗口长度默认为EEG.srate的整数倍。在分段数据上，spectopo内部会用每个试次分别计算再平均，这个平滑效果对FOOOF非常友好。千万别图省事直接把数据拼成一长串再调pwelch，试次之间的边缘跳变会引入高频能量，把beta/gamma频段的峰搞出很多假阳性。
3. FOOOF拟合：把频谱拆成可解释的成分
3.1 FOOOF究竟在做什么
FOOOF的核心思想不复杂，它假设观测到的功率谱由两部分组成：非周期背景（apériodic background）加上若干个周期性的高斯峰（oscillatory peaks）。
数学形式大致是：
PSD(f) = b + log(1/f^exponent) + sum(高斯峰)拟合时，FOOOF先估计非周期背景的参数（offset和exponent），然后从残差里逐次迭代找出符合阈值的高斯峰，记录每个峰的中心频率、高度和带宽。这一步看着简单，但真正做起来，参数设置对结果影响非常大。
3.2 MATLAB调用FOOOF的完整代码
MATLAB版FOOOF官方仓库是TDonoghue/fooof_matlab，clone下来后把仓库路径加入MATLAB搜索路径即可。
%% 添加FOOOF路径 addpath(genpath('D:\toolbox\fooof_matlab')); %% 设置拟合参数 settings = struct(); settings.peak_width_limits = [1 8]; % 峰的半高全宽范围（Hz） settings.max_n_peaks = 6; % 最多检测6个峰 settings.min_peak_height = 0.1; % 峰高阈值，过低容易把噪声当峰 settings.peak_threshold = 2.0; % 峰显著性阈值，相对于残差标准差 settings.aperiodic_mode = 'fixed'; % 固定模式：背景为1/f %% 拟合范围 fit_range = [1 40]; % 避开50Hz工频 %% 循环处理每个电极 channel_results = struct(); for ch = 1:size(psd_lin, 1) freqs_tmp = freqs(freqs >= fit_range(1) & freqs <= fit_range(2)); psd_tmp = psd_lin(ch, freqs >= fit_range(1) & freqs <= fit_range(2)); current_result = fooof(freqs_tmp, psd_tmp, fit_range, settings); channel_results(ch).aperiodic_params = current_result.aperiodic_params; channel_results(ch).peak_params = current_result.peak_params; channel_results(ch).fooofed_spectrum = current_result.fooofed_spectrum; channel_results(ch).r_squared = current_result.r_squared; end几点说明：
peak_width_limits设[1 8]，是经验值。EEG的振荡峰半高宽一般在2~6 Hz之间，太窄的峰（比如0.5 Hz）基本是伪迹；太宽的峰（比如20 Hz）可能把两个邻近峰糊在一起。
max_n_peaks我设6，是因为在1~40 Hz范围内，典型能看到的节律就是theta（4~7）、alpha（8~12）、beta（13~30）这几个谱段，超过6个峰基本是过拟合。
aperiodic_mode选'fixed'还是'knee'取决于数据。如果功率谱在低频段有明显的“拐点”（即低频不完全是线性下降，而是变平），可以用'knee'，它会多拟合一个膝盖参数；但如果所有被试的曲线形状差异不大，用'fixed'更稳定，组间比较时参数也更少、更好解释。我一般先画总平均功率谱，看到明显拐点才用knee。
3.3 参数设置的实验依据
FOOOF结果好不好，八成取决于三个设置：拟合频率范围、peak高度阈值、peak宽度限制。
拟合频率范围我坚持1~40 Hz，避开45 Hz以上有两个考虑：一是50Hz工频即使被陷波滤波拉低了，残余能量仍可能在拟合时被识别成高峰；二是多数认知实验关心的节律都在40 Hz以下，范围设窄一点能减少计算量。
peak高度阈值min_peak_height需要根据数据量级调整。如果用的是μV²/Hz线性单位，典型alpha峰高度大约在0.5~3之间，阈值0.1能过滤掉大部分随机波动；但如果你的数据是微幅级别（比如小动物EEG），功率整体较小，阈值可能要降到0.02~0.05。一个实用技巧：先不设阈值跑一遍，画出原始谱和拟合谱的对比图，看看那些“你认为应该被识别”的峰有没有被忽略。如果没有，说明阈值合适；如果alpha峰都没被识别，就把阈值调低一个数量级。
我踩过的坑是，把peak_threshold设得过大，比如3.0，结果alpha峰虽然肉眼可见，但因为残差标准差较大，FOOOF认为它不显著，直接不报这个峰。后来我改用peak_threshold=2.0，同时配合min_peak_height，效果稳定很多。
4. 绘制拟合结果图：从单条曲线到论文级成图
4.1 单电极、单被试的频谱拟合对比图
这是最基础的一类图，也是最直观的：一条黑色线是原始功率谱，一条红色虚线是FOOOF拟合谱，两条线重合度高，说明拟合效果好；差距大，说明参数没调对。
%% 以Cz电极为例 ch_fit = strcmp({EEG.chanlocs.labels}, 'Cz'); figure('Color', 'w', 'Position', [100 100 800 450]); plot(freqs_tmp, 10*log10(psd_tmp), 'k', 'LineWidth', 1.2); hold on; plot(freqs_tmp, 10*log10(channel_results(ch_fit).fooofed_spectrum), ... 'r--', 'LineWidth', 1.8); xlabel('Frequency (Hz)', 'FontSize', 12); ylabel('Power (dB)', 'FontSize', 12); legend({'Original PSD', 'FOOOF fit'}, 'Location', 'northeast', 'FontSize', 11); set(gca, 'FontSize', 11, 'box', 'off'); xlim(fit_range);注意这里绘图时统一用dB，因为线性功率的数值范围太大，画出来低频段会把高频段压扁，什么都看不清。dB转换后的图，视觉上才和EEGLAB里常见的功率谱图一致。
这个图建议作为第一个检查项——每当要批量处理被试之前，先随机挑一个被试、两三个电极跑一遍，把拟合谱和原始谱叠在一起看，确认没问题再上批量。批量跑完以后，最好输出每个电极的拟合R²，低于0.9的电极单独标记出来重新检查。
4.2 群体水平的参数可视化
FOOOF的核心产出是参数。最常用的对比指标是每个振荡峰的高度和中心频率。比如想比较两组被试的alpha峰高度差异，可以画出柱状图叠加误差棒和散点：
%% 假设g1和g2是两个组的alpha peak height g1 = [1.8 2.1 1.5 2.3 1.9 2.0]; g2 = [1.2 1.5 1.0 1.8 1.4 1.6]; g1_mean = mean(g1); g1_sem = std(g1) / sqrt(length(g1)); g2_mean = mean(g2); g2_sem = std(g2) / sqrt(length(g2)); figure('Color', 'w', 'Position', [100 100 600 450]); bar([g1_mean, g2_mean], 0.6, 'FaceColor', [0.7 0.7 0.7]); hold on; errorbar([1 2], [g1_mean, g2_mean], [g1_sem, g2_sem], ... 'k', 'LineStyle', 'none', 'LineWidth', 1.5); % 叠加散点，让读者看到个体分布 scatter(ones(size(g1)) + randn(size(g1))*0.05, g1, 40, 'k', 'filled', 'jitter', 'on', 'jitterAmount', 0.05); scatter(2*ones(size(g2)) + randn(size(g2))*0.05, g2, 40, 'k', 'filled', 'jitter', 'on', 'jitterAmount', 0.05); set(gca, 'XTick', [1 2], 'XTickLabel', {'Group1', 'Group2'}, 'FontSize', 12, 'box', 'off'); ylabel('Alpha Peak Height (μV²/Hz)', 'FontSize', 12);这里用散点叠加是近年神经影像论文的主流画法，比单纯柱状图+误差棒传递的信息多得多。如果样本量比较大（比如n>30），散点可以用半透明的“小提琴图”或者“雨云图”，但MATLAB内置函数没有现成的雨云图，我用的是第三方函数raincloud_plot，效果也不错。
4.3 地形图呈现FOOOF参数的空间分布
参数算出来后，除了看单个电极，还应该看全脑分布。用EEGLAB的topoplot可以把每个电极的alpha峰值高度画成地形图，一眼看出是顶枕部高、额部低，还是全脑一致。
%% 假设已经算出所有电极的alpha峰值高度：alpha_peak_heights（通道×1） figure('Color', 'w', 'Position', [100 100 500 420]); topoplot(alpha_peak_heights, EEG.chanlocs, ... 'maplimits', 'maxmin', ... 'electrodes', 'on', ... 'shading', 'interp'); colorbar; title('Alpha Peak Height Topography', 'FontSize', 12);这个图的细节在于maplimits。如果两个条件或两组被试要放在同一scale下比较，不要用maxmin，而是取所有数据的最小值和最大值作为统一的maplimits，否则每一张图各自归一化，视觉上会掩盖真实的组间差异。
我还习惯把地形图叠加在单个被试的拟合对比图旁边，这样单被试报告里既有局部频谱细节，又有全脑分布概览。论文里做“figure 1”或者“supplementary”都很合适。
5. 常见问题与排查技巧实录
5.1 FOOOF拟合效果差、R²偏低怎么办
先别急着调FOOOF参数，先怀疑PSD本身。我遇到过最多次的情况是：spectopo输入的freqrange上限超过了Nyquist频率，或者数据里有大量未剔除的噪声段，导致某个电极的PSD完全不像脑电。
排查顺序： 1. 看原始PSD曲线有没有异常尖峰 → 有则查工频干扰、肌电伪迹 2. 看R²分布 → 多个电极R²<0.8，说明预处理阶段有问题 3. 看每个电极检测到的peak数量 → 全部都是0，说明阈值太高 4. 看拟合谱是否平滑跟随原始谱 → 不跟，说明PSD输入单位错了单位错误是最隐蔽的。我建议在FOOOF之前，先用一个简单测试：把psd_lin取log，画出来，看低频段是否是明显的向下倾斜直线。如果不是，说明你的PSD可能经过了一次log变换，需要回退。
5.2 批量处理几十个被试，速度太慢怎么优化
FOOOF的拟合过程是逐电极、逐被试循环的，每个电极要迭代搜索高斯峰，电极大（比如128导）的时候确实慢。我实测下来，64导、1~40 Hz、6个peak上限，单个电极大约0.3~0.8秒，一个被试一两分钟，20个被试一小时左右。如果觉得慢，两个办法：
第一个，用parfor并行处理电极。前提是每个电极的拟合相互独立，这正好满足。把for ch替换成parfor ch，再开个并行池：
parpool('local', 4); % 按CPU核心数调整注意，parfor里不要频繁写入同一个struct数组，最好是把结果存成cell，循环结束后再合并。否则传输开销会吃掉并行收益。
第二个办法，降低峰值搜索上限。如果研究只关心alpha峰，max_n_peaks设成3就够，别让FOOOF把精力浪费在找不存在的beta/gamma峰上。
5.3 不同被试的peak param对不齐，没法做统计怎么办
这是FOOOF入坑之后最常见的苦恼：被试A的alpha峰中心频率是10.2 Hz，被试B是8.5 Hz，两组数据直接做t检验，看起来像是同一个峰，但实际生理意义可能不完全一致。
我的处理方法是“band-constrained peak extraction”。在FOOOF输出peak_params之后，按照中心频率所在频带进行归类：
%% 把peak按频带归类 alpha_peaks = peak_params(peak_params(:,1) >= 8 & peak_params(:,1) <= 13, :); theta_peaks = peak_params(peak_params(:,1) >= 4 & peak_params(:,1) <= 7, :);如果某个被试的某个频带里检出了多个峰，取最高的那个；如果一个都没检出，记为NaN。这样后续统计时，组间比较的是“该频带最显著峰的高度”，而不是笼统的“所有峰的高度”，结果解释起来更干净。
这个方法唯一的风险是，如果某组被试的alpha中心频率普遍超出了8~13 Hz范围，那就会被误记为缺失。所以正式分析前，一定要画一张所有被试peak中心频率的分布直方图，确认频带边界设置合理再往下走。
5.4 一个容易被忽略的可重复性问题
FOOOF的输出依赖初始化状态，同一份数据跑两次，结果理论上是一致的，但我在实践里发现，当某个峰的显著性接近阈值时，拟合结果可能出现细微的差别。这通常发生在peak_threshold在临界值附近的时候。
为了可重复性，我在批量处理前固定随机种子：
rng(42);然后在保存结果时，把FOOOF的版本号、MATLAB版本、settings结构体、拟合范围全部存进一个JSON或MAT文件里，和输出的CSV放一起。这样审稿人问起来，我能精确说明每一行数据是怎么得到的。
最后再分享一个小技巧：FOOOF拟合出来的背景参数（exponent）本身就是一个很有意思的指标，它反映的是棘波活动、兴奋/抑制平衡等生理信息。不要只盯着峰看，把1/f斜率也纳入分析，很多时候组间差异恰恰藏在这个背景趋势里。我第一次跑完整套流程时，发现两组的peak height没有显著差异，但exponent差异很显著，这个发现直接改变了文章的分析方向。
本文还有配套的精品资源，点击获取
