随机信号参数建模法

在进行语音信号编码传输时，往往对语音信号信源进行建模，然后对模型参数进行编码，只传送编码后的模型参数。

这样在解码端获得模型参数重建信源模型后，即可获得重建后的语音信号。

信号的现代建模法是建立在具有最大的不确定性基础上的预测。其广义定义为：随机信号x(n)由白噪w(n)激励某一确定系统的响应。确定了白噪就，研究随机信号就等价于研究产生随机信号的系统。

对平稳随机信号，三种常用的线性模型为：

• AR模型（自回归模型，Auto-regression model）
• MA模型 （滑动平均模型 Moving average model）
• ARMA模型 （自回归滑动平均模型 Auto-regression-Moving average model）

着重介绍AR模型的基本原理和方法

三种参数模型

MA模型

随机信号x(n)由当前激励w(n)和若干次过去的激励w(n-k)线性组合产生：
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)$$
系统函数为
$$H(z)=\frac{X(z)}{W(z)}=\sum _{k=0}^q{z}^{-k}$$
其中q表示阶数。系统函数只有零点没有极点，为全零点模型，一定是稳定的。全零点用MA(q)表示。

AR模型

随机信号由本身的若干次过去值x(n-k)和当前激励w(n)线性组合产生：
$$x(n)=w(n)-\sum _{k=1}^P{a}_{k}x(n-k)$$
系统函数为:
$$H(z)=\frac{1}{1+{\sum }_{k=1}^P{a}_k{z}^{-k}}$$
其中p为阶数。这是个全极点模型，只有极点没有零点。系统的稳定性取决于极点的分布位置，用AR§表示。

ARMA表示

是AR模型和MA模型的结合,由数个x(n-k)和w(n-k)线性组合形成：
$$x(n)=\sum _{k=0}^q{b}_{k}w(n-k)-\sum _{k=1}^p{a}_{k}x(n-k)$$
系统函数为
$$H(z)=\frac{{\sum }_{k=0}^q{b}_k{z}^{-k}}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$
既有零点又有极点。稳定性由零极点分布确定，表示为ARMA(p,q).

AR模型参数估计

使用AR模型的场合较多，这是因为计算工作比较容易。实际上三种模型都可以相互转化，选择哪一种模型主要取决于计算量和节约。选定模型后下一步就是确定模型参数了。

AR模型和自相关函数

对之前的ZR模型x(n)表达式等式两边同乘x(n-m),然后求均值
$$E[x(n)x(n-m)]=E[w(n)x(n-m)-\sum _{k=1}^P{a}_{k}x(n-k)x(n-m)]$$
有自相关函数:
$$R_{xx}(m)=E[x(n)x(n+m)]=E[x(k-m)x(k)]=E[x(n)x(n-m)]=R_{xx}(-m)$$
呈偶对成，所以原式可化为：
$$R_{xx}(m)=R_{xw}(m)-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)$$
系统的单位脉冲响应h(n)是因果的，有下列推导：
$$\begin{gathered} R_{xw}(m)=E[x(n)w(n+m)] \\ x(N)=w(n)\ast h(n)=\sum _{k=0}^{\infty }h(k)w(n-k) \\ {R}_{xw}(m)=E[\sum _{k=0}^{\infty }h(k)w(n-k)w(n+m)]=h(k)E[\sum _{k=0}^{\infty }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{gathered}$$
最后可得
$$R_{xw}(m)=\{\begin{array}{ll}0& m>0\\ {\sigma }_w^{2}h(-m)& m\le 0\end{array}$$

带入自相关函数可得
$$R_{xx}(m)=\{\begin{array}{ll}-{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)& m>0\\ -{\sum }_{k=1}^p{a}_k{R}_{xx}(m-k)+h(0){\sigma }_w^2& m=0\\ R_{xx}(-m)& m<0\end{array}$$
其中由系统函数可得h(0)=1。

可得AR模型输出信号自相关函数由递推性质：
$$R_{xx}(m)=-\sum _{k=1}^p{a}_k{R}_{xx}(m-k)m>0$$
这也就是Yule-Walker（Y-W）方程，变化后渴求输入白噪方差：
$$\begin{gathered} 0=R_{xx}(m)+\sum _{k=1}^p{a}_k{R}_{xx}(m-k)m>0 \\ sigm{a}_w^2=R_{xx}(0)+\sum _{k=1}^p{a}_k{R}_{xx}(-k)m=0 \end{gathered}$$
写成矩阵形式，可得单一的正规矩阵方程：
$$\left[\begin{array}{cccc}R(0)& R(-1)& ...& R(-p)\\ R(1)& R(0)& ...& R(-p+1)\\ ...& ...& ...& ...\\ R(p)& R(p-1)& ...& 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]$$
第一个矩阵的元素都是自相关函数。由于自相关函数是偶函数，对角线上的元素相同，存在高效算法。

只要已知输出平稳随机信号的自相关函数，就可的AR模型参数，且不需要很多观测数据。

实例

x(n)表达式中的系数即为a1,a2,a3.代入矩阵表达式可得：

$$\left[\begin{array}{cccc}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\\ -\frac{14}{24}\\ -\frac{9}{24}\\ \frac{1}{24}\end{array}\right]=\left[\begin{array}{c}{\sigma }_w^2\\ 0\\ ...\\ 0\end{array}\right]$$

转化可得：
$$\left[\begin{array}{cccc}a(1)& a(2)& a(3)& a(4)\\ a(2)& a(1)+a(3)& a(4)& 0\\ a(3)& a(2)+a(4)& a(1)& 0\\ a(4)& a(3)& a(2)& a(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\end{array}\right]$$
利用matlab解线性方程组：

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

可得R（0)-R（3）：

Rxx =

    4.9377
    4.3287
    4.1964
    3.8654

通过推导公式可得R(4)R(5)：
$$\begin{gathered} R(4)=-\sum _{k=1}^3{a}_k{R}_{xx}(4-k)=3.6418 \\ R(4)=-\sum _{k=1}^3{a}_k{R}_{xx}(5-k)=3.4027 \end{gathered}$$
将R(0)-R(3)带入矩阵表达式
$$\left[\begin{array}{cccc}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\\ {\hat{a}}_1\\ {\hat{a}}_2\\ {\hat{a}}_3\end{array}\right]=\left[\begin{array}{c}{\hat{\sigma }}_w^2\\ 0\\ ...\\ 0\end{array}\right]$$
通过计算可得
$$\begin{gathered} {\hat{a}}_1=-\frac{14}{24} \\ {\hat{a}}_2=-\frac{9}{24} \\ {\hat{a}}_3=-\frac{1}{24} \\ {\hat{\sigma }}_w^2=1 \end{gathered}$$
和题目中参数一样，没有失真，这是因为AR模型已知。

已有一定长的观察值，则可以求自相关函数：

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 = zeros( 1, length(xn) );
for m = 0 : length(xn) - 1
    i = 1;
    while i+m <= length(xn)
        Rxx(m+1) = Rxx(m+1) + xn(i) * xn(i + m);
        i = i + 1;
    end
    Rxx(m+1) = Rxx(m+1) / 32;
end
Rxx
%结果
Rxx =

  1 至 10 列

    1.9271    1.6618    1.5381    1.3545    1.1349    0.9060    0.8673    0.7520    0.7637    0.8058

  11 至 20 列

    0.8497    0.8761    0.9608    0.8859    0.7868    0.7445    0.6830    0.5808    0.5622    0.5134

  21 至 30 列

    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

由此可得Rxx，取前四位，即R(0)-R(3)，带入矩阵，用同样的方法可得估计值：
$$\begin{gathered} {\hat{a}}_1=-0.6984 \\ {\hat{a}}_2=-0.2748 \\ {\hat{a}}_3=0.0915 \\ {\hat{\sigma }}_w^2=0.4678 \end{gathered}$$
与真实AR模型有以下误差：
$$\begin{gathered} e_1=0.1151 \\ {e}_2=0.1001 \\ {e}_3=0.0498 \\ {e}_{\sigma }=0.5322 \end{gathered}$$
这是因为只有一部分观测数值，不能完全理想的还原AR模型,且模拟的白噪声长度也只有32点。

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

为了获得更精确的估计值，就要建立更高阶的AR模型，但同时观测值求解的计算量会增加很多。因此把AR模型和预测系统联系起来，用另一个方法估计AR模型参数。

从AR模型表达式可知其当前输出值与过去输出值有关。可以通过预测的方法利用信号的前后相关性估计未来的信号值，即前向预测，可以表示为：
$$\hat{x}(n)=-\sum _{k=1}^m{a}_m(k)x(n-k)$$
m为预测器的阶数，k=1,2,3…m。

预测结果与真实值之间的误差，即前向预测误差e(n):
$$e(n)=x(n)-\hat{x}(n)=x(n)+\sum _{k=1}^m{a}_m(k)x(n-k)$$
将e(n)看作系统的输出，x(n)看作输入，得到系统函数:
$$\frac{E(z)}{X(z)}=1+\sum _{k=1}^m{a}_m(k)z^{-k}$$
这里假设m=p,即两个系统——预测器和AR模型阶数相同，且预测系数和AR模型参数相同。

有w(n)=e(n)。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7OnSGXlF-1592722510260)(E:\大三下\数据压缩\作业七\两个系统.png)]

求预测误差均方值：
$$\begin{gathered} E[e^2(n)]=E[(x(n)+\sum _{k=1}^m{a}_m(k)x(n-k){)}^2] \\ =R_{xx}(0)+2[\sum _{k=1}^m{a}_m(k)R_{xx}(k)+\sum _{k=1}^m\sum _{l=1}^m{a}_m(l)a_m(k)R_{xx}(l-k)] \end{gathered}$$
要是均方误差最小，即对等式右侧求偏导使其等于零，得到m个等式：
$$R_{xx}(l)=-\sum _{k=1}^m{a}_m(k)+R_{xx}(l-k)l=1,2,...,m$$
带入均方误差表达式，可得：
$$\begin{gathered} R_{xx}(l)=R_{xx}(0)+\sum _{k=1}^m{a}_m(k)R_{xx}(k) \\ 或E_p[e^2(n)]=R_{xx}(0)+\sum _{k=1}^{m}ak{R}_{xx}(k) \\ {a}_k=a_m(k)m=p \end{gathered}$$
同时由于e(n)=w(n)，最小均方误差即等于白噪声方差。

如此，在计算AR模型参数的时候使用L-D算法，即使用Y-W方程式（递推式）。L-D递推算法是模型阶数逐渐增大的一种算法。

首先计算m=1时的预测系数和预测均方误差（白噪声），然后计算m=2的，一直计算到m=p。每一阶的参数都是从低一阶的参数推导得到的。这样在减小计算量的同时也可以同时寻找最佳阶数，满足精度要求即可停止递推，选定阶数。

大致流程图如下：

推导式为：
$$\{\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=[1-a_m^2(m)]E_{m-1}=R(0){\prod }_{k=1}^m[1-a_k^2(k)]\end{array}$$
其中am(m)称为反射系数。通过初始化E0=R(0)，a(0)=1，阶数p，就可以按照流程图进行AR模型参数估计。

L-D算法的优点时计算速度快，所得AR模型必定稳定，且均方误差随着阶数上升而下降。其缺点时在求自相关系数时，假设除了观测值之外的数据都为0，会引入较大误差。

实例

首先和之前一样，先根据观测值求自相关序列Rxx。

随后根据原理进行初始化，E0=Rxx(0)=1.9271,a0=1，阶数p.

然后根据递推式计算m=1，2，3，当p=3时所得均方误差和前一例题一致，说明这一阶数比较合理。

matlab中有专门实现L-D算法的函数aryule()。

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 = zeros(1, length(xn));
for m = 0 : length(xn) - 1
    i = 1;
    while i+m <= length(xn)
        Rxx(m+1) = Rxx(m+1) + xn(i) * xn(i + m);
        i = i + 1;
    end
    Rxx(m+1) = Rxx(m+1) / 32;
end
p = 4;
E = zeros(1, p); E(1) = Rxx(1);
a = zeros(p, 3); a(1, 1) = 1;
for m = 1:3
    a(m + 1, m) = Rxx(m + 1);
    for k = 1:(m - 1)
        a(m + 1, m) = a(m + 1, m) + a(m, k) * Rxx(m + 1 - k);
    end
    a(m + 1, m) = - a(m + 1, m) / E(m);
    for k = 1:(m - 1)
       a(m + 1, k) = a(m, k) + a(m + 1, m) * a(m, m - k);
    end 
    E(m + 1) = (1 - a(m + 1, m).^2) * E(m); 
end
%结果
a =

    1.0000         0         0 %m=0
   -0.8623         0         0 %m=1
   -0.6789   -0.2127         0 %m=2
   -0.6984   -0.2748    0.0915 %m=3=p


E =

    1.9271    0.4941    0.4718    0.4678 %m=0,1,2,3

也可以直接使用matlab提供的函数：

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];
[a,E]=aryule(xn,3)
%结果
a =

    1.0000   -0.6984   -0.2748    0.0915


E =

    0.4678

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

为了克服L-D算法的缺点，1968年Burg提出了Burg算法，基本思想是对数据同时机型前项和后向预测，两者的均方误差值和最小为估计准则，进而通过L-F算法递推得到AR模型参数。

其优点为计算效率高，所得AR模型稳定，但递推使用的仍然为L-D算法，没有解决根本缺点。

matlab中有专门的函数arburg()，使用方法如下：

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];
[a,E]=arburg(xn,3)
%结果
a =

    1.0000   -0.6982   -0.2626    0.0739


E =

    0.4567
%高阶计算
[a,E]=arburg(xn,12)
%结果
a =

  1 至 10 列

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

  11 至 13 列

    0.0382    0.1628   -0.2501


E =

    0.3237

计算速度比原来的方法快了很多。

1980年Marple在前人基础上提出了一种高效算法——Marple算法或称不受约束的最小二乘法（LS）。基本思想为：让每一个预测系数的确定直接与前向后向预测的总平方误差最小，此时系数不能从低一阶的系数推导得到了，脱离了L-D算法的大框架。实践表明这种算法更加优越，但不能保证AR模型的稳定性。

最终预测误差

FPE（final predidyion error）准则定义：给定观测长度N，从某个过程的一次观测数据中估计到了预测系数，然后用该预测系数构成的系统处理另一次观察数据，则有预测均方误差，该误差在某个阶数时为最小，其表达式为：
$$FPE(p)={\hat{\sigma }}_{wp}^2(\frac{N+P+1}{N-P-1})$$
随着阶数p的增加，方差减小，而括号内的值增加。因此可以找到最佳的p使FPE最小。