先说我最近遇到的一件小事。某次在现场处理一台旋转设备的振动数据，采样率设了50kHz，转速大概3000转/分。我随手用FFT画了幅频谱图，峰值频率看着挺对，等换了另一段时间窗口重算，幅值变了，低频段噪声底也跟着上下跳。起初以为设备状态有波动，折腾半天发现不是设备的问题，根子出在“谱估计方法”这个环节上。
功率谱密度的估计，说白了就是把一段信号从时域转到频域，搞清楚它的能量在不同频率上是怎样分布的。但这件看似简单的事，背后埋着方差、偏差、分辨率、泄漏这一整套博弈。你拿不同方法算，或者同一种方法换了参数，出来的谱可能完全不是一幅图。这篇东西我会从经典周期图法讲到Welch法，再讲到AR模型和MUSIC这类现代方法，最后给出一套可以直接跑的实操流程和调试心得。适合正在做信号分析、故障诊断、音频处理，或者刚入门想搞懂谱估计水有多深的朋友。
1. 先搞清楚谱估计到底在解决什么问题
1.1 一条信号里藏的信息，时域真的看不全
很多朋友一开始接触谱估计，是从傅里叶变换开始的：把信号x(t)做个FFT，然后画出幅度谱。但“功率谱密度”和“幅度谱”是完全不同的两个概念。幅度谱反映的是各个频率分量的幅度大小，而功率谱密度反映的是信号功率在频率轴上的密度分布，它的物理单位常见的是W/Hz，或者工程上喜欢用g²/Hz、V²/Hz这类量纲。
我有段时间做结构振动测试，传感器输出的是一个时域波形，光看波形你只能知道振动大小大概什么样，但分不清这个振动是由哪个部件引起的。轴承外圈故障的振动特征频率、齿轮啮合频率、转子不平衡导致的基频，这些信息全部藏在频域里。功率谱密度估计就是干这个用的：给定有限长的观测数据，我们想尽可能准确地复现出信号在频域上的能量分布。
注意我用的词是“复现”，不是“还原”。真正的功率谱密度定义是针对无限长随机信号做自相关后再傅里叶变换得到的理论值，而我们手里的数据永远有限，所以任何方法说到底都是在做“估计”。理解这一点，你就能接受后面所有方法的妥协和取舍。
1.2 “估计”二字的真实分量
既然是估计，就会有偏差、有方差。这是理解整个谱估计的关键门槛。举个例子，你用一段N点数据做DFT，得到的是真实频谱的某种近似。经典周期图法确实很简单——直接对N点做FFT再取模平方——但它的估计方差不会随着数据长度N增加而像直觉中那样有效收敛。或者说，它的一致性问题很麻烦：数据越多，谱图的细节越丰富，但每个频率点的估计值依然在真实谱附近剧烈波动。
这就导致了第一个经典矛盾：想提高频率分辨率就得用更长的数据，可数据一长，周期图法在单个频率点上的波动并不会因平均而消除，噪声底在视觉上反而更乱了。你用一万点数据算出来的噪声基底，可能比用一千点数据算出来的还要毛糙。
所以，谱估计发展史基本上就是在“分辨率”和“方差稳定”之间找平衡的历程。经典方法靠加窗和分段平均来压低方差，代价是牺牲分辨率；现代方法靠参数模型或子空间分解来突破某些限制，代价是计算复杂度和对先验信息的依赖。没有一种方法是万能的，关键在于你手里的数据长什么样、任务最看重什么。
2. 经典方法：周期图与BT法的真实表现
2.1 周期图法：最直观，但毛病最明显
周期图法的做法，就是把长度为N的观测序列直接做离散傅里叶变换，然后取模平方再除以N。公式上可以写成：
P_xx(k) = (1/N) · |X(k)|²
其中X(k)是N点DFT的结果。我刚开始接触谱分析时，觉得这个办法理所当然：傅里叶变换把信号分解成正弦分量，模平方就是能量，除以N就是平均功率。实测下来，用一两个正弦信号做测试，周期图法的谱峰位置确实很准，但问题藏在细节里。
第一个麻烦是频谱泄漏。如果信号频率不是FFT频率分辨率的整数倍，能量就会“漏”到相邻谱线上去，峰值变矮、附近多出很多旁瓣。第二个麻烦是方差大。这在前面提过，单次周期图的方差大约跟信号功率的平方一个量级，这意味着谱图的随机起伏非常明显，两个频率相近、幅度也接近的分量，很可能淹没在这种波动里。实测对比下来，同样的数据，用周期图法和用Welch法得到的谱图，视觉差异会让你怀疑到底谁算错了。
很多人误以为加个窗能解决泄漏，但加窗只是抑制旁瓣，代价是主瓣变宽，频率分辨率变差。这件事我在第三节会展开。
2.2 维纳-辛钦定理与BT法：换个路径照样有坑
既然直接对数据做FFT不省心，那就绕道走。维纳-辛钦定理说，平稳随机信号的自相关函数和它的功率谱密度是一对傅里叶变换对。基于这个原理，BT法（Blackman-Tukey法）先估计自相关函数，再对自相关函数做傅里叶变换得到谱。
思路听起来很顺，实际操作会有新的问题。自相关函数估计本身就需要有限长度数据，而且延迟点数越大，自相关估计的样本数越少，可靠性越差。所以BT法通常要给自相关函数加一个滞后窗，只保留前面一部分延迟，这其实就是在做一种加权平均。
我早年间做过一段语音信号分析，用BT法跟周期图法对比，BT法能得到相对平滑的谱，但滞后窗取得太窄时谱峰被显著展宽，两个靠得很近的谐波直接糊成一坨。这个体验让我记住一件事：谱估计方法没有银弹，每种平滑手段都在拿分辨率做交易。事实上现在许多工程实践里，Welch法已经取代了BT法成为更常用的经典方法，因为它的参数语义更直观。
2.3 窗函数：一面是泄漏，另一面是分辨率
窗函数这个话题值得单独拎出来说。我们在做DFT时，隐含的操作其实是用一个矩形窗去截断无限长信号。矩形窗的频谱主瓣窄、旁瓣高，能量泄漏严重。换用汉宁窗或汉明窗，旁瓣能压下去二三十个dB，但主瓣宽度会翻倍。这意味着原本可以区分两个间隔为Δf的分量，加窗后可能需要2Δf甚至更宽才能分开，这就是谱估计里常说的“分辨率代价”。
我个人的选择习惯：如果目标是找强信号附近的小分量，用汉宁窗或平顶窗；如果目标是保留幅度精度，平顶窗更好；如果目标仅仅是找谱峰频率位置，汉宁窗的旁瓣抑制够用且主瓣不算太宽。加窗之后，记得谱的幅度要做相应的幅度修正或能量修正，否则你会觉得幅值莫名变低了。这个细节很多人忽略，导致谱图幅值对不上。
3. Welch法：工程实践的事实标准
3.1 分段、加窗、平均：用方差换稳定
Welch法的核心思想很简单：把长度为N的数据先分成若干段，每段做周期图，然后把所有段的周期图取平均。理论依据是，各段数据近似独立，平均之后方差大约降到原来的1/K，K是段数。就这么一个朴素的平均操作，让谱图在视觉上稳定得多，所以它成了几乎所有商用软件里“频谱分析”的默认选项。
但请注意，天下没有白拿的好处。分段意味着每段长度变短，段的DFT分辨率就变成fs/L，L是段长。于是Welch法的实际分辨率比全数据周期图法低，但换来了方差大幅降低。也就是说，Welch法把经典周期图法“高分辨率、高方差”的毛病，扭成了“可调分辨率、低方差”。
我经常跟朋友打比方：周期图法像拿一台高倍望远镜看远方，能看清细节但手一抖画面就糊；Welch法像把镜头换成广角，稳定是稳了，但细节不如高倍镜。关键是你想要哪个。
3.2 关键参数怎么选：nperseg、noverlap与窗函数
在Python里调用scipy.signal.welch大概是做谱估计最省事的路径。函数里有几个参数直接影响结果形态：nperseg指定每段长度，noverlap指定前后两段重叠的样本数，window指定窗函数。
重叠率是个值得说的点。按理说分段越多方差越低，但如果段与段完全不重叠，数据利用率太低。实际中常用50%重叠，此时相邻段部分相关，等效独立段数增加得没那么理想，不过方差改善依然明显。有些场景会用到75%重叠，计算量上去了但谱图更平滑。我实测过50%和75%的差异，对大多数振动和音频信号来说，50%重叠已经能给出稳定结果，75%的边际收益不大。
窗函数的选法配合段长一起考虑：nperseg越小，主瓣越宽，相邻谱峰越容易混在一起；窗函数旁瓣抑制能力又决定小信号会不会被大信号的旁瓣淹没。原则很简单：先根据你关心的最小频率间隔定出所需的段长，再选一个旁瓣足够低的窗。
3.3 频率分辨率的账，得自己先算清楚
很多朋友用Welch法时只关心“图好不好看”，却不验算频率分辨率。频率分辨率通常定义为Δf = fs / L，L就是nperseg。假设采样率是50kHz，你关心两个间隔20Hz的分量，那L至少要2500点，再考虑加窗展宽，实际需要更多的点数。如果你想一个段长既满足分辨率又控制方差，就得看总数据长度够不够。
我踩过这样一个坑：手头数据只有1秒长，采样率50kHz，也就是5万点。当时图省事，nperseg填了8192，分辨率大约6.1Hz，看起来挺不错。但问题是我关心的一个特征频率是12Hz附近，加窗展宽之后，跟直流分量附近的低频噪声糊在一起了。后来把nperseg改成16384，分辨率到3Hz左右，情况才缓解。这就是“先算账再选参”的意义：把分辨率目标换算成最小段长，再去平衡重叠率和窗函数。
4. 现代谱估计：AR模型与MUSIC的取舍
4.1 AR模型谱估计：适合短数据的参数化思路
当数据段很短，经典方法的分辨率跟不上时，参数化方法就登场了。AR模型谱估计的思路是：把信号建模成一个白噪声激励下的全极点滤波器输出，然后估计滤波器系数，再根据系数算出功率谱。
具体估计系数时常用的有Yule-Walker方程和Levinson-Durbin递推，它们在Python里有现成工具，比如statsmodels 或者 spectrum库。AR谱估计最大优势是在数据短的情况下能给出比周期图更平滑、峰值更尖锐的谱。这在雷达、生物医学信号处理里很常见——单次心动周期就那么零点几秒，还要从中提取频谱特征，经典方法很难办。
但AR模型不是拿来就能用。模型阶数p需要选择，选低了谱峰被平滑掉，选高了会出现虚假峰。理论上可以用AIC、BIC这类准则来找合理阶数，实际中我建议先跑几个阶数做对比，看看谱峰数量是否稳定。如果阶数一变动、谱峰就跳来跳去，那说明这个数据可能不适合用AR模型。
4.2 MUSIC算法：子空间方法为什么能“超分辨”
MUSIC算法是子空间类方法的代表，它的基本认识是：信号由若干正弦分量和噪声组成，协方差矩阵的特征向量可以分成信号子空间与噪声子空间。因为信号和噪声子空间正交，我们可以构造一个“空间谱”函数，在信号频率处会出现尖锐峰值。
这方法厉害的地方在于，它的频率分辨率理论上不受数据长度限制，两个频率很接近的正弦也能分开。我第一次在仿真里看到MUSIC把相隔5Hz的两个分量干净利落地分开时，确实被震撼了一下。但震撼归震撼，MUSIC需要你事先知道信号源的个数，这个数错了结果就完全不可信。估计源个数本身也是个难题，工程上常常用特征值分布或者信息论准则去猜。
我做设备故障诊断时试过用MUSIC分析轴承振动信号，发现它擅长分辨相邻谐波，但对宽带噪声的谱形描述并不好。你问它“哪个频率有峰值”它很强，你问它“这个频段的能量有多高”它就有点力不从心。所以它更适合做“谱线检测”，不太适合做“功率密度计量”。
4.3 选型判断：经典和现代，不是谁取代谁
聊到这儿，你应该能感受到，经典方法和现代方法其实是不同目标下的产物。经典方法基于非参数化的傅里叶分析路径，通用性极强，噪声特性可控，适合大多数测量任务和未知信号的分析；现代方法利用模型假设换取更高的分辨率和短数据适应性，但一旦模型假设与真实信号不符，结果会错得离谱。
我的判断标准很实际：信号信噪比不错、数据量充足、需要定量功率谱密度，首推Welch法；数据很短、信噪比高、只关注频率位置，可以上MUSIC；如果数据平稳性存疑，或者你需要平滑谱形又不挑分辨率，AR模型可以考虑。把这些方法当成工具库，而不是互相替代的关系，是做好谱估计的第一步。
5. 一次完整的谱估计实操：参数怎么定
5.1 构造测试信号：先把方法验证明白
纸上谈兵没意思，我习惯先构造一个已知答案的信号来校准方法。比如采样率fs=5000Hz，信号时长2秒，里面放两个频率分量：一个250Hz幅度0.8，一个1250Hz幅度0.4，再加一点点白噪声。这样我们心里清楚真实谱应该在哪两个位置有峰，谱峰值比应该接近4:1。
验证时先跑周期图法，看两个峰位置是否准确。因为信号频率可能不是FFT分辨率的整数倍，你会看到峰值幅度偏低且旁边有旁瓣。这时候做个加窗再对比，就知道泄漏长什么样了。接下来把Welch法接上，用nperseg=1024、noverlap=768、window='hann'跑一遍，两个峰的位置依然准确，而且谱线明显平滑很多。
5.2 Python实现与结果解读
下面这段代码可以直接拿来跑，配合注释理解每一步的作用：
import numpy as np from scipy import signal import matplotlib.pyplot as plt fs = 5000 N = fs * 2 # 2秒钟数据 t = np.arange(N) / fs rng = np.random.default_rng(42) x = 0.8 * np.sin(2 * np.pi * 250 * t) + 0.4 * np.sin(2 * np.pi * 1250 * t) x = x + 0.05 * rng.standard_normal(N) # 周期图法 f1, P1 = signal.periodogram(x, fs, window='hann', nfft=4096) # Welch法 f2, P2 = signal.welch(x, fs, window='hann', nperseg=1024, noverlap=768) fig, ax = plt.subplots(1, 2, figsize=(12, 4)) ax[0].semilogy(f1, P1) ax[0].set_title('Periodogram') ax[1].semilogy(f2, P2) ax[1].set_title('Welch') plt.show()跑完你会看到：周期图法的噪声基底毛糙，Welch法的噪声基底平坦得多。两个方法在峰值频率附近都能找到位置，但Welch的峰值更“干净”。另外注意，由于我加的是汉宁窗，周期图法直接得到的峰值幅度比真实值低一些，这是幅度修正没做导致的，不影响频率定位，但如果你做的是定量幅值测量，就得乘上窗函数的相干增益修正系数。
5.3 实操中的参数调整经验
在真实项目里，我一般按这个顺序来调参数：先看总数据长度够不够，再定nperseg使分辨率满足需求；接着选一个旁瓣性能合适的窗；最后设重叠率50%起步，视谱线平滑度决定是否增大。
如果谱线总体平滑但就是有几个异常的尖刺，优先怀疑是不是窗旁瓣泄漏加上微弱干扰；如果谱图太平滑、连本该有的峰都被压扁了，多半是nperseg太小；如果波形有趋势项或直流偏置，记得先做去均值或高通滤波，否则0Hz附近的巨大能量会把低频段拖得很丑。这个细节我在下一节会专门讲。
有一回我处理加速度传感器信号，频谱图上在10Hz以下出现一个超级大的峰，换了多少个窗都消不掉。后来才想起来传感器安装时带了零点几伏的偏置电压，直接进ADC了。把直流分量去掉之后，低频段立刻正常。所以说，谱估计的很多“怪象”根本不是估计方法的问题，是前端数据质量的问题。
6. 常见坑位与排查心得
6.1 谱线“飘移”与泄漏：先检查窗和段长
现象：同一组数据，不同时间段截出来算Welch，峰值频率在不同帧之间小幅度跳动。这其实不是谱估计方法算错了，而是信号本身可能有轻微的非平稳性，或者信噪比太低，噪声让峰值估计在几个频点之间随机摆动。
排查顺序：先看总谱图，如果峰值底下的噪声底很厚，说明SNR低，考虑增加数据量、提高段长；如果峰值本身占了好几根谱线宽度，那就是主瓣展宽，不是飘移，是分辨率不足。如果换了更长的nperseg之后峰值变尖锐，问题就定位清楚了。
还有一类情况是谱峰位置准确但旁边出现了对称的“裙边”，这强烈提示窗函数旁瓣不够低。换成旁瓣衰减更快的窗（如Blackman-Harris窗）通常能缓解。实测中窗函数从汉宁换成Blackman-Harris，旁瓣能差出好几倍，代价是主瓣又宽了一点，务必平衡。
6.2 负频率、直流分量与去均值
做实信号FFT时，频谱是关于0Hz对称的。很多人画谱图只画正频率部分，这没问题，但如果你用单边谱，记得某些库默认没有对负频率做合并，幅值处理会和双边谱不一致。具体来说，单边谱需要把除0Hz和Nyquist频率外的正频率幅值乘2，才能等效还原真实功率。
直流分量是另一个高频坑。信号里只要有个均值不为零的偏置，谱图在0Hz处就会出现一个巨大峰值，而泄漏效应会拖着低频段一片虚高。处理办法很简单：先减均值。scipy.signal.detrend可以做线性趋势去除，我在处理振动数据时基本必做。去均值之后再看低频，很多“莫名其妙”的谱峰就消失了。
还有一点，重采样或重采样到整数倍频率后，如果数据长度没有对齐，FFT可能引入频率误差。最好先确认fs是不是恒定值，存在重采样误差时优先做时钟校正。
6.3 我踩过几次坑之后的小结
用Welch法做振动谱分析时间长了，我总结出几个原则：第一，永远先验算分辨率，估算目标频率间隔再定段长；第二，画图时同时把周期图和Welch放在一起对比，能帮你快速判断哪些特征是真的谱峰，哪些是估计抖动；第三，参数调完看谱之前，先看看原始时域信号有没有毛刺、断点、溢出，数据质量不过关时，再精细的谱估计都是白搭。
功率谱密度估计是那种看着简单、深入之后发现到处是选择的领域。每一套参数背后都对应一组权衡，没有绝对正确，只有在该场景下合不合适。我个人的体会是，与其追求某个方法的“最强性能”，不如先吃透经典方法在分辨率与方差之间的取舍逻辑。把Welch用熟、用到能预判它的表现，再去尝试AR模型和MUSIC，路会顺很多。最后分享一个小技巧：每次做谱分析前，先用一个已知的合成信号做方法校准。这样当你面对真实数据时，再也不担心是被信号特征骗了，还是被估计方法骗了。
