随机信号的参数建模法--AR模型及Matlab实现_matlab-ar函数

一、随机信号参数模型

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

经典信号建模法（classical modeling method for signal）前面已经指出，医学信号处理的目的是提取包含于随机信号中的确定性成分，以便在一定的准确性（最小二乘意义）上进行预测。这就是建立各种各样的确定性数学模型，包括代数、微分、积分、差分方程模型。这是经典的信号建模方法。
信号的现代建模方法（Modern modeling method for signal）是建立在具有最大的不确定性基础上的预测。提出了众多的数学模型( mathematical models)。根据 Wold 的证明：任何平稳的 ARMA（自回归移动平均）模型或 MA 模型均可用无限阶或阶数足够的 AR 模型去近似。因此本节着重介绍 AR 模型的基本原理和方法。
对平稳随机信号，三种常用的线性模型分别是 AR 模型（自回归模型 Auto-regression model），MA 模型（滑动平均模型 Moving average model）和 ARMA 模型（自回归滑移平均模型 Auto-regression-Moving average model）。

这里对AR模型展开讨论。

二、AR模型

随机信号$x(n)$由本身的若干次过去值 $x(n-k)$ 和当前的激励值$w(n)$ 线性组合产生：

该模型的系统函数是：

$p$是系统阶数，系统函数中只有极点，无零点，也称为全极点模型，系统由于极点的原因，要考虑到系统的稳定性，因而要注意极点的分布位置，用$AR(p)$来表示。

三、AR模型参数的估计

随机信号的建模法最近在生物医学信号处理中应用相当普遍，在自发脑电、诱发脑电、肌电、心电、胃电等方面都有人尝试应用模型法进行研究。应用较多的是 AR 模型，因为建立这种的模型计算工作比较容易。作为数学逼进，三种模型都可以互相转换。实际中选用哪一种模型就要考虑到节约和计算量，选定模型后，剩下的任务就是用适当的算法估计模型参数$(a_k、b_k、p、q)$，以便用模型对随机信号进行预测。
下面我们就以 AR 模型为例进行参数的估计。

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

根据图 7.5 和式(7-24)，对该式两边同时乘以 $x(n-m)$，然后求均值：

因为自相关函数：

所以自相关函数呈现偶对称，(7-27)式化为：

系统的单位脉冲响应 是因果的，所以输出的平稳随机信号和输入的白噪声之间的互相
关函数有下列推导：

所以

带入式(7-28)得到：

由于$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(0)=1$,显然，AR 模型输出信号的自相关函数具有递推的性质，即:

上式就是著名的 Yule－Walker（Y-W）方程，将上式变换：

从（7-33）求得输入的白噪声方差为：

将(7-35)和(7-36)结合，把该式的下标简化并写成矩阵的形式，可以写成单一的正规矩阵方程：

式（7-37）方程组的系数都是自相关矩阵 ,由于自相关函数是偶对称函数：$R_{xx}(m)=R_{xx}(-m)$，因而自相关矩阵是对称矩阵，与主对角线平行的斜对角线上的元素都是相同的，是 的维托毕兹（Toeplitz）矩阵，所以存在高效算法，其中应用广泛的有 Levinson－Durbin（L-D）算法。Yule－Walker（Y-W）方程表明，只要已知输出平稳随机信号的自相关函数，就能求出 AR 模型中的参数 { a k } {\{a_k }\} {ak​}，并且需要的观测数据较少。

例1 已知自回归信号模型AR(3)为：

上述模型的Matlab实现如下，具体代码含义已标注。用fsolve解函数，其中序列x_T=[$a_1$ $a_2$ $a_3$ ${\sigma }_w^2$]

clc;
a=[-14/24 -9/24 1/24];%模型参数{ak}
b=[1;zeros(3,1)];
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];
Rx_T=A\b  %题a.真实的自相关序列
R_T=[Rx_T(1),Rx_T(2),Rx_T(3),Rx_T(4);
    Rx_T(2),Rx_T(1),Rx_T(2),Rx_T(3);
    Rx_T(3),Rx_T(2),Rx_T(1),Rx_T(2);
    Rx_T(4),Rx_T(3),Rx_T(2),Rx_T(1)];
%题b.参数ak及方差
x_T=fsolve(@(x_T)R_T*[1 x_T(1) x_T(2) x_T(3)].'-[x_T(4) 0 0 0].',rand(1,4))

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];
%自相关估计值 
Rx=xcorr(x_n)./length(x_n);
Rx=Rx(length(x_n):end)
R=[Rx(1),Rx(2),Rx(3),Rx(4);
    Rx(2),Rx(1),Rx(2),Rx(3);
    Rx(3),Rx(2),Rx(1),Rx(2);
    Rx(4),Rx(3),Rx(2),Rx(1)];
%题c.估计的参数ak及方差
x=fsolve(@(x)R*[1 x(1) x(2) x(3)].'-[x(4) 0 0 0].',rand(1,4))

%误差
e=abs(x-x_T)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

结果：

Rx_T =

    4.9377
    4.3287
    4.1964
    3.8654

x_T =

   -0.5833   -0.3750    0.0417    1.0000
   
Rx =

  1 至 14 列

    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

  15 至 28 列

    0.7868    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

  29 至 32 列

   -0.0015   -0.0089   -0.0143   -0.0083

x =

   -0.6984   -0.2748    0.0915    0.4678

e =

    0.1150    0.1002    0.0498    0.5322
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

b结果分析：

$a_1$=-0.5833 (即-14/24) $a_2$=-0.3750 (即-9/24)$a_3$= 0.0417 (即1/24) ${\sigma }_w^2$=1

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

c结果分析：

求出估计值为$a_1$= -0.6984 $a_2$=-0.2748 $a_3$= 0.0915 ${\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. Y－W 方程的解法——L-D 算法

在求解【例 7-2】时要得到更精确的估计值，就要建立更高阶的 AR 模型，直接用观测值的自相关序列来求解 Y-W 方程计算量太大。因此把 AR 模型和预测系统联系起来，换个方法来估计 AR 参数。

2-1 前向预测器

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

若序列的模型已知而用过去观测的数据来推求现在和将来的数据称为前向预测器，表示为：

式中${a(k)}$，k＝1，2，…，m，代表 m 阶预测器的预测系数，负号是为了与技术文献保持一致。显然预测出来的结果与真实的结果存在预测误差或前向预测误差，设误差为$e(n)$:

把$e(n)$看成是系统的输出，$x(n)$看成是系统的输入，得到系统函数：

假如$m＝p$，且预测系数和 AR 模型参数相同，比较式(7-41)和式(7-24)，把预测误差系统框图和 AR 模型框图给出，如图 7.6 所示，即有 $w(n)=e(n)$，即前向预测误差系统中的输入为 $x(n)$，输出为预测误差$e(n)$等于白噪声。也就是说前向预测误差系统对观测信号起了白化的作用。由于 AR 模型和前向预测误差系统有着密切的关系，两者的系统函数互为倒数，所以求 AR 模型参数就可以通过求预测误差系统的预测系数来实现。

对(7-40)式求预测误差均方值：

要使得均方误差最小，将上式右边对预测系数求偏导并且等于零，得到 m 个等式：

将(7-43)代入(7-42)求得最小均方误差：

或

也就是 p 阶预测器的预测系数等于 p 阶 AR 模型的参数，由于 $w(n)=e(n)$，所以最小均方预测误差等于白噪声方差，即 $E_p[e^2(n)]={\sigma }_w^2$。

2-2 建立更高阶的AR模型

有了上面的知识后，我们回来看怎样估计 AR 模型参数，也即要估计参数{$a_k,p,{\sigma }_w^2$}，这里介绍应用广泛的L-D算法。L-D 算法的基本思想就是根据Y-W方程式或式(7-43)、(7-44)、(7-45)，自相关序列具有递推的性质，L-D 递推算法是模型阶数逐渐加大的一种算法，先计算阶次 m＝1 时的预测系数{$a_m(k)$}=$a_1(1)$和${\sigma }_{w1}^2$ ，然后计算 m＝2 时的{$a_m(k)$}=$a_2(1)$和${\sigma }_{w2}^2$ ,一直计算到 m＝p 阶时的$a_p(1),a_p(2)...,a_p(p)$ 以及${\sigma }_{wp}^2$ 。这种递
推算法的特点是，每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。

经过推导，得到预测系数和均方误差估计的通式：

流程图：

L-D 算法的优点就是计算速度快，求得的 AR 模型必定稳定，且均方预测误差随着阶次的增加而减小（见(7-50)的(3)式）。L-D 算法的缺点，由于在求自相关序列时，是假设除了观测值之外的数据都为零，必然会引入较大误差。

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

可以套用公式(7-50)进行计算，也可以在Matlab 使用专门的函数实现 L-D 算法的 AR 模型参数估计：[a E]＝aryule（x，p），输入 x表示观测信号，输入 p 表示要求的阶数，输出 a 表示估计的模型参数，e 表示噪声信号的方差估计。例如本题用该函数计算结果为：

clc;
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];
[a E] = aryule(x_n, 3)
1
2
3

a =
    1.0000   -0.6984   -0.2748    0.0915
E =
    0.4678
1
2
3
4

用更高阶的AR模型估计，p=12：

a =
    1.0000   -0.6703   -0.3254   -0.0793    0.1407    0.3676   -0.2451    0.0483   -0.0912   -0.0522    0.0515    0.0186   -0.0955
E =
    0.3783
1
2
3
4

由此发现，阶数越高，均方误差越小 。

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

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

clc;
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];
[a E] = arburg(x_n, 3)
1
2
3

结果：与L-D算法结果略有不同，$a_3$估计得更精确。

a =
    1.0000   -0.6982   -0.2626    0.0739
E =
    0.4567
1
2
3
4

阶数12：

a =
    1.0000   -0.6495   -0.3066   -0.0934    0.0987    0.4076   -0.1786   -0.0126   -0.0805   -0.0899    0.0382    0.1628   -0.2501
E =
    0.3237
1
2
3
4

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

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

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

上式中估计的方差随着阶数的增加而减小，而括号内的值随着 的增加而增加，因而能找到一最佳的$p_{opt}$，使得 FPE 最小。