引言

为随机信号建立参数模型是研究随机信号的一种基本方法，其含义是认为随机信号x（n）x（n）x（n）是由白噪w（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）由本身的若干次过去值x（n−k）x（n-k）x（n−k）和当前的激励值w（n）w（n）w（n）线性组合产生：

x(n)=w(n)−∑k=1pakx(n−k)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)=11+∑k=1pakz−kH(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（p）来表示。

AR 模型参数的估计

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

证明过程

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

已知随机信号 x(n)=w(n)−∑k=1pakx(n−k)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）x（n−m），然后求均值：

E[x(n)x(n−m)]=E[w(n)x(n−m)−∑k=1pakx(n−k)x(n−m)](1.1)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)

又因为自相关函数呈现偶对称：Rxx(m)=Rxx(−m)R_{x x}(m)=R_{x x}(-m)Rxx​(m)=Rxx​(−m)

所以：

Rxx(m)=Rxw(m)−∑k=1pakRxx(m−k)(1.2)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. 求互相关函数RxwR_{x w}Rxw​:

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)=σw2h(−m)(1.3)\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)w(n)是具有方差σw2\sigma_{w}^{2}σw2​ 的平稳白噪声

所以：

Rxw(m)={0m>0σw2h(−m)m≤0(1.4)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)得：

Rxx(m)={−∑k=1pakRxx(m−k)m>0−∑k=1pakRxx(m−k)+h(0)σw2m=0Rxx(−m)m<0(1.5)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)=11+∑k=1pakz−kh(n)+∑k=1pakh(n−k)=δ(n)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)=1h(0)=1h(0)=1，所以：

Rxx(m)={−∑k=1pakRxx(m−k)m>0−∑k=1pakRxx(m−k)+σw2m=0Rxx(−m)m<0(1.6)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 模型参数和自相关函数的关系为：

[R(0)R(−1)⋯R(−p)R(1)R(0)⋯R(−p+1)⋮⋮⋯⋮R(p)R(p−1)⋯R(0)][1a1⋮ap]=[σw20⋮0](1.7)\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}[R]p+1​，由于自相关函数是偶对称函数：Re（m）=Re（−m）Re（m）=Re（-m）Re（m）=Re（−m），因而自相关矩阵是对称矩阵，与主对角线平行的斜对角线上的元素都是相同的，是（p+1）×（p+1）（p+1）×（p+1）（p+1）×（p+1）的维托毕兹（Toeplitz）矩阵，所以存在高效算法，其中应用广泛的有Levinson-Durbin（L-D）算法。Yule-Walker（Y-W）方程表明，只要已知输出平稳随机信号的自相关函数，就能求出AR模型中的参数{ak}\left\{a_{k}\right\}{ak​}，并且需要的观测数据较少。

【例1】已知自回归信号模型AR（3）为：

x(n)=1424x(n−1)+924x(n−2)−124x(n−3)+w(n)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）是方差σw2=1{\sigma}_{w}^{2}=1σw2​=1的平稳白噪声，求：

1.自相关序列Rxx(m)R_{x x}(m)Rxx​(m)，m=0，1，2，3，4，5；

2.用上一问求出的自相关序列来估计AR（3）的参数{a^k}\left\{\hat{a}_{k}\right\}{a^k​}，以及输入白噪声的方差σ^w2\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）的参数{a^k}\left\{\hat{a}_{k}\right\}{a^k​}以及输入白噪声的方差σ^w2\hat{\sigma}_{w}^{2}σ^w2​。

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

[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]\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解线性方程组得： R(0)=4.9377R(1)=4.3287R(2)=4.1964R(3)=3.8654\mathrm{R}(0)=4.9377 \quad \mathrm{R}(1)=4.3287 \quad \mathrm{R}(2)=4.1964 \quad \mathrm{R}(3)=3.8654R(0)=4.9377R(1)=4.3287R(2)=4.1964R(3)=3.8654

syms R0 R1 R2 R3A = [1-14/24-9/241/24 ];B = [1000];R = [ R0 R1 R2 R3R1 R0 R1 R2R2 R1 R0 R1R3 R2 R1 R0];[R0 ,R1, R2, R3]=solve(R*A==B);Rxx=vpa([R0 ,R1, R2, R3],5)

利用式（1.6） Rxx(m)=−∑k=1pakRxx(m−k)m>0,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, 可以求出 R(4),R(5)⋯\mathrm{R}(4), \mathrm{R}(5) \cdotsR(4),R(5)⋯

for m = 5 : 6Rxx(m) = 0;for k = 1 : 3Rxx(m) = Rxx(m) - A(k+1) * Rxx(m - k);endendRxx=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 0Rxx(2) Rxx(1) Rxx(2) Rxx(3)Rxx(3) Rxx(2) Rxx(1) Rxx(2)Rxx(4) Rxx(3) Rxx(2) Rxx(1)];B = [ 1000];A = R\BR(1,:) = [Rxx(1) Rxx(2) Rxx(3) Rxx(4)];sigma = R(1,:) * A

结果

A =1.0000-0.5833-0.37510.0417sigma =1.0000

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

离散序列的自相关公式为：

Rxx(m)=1N−m∑n=1N−mXnXn+m,m=0,1,……MR_{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.29411.0475 0.43190.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.786816 至 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.008931 至 32 列-0.0143 -0.0083

把头 4 个相关序列值代入公式(1.7)，利用第二问的代码求得估计值：

a^1=−0.6984a^2=−0.2748a^3=0.0915,σ^w2=0.4678\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 模型参数误差为：e1=0.1151,e2=0.1002,e3=0.0498e_{1}=0.1151, \quad e_{2}=0.1002, \quad e_{3}=0.0498e1​=0.1151,e2​=0.1002,e3​=0.0498，原因在于我们只有一部分的观测数据，使得自相关序列值与理想的完全不同，输入信号的方差误差比较大：eσ=0.5322e_{\sigma}=0.5322eσ​=0.5322。给出一段观测值求 AR 模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

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

AR模型和差分预测的关系

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

前向预测误差系统：

e(n)=x(n)−x^(n)=x(n)+∑k=1mam(k)x(n−k)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)−∑k=1pakx(n−k)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)w(n)=e(n)，且预测系数和 AR 模型参数相同，则两者的系统函数互为倒数， 如图 2 所示：

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

所以求 AR 模型参数就可以通过求预测误差系统的预测系数来实现，对(1.6)式改写：

Rxx(l)=−∑k=1mam(k)Rxx(l−k)l=1,2,⋯m(2.3)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)

Em[e2(n)]=Rxx(0)+∑k=1mam(k)Rxx(k)或Ep[e2(n)]=Rxx(0)+∑k=1pakRxx(k)其中ak=am(k)m=pm是预测器的阶数(2.4)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$ 阶 ARARAR 模型的参数，由于 e(n)=w(n)e(n)=w(n)e(n)=w(n), 所以最小均方预测误差等于白噪声方差，即 Ep[e2(n)]=σw2E_{p}\left[e^{2}(n)\right]=\sigma_{w}^{2}Ep​[e2(n)]=σw2​。

L-D算法基本思想

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

推导过程

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

m=1:

{R(1)=−a1(1)R(0)(1)E1=R(0)+a1(1)R(1)(2)⟶σw12=E1=R(0)[1−a12(1)](2.5)\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：

{R(1)=−a2(1)R(0)−a2(2)R(1)R(2)=−a2(1)R(1)−a2(2)R(0)\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）代入上式得到：

{a2(1)=−R(1)−a2(2)R(1)R(0)=a1(1)+a2(2)a1(1)(1)a2(2)=−R(2)+a2(1)R(1)R(0)=−R(2)+[a1(1)+a2(2)a1(1)]R(1)R(0)(2)⇒a2(2)=−R(2)+a1(1)R(1)R(0)+a1(1)R(1)=−R(2)+a1(1)R(1)E1(3)(2.6)\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)代人:

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)\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)

这样递推下去可以得到预测系数和均方误差估计的通式。

结论

预测系数和均方误差估计的通式：

{am(k)=am−1(k)+am(m)am−1(m−k)(1)am(m)=−R(m)+∑k=1m−1am−1(k)R(m−k)Em−1(2)Em=σwm2=[1−am2(m)]Em−1=R(0)∏k=1m[1−ak2(k)](3)(2.8)\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)

其中 am(m)a_{m}(m)am​(m) 称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值E0=R(0),a0(0)=1,E_{0}=R(0), \quad a_{0}(0)=1, \quadE0​=R(0),a0​(0)=1, 以及 AR\mathrm{AR}AR 模型的阶数 p,p,p,。

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

图3 L-D 算法流程图

L-D 算法的优缺点

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

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

例题

【例2】已知自回归信号模型AR（3）为：

x(n)=1424x(n−1)+924x(n−2)−124x(n−3)+w(n)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）是方差σw2=1{\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）的参数{a^k}\left\{\hat{a}_{k}\right\}{a^k​}以及输入白噪声的方差σ^w2\hat{\sigma}_{w}^{2}σ^w2​。

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

步骤一：计算自相关函数（同例1.3）Rxx(m)=[1.92711.66181.53811.3545⋯]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⋯]​

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

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

m=1:{a1(1)=−R(1)E0=−1.66181.9271=−0.8623E1=R(0)[1−a12(1)]=0.4942m=2:{a2(2)=−R(2)+a1(1)R(1)E1=−0.2127a2(1)=a1(1)[1+a2(2)]=−0.6789E2=E1[1−a22(2)]=0.4718m=3:{a3(3)=−R(3)+a2(1)R(2)+a2(2)R(1)E2=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\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=3p=3p=3 时, 预估到的 AR\mathrm{AR}AR 模型参数为 :a^1=−0.6983,a^2=−0.2748,a^3=0.0914: \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。估计的输入信号的方差为: σ^w2=E3=0.4679\hat{\sigma}_{w}^{2}=E_{3}=0.4679σ^w2​=E3​=0.4679 。和例题 111 的 c\mathrm{c}c 结果一致，误差分析也一样，假设满足精度，停止递推。

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

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

结果

a =1.0000 -0.6984 -0.2748 0.0915E =0.4678

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

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

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

Burg 算法

基本思想:

对观测的数据进行前向和后向预测，然后让两者的均方误差之和为最小作为估计的准则估计处反射系数am(m)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.1759E =0.4426

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

Marple 算法

基本思想:

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

关键在于阶数的选择：

最终预测误差（FPE:final predidyion error）准则定义：给定观测长度为N，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差，该误差在某个阶数p时为最小。其表达式为：

FPE(p)=σ^wp2(N+P+1N−P−1)(3.1)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 的范围内较好。