Levinson-Durbin 算法

Levinson-Durbin 算法
为了使读者对Levinson-Durbin 算法的应用背景有所了解，在此我们先简单地描述一下基于最小均方误差的线性预测编码（LPC ）算法。

假设)(n y 是一实数据列，+∞<<∞-n ，我们可以用过去时刻的N 个数据来预测当前时刻的数据)(n y ， 即：
+∞<<∞--=∑=n - , )()()(ˆ1
N
k N k n y k a n y
(A.1)
这里)(k a N 即为预测系数。

定义预测误差)(n e 为
)(ˆ)()(n y
n y n e -= (A.2)
我们将采用最小均方误差准则来选择)(k a N 的值，使得式(A.3)总误差N E 最小。

∑∑∑∞
-∞==∞
-∞=⎥⎦
⎤⎢⎣⎡-+==n N
k N n N k n y k a n y n e E 2
12
)()()()(
(A.3)
这种优化参数)(k a N 的方法导致了求解如下的正则方程组
N l l R l k R k a
yy N
k yy N
,,2, , )()()(1
=-=-∑=
(A.4)
这里的)(k R yy 是序列)(n y 的自相关系数。

式（A.4）可写成如下的矩阵形式：
N N N r a R -=⋅
（A.5）
其中
[]T
N N N N N a a a )(,),2(),1( =a
(A.6) []
T
yy yy yy N N R R R )(,),2(),1( =r
(A.7)
注意到)(k R yy 具有)()(j i R j i R yy yy -=-的性质，式(A.5)中的N R 可写成如下形式
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡----=)0()2()1()2()0()1()1()1()
0(yy yy yy yy yy yy yy yy yy N R N R N R N R R R N R R R R (A.8)
Levinson-Durbin 算法是求解正则方程组中的预测系数N a 的有效算法。

这种算法利用了自相关矩阵中特殊的对称性。

注意到)(),(j i R j i R yy N -=，即对角线上的元素都相等，所以这个自相关矩阵是Toeplitz 矩阵。

Levinson-Durbin 算法利用了Toeplitz 矩阵的特点来进行迭代计算。

即首先由一阶预测器(1=N )开始，计算预测系数)1(1a 。

然后增加阶数，利用低阶的结果得到下一个高阶的计算结果。

根据(A.4)式求解得到的一阶预测器的预测系数)1(1a 是：
)
0()1()1(1yy yy R R a -
= (A.9)
其最小均方误差是：
[]
212
111)1(1)0()0()1()1()1(2)0(a R R a R a R E yy yy yy yy -=+⋅+=
(A.10)
这里11)1(ρ=a 是格形滤波器的第一反射系数。

下一步是求解二阶预测器的系数)1(2a 和)2(2a ，并将结果用)1(1a 表示，根据式(A.5)得到的两个方程是:
)1()1()2()0()1(22yy yy yy R R a R a -=+ )2()0()2()1()1(22yy yy yy R R a R a -=+
（A.11）
通过用(A.9)的解来消去)1(yy R ，我们得到解：
[
]
1
121
12)
1()1()2()
1(1)0()1()1()2()2(E R a R a R R a R a yy yy yy yy yy ⋅+=
-⋅+-
=
)1()2()1()1(1212a a a a +=
（A.12）
这样我们得到了二阶预测器的预测系数，我们再次注意到22)2(ρ=a 是格形滤波器中的第
二反射系数。

据此类推，我们可以用)1(-m 阶预测器的预测系数来表示m 阶预测器的系数。

这样，我们可将m 阶预测系数矢量m a 写成两矢量的和，也就是
⎥⎦⎤
⎢⎣⎡+⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=--m m m m m m m m a a a ρ110)()2()1(d a
a
(A.13)
这里1-m a 矢量是第)1(-m 阶预测器的预测系数，)1(-m 维的矢量1-m d 和标量m ρ是待定的。

我们将m m ⨯自相关矩阵m R 分区如下：
⎥⎥⎦
⎤⎢⎢⎣⎡=---)0(1
1
1
yy bt
m b
m m m R r r R R
(A.14)
这里[]
T
b
m yy yy yy bt
m R m R m R )
()1(,),2(),1(11--=--=r r ，b
m 1-r 的上标
b 表示
[]
T
yy yy yy m m R R R )1(,),2(),1(1-=- r 元素的倒序排列。

根据式(A.13)和(A.14)，式(A.5)可以写成如下形式
⎥⎦⎤⎢⎣⎡-=⎪⎪⎭⎫ ⎝⎛⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣
⎡⎥⎥⎦⎤⎢⎢⎣⎡------)(0)0(11111
1m R R yy m m m m yy bt m b
m m r d a r r R ρ
(A.15)
这是Levinson-Durbin 算法中的关键一步，从式(A.15)中我们得到两个方程
111111-------=++m b
m m m m m m r r d R a R ρ
(A.16) )()0(1111m R R yy yy m m bt m m bt m -=++----ρd r a r
(A.17)
由于111----=m m m r a R ， 由式(A.16)得到
b m m m m 1111-----=r R d ρ
(A.18)
又由于b
m 1-r 仅是1-m r 倒序排列，且1-m R 是Toeplitz ，因此可得
b
m b m m r a R -=-1
(A.19)
即：
b m m b m r R a ⋅-=--11
(A.20)
因此式(A.18)可写成
⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎣⎡--==-----)1()2()1(11
11
1m m m m b m m m a m a m a ρρa d (A.21)
现在可用式(A.17)这个方程来求解m ρ，如果我们用式(A.21)来消去式(A.17)中的1-m d ，可得
[])()0(1111m R R yy m bt
m b m bt m yy m -=++----a r a r ρ （A.22）
注意到1111----⋅=⋅m t
m b m bt m a r a r ，由式(A.22)可得
1
111)0()(----++-
=m t m yy m bt
m yy m R m R a r a r ρ (A.23)
因此，通过用式(A.21)和(A.23)的结果，替换式(A.13)中的值，我们得到了求解预测器系数的Levinson-Durbin 迭代算法。

1
1
11
111)()0()()(-------+-
=++-
==m m bt
m yy m t m yy m bt
m yy m m E m R R m R m a a r a r a r ρ (A.24)
)
()()( )()()(1111k m a m a k a k m a k a k a m m m m m m m -+=-+=----ρ
p
m m k ,,2,11
,,2,1 =-=
(A.25)
最后我们来确定最小均方误差m E 的表达式。

对于m 阶预测器我们有
[][]
)()( )()()()()()()0( )
()()()()0( )
()()0(1111
111
111
111----=--=-=--=++=⎥⎦
⎤
⎢⎣⎡-+++=⋅-++=+=∑∑∑∑m bt
m yy m m m k yy m yy m m k yy m yy m
k yy m m m yy m
k yy m yy m m R m a E k R k m a m R m a k R k a R k R k m a m a k a R k R k a R E a r
(A.26)
利用公式（A.24）可得
[][]
2
1211)(1m m m m m E m a E E ρ-=-=--
(A.27)
这里)0(0yy R E =。

由于N m E m ,,1,0 , 0 =>，根据式（A.27）知反射系数满足1≤m ρ，这样预测器
最小均方误差序列满足条件：
N E E E E ≥≥≥≥ 210
（A.28）
至此，求解N N N r a R -=线性方程组的Levinson-Durbin 算法就介绍完了。