数据压缩学习笔记（四）语音参数编码与随机信号的AR模型参数建模

语音信号的特征

浊音（Voiced sounds）：声带振动，引起声门的打开和关闭，从而发送压力变化的脉冲到声道，在声道形成声波。
清音（Unvoiced sounds）：清音是声门保持打开并将气体压进一个收缩声道的结果。清音样本表现出很少的相关性，是随机的或接近随机的。
爆破音（Plosive sounds）：声道关闭之后产生的压缩空气然后突然打开声道所发出的音。

因此语音信号的生成可以被归为两大类：

在话音生成模型中，声道被等效成一个随时间变化的滤波器，它由白噪声—清音话音段激励，或者由脉冲串——浊音段激励。
因此语音编码需要编码的对象就是发声的种类和滤波器的参数。可见，对语音信号进行建模并获得参数的过程非常重要。

参数建模法

为随机信号建立参数模型是研究随机信号的一种基本方法，其含义是认为随机信号x(n)
是由白噪w(n)激励某一确定系统的响应。

参数模型可分为三种：

1. MA模型
   
   
   q表示系统阶数，系统函数只有零点，没有极点，所以该系统一定是稳定的系统，也称为全零点模型，用 MA（q）来表示。

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

3. ARMA模型
   
   
   它既有零点又有极点，所以也称极零点模型，要考虑极零点的分布位置，保证系统的稳定，用 ARMR（ p ，q）表示。

根据 Wold 的证明：任何平稳的 ARMA（自回归移动平均）模型或 MA 模型均可用无限阶或阶数足够的 AR 模型去近似。
在语音的参数编码中，声道模型一般使用AR模型。

AR 模型参数的估计

AR模型的模型为：

对上式两边同时乘以x(n-m)并求均值，有

因为自相关函数偶对称，所以上式化为：

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



将系统函数H(z)转换至时域，有：
因而h(0) = 1
带入(7-33)，可得Yule－Walker（Y-W）方程：


将上述二式下标简化，写成矩阵形式：

由于自相关函数是偶对称函数，因而自相关矩阵是对称矩阵，与主对角线平行的斜对角线的元素都是相同的。根据输入信号的自相关函数求解参数可用Levinson－Durbin（L-D）算法。

Levinson－Durbin（L-D）算法

由AR 模型的时域表达式：

可知当前输出值与它过去的输出值有关。

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

其预测误差为：

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

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


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

求得最小均方误差：

p 阶预测器的预测系数等于 p 阶 AR 模型的参数，由于e(n)＝w(n)，所以最小均方预测误差等于白噪声方差

L-D 递推算法是模型阶数逐渐加大的一种算法，先计算阶次 m＝1 时的预测系数，一直计算到m=p阶的预测系数。这种递推算法的特点是，每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。

按照上述推导中m个使偏导为0的等式和最小均方误差等式可得预测系数和均方误差估计的通式：

其中am(M) 为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值E0 = R(0), a0(0) = 1以及阶数p就可以对参数进行估计。

L-D 递推算法流程为：

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

对于上述算法，可用matlab中自带函数：[a E] = aryule（x，p）实现，其中x为输入数据，p为阶数。输出a为估计的模型参数，E为噪声信号的方差估计。

L-D算法的改进

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

Matlab 中有自带函数实现 Burg 算法的 AR 模型参数估计：[a E]＝arburg（x，p）

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

AR模型预测阶数的选择

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

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

AR模型预测在matlab中的实现
对如下模型生成的32点数据使用L-D算法和Brug算法分别进行3阶和12阶的估计：

数据：

data = [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];

matlab代码：

[a_LD3,E_LD3] = aryule(data,3)
[a_LD12,E_LD12] = aryule(data,12)
[a_Burg3,E_Burg3] = aryule(data,3)
[a_Burg12,E_Burg12] = aryule(data,12)

结果：

a_LD3 =
    1.0000   -0.6984   -0.2748    0.0915

E_LD3 =
    0.4678


a_LD12 =
  Columns 1 through 9
    1.0000   -0.6703   -0.3254   -0.0793    0.1407    0.3676   -0.2451    0.0483   -0.0912
  Columns 10 through 13
   -0.0522    0.0515    0.0186   -0.0955

E_LD12 =
    0.3783


a_Burg3 =
    1.0000   -0.6984   -0.2748    0.0915

E_Burg3 =
    0.4678


a_Burg12 =
  Columns 1 through 9
    1.0000   -0.6703   -0.3254   -0.0793    0.1407    0.3676   -0.2451    0.0483   -0.0912
  Columns 10 through 13
   -0.0522    0.0515    0.0186   -0.0955

E_Burg12 =
    0.3783

以三阶L-D算法输出为例，其结果a_LD3含义为：a1 = -0.06984；a2 = -0.2748；a3 = 0.0915。
对于两种算法而言，提高预测阶数均可降低均方误差。

使用最终预测误差（FPE：final predidyion error）准则确定阶数：

使用上述模型再生成两组数据：

dat1 = [-0.98976,0.76223,0.36297,2.0177,2.4194,1.4687,0.388,0.60336,0.10567,-0.57189,0.17865,1.3738,0.34575,-0.13732,-0.25403,0.44894,-0.68185,-1.4201,-1.2227,-1.2827,-0.66225,-1.4119,-1.1682,-1.618,-1.4024,0.15899,-0.97222,-1.7964,-0.94975,-2.0907,-1.465,-2.2266];
dat2 = [0.5354,0.8652,0.50178,-1.4595,-0.56668,0.69417,1.2716,-0.55468,0.045716,-0.91598,-1.5186,-2.4656,-1.6808,-2.2711,-1.7966,-2.1975,-2.3259,-1.735,-1.0645,0.93749,-1.1373,-1.2901,-0.1803,-0.93138,-1.9384,-1.157,0.19024,0.46577,2.3486,2.0413,3.9166,2.6125];

计算得到不同阶数时的预测均方误差（1-30阶）：

w = [1.1048 1.1077 1.1087 1.1096 1.0958 1.0983 1.1065 1.1018 1.1094 1.1218 1.1312 1.1319 1.1324 1.0847 1.0085 1.0125 1.0007 0.9875 0.9881 0.9861 1.0178 1.0342 1.0188 1.0289 1.0220 1.0193 1.0268 1.0327 1.0299 1.0088]

根据公式计算各阶FPE：

fpe = [1.2521 1.3368 1.4255 1.5206 1.6015 1.7134 1.8442 1.9640 2.1179 2.2971 2.4886 2.6807 2.8940 2.9990 3.0254 3.3075 3.5740 3.8739 4.2819 4.7513 5.4960 6.3200 7.1313 8.3784 9.8790 12.0279 15.4020 20.9979 31.9284 63.5575]

找到最小的FPE值，确定最佳预测阶数为1阶。

Matlab代码：

close all; clear; clc;

N = 32;
dat1 = [-0.98976,0.76223,0.36297,2.0177,2.4194,1.4687,0.388,0.60336,0.10567,-0.57189,0.17865,1.3738,0.34575,-0.13732,-0.25403,0.44894,-0.68185,-1.4201,-1.2227,-1.2827,-0.66225,-1.4119,-1.1682,-1.618,-1.4024,0.15899,-0.97222,-1.7964,-0.94975,-2.0907,-1.465,-2.2266];
dat2 = [0.5354,0.8652,0.50178,-1.4595,-0.56668,0.69417,1.2716,-0.55468,0.045716,-0.91598,-1.5186,-2.4656,-1.6808,-2.2711,-1.7966,-2.1975,-2.3259,-1.735,-1.0645,0.93749,-1.1373,-1.2901,-0.1803,-0.93138,-1.9384,-1.157,0.19024,0.46577,2.3486,2.0413,3.9166,2.6125];

w = zeros(1,N-2);
fpe = zeros(1,N-2);
for p=1:(N-2)
    [a,E] = aryule(dat1,p);
    w(p) = Predictor_error(dat2,a);
    fpe(p) = w(p)*((N+p+1)/(N-p-1));
end

w
fpe
p_opt = find(fpe==min(fpe))

其中Predictor_error()函数为：

function [ E ] = Predictor_error( data,a )

a_temp = a(2:length(a));

N_dat = length(data);
N_a = length(a_temp);
N_all = N_a + N_dat;
data_temp = [zeros(1,N_a) reshape(data,[1,N_dat])];
a_temp = flip(reshape(a_temp,[1,N_a]));

errors = zeros(1,N_dat);
for i=N_a+1:N_all
    p = data_temp(i-N_a:i-1)*(a_temp');
    errors(i-N_a) = data(i-N_a) + p;
end

E = mean(errors.^2);
end