线性预测之Levinson-Durbin算法

95 阅读 0 评论 63 点赞

我是靠谱客的博主缓慢皮皮虾，最近开发中收集的这篇文章主要介绍线性预测之Levinson-Durbin算法，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

自相关法中求解预测系数方程组为

$sum_{k=1}^{p}a_{k}R[|i-k|]=R[i] , 1leqslant ileqslant p$

方程的矩阵形式表示为：

$Ra=r$

式中R是正定对称的拓普利兹矩阵，且第 $left ( i,j right )$ 个元素为 $R[i-j]$ ，a和r分别是元素为 $a_{i}$ 和 $r[i]=R[i]$ 的列向量。

上式中方程组可做如下变换：

$R[i]-sum_{k=1}^{p}a_{k}R[|i-k|]=0,i=1,2,cdots ,p$

写成矩阵形式为

$begin{bmatrix} R[1]& R[0] & R[1]& cdots &R[p-1] \ R[2]& R[1]& R[0]& cdots &R[P-2] \ vdots & vdots & vdots & cdots &vdots \ R[p]& R[p-1]&R[P-2] & cdots & R[0] end{bmatrix},begin{bmatrix} 1\ -alpha _{1}^{(p)}\ -alpha _{2}^{(p)}\ cdots \ -alpha _{p}^{(p)} end{bmatrix}=begin{bmatrix} 0\ 0\ vdots \ 0\ end{bmatrix}$

根据最小均方误差公式，对于一个p阶最佳预测器，其最小均方预测误差为

$R[0]-sum_{k=1}^{p}a_{k}R[k]=varepsilon ^{(p)}$

写成矩阵的形式为

$[R[0], R[1], R[2]cdots R[p]],begin{bmatrix} 1\ -alpha _{1}^{(p)}\ -alpha _{2}^{(p)}\ cdots \ -alpha _{p}^{(p)} end{bmatrix}=[varepsilon ^{(p)}]$

把上面两个矩阵合在一起，写成新p+1个方程，方程组满足p各未知预测器系数和对应的未知均方预测误差 $varepsilon ^{(p)}$ ,新的矩阵方程为

$begin{bmatrix} R[0]& R[1] & R[2]& cdots &R[p] \ R[1]& R[0] & R[1]& cdots &R[p-1] \ R[2]& R[1]& R[0]& cdots &R[P-2] \ vdots & vdots & vdots & cdots &vdots \ R[p]& R[p-1]&R[P-2] & cdots & R[0] end{bmatrix},begin{bmatrix} 1\ -alpha _{1}^{(p)}\ -alpha _{2}^{(p)}\ cdots \ -alpha _{p}^{(p)} end{bmatrix}=begin{bmatrix} varepsilon ^{(p)}\ 0\ 0\ vdots \ 0\ end{bmatrix}$

上式中构造出的 $(p+1)times (p+1)$ 阶矩阵也是拓普利兹矩阵。该方程组用Levinson-Durbin算法递归求解。该算法通过在每次迭代中顺序地结合一个新的相关值来实现，并且根据新的相关值和已获得的预测器就能解出下一个高一阶的预测器。

对于任意阶数i，上式中的方程组都可以以矩阵形式表示：

$R^{i}alpha ^{i}=e^{i}$

我们希望第i阶的解是如何由第i-1阶的解导出的。换而言之是给定 $a^{(i-1)}$ ，即 $R^{(i-1)}alpha ^{(i-1)}=e^{(i-1)}$ 的解，我们希望导出 $R^{(i)}alpha ^{(i)}=e^{(i)}$ 的解。

1、先将矩阵方程 $R^{(i-1)}alpha ^{(i-1)}=e^{(i-1)}$ 以扩展形式写为

$begin{bmatrix} R[0]& R[1] & R[2]& cdots &R[i-1] \ R[1]& R[0] & R[1]& cdots &R[i-2] \ R[2]& R[1]& R[0]& cdots &R[i-3] \ vdots & vdots & vdots & cdots &vdots \ R[i-1]& R[i-2]&R[i-3] & cdots & R[0] end{bmatrix},begin{bmatrix} 1\ -alpha _{1}^{(i-1)}\ -alpha _{2}^{(i-1)}\ cdots \ -alpha _{i-1}^{(i-1)} end{bmatrix}=begin{bmatrix} varepsilon ^{(i-1)}\ 0\ 0\ vdots \ 0\ end{bmatrix}$

2、将数0附加到向量 $a^{(i-1)}$ 中，并与矩阵 $R^{(i)}$ 相乘，这时得到新的一组i+1个方程：

$begin{bmatrix} R[0]& R[1] & R[2]& cdots &R[i] \ R[1]& R[0] & R[1]& cdots &R[i-1] \ R[2]& R[1]& R[0]& cdots &R[i-2] \ vdots & vdots & vdots & cdots &vdots \ R[i-1]& R[i-2]&R[i-3] & cdots & R[1] \ R[i]& R[i-1]&R[i-2] & cdots & R[0] end{bmatrix},begin{bmatrix} 1\ -alpha _{1}^{(i-1)}\ -alpha _{2}^{(i-1)}\ cdots \ -alpha _{i-1}^{(i-1)}\ 0\ end{bmatrix}=begin{bmatrix} varepsilon ^{(i-1)}\ 0\ 0\ vdots \ 0\ gamma ^{(i-1)}\ end{bmatrix}$

上式成立的条件是：

$gamma ^{(i-1)}=R[i]-sum_{j=1}^{i-1}alpha _{j}^{(i-1)}R[i-j]$

3、根据拓普利兹矩阵 $R^{(i)}$ 的对称性，方程组可以以相反顺序写出(第一个方程写在最后面，最后一个方程写在最前面，以此类推)，可得

$begin{bmatrix} R[0]& R[1] & R[2]& cdots &R[i] \ R[1]& R[0] & R[1]& cdots &R[i-1] \ R[2]& R[1]& R[0]& cdots &R[i-2] \ vdots & vdots & vdots & cdots &vdots \ R[i-1]& R[i-2]&R[i-3] & cdots & R[1] \ R[i]& R[i-1]&R[i-2] & cdots & R[0] end{bmatrix},begin{bmatrix} 0\ -alpha _{i-1}^{(i-1)}\ -alpha _{i-2}^{(i-1)}\ cdots \ -alpha _{1}^{(i-1)}\ 1\ end{bmatrix}=begin{bmatrix} gamma ^{(i-1)}\ 0\ 0\ vdots \ 0\ varepsilon ^{(i-1)}\ end{bmatrix}$

4、将2和3两步中的公式合并可得

$R^{i}cdot begin{bmatrix} begin{bmatrix} 1\ -alpha _{1}^{(i-1)}\ -alpha_{2}^{(i-1)}\ vdots \ -alpha _{i-1}^{(i-1)}\ 0 end{bmatrix}& - &k _{i} begin{bmatrix} 0\ -alpha _{i-1}^{(i-1)}\ -alpha_{i-2}^{(i-1)} \ vdots \ -alpha _{1}^{(i-1)}\ 1 end{bmatrix} end{bmatrix}= begin{bmatrix} begin{bmatrix} varepsilon ^{(i-1)}\ 0\ 0\ vdots \ 0\ gamma ^{(i-1)} end{bmatrix}& - &k _{i} begin{bmatrix} gamma ^{(i-1)}\ 0\ 0\ vdots \ 0\ varepsilon ^{(i-1)} end{bmatrix} end{bmatrix}$

该式逼近了表达式 $R^{(i)}alpha ^{(i)}=e^{(i)}$ ，然后选择 $gamma ^{(i-1)}$ 使方程右边的向量只有一个非零元素，此时新参数 $k_{i}$ 必须为

$k_{i}=frac{gamma ^{(i-1)}}{varepsilon ^{i-1}}=frac{R[i]-sum_{j=1}^{i-1}alpha _{j}^{(i-1)}R[i-j]} {varepsilon ^{i-1}}$

这样就确保方程右边的向量的最后一个元素为0，则此时第一个元素为

$varepsilon ^{i}=varepsilon ^{i-1}-k_{i}gamma ^{(i-1)}=varepsilon ^{i-1}(1-k_{i}^{2})$

5、选定 $gamma ^{(i-1)}$ 后，第i阶的预测系数向量为

$begin{bmatrix} 1\ -alpha _{1}^{(i)} \ -alpha _{2}^{(i)}\ vdots \ -alpha _{i-1}^{(i)}\ -alpha _{i}^{(i)} end{bmatrix}=begin{bmatrix} 1\ -alpha _{1}^{(i-1)}\ -alpha_{2}^{(i-1)}\ vdots \ -alpha _{i-1}^{(i-1)}\ 0 end{bmatrix}& - &k _{i} begin{bmatrix} 0\ -alpha _{i-1}^{(i-1)}\ -alpha_{i-2}^{(i-1)} \ vdots \ -alpha _{1}^{(i-1)}\ 1 end{bmatrix}$

6、可得系数更新方程组为

$begin{matrix} alpha _{j}^{(i)}=alpha _{j}^{(i-1)} -k_{i}alpha _{i-j}^{(i-1)} ,j=1,2,cdots ,i-1 &\alpha _{i}^{(i)}=k_{i} & end{matrix}$

因此，对于某个特定的阶数p，最佳预测系数为 $alpha _{j}=alpha _{i}^{(i)}=k_{i},j=1,2,cdots ,p$

以上就是Levinson-Durbin算法的核心步骤，其算法的伪代码如下