随机信号的参数建模法及MATLAB实现

随机信号的参数建模法及MATLAB实现

为随机信号建立参数模型是研究随机信号的一种基本方法。在对语音信号进行编码时，往往通过分析不同种类语音信号的特点及产生，用数学模型表示信源，而编码器根据输入信号计算模型参数，然后对模型参数进行编码，也就是说，只需要对编码后的参数进行传送（而不需要传送语音信号本身），解码器通过收到的模型参数，直接利用相同的数学模型即可重建出语音信号，大大减小了传送的数据量。本文将主要介绍AR模型的参数估计以及L-D算法。

一. 概述

在随机信号的参数模型中，我们认为随机信号$x(n)$是由白噪声$w(n)$激励某一确定系统的响应，只要$w(n)$的参数确定了，研究随机信号就可以转化成研究产生随机信号的系统。

图1-1 随机信号的参数模型

对于平稳随机信号，主要有三种常用的线性模型：AR（Auto-Regression，自回归）模型、MA（Moving Average，滑动平均）模型和ARMA（Auto-Regression-Moving Average，自回归滑动平均）模型。

1. MA模型

随机信号$x(n)$由当前的激励$w(n)$和若干次过去的激励$w(n-k)$线性组合产生：
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)\text{\hspace{1em}(1-1)}$$
两边进行$z$变换，可得系统函数：
$$H(z)=\frac{X(z)}{W(z)}=\sum _{k=0}^q{b}_k{z}^{-k}\text{\hspace{1em}(1-2)}$$
其中，$q$表示系统阶数。

可以看出，MA模型的系统函数只有零点，没有极点，所以该系统一定是稳定的系统，也称为全零点模型，表示为$\mathrm{M}\mathrm{A}(q)$。

2. AR模型

随机信号$x(n)$由本身的若干次过去值$x(n-k)$和当前的激励值$w(n)$线性组合产生：
$$x(n)=w(n)-\sum _{k=1}^p{a}_{k}x(n-k)\text{\hspace{1em}(1-3)}$$
该模型的系统函数为：
$$H(z)=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}\text{\hspace{1em}(1-4)}$$
其中，$p$为系统阶数。

该模型系统函数中只有极点，无零点，也称为全极点模型。系统由于极点的原因，要考虑到系统的稳定性，因而要注意极点的分布位置。可表示为$\mathrm{A}\mathrm{R}(p)$。

3. ARMA模型

ARMA模型是AR模型和MA模型的结合：
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)-\sum _{k=1}^p{a}_{k}x(n-k)\text{\hspace{1em}(1-5)}$$
其系统函数为：
$$H(z)=\frac{{\sum }_{k=0}^q{b}_k{z}^{-k}}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}\text{\hspace{1em}(1-6)}$$
可以看出，它既有零点又有极点，所以也称极零点模型，要考虑极零点的分布位置，保证系统的稳定。可表示为$\mathrm{A}\mathrm{R}\mathrm{M}\mathrm{R}(p,q)$。

二. AR模型的参数估计

以上三种模型实际上都可以进行相互转换，实际中选用哪一种模型就要考虑到节约和计算量，剩下的任务就是用适当的算法估计模型参数（$a_k,b_k,p,q$），以便用模型对随机信号进行预测。

实际应用中，使用较多的是AR模型，下面我们对AR模型进行参数估计，其它两种模型可以进行类比计算。

1. AR模型参数与自相关函数的关系

推导过程

对式$(1-3)$两边同时乘以$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]\text{\hspace{1em}(2-1)}$$
根据自相关函数的偶对称特性，可进一步化简为：
$$R_{xx}(m)=R_{xw}(m)-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)\text{\hspace{1em}(2-2)}$$
根据
$$x(n)=w(n)\ast h(n)=\sum _{k=0}^{\infty }h(k)w(n-k)\text{\hspace{1em}(2-3)}$$
有：
$$\begin{array}{rl}{R}_{xw}(m)& =E[x(n)w(n+m)]\\ & =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_{ww}(m+k)\\ & =\sum _{k=0}^{\infty }h(k){\sigma }_w^2\delta (m+k)\\ & ={\sigma }_w^{2}h(-m)\end{array}\text{\hspace{1em}(2-4)}$$
由于系统的单位脉冲响应$h(n)$因果，故
$$R_{xw}(m)=\{\begin{array}{cc}{\sigma }_w^{2}h(-m)& ,m\le 0\\ 0& ,m>0\end{array}\text{\hspace{1em}(2-5)}$$
带入式$(2-2)$得：
$$R_{xx}(m)=\{\begin{array}{ll}{R}_{xx}(-m)& ,m<0\\ -{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)+h(0){\sigma }_w^2& ,m=0\\ -{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)& ,m>0\end{array}\text{\hspace{1em}(2-6)}$$
将式$(1-4)$变换到时域得：$h(n)+{\sum }_{k=1}^p{a}_{k}h(n-k)=\delta (n)$，故当$n=0$时有$h(0)=\delta (0)=1$。

显然，结合式$(2-2)$和式$(2-5)$可得出AR模型输出信号的自相关函数的递推性质：
$$R_{xx}(m)=-\sum _{k=1}^p{a}_k{R}_{xx}(m-k),m>0\text{\hspace{1em}(2-7)}$$
这就是著名的Yule－Walker（Y-W）方程，并将Y-W方程进行变换：
$$R_{xx}(m)+\sum _{k=1}^p{a}_k{R}_{xx}(m-k)=0,m>0\text{\hspace{1em}(2-8)}$$
由$(2-6)$可求得输入的白噪声方差为
$$R_{xx}(0)+\sum _{k=1}^p{a}_k{R}_{xx}(-k)={\sigma }_w^2,m=0\text{\hspace{1em}(2-9)}$$
再将式$(2-8)$和$(2-9)$结合，把该式的下标简化并写成矩阵的形式，可写成单一的正规矩阵方程：
$$\left[\begin{array}{cccc}R(0)& R(-1)& \cdots & R(-p)\\ R(1)& R(0)& \cdots & R(-p+1)\\ ⋮& ⋮& \ddots & ⋮\\ R(p)& R(p-1)& \cdots & R(0)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}{\sigma }_w^2\\ 0\\ ⋮\\ 0\end{array}\right]\text{\hspace{1em}(2-10)}$$
其中方程组的系数组成自相关矩阵$[\mathbf{R}{]}_{p+1}$，由于自相关函数是偶对称函数，即$R_{xx}(m)=R_{xx}(-m)$，因而自相关矩阵为对称矩阵，与主对角线平行的斜对角线上的元素都是相同的，故存在高效算法，其中广泛应用的有Levinson－Durbin（L-D）算法。Y-W方程表明，只要已知输出平稳随机信号的自相关函数，就能求出AR模型中的参数 { a k } \left\{a_k\right\} {ak​}，并且需要的观测数据较少。

实例

【例1】已知自回归信号模型$\mathrm{A}\mathrm{R}(3)$为：
$$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)\text{\hspace{1em}(2-11)}$$
其中$w(n)$是方差${\sigma }_w^2=1$的平稳白噪声，求：

1. 自相关序列$R_{xx}(m),m=0,1,2,3,4,5$；
2. 用上一问求出的自相关序列来估计$\mathrm{A}\mathrm{R}(3)$的参数 { a ^ k } \left\{\hat a_k\right\} {a^k​}，以及输入白噪声的方差${\hat{\sigma }}_w^2$；
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]
   并用观测值的自相关序列直接来估计$\mathrm{A}\mathrm{R}(3)$的参数 { a ^ k } \left\{\hat a_k\right\} {a^k​}以及输入白噪声的方差${\hat{\sigma }}_w^2$。

【解】

1. 由式$(2-11)$可知模型参数 { a k } \left\{a_k\right\} {ak​}，即$a_1=-\frac{14}{24},a_2=-\frac{9}{24},a_3=\frac{1}{24}$，我们只需要代入$(2-10)$并利用其偶对称特性：
   $$\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\\ a_1\\ a_2\\ a_3\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]\text{\hspace{1em}(2-12)}$$
   并重新整理，即可求出$R_{xx}(m)$：
   $$\left[\begin{array}{cccc}1& a_1& a_2& a_3\\ a_1& 1+a_2& a_3& 0\\ a_2& a_1+a_3& 1& 0\\ a_3& a_2& a_1& 1\end{array}\right]\left[\begin{array}{c}R(0)\\ R(1)\\ R(2)\\ R(3)\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]\text{\hspace{1em}(2-13)}$$
   对应的MATLAB程序为：

   a = [-14 / 24, -9 / 24, 1 / 24];
   A = [1, a(1), a(2), a(3);...
        a(1), 1 + a(2), a(3), 0;...
        a(2), a(1)+a(3), 1, 0;...
        a(3), a(2), a(1), 1];
   b = [1; 0; 0; 0];
   Rxx = A \ b
   

   运行结果：

   Rxx =
   
       4.9377
       4.3287
       4.1964
       3.8654
   

   对于$R(4),R(5),\cdots$，可以根据式$(2-7)$：$R_{xx}(m)=-{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k),m>0$求出：

   for m = 5 : 6
       Rxx(m) = 0;
       for k = 1 : 3
           Rxx(m) = Rxx(m) - a(k) * Rxx(m - k);
       end
   end
   Rxx
   

   运行结果：

   Rxx =
   
       4.9377
       4.3287
       4.1964
       3.8654
       3.6481
       3.4027
   

2. 现在已知$R(0),R(1),R(2),R(3)$，则只需要代入到$(2-12)$，即可反求出${\hat{a}}_1,{\hat{a}}_2,{\hat{a}}_3$。相应代码如下：

   R = [Rxx(1) Rxx(2) Rxx(3) Rxx(4);...
        Rxx(2) Rxx(1) Rxx(2) Rxx(3);...
        Rxx(3) Rxx(2) Rxx(1) Rxx(2);...
        Rxx(4) Rxx(3) Rxx(2) Rxx(1)];
   R \ b
   

   结果为：

   ans =
   
       1.0000
      -0.5833
      -0.3750
       0.0417
   

   其中后三个元素为参数 { a ^ k } \left\{\hat a_k\right\} {a^k​}。可以看出对AR模型参数是无失真的估计，因为已知AR模型，我们可以得到完全的输出观测值，因而求得的自相关函数没有失真，当然也就可以不失真地估计。

3. 根据
   $$R_{xx}(m)=\frac{1}{n}\sum _{i=1}^n{x}_i{x}_{i+m}\text{\hspace{1em}(2-14)}$$
   可求出自相关序列$R_{xx}(m),m=0,1,\cdots ,31$。

   代码如下：

   clear; clc;
   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];
   
   Rxx_obs = xcorr(xn) ./ length(xn);
   Rxx_obs = Rxx_obs(length(xn) : end)
   

   结果如下：

 Rxx_obs =
 
   列 1 至 8
 
     1.9271    1.6618    1.5381    1.3545    1.1349    0.9060    0.8673    0.7520
 
   列 9 至 16
 
     0.7637    0.8058    0.8497    0.8761    0.9608    0.8859    0.7868    0.7445
 
   列 17 至 24
 
     0.6830    0.5808    0.5622    0.5134    0.4301    0.3998    0.3050    0.2550
 
   列 25 至 32
 
     0.1997    0.1282    0.0637    0.0329   -0.0015   -0.0089   -0.0143   -0.0083

将$R_{xx}(0),R_{xx}(1),R_{xx}(2),R_{xx}(3)$代入$(2-12)$可求得估计值：${\hat{a}}_1=-0.6984,{\hat{a}}_2=-0.2748,{\hat{a}}_3=0.0915,{\hat{a}}_4=0.4678$，与真实AR模型参数误差为：$e_1=0.1151,e_2=0.1002,e_3=0.0498$，原因在于我们只有一部分的观测数据，使得自相关序列值有偏差。输入信号的方差误差比较大：$e_{\sigma }＝0.5322$，造成的原因比较多，如计算机仿真的白噪声由于只有32点长，其方差不可能刚好等于1。给出一段观测值求AR模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法减小计算量并提高精度。

2. Yule－Walker方程的解法——Levinson－Durbin算法

准备工作

在上面的例子中，为了得到更精确的估计值并减小计算量，就要建立更高阶的AR模型。因此把 AR 模型和预测系统联系起来，换个方法来估计 AR 参数。

由式$(1-3)$：$x(n)=w(n)-{\sum }_{k=1}^p{a}_{k}x(n-k)$知，模型的当前输出值与它过去的输出值有关。预测是推断一个给定序列的未来值，即利用信号前后的相关性来估计未来的信号值。

若序列的模型已知，而用过去观测的数据来推求现在和将来的数据的模型称为前向预测器，表示为：
$$\hat{x}(n)=-\sum _{k=1}^m{a}_m(k)x(n-k)\text{\hspace{1em}(2-15)}$$
其中 { a ( k ) } , k = 1 , 2 , ⋯   , m \left\{a(k)\right\},k=1,2,\cdots,m {a(k)},k=1,2,⋯,m表示$m$阶预测器的预测系数，负号是为了与技术文献保持一致。显然预测出来的结果与真实的结果存在预测误差（或前向预测误差），设误差为$e(n)$：
$$e(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^m{a}_m(k)x(n-k)\text{\hspace{1em}(2-16)}$$
把$e(n)$看成是系统的输出，$x(n)$看成是系统的输入，得到系统函数
$$\frac{E(z)}{X(z)}=1+\sum _{k=1}^m{a}_m(k)z^{-k}\text{\hspace{1em}(2-17)}$$
设$m=p$，且预测系数与AR模型参数相同，比较$(1-4)$：$H(z)=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$与上式，可以看出系统函数互为倒数。结合图2-1，则有$w(n)=e(n)$，即前向预测误差系统中的输入为$x(n)$，输出为预测误差$e(n)$等于白噪声。也就是说前向预测误差系统对观测信号起了白化的作用。由于AR模型和前向预测误差系统函数互为倒数，所以求AR模型参数就可以通过求预测误差系统的预测系数来实现。

图2-1 左：前向预测误差系统；右：AR模型

对$(2-16)$求预测误差均方值：
$$\begin{array}{rl}E\left[e^2(n)\right]& =E\left[{\left(x(n)+\sum _{k=1}^m{a}_m(k)x(n-k)\right)}^2\right]\\ & =R_{xx}(0)+2\left[\sum _{k=1}^m{a}_m(k)R_{xx}(k)\right]+\sum _{k=1}^m\sum _{l=1}^m{a}_m(l)a_m(k)R_{xx}(l-k)\end{array}\text{\hspace{1em}(2-18)}$$
要使均方误差最小，将上式右边对预测系数求偏导并且等于零，得到$m$个等式：
$$R_{xx}(l)=-\sum _{k=1}^m{a}_m(k)R_{xx}(l-k)l=1,2,\cdots m\text{\hspace{1em}(2-19)}$$
并代入$(2-18)$，可得最小均方误差：
$$E_m\left[e^2(n)\right]=R_{xx}(0)+\sum _{k=1}^m{a}_m(k)R_{xx}(k)\text{\hspace{1em}(2-20)}$$
由于
$$E_p\left[e^2(n)\right]=R_{xx}(0)+\sum _{k=1}^p{a}_k{R}_{xx}(k)\text{\hspace{1em}(2-21)}$$
则有
$$a_k=a_m(k),m=p\text{\hspace{1em}(2-22)}$$
也就是$p$阶预测器的预测系数等于$p$阶AR模型的参数，由于$e(n)=w(n)$，所以最小均方预测误差等于白噪声方差，即$E_p\left[e^2(n)\right]={\sigma }_w^2$。

使用L-D算法估计AR模型参数

有了上面的知识后，我们再看如何使用L-D模型估计AR模型参数，即如何估计参数$\left\{a_k,p,{\sigma }_w^2\right\}$。L-D算法的基本思想基于Y-W方程式或$(2-19),(2-20),(2-21)$和自相关序列具有递推的性质。L-D递推算法是模型阶数逐渐加大的一种算法，先计算阶次$m＝1$时的预测系数 { a m ( k ) } = a 1 ( 1 ) \left\{ a_m(k)\right\}=a_1(1) {am​(k)}=a1​(1)和${\sigma }_{w1}^2$，再计算$m＝2$时的预测系数 { a m ( k ) } = a 2 ( 1 ) ,   a 2 ( 2 ) \left\{ a_m(k)\right\}=a_2(1),\ a_2(2) {am​(k)}=a2​(1), a2​(2)和${\sigma }_{w2}^2$，一直计算到$m=p$阶时的$a_p(1),a_p(2),\cdots ,a_p(p)$及${\sigma }_{wp}^2$。这种递推算法的特点是，每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。

经过一系列推到，最终可以得到：
$$\{\begin{array}{l}{a}_m(k)=a_{m-1}(k)+a_m(m)a_{m-1}(m-k)\\ a_m(m)=-\frac{R(m)+{\sum }_{k=1}^{m-1}{a}_{m-1}(k)R(m-k)}{E_{m-1}}\\ E_m={\sigma }_{wm}^2=\left[1-a_m^2(m)\right]E_{m-1}=R(0){\prod }_{k=1}^m\left[1-a_k^2(k)\right]\end{array}\text{\hspace{1em}(2-23)}$$
其中，$a_m(m)$称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值$E_0=R(0),a_0(0)=1$以及AR模型的阶数$p$，就可以按照图2-2所示的流程图进行估计。

图2-2 Levinson－Durbin算法流程图

L-D算法的特点

L-D算法的优点就是计算速度快，求得的AR模型必定稳定，且均方预测误差随着阶次的增加而减小。

由于使用L-D算法在求自相关序列时，我们假设除了观测值之外的数据都为零，因此必然会引入较大误差。

实例

【例2】已知自回归信号模型$\mathrm{A}\mathrm{R}(3)$为：
$$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)\text{\hspace{1em}(2-11)}$$
其中$w(n)$是方差${\sigma }_w^2=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]，用L-D算法来估计$\mathrm{A}\mathrm{R}(3)$的参数 { a ^ k } \left\{\hat a_k\right\} {a^k​}以及输入白噪声的方差${\hat{\sigma }}_w^2$。

【解】

• 步骤1：求$R_{xx}(m),m=0,1,\cdots ,31$，这部分计算的方法与例1中第3问完全相同，这里不再赘述。

• 步骤2：初始化：$E_0=R_{xx}(0)=1.9271,a_0=1$

• 步骤3：根据$(2-23)$有：
  $$\begin{array}{rl}& m=1:\{\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}\\ & \\ & m=2:\{\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}\\ & \\ & m=3:\{\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}\end{array}$$
  因而当$p=3$时，估计得到的AR模型参数为：${\hat{a}}_1=-0.6983,{\hat{a}}_2=-0.2748,{\hat{a}}_3=0.0914$，估计的输入信号方差为${\sigma }_w^2=E_3=0.4679$，这与例1的第3问结果一致。

在MATLAB中，有专门的函数[a E] = aryule(x, p)来实现L-D算法的AR模型参数估计，其中x表示观测信号，p表示阶数，a表示估计的模型参数，E表示噪声信号的方差估计。本例的代码为：

[a E] = aryule(xn, 3)

运行结果为：

a =

 1.0000   -0.6984   -0.2748    0.0915


E =

 0.4678

我们也可以使用更高阶的AR模型来进行估计，如$p=12$：

a =

  列 1 至 8

    1.0000   -0.6703   -0.3254   -0.0793    0.1407    0.3676   -0.2451    0.0483

  列 9 至 13

   -0.0912   -0.0522    0.0515    0.0186   -0.0955


E =

    0.3783

可见阶数越高，均方误差就越小（$0.3783<0.4678$）。

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

Burg算法

为了克服L-D算法导致的误差，1968年Burg提出了Burg算法，其基本思想是对观测的数据进行前向和后向预测，然后让两者的均方误差之和为最小作为估计的准则估计处反射系数，进而通过L-D算法的递推公式求出AR模型参数。Burg算法的优点是，求得的AR模型是稳定的，较高的计算效率，但递推还是用的L-D算法，因此仍然存在明显的缺点。

MATLAB中也有专门的函数来实现Burg算法的AR模型参数估计：[a E] = arburg(x, p)，例如上面例2在$p=3$时和$p=12$时的计算结果分别为

a =

    1.0000   -0.6982   -0.2626    0.0739


E =

    0.4567

和

a =

  列 1 至 8

    1.0000   -0.6495   -0.3066   -0.0934    0.0987    0.4076   -0.1786   -0.0126

  列 9 至 13

   -0.0805   -0.0899    0.0382    0.1628   -0.2501


E =

    0.3237

结果与L-D算法略有不同，Burg算法略微更精确些。

Marple算法

1980年Marple在前人的基础上提出一种高效算法，Marple算法或者称不受约束的最小二乘法（LS），该算法的思想是，让每一个预测系数的确定直接与前向、后向预测的总的平方误差最小，这样预测系数就不能由低一阶的系数递推确定了，所以不能用L-D算法求解，实践表明该算法比L-D、Burg算法优越。由于该算法是从整体上选择所有的模型参数达到总的均方误差最小，与自适应算法类似，不足是该算法不能保证AR模型的稳定性。

AR模型阶数的选择

AR模型的阶数选择对模型效果有有较大的影响，因而如何选择阶数很重要。因此国内外学者在这方面都做了许多研究工作，其中基于均方误差最小的最终预测误差（FPE）准则是确定AR模型阶次比较有效的准则。

最终预测误差准则定义如下：给定观测长度为$N$，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差，该误差在某个阶数$p$时为最小，其表达式为：
$$\mathrm{F}\mathrm{P}\mathrm{E}(p)={\hat{\sigma }}_{wp}^2\left(\frac{N+p+1}{n-p-1}\right)\text{\hspace{1em}(7-24)}$$
上式中估计的方差随着阶数的增加而减小，而括号内的值随着$p$的增加而增加，因而一定存在一个最佳的$p$，使得FPE最小。

我们还可以使用试验的方法来进行阶数的选择。如果某一个$p$值满足预先规定的某些要求，则认为该$p$值是最优的。在短数据情况下，根据经验，AR模型的阶次选在$\frac{N}{3}\sim \frac{N}{2}$的范围内较好。