做脑电数据分析的，迟早会撞上“谱分析”这堵墙。无论是看静息态alpha节律有没有增强，还是算事件相关任务里的theta/beta能量变化，频谱分析都是第一个绕不开的工具。我最早接触脑电时，拿到一段波形图就懵了，密密麻麻的曲线根本看不出门道，直到把数据从时域变换到频域，那些隐含的节律才真正浮出水面。这篇就拆解谱分析这条主线，从FFT、功率谱密度讲到Welch平均、短时傅里叶变换，把常用的参数怎么选、数据怎么算、结果怎么解释一次性理清楚。适合刚接触脑电数据处理、想尽快看懂频谱图并独立跑通分析流程的读者。
1. 为什么脑电要先做谱分析——从时域看不出名堂的波形说起
1.1 时域信号里藏的信息，靠眼睛很难看全
脑电本质是一串随时间变化的电压信号，采集下来就是典型的时间序列。但人眼在时域上能读到的信息非常有限，我们都知道alpha波大概8到12赫兹，睁眼闭眼时幅度会有变化，可是这事儿靠盯着波形图肉眼判断一点都不靠谱：信号里的噪声、肌电伪迹、基线漂移叠在一起，某一段波形看着“抖得厉害”，你很难说清楚这到底是alpha节律变强了，还是眨眼伪迹带来的低频扰动。
这就是谱分析存在的根本原因。它做的事情，说穿了就是把“随时间变化的电压”拆成“不同频率成分的叠加”，告诉我们每一段信号里，4赫兹的theta有多少能量，10赫兹的alpha有多少能量，40赫兹的gamma又有多少。从频域去看脑电，等于换了一副眼镜，很多在时域里模模糊糊的特征会变得非常清晰。
我在实际项目里最常用谱分析的场景有三个：静息态闭眼和睁眼两种条件下的alpha功率对比、睡眠分期里的delta/theta/spindle检测、以及事件相关任务中特定频段能量的动态变化。这三种场景基本覆盖了目前脑电研究的主流需求，也解释了为什么谱分析会成为脑电数据处理流程里默认的第一步。
1.2 脑电频段划分：不是拍脑袋的，背后有生理意义
做谱分析之前，先得把脑电的频段划分搞清楚。标准的分法是把频谱分成几个区间，delta（1到4赫兹）主要是深睡眠状态的优势节律，theta（4到8赫兹）和记忆编码、情绪加工关系密切，alpha（8到13赫兹）在放松闭眼时最明显，beta（13到30赫兹）与主动思考、运动准备相关，gamma（30到45赫兹甚至更高）通常和跨脑区的信息整合绑定。
不同文献的边界会略有浮动，比如有的写alpha是8到12赫兹，有的写8到13赫兹，这都不算错，关键是你在项目里得固定一套标准，不能这篇用A标准、那篇用B标准。频段选定之后，所有受试者、所有条件都用同一套边界，否则结果根本没有可比性。
我在分析里习惯把频段边界做成一个字典存着，每条数据算完直接查表取值，这样既保证一致性，又方便后期调整。需要注意的是，低频段（尤其是1赫兹以下）往往混入大量漂移噪声，如果不做高通滤波就去算谱，delta频段的能量会被严重高估，后面会专门讲这个问题。
1.3 谱分析能回答哪些实际问题
谱分析最直接的产品就是频谱图，横轴是频率，纵轴是功率或幅值。这张图可以回答几个层面的问题。第一，这段脑电的主导节律是什么？比如静息态下在10赫兹附近看到一个明显的峰，基本可以判断这个人的alpha节律很突出。第二，不同条件下的节律强度差别有多大？闭眼时alpha峰明显高于睁眼，这个在频谱图上几乎是一眼可见的差异。第三，某种疾病或干预手段是否改变了特定频段的能量？比如用tACS在10赫兹做刺激，理论上应该提高alpha功率，频谱分析就是检验这种效应的标准方法。
除了静态频谱，谱分析还有时间维度上的扩展，也就是时频分析。一般脑电实验里，刺激出现之后大脑的反应并不是平稳的，theta可能在刺激后300毫秒短暂增强，alpha可能随后被抑制，这种动态变化单靠一张平均频谱图根本看不出来，必须把时间切成小窗、逐段计算频谱，这就是后面要讲的短时傅里叶变换。
2. 谱分析的数学地基：FFT与功率谱密度
2.1 FFT到底干了什么活
说起频谱分析，FFT（快速傅里叶变换）是躲不开的。傅里叶变换的核心思想是：任何一段信号，都可以看作若干不同频率、不同振幅、不同相位的正弦波的叠加。FFT就是高效计算这种分解的算法，它把一个长度为N的时间序列变换成同样长度的复数数组，每个复数对应一个频率成分的振幅和相位。
我在给朋友解释FFT时常用一个类比：把一杯咖啡看成由水和咖啡豆萃取物混合而成，FFT就是反推这杯咖啡里各成分的比例。时域信号是什么味道只有喝下去才知道，频域信号直接告诉你成分表。对于脑电来说，FFT的输入是一段电压数据，输出就是每个频率上的能量分布。
实际操作中，需要注意输出数组的下半部分（N/2之后的点）是正频率的镜像，没有额外的物理意义。取频谱时只用前半部分，也就是0到奈奎斯特频率（采样率的一半）这个区间。比如采样率是500赫兹，那么FFT能分析的频率范围就是0到250赫兹，超出这个范围的高频成分会在采样时混叠进低频区间，这也是为什么采集脑电时抗混叠滤波器必不可少。
2.2 功率谱密度（PSD）的两种口径
谱分析里绕不开的概念是PSD，即功率谱密度。两个最容易混淆的版本是幅值谱和功率谱。幅值谱直接看FFT结果的绝对值，单位是微伏；功率谱则是幅值的平方，单位是微伏的平方。两者在图形上形状一致、纵轴刻度不同，但解读时不能混为一谈。
更严谨的医学工程实践里，用的通常是功率谱密度，它额外考虑了频率分辨率的影响。计算时要把每个频点的功率除以频率分辨率（也就是1/信号长度），这样得到的数值才有跨不同参数设置的可比性。如果换了窗长，频率分辨率变了，未归一化的功率值直接比较是会出偏差的，这是一个非常容易掉进去的坑。
我在实际分析中使用的是scipy里signal.welch这个函数，它默认返回的就是功率谱密度的估计值，单位是V^2/Hz，如果输入信号单位是微伏，算出来就是μV^2/Hz。这个函数内部自动处理了窗函数归一化和频点功率的缩放，比自己手写FFT再平方要省心不少。
2.3 参数计算：采样率、频率分辨率、分段长度怎么定
谱分析的参数之间是互相牵扯的，核心关系式是：
频率分辨率Δf = fs / N
其中fs是采样率，N是参与FFT的数据点数。举个例子，采样率500赫兹，取5秒的数据，N=2500，那么频率分辨率就是0.2赫兹，意味着频谱图上相邻两个频点之间隔0.2赫兹。如果只取1秒数据，N=500，Δf变成1赫兹，低频段的细节就看不清了。
另一个关键参数是最大可分析频率，也就是奈奎斯特频率fs/2。采样率500赫兹时，能看到的最高频率是250赫兹，这远远覆盖了脑电的主要频段，gamma频段上限一般也就到45到50赫兹左右。
实际操作中，我一般会让分段长度至少在2到4秒以上，这样低频段的频率分辨率足够细，能分辨0.25到0.5赫兹的差异。如果数据本身很短，比如事件相关电位里只有刺激前500毫秒的基线，这段数据的频谱分辨率就只有2赫兹，delta和theta的分界会变得很模糊。遇到这种情况，一个折中方法是把相邻的多个trial拼在一起算频谱，来弥补单段时长不足的问题。
3. 从周期图到Welch法：为什么不能直接对整段数据做FFT
3.1 周期性假设与频谱泄漏
很多人第一次接触谱分析时有个疑问：数据多长，就直接做多长的FFT不行吗？理论上可以，但实际效果往往很差。根本原因在于，FFT隐含了一个假设：这段信号是周期性重复的。可脑电信号本身就不满足这个条件——段首和段尾的幅值、相位都对不上，强行拼接会人为制造一个跳变边缘，这在频谱上会表现为能量向周围频点扩散，也就是“频谱泄漏”。
频谱泄漏的直观后果是，本来非常集中的alpha峰变得又胖又矮，能量散布到相邻频点上，峰值幅度被低估，频谱图看着脏兮兮的。解决频谱泄漏的经典办法是加窗函数，让数据在边缘处平滑衰减到零，从而消除段间的跳变。常用的窗有Hamming窗、Hanning窗和Blackman窗，它们在主瓣宽度和旁瓣衰减之间各有取舍。
我在处理脑电时默认用Hamming窗，它主瓣宽度适中，旁瓣衰减也有保证，做Welch频谱估计时是scipy.welch的默认选项之一。窗函数加在每个分段上，计算完之后在平均时会做归一化补偿，所以不会造成功率值的系统性低估。
3.2 Welch平均法：用分段平均换稳定估计
直接对整段数据做一次FFT，得到的频谱方差非常大，因为相当于只用了一个样本去估计每个频率上的功率。脑电信号本身不平稳，频谱的随机波动就更明显。Welch法的思路很朴素：把长信号切成若干小段，每段分别加窗做FFT，再把所有段的功率谱平均起来。由于噪声和随机波动会随着平均段数增加而相互抵消，最终得到的频谱曲线会平滑稳定得多。
我在处理一个5分钟的静息态记录时，采样率1000赫兹，把数据切成4秒的段，一共可以分出75段，Welch平均之后频谱曲线的标准误大约只有单段估计的约1/8.7，稳定性提升非常可观。代价是频率分辨率下降了，从整段直接做FFT的约0.003赫兹降到0.25赫兹，但对脑电频段划分来说完全够用。
3.3 重叠率与窗函数的选择逻辑
Welch法的两条关键参数：分段长度和段间重叠率。分段长度决定频率分辨率，重叠率决定参与平均的分段数量。常见设置是重叠50%，也就是每次滑动半个窗长的距离。这样做的好处是，在不缩短分段长度的情况下增加平均次数，让频谱估计更平滑。
有的流程里会用75%甚至90%的重叠率，平均段数更多、曲线更光滑，但相邻分段之间的独立性下降了，每一段的信息重复度很高，实际提升并没有表面看起来那么大。我一般折中选50%到75%，数据量充足的情况下不会刻意追求高重叠率。
窗函数的选择也有讲究。Hamming窗计算速度快、性能均衡，最适合常规脑电频谱估计。如果主要关心的是alpha这样相对平滑的宽频峰，参数差异其实不大；但如果要分辨两个靠得很近的频率成分，需要主瓣更窄的矩形窗或Kaiser窗，不过旁瓣泄漏会变严重。对绝大多数脑电研究场景，Hamming窗是安全和稳妥的选择。
4. 时频分析：给频谱加上时间轴，短时傅里叶变换解析
4.1 为什么静态频谱不够用
静态频谱分析对平稳信号很有效，但脑电本质上非平稳。事件相关实验里，刺激出现后各频段的功率会随时间动态变化：theta在编码阶段短暂增强，alpha在注意加工时被抑制，beta在运动准备中有特定的能量起伏。这些时间上的变化如果全部平均到一个频谱里，动态信息就全部丢失了。
时频分析就是为了同时保留时间和频率两个维度的信息。实现方式有很多种，短时傅里叶变换（STFT）是最直观的一种：把信号切成一个一个短时间窗，对每个窗内的小段数据做FFT，然后把结果按时间顺序排列，就得到一张二维时频图，横轴时间、纵轴频率、颜色深浅表示功率强度。
我在研究视觉工作记忆任务时，最关心的就是刺激保持阶段theta频段有没有持续增强。如果用静态频谱，只能看到整个实验期间的theta平均值，变化细节全被抹平了。换成时频图之后，theta增强的时间范围、峰值位置、以及和被试行为表现的相关性都变得一目了然，这就是时频分析带来的增量信息。
4.2 STFT的两个关键参数和取舍逻辑
做STFT需要设置两个参数：窗函数长度和窗间步长。窗长决定频率分辨率，步长决定时间分辨率，两者之间存在此消彼长的关系。窗越长，频率分辨率越细，但时间定位越模糊；窗越短，时间定位越准，但频率分辨率变差。我必须强调，这是时频分析最核心的取舍，理解了这个，STFT就算学会了。
以500赫兹采样率为例，用250毫秒窗长（125个点）做FFT，频率分辨率是4赫兹，可以清晰区分delta、theta、alpha这些频段的差异，但对时间上快速变化的成分定位能力有限。如果改用100毫秒窗长（50个点），频率分辨率变成10赫兹，delta和theta基本就分不清了，但时间上能辨别的变化幅度降到100毫秒的量级。
实际做脑电时频分析时，我默认用256毫秒或者512毫秒的窗长，步长设为50毫秒左右，这样出来的时频图在频率和时间两个维度上都比较均衡。如果实验任务是快速呈现刺激、需要精确到几十毫秒的时间差异，我才会把窗长缩短，代价是低频段分辨率变差。
4.3 用STFT看事件相关振荡变化
STFT在脑电研究中的一个典型应用是事件相关频谱扰动（ERSP）和事件相关去同步/同步（ERD/ERS）。这两个指标都是拿实验条件时期的时频功率减去基线时期（通常是刺激前500毫秒左右）的时频功率，看能量是增加了（ERS）还是减少了（ERD）。
我在实际操作中会先把每个trial的时频功率算出来，然后对trial维度求平均，再对每个时间-频率点做基线校正。基线校正有两种方式：减法和除法。减法得到的是绝对功率变化，除法得到的是相对百分比变化，后者在跨被试比较时更常用。ERD/ERS通常用百分比来表示，比如alpha频段在刺激后下降30%，就意味着该频段的功率比基线期少了三成。
需要特别提醒的是，基线校正能不能做好，直接决定了时频图的质量。如果基线期的数据本身就包含任务相关的准备活动，基线功率被人为抬高，后面对比时就会出现假性抑制，这在心理生理学文献中已有讨论。我的经验是把基线窗口尽量缩短，放在刺激前350到150毫秒之间，避开早期的准备效应。
5. 实操流程梳理：从原始脑电到可以放进论文的频谱图
5.1 预处理：不能偷懒的一步
频谱分析对数据质量极其敏感，垃圾进垃圾出这句话在时频分析里体现得淋漓尽致。我自己的标准流程是先做带通滤波，频率范围设成0.5到45赫兹或者0.5到50赫兹，把工频干扰和高频肌电噪声先压下去。高通0.5赫兹的目的就是去除基线漂移，如果这里偷懒，低频段的delta功率会被漂移严重污染，结果完全没法用。
滤波之后做坏段剔除和独立成分分析去伪迹。眨眼伪迹主要表现为额叶区低频高幅的成分，在频谱上会污染delta和theta频段；水平眼动主要影响额叶和颞叶的低频段。这些伪迹如果不处理，会在时频图里形成大片假阳性增强区域。我的惯例是先把明显的坏段标记出来删掉，再用ICA跑一遍，把和眼电高度相关、位于额叶的成分剔除。
预处理阶段的另一个常见疏漏是重参考。头皮脑电在没有参考的情况下的绝对功率意义有限，不同参考方式会影响各频段的功率分布。我在项目里统一用全脑平均参考或者双侧乳突平均参考，并在方法部分明确写清楚，保证结果可复现。
5.2 分段、去伪迹和基线校正
预处理完成之后，正式进入谱分析流程。第一步是把连续数据切成epoch，如果是事件相关任务，以刺激出现时刻为0点，往前取500毫秒基线，往后取1000到1500毫秒分析窗口。如果是静息态，直接按4秒一段切分，段与段之间可以留50%重叠，然后把所有段的频谱平均起来。
每段数据在进入频谱计算之前，还要再检查一遍有没有残留的伪迹。我的做法是设定一个幅值阈值，比如超过±100微伏的段直接剔除，因为这种段通常包含剧烈的肌电收缩或电极松动。阈值设多少合理没有统一标准，需要结合自己的数据质量调整，但不管设多少都要在方法里写清楚，保证透明度。
基线校正这一步，如果是静息态就不需要，如果是事件相关的时频分析，必须做。基线窗口选刺激前300到150毫秒，计算基线期的平均功率，然后用除法做归一化，得到相对变化的百分比。这一步骤在时频分析里属于常规操作，但要留意基线窗口的稳健性，最好多试几个窗口看看结果对基线选择敏不敏感。
5.3 参数选择与完整代码示例
我整理了一个可以直接跑通的Python示例，基于mne读取脑电数据，然后完成Welch频谱估计和STFT时频分析。下面的代码使用的是公开的示例数据，核心逻辑可以直接迁移到自己的数据上。
import numpy as np import matplotlib.pyplot as plt from scipy import signal as sig import mne # 读取数据，假设已经有了原始数据data，采样率sfreq sfreq = 500 data, times = mne.io.read_raw_fif('sample_raw.fif').load_data()[:] # 带通滤波 raw = mne.io.read_raw_fif('sample_raw.fif', preload=True) raw.filter(1, 45) # 提取数据为数组形式 epochs_data = raw.get_data() # shape: (n_channels, n_times) # 用Welch法计算静息态功率谱 channel = 0 # 以第一个通道为例 freqs, psd = sig.welch(epochs_data[channel], fs=sfreq, nperseg=2 * sfreq, noverlap=sfreq, window='hamming') # freqs单位Hz，psd单位μV^2/Hz # 绘制功率谱图 plt.figure(figsize=(8, 4)) plt.plot(freqs, psd) plt.xlim([1, 45]) plt.xlabel('Frequency (Hz)') plt.ylabel('Power Spectral Density (μV^2/Hz)') plt.title('Resting-State EEG Spectrum') plt.grid(alpha=0.3) plt.show()需要说明的一点是，上面代码里nperseg=2*sfreq表示分段长度2秒，noverlap=sfreq表示50%重叠。实际项目中，分段时间要根据你的实验设计调整。静息态分析用2秒到4秒的分段都可以，分段越长时间窗内的频率点越密，低频分辨越细。
STFT部分的代码思路也列一下，方便对照理解：
# 短时傅里叶变换参数 window_len = 256 # 采样点，等于512毫秒（500Hz采样率） step = 25 # 步长50毫秒 f, t, Zxx = sig.stft(epochs_data[channel], fs=sfreq, nperseg=window_len, noverlap=window_len - step, window='hamming', boundary=None) # Zxx是复数矩阵，功率谱是幅值的平方 power = np.abs(Zxx) ** 2 # 绘制时频图 plt.figure(figsize=(10, 5)) plt.imshow(power, aspect='auto', origin='lower', extent=[t[0], t[-1], f[0], f[-1]], cmap='jet', vmin=0, vmax=np.percentile(power, 95)) plt.colorbar(label='Power (μV^2)') plt.xlabel('Time (s)') plt.ylabel('Frequency (Hz)') plt.title('Time-Frequency Representation') plt.show()STFT输出里的f和t分别是频率轴和时间轴，Zxx是复数矩阵，每个元素代表对应时间窗内对应频率的复振幅。做后续分析时通常取幅值平方作为功率值，再对trial进行平均和基线校正。
6. 踩坑记录：谱分析中我遇过的典型问题与排查思路
6.1 频谱泄漏导致的谱峰变宽和旁瓣
最早期做频谱分析时，我直接拿整段数据不加窗做了FFT，结果alpha峰非常宽，旁边还多了很多不规则的毛刺。后来才意识到这就是频谱泄漏。解决办法就是给每个分段加窗，我用的是Hamming窗，加窗之后alpha峰明显收窄变高，整个频谱也干净了很多。
如果你在频谱图上看到某频段出现异常的宽峰，先不要怀疑数据有病理问题，检查一下有没有加窗、以及分段长度是否太短。还有一个容易被忽略的点：如果采集时存在基线漂移，即使加了窗，低频段也可能出现一个很高的斜坡状能量，解决途径是增加高通滤波的截止频率。
6.2 工频干扰抬高了50赫兹附近的能量
50赫兹工频干扰是做脑电的人躲不开的老朋友。即使采集系统声称有陷波滤波，数据里还是可能残留明显的50赫兹及其谐波。频谱图上最典型的表现是50Hz处出现一个尖锐的窄峰。如果这个峰的能量过高，处理方式是在预处理阶段额外加一个48到52赫兹的陷波滤波器。
但我要提醒一点：如果实验任务本身就涉及高频gamma频段，陷波滤波可能会损失掉部分有效信号，因为gamma频段恰好和工频有重叠。这种情况下更好的选择是优化采集环境、做好屏蔽和接地，而不是在数据后处理时亡羊补牢。实在不行再使用陷波，并在方法部分如实说明。
6.3 边缘效应和短数据段的频谱误差
STFT在做时频分析时会引入边缘效应，主要原因是数据段两端的数据不完整，窗函数滑动到边缘时有效数据点变少，功率估计会有误差。我的习惯是在STFT之后把时频图两端各裁剪掉半个窗长的时间范围，只保留中间可靠区域，避免边缘的假性效应被误读。
另一个常见场景是事件相关到刺激后效应比较晚，分段时间不够长，导致感兴趣的频段变化还没结束，STFT已经跑到数据末尾了。这种问题要从实验设计阶段规避，把epoch时间窗设置得长一些（比如刺激后1500毫秒），给时频分析留足剩余空间。
6.4 统计检验时频功率时的多重比较问题
时频分析结束后，常规做法是比较不同条件下的功率变化，这就涉及统计检验。时频图上的时间-频率点数量非常多，如果逐点做t检验或者ANOVA，会触发严重的多重比较问题，假阳性率非常高。最常见的解决办法是进行聚类置换检验，把相邻的显著时间-频率点聚成一个簇，再对簇层面做置换检验控制总体错误率。
这一块我在mne-python里用mne.stats.permutation_cluster_test直接实现，比逐点检验稳健很多。另外一个务实的建议是，在做全脑全频段搜索之前，最好基于研究假设预先锁定目标频段和时间窗口，只在假设的区域内做检验，既提升统计效力，也降低解读难度。别试图把所有显著的时频色块都“发现”一遍，这往往是假阳性工厂。
写在最后的几点心得
谱分析这条路，我走了很久才真正理解那些参数不是孤立的，而是彼此牵制的一整套体系。采样率决定了能看到多高的频率，数据时长决定了能分辨多细的频率，窗函数和重叠率决定了频谱的稳定程度，而STFT里的窗长和步长决定了时频图上时间和频率精度的天平怎么偏。做分析前先把这些关系在脑子里过一遍，比拿到数据就埋头跑代码要高效得多。
如果只让我说一个建议，那就是：永远把预处理放在第一位。滤波、去伪迹、坏段剔除这些步骤做得越扎实，后面谱分析的结果就越可靠。很多时频图上说不清道不明的“异常色块”，追根溯源都来自预处理阶段偷懒留下的伪迹。数据准备干净了，频谱和时频分析都只是按部就班的计算而已。
最后再分享一个小技巧：当你不确定某个频段功率变化到底可不可信时，把原始信号画出来，结合你的处理步骤推理一遍。假如基线期有眼球运动残留，那theta频段的增强很可能就是伪迹在捣鬼。把时域信号和频域结果对照着看，是我排查假阳性最顺手的方式，你试过之后大概率也会觉得它很实用。
