做脑电分析这几年我越来越觉得频域这关是绕不过去的坎。时域波形图谁都会看但一到“这个被试的 alpha 波到底比那个被试高多少”“睡眠分期里 delta 功率怎么变化”“抑郁症个体是不是 frontal 区的 theta 功率异常”这类问题全部都要回到功率谱密度PSD上找答案。可以说只要是跟脑电节律、神经振荡、频段能量有关的分析PSD 就是地基中的地基。这篇是脑电分析系列的第 7 篇我把功率谱密度的计算方法和常用方案一次性讲透。不堆公式吓唬人但关键原理必须说清不搞“抄个函数就完事”的糊弄文但代码和参数会给你能直接用的版本。适合刚接触脑电频域分析的研究生、刚接手 EEG 数据分析的工程师以及想把 Welch、AR 模型、多窗口法这些方法区别搞明白的从业者。看完你至少能回答三个问题PSD 到底在算什么不同方法差在哪我自己处理脑电数据时该选哪种。1. 频域分析为什么绕不开 PSD1.1 从时域波形里看不到的那些信息脑电信号本质上是一堆不同频率正弦波的叠加。时域图上你看到的是一条连续起伏的曲线但如果把这条曲线拆开来看里面其实藏着 1 Hz 的慢波、8-12 Hz 的 alpha 节律、20-30 Hz 的 beta 成分以及更高频的 gamma。我们的大脑在清醒、困倦、睡眠、认知任务等不同状态下这些频率成分的能量会发生变化而这些变化在时域波形里往往是“看不出来”的——看起来差不多的噪声样波形其频谱结构可能完全不同。频域分析的核心目标就是把时域信号里各个频率成分的相对大小找出来。比如同样是闭眼静息态的两段脑电一段 alpha 波明显占优另一段却几乎没有 alpha 节律二者的时域波形都是杂乱无章的曲线但频率结构截然不同。PSD 解决的就是这种“看不出但实际存在”的问题。1.2 功率谱密度到底在算什么功率谱密度描述的是信号功率在不同频率上的分布情况。脑子里可以想一个画面总功率是固定的PSD 告诉我们这些功率是按什么比例分配给不同频率的。在这个频率上功率高说明该频率对应的神经振荡活动强那个频率功率低说明这个频段的活动弱。从数学上看脑电信号属于随机信号它没有像 sin(t) 那样的确定性傅里叶变换所以要用统计意义上的“功率谱”来描述。我这里只点出最核心的关系对一段信号 x(t)它的自相关函数与功率谱密度是一对傅里叶变换对这就是维纳-辛钦定理。另个一常见的理解角度是PSD 可以写成对信号加窗后做傅里叶变换再取模平方的某种平均。后者更直观也是工程实现中最常用的思路。实际计算中最常见的方法是周期图法及其改进版。脑电分析里Welch 法出现频率最高因为它在计算稳定性和结果可读性之间找到了很好的平衡点。后面我会对这几种方法逐个展开。1.3 脑电频段的划分是PSD分析的落脚点做 PSD 不是为了画一条漂亮的曲线而是为了把结果对应到脑电的经典频段上去。目前广泛使用的频段划分大致为delta0.5-4 Hz深度睡眠、某些疾病状态、婴幼儿发育期占优theta4-8 Hz困倦、冥想、工作记忆加工、海马节律alpha8-13 Hz闭眼静息、放松状态、抑制控制相关beta13-30 Hz警觉、运动执行、认知活跃gamma30-45 Hz 或更高高级认知整合、跨脑区信息绑定有了频段划分PSD 曲线就变成了可解释的指标——把 alpha 频段范围内的功率积分起来就得到 alpha 功率把 alpha 功率除以总功率就得到 alpha 相对功率。这种频带功率指标比直接描述整条曲线要实用得多。所以在动手计算 PSD 之前建议先问自己一个问题我是想看整段信号的频谱全貌还是只需要某个频带的具体功率值这个问题的答案直接决定你后续的方法选择和参数设置。2. PSD 计算的五大核心方法与对比2.1 周期图法最朴素但方差很大周期图法是最直接的思路对整段信号直接做傅里叶变换然后把幅值取平方再除以信号长度。对应到代码层面SciPy 里的signal.periodogram函数就是干这个事的。但周期图法有个致命伤——估计方差大。单独一次周期图得到的 PSD 曲线会抖得厉害相邻频率点的功率值忽高忽低基本没法直接用来做组间比较。为什么会这样呢因为瞬时频谱的波动性太高一次观测对随机过程来说样本量不够估计结果就不稳。在脑电分析中如果不分段、直接用整段 5 分钟的静息态数据算一个周期图你会发现曲线参差不齐但大的节律峰比如 alpha 峰还是能看出来的。所以周期图可以作为“快速看一眼频谱长什么样”的工具但不建议作为正式统计的指标来源。2.2 Welch 平均法脑电分析的首选方案Welch 法是周期图法的高级版。核心思想很简单把一段长信号切成一堆小段每段分别做周期图然后把所有小段的功率谱取平均。这样做方差显著降低代价是频率分辨率下降——因为你每段长度短了能分辨的最小频率间隔就变大了。Welch 法有两个关键参数窗函数长度和重叠率。窗函数长度决定频率分辨率比如采样率 500 Hz、窗长 2 秒的信号每个窗内有 1000 个点做 FFT 后频率分辨率为 0.5 Hz。重叠率影响估计的平滑程度和计算量脑电分析里 50% 重叠是默认配置有些场景会用 75% 来进一步平滑。为什么脑电分析特别适合 Welch原因有三一是脑电信号长通常至少几分钟分段后每段基本可以视为近似平稳二是 Welch 法硬件计算成本低处理几十个通道的长时间数据毫无压力三是它对噪声的抗性比单次周期图好很多。2.3 Blackman-Tukey 法从自相关的角度切一刀Blackman-Tukey 法走的是维纳-辛钦定理这条路先估计信号的自相关函数再对自相关序列加窗后做傅里叶变换得到功率谱。逻辑上它和周期图是等价的区别在于加窗的位置不同——周期图在数据上加窗BT 法在自相关函数上加窗。BT 法在脑电分析中的使用率不高原因在于计算自相关函数的开销较大而且为了减小方差必须对自相关函数截短这会引起功率谱估计的偏差和频率泄漏。但在一些特定场景——比如对短数据段、或需要从谱分解里提取参数时——BT 法仍有它的用处。日常脑电流水线分析不推荐首选 BT 法了解原理能帮你更好地理解那个 PSD 公式为什么能成立就够了。2.4 AR 模型法参数化估计的粗犷版AR 模型法自回归模型Auto-Regressive的处理思路换了个赛道。它不再直接对数据做傅里叶变换而是先假设脑电信号符合一个由自身历史值线性组合生成的模型当前时刻的值等于前 p 个时刻值的加权和再加上一个白噪声激励。用这个拟合出的模型参数可以推导出信号的功率谱表达式。如何确定阶数 p这是 AR 法的核心问题。阶数太低谱峰会被抹平阶数太高会出来一堆假峰。常用的选择准则包括 AIC赤池信息准则和 BIC贝叶斯信息准则。脑电中常见的选择范围在 10 到 50 之间具体和信号的复杂度有关。AR 模型法的优势在于对短数据段的谱估计效果好谱峰尖锐且平滑。因此它在睡眠脑电分期、某些事件相关频谱分析中仍有应用。缺点是模型阶数选择要靠经验错了结果偏差很大而且如果信号本身不符合 AR 模型假设估计出来的谱会有系统性偏误。现在大多数通用脑电分析流程已转向 Welch 法AR 模型更适合有特定需求的研究者。2.5 多窗口法MTM与子空间法骨灰级但也别忽略多窗口法Multitaper MethodMTM用一组正交的 Slepian 塔珀窗DPSS 窗分别对信号加窗得到多个独立的周期图再求平均。它在低频段和窄带信号的功率估计上优势明显谱估计方差小、能量泄漏少。缺点是实现复杂、计算慢参数时间带宽积 NW 和塔珀数需要调。脑电研究里想精细分析低频慢波比如 0.1-1 Hz 的皮层慢电位时MTM 的表现会比 Welch 更好。子空间类方法如 MUSIC、ESPRIT原本是阵列信号处理和雷达领域的主力通过将信号分解为信号子空间和噪声子空间来估计频率成分。它们的频率分辨率极高可以分辨靠得很近的谱峰但计算复杂度高、对模型阶数敏感。在脑电分析中用得极少主要出现在特殊科研场景里比如需要精确提取特定频率振荡的相位和幅值。为了直观对比我做了一个方法性能速查表方法频率分辨率方差水平计算成本脑电适用场景缺点周期图法高整段信号很高低快速预览频谱形态曲线抖容易出假峰Welch 平均法中受窗长限制低低常规脑电频谱分析首选高频细节相对弱化Blackman-Tukey中中中自相关域分析场景堆计算使用率低AR 模型法高模型分辨率好低中短数据段、稳态谱估计阶数选择敏感怕模型假设不成立MTM 多窗口法高低高低频精细谱分析、慢波研究调参复杂运行慢子空间法极高中高少数特殊科研场景阶数难定脑电里很少用3. 脑电 PSD 实操全流程与代码实现3.1 预处理环节容易被忽视的影响因素预处理决定了 PSD 的质量。这里面有个关键认知功率谱分析对数据的平稳性和干净程度极其敏感。如果信号里有大幅的眼动伪迹低频段尤其是 delta 和 theta的功率会被严重抬高如果肌肉紧张度高高频段beta 和 gamma会被肌电污染。这也是为什么先做 ICA 去眼电/肌电、再做 PSD 是标准套路的原因。具体到预处理流程我建议至少做到这几步重参考常用平均参考或双侧乳突参考带通滤波0.5-45 Hz 或 1-40 Hz取决于你关注的频段去坏段人工或算法剔除幅值突变明显的时段ICA 去伪迹去掉眨眼、眼动、心电、肌电成分去均值、去线性趋势消除直流偏置和基线漂移的影响去均值这一步极其重要。如果不去均值直流分量会在 0 Hz 附近产生一个巨大的能量峰还会通过泄漏影响低频段。去线性趋势则是为了抑制低频漂移这种漂移在原始脑电中非常常见会大幅抬高 delta 频段的功率。3.2 参数选择窗长、窗函数、重叠率怎么定在脑电分析中Welch 法的参数设置值得专门讲。我自己的经验是先明确你要的频率分辨率再反推窗长。频率和时间的分辨率之间存在刚性约束纱窗效应在脑电里同样存在。窗越长频率上看得越细窗越短时间上定位越准。做静息态或睡眠分析时关心的是某个时间段内的平均功率这时可以放宽时间分辨率选择 2 秒或 4 秒窗长既能把 alpha在 10 Hz 附近和 theta6 Hz 附近这类节律峰解析清楚又能保持足够的平滑度。窗函数建议选汉宁窗Hann。对比矩形窗汉宁窗的旁瓣更低能显著减少频域泄漏。注意加窗会让频谱幅值衰减所以用 Welch 法算出来的 PSD 需要做归一化校正。SciPy 的scipy.signal.welch内部已经做了正确归一化直接使用即可。重叠率方面50% 的默认值是稳妥的选择。如果你的数据段长度有限、希望曲线再平滑一点可以提高到 75%。但重叠率过高会让相邻窗高度相关实际增加的信息量有限反而增加计算量一般不建议超过 90%。3.3 用 Python 完成 PSD 计算与频段功率提取下面给出一个可以直接套用的示例代码用的是 MNE 库处理数据、SciPy 计算 PSD。假设你已经有清洗后的连续脑电数据raw并且已经设定好电极位置和采样率。import numpy as np from scipy import signal import mne # 读取清洗后的数据 raw mne.io.read_raw_fif(preprocessed_data.fif, preloadTrue) sfreq raw.info[sfreq] # 采样率通常 250/500/1000 Hz # 提取数据矩阵形状为 (n_channels, n_times) data raw.get_data(pickseeg) * 1e6 # 转为微伏单位 # Welch 法计算 PSD freqs, psd signal.welch( data, fssfreq, npersegint(2 * sfreq), # 窗长 2 秒 noverlapint(1 * sfreq), # 50% 重叠 nfftint(2 * sfreq), windowhann, axis-1, ) # freq: (n_freqs,), psd: (n_channels, n_freqs)单位是 uV^2/Hz # 定义频段 bands { delta: (0.5, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30), gamma: (30, 45), } # 计算各频段的绝对功率与相对功率 band_power {} for band_name, (f_lo, f_hi) in bands.items(): mask (freqs f_lo) (freqs f_hi) band_power[band_name] psd[:, mask].sum(axis-1) # 积分近似 total_power sum(band_power.values()) relative_power {name: power / total_power for name, power in band_power.items()} # 打印示例前 5 个通道的 alpha 相对功率 print(Alpha relative power (first 5 channels):, relative_power[alpha][:5])这段代码跑出来的psd是二维数组每一行对应一个通道。如果你要画某个电极比如 Cz的 PSD 曲线只需要把对应行的数据与freqs一起画出来纵轴用对数刻度即可。如果用的是 MNE 自身的功能也有现成接口spectrum raw.compute_psd(methodwelch, fmin0.5, fmax45, n_fftint(2 * sfreq)) spectrum.plot() band_power_mne spectrum.get_band(pick_bands{alpha: (8, 13)})用compute_psd会更省心因为它会直接处理raw对象里的元数据输出也是 MNE 的频谱容器类型方便后续批量统计和可视化。3.4 结果怎么解读才是对的算出频段功率之后别忘了做“归一化再比较”这一步。绝对功率在不同人之间差异巨大受头皮厚度、电极阻抗、颅骨传导特性影响组间比较通常用相对功率某频带功率占总功率的比例或功率比值如 theta/beta ratio来抵消个体差异。另一方面单次测量得到的 PSD 只是该时段内的大脑状态快照。静息态 alpha 功率在闭眼时明显高于睁眼睡眠和清醒状态下频谱差异更大。所以解读结果前一定要问清记录条件被试是睁眼还是闭眼是否执行了任务数据采集时长是多少睡眠分期的信息有没有带上基线校正也是常见操作。如果是任务态研究通常会减去任务前静息期的 PSD得到“事件相关频谱扰动”或“任务相关功率变化”。这个过程本质上就是逐频率点相减或相除操作上非常简单但解释上要谨慎——它反映的是从基线到任务状态的功率变化不是绝对的功率高低。4. 实操中的典型问题与避坑实录4.1 频谱泄漏和窗函数选择频谱泄漏是指本来集中在一个频率上的能量因为截断效应“漏”到相邻的频率点上去了表现为 PSD 曲线在谱峰附近出现拖尾。泄漏严重时会淹没有价值的低幅峰也会让相邻频段的功率估计互相污染。解决方向有两个一是用非矩形窗推荐汉宁窗二是保证窗长是信号周期内容的合理倍数。脑电信号频率成分复杂不可能找到精确的倍数关系所以主要靠窗函数来压制泄漏。如果你发现曲线高频段有种“波浪状”的起伏多半是泄漏在作怪更换窗函数往往能改善。另外在某些实时分析系统中为了省资源会用矩形窗。这个做法对频谱形态要求不高的场景能接受但你要做正式的频带功率统计就别省这一步。4.2 基线漂移对低频功率的污染低频段0.5-2 Hz的功率统计最让人头疼的问题就是基线漂移。电极与头皮接触阻抗的慢变化、放大器的直流漂移、受试者的缓慢动头都会引起低频大幅摆动。这种漂移反映在 PSD 上就是低频端能量急剧抬升直接污染 delta 段。处理手段包括高通滤波0.5 Hz 或 1 Hz 以上、去趋势线性或多项式、以及在预处理中剔除长时段漂移明显的片段。注意高通滤波的截止频率不能设太高否则会真实削减 delta 和 theta 的生理成分。0.5 Hz 高通是脑电分析里最常见的折中方案。4.3 伪迹对 PSD 的影响远比想象中大眼电伪迹的能量集中在低频段主要在 0-4 Hz会严重高估 delta、theta 功率肌电伪迹则散布在 20 Hz 以上的频段会让 beta 和 gamma 功率虚高。也就是说伪迹对 PSD 的污染是有频率特异性的不清理干净最后统计结果很可能只是“伪迹功率组间差异”的翻版。ICA 是清理脑电伪迹的常用方式但也有它的局限——如果某个独立成分同时包含眼电和脑电活动分裂不彻底就会在分离过程中损失部分脑电信号。我的建议是预处理阶段宁可在时间上多剔除一些疑似伪迹段也不要留着硬扛。有问题的时窗数据留着不进 PSD 计算比后期任何补救都有效。4.4 绝对功率和相对功率的选择焦虑这个问题几乎每次组会都会被问一遍到底用绝对功率还是相对功率我个人的经验是单被试纵向比较比如同一个患者治疗前后的对比可以用绝对功率因为记录条件基本一致跨被试横向比较务必用相对功率或比例指标。绝对功率的个人差异太大了头皮形态、颅骨电阻抗、电极摆放细微差异都会带来系统性的水平差异直接比较很容易被个体间的“基准”差异掩盖了真正的神经活动差异。不过相对功率也有坑。由于各频段功率总和是固定的一个频段相对功率上升必然伴随至少一个其他频段相对功率下降。在解读结果时,需要小心“频段补偿效应”——比如额叶 theta 相对功率升高可能会造成 alpha 相对功率看似下降但这并不意味着 alpha 的绝对功率真的减少了。忠告是重要发现最好同时报告绝对功率和相对功率两种指标并注明所用的换算方式。4.5 频段边界和参数一致性问题另一个隐蔽的坑是频段边界定义不统一。同样是 theta 频段有的文献定义为 4-8 Hz有的用 4-7 Hzbeta 有 13-30 Hz也有 15-30 Hz。如果你做文献复现务必把原文的频段定义copy到自己的脚本里别想当然地挑一个“看起来差不多”的。频段差 1 Hz在小样本研究中完全可能影响显著性。另外分析参数窗长、重叠率、窗函数、滤波范围在组间比较时必须在所有被试上保持一致。哪怕只是把某个被试的数据处理流程从“2 秒窗”换成“4 秒窗”都会引入本来不存在的组间差异。做批量分析时把全部参数写进一个配置文件或字典里这是个普通但很管用的习惯。现象可能原因排查与解决方案PSD 低频端异常上翘基线漂移未去除加高通滤波、去趋势剔除漂移段高频段波浪状起伏频谱泄漏换汉宁窗适当增加 FFT 点数谱线抖动剧烈段数太少或未平均减小窗长增加段数或提高重叠率alpha 峰不明显睁眼状态记录确认记录条件必要时分组对比同一被试两次 PSD 差异大预处理参数不一致固定全部参数用配置脚本统一跑4.6 从相反视角认识 PSD 的局限性PSD 也不是万能的。它假设信号在一定时长内是平稳的而脑电实际上充满了瞬态变化。对于任务态、事件相关态的动态变化PSD 往往会“平均掉”关键的时间演变信息。这个时候更适合用时频分析比如短时傅里叶变换STFT、小波变换或希尔伯特黄变换。日常分析里我总结出来最简单的判别标准如果问题是“大脑在这段时间内的平均振荡状态是什么样”用 PSD 就够如果问题是“某个刺激出现后 300 到 500 毫秒内 alpha 功率怎么变化”请去用时频分析。还有一点关于参考电极的提醒。不同参考策略平均参考、双侧乳突参考、CZ 参考对 PSD 的影响是全局性的尤其影响顶区和枕区 alpha 的功率幅度。同一个数据集用不同参考处理alpha 功率的绝对值可能相差数倍。所以做跨研究比较时参考方式必须匹配做组内比较时至少保持同一实验使用同一种参考。5. 从 PSD 到常用衍生指标5.1 峰值频率与频带中心频率在某些疾病研究里PSD 的形态参数比能量参数更敏感。比如帕金森病患者服用多巴胺药物后beta 频段13-30 Hz的峰值功率会下降且峰值频率是否发生偏移也有临床意义。提取峰值频率的步骤是在目标频段内找 PSD 最大值对应的频率点。做这个操作时要注意数据平滑——如果直接用原始 PSD 找峰值容易被噪声干扰导致频率点反复横跳。建议先对 PSD 做一次平滑Savitzky-Golay 滤波或移动平均再定位峰值。峰值频率这个指标在睡眠纺锤波、alpha 个人频率差异等研究中都是核心输出。不过它对计算参数也比较敏感窗长偏短时频率分辨率不足峰值定位偏差会比较大所以要保证该频段内有足够的频率分辨率。5.2 谱陡度spectral slope近年来脑电领域对宽频功率broadband power和谱陡度非常感兴趣。这套思路认为 PSD 在对数-对数坐标上服从一个近似幂律的形式即功率与频率的某个幂次成正比。对数坐标下这条关系表现为一条直线直线的斜率叫谱陡度spectral slope截距则反映整体的宽频功率水平。谱陡度与神经元放电活动、兴奋抑制平衡有关也受年龄、药物、疾病状态影响。计算谱陡度的过程比想象中麻烦要去除周期性的振荡峰只拟合非振荡背景再对拟合结果结合线性回归提取斜率和截距。现在有一些开源工具如 speclib、FOOOF 算法专门做这件事大大降低了计算门槛。如果你们团队的方向偏向神经机制研究这个指标值得纳入分析计划。5.3 功能连接里的 PSD 参与在静息态脑电功能连接分析中PSD 常作为加权项参与计算。比如加权相位滞后指数wPLI或加权相干Weighted Coherence都需要先估计每对导联之间的互功率谱cross-spectrum和各自的 PSD。这个场景下PSD 不仅是结论更是中间产物。此时对 PSD 的准确度要求更高因为相位和幅度的偏差都会传递到连接指标里。如果做过连接分析你会遇到一个实际问题到底该用单次周期图、Welch 还是 MTM我的经验是Wiener 型连接指标比如功率谱相干用 Welch 的平均谱比较稳定而基于相位连接(如 PLV、wPLI)的则对相位一致性更敏感用 MTM 或分段处理后的相位信息效果更稳。心理物理上这种差异不明显但实际计算时你会发现结果可能显著不同。5.4 任务态与静息态的 PSD 使用差异静息态分析关注的是“基线状态下大脑各个频带的功率格局”一般把整段数据计算得到一条 PSD再提取各频带功率。任务态分析则复杂得多一是时间锁定的瞬态功率变化需要短窗分析二是任务引起的功率变化往往表现为事件相关同步化ERS和事件相关去同步化ERD需要对比任务段和基线段的功率。经典做法是对每个 trial取刺激前 500 ms 到刺激后 1000 ms 之间做时频分解然后在频域上对比基线。注意了只用 PSD 做这种对比是有风险的——它只能告诉你功率变了多少不能告诉你这个变化锁定的时间位置。结合 STFT 短窗滑动或小波变换才能既保留时间信息又得到频率信息。这也是很多人把时频分析和 PSD 搭配使用的原因。6. 写在最后的实操体会做脑电频域分析这几年我踩过的最大的坑不是算法不会写而是参数设置和预处理不一致导致结果不可复现。同一个数据集换个人来跑因为滤波器截止频率不同、窗长不同、重叠率不同最后的频带功率可能完全不同。所以我现在做任何批处理之前第一件事就是写一个元数据文件记录所有分析参数包括版本号。这个习惯救了我很多次。另外我强烈建议把 Welch 作为常规分析的首选但不要排斥其他方法。早期的研究工作里用周期图法被审稿人质疑方差问题后来换用 Welch 后结果稳定了一大截。等做低频慢波研究时又发现 Welch 的低频端精度不够换成 MTM 才看清楚 1 Hz 以下的谱结构。不同方法适合不同问题方法本身没有绝对的优劣适合你的数据和科学问题最重要。如果你还拿不准手里的数据该用哪种方法给一条最简单但不踩坑的建议先画一张不重叠、无平滑的周期图看看图谱形态再画一张 Welch 对比一下两张图放在一起你基本就能判断信号里哪些是稳定的振荡峰哪些是噪声和伪迹。多对比几次你对 PSD 的理解会超过很多照猫画虎的脚本本身。