随机信号的参数建模法

最新推荐文章于 2021-06-05 17:42:39 发布

MaaaaaCro

最新推荐文章于 2021-06-05 17:42:39 发布

阅读量461

点赞数

分类专栏：数据压缩文章标签： matlab

本文链接：https://blog.csdn.net/weixin_44949552/article/details/106878638

版权

数据压缩专栏收录该内容

17 篇文章 1 订阅

订阅专栏

随机信号的参数建模法

在对随机信号的研究和处理中，根据其参数建立相应的参数模型是一种重要的方法。随机信号的参数建模法可以应用在语音信号的编码过程中：对不同的物理过程产生的声音信号建立数学模型，根据语声模型提取语音信号的特征参数，对参数进行编码传输，接收端根据特征参数重建语音信号。

随机信号的参数建模法将随机信号 $x (n)$ 看做是由白噪声 $w (n)$ 激励某一确定系统的响应：

随机信号的参数建模法
对于平稳随机信号，有三种常用的线性模型： AR 模型（自回归模型 Auto-regression model），MA 模型（滑动平均模型 Moving average model）和 ARMA 模型（自回归滑移平均模型 Auto-regression-Moving average model）。

AR模型

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

$x(n)=w(n)-\sum_{k=0}^pa_kx(n-k)$

其中p表示系统阶数。进行z变换后，可以得到模型的系统函数：

$H(Z)=\frac{1}{1+\sum_{k=0}^pa_kz^{k}}$
AR模型系统函数的特点是只有极点，没有零点，为全极点模型，要建立稳定的系统，则要注意极点的分布，用 $A R (p)$ 来表示。

MA模型

随机信号 $x (n)$ 由当前的激励 $w (n)$ 和若干次过去的激励 $w (n - k)$ 线性组合产生：
$x(n)=\sum_{k=0}^qb_kw(n-k)$

其中q表示系统阶数。进行z变换后，可以得到模型的系统函数：
$H(Z)=\frac{X(Z)}{W(Z)}=\sum_{k=0}^qb_kz^{-k}$

MA模型系统函数的特点是只有零点，没有极点，为全零点模型，故系统稳定。

ARMA模型

ARMA模型综合了AR模型与MA模型：

$x(n)=\sum_{k=0}^qb_kw(n-k)-\sum_{k=1}^pa_kx(n-k)$

模型的系统函数为；

$H(Z)=\frac{\sum_{k=0}^qb_kz^{-k}}{1+\sum_{k=1}^pa_kz^{k}}$

系统既有零点，又有极点，考虑到系统的稳定性，要注意零极点的分布，用 $A R M A (p ， q)$ 表示。

AR模型的参数估计

实际应用中，三种模型选用哪一种，就要考虑到节约和计算量。选定模型后，剩下的任务就是用适当的算法估计模型参数 $a_k，b_k，p，q)$ ，以便用模型对随机信号进行预测。在信号处理中，AR模型的应用较为广泛，我们以AR模型为例进行参数估计：

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

根据AR模型的公式 $x(n)=w(n)-\sum_{k=0}^pa_kx(n-k)$ ，两边同乘以 $x (n - m)$ ,再求均值得到：
$E[x(n)x(n-m)]=E[w(n)x(n-m)-\sum_{k=0}^pa_kx(n-k)x(n-m)]$
根据自相关函数的对称性：
$R_{xx}(m)=R_{xx}(-m)$
将上述式子化为：
$R_{xx}(m)=R_{xw}(m)-\sum_{k=0}^pa_kR_{xx}(m-k)$
由于系统的单位冲激响应 $h (n)$ 是因果的，所以输出的平稳随机信号和输入的白噪声之间的互相关函数有下列推导：

$x(n)=w(n)*h(n)=\sum_0^∞h(k)w(n-k)$
$R_{xw}(m)=E[\sum_0^∞h(k)w(n-k)w(n+m)]=\sum_0^∞h(k)E[w(n-k)w(n+m)]\\=\sum_0^∞h(k)R_{ww}(n+m)=\sum_0^∞h(k)σ^2_wδ_w(m+k)\\=σ^2_wh(-m)$
根据因果性可知：
$R_{xw}(m)= \bigg \{ \begin{aligned} σ^2_wh(-m) & \ & \ m>0\\ 0 & \ & \ m\leq0 \end{aligned}$

代入到 $R_{xx}(m)=R_{xw}(m)-\sum_{k=0}^pa_kR_{xx}(m-k)$ 式子中，便得到：
$R_{xx}(m)= \left\{ \begin{aligned}-\sum_{k=1}^pa_kR_{xx}(m-k) & \ & \ m>0\\\ -\sum_{k=1}^pa_kR_{xx}(m-k) +σ^2_wh(0) & \ & \ m=0\\R_{xx}(-m) & \ & \ m\leq0\end{aligned} \right.$

将 $H (Z)$ 进行逆Z变换，转换到时域有： $h(n)+\sum_{k=1}^pa_kh(n-k)=δ(n)$ ,故 $h (0) = 1$
显然，AR模型具有递推特性：
$R_{xx}(m)=-\sum_{k=1}^pa_kR_{xx}(m-k) \quad \quad m>0$
这就是著名的Yule－Walker方程，将上式变换：
$R_{xx}(m)+\sum_{k=1}^pa_kR_{xx}(m-k)=0 \quad \quad m>0$
输入的白噪声方差为：
$σ^2_w=R_{xx}(0)+ -\sum_{k=1}^pa_kR_{xx}(-k)$

将Y-W方程和 $σ^2_w$ 结合，可以写成矩阵形式，得到矩阵方程：
$\begin{bmatrix} R(0) & R(-1)&...&\ R(-p)\\ R(1)&R(0) & ... &R(-p+1)\\ ...&...& &...\\ R(p) & R(p-1)&...&\ R(0)\\ \end{bmatrix} \begin{bmatrix}1\\a_1\\...\\a_p \end{bmatrix} = \begin{bmatrix}σ^2_w\\0\\...\\0\end{bmatrix}$

上述矩阵方程组的系数均为自相关函数，由于自相关函数具有对称性，故其组成的系数矩阵为对称矩阵，存在着更高效的算法，其中常用的有Levinson－Durbin（L-D）算法。Yule－Walker（Y-W）方程表明，只要已知输出平稳随机信号的自相关函数，就能求出 AR 模型中的参数 $a_k$ ，并且所需要的数据较少。

实例分析：

已知自回归信号模型 AR（3）为： $\frac{14}{24}x(n-1)+ \frac{9}{24}x(n-2)- \frac{1}{24}x(n- 3)+w(n)$ 式中 $w (n)$ 是具有方差 $σ^2_w=1$ 的平稳白噪声，求:
1 .自相关序列 $R_{xx}(m)$ ，m＝0，1，2，3，4，5。
2 . 用 1.中求出的自相关序列来估计 AR（3）的参数 $a'_k$ ，以及输入白噪声的方差 $σ'^2_w$ 大小。
3 . 利用给出的 AR 模型，用计算机仿真给出 32 点观测值 $x(n)=[\quad0.4282\quad 1.1454 \quad 1.5597 \quad 1.8994 \quad 1.6854 \quad 2.3075 \quad 2.4679 \quad1.9790 \quad 1.6063 \quad 1.2804 \quad -0.2083 \quad 0.0577 \quad 0.0206 \quad 0.3572 \quad 1.6572 \quad 0.7488 \quad 1.6666 \quad1.9830 \quad 2.6914 \quad 1.2521 \quad 1.8691 \quad 1.6855 \quad0.6242 \quad0.1763 \quad 1.3490 \quad 0.6955 \quad 1.2941 \quad 1.0475\quad 0.4319 \quad 0.0312 \quad 0.5802\quad -0.6177\quad]$ ,用观测值的自相关序列直接来估计 AR(3)的参数 $a'_k$ 以及输入白噪声的 $σ'^2_w$ 。

解：
1、已知模型参数 $a_k$ ， $a_1=-\frac{14}{24}$ , $a_2=-\frac{9}{24}$ , $a_3=\frac{1}{24}$ ,代入矩阵中得到：
$\begin{bmatrix} 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{bmatrix} \begin{bmatrix}1\\-\frac{14}{24}\\-\frac{9}{24}\\\frac{1}{24} \end{bmatrix} = \begin{bmatrix}1\\0\\0\\0\end{bmatrix}$
调用MATLAB解方程：

a=[-14/24 -9/24 1/24];
R=[ 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];
w=[1;0;0;0];
Rxx=R\w

解得结果： $R (0) = 4.9377$ , $R (1) = 4.3287$ , $R (2) = 4.1964$ ， $R (3) = 3.8654$

根据递推性质，可以的得出：
$R_{xx}(4)=-\sum_{k=1}^3a_kR_{xx}(4-k)=3.6481 \\ R_{xx}(5)=-\sum_{k=1}^3a_kR_{xx}(5-k)=3.4027$

2、已知自相关序列值，代入方程：
$\begin{bmatrix} 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{bmatrix} \begin{bmatrix}1\\a_1\\a_2\\a_3 \end{bmatrix}= \begin{bmatrix}σ^2_w\\0\\0\\0\end{bmatrix}$

解得：
$a_1=-\frac{14}{24}$ , $a_2=-\frac{9}{24}$ , $a_3=\frac{1}{24}$ , $σ^2_w=1$

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

3、利用给出的 32 点观测值，先求自相关序列，根据自相关的定义：
$R_{xx}=\frac{1}{n}\sum_{i=1}^nx_ix_{i+m}$
用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];

Rxx=xcorr(xn) ./ length(xn)

由于自相关函数的对称性，在得到的结果的中心处为 $R_{xx}(0)$ 的值

.........
//第31列为m=0时Rxx的值
23 至 33 列

    0.8058    0.7637    0.7520    0.8673    0.9060    1.1349    1.3545    1.5381    1.6618    1.9271    1.6618

34 至 44 列

    1.5381    1.3545    1.1349    0.9060    0.8673    0.7520    0.7637    0.8058    0.8497    0.8761    0.9608
...........

这里我们只用到 $R_{xx}(m)$ 当 m=0,1,2,3 的值
计算出 $R_{xx}(m)=[1.9271,1.6618 ,1.5381 ,1.3545 ]\quad\quad m=0,1,2,3$ 代入方程中，得到估计值： $a'_1=-0.6984$ , $a'_2=-0.2748$ , $a'_3=0.0915$ , $σ'^2_w=0.4678$
估计值与真实AR模型参数的误差为： $e_1=0.1151$ , $e_2=0.1002$ , $e_3=0.0498$ , $e_σ=0.5322$ 。造成误差的原因在于只观测了32点数据，数据量不足，计算机仿真的32点白噪声方差也不等于1。给出一段观测值求 AR 模型参数这样直接解方程组，当阶数越高时直接解方程组计算就越复杂，因而要用特殊的算法使得计算量减小且精确度高。

2、 Y－W 方程的解法——L-D 算法

使用直接观测的数据估计AR模型的参数，所需要的阶数需要很大，才能获得精度较高的估计值，阶数增大必然造成计算量的增加，因此把 AR 模型和预测系统联系起来，换个方法来估计 AR 参数。

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

若序列的模型已知而用过去观测的数据来推求现在和将来的数据称为前向预测器，表示为： $\hat{x}(n)=-\sum_{k=1}^ma_m(k)x(n-k)$
其中 $a_k$ 代表 m 阶预测器的预测系数。设预测误差为 $e (n)$ :
$e(n)=x(n)-\hat{x}(n)=x(n)+\sum_{k=1}^ma_m(k)x(n-k)$
把 $e (n)$ 看成是系统的输出， $x (n)$ 看成是系统的输入，得到系统函数：
$\frac{E(Z)}{X(Z)}=1+\sum_{k=1}^ma_m(k)z^{-k}$
假如 m＝p，且预测系数和 AR 模型参数相同,把预测误差系统框
图和 AR 模型框图给出有 $w (n) = e (n)$ ,即前向预测误差系统中的输入
为 $x (n)$ ，输出为预测误差 $e (n)$ 等于白噪声。也就是说前向预测误差系统对观测信号起了白化的作用。由于 AR 模型和前向预测误差系统有着密切的关系，两者的系统函数互为倒数，所以求 AR 模型参数就可以通过求预测误差系统的预测系数来实现。
在这里插入图片描述
求预测误差的均方值：
$E[e^2(n)]=e[(x(n)+\sum_{k=1}^ma_m(k)x(n-k))^2]\\=R_{xx}(0)+2[\sum_{k=1}^ma_m(k)R_{xx}(k)]+\sum_{k=1}^m\sum_{l=1}^ma_m(l)a_m(k)R_{xx}(l-k)$

要使得均方误差最小，将上式右边对预测系数求偏导并且等于零，得到 m 个等式：
$R_{xx}(l)=-\sum_{k=1}^ma_m(k)R_{xx}(l-k)$
求得最小均方误差为：
$E_m[e^2(n)]=R_{xx}(0)+\sum_{k=1}^ma_m(k)R_{xx}(k)$
$E_p[e^2(n)]=R_{xx}(0)+\sum_{k=1}^pa(k)R_{xx}(k)$
所以对照上述二式得知：
$a_k=a_m(k) \quad\quad m=p$ 也就是 p 阶预测器的预测系数等于 p 阶 AR 模型的参数,由于 $e (n) = w (n)$ ，所以二最小均方预测误差等于白噪声方差 $σ^2_w$

L-D 算法
L-D算法的基本思想就是根据Y-W 方程，自相关序列具有递推的性质，，L-D 递推算法是模型阶数逐渐加大的一种算法，先计算阶次 m＝1 时的预测系数 $a_1(1)$ 和 $σ^2_{w1}$ ,再计算m＝2 时的 $a_2(1)$ 、 $a_2(2)$ 和 $σ^2_{w2}$ ，一直计算到 m＝p 阶时的预测系数。
每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作
量又便于寻找最佳的阶数值，满足精度时就停止递推。

下面进行对 $R_{xx}(l)=-\sum_{k=1}^ma_m(k)R_{xx}(l-k)$ 和 $E_p[e^2(n)]=R_{xx}(0)+\sum_{k=1}^pa(k)R_{xx}(k)$ 式子的递推 ,取m=1：
$\left\{ \begin{aligned}R(1)=-a_1(1)R(0) \quad (1)\\E_1=R(0)+a_1(1)R(1) \quad (2)\end{aligned} \right.$
$σ^2_{w1}=E_1$
m=2时：
$\left\{ \begin{aligned}R(1)=-a_2(1)R(0)-a_2(2)R(1) \\R(2)=-a_2(1)R(1)-a_2(2)R(0) \end{aligned} \right.$
估计的方差：
$E_2=R(0)+a_2(1)R(1)+a_2(2)R(2)$
将m=1时的 $R (1)$ , $E (1)$ 代入，可以解得m=2时预测系数。一直递推下去，可以得到预测系数和预测均方误差估计的通式：
$\left\{ \begin{aligned}a_m(k)= a_{m-1}(k)+a_m(m)a_{m-1}(m-k)\quad (1)\\a_m(m)=-\frac{R(m)+\sum_{k=1}^ma_{m-1}(k)R(m-k)}{E_{m-1}} \quad (2)\\E_m=σ^2_{wm}=[1-a_m^2(m)]E_{m-1}=R(0)\prod_{k=1}^m[1-a_k^2(k)]\quad (3)\end{aligned} \right.$

其中 $a_m(m)$ 称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值 $E_0=R(0)$ 、 $a_0(0)=1$ ，以及AR模型的阶数p，就可以依照以下的流程图进行估计：

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

实例分析：

已知自回归信号模型 AR（3）为： $\frac{14}{24}x(n-1)+ \frac{9}{24}x(n-2)- \frac{1}{24}x(n- 3)+w(n)$ 式中 $w (n)$ 是具有方差 $σ^2_w=1$ 的平稳白噪声，利用给出的 AR 模型，用计算机仿真给出 32 点观测值 $x(n)=[\quad0.4282\quad 1.1454 \quad 1.5597 \quad 1.8994 \quad 1.6854 \quad 2.3075 \quad 2.4679 \quad1.9790 \quad 1.6063 \quad 1.2804 \quad -0.2083 \quad 0.0577 \quad 0.0206 \quad 0.3572 \quad 1.6572 \quad 0.7488 \quad 1.6666 \quad1.9830 \quad 2.6914 \quad 1.2521 \quad 1.8691 \quad 1.6855 \quad0.6242 \quad0.1763 \quad 1.3490 \quad 0.6955 \quad 1.2941 \quad 1.0475\quad 0.4319 \quad 0.0312 \quad 0.5802\quad -0.6177\quad]$ ,用 L-D 算法来估计 AR(3)的参数 $\hat{a}_k$ 以及输入白噪声的 $\hat{σ}^2_w$

利用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];
p = 4;
Rxx = zeros( 1, length(xn) );

for i = 0 : length(xn) - 1 
    j = 1; 
    while j + i <= length(xn)
        Rxx(i+1) = Rxx(i+1) + xn(j) * xn(j + i );
        j = j + 1;
    end
    Rxx(i+1) =  Rxx(i+1) / length(xn);
end


E = zeros(1, p); E(1) = Rxx(1);
a = zeros(p, 3); a(1 , 1) = 1;
ak=zeros(1,3);

for i = 1:3
    a( i + 1, i ) = Rxx(i + 1);
    for k = 1:( i - 1 )
        a( i + 1, i ) = a( i + 1, i ) + a( i, k ) * Rxx( i + 1 - k );
    end
    a( i + 1, i ) = - a( i + 1, i ) / E(i);
    
    for k = 1:(i - 1)
       a( i + 1, k ) = a( i, k ) + a( i + 1, i ) * a( i, i - k );
    end 
        
    E( i + 1 ) = ( 1 - a( i + 1, i ).^2 ) * E(i); 
end

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

求得结果： $a'_1=-0.6984$ , $a'_2=-0.2748$ , $a'_3=0.0915$ , $σ'^2_w=0.4678$ ，与前文中用解方程所得结果一致。

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

L-D 算法也存在缺点：由于在求自相关序列时，是假设除了观测值之外的数据都为零，必然会引入较大误差。为了克服L-D 算法的缺点，可以选择其他优化的算法。

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

Marple 算法
Marple算法的思想是，让每一个预测系数的确定直接与前向、后向预测的总的平方误差最小，这样预测系数就不能由低一阶的系数递推确定了，所以不能用 L-D 算法求解，实践表明该算法比 L-D、Burg 算法优越。由于该算法是从整体上选择所有的模型参数达到总的均方误差最小，与自适应算法类似，不足是该算法不能保证 AR 模型的稳定性。

阶数选择
AR 模型的阶数选择不同得到的模型不同，效果相差较大，因而如何选择阶数很重要。
其中基于均方误差最小的最终预测误差（FPE：final predidyion error）准则是确定 AR 模型阶次比较有效的准则。

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

$FPE(p)=\hat{σ}^2_{wm}(\frac{N+P+1}{N-P-1})$
上式中估计的方差随着阶数的增加而减小，而括号内的值随着的增加而增加，因而能找到一最佳的阶数p，使得 FPE 最小。

MaaaaaCro

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
随机信号的参数建模法

随机信号的参数建模法在对随机信号的研究和处理中，根据其参数建立相应的参数模型是一种重要的方法。随机信号的参数建模法可以应用在语音信号的编码过程中：对不同的物理过程产生的声音信号建立数学模型，根据语声模型提取语音信号的特征参数，对参数进行编码传输，接收端根据特征参数重建语音信号。随机信号的参数建模法将随机信号 x(n)x(n)x(n) 看做是由白噪声 w(n)w(n)w(n) 激励某一确定系统的响应：对于平稳随机信号，有三种常用的线性模型： AR 模型（自回归模型 Auto-regression m
复制链接

扫一扫