数据压缩（十四）——AR模型的参数估计阅读

一、AR模型

  $AR$模型是随机信号参数模型的一种。

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

二、AR模型参数的估计

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

2.1.1 推导过程

2.1.2 举例说明

【例】已知AR模型AR(3)为：$$x(n)=\frac{14}{24}x(n-1)+\frac{9}{24}x(n-2)-\frac{1}{24}x(n-3)+w(n)$$  式中$w(n)$是具有方差${\sigma }_w^2=1$的平稳白噪声。
【1】求自相关序列$R_{xx}(m),m=0,1,2,3,4,5$。

【2】用a求出的自相关序列来估计AR(3)的参数{${\hat{a}}_k$}，以及输入白噪声的方差${\hat{\sigma }}_w^2$大小。

  与[1]相比，只是已知条件相反而已。

【3】利用给出的AR模型，用计算机仿真给出32点观测值$$x(n)=\left[\begin{array}{cccccccccccccccccccccccccccccccc}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\end{array}\right]$$用观测值的自相关序列直接来估计AR(3)的参数{${\hat{a}}_k$}以及输入的白噪声${\sigma }_w^2$

  把头4个相关序列值带入前两问一直在使用的矩阵方程，求得估计值:$${\hat{a}}_1=-0.6984,{\hat{a}}_2=-0.2748,{\hat{a}}_3=0.0915,{\hat{\sigma }}_w^2=0.4678$$
  可以发现，计算所得的估计值和真实的AR模型参数的误差为：$e_1=0.1151,e_2=0.1002,e_3=0.0498$。其原因在于：我们只有一部分的观测数据，使自相关序列值和理想的值完全不同。
  输入信号的方差误差比较大：$e_{\sigma }=0.5322$，其原因比较多：其中一个是计算机仿真的白噪声只有32点长，32点序列的方差不可能刚好等于1。
  给出一段观测值求AR模型参数这样解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

2.2 Y-W方程的解法——L-D算法

2.2.1 先验知识

  在求解时要得到更精确的估计值，就要建立更高阶的AR模型，直接用观测值的自相关序列来求解Y-W方程的计算量太大。因此把AR模型和预测系统联系起来，换个方法来估计AR参数。
  从AR模型的时域表达式：$x(n)=w(n)-{\Sigma }_{k=1}^p{a}_{k}x(n-k)$，我们知道模型的当前输出值与他过去的输出值有关。预测是推断一个给定序列的未来值，即利用信号前后的相关性来估计未来的信号值。
  若序列的模型已知而用过去观测的数据来推求现在和将来的数据称为前向预测器，表示$$\hat{x}(n)=-{\Sigma }_{k=1}^m{a}_m(k)x(n-k)$$
  显然预测出来的结果与真实的结果存在预测误差或前向预测误差，设误差为$e(n)$:
$$e(n)=x(n)-\hat{x}(n)=x(n)+{\Sigma }_{k=1}^m{a}_m(k)x(n-k)$$
  把$e(n)$看成是系统的输出，$x(n)$看成是系统的输入，得到系统函数：$$\frac{E(z)}{X(z)}=1+{\Sigma }_{k=1}^m{a}_m(k)z^{-k}$$
  假如$m=p$,且预测系数和AR模型参数相同，比较$\frac{E(z)}{X(z)}=1+{\Sigma }_{k=1}^m{a}_m(k)z^{-k}$和$x(n)=w(n)-{\Sigma }_{k=1}^p{a}_{k}x(n-k)$，把预测误差系统框图和AR模型框图给出

  即有$w(n)=e(n)$,即前向预测误差系统中的输入为$x(n)$,输出预测误差$e(n)$等于白噪声。也就是说前向预测误差系统对观测信号起了白化的作用。由于AR模型和前向预测误差系统有着密切的关系，两者的系统函数互为倒数，所以求AR模型参数就可以通过求预测误差系统的预测系数来实现。

2.2.2 L-D算法

  但是按照通式计算太为繁琐，所以直接按照递推的思路绘制流程图，使用代码来帮助解决问题。流程图如下：

L-D算法的优点：

• 计算速度快
• 求得的AR模型必定稳定
• 且均方预测误差随着阶次的增加而减小

L-D算法的缺点：

• 有较大误差【因为在求自相关序列时，假设除了观测值之外的数据均为0】

2.2.3 L-D算法Matlab实现

值得注意的是：matlab从1开始计数，所以需要对模型的参数进行一定的调整，将0变为1，并且往后所有的下标都向后挪动一位。

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];

R_xx = zeros( 1, length(x_n) );

for m  = 0 : length(x_n) - 1 
    j = 1; 
    while j + m <= length(x_n)
        R_xx(m+1) = R_xx(m+1) + x_n(j) * x_n(j + m );
        j = j + 1;
    end
    R_xx(m+1) =  R_xx(m+1) / 32;
end

%初始化
p = 4;
E = zeros(1, p); E(1) = R_xx(1);
a = zeros(p, 3); a(1 , 1) = 1;
a_final=zeros(1,3);

for m = 1:3
    %开始计算a(m+1,m)[先计算是因为后面的式子里用到了这个值]
    a( m + 1, m ) = R_xx(m + 1);
    for k = 1:( m - 1 )
        a( m + 1, m ) = a( m + 1, m ) + a( m, k ) * R_xx( m + 1 - k );
    end
    a( m + 1, m ) = - a( m + 1, m ) / E(m);
    
    %开始计算a(m+1,k)
    for k = 1:(m - 1)
       a( m + 1, k ) = a( m, k ) + a( m + 1, m ) * a( m, m - k );
    end 
    
    %开始计算E
    E( m + 1 ) = ( 1 - a( m + 1, m ).^2 ) * E(m); 
end

for k=1:3
    a_final(k)=a(4,k);
end

  最后输出的结果为：

  同时可以调用matlab本身的函数来进行参数估计。

  或者采用更高阶的AR模型来估计：

  可以发现，阶数越高，误差越小。

2.3 对L-D算法的改进

2.3.1 Burg算法

  为了克服L-D算法的误差，Burg提出了Burg算法。其基本思想是对观测的数据进行前向和后向预测，然后让两者的均方误差之和作为估计的准则估计出反射系数，进而通过L-D算法的递推公式求出AR模型参数。
  Burg算法的优点：

• 求得的AR模型稳定
• 有较高的计算效率

  但是由于Burg的递推依然采用的是L-D算法，所以存在明显缺点。

  结果和L-D算法略有不同，${\hat{a}}_3$估计的更精确一点。

2.3.2 Marple算法

2.4 AR模型阶数选择

  基于均方误差最小的最终预测误差（FPE:final predidyion error）准则是确定AR模型阶次比较有效的准则。
  其定义是：给定观测长度N，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差。该误差在某个阶数p时为最小，其表达式为：$$FPE(p)={\hat{\sigma }}_{wp}^2(\frac{N+P+1}{N-P-1})$$  上式中估计的方差随着阶数的增加而减小，而括号内的值随着p的增加而增加，因而能找到一个最佳的$p_{opt}$，使FPE最小。
  除了用这种方法进行阶数的选择，还可以采用试验的方法：如果某一个p值满足预先规定的某些要求，则认为该p值是最优的。在短数据情况下，根据经验，AR模型的阶次选在N/3~N/2的范围内较好。