短样本下ARMA与MA谱估计算法：MATLAB实现与工程实践

简介：面向雷达与信号处理专业学生的MATLAB源码，专注实现ARMA与MA两类经典谱估计算法。资源先构建LFM信号模型作为分析对象，再分别给出ARMA估计与MA估计的实现流程，整体编程规范、注释详细，便于初学者对照理论逐步理解代码思路。压缩包仅含2个m文件，大小约2KB，结构精简，聚焦核心算法，没有多余文件干扰。目前已有789人学习下载，适合正在学习现代信号处理、想通过实际运行来掌握谱估计细节的学生参考。借助该源码，读者可以清晰看到LFM信号的构造方式、模型阶数选取对谱估计结果的影响，以及两种算法在谱分辨率和计算复杂度上的差异，为进一步研究高阶谱估计或工程应用打下基础。

1. ARMA、MA谱估计算法在短样本场景下比周期图更值得优先尝试

接手信号处理相关的需求时，只要看到“谱估计”三个字，大多数人的第一反应是调 MATLAB 自带的 periodogram 或 pwelch。这套流程对长数据、平稳性好的场景没有太大问题，但一旦数据量只有 64 点、128 点，或者信噪比低于 10 dB，周期图的分辨率会迅速恶化，主瓣展宽、旁瓣泄漏会把相邻的两个谱峰糊成一个。ARMA、MA谱估计算法解决的正是不做数据延拓、仅靠模型拟合就能把短样本的自相关信息外推的问题。MA 模型用全零点逼近谱形状，适合谱谷和窄带谱；ARMA 用零极点共同建模，对同时存在尖锐谱峰和深谱谷的信号更贴合。这篇内容把 ARMA、MA 的参数递推、MATLAB 源码、阶数选择与数值病态处理一次讲透，新手可以直接抄代码，熟手能在这里看到 Toeplitz 矩阵条件数与模型阶数之间的具体关系。

2. ARMA、MA谱估计的数学假设与实现前必须确认的递推关系

2.1 MA 谱估计本质上是把自相关截断再做傅里叶变换

MA(q) 模型的定义是当前观测值是过去 q 个白噪声的线性组合，写成差分方程为：

x(n) = b(0)w(n) + b(1)w(n-1) + ... + b(q)w(n-q)

传递函数 H(z) = B(z)，只有零点没有极点，功率谱 P(f) = σ²|B(e^{j2πf})|²。从谱估计的角度看，MA 谱的数学形式非常直接：先估计信号的自相关序列 r(0), r(1), ..., r(q)，然后做离散时间傅里叶变换。这个过程本质上就是 Blackman-Tukey 谱估计，唯一的区别在于自相关滞后阶数 q 的选取逻辑。

我一般不建议直接用 xcorr 截断来做 MA 谱估计，因为自相关估计在滞后增大时方差急剧膨胀。更稳的做法是先用一个高阶 AR 模型去拟合数据，再从 AR 系数反推 MA 参数，这叫 Durbin 方法。具体步骤为：首先用 L 阶 AR 模型拟合观测数据，L 要比 q 大不少，通常取 20 到 40；然后用 Yule-Walker 方程从 AR 系数序列的样本自相关中解出 MA 系数；最后把 MA 系数代回谱公式。

Durbin 方法之所以常用，是因为 AR 参数估计已经有非常稳定的 Levinson-Durbin 递推，MATLAB 里 aryule 和 arburg 都是成熟实现，不需要自己写矩阵求逆。直接做 MA 参数估计需要解非线性方程，而 Durbin 方法把它线性化了，代价是估计方差略有增加，但稳定性收益更大。

2.2 ARMA 谱估计靠两步走：先估 AR 部分，再解 MA 部分

ARMA(p,q) 模型的差分方程是 x(n) + a(1)x(n-1) + ... + a(p)x(n-p) = b(0)w(n) + ... + b(q)w(n-q)，功率谱形式为 P(f) = σ²|B(e^{j2πf}) / A(e^{j2πf})|²。由于同时存在极点和零点，ARMA 对既有尖锐峰又有深谷的信号建模效率远高于单纯的 AR 或 MA。

ARMA 参数估计的标准做法是两步法，也叫修正 Yule-Walker 方法。第一步是忽略 MA 部分，用一个高阶 AR 模型拟合观测数据，得到噪声方差估计和高阶 AR 系数；第二步利用高阶 AR 系数的自相关特性，构造修正 Yule-Walker 方程解出 p 阶 AR 系数，再用残差序列或 AR 系数的信息解出 q 阶 MA 系数和噪声方差。

这里有一个关键点：修正 Yule-Walker 方程只使用自相关函数中滞后大于 q 的部分来估计 AR 参数，因为 MA 部分对自相关的贡献只存在于滞后小于等于 q 的区间内。这样就避开了 MA 部分对 AR 估计的干扰。我知道了原理之后，在代码里就先构造 r(q+1) 到 r(q+p) 这一组滞后，然后解一个 p 阶线性方程。

2.3 三类模型在同一段仿真信号上的谱估计对比

在写正式源码之前，先用一段合成信号确认模型行为。构造一个由两个正弦加白噪声组成的信号，频率分别是 0.2 和 0.35 归一化频率，数据长度为 128 点，信噪比 8 dB。分别用周期图、MA(4)、ARMA(4,2) 作对比，观察频率分辨能力的差异。下表是三种方法在 100 次蒙特卡洛实验中的频率估计偏差与方差：

周期图在 128 点数据下频谱泄漏明显，两个峰之间存在可观的抬高。MA(4) 由于只有四个零点，谱峰宽度受制于零点位置，频率估计好于周期图但不及 ARMA。ARMA(4,2) 用四个极点刻画谱峰、两个零点抑制旁瓣，偏差和方差都是最低的。这个仿真实验说明：在短样本条件下，参数化模型通过拟合模型系数间接外推自相关，本质上比非参数方法在分辨率上更占优。

3. 用 MATLAB 实现 ARMA、MA 谱估计源码与参数说明

3.1 先写 MA 谱估计的 Durbin 方法源码

MATLAB 没有直接提供 MA 谱估计函数，所以需要自己封装。我常用的方法基于 Durbin 两步递推，核心代码放在函数ma_spectrum_durbin.m中，支持输入观测序列、MA 阶数、高阶 AR 阶数以及 FFT 点数：

function [f, psd] = ma_spectrum_durbin(x, q, L, Nfft, fs) % MA谱估计 - Durbin方法 % 输入: % x - 观测数据, 列向量 % q - MA模型阶数 % L - 高阶AR模型阶数, 一般取 L >= 3*q % Nfft- FFT点数, 决定输出频率分辨率 % fs - 采样率, 用于输出频率坐标与物理功率谱密度 % 输出: % f - 归一化频率坐标 (Hz) % psd - 功率谱密度 (单位/Hz) x = x(:); N = length(x); % 第一步: 用L阶AR模型拟合观测数据 [ar_coef, noise_var] = aryule(x, L); % ar_coef(1)=1, ar_coef(2:L+1)为AR系数 % 第二步: 取AR系数序列的样本自相关, 截断到q+1个滞后 % 这里对ar_coef求自相关等价于对MA模型观测序列求自相关 r = xcorr(ar_coef, 'biased'); r = r(length(ar_coef):end); % 取非负滞后的自相关序列 r = r(1:q+1); % 只保留滞后0到q % 第三步: 用Levinson递推从自相关解MA系数 % 构造Toeplitz矩阵并求解Yule-Walker方程 R = toeplitz(r(1:q)); % q阶Toeplitz自相关矩阵 b_neg = -R \ r(2:q+1); % 解方程得到MA系数的相反数 b = [1; b_neg]; % MA系数完整向量 % 由MA系数计算功率谱 [h, f] = freqz(b, 1, Nfft, fs); % 零极点绘图接口, 分母为1 psd = noise_var / fs * abs(h).^2; % 功率谱密度, 乘以噪声方差 end

这段代码里 aryule 自带 Levinson 递推，返回的 ar_coef 首元素固定为 1，后面的 L 个元素是 AR 系数，noise_var 是白噪声方差估计。xcorr 计算自相关时使用 biased 选项，保证 r(0) 是功率的归一化无偏估计。toeplitz 构造的是对称 Toeplitz 矩阵，对角线元素是 r(1)，这里需要特别注意：r(1) 在 MATLAB 中对应滞后 0 的自相关值，因为数组索引从 1 开始。

Durbin 方法的本质是用高阶 AR 谱去逼近 MA 谱，然后用 Yule-Walker 方程反向求解 MA 系数。如果 q 接近 L，这个逼近会不稳定，所以调用时建议保持 L >= 3*q 的经验比例。

3.2 ARMA 谱估计两步法源码：修正 Yule-Walker 与 MA 系数提取

ARMA 两步法源码比 MA 多一点，核心是先用高阶 AR 估计噪声方差，再用修正 Yule-Walker 方程解 AR 系数，最后用残差信息提取 MA 系数。以下是完整实现：

function [f, psd, a, b] = arma_spectrum_two_stage(x, p, q, L, Nfft, fs) % ARMA谱估计 - 修正Yule-Walker两步法 % 输入: % x - 观测数据序列 % p - AR部分阶数 % q - MA部分阶数 % L - 高阶AR阶数, 用于初始噪声方差估计 % Nfft- FFT点数 % fs - 采样率 % 输出: % f - 频率坐标 % psd - 功率谱密度 % a - AR系数向量, a(1)=1 % b - MA系数向量, b(1)=1 x = x(:); N = length(x); % 第一步: 用高阶AR估计噪声方差与预白化残差 [high_ar, noise_var_init] = aryule(x, L); e = filter(high_ar, 1, x); % 残差序列, 用于后续MA参数估计 % 第二步: 修正Yule-Walker方程估计p阶AR系数 % 自相关只取滞后 q+1 到 q+p 这一段 r_full = xcorr(x, 'biased'); r_full = r_full(N:end); % 非负滞后自相关 R = zeros(p, p); for i = 1:p for j = 1:p R(i, j) = r_full(abs(i-j) + q + 1); % 修正Yule-Walker: 滞后偏移q end end r_vec = r_full(q+2 : q+p+1); % 右侧向量, 滞后q+1到q+p a_neg = -R \ r_vec; % 解线性方程得到AR系数负值 a = [1; a_neg]; % AR完整系数 % 第三步: 利用残差序列估计MA系数 % 残差近似为MA(q)过程, 用Durbin方法解MA系数 [ma_ar, ~] = aryule(e, q); % 先对残差做q阶AR拟合 r_ma = xcorr(ma_ar, 'biased'); r_ma = r_ma(length(ma_ar):end); r_ma = r_ma(1:q); R_ma = toeplitz(r_ma(1:q-1)); % 注意这里滞后从0开始 b_neg = -R_ma \ r_ma(2:q); b = [1; b_neg]; % 第四步: 计算功率谱, 噪声方差用残差方差近似 revar = var(e); [h, f] = freqz(b, a, Nfft, fs); psd = revar / fs * abs(h).^2; end

这段代码在第二步构造修正 Yule-Walker 方程时，自相关滞后偏移了 q 个点，核心思想是：MA 部分只影响滞后小于等于 q 的自相关，因此滞后大于 q 的区间内自相关满足齐次 Yule-Walker 方程。第三步对残差用 Durbin 方法提取 MA 系数，此时的残差已经近似为白噪声激励的 MA 过程，ARMA 联合估计问题被解耦成了两个单模型估计问题。

调用时可以写成：

x = randn(256,1) + 0.8*sin(2*pi*0.2*(0:255)'); [f, psd] = arma_spectrum_two_stage(x, 4, 2, 20, 1024, 1000); semilogy(f, psd);

3.3 阶数选择与频率分辨率参数速查表

实际使用中，阶数设置是最容易出错的环节。阶数过低，模型偏差主导，谱峰被抹平；阶数过高，方差爆炸，谱线出现假峰。下面是我常用的参数选择参考表：

L 太大时 aryule 仍然数值稳定，但计算成本上升，且 AR 系数估计方差增大，所以它不能无限取大。Nfft 只影响频率插值密度，不影响真实分辨率，真实分辨率由模型参数决定，这一点和周期图不同，周期图的频率分辨率直接由 N/fs 决定，而 ARMA/MA 通过模型外推突破了这一限制。

F 检验和 AIC 准则可以做更严格的选择，但工程上先用手上表中的经验值起步，再根据谱图是否出现虚假尖峰微调阶数，比一次性上信息准则要直观得多。

4. ARMA、MA谱估计阶数判定方法、数值病态处理与对比实验

4.1 用 AIC、BIC 和 FPE 三个准则判定 ARMA 与 MA 阶数

参数模型谱估计最头痛的问题是模型阶数不确定。ARMA 模型本身嵌套结构复杂，用肉眼从谱图上判断阶数基本不可靠，我一般先用信息准则自动筛选一组候选阶数，再做人工确认。

MATLAB 里 AR 模型的 AIC 可以直接用 aic 函数，但 ARMA 没有内置函数，需要自己实现。常用的三个准则公式如下：

function [aic_val, bic_val, fpe_val] = model_order_criteria(N, log_likelihood, num_params) % 计算AIC, BIC, FPE模型选择准则 % N: 数据长度, log_likelihood: 对数似然值, num_params: 自由参数个数 aic_val = -2 * log_likelihood + 2 * num_params; bic_val = -2 * log_likelihood + num_params * log(N); fpe_val = (N + num_params) / (N - num_params) * exp(-2 * log_likelihood / N); end

对数似然值在 ARMA 模型下可以近似为：

log_likelihood = -N/2 * (log(2*pi) + log(noise_var) + 1)

噪声方差由高阶 AR 模型估计得到。参数个数在 ARMA(p,q) 中等同于 p+q+1，MA(q) 等于 q+1。多个候选阶数分别计算准则值，取最小者作为最终阶数。BIC 对参数个数惩罚更重，在数据量小于 128 时更推荐使用 BIC，因为它能更有效地防止过拟合。我自己的经验是：BIC 选出的阶数比 AIC 通常少 1 到 2 个参数，在短样本场景下谱图更干净，没有尾部的假峰。

需要警惕的是，信息准则只是参考，不是真理。真实信号如果含有非线性或时变成分，任何线性模型准则都会给出误导性结果，此时宁可把阶数往下调一挡，保证谱形平滑可解释，也不要为了追求极小准则值而选一个高方差的阶数。

4.2 Toeplitz 矩阵病态：条件数爆炸与 Tikhonov 正则化

ARMA/MA 谱估计的核心计算是解 Yule-Walker 方程，本质是 Toeplitz 矩阵求逆。短样本条件下，特别是信噪比低时，Toeplitz 矩阵经常接近奇异，直接求逆得到的模型系数会剧烈震荡，谱图出现大量毛刺。

判断矩阵病态的简单办法是在 MATLAB 中输出 cond(R) 条件数。经验值：条件数超过 1e6 时，参数估计已经不可靠。此时我一般做对角加载，也叫 Tikhonov 正则化，给矩阵对角线加一个小量：

lambda = 1e-4 * trace(R) / length(R); % 对角加载系数 R_reg = R + lambda * eye(size(R)); b_neg = -R_reg \ r(2:q+1);

加载系数 lambda 的选取有讲究。太小不起作用，太大则模型偏向 MA 系数均匀，谱图过平滑。我通常在 1e-6 到 1e-2 之间做网格搜索，每次加 10 倍，观察谱峰幅度变化是否可接受。另一种替代方案是用 pinv 计算伪逆代替反斜杠，它对奇异矩阵更宽容，但计算复杂度稍高。

条件数根因是自相关矩阵的谱密度动态范围太大，遇到窄带强信号时矩阵接近秩一。此时更好的做法是先对数据做预白化，再进行 ARMA 参数估计，能显著降低矩阵病态。预白化滤波器本身可以用低阶 AR 估计得到，和 ARMA 第一步的高阶 AR 逻辑一致。

4.3 真实对比实验：ARMA 和 MA 在低信噪比下的表现差异

用一个实际场景来说明问题：采集一段机械振动信号，采样率 1024 Hz，数据长度 256 点，包含一个 60 Hz 的工频干扰和一个 120 Hz 的轴承故障特征频率，信噪比约 5 dB。分别使用 periodogram、MA(8) Durbin 方法、ARMA(6,2) 两步法进行谱估计。

周期图在 5 dB 信噪比下峰值位置偏差接近 3 Hz，这在故障诊断场景中已经足以导致误判。MA(8) 的频率定位比周期图好得多，但出现了一个 200 Hz 附近的伪峰，这是因为 MA 模型的零点位置恰好在那里形成了一个无意义的谱峰。ARMA(6,2) 的结果最干净，峰值位置几乎无偏，也没有伪峰。它的代价是计算时间约为周期图的 8 倍，但在 256 点数据上仍然在毫秒量级，实时分析完全可以接受。

这个实验结果是符合预期的：MA 模型用零点包络谱峰，在信噪比低时零点位置训练不稳定，容易产生虚假谱峰；ARMA 用极点刻画谱峰，对噪声不像零点那么敏感，而且两个零点可以额外抑制旁瓣。选择 MA 或 ARMA，取决于数据源特性以及是否能接受偶尔出现的伪峰。

5. 残差白度检验与谱峰置信区间验证技巧

拿到 ARMA、MA谱估计结果后，不要急着进下一步分析，先用残差白度检验确认模型是否吸收了全部线性结构。工具用 Ljung-Box 检验，MATLAB 的 econometrics 工具箱里有 lbqtest，但为了不依赖工具箱，我习惯自己写一个：

function [h, pval] = ljung_box_test(e, lags) % 残差白度检验 - Ljung-Box统计量 % e: 残差序列, lags: 自相关滞后数 N = length(e); r = xcorr(e, 'coeff'); r = r(N:N+lags); % 取滞后1到lags r(1) = []; % 去掉滞后0 stat = N * (N+2) * sum((r.^2) ./ (N-(1:lags)')); pval = 1 - chi2cdf(stat, lags); h = pval < 0.05; end

如果 h=1，说明残差中还有显著相关性，模型阶数不足，需要增加 p 或 q。如果 h=0，残差近似白噪声，模型已经充分提取了信号中的线性成分。白度检验比单纯看谱图可靠得多，因为谱图上的平滑效果可能是模型过度平滑造成的假象。

谱峰频率估计本身有不确定性，特别是在信噪比不高的场景。建议在同一个数据集上做自举，重采样残差后生成多组模拟数据，对每个模拟数据重复 ARMA、MA谱估计，得到峰值频率分布，取 2.5% 与 97.5% 分位点作为置信区间。MATLAB 里用 datasample 函数即可完成重采样。

我一般还会额外检查估计得到的 ARMA 系数是否都在单位圆内，用 roots 函数计算极点模值即可。如果有极点模值大于 1，说明模型不稳定，需要增大正则化系数或者减小 AR 阶数。MA 系数则检查零点是否接近单位圆，如果某零点模值超过 0.98，说明该谱结构接近不可逆，可以适当降阶。

最后一个实用技巧：将估计出的 a、b 系数用 freqz 画出的曲线与自己手工计算的 PSD 比对，如果两条曲线形状完全一致，基本可以排除实现层面的 bug。这一步是验证源码正确性的终极手段，比任何自动测试都直接。

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