如果你做过信号处理，大概率遇到过这种让人烦躁的场景：手里数据不算长，可能是一次雷达快拍，也可能是传感器采到的一段冲击响应，或者是一段只有几百个采样点的生物电记录。你想知道里面到底有哪几个频率分量，却怎么折腾都不够清楚——直接上FFT，谱峰又宽又糊，两个靠得很近的频率成分黏在一起，换窗函数、补零、分段平均，折腾一圈下来还是分不开。这时候，就该轮到现代频谱估计的几种方法出场了。
我最早接触这个概念，是因为某次处理一段振动监测数据，采样率不高、记录时间很短，经典谱图上看不出两个相邻模态频率的边界。后来换成现代谱估计方法，用同样长的数据硬是把两个峰给劈开了。那一刻我才意识到，FFT不是唯一的答案，甚至不是最聪明的答案。这篇文章想跟你聊清楚：现代频谱估计到底在解决什么问题，主流的几条技术路线分别是什么，各自的性格脾气怎么样，以及落到工程实践里怎么选、怎么调、怎么避坑。适合正在做雷达、声呐、语音、振动诊断、生物电信号处理的人，也适合刚学完经典数字信号处理、想往更高分辨率方向进阶的同学。
1. 为什么还要搞频谱估计：经典方法到底卡在哪
1.1 分辨率瓶颈：瑞利限与“谱线看不清楚”的真相
经典频谱估计的起点是周期图，说白了就是把一段有限长度数据跟不同的复指数做相关，看哪个频率成分跟数据的“步调”最一致。这个方法简单直观，但有一个硬伤：频率分辨率被观测时间死死卡住。理论上的分辨率极限大约是 Δf ≈ 1/T，T 是观测时长。这不是算法不够努力，而是有限长数据本身携带的信息量有限，两个频率在观测窗口内只差了一小段相位变化，傅里叶核相关运算压根反应不过来。
举个例子你就明白了。采样率 1000 Hz，采 0.1 秒，那么 T=0.1 秒，分辨率极限大约 10 Hz。如果信号里同时有 502 Hz 和 505 Hz 两个分量，间隔只有 3 Hz，比 10 Hz 小得多，FFT 谱图上只会有一个鼓包，分不出里面是几个峰。你说那我补零到 10000 点？补零只是让频谱曲线的采样网格变细了，谱峰本身的物理宽度没变，两个峰该糊在一起还是糊在一起。这就跟拿放大镜看一张模糊的照片一样，放大镜度数再高，照片本身的细节不够，也看不出更多东西。
1.2 加窗与平均：治标不治本的招数
既然谱峰太宽，很多人第一反应是加窗。窗函数能压低旁瓣，但代价是主瓣变宽，两个本来就靠近的峰反而更难分开。汉宁窗、海明窗这些我都试过，降低泄漏的效果确实有，但分辨率并没有变好。还有一招是 Welch 平均，把一段长数据切成若干小段，分别做傅里叶变换再平均，谱线的方差确实小了不少，但每一段更短了，分辨率进一步恶化。这是经典谱估计内部一个绕不过去的跷跷板：方差和分辨率互相拉扯，你不可能两头都占。
现代谱估计的思路完全不同。它不再把数据当成“有限长确定性波形”去硬做相关，而是先把数据看成某个随机过程的一段实现，然后假设一个模型，用数据把模型参数估计出来，最后从模型里把谱抽出来。换句话说，经典方法是在问“这段数据里有哪些频率成分”，现代方法是在问“这个数据最可能由什么样的模型产生，这个模型的谱长什么样”。后者利用先验的结构信息，所以能在更短的数据条件下获得更尖锐的谱峰。当然它也不是魔法，它受克拉美-罗界这类统计极限约束，但确实比 FFT 更懂得“榨干”数据里的信息。
2. 现代频谱估计的几条路线：从模型、子空间到自适应滤波
2.1 参数化模型法：先假设一个模型，再让数据说话
这一类方法的核心是自回归模型，也就是 AR 模型。它假设观测信号 x[n] 可以由自己过去 p 个值的线性组合再加上一个白噪声激励来产生：x[n] = -∑ a_k x[n-k] + w[n]。这个模型的功率谱非常干脆，就是输出白噪声的功率除以一个多项式模平方。模型的极点位置直接对应谱峰，极点越接近单位圆上的某个角度，那个频率处的峰就越尖锐、越窄。
AR 模型谱估计有几种经典算法，性格差异挺大。最早的 Yule-Walker 方法最直白：先估计自相关序列，再解一组托普利茨方程得到 AR 系数。它算法稳定、计算量最小，但短数据下自相关尾部估计不靠谱，分辨率会打折扣。Burg 算法是实践里用得最多的，它不先估自相关，而是直接递推地最小化前向、后向预测误差的平方和，再用 Levinson 递推更新系数。短数据下分辨率非常好，我很多次用它救活了看起来没救的数据。但它有个小脾气：对相位比较敏感，信噪比偏低的时候谱峰位置可能有一点点偏移。还有一个是修正协方差法，同时最小化前后向预测误差且不强迫满足 Levinson 递推，对正弦型信号特别友好，代价是计算量偏大、内存占用高。
AR 方法的优点是计算效率高、物理意义直观，非常适合那种“白噪声背景 + 几根窄带谱线”的场景。但反过来，如果你的信号本来就是宽带的、平坦的谱，AR 模型就不合适了，它天生擅长画尖峰，不擅长画平顶。
2.2 特征分解与子空间方法：把噪声“分离”出去
子空间类方法走的是另一条路，它不假设信号是某个低阶 AR 过程，而是把观测向量构造成一个协方差矩阵，然后对这个矩阵做特征分解。核心思想很漂亮：如果信号由 m 个复指数成分叠加噪声组成，那么协方差矩阵的大特征值对应的特征向量张成“信号子空间”，其余小特征值对应的特征向量张成“噪声子空间”，而这两个子空间是正交的。于是你可以构造一个谱函数，在每个候选频率处用信号的导引向量去碰噪声子空间，碰得越“正交”，分母越接近零，谱峰就越尖锐。
MUSIC 算法就是这一类里的代表。它把上面的思想直接做成了谱扫描，得到的是一个“伪谱”——谱峰的高度并不代表那个频率的实际功率，但峰的位置极其精确，尤其是在信噪比不错、数据足够的情况下。它的分辨率远超经典方法，甚至可以突破到接近克拉美-罗界的水平。ESPRIT 算法则更进一步，干脆不做谱扫描，它利用均匀采样或者均匀阵列的旋转不变性，直接从子空间里解出频率，计算量更低，还不用逐点扫描。root-MUSIC 是把 MUSIC 的谱峰搜索问题变成多项式求根，精度高、速度快，是我自己比较偏爱的一种变体。
不过子空间方法有一个绕不开的前提：你必须提前知道信号源数 m。源数估错，谱图立刻给你脸色看，要么多出一堆假峰，要么真实峰直接消失。处理相干源（比如多径、同频干扰）的时候，标准 MUSIC 还会失效，因为相干信号会把协方差矩阵的秩“压扁”，信号子空间和噪声子空间不再完全正交。实际里需要做空间平滑或者前后向平滑来“恢复秩”，这属于子空间方法的常规组装套餐。
2.3 Capon 最小方差法：介于经典与子空间之间的稳健派
Capon 方法，也叫最小方差无失真响应谱估计，最初是从自适应阵列信号处理里搬过来的思想。它不假设信号由少数几个离散频率成分构成，而是逐个频率点做自适应滤波：在保持该频率方向增益刚好为 1 的前提下，让滤波器输出功率最小。这样做的效果是，能量会被自动引导到“避开”其他频率成分的干扰，最终得到的谱比经典谱尖锐，又没有 MUSIC 那么娇气。
Capon 谱估计的唯一公式非常简洁：P(f) = 1 / (a^H(f) R^{-1} a(f))，其中 a(f) 是频率导引向量，R 是数据的协方差矩阵。从公式就能看出来，它需要求协方差矩阵的逆，所以协方差估计的质量直接决定谱图质量。快照太少，R 病态，谱就发毛；快照足够，它在低信噪比下比 MUSIC 稳健得多。如果你手里的数据信噪比一般，又不能用太多阶数去凑 AR 模型，Capon 是一个很不错的中间选择。
3. 算法落地：公式、参数选择与参考实现要点
3.1 核心公式速览：每个方法到底在算什么东西
把前面提到的三类方法的核心公式放在一起对比，你就能更清楚地看到它们的共同点和差异。
AR 模型功率谱：
P(f) = σ² / |1 + Σ_{k=1}^{p} a_k e^{-j2πfk}|²
这里 σ² 是白噪声激励的方差，a_k 是估计出来的 AR 系数，p 是模型阶数。这个公式说明：分母越小，谱值越大，所以 AR 谱天然擅长画尖峰。
MUSIC 伪谱：
P(f) = 1 / (a^H(f) E_N E_N^H a(f))
E_N 是噪声子空间的特征向量矩阵，a(f) 是频率 f 处的导引向量。当 f 落在真实频率上时，a(f) 与噪声子空间正交，分母近乎零，谱值飙高。
Capon 谱：
P(f) = 1 / (a^H(f) R^{-1} a(f))
R 是数据协方差矩阵，这个公式体现了在约束增益为 1 的前提下最小化输出功率。每个方法的谱峰位置解释是一致的：分母越小、越接近零，这个频率越像信号的真实分量。理解了这一点，你就抓住了现代谱估计的核心手感——谱峰位置就是答案，谱峰高度则不一定有意义。
3.2 阶数和源数怎么定：AIC/MDL 准则与工程扫描习惯
参数选择是现代谱估计里最容易翻车的环节。AR 方法的阶数 p 如果太小，谱被抹得太平，该有的峰出不来；如果太大，谱上开始长毛刺，全是假峰。我个人的习惯是先做一个阶数扫描：从 p=1 到 p=N/3 左右，把每个阶数下的谱峰位置记下来，看主峰在哪个阶数区间稳定不动。稳定区间里的最低阶数，通常就是比较稳妥的选择。
信息论准则可以做更自动化的判断。对 AR 阶数选择，常见的有 AIC：AIC(p) = N ln(σ_p²) + 2p，以及 MDL：MDL(p) = N ln(σ_p²) + 0.5 p ln(N)。N 是数据长度，σ_p² 是 p 阶模型对应的预测残差方差。AIC 偏向稍大的阶数，MDL 更保守一些。工程实践里我更推荐“扫描 + 稳定性判据”的组合，因为信息论准则给的是一个“建议值”，真正到具体信号上，谱峰稳定了没有、有没有假峰，只有扫过之后才心里有底。
至于子空间方法需要的源数 m，可以从协方差矩阵的特征值入手。特征值按从大到小排列，前面几个明显偏大的对应信号，后面一串平缓的就对应噪声。判断“明显偏大”可以用门槛，也可以看特征值对数图的拐点。更正规的做法是套用 AIC/MDL 类准则，核心思想都是“似然残差 + 参数数量惩罚项”，取最小值对应的源数。这里给你一个实操提醒：处理实值单通道信号时，一个实正弦在协方差里会占用两个共轭对称的分量，所以源数要按实际正弦数量的两倍来给，否则 MUSIC 会找不到峰或者多出峰。
3.3 一段可以跑的参考实现：自己写一个 MUSIC 谱
理论说多了容易飘，直接上一段可以对照实现的 Python 示意代码。这个例子用滑窗构造数据矩阵，再对协方差矩阵做特征分解，最后扫描频点得到 MUSIC 谱。
import numpy as np def music_spectrum(x, n_sources, m=None, nfft=1024, fs=1.0): N = len(x) if m is None: m = N // 2 # 用滑窗构造观测向量矩阵 X = np.array([x[i:i+m] for i in range(N - m)]) # 协方差矩阵 R = X.conj().T @ X / (N - m) # 特征分解，默认按特征值升序排列 evals, evecs = np.linalg.eigh(R) # 取最小的 (m - n_sources) 个特征向量作为噪声子空间 En = evecs[:, :m - n_sources] freqs = np.linspace(0.0, fs / 2.0, nfft) spec = np.zeros(nfft) for i, f in enumerate(freqs): a = np.exp(-2j * np.pi * f * np.arange(m) / fs) spec[i] = 1.0 / (a.conj() @ En @ En.conj().T @ a).real return freqs, spec有几个点必须提醒你。第一，m 取多少很关键，我建议取 N/3 到 N/2 之间。m 太小，子空间维度不够，分辨率上不去；m 太大，滑窗快照数量 N-m 太少，协方差矩阵估不准，谱会发毛。第二，实信号先做希尔伯特变换转成解析信号再跑 MUSIC，或者源数按双倍填，不然很容易得到一张莫名其妙的谱图。第三，这个代码是教学向的时间平滑版本，对付单次记录够用；要是阵列数据，直接把 X 换成各阵元的快照即可。
如果你手头有更顺手的现成工具，MATLAB 里的 pburg、pmusic、peig、capon 这些函数都能直接调，很多细节已经给你处理好了。重点不是背函数名，而是理解每个参数背后对应的物理量——阶数、源数、窗长、快照数，这些才是决定谱图质量的真正变量。
3.4 别忘了前处理：现代方法对数据卫生更挑剔
经典傅里叶方法皮实，一段带趋势、带直流偏置、甚至带点非平稳漂移的数据，照样能画出个大概。现代谱估计就不行了，AR 模型和子空间方法对模型失配非常敏感。数据里的线性趋势会被当成一个极低频的大能量成分，直流分量会在零频附近造成强干扰，带外噪声会把协方差矩阵污染掉。所以我的标准流程是：先去均值、去趋势，再做带通滤波把关心的频段单独切出来，最后才交给现代谱估计处理。这一步看着不起眼，往往能直接决定一个方法“好用”还是“翻车”。
4. 实战对比：几种方法在典型场景下怎么选
4.1 横向对比表：先建立整体印象
把前面提到的几种方法放在一张表里看，各自的性格会更清楚。
| 方法 | 分辨率 | 计算量 | 短数据表现 | 主要短板 | 适合场景 |
|---|---|---|---|---|---|
| Yule-Walker | 中等 | 低 | 一般 | 短数据分辨率受限 | 实时粗估、嵌入式 |
| Burg | 高 | 低 | 好 | 低信噪比时谱峰可能偏移 | 短数据窄带分析 |
| 修正协方差 | 高 | 偏高 | 好 | 计算量偏大 | 正弦/窄带混合信号 |
| Capon | 较高 | 中高 | 较好 | 需要足够快照求逆 | 低信噪比稳健谱图 |
| MUSIC | 极高 | 高 | 不错 | 需已知源数、怕相干源 | 高信噪比精细频率定位 |
| ESPRIT | 极高 | 中 | 不错 | 需均匀阵列或均匀采样 | 阵列测向、频率跟踪 |
4.2 不同场景的选型建议
如果你手里只有一次短记录，比如 128 点或者 256 点，主要目标是看有几个主峰，Burg 是最稳的起手式，计算量小、谱峰尖锐，基本不会给你出大幺蛾子。想进一步把峰的位置定得更准，在信噪比不错的前提下切到 MUSIC 或 root-MUSIC，效果会更好。
如果信噪比不高，MUSIC 的谱会开始发飘，这时候 Capon 更值得信任。我在某次实验数据里对比过：同样一段叠加了强噪声的信号，MUSIC 谱出现了明显的假峰，而 Capon 虽然峰宽了一点，但位置和数量都稳。反过来，如果两个频率间隔特别小，Capon 的分辨率又不够用了，还是得靠 MUSIC 这类子空间方法把“缝隙”撬开。
如果是实时嵌入式系统，算力有限，AR 模型几乎是一眼选中。Levinson 递推的复杂度是 O(p²)，p 通常几十阶，在微控制器上也能跑。AR 系数算出来之后，你甚至不用算整个谱，只选几个关心的频点计算 P(f) 就行，省得整段频谱都扫一遍。至于 ESPRIT，它在均匀阵列测向里是真香，不需要逐点扫描，直接出精度很高的频率或角度估计，非常适合那种连续跟踪多个目标的应用。
还有个经验是：别把宝押在一个方法上。我习惯先用经典谱，比如 Welch，快速看一遍全貌和大致频带；再用 Burg 或 MUSIC 去精修局部峰。两轮结果对照着看，某个峰在两个方法里都出现、位置还几乎一致，那基本可以放心；如果只有一种方法出峰、另一种没有，那就要警惕它是模型失配带来的赝品。
5. 常见问题与排查技巧实录
5.1 谱峰分裂与伪峰：最常见的翻车现场
现代谱估计最让人头疼的就是谱峰分裂和伪峰。谱峰分裂，一个本该是单峰的地方被劈成两三个挨着的尖峰；伪峰，谱图里无中生有冒出一堆看起来很“像样”的峰。这两种情况大多逃不开几个原因：AR 阶数选过大、MUSIC 的源数给多、数据里没去掉趋势项、或者数据本身含非平稳成分。
排查的时候我建议按顺序来。先看数据前处理有没有做到位，去均值去趋势有没有漏掉；再做阶数或源数的扫描，画一张“参数 vs 主峰位置”的图，看峰的稳定性；最后把数据切短再试一次，如果谱峰位置随数据段变动很大，说明信号本身或模型假设可能有更根本的问题。伪峰的问题多出在“模型阶数贪心”上，记住我的原则：能解释峰的最小参数模型，才是好模型。
5.2 谱峰位置偏移：看着像，其实不准
谱峰是出来了，但位置比真实频率差了几个赫兹，这种情况比伪峰更隐蔽。Burg 算法在低信噪比下有轻微的系统性偏差，MUSIC 在源数估计错误时会对峰位置造成拉扯，实信号直接跑 MUSIC 但没有把源数翻倍时，也会因为模型不匹配导致峰位置移动。
解决办法一个是把数据“复数化”，也就是用希尔伯特变换构造解析信号，让实信号的正负频率不再互相干扰，子空间方法的表现会稳定不少。另一个是多快照平均，如果实验条件允许重复采集，多组数据加噪协方差平均能显著压低偏差。干活的时候我还会顺手做个正弦拟合校验：估出来的频率作为初始值，用最小二乘法对原始信号拟合正弦，对比拟合残差。残差小说明这个频率靠得住，残差突然变大就说明谱峰位置可能被带偏了。
5.3 相干源面前，MUSIC 也会“宕机”
相干源指的是多个信号成分之间存在完全的相位关系，比如多径传播到达同一个接收阵列，或者同频干扰和真实信号完全相关。协方差矩阵的秩会因此降低，信号子空间和噪声子空间的维度关系被破坏，MUSIC 谱峰直接消失。这是子空间方法的通病，不是你代码写错了。
基础解法是空间平滑：把阵列或数据窗切成若干重叠的子阵，对每个子阵求协方差再取平均，能把相关源的秩“恢复”一部分。单通道时间序列里也可以用类似的滑动窗平滑，只是效果不如阵列里那么彻底。如果平滑之后还是不行，干脆换 Capon 或者 AR 类方法，它们对付相关源不像子空间方法那么脆弱。实践中我的态度是：数据条件不够好就别硬上 MUSIC，方法选型是对数据复杂度的尊重。
5.4 我的几条实操心得
这几年在现代谱估计上踩过的坑不少，也慢慢攒了一些自己的规矩。第一，拿到任何数据先做“快速检查”——长度多少、有没有趋势、有没有明显非平稳、目标频段在哪。这几件事花不了五分钟，但能帮你避开后面一整天的玄学调参。第二，永远对谱峰高度保持怀疑态度，MUSIC 谱和 AR 谱的峰高不等于功率，想对比不同频率成分的强弱，老老实实回频域做功率积分。第三，多方法做交叉验证，这不丢人，反而是最靠谱的质量控制。我通常会拿 Welch 看轮廓，Burg 或 Capon 看细节，MUSIC 定位置，三个结果对得上才写进报告里。
最后再分享一个实用小技巧：现代谱估计得到峰值频率之后，把峰值频率带代回原始数据做一次最小二乘正弦拟合，用拟合的残差和幅度进一步验证。比如你怀疑某条谱峰是假的，拟合一下就会发现残差明显偏大。这个方法补上了现代谱估计最大的一个短板——它本身不给出可靠的幅度和可信度。稍加延伸，还可以把多个候选频率同时放进一个多正弦模型里做联合拟合，能进一步确认整体频率集合是否自洽。这个习惯帮我排掉过不少假峰，比调参数管用得多。
