随机信号的参数建模法( AR 模型)

引言

为随机信号建立参数模型是研究随机信号的一种基本方法，其含义是认为随机信号$x（n）$x（n）是由白噪$w（n）$w（n）激励某一确定系统的响应（如图7.5）。只要白噪的参数确定了，研究随机信号就可以转化成研究产生随机信号的系统。

图 1 随机信号的参数模型

对平稳随机信号，三种常用的线性模型分别是 AR 模型（自回归模型 Auto-regression model），MA 模型（滑动平均模型 Moving average model）和 ARMA 模型（自回归滑移平均模型 Auto-regression-Moving average model）。 *根据$ Wold $的证明：任何平稳的 ARMA（自回归移动平均）模型或 MA 模型均可用无限阶或阶数足够的 AR 模型去近似。*因此介绍 AR 模型的基本原理和方法。

AR 模型概述

随机信号$x（n）$x（n）由本身的若干次过去值$x（n-k）$x（n−k）和当前的激励值$w（n）$w（n）线性组合产生：
$x(n)=w(n)-\sum_{k=1}^{p} a_{k} x(n-k)$x(n)=w(n)−k=1∑p​ak​x(n−k)
该模型的系统函数是：
$H(z)=\frac{1}{1+\sum_{k=1}^{p} a_{k} z^{-k}}$H(z)=1+∑k=1p​ak​z−k1​
$p $是系统阶数，系统函数中只有极点，无零点，也称为全极点模型，系统由于极点的原因， 要考虑到系统的稳定性，因而要注意极点的分布位置，用 $AR（ p ）$AR（p）来表示。

AR 模型参数的估计

AR 模型参数和自相关函数的关系

证明过程

1. 求自相关函数$R_{x x}(m)$Rxx​(m)

已知随机信号 $x(n)=w(n)-\sum_{k=1}^{p} a_{k} x(n-k)$x(n)=w(n)−∑k=1p​ak​x(n−k)，对该式两边同时乘以$x（n-m）$x（n−m），然后求均值：
$E[x(n) x(n-m)]=E\left[w(n) x(n-m)-\sum_{k=1}^{p} a_{k} x(n-k) x(n-m)\right]\tag{1.1}$E[x(n)x(n−m)]=E[w(n)x(n−m)−k=1∑p​ak​x(n−k)x(n−m)](1.1)
又因为自相关函数呈现偶对称：$R_{x x}(m)=R_{x x}(-m)$Rxx​(m)=Rxx​(−m)

所以：
$R_{x x}(m)=R_{x w}(m)-\sum_{k=1}^{p} a_{k} R_{x x}(m-k) \tag{1.2}$Rxx​(m)=Rxw​(m)−k=1∑p​ak​Rxx​(m−k)(1.2)
2. 求互相关函数$R_{x w}$Rxw​:
$\begin{aligned} R_{x w}(m)&=E[x(n) w(n+m)] \\&=E\left[w(n) * h(n) w(n+m)\right] \\&=E\left[\sum_{k=0}^{\infty} h(k) w(n-k) w(n+m)\right] \\&=\sum_{k=0}^{\infty} h(k) E[w(n-k) w(n+m)] \\&=\sum_{k=0}^{\infty} h(k) R_{w w}(m+k) \\&=\sum_{k=0}^{\infty} h(k) \sigma_{w}^{2} \delta(m+k) \\&=\sigma_{w}^{2} h(-m) \end{aligned}\tag{1.3}$Rxw​(m)​=E[x(n)w(n+m)]=E[w(n)∗h(n)w(n+m)]=E[k=0∑∞​h(k)w(n−k)w(n+m)]=k=0∑∞​h(k)E[w(n−k)w(n+m)]=k=0∑∞​h(k)Rww​(m+k)=k=0∑∞​h(k)σw2​δ(m+k)=σw2​h(−m)​(1.3)

因为系统的单位脉冲响应是因果的, $w(n)$w(n)是具有方差$\sigma_{w}^{2}$σw2​ 的平稳白噪声

所以：
$R_{x w}(m)=\left\{\begin{array}{cc} 0 & m>0 \\ \sigma_{w}^{2} h(-m) & m \leq 0 \end{array}\right.\tag{1.4}$Rxw​(m)={0σw2​h(−m)​m>0m≤0​(1.4)
3. 结合(1.2)和(1.4)得：
$R_{x x}(m)=\left\{\begin{array}{cc} -\sum_{k=1}^{p} a_{k} R_{x x}(m-k) & m>0 \\ -\sum_{k=1}^{p} a_{k} R_{x x}(m-k)+h(0) \sigma_{w}^{2} & m=0 \\ R_{x x}(-m) & m<0 \end{array}\right.\tag{1.5}$Rxx​(m)=⎩⎨⎧​−∑k=1p​ak​Rxx​(m−k)−∑k=1p​ak​Rxx​(m−k)+h(0)σw2​Rxx​(−m)​m>0m=0m<0​(1.5)
又因为：
$H(z)=\frac{1}{1+\sum_{k=1}^{p} a_{k} z^{-k}} \\h(n)+\sum_{k=1}^{p} a_{k} h(n-k)=\delta(n)$H(z)=1+∑k=1p​ak​z−k1​h(n)+k=1∑p​ak​h(n−k)=δ(n)
即$h(0)=1$h(0)=1，所以：
$R_{x x}(m)=\left\{\begin{array}{cc} -\sum_{k=1}^{p} a_{k} R_{x x}(m-k) & m>0 \\ -\sum_{k=1}^{p} a_{k} R_{x x}(m-k)+ \sigma_{w}^{2} & m=0 \\ R_{x x}(-m) & m<0 \end{array}\right.\tag{1.6}$Rxx​(m)=⎩⎨⎧​−∑k=1p​ak​Rxx​(m−k)−∑k=1p​ak​Rxx​(m−k)+σw2​Rxx​(−m)​m>0m=0m<0​(1.6)
可以得知：AR 模型输出信号的自相关函数具有递推的性质，自相关函数是偶对称函数。

结论及例题

AR 模型参数和自相关函数的关系为：
$\left[\begin{array}{cccc} R(0) & R(-1) & \cdots & R(-p) \\ R(1) & R(0) & \cdots & R(-p+1) \\ \vdots & \vdots & \cdots & \vdots \\ R(p) & R(p-1) & \cdots & R(0) \end{array}\right]\left[\begin{array}{c} 1 \\ a_{1} \\ \vdots \\ a_{p} \end{array}\right]=\left[\begin{array}{c} \sigma_{w}^{2} \\ 0 \\ \vdots \\ 0 \end{array}\right]\tag{1.7}$⎣⎢⎢⎢⎡​R(0)R(1)⋮R(p)​R(−1)R(0)⋮R(p−1)​⋯⋯⋯⋯​R(−p)R(−p+1)⋮R(0)​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​1a1​⋮ap​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​σw2​0⋮0​⎦⎥⎥⎥⎤​(1.7)
式（1-7）方程组的系数都是自相关矩阵$[R]_{p+1}$[R]p+1​，由于自相关函数是偶对称函数：$Re（m）=Re（-m）$Re（m）=Re（−m），因而自相关矩阵是对称矩阵，与主对角线平行的斜对角线上的元素都是相同的，是$（p+1）×（p+1）$（p+1）×（p+1）的维托毕兹（Toeplitz）矩阵，所以存在高效算法，其中应用广泛的有Levinson-Durbin（L-D）算法。Yule-Walker（Y-W）方程表明，只要已知输出平稳随机信号的自相关函数，就能求出AR模型中的参数$\left\{a_{k}\right\}${ak​}，并且需要的观测数据较少。

【例1】已知自回归信号模型AR（3）为：
$x(n)=\frac{14}{24} x(n-1)+\frac{9}{24} x(n-2)-\frac{1}{24} x(n-3)+w(n)$x(n)=2414​x(n−1)+249​x(n−2)−241​x(n−3)+w(n)
其中w（n）是方差${\sigma}_{w}^{2}=1$σw2​=1的平稳白噪声，求：
1.自相关序列$R_{x x}(m)$Rxx​(m)，m=0，1，2，3，4，5；
2.用上一问求出的自相关序列来估计AR（3）的参数$\left\{\hat{a}_{k}\right\}${a^k​}，以及输入白噪声的方差$\hat{\sigma}_{w}^{2}$σ^w2​大小；
3.利用给出的AR模型，用计算机仿真给出32点观测值x（n）＝[0.4282 1.1454
1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242 0.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.5802 -0.6177]，用观测值的自相关序列直接来估计AR（3）的参数$\left\{\hat{a}_{k}\right\}${a^k​}以及输入白噪声的方差$\hat{\sigma}_{w}^{2}$σ^w2​。

1. $a_{1}=-\frac{14}{24}, a_{2}=-\frac{9}{24}, a_{3}=\frac{1}{24},{\sigma}_{w}^{2}=1$a1​=−2414​,a2​=−249​,a3​=241​,σw2​=1代入公式(1.7)得：

$\left[\begin{array}{llll}R(0) & R(1) & R(2) & R(3) \\ R(1) & R(0) & R(1) & R(2) \\ R(2) & R(1) & R(0) & R(1) \\ R(3) & R(2) & R(1) & R(0)\end{array}\right]\left[\begin{array}{c}1 \\ -14 / 24 \\ -9 / 24 \\ 1 / 24\end{array}\right]=\left[\begin{array}{l}1 \\ 0 \\ 0 \\ 0\end{array}\right]$⎣⎢⎢⎡​R(0)R(1)R(2)R(3)​R(1)R(0)R(1)R(2)​R(2)R(1)R(0)R(1)​R(3)R(2)R(1)R(0)​⎦⎥⎥⎤​⎣⎢⎢⎡​1−14/24−9/241/24​⎦⎥⎥⎤​=⎣⎢⎢⎡​1000​⎦⎥⎥⎤​

• 利用matlab解线性方程组得： $\mathrm{R}(0)=4.9377 \quad \mathrm{R}(1)=4.3287 \quad \mathrm{R}(2)=4.1964 \quad \mathrm{R}(3)=3.8654$R(0)=4.9377R(1)=4.3287R(2)=4.1964R(3)=3.8654

syms R0 R1 R2 R3
A = [1
     -14/24
     -9/24
     1/24 ];
B = [1
     0
     0
     0];
R = [ R0 R1 R2 R3
        R1 R0 R1 R2
        R2 R1 R0 R1
        R3 R2 R1 R0];
[R0 ,R1, R2, R3]=solve(R*A==B);
Rxx=vpa([R0 ,R1, R2, R3],5)

• 利用式（1.6） $R_{x x}(m)=-\sum_{k=1}^{p} a_{k} R_{x x}(m-k) \quad m>0,$Rxx​(m)=−∑k=1p​ak​Rxx​(m−k)m>0, 可以求出 $\mathrm{R}(4), \mathrm{R}(5) \cdots$R(4),R(5)⋯

for m = 5 : 6
    Rxx(m) = 0;
    for k = 1 : 3
        Rxx(m) = Rxx(m) - A(k+1) * Rxx(m - k);
    end
end
Rxx=vpa(Rxx,5)

• 结果

Rxx =

	[ 4.9377, 4.3287, 4.1964, 3.8654, 3.6481, 3.4027]

• 代码

Rxx =[ 4.9377, 4.3287, 4.1964, 3.8654, 3.6481, 3.4027];
R = [1 0 0 0
     Rxx(2) Rxx(1) Rxx(2) Rxx(3)
     Rxx(3) Rxx(2) Rxx(1) Rxx(2)
     Rxx(4) Rxx(3) Rxx(2) Rxx(1)];
B = [ 1
      0
      0
      0];
A = R\B
R(1,:) = [Rxx(1) Rxx(2) Rxx(3) Rxx(4)];
sigma = R(1,:) * A 

• 结果

A =

    1.0000
   -0.5833
   -0.3751
    0.0417


sigma =

    1.0000

可以发现对 AR 模型参数是无失真的估计，因为已知 AR 模型，我们可以得到完全的输出观 测值，因而求得的自相关函数没有失真，当然也就可以不失真的估计。

3. 离散序列的自相关公式为：
   $R_{xx}(m)=\frac{1}{N-m} \sum_{n=1}^{N-m} X_{n} X_{n+m}, m=0,1, \ldots \ldots M$Rxx​(m)=N−m1​n=1∑N−m​Xn​Xn+m​,m=0,1,……M

根据公式可以求得前32点的自相关序列：

xn = [0.4282    1.1454    1.5597    1.8994    1.6854    2.3075    2.4679    1.9790   
      1.6063    1.2804   -0.2083    0.0577    0.0206    0.3572    1.6572    0.7488    
      1.6666    1.9830    2.6914    1.2521    1.8691    1.6855    0.6242    0.1763    
      1.3490    0.6955    1.2941     1.0475    0.4319     0.0312    0.5802   -0.6177];
N=length(xn(:));

% xcorr可以直接计算自相关函数
xn = reshape (xn',1,[]);
Rxx=xcorr(xn)/N;
Rxx(N : end)

% for m = 0:N-1
%     Rx_1 = zeros (1,N);
%     Rx_2 = zeros (1,N);
%     Rx_1(1:N - m) = xn (1:N - m);
%     Rx_2(1:N - m) = xn (1 + m:N);
%     Rxx(m+1) = sum( Rx_1 .* Rx_2 ) / (N-m) ;
% end
% Rxx

ans =

  1 至 15 列

    1.9271    1.6618    1.5381    1.3545    1.1349    0.9060    0.8673    0.7520    0.7637    0.8058    0.8497    0.8761    0.9608    0.8859    0.7868

  16 至 30 列

    0.7445    0.6830    0.5808    0.5622    0.5134    0.4301    0.3998    0.3050    0.2550    0.1997    0.1282    0.0637    0.0329   -0.0015   -0.0089

  31 至 32 列

   -0.0143   -0.0083

把头 4 个相关序列值代入公式(1.7)，利用第二问的代码求得估计值：
$\hat{a}_{1}=-0.6984 \quad \hat{a}_{2}=-0.2748 \quad \hat{a}_{3}=0.0915, \quad \hat{\sigma}_{w}^{2}=0.4678$a^1​=−0.6984a^2​=−0.2748a^3​=0.0915,σ^w2​=0.4678
与真实 AR 模型参数误差为：$e_{1}=0.1151, \quad e_{2}=0.1002, \quad e_{3}=0.0498$e1​=0.1151,e2​=0.1002,e3​=0.0498，原因在于我们只有一部分的观测数据，使得自相关序列值与理想的完全不同，输入信号的方差误差比较大：$e_{\sigma}=0.5322$eσ​=0.5322。给出一段观测值求 AR 模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

通过L-D算法来估计$AR$AR模型的参数$\left\{a_{k}, p, \sigma_{w}^{2}\right\}${ak​,p,σw2​}

AR模型和差分预测的关系

AR模型可以通过差分预测的的方式计算得到，但二者本质没有关系，只是形式恰好相同

前向预测误差系统：
$e(n)=x(n)-\hat{x}(n)=x(n)+\sum_{k=1}^{m} a_{m}(k) x(n-k)$e(n)=x(n)−x^(n)=x(n)+k=1∑m​am​(k)x(n−k)
平稳随机信号的AR模型：
$x(n)=w(n)-\sum_{k=1}^{p} a_{k} x(n-k)$x(n)=w(n)−k=1∑p​ak​x(n−k)
假如 m＝p，$w(n)=e(n)$w(n)=e(n)，且预测系数和 AR 模型参数相同，则两者的系统函数互为倒数， 如图 2 所示：

图 2 预测误差系统和 AR 模型

所以求 AR 模型参数就可以通过求预测误差系统的预测系数来实现，对(1.6)式改写：
$R_{x x}(l)=-\sum_{k=1}^{m} a_{m}(k) R_{x x}(l-k) \quad l=1,2, \cdots m \tag{2.3}$Rxx​(l)=−k=1∑m​am​(k)Rxx​(l−k)l=1,2,⋯m(2.3)

$E_{m}\left[e^{2}(n)\right]=R_{x x}(0)+\sum_{k=1}^{m} a_{m}(k) R_{x x}(k) \quad或\quad E_{p}\left[e^{2}(n)\right]=R_{x x}(0)+\sum_{k=1}^{p} a_{k} R_{x x}(k) \\其中a_{k}=a_{m}(k) \quad m=p \quad m是预测器的阶数 \tag{2.4}$Em​[e2(n)]=Rxx​(0)+k=1∑m​am​(k)Rxx​(k)或Ep​[e2(n)]=Rxx​(0)+k=1∑p​ak​Rxx​(k)其中ak​=am​(k)m=pm是预测器的阶数(2.4)

也就是说：自相关序列具有递推的性质；$ p$ 阶预测器的预测系数等于$ p$ 阶 $AR$AR 模型的参数，由于 $e(n)=w(n)$e(n)=w(n), 所以最小均方预测误差等于白噪声方差，即 $E_{p}\left[e^{2}(n)\right]=\sigma_{w}^{2}$Ep​[e2(n)]=σw2​。

L-D算法基本思想

L-D 递推算法是模型阶数逐渐加大的一种算法：利用式(2.3)、(2.4)自相关序列具有递推的性质，先计算阶次 $\mathrm{m}=1$m=1 时的预测系数 $\left\{a_{m}(k)\right\}=a_{1}(1)${am​(k)}=a1​(1) 和 $\sigma_{w 1}^{2},$σw12​, 然后计算 $\mathrm{m}=2$m=2 时的 $\left\{a_{m}(k)\right\}=a_{2}(1),a_{2}(2)${am​(k)}=a2​(1),a2​(2) 以及 $\sigma_{w 2}^{2},$σw22​, 一直计算到 $\mathrm{m}=\mathrm{p}$m=p 阶时的 $a_{p}(1), a_{p}(2), \cdots, a_{p}(p)$ap​(1),ap​(2),⋯,ap​(p) 以及 $\sigma_{w p}^{2}$σwp2​ 。这种递推算法的特点是，每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。

推导过程

结合式(2.3)、(2.4)，递推过程如下：

m=1:
$\left\{\begin{array}{c} R(1)=-a_{1}(1) R(0) \quad (1)\\ E_{1}=R(0)+a_{1}(1) R(1) \quad (2) \end{array} \quad \longrightarrow \quad \sigma_{w 1}^{2}=E_{1}=R(0)\left[1-a_{1}^{2}(1)\right]\right.\tag{2.5}${R(1)=−a1​(1)R(0)(1)E1​=R(0)+a1​(1)R(1)(2)​⟶σw12​=E1​=R(0)[1−a12​(1)](2.5)

m=2：
$\left\{\begin{array}{l} R(1)=-a_{2}(1) R(0)-a_{2}(2) R(1) \\ R(2)=-a_{2}(1) R(1)-a_{2}(2) R(0) \end{array}\right.${R(1)=−a2​(1)R(0)−a2​(2)R(1)R(2)=−a2​(1)R(1)−a2​(2)R(0)​

把（2.5）的（1）代入上式得到：

$\left\{\begin{array}{c} a_{2}(1)=\frac{-R(1)-a_{2}(2) R(1)}{R(0)}=a_{1}(1)+a_{2}(2) a_{1}(1) \quad\quad\quad (1)\\ a_{2}(2)=-\frac{R(2)+a_{2}(1) R(1)}{R(0)}=-\frac{R(2)+\left[a_{1}(1)+a_{2}(2) a_{1}(1)\right] R(1)}{R(0)} \quad (2) \end{array}\right. \\\Rightarrow a_{2}(2)=-\frac{R(2)+a_{1}(1) R(1)}{R(0)+a_{1}(1) R(1)}=-\frac{R(2)+a_{1}(1) R(1)}{E_{1}} \quad (3) \tag{2.6}${a2​(1)=R(0)−R(1)−a2​(2)R(1)​=a1​(1)+a2​(2)a1​(1)(1)a2​(2)=−R(0)R(2)+a2​(1)R(1)​=−R(0)R(2)+[a1​(1)+a2​(2)a1​(1)]R(1)​(2)​⇒a2​(2)=−R(0)+a1​(1)R(1)R(2)+a1​(1)R(1)​=−E1​R(2)+a1​(1)R(1)​(3)(2.6)

根据(2.4)，估计的方差为：
KaTeX parse error: Unknown column alignment: 1 at position 16: \begin{array}{1̲} E_{2}&=R(0)+a…
把(2.6)的（1）(2)代人:
$\begin{array}{l} E_{2}&=R(0)+a_{2}(1)\left[1-a_{2}(2)\right] R(1)-a_{2}^{2}(2) R(0) \\ &=R(0)+a_{1}(1)\left[1+a_{2}(2)\right]\left[1-a_{2}(2)\right] R(1)-a_{2}^{2}(2) R(0) \\ &=\left[1-a_{2}^{2}(2)\right]\left[R(0)+a_{1}(1) R(1)\right]=\left[1-a_{2}^{2}(2)\right] E_{1} \end{array} \tag{2.7}$E2​​=R(0)+a2​(1)[1−a2​(2)]R(1)−a22​(2)R(0)=R(0)+a1​(1)[1+a2​(2)][1−a2​(2)]R(1)−a22​(2)R(0)=[1−a22​(2)][R(0)+a1​(1)R(1)]=[1−a22​(2)]E1​​(2.7)
这样递推下去可以得到预测系数和均方误差估计的通式。

结论

预测系数和均方误差估计的通式：
$\left\{\begin{array}{c} a_{m}(k)=a_{m-1}(k)+a_{m}(m) a_{m-1}(m-k) \quad (1) \\ a_{m}(m)=-\frac{R(m)+\sum_{k=1}^{m-1} a_{m-1}(k) R(m-k)}{E_{m-1}} \quad (2) \\ E_{m}=\sigma_{w m}^{2}=\left[1-a_{m}^{2}(m)\right] E_{m-1}=R(0) \prod_{k=1}^{m}\left[1-a_{k}^{2}(k) \right] \quad (3) \end{array}\right. \tag{2.8}$⎩⎪⎨⎪⎧​am​(k)=am−1​(k)+am​(m)am−1​(m−k)(1)am​(m)=−Em−1​R(m)+∑k=1m−1​am−1​(k)R(m−k)​(2)Em​=σwm2​=[1−am2​(m)]Em−1​=R(0)∏k=1m​[1−ak2​(k)](3)​(2.8)
其中 $a_{m}(m)$am​(m) 称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值$E_{0}=R(0), \quad a_{0}(0)=1, \quad$E0​=R(0),a0​(0)=1, 以及 $\mathrm{AR}$AR 模型的阶数 $p,$p,。

参数估计流程图，如图3所示。

图3 L-D 算法流程图

L-D 算法的优缺点

• 优点：计算速度快，求得的 AR 模型必定稳定，且均方预测误差随着阶次的增加而减小（见(2.8)的(3)式）。

• 缺点，由于在求自相关序列时，是假设除了观测值之外的数据都为零，必然会引入较大误差。

例题

【例2】已知自回归信号模型AR（3）为：
$x(n)=\frac{14}{24} x(n-1)+\frac{9}{24} x(n-2)-\frac{1}{24} x(n-3)+w(n)$x(n)=2414​x(n−1)+249​x(n−2)−241​x(n−3)+w(n)
其中w（n）是方差${\sigma}_{w}^{2}=1$σw2​=1的平稳白噪声，利用给出的AR模型，用计算机仿真给出32点观测值x（n）＝[0.4282 1.1454
1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242 0.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.5802 -0.6177]，用观测值的自相关序列直接来估计AR（3）的参数$\left\{\hat{a}_{k}\right\}${a^k​}以及输入白噪声的方差$\hat{\sigma}_{w}^{2}$σ^w2​。

按照图3 算法流程图 的思路进行计算

• 步骤一：计算自相关函数（同例1.3）$R_{x x}(m)=\left[\begin{array}{llll}1.9271 & 1.6618 & 1.5381 & 1.3545 \cdots ]\end{array}\right.$Rxx​(m)=[1.9271​1.6618​1.5381​1.3545⋯]​

• 步骤二：初始化: $E_{0}=R_{x x}(0)=1.9271, \quad a_{0}=1$E0​=Rxx​(0)=1.9271,a0​=1

• 步骤三：根据公式2.8进行递推，直到满足精度为止。

$\begin{array}{l} m=1:\left\{\begin{array}{l} a_{1}(1)=-\frac{R(1)}{E_{0}}=-\frac{1.6618}{1.9271}=-0.8623 \\ E_{1}=R(0)\left[1-a_{1}^{2}(1)\right]=0.4942 \end{array}\right. \\ \begin{array}{l} m=2:\left\{\begin{array}{l} a_{2}(2)=-\frac{R(2)+a_{1}(1) R(1)}{E_{1}}=-0.2127 \\ a_{2}(1)=a_{1}(1)\left[1+a_{2}(2)\right]=-0.6789 \\ E_{2}=E_{1}\left[1-a_{2}^{2}(2)\right]=0.4718 \end{array}\right. \\ m=3:\left\{\begin{array}{l} a_{3}(3)=-\frac{R(3)+a_{2}(1) R(2)+a_{2}(2) R(1)}{E_{2}}=0.0914 \\ a_{3}(1)=a_{2}(1)+a_{3}(3) a_{2}(2)=-0.6983 \\ a_{3}(2)=a_{2}(2)+a_{3}(3) a_{2}(1)=-0.2748 \\ E_{3}=E_{2}\left[1-a_{3}^{2}(3)\right]=0.4679 \end{array}\right. \end{array} \end{array}$m=1:{a1​(1)=−E0​R(1)​=−1.92711.6618​=−0.8623E1​=R(0)[1−a12​(1)]=0.4942​m=2:⎩⎨⎧​a2​(2)=−E1​R(2)+a1​(1)R(1)​=−0.2127a2​(1)=a1​(1)[1+a2​(2)]=−0.6789E2​=E1​[1−a22​(2)]=0.4718​m=3:⎩⎪⎪⎨⎪⎪⎧​a3​(3)=−E2​R(3)+a2​(1)R(2)+a2​(2)R(1)​=0.0914a3​(1)=a2​(1)+a3​(3)a2​(2)=−0.6983a3​(2)=a2​(2)+a3​(3)a2​(1)=−0.2748E3​=E2​[1−a32​(3)]=0.4679​​​

因而当 $p=3$p=3 时, 预估到的 $\mathrm{AR}$AR 模型参数为 $: \hat{a}_{1}=-0.6983, \hat{a}_{2}=-0.2748, \hat{a}_{3}=0.0914$:a^1​=−0.6983,a^2​=−0.2748,a^3​=0.0914。估计的输入信号的方差为: $\hat{\sigma}_{w}^{2}=E_{3}=0.4679$σ^w2​=E3​=0.4679 。和例题 $1$1 的 $\mathrm{c}$c 结果一致，误差分析也一样，假设满足精度，停止递推。

当要计算的阶数比较高时，可以利用递推程序来实现。Matlab有专门的函数实现 L-D 算法的 AR 模型参数估计: [a E]=aryule（x，p），输入 x表示观测信号，输入 p表示要求的阶数，输出 a 表示佔计的模型参数，e 表示噪声信号的方差估让。

[a, E] = aryule(x_n,3);
a
E

• 结果

a =

    1.0000   -0.6984   -0.2748    0.0915


E =

    0.4678

假如使用更高阶的 AR 模型来估计，均方误差更小，结果越精确。

AR 模型参数估计的各种算法的比较和阶数的选择

L-D 算法由于在求自相关序列时，是假设除了观测值之外的数据都为零，必然会引入较大误差；阶数$p$p越高，精度越大，但算法复杂性也随之上升。为了解决以上问题，分别对 $AR$AR 模型参数有着不同的优化。

Burg 算法

基本思想:

对观测的数据进行前向和后向预测，然后让两者的均方误差之和为最小作为估计的准则估计处反射系数$a_{m}(m)$am​(m),，进而通过 L-D 算法的递推公式求出 AR 模型参数。相比于L-D 算法只进行前向预测而言，精准度和可靠性更高，但本质没有区别。

代码实现:

[a, E] = arburg(x_n,p);
a =

    1.0000   -0.6852   -0.3088   -0.0489    0.1759


E =

    0.4426

可以看到，$\hat{a}_{3}$a^3​ 佔计得更精确些。

Marple 算法

基本思想:

与自适应算法类似，每一个预测系数的确定直接与前向、后向预测的总的平方误差最小，而不是由低一阶的系数递推确定。由于该算法是从整体上选择所有的模型参数达到总的均方误差最小，所以比 L-D、Burg 算法优越，但无法保证 AR 模型的稳定性。

关键在于阶数的选择：

最终预测误差（FPE:final predidyion error）准则定义：给定观测长度为N，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差，该误差在某个阶数p时为最小。其表达式为：
$F P E(p)=\hat{\sigma}_{w p}^{2}\left(\frac{N+P+1}{N-P-1}\right) \tag{3.1}$FPE(p)=σ^wp2​(N−P−1N+P+1​)(3.1)
除了用式(3.1)进行阶数的选择，还可以采用试验的方法，在短数据情况下，根据经验，AR 模型的阶次选在N/3～N/2 的范围内较好。