简介：这是一份信号处理领域的现代谱估计实验报告，面向通信、声学、电子工程等方向的学习者与工程技术人员，解决经典谱估计在低信噪比下频率分辨率不足、方差性能欠佳的问题。包内含1个doc文档，压缩包约154KB，系统梳理了参数模型法、AR模型、Levinson-Durbin递推算法，并结合自相关法、Burg法、协方差法、改进协方差法开展对比实验。文档记录了信号生成与功率谱估计的完整编程步骤，给出了不同信噪比和阶次下的功率谱图、结果分析与结论，可帮助读者直观理解各方法在噪声抑制、频率分辨率和稳定性方面的差异，适合作为课程实验、期末复习或算法选型时的参考资料。已有153人学习下载。

1. 噪声中正弦信号的现代法频谱分析：在短窗与低信噪比下把频率抠出来

做振动诊断、结构健康监测或者无线电测向的人，大概率都遇到过同一个窘境：一段数据里明明有正弦分量，拿FFT画出来却是一片鼓包，或者峰值断层皮，频率读不准。经典傅里叶谱把分辨率压在“数据长度”这堵墙上，一旦信噪比掉到0 dB附近，周期图上的谱峰就开始和旁瓣糊成一体。所谓噪声中正弦信号的现代法频谱分析，指的是用随机信号的参数模型或特征分解手段——常见的是AR模型、Capon最小方差法、MUSIC和ESPRIT——替代傅里叶周期图来估计这些埋在白噪声里的正弦频率。它最大的价值在于：不依赖长数据窗，频率相隔几个赫兹也能分开，并且在低信噪比下仍然给出可读的谱峰。适合的是手里只有几十毫秒到几秒数据、却要做高分辨频率估计的工程场景。

2. 现代法频谱分析的三条路线：AR模型、Capon与子空间方法的选型依据

2.1 经典傅里叶方法失效在哪

在做现代法之前，先把“现代”这两个字的落点说清楚。FFT类的经典法本质上是把有限长观测数据窗与无限长正弦信号相乘，再假设窗外是周期延拓。这引入了两个先天缺陷：一是矩形窗的频谱旁瓣会掩盖邻近小信号，二是频率分辨率受限于瑞利判据，1/T的数据长度摆在那里，两个频率差小于1/T时，周期图的峰值必然连成一个包。就算加窗、补零、加平均，也只是平滑了方差，并不能突破分辨率物理极限。

噪声中的正弦信号，恰好踩在经典法的痛点上。比如一段0.5秒的数据，100 Hz和105 Hz这两个分量，FFT理论分辨极限是2 Hz，能勉强分清，但实际信噪比一低，噪声方差摊在每一个频点上，峰值定位的抖动可能达到好几个赫兹。现代法频谱分析换了一条路：不再把数据当成一段截断的确定信号，而是建立一个随机信号模型，比如“白噪声激励一个全极点系统”的AR模型，或者“K个正弦加噪声”的子空间模型，用模型参数去反推正玄频率。模型的自由度比数据样本数小得多，相当于用先验结构换来了分辨力。

2.2 主流现代法频谱分析算法与适用场景比较

把市面上真正能落地的现代法频谱分析工具排一下，主要就是三类。

AR模型适合“我只知道数据平稳，不清楚里面有几个正弦”的情况，因为谱自然连续，读出峰个数就行。Capon的定位更像“频率跟踪器”，适合慢变非平稳信号，但它得到的谱叫做伪谱，峰值不代表功率。MUSIC是线谱估计的杀手锏，两个频率差只有FFT极限的几分之一也能分开，但代价是要先估计信源个数，这个前置参数一旦猜错，谱就是乱码。

2.3 选型口诀与最少动手准备

我的习惯是：先拿AR模型做探索，因为它的输入只有一个阶数p，谱画出来直观；如果AR谱已经有分离峰，就满足需求，不用上MUSIC。如果频率靠得太近，AR谱裂不开，再切到MUSIC，把p换成信源个数。Capon放在中间一层，用来交叉验证AR峰的可信度——如果AR和Capon在同一个频率都出现峰值，基本可以断定不是伪峰。

动手前只需要两个工具：一段采样率已知的离散信号、一个NumPy。所有算法在百量级数据点上就能跑，不需要GPU，不需要专业软件。下面两章分别把AR模型和MUSIC的完整实现写出来，代码可以直接复制到本地跑通。

3. 用AR模型做现代法频谱分析：协方差法与Burg的完整Python实现

3.1 从自回归模型到频谱：核心推导

AR模型假设观测序列x[n]满足

x[n] = -∑ a[k] x[n-k] + e[n]

其中e[n]是均值为零、方差为σ²的白噪声。把这个方程看成一个全极点滤波器，输入是白噪声，输出就是x[n]，那么输出频谱为

P(f) = σ² / |1 + ∑ a[k] exp(-j2πfk)|²

所以只要从数据里估计出系数a和激励方差σ²，然后对f逐点算分母，就能得到平滑连续的现代法频谱分析结果。这里有一个容易误解的点：AR谱不是对数据直接做变换，而是对“模型”做变换，谱的平滑度由模型阶数决定，和数据点数不是一回事。

估计a的办法有两类。一类是解Yule-Walker方程，前提是假设x是平稳随机过程，用自相关函数的Toeplitz结构求解；另一类是直接从最小二乘残差出发，把每个样本用前面p个样本回归，这是协方差法。Burg算法则是折中：它用前向和后向预测误差的均方和最小来递推反射系数，对短数据段的谱估计效果通常比Yule-Walker好。

3.2 协方差法的最小二乘实现

协方差法的思路最直观：把x[p:]中的每个点当作y，把前p个点当作特征，做一次普通的最小二乘回归。

import numpy as np def ar_covariance(x, order): """用协方差法估计AR系数，返回包含首项1的系数数组""" n = len(x) # 构造设计矩阵：每行是某个样本点之前的order个样本，逆序排列 X = np.empty((n - order, order)) for k in range(order): X[:, k] = x[order - 1 - k: n - 1 - k] y = x[order:] # 最小二乘解：a = argmin || Xa - y ||^2 coef = np.linalg.lstsq(X, y, rcond=None)[0] # 按AR标准形式补上1，注意正负号 return np.r_[1.0, -coef] def ar_spectrum(a, sigma2, fs, nfft=2048): """根据AR系数和激励方差计算功率谱密度，频率轴从0到fs/2""" freqs = np.linspace(0, fs / 2, nfft) w = 2 * np.pi * freqs / fs exps = np.exp(-1j * np.outer(w, np.arange(len(a)))) denom = np.abs(exps @ a) ** 2 psd = sigma2 / denom return freqs, psd # 示例：0.5秒数据，其中包含100Hz和105Hz正弦，叠加白噪声 fs = 1024 t = np.arange(0, 0.5, 1/fs) x = np.sin(2*np.pi*100*t) + 0.8*np.sin(2*np.pi*105*t + 0.3) x += 0.5 * np.random.randn(len(t)) # 信噪比大约10dB order = 8 a = ar_covariance(x, order) # 激励方差由残差估计：用拟合误差的均方值近似 X_pred = np.empty((len(x) - order, order)) for k in range(order): X_pred[:, k] = x[order - 1 - k: len(x) - 1 - k] residual = x[order:] - X_pred @ a[1:] sigma2 = np.mean(residual**2) freqs, psd = ar_spectrum(a, sigma2, fs, 4096) # 打印谱峰所在频率 peak_freq = freqs[np.argmax(psd)] print(f"主峰频率: {peak_freq:.2f} Hz")

这段代码里，设计矩阵X的构造方式是关键：第k列取的是每个样本点前k+1个样本，行数等于样本数减阶数。np.linalg.lstsq用最小二乘求解回归系数，得到的系数以“x[n] = -∑a[k]x[n-k]”的标准形式存储，所以在返回时取了负号。ar_spectrum里用exp(-1jwn)构造频率响应，分母越小，谱值越大，谱峰位置对应系统极点。

参数说明：order是AR模型阶数，理论上要大于等于正弦数的两倍，实际工程经验是取数据点数的三分之一到一半，但不要超过数据点数的一半，否则矩阵病态。nfft只是频率网格的密度，不影响谱的分辨率。这里order取了8，恰好能容纳两个正弦的四个极点。

3.3 Burg算法的实现与参数选择

Burg算法的实现稍微绕一点，它通过递推估计反射系数，再转换为AR系数。好处是前向和后向预测误差一起最小化，短数据段的稳健性比协方差法好。

def ar_burg(x, order): """Burg法估计AR系数，返回(a, sigma2)""" x = x.astype(float) n = len(x) # 初始化前向误差ef和后向误差eb为原始信号 ef = x.copy() eb = x.copy() a = np.array([1.0]) for m in range(order): # 当前阶的反射系数：最小化前向和后向误差均方和 ef_next = ef[1:] eb_prev = eb[:-1] num = -2.0 * np.sum(ef_next * eb_prev) den = np.sum(ef_next**2) + np.sum(eb_prev**2) k = num / den # 反射系数，范围在-1到1之间 # 更新格型滤波器的前向、后向误差 new_ef = ef_next + k * eb_prev new_eb = eb_prev + k * ef_next ef = new_ef eb = new_eb # Levinson递推：由反射系数更新AR系数 a = np.concatenate([a, [0]]) a = a + k * a[::-1] # 激励方差由最终的前向误差均方值近似 sigma2 = np.mean(ef**2) return a, sigma2 a_burg, sigma2_burg = ar_burg(x, 8) freqs_burg, psd_burg = ar_spectrum(a_burg, sigma2_burg, fs, 4096) peak_freq_burg = freqs_burg[np.argmax(psd_burg)] print(f"Burg主峰频率: {peak_freq_burg:.2f} Hz")

Burg迭代里每一次m更新：先算出反射系数k，然后更新误差序列，最后用Levinson递推把AR系数累加进来。这段代码里的a + k * a[::-1]是Levinson递推在实信号下的简洁写法，它保证每一阶都满足反射系数绝对值小于1，因此AR模型稳定性天然有保证。sigma2取最终前向误差的均方，相当于模型无法解释的部分，低频段的基线就靠它托底。

3.4 频轴与峰值读数

AR谱画出来之后，读数步骤有一点容易错。推荐的做法是：直接对psd做峰值检测，找局部极大值而不是全局最大值，因为工程里常常有多个正弦。峰值检测时建议对psd做10 * log10(psd)“Log化”，这样让小峰更容易被峰检测器识别；如果两个相邻峰距离在2个频率网格以内，把nfft调到8192再画一次，不要直接认为模型失效。

AR谱的纵轴有真实物理意义，单位是W/Hz或者V²/Hz，半功率带宽对应极点离单位圆的距离，可以粗略换算成阻尼比。这一点比MUSIC伪谱友好得多，所以我总是先看AR谱，而不是直接上子空间法。

4. 用MUSIC做线谱估计：特征分解与频率搜索的实操

4.1 MUSIC原理与数据矩阵构造

MUSIC把观测信号建模成K个正弦加白噪声，即

x[n] = ∑ A_i exp(j2πf_i n) + w[n]

构造一个长度为L的滑窗，把数据截成多段，每段组成一个列向量，它们的协方差矩阵R = E[s sᴴ]可以做特征分解。理想情况下，R的前K个大特征值对应信号子空间，剩下的L-K个小特征值对应噪声子空间。因为信号子空间和噪声子空间正交，所以扫描频率f时，构造导向矢量a(f)，如果a(f)与噪声子空间近似正交，那么a(f)ᴴ E_noise E_noiseᴴ a(f)接近零，其倒数会出现尖峰。

工程实现的关键在R的估计。数据段数太少，R的秩就凑不齐；段数太多，每段长度L太短，频率分辨力又不够。一般取L在数据点数的1/4到1/3之间，重叠率50%来切段。

4.2 Python实现与关键参数

def music_frequency(x, fs, num_sources, L=None, nfft=8192, min_freq=0, max_freq=None): """MUSIC谱估计，返回频率轴和功率谱""" n = len(x) if max_freq is None: max_freq = fs / 2 if L is None: L = n // 3 # 重叠50%切段，构造自相关矩阵 seg_stride = L // 2 segments = [] for start in range(0, n - L + 1, seg_stride): segments.append(x[start:start + L]) R = np.zeros((L, L), dtype=complex) for seg in segments: s = np.asarray(seg, dtype=complex) R += np.outer(s, s.conj()) R /= len(segments) # 特征分解：特征值从大到小排列 evals, evecs = np.linalg.eigh(R) order = np.argsort(evals)[::-1] evecs = evecs[:, order] # 取噪声子空间：去掉前num_sources个特征向量 noise_eigen = evecs[:, num_sources:] # 扫描频率，计算MUSIC谱 freqs = np.linspace(min_freq, max_freq, nfft) music_psd = np.zeros(nfft) for i, f in enumerate(freqs): # 导向矢量：按实际采样率换算相位增量 steering = np.exp(-2j * np.pi * f * np.arange(L) / fs) # 与噪声子空间做投影，取倒数形成谱峰 projection = steering @ noise_eigen @ noise_eigen.conj().T @ steering.conj() music_psd[i] = 1.0 / np.abs(projection) return freqs, music_psd # 使用示例：两个相距更近的正弦 x2 = np.sin(2*np.pi*100.5*t) + np.sin(2*np.pi*102.0*t + 0.8) x2 += 0.8 * np.random.randn(len(t)) freqs_m, psd_m = music_frequency(x2, fs, num_sources=2, L=180, nfft=8192) # 取出前三个峰 threshold = 0.3 * psd_m.max() peak_idx = np.where((psd_m[1:-1] > psd_m[:-2]) & (psd_m[1:-1] > psd_m[2:]) & (psd_m[1:-1] > threshold))[0] + 1 peaks = freqs_m[peak_idx][:3] print("MUSIC检测频率:", peaks)

这段代码里的R构造用的是“快拍平均”思想，重叠率50%是折中的选择，即保证段数量够，又不让相邻段过度相关。noise_eigen取的是除了前num_sources个特征向量之外的剩余部分；如果num_sources估大，噪声子空间里会被掺入信号成分，伪峰会冒出来。导向矢量用exp(-2jπf n/fs)而不是exp(2jπf n/fs)，两个公式互为共轭，投影的模值不变，不会影响结果，但保持统一能避免调试时晕头转向。

从代码里可以看到，MUSIC谱本质上是“伪谱”，纵轴是导向矢量到噪声子空间距离的倒数，不是真实功率，所以不能用峰高比较信号强度。如果你想估计每个正弦的真实幅度，需要在MUSIC找到频率之后，再用最小二乘拟合幅度和相位，一步到位。

4.3 频率搜索的网格与误差控制

MUSIC的频率搜索步长直接决定定位精度和计算量的平衡。nfft=8192在0到512 Hz范围内对应约0.0625 Hz的网格，够大多数工程需求。但注意，MUSIC的极限分辨力不依赖网格，而是依赖数据长度和信噪比，理论上能达到Cramer-Rao下界；网格太小只是让峰值位置更细，不会提升真实分辨力。

如果峰值一直出现在频率轴的尽头，要检查max_freq设置是否越过奈奎斯特频率。另一个常见问题是在低信噪比下，num_sources设置成2时，真实峰出现在第1和第200个网格，但中间多出很多毛刺；解决方法是先做一个AR谱看峰数，把MUSIC的num_sources设为AR谱可分辨峰数再加0到1。毛刺还可以用谱平滑来压制：对music_psd做滑动平均滤波，窗口3个网格就够，太大反而把相邻峰抹平。

5. 现代法频谱分析的避坑清单：阶数与伪峰问题的真实记录

5.1 阶数选不对，谱峰裂成两片

现象：同一个信号，AR模型阶数取8时一个峰清晰利落，取16时同一个峰变成双峰，频率读数偏移好几赫兹。

原因：AR模型阶数过高，等于给白噪声也建模了。噪声被当成确定性分量，模型把每一个噪声样本都用一个极点去拟合，结果就是谱峰分裂和虚假极点出现。这在低信噪比时尤其严重，因为噪声的高频成分会诱导极点分布变散。

解决：限定order不超过数据点数的一半，这是硬上限。更精细的做法是用信息准则，比如AIC或BIC，order从1扫到min(50, n//3)，取使准则最小的阶数。工程上我喜欢再保守一点：用BIC选出来的阶数减1到2，因为BIC在样本量不大时仍然偏乐观，减几阶能让谱更平滑。

5.2 MUSIC伪峰的四重来源

现象：MUSIC谱里出现好几个峰，但没有一个对应真实正弦频率；或者真实频率处反而凹陷。

原因：最常见的四个原因分别是——num_sources估计过大导致信号特征向量混入噪声子空间；数据段L太短导致协方差矩阵的秩不足，噪声特征值分布不均；信号相关性太强，两个正弦幅度差别很大时，小幅值正弦的特征向量接近正交于信号子空间，解不出来；谱扫描时nfft过高，随机毛刺的峰被放大显示。

解决：先用AR谱估计峰数；L取数据长度的1/4到1/2并至少保证L > 2*num_sources；两个正弦的幅度差超过20 dB时，建议把数据做预白化处理再进MUSIC；扫描频率范围收缩到目标频带，不要全频带瞎扫。伪峰还有一个隐蔽来源是频率恰好落在某一个特征向量的零点上，这在信号带有直流和工频干扰时经常出现，解决办法是先把数据减均值、做带通滤波。

5.3 频率归一化与Nyquist边界

现象：写代码时直接用freqs = np.linspace(0, 1, nfft)扫描归一化频率，乍一看谱峰在0.38处，但换算成Hz是多少分不清；或者代码里fs写错成1000，实际采样率是1024，整个谱轴全部偏移4%。

原因：MUSIC谱的导向矢量用的是相位增量2πf/fs，这里f的单位是Hz；AR谱的w单位是rad/sample归一化角频率。两个混在一起用，不加换算直接把两个峰值读数相减，频率全错。Nyquist边界的问题则是扫描范围超过fs/2，导向矢量出现混叠，真实频率f=400 Hz的信号会在fs-f=124 Hz处出现镜像峰。

解决：定一套规则固定下来：所有现代法频谱分析代码里，内部统一使用Hz，导向矢量里除以fs；对外显示频率时，先检查峰值是否都小于fs/2，再打印。fs变量在项目开头定义一个全局常量，不要散落各处。AR谱的w = 2πf/fs也要换算成f再显示，我在3.2节的ar_spectrum函数里已经做了这个换算。

5.4 短段数据下的Burg稳定性

现象：数据只有64个点，Burg法算出的AR谱在部分频段是负功率，或者整个谱在低频段剧烈振荡。

原因：Burg算法里每一步的反射系数是通过分子分母求除法得到的，当数据段太短、预测残差能量很小时，除法数值不稳定。反射系数理论上被约束在[-1,1]，但浮点误差或数据非平稳会把它推出界外，导致模型失稳。

解决：短数据段优先用协方差法，它的最小二乘求解不涉及迭代递推，数值上更皮实；必须用Burg时，对反射系数加一个clip，让它不超过0.99。我在代码示例里没有加这个约束，真实工程中可以在k计算后加一行k = np.clip(k, -0.99, 0.99)，损失一点精度换稳定性，值得。

5.5 幅度估计的偏差

现象：AR谱的峰高和真实正弦功率差十几倍，而Capon谱的峰高也不等于功率。

原因：AR谱是模型谱，谱值表示“系统在这一点上的增益”，如果两个正弦频率靠得很近，极点之间相互影响，峰高天然被压低或抬高。Capon伪谱更是只反映滤波器输出功率比，和真实功率谱相差一个比例系数。

解决：不要用现代法频谱分析的峰高来做幅度估计。正确做法是：先用AR或MUSIC得到频率，再用最小二乘拟合x[n] = A sin(2πf n + φ)，把幅度和相位单独估出来，这样幅度精度能够逼近理论最优。

6. 现代法频谱分析的验证闭环：CRB对照与两个实用技巧

6.1 用Cramer-Rao下界判断估计好不好

所有频率估计的方差有一个理论下限，叫克拉美-罗下界（CRB）。单正弦加白噪声、数据长度为N时，频率估计的CRB近似为

var(f_est) ≥ 12 σ² / ((2π)² A² N³) × fs²

其中A是正弦幅度，σ²是白噪声方差。你的估计方差如果落到这个下界的10倍以内，算合格；如果差几个数量级，回去查阶数和预处理。这段代码可以对照实验结果：

def estimate_var_and_crb(peaks, true_freq, sig_amp, noise_var, fs, N): # 多次MUSIC测得的频率样本方差 freq_var = np.var(peaks) # 理论CRB，单位Hz的平方 crb = 12 * noise_var / ((2*np.pi)**2 * sig_amp**2 * N**3) * fs**2 return freq_var, crb # 模拟50次独立试验，统计频率估计稳定性 freq_estimates = [] for _ in range(50): x_trial = np.sin(2*np.pi*100*t) + 0.5*np.random.randn(len(t)) f_m, psd_m = music_frequency(x_trial, fs, num_sources=1, L=120) peak_idx = np.argmax(psd_m[:4096]) freq_estimates.append(f_m[peak_idx]) var_est, crb = estimate_var_and_crb(freq_estimates, 100, 1.0, 0.25, fs, len(t)) print(f"估计方差: {var_est:.6f}, 理论下界: {crb:.6f}")

CRB的意义在于：它告诉你当前数据量下，频率估计能有多准，不要把时间花在无意义的参数穷举上。如果方差离CRB很远，优先检查预处理和数据长度，而不是降阶或加网格。

6.2 抛物线插值细化频率与两个送分技巧

MUSIC谱峰在网格上是一个离散点，直接用argmax读数误差最多半个网格。一个便宜好用的改进是抛物线插值：取峰值点和左右两个相邻点的(f, psd)，拟合一条抛物线，顶点横坐标就是细化后的频率。这个技巧能零成本地把频率读数精度提高三到五倍。

def refine_peak(freqs, psd, idx): if idx == 0 or idx == len(freqs) - 1: return freqs[idx] f0, f1, f2 = freqs[idx-1], freqs[idx], freqs[idx+1] p0, p1, p2 = psd[idx-1], psd[idx], psd[idx+1] denom = (f0 - f1) * (f0 - f2) * (f1 - f2) a = (f0*(p2-p1) + f1*(p0-p2) + f2*(p1-p0)) / denom b = (f0**2*(p1-p2) + f1**2*(p2-p0) + f2**2*(p0-p1)) / denom return -b / (2*a) if a != 0 else f1

第二个技巧是留一法交叉验证：把数据分成前后两半，分别做现代法频谱分析，两个谱峰位置如果一致，当作可信；不一致，说明数据里有非平稳成分或瞬态干扰，先做分段处理再估计。第三个技巧是把现代法谱峰和FFT谱峰放在同一张图上看，FFT只能在粗尺度上验证，现代法在细尺度上读数，两者互相背书。

做完这些验证，你就不会在交付时被老板问“这个频率测得到底准不准”问得发慌。我的工作习惯是无论用什么算法，最后总要在同样信噪比和信号形状的仿真数据上过一遍方差，确信结果贴着CRB再拿去处理实测数据。这套流程让我在齿轮箱边频识别和电网谐波分析上少走了很多弯路，希望帮到你。

本文还有配套的精品资源，点击获取