脑电信号频谱分析实战：从功率谱密度到Welch参数调优

写这个系列的第一篇之前，我先说一个后台被问过很多次的问题：手头有一段脑电数据，到底应该先看时域波形，还是直接看频谱？我的答案一直很固定——时域波形只适合判断有没有坏段、有没有漂移，用它来判断“这个被试当下处于什么状态”“某个实验条件到底改变了什么”，效率非常低。脑电信号最核心的价值恰恰藏在频率成分里：alpha 波有没有变强、theta 能量是不是升高了、beta 活动在哪个电极最明显，这些都要靠谱分析来量化。这个系列里，我计划按“谱分析 → 时频联合分析 → 功能连接”的顺序往下写，所以第一期先把最基础的功率谱密度估计讲透。这篇文章不打算堆数学，而是从实操角度告诉你：谱分析能找到什么、找不到什么，用什么方法估计，参数到底该怎么定，以及哪些坑我替你先踩过了。

1. 动手之前先把问题对齐：谱分析到底向你报告什么

1.1 谱分析不是“把一个波形变成另一个图”那么简单

很多人第一次接触频谱，是从 FFT 开始的。把一段脑电信号做傅里叶变换，拿到幅值谱，画出来，然后对着某个峰值说“这是 alpha”——这个流程本身没错，但如果你只做到这一步，很容易在后续解释中犯错误。因为谱分析本质上是一个统计估计过程，它回答的问题是：在一段给定的时间里，信号能量在频率方向上是怎么分布的？

这句话里有两个关键词值得画重点。

第一个是“一段给定时间”。做普通谱分析时，我们把时间信息抹掉了。它假设被分析的这一段数据在统计意义上是平稳的，即频率成分不随时间剧烈变化。静息态脑电基本满足这个假设，但事件相关任务、睡眠纺锤波、癫痫发作期这类非平稳信号，直接用功率谱去概括整体状态，就会丢掉关键信息。这个问题我会在系列第二篇讲时频分析时重点展开。

第二个是“能量分布”。功率谱的核心指标不是“哪个频率出现了”，而是“哪个频率携带了多少能量”。它和幅值谱的区别在于：幅值谱看的是单个频率分量的幅度，功率谱看的是该分量在单位频率上的平均能量。脑电图分析里绝大多数场景都使用功率谱，因为脑电信号的能量更符合我们对节律强度的直觉——alpha 峰高，代表枕区 alpha 节律的功率更强。

1.2 谱分析能帮你回答哪几类问题

根据我自己的经验，脑电谱分析最常被用来解决四类问题：

• 静息态特征刻画：比如给被试做 5 分钟闭眼静息记录，比较不同人群的 alpha 峰值频率、alpha 功率、theta/beta 比值。
• 状态变化检测：比如比较睁眼和闭眼、清醒和困倦状态下 alpha 功率的差异；比较任务前后不同频带功率的变化。
• 事件相关同步/去同步（ERS/ERD）的预备研究：虽然真正的 ERS/ERD 需要在时间维度上滑动窗口，但你先用整段谱分析确认该频带在这个任务里确实有明显表现，再用时频方法去追踪动态变化，会更有把握。
• 数据质量筛查：查看全频段功率谱，可以快速发现 50 Hz 工频干扰、高频肌电污染、超低频漂移等问题。它比肉眼盯原始波形可靠得多。

所以，我建议你在读入数据后做的第一件正经事，不是画时域图，而是先画一张全通道的平均功率谱。这一步能让你对数据“健康程度”获得一个非常直观的印象，之后再做预处理、坏段剔除，你的决策都有依据。

2. 频谱估计的底层逻辑：从 DFT 到功率谱，为什么不能直接拿 FFT 结果去解释

2.1 FFT 只是数学工具，它给的结果并不能直接用

如果你在 MATLAB 或 Python 里对一段脑电波形直接调用 FFT，然后把每个点的幅值画出来，你大概率会得到一条毛刺很多、谱峰很宽、看起来并不干净的曲线。这不是你数据有问题，而是因为你忘了：原始 FFT 的谱线数量和信号长度绑定，频率轴的间隔等于 fs / N，其中 fs 是采样率，N 是参与 FFT 的点数。

假设采样率是 1000 Hz，你对 1 秒的数据做 FFT，频率分辨率只有 1 Hz；对 10 秒的数据做 FFT，分辨率提升到 0.1 Hz。但脑电信号里 alpha 峰和 theta 峰可能只差 3~4 Hz，1 Hz 分辨率虽然能把两者粗略分开，却很难准确估计峰值频率；如果两个节律靠得更近，比如 9.8 Hz 和 10.4 Hz 的子峰，1 Hz 分辨率根本分辨不出来。

这就是频谱估计的第一个核心概念：频率分辨率和数据长度成反比。想要 0.1 Hz 的分辨率，至少要取 10 秒连续数据。很多人调试半天参数，最后发现谱峰还是糊成一团，根源往往不是算法，而是给算法的数据长度不够。

2.2 从幅值谱到功率谱密度：单位是个大坑

做过一次 FFT 后，标准做法是对每个频率点的复数幅值取模平方，得到功率，再根据你的分析目标决定要不要除以频率单位。脑电数据分析中经常会遇到“功率谱密度（PSD）”和“功率谱”混用的局面。PSD 的单位是 μV²/Hz，反映的是单位频率宽度上的功率；而窄带功率的单位是 μV²，反映某个频带宽度内的总功率。

这两个单位在不同的软件包里默认输出不一样。EEGLAB 的频谱函数经常会给出 μV²/Hz，MNE-Python 的默认情况也需要留意。如果你投稿时只写“alpha 功率显著升高”，却不注明单位，那审稿人大概率会质疑你。我自己习惯统一用 μV²/Hz，并且在方法部分专门写清楚：计算时使用 Welch 法，Hanning 窗，窗长 4 秒，50% 重叠，频率分辨率约 0.5 Hz。

有一点必须强调：对单段数据做 FFT 得到的谱，是“周期图”（periodogram），它并不是一致估计。什么意思？就是随着数据长度越来越大，谱估计的方差不会明显下降。纯粹用周期图去看脑电，你会看到波动剧烈的谱线，真正的谱峰会淹没在噪声里。所以工程和科研中真正可靠的做法，要么是对多段数据取平均（Welch 法），要么是多锥度法，要么是其他平滑手段。这一点下一节我会展开讲。

3. 数据预处理和分段：谱的质量往往从这里决定

3.1 伪迹在频域里长什么样

我在处理学生数据时，最常看到的一种错误是：预处理没做干净，直接把坏数据丢进 FFT，然后得到一个看起来“很漂亮”但实际上充满伪迹的谱图。脑电里几类最常见伪迹在频域里有非常明显的痕迹：

• 工频干扰：50 Hz（国内）或 60 Hz（部分地区）处出现极窄的尖峰，而且还会在 100 Hz、150 Hz 等谐波位置出现小峰。
• 肌电伪迹：肌肉活动主要以 20~200 Hz 甚至更高频为主，表现为宽频带的高能量基底，额头、颞肌附近的电极尤其明显。
• 眼动和眨眼：主要能量集中在 4 Hz 以下，会抬高 delta 频段功率，严重时让低频段出现一个向下倾斜的“悬崖”。
• 电极滑动或接触不良：表现为全频段不规则的功率升高，通常伴随原始波形中的陡峭偏移或方波样变化。

所以在做谱分析前，我强烈建议你先按完整流程走一遍：定位坏导、剔除明显漂移和肌肉段、用 ICA 去掉眨眼和心电成分、必要时做 0.5~40 Hz 带通滤波。这里要特别提醒：滤波会改变谱的边缘区域，所以如果分析频带包含 theta（4~8 Hz），高通截止频率不要设在 1 Hz 以上；如果关注 gamma，带通滤波的低通截止频率至少要留出 10 Hz 以上的余量，同时要留意陷波滤波带来的频点凹陷。

3.2 参考方式对谱的影响比很多人以为的大

很多初学者对“参考”不敏感，总觉得影响不大。实际上，参考电极的选择会改变每个导联记录的瞬时电压，并间接影响功率谱的幅值和空间分布。常见做法是把数据重参考为平均参考（average reference），或者在顶区做双极导联。但无论选哪种，你都要清楚：参考改变的是脑电的绝对幅值和空间分布，不会改变某个频率“有没有峰”这一基本事实。换句话说，参考方式会影响 alpha 功率具体多大、哪个电极高，但不会把 alpha 峰整消失，也不会凭空造出一个 alpha 峰。

这里最需要注意的坑是：如果你分析的是整段平均功率谱，并且对多个被试做了不同参考方式，后果就是组间比较全部失真。所以，项目一开始就要统一重参考策略，最好在预处理流程里固定下来，别做到一半再改。

3.3 分段窗长、去趋势与数据重叠的关键逻辑

谱分析里的“分段”有两层意思。一层是 Welch 法把长数据切成若干小段，分别做 FFT 再取平均；另一层是指你在预处理时按事件或按时间窗截取数据。

在被试水平上做功率谱估计时，数据长度必须足够覆盖你关心的最低频率。我常用的经验公式是：如果你要可靠估计 1 Hz 以上的频带，每一段被分析的数据至少要有 2 秒；如果要区分 0.5 Hz 以下的慢波，至少要 4~6 秒。分段前最好把线性趋势去掉，否则直流漂移会在低频段制造大量虚假功率，尤其是 1 Hz 以下区域。去趋势可以直接用 detrend 函数，也可以做 0.1 Hz 或 0.5 Hz 高通滤波。但不要同时对同一段数据既做高通滤波又做 detrend，那样会过度压低低频，造成 theta/delta 功率被低估。

另一个容易忽略的问题是：Welch 法里的“重叠”是不是越多越好？答案是否定的。重叠率提高可以增加段数，从而降低方差，但段与段之间高度相关时，实际带来的自由度提升没有想象中大。一般 50% 重叠是工程上的经典选择，75% 重叠改善有限，计算量却显著上升。我实测过不同重叠率下 alpha 峰的稳定程度，50% 和 75% 几乎没有肉眼可见差别，所以我一般就固定在 50%。

4. 主流的三种谱估计方法：周期图、Welch、多锥度的取舍

4.1 三种算法到底在做什么

周期图法是最原始的思路：把一段 N 点信号做 FFT，然后取模平方。它干净、直接、没有任何额外参数，但方差很大。随性波动的脑电信号本身就是随机过程，单次周期图的每个频点估计值可能偏离真实功率很多。

Welch 法本质上是对周期图求平均。它把总长度信号分成多段，允许段与段之间有重叠，每段加窗后再做 FFT，最后把所有段的功率谱平均。这个平均操作直接把方差压下去了，代价是频率分辨率变差。如果你用 4 秒窗，频率分辨率大约是 0.25 Hz（实际与窗函数有关）；而如果你拿整段 60 秒数据做一次 FFT，分辨率可以达到 0.017 Hz。所以 Welch 是在“分辨率”和“稳定性”之间做权衡，绝大多数脑电静息态分析都会选择 Welch 法。

多锥度法（multitaper）是另一种策略：它对同一段数据使用多个正交的锥度窗函数（通常使用 Slepian 锥度，也叫 DPSS），每个锥度给出一个谱估计，最后再把这些估计平均起来。这么做的好处是，它能在保持较好频率定位精度的同时抑制谱估计方差，而且可以通过锥度参数控制频域平滑程度。多锥度法非常适合短数据、包含明显窄带节律（如 alpha 峰）或需要高精度频率定位的场景。它的缺点是参数理解门槛稍微高一点，常见的参数包括时间带宽积（time-bandwidth product）和锥度数量。

4.2 我的推荐参数表

在实际处理脑电数据时，我一般会根据分析目的做如下选择：

强调一个很常见的误解：多锥度法并不是“更高级、一定更好”。它引入的频域平滑会把原本很窄的 alpha 峰拉宽一点，如果你的分析目标是“把 alpha 和 beta 边界切清楚”，这个方法反而会让边界变得模糊。所以不要盲目追新，先想清楚你需要的到底是频率精度还是估计稳定性。

4.3 自由度与统计检验的关系

做组分析时，很多人会忽略 Welch 法里“段数”带来的自由度问题。一段 60 秒的数据，用 4 秒窗、50% 重叠，大概能得到约 29 个分段。这些分段不是完全独立，但因为相邻分段有重叠，实际独立段数要少于 29。做 t 检验或方差分析时，如果你把每个频点单独拿出来检验，多重比较问题会被放大；如果你把 alpha 频带（8~13 Hz）的平均功率作为一个指标，相当于对多个频点做了降维，统计稳定性会好很多。

这也是为什么我在实际项目中，很少直接报告“某个频点上有显著性差异”，而是报告“某个频带的平均功率存在差异”。频带平均既符合神经生理学原理，又能在一定程度上避免多重比较和自由度过低带来的假阳性。

5. 窗函数、重叠率和 FFT 点数：调参的规则与边界

5.1 频率分辨率其实是个“窗户”问题

窗函数这个词听起来很抽象，实际上可以这样理解：你把一段长脑电信号切成 4 秒的小段，相当于用一个“矩形窗户”把无限长的信号截了一段出来。但矩形窗在频率域里的表现很差，旁瓣很高，会让一个强 alpha 峰的功率“泄漏”到相邻频点，形成很宽的谱拖尾。为了抑制这种泄漏，我们改用汉宁窗、汉明窗或者布莱克曼窗，它们的特点是中间高、两边低，段中心和段边缘样本的权重不同，从而让截断更平滑。

代价也很明确：加窗会让主瓣变宽，频率分辨率变差。矩形窗理论上可以把谱峰压得最窄，但它旁瓣高；汉宁窗主瓣宽度大约是矩形窗的两倍，但旁瓣大幅下降。对脑电这种多频率混合的信号来说，泄漏带来的误差通常比主瓣变宽更麻烦，所以我默认首选汉宁窗，特殊场景用 Slepian 窗口。

5.2 窗长、重叠率、FFT 点数三者联动怎么调

窗长决定了频率分辨率，FFT 点数决定了谱线“画得多细”。一个常见误区是：为了把谱线画得更高清，把 FFT 补零到很大点数，比如 8192 点。但补零并不会提高真实频率分辨率，它只是插值，让曲线看起来更平滑。真实的频率分辨率仍然由窗长决定。

所以，正确的调参优先级应该是：

1. 先确定分析频段，计算最低关心频率，从而初步确定最小窗长。比如最低关心 1 Hz，窗长至少 2 秒。
2. 再确定你要区分的最小频率间隔。如果想区分 9 Hz 和 10 Hz 的两个邻近峰，频率分辨率至少要优于 1 Hz，所以窗长应大于 1 秒；如果想精细测量 alpha 峰，让 4 秒甚至 6 秒更稳妥。
3. 根据数据总长度选择重叠率。数据很短时，尽量提高重叠率以增加段数，比如 1 分钟数据搭配 4 秒窗、50% 重叠；数据很长时，重叠率不需要刻意调高。
4. 最后设置 FFT 点数。通常设为大于窗长点数的最小二次幂即可，或直接保留与窗长相同点数，不会影响结果。

5.3 我自己踩过的窗函数坑

曾经有一次分析静息态数据，随手用了默认参数，结果所有被试的 alpha 峰都宽得离谱，峰值频率还明显右移。排查了很久，发现是窗长设成了 1 秒，频率分辨率只有 1 Hz，alpha 峰 8~13 Hz 本来就只占几个频点，再叠加窗函数展宽，谱峰自然糊成一片。换成 4 秒窗之后，alpha 峰干净得多。这件事让我以后形成习惯：任何分析流程上线前，先人工抽 3~5 个被试画原始谱图，看一眼峰形是否合理，再进入批量处理。

6. 全相位数字谱分析：一个来自工程领域的补充方案

6.1 全相位思路和传统 FFT 的差别

“全相位数字谱分析方法”这几年在信号处理圈里讨论得不少，它的想法是：传统 FFT 只取一段 N 点信号，起始相位不同会直接影响 FFT 结果的相位和泄漏特性；而全相位方法会把所有可能包含某个采样点的 N 点截断序列都考虑进来，并把它们的相位对齐后再做加权平均，最后只取中心样本对应的谱值。这样做的好处是相位估计非常稳定，同时抑制频谱泄漏的能力比普通加窗 FFT 更好。

如果你要自己实现一个最简版本，流程大概是这样：对长度为 2N-1 的数据，设计一个长度为 N 的窗函数（比如汉宁窗），用这个窗和自身反折卷积得到全相位窗，把 2N-1 个点加权后按周期延拓求和，再取 N 点做 FFT，最后乘以校准系数得到谱。工程实现上，MATLAB 里可以直接按这个逻辑写，Python 也可以用 numpy 实现。由于它不是数学软件内置的标准函数，我建议不要在生产管线里随便替换 Welch 法，而应先把它作为“对照方法”用来验证窄带峰的位置和相位。

6.2 这个方法的适用边界

在脑电场景里，全相位数字谱分析的真实优势主要体现在两类情况：一是你对某个窄带节律的相位非常敏感，比如做锁相分析、稳态诱发电位；二是数据里有很强的高频干扰，你想更干净地分离窄带峰。但它也有明显短板：它本质上没有做多段平均，对单段噪声数据的抑制能力并不比周期图强多少；如果你拿它估计静息态全频带功率，曲线仍然会很毛糙。所以我的建议是：日常全频带功率比较仍用 Welch 或者多锥度，全相位谱可以放在特定问题里做交叉验证，不要因为它名字听起来“高级”就盲目改用。

7. 读谱阶段的常见误区：频带、伪迹、参考与体积传导

7.1 频带划分不能只背数字

经典的频带划分是 delta（1~4 Hz）、theta（4~8 Hz）、alpha（8~13 Hz）、beta（13~30 Hz）、gamma（30 Hz 以上）。但这只是通用模板，个体差异很大。同一个被试在不同记录日的 alpha 峰值可能漂移 1~2 Hz；同一个被试在睁眼和闭眼条件下，alpha 峰值也可能不同。如果机械地把 alpha 固定在 8~13 Hz，一部分被试的 alpha 峰可能在 7.5 Hz 或 13.5 Hz，你的频带平均功率就会把他一半的节律能量切到 theta 或 beta 里去。

我常用的做法是：先画个体功率谱，确定可视化峰值；再做频带划分时考虑个体峰值，比如以峰值频率 ±2 Hz 作为个体的 alpha 频带；最后做组分析时，同时报告固定频带和个体化频带两种结果。个体化频带对效应量的提升通常非常明显。

7.2 体积传导和参考会让频谱解读翻车

脑电记录在头皮上，每一个电极收集到的都是全脑大量神经元同步活动在空间上叠加后的结果，这叫做体积传导效应。因此，当你看到枕区电极（如 O1/O2/Oz）上 alpha 功率最高时，你说的不应该是“枕叶产生了 alpha”，而应该说“这些传感器上方记录到了最强的 alpha 功率”，或者最多谨慎地说“alpha 的分布集中在枕区”。想讨论真正的源位置，需要做偶极子定位或者源重建。

另一个很容易被忽略的问题是分析单电极功率谱时的参考污染。使用某一只电极做参考时，参考电极本身也可能混入很强的低频活动，导致全脑功率谱的低频段都被抬高。我习惯做平均参考，但即使这样，也要在文章中写清楚参考方案，因为参考不同，功率谱的绝对数值在不同实验室之间很难直接对比。

7.3 不要把所有低频或高频分量都当成神经信号

读谱时我会先问自己一个问题：“这个峰是脑来源，还是来源很明显？”如果 alpha 峰在枕区、睁眼时降低，那大概率是脑节律；如果 50 Hz 附近有尖峰且各电极普遍存在，那是工频；如果 gamma 频段功率在所有前额电极普遍偏高，那基本是肌电。低频段如果出现一个非常锐利的窄峰，还要考虑是否来自参考电极的微弱漂移或电极线运动。

为了减少误判，我建议在正式分析前做一个“伪迹谱对照”：选一段明显包含眨眼和明显包含肌肉伪迹的数据，分别算功率谱并保存为模板，之后看到类似形态时就能第一时间对齐。

8. 实操工作流：一个静息态 alpha 谱分析的完整示例与问题速查

8.1 最小可复现流程

下面我按自己最常用的一套流程写一个示例。假设数据是 64 导脑电、采样率 1000 Hz、时长 60 秒的闭眼静息数据。处理目标是计算 O1/O2/Oz 的平均 alpha 功率，并对比两个条件。

第一步，读取数据并做基本预处理。使用 MNE-Python 的话，大致的操作链路是：导入数据 → 定位坏导 → 0.5~40 Hz 带通滤波 → 运行 ICA 去掉眨眼和心电成分 → 参考改为平均参考 → 剔除明显伪迹段。这一步不要偷懒，因为谱分析对低频漂移和高频肌电都很敏感。

第二步，分段与计算。取每段 4 秒、50% 重叠，对每一段加汉宁窗后做 FFT，再平均所有段的功率谱。MNE 里可以直接这样写：

import mne from mne.time_frequency import psd_welch epochs = mne.make_fixed_length_epochs(raw, duration=4.0, overlap=2.0) psds, freqs = psd_welch( epochs, fmin=1.0, fmax=40.0, n_fft=4096, n_overlap=int(2.0 * raw.info["sfreq"]), window="hann", average="mean" )

第三步，提取枕区电极 alpha 频带平均功率。写一个很简单的频带平均逻辑：把 8~13 Hz 范围内所有频点的功率取平均。由于 alpha 个体峰值有漂移，我建议同时用可视化方法确认 8~13 Hz 内确实有峰，再决定是否做个体化频带。

import numpy as np alpha_idx = np.logical_and(freqs >= 8, freqs <= 13) alpha_power = psds[:, [ch_names.index("O1"), ch_names.index("O2"), ch_names.index("Oz")], alpha_idx].mean(axis=(1, 2))

第四步，根据统计设计做组间或组内比较。由于功率谱是非负、右偏的，一般先做以 10 为底的对数转换或除以平均功率做相对值，再进入 t 检验或方差分析。千万别直接拿原始 μV²/Hz 去做正态性很强的统计，除非你的数据已经显示正态。

8.2 我常用的问题排查速查表

这套工作流看起来很简单，但实战中最大的成本往往不是代码，而是“判断”。什么时候该剔除一个段，什么时候该保留；一个信号峰到底是真实节律还是参考电极干扰，都需要反复看图积累经验。我的建议是：头几次做谱分析，把所有中间结果图都画出来，一张一张过。不要直接跳到统计结果。

最后，如果你后续要分析任务态脑电或者需要追踪某一频带功率随时间的变化，谱分析只是起点，必须过渡到时间分辨的时频分析。这也是我接下来系列文章要重点写的方向。