文件类型：未知 文件大小：字节
1
第第第第五五五五章章章章 预测预测预测预测编码编码编码编码
对信号样值进行预测的概念在第2章讨论差值量化时已经提出.在那里,预测器只是一
个能输出预测信号的功能块,并指出对差值信号进行量化编码的比特率要低于对原信号的量
化编码的比特率.
本章将讨论预测器的工作原理,其中主要讨论线性预测的原理,算法以及它在语音,图
像信号编码中的应用.
由信号的取样值构成的信号序列中样值间存在相关性.利用相关性一个成功可追溯到
Wiener在1949年发表的论文"平稳时间序列的外推,内插,平滑及其在工程中的应用".由
于Wiener的开拓性工作,使滤波理论发展到一个新的阶段,即用统计学的观点来描述通信与
控制问题.整个滤波理论应包括对时间序列的滤波,预测与平滑三个方面.这些工作都是以
信号样值之间的相关性为依据的.如果在n时刻从噪声污染中复现出信号( )x n的近似值为
^( )x n,这便是滤波问题;如果在n时刻输出^( )x n k+,即对n + k时刻信号值的估计,则这是
预测问题;如果在n时刻输出^( )x n k-,即对n - k时刻信号的修正,则称作平滑问题.这里,
我们关心只是预测问题,但必须搞清楚时间序列分析的一些基本问题,诸如时间序列的概率
模型,时域估计等.
应用线性预测来进行数据压缩编码的一个种实现就是图像与语音信号的DPCM编码方
法.这个方法的侧重点在于寻找一个预测器,该预测器能够减小差值信号的方差,从而也能
够减小量化误差或降低数码率.在语音的线性预测编码中则另辟蹊径,用准周期脉冲与随机
噪声去激励一个线性时变系统所产生的输出作为语音主号.不同的语音所对应的线性时变系
统的参数也不同.应用线性预测的方法对这些时变的参数进行估值与编码来代替对语音波形
的直接编码,可以达到数据压缩的目的.这种能够用准周期脉冲与随机噪声合成语音信号的
时变参数系统,叫作线性预测声码器LPC(Linear Prediction Vocoder).目前用2400 b/s的速率
传输LPC参数已可获得较满意的话音质量,比较PCM编码所需64000 b/s的数据率,其压缩
比约为27.因此,这是一个很有效的数据压缩方法,在5.4节后面我们会更详细的讨论.5.5
节讨论图像线性预测编码的方法,主要讨论帧内与帧间的DPCM概念.
在讨论线性预测方法之前,先学习与介始随机过程的线性模型与线性预测的基本原理.
这是很重要的理论基础.上面已提到,它不仅应用于线性预测,而且也应用于整滤波理论与
2
现代谱分析研究.
5.1 时间序列的概率模型时间序列的概率模型时间序列的概率模型时间序列的概率模型
本节讲述时间序列的几种概率模型,它们总称为随机模型.根据信息率的观点,载有信
息的信号总是含有随机(也叫不确定)因素的.事实上,如果尚未听到声音就知道对方想说什
么,则从信息论意义上廛,这个语音信号没有意义,可以完全取消.但从其他方面考虑,比
如欣赏音乐,还是有意义的.
从数学上说,随机过程可定义为随机变量的一个集合{ ( ), }X t t T∈,其中T为定义随机过
程的时间点的集合.如果T是连续的(即,t-∞ <
= =
-
- <
∑
∑
(5.1.15)
由式(5.1.15)可见,自相关系数ρ是k的偶函数,且在k q=处被"截尾".这是MA过程的一
个特点.对于01β=的MA(1)过程,自相关系统为:
1
2
1
1 0
( ) 1
1
0
k
k k
β
ρ
β
=
= = ±
+
其他
(5.1.16)
3. MA过程的可逆性过程的可逆性过程的可逆性过程的可逆性
7
为了说明可逆性问题,先考察下面两个1阶的MA过程.
过程A:
1t t tX Z Zβ-= +
过程B:
1
1
t t tX Z Z
β-= +
由式(5.1.16)很容易看出,这两个不同的过程具有完全一样的自相关系数.这就说明,我们不
能从给定的自相关系数惟一地识另一个MA过程.
那么,{ }iβ满足什么条件下能使过程可逆呢 对此,让我们作进一步考察.对过程A和
B做如下改写:
过程A:
2
1 1 2 1 2( )t t t t t t t t tX Z Z Z Z Z Z Z Zβ β β β β- - - - -= - = - - = = - + -
过程B:
1 12
1 1
t t t tX Z Z Z
β β- -= - + -
如果| | 1β<,则A的级数收敛而B的级数不收敛,故称过程A是可逆的,而过程B不可逆的.
也就是说,条件| | 1β<使自相关系统能惟一对应一个MA过程.
现在来讨论一般情况下的可逆性条件.为了讨论方便,引入"后移算子"B,其定义为:
j
t t jB X X-= (5.1.17)
其中j为整数.于是,式(5.1.11)可写成,
0 1( ) ( )q
t q t tX B B Z B Zβ β β θ= + + + = (5.1.18)
仿照上面的改写方法,式(5.1.18)又可写成:
1( )t tZ B Xθ-= (5.1.19)
其中( )Bθ为B的q次多项式.若将B看作复变量而不再作算子,且设多项式,
0 1( )q
qB B Bθ β β β= + + + (5.1.20)
的根为1 2, , ,qb b b ,则多项式( )Bθ可以分解因式,并有:
8
1
1
1
1
( )
( )
q
i
q
ii
i
i
k
B
b B
b B
θ-
=
=
= =
-
-
∑
∏
(5.1.21)
其中ik是由1 2, , ,qb b b 确定的常数.
将式(5.1.21)代入式(5.1.19),并将1/( )ib B-展成,
2
2
1
1
i i i
B B
b b b
+ + +
得到
2
1
1
q
i
tt
ii i i
kB B
Z X
b b b=
= + + +
∑ (5.1.22)
现在应该恢复算子B的功能,可见式(5.1.22)中每个i值都对应1个无穷级数,q个无穷级数之
和代表了前面模型A与B的情况,它的可逆条件对应于无穷级数
2
2
1
i i
B B
b b
+ + +
的收敛条件,
即希望
1
1
ib
<或| |ib.于是有结论:
一个q阶的MA过程的可逆条件为:该过程的后移算子多项式( )Bθ的根,全部在B的平
均单位圆之外.
4.1.5 自回归过程自回归过程自回归过程自回归过程(AR模型模型模型模型)
假设{ }tZ是一个均值为0,方差为2
zσ的纯随机过程.如果有一个过程{ }tX满足:
1 1t t q t q tX X X Zα α- -= + + + (5.1.23)
则称该过程为p阶回归过程,记作( )AR p或称AR(Auto-Regressive Process)模型.这个模型很
象多元回归模型,但不是对于独立变量,而是对于tX的过去值的回归,故称之为自回归自回归自回归自回归.AR
模式在线性预测编码中起主角作用的一个线性模式,因此对它的讨论要深入一点.
1. AR(1)模型模型模型模型(Markov过程过程过程过程)
1阶自回归过程简写为AR(1)过程(模型),也叫Markov过程,即每个当前信号tX只与前
9
一时刻的信号值1tX-有关,其表达式为:
1t t tX X Zα-= + (5.1.24)
应用后移算子B(定义见式(5.1.17)),上式改写为:
(1 )t tB X Zα- =
因而有:
2 2
2
1 2
0
(1 )
1
t
tt
s
t t t t s
s
Z
X B B Z
B
Z Z Z Z
α α
α
α α α
∞
- - -
=
= = + + +
-
= + + + =∑
(5.1.25)
这就是MA(∞)的表达式.也就是说,1阶自回归过程可用无限多阶MA过程来表示.
随机变量tX的均值(1阶矩)为:
2[ ] [ ](1 )
1
z
t tE X E Z
μ
α α
α
= + + + =
-
| | 1α< (5.1.26)
在对Z作零均值假设条件下,可见[ ] 0tE X=.
2阶矩量有方差与自协方差,现分别计算如下:
2def
2 4 2
2
0 0
var[ ] var[ ](1 ) | | 1
1
cov[ , ] [ ]
z
t tx
s j
t t k t t k t s t k j
s j
X Z
X X E X X E Z Z
α
α α α α
α
α α
∞ ∞
+ + - + -
= =
= + + + = <
-
= =
==
∑ ∑
(5.1.27)
当j s k= +时,该两项相乘的均值为tZ的方差2
zσ;当j s k≠ +时,各项乘积的均值为0(因不
相关),故,
2
22
2
0
cov[ , ] ( )
1
k
s k s kz
t t k zz
s
X X k
α σ
γ σ α α α σ
α
∞
+
+
=
= = = =
-
∑ | | 1α< (5.1.28)
由此可见,1阶AR过程(Markov过程)的1阶矩与2阶矩都为常数,故该过程是2阶平
稳的,其条件为| | 1α
10
部分.由图可见,当α值减小时,自相关系数衰减加剧,当α为负值时,( )kρ变成交变状态,
即出负相关.负相关意味着两个随机变量中一个取较大值时另一个却取较小值.
从图中还可以看到AR模型的所谓"拖尾效应",即当k p≥(这里1p=)时,( )kρ并不为
0.虽然AR(1)模型的tX只与1tX-有关,但由于相关的传递性,1tX-又与2tX-有关,……,这
与MA模型的自相关系数的"截尾"特性相对应,构成了这两种模型的不同特点.
图图图图5.1.1 AR(1)过程的自相关函数例子过程的自相关函数例子过程的自相关函数例子过程的自相关函数例子
2. AR(p)过程过程过程过程
AR(p)过程的表达式见式(5.1.23),应用算子B仍可改写为:
1(1 )p
p t tB B X Zα α- - = (5.1.30)
我们所关心的问题,仍然是过程的平稳及过程中各样值之间的相关性.相关性问题将在5.4
节详细讨论,这里只讨论AR(p)过程的平稳条件.为此,令
1( ) (1 )p
pB B B α α= - - (5.1.31)
于是,式(5.1.30)可改写为:
( )t tB X Z = (5.1.32)
或者写成tX的表达式:
11
1( )t tX B Z -= (5.1.33)
所谓过程平稳的条件,就是要求表达式(5.1.33)有是限的,并且不是同时t的函数.由于式(5.1.31)
是一个以iα为系数的B的p次多项式,我们设该多项式有p个异实根1 2, , ,pb b b ,并进行类
似于式(5.1.21)与式(5.1.22)的讨论,即对( )B 作因式分解:
1
( ) ( )
p
i
i
B b B
=
= -∏ (5.1.34)
然后写出1( )B -为:
1
1
( )
p
i
ii
k
B
b B
-
=
=
-
∑ (5.1.35)
于是,式(5.1.33)可展开为:
2
2
1
1
p
i
tt
ii i i
kB B
X Z
b b b=
= + + +
∑ (5.1.36)
由此可见,过程tX平稳的条件,也就是上式中括号内无穷级数收敛的条件.显然,该条件就
是
1
1
ib
=的1 2,α α取值范围为:
2 1 2 1 2( ) 1,( ) 1, 1 1α α α α α+ < - < - < 时,用上述方法来确定AR过程的平稳域,会遇到计算上的困难,要寻
求更简单的判别方法.下面介绍常用一个方法.
3. Schur-Cohn准则准则准则准则
这里介绍的是一个比较简单的平稳性判别准则.由于从数据压缩应用出发来确定的模型
系数必须使系统能稳定工作,因此,我们主要关心的就是AR(p)模型的系数iα与系统稳定性
关系.
AR(p)模型所对应系数稳定的条件的充分必要条件就系统的特征多项式
1
1 1 0( )p p
p pF a a a aω ω ω ω-
-= + + + (5.1.38)
的根在复平面的单位圆内.这个结论可从( )B 的根的条件来推证.这里应该注意式(5.1.38)
与( )B 的区别,即i p iaα-=.
Schur-Cohn准则的关键是构造所谓的S-C行列式,它的定义为:
13
01 1
1 02
1 0
0 1 1
10 1
10
0
0
0
0
p p p k
p p k
kp
k
pk
p pk
p k p
a a a a
a a a a
a a a
a a a a
a a a a
a a a
- - +
- +
-
-
--
- +
=
,1, 2, ,k p= (5.1.39)
其中,ia为ia的共轭复数,且令01 = -.
S-C准则准则准则准则:( )Fω为稳定系统的特征多项式的充要条件是:它的S-C行行形式满足:
0
0k
k
k
为奇数
为偶数
(5.1.40)
显然,满足上述条件,意味着( ) 0Fω=的所有根都在单位圆内;否则,至少有一个根不
在单位圆内,但不能判定它们在圆外还是在圆上.
我们不去证明准则的正确性,只是举一个用它对AR(2)过程所对应的系统的稳定性判断
的例子,来说明它的应用.
AR(2)模型所对系统的特征多项式为:
2
2 1 0( )F a a aω ω ω= + +
S-C行列式为:
0 2
1
2 0
a a
a a
=
因为式(5.1.37)中的ia为实数,故这里i ia a=.然后用0
2 1 1 2, , 1a a aα α= - = - =代入上式,可得:
2
1 21α = -
由条件(5.1.40)应有10 <.因此,要求21 1α- +,不等式两边消去2
2(1 )α+后,
有:
2 2
2 1(1 )α α- >
得到稳定性条件为:2 11α α± 乘等式的两边,再各自求均值,便得到自协方差函数的关系式:
1 2( ) ( 1) ( 2) ( )pk k k k pγ α γ α γ α γ= - + - + + - (5.2.21)
相应的自相关函数的表达式为:
1 2( ) ( 1) ( 2) ( )pk k k k pρ α ρ α ρ α ρ= - + - + + - 1, 2 ,k p= (5.2.22)
这便是由p个方程组成的线性方程组,称之为Yule-Walker方程.这是一个在时间序列分析中
非常有用的方程,它反映了线性模型参数{ }iα与自相关系数{ }kρ之间的关系.写成矢量与矩
阵形式为:
=ρ αρ αρ αρ αR (5.2.23)
式中,
( (1), (2), , ( ))Tpρ ρ ρ= ρρρρ (5.2.24a)
1 2( , , , )T
pα α α= αααα (5.2.24b)
1 (1) ( 1)
(1) 1 ( 2)
( 1) ( 2) 1
p
p
p p
ρ ρ
ρ ρ
ρ ρ
-
- =
- -
R (5.2.24c)
可以看以,R是一个Toeplitz矩阵.
21
2. AR模型参数的估计模型参数的估计模型参数的估计模型参数的估计
若取得了信号的N个样值1 2, ,Nx x x ,它们是随机序列tX的一段样本值,则称N为样本
长度.为根据该段样本数据来估计tX,先要对自协方差( )kγ与自相关函数( )kρ作出估计.
我们采用:
1
1
^ ^( ) ( )
N k
l l k
l
k k x x
N
γ γ
-
+
=
= - =∑ 0,1, 2, , 1k N= - (5.2.25)
及
^( )
^ ^( ) ( )
^(0)
k
k k
γ
ρ ρ
γ
= - = 0,1, 2, , 1k N= - (5.2.26)
也有采用如下定义的:
* *
1
1
^ ^( ) ( )
N k
l l k
l
k k x x
N k
γ γ
-
+
=
= - =
-
∑ 0,1, 2, , 1k N= - (5.2.27)
从式(5.2.25)和(5.2.27)可以发现,
*^ ^( ) ( )
N k
k k
N
γ γ
-
=
因此,可见它们是渐近相等的,即当N→ ∞时有*^ ^( ) ( )k kγ γ=.
有了样本的自相关系数,便可由式(5.2.23)解出:
1^^^-=α ρα ρα ρα ρR (5.2.28)
展开得:
1
1
2
^^1 (1) ( 1) (1)
^^(1) 1 ( 2) (2)
^^( 1) ( 2) 1 ( )p
p
p
p p p
αρ ρ ρ
αρ ρ ρ
αρ ρ ρ
--
- =
- -
(5.2.29c)
常常称^αααα为αααα的Yule-Walker的估计.
为了得到残差能量的估计,将式(5.2.20)移项,两边平方后再取均值,得:
2 2
1 , 1
[ ] (0) 2 ( ) ( )
p p
z t j i j
j i j
E Z j j iσ γ α γ α α γ
= =
= = - + -∑ ∑
为了化简上式,可先把(5.2.23)按式(5.2.24)展开,并对第j个方程两边乘以jα;然后将方程组
的左边与右边分别加起来,这样就可得到:
22
1 , 1
( ) ( )
p p
j i j
j i j
j j iα γ α α γ
= =
= -∑ ∑ (5.2.30)
将上式代入2
zσ表达式,得:
2
1 , 1
(0) ( ) (0) ( )
p p
z j i j
j i j
j j iσ γ α γ γ α α γ
= =
= - = - -∑ ∑ (5.2.31)
将估计值^( )j iγ-与^
iα代入(5.2.31)便得到2
zσ的估计值2^
zσ:
2
1 , 1
^ ^ ^ ^ ^ ^ ^^(0) ( ) (0) ( )
p p
z j i j
j i j
j j iσ γ α γ γ α α γ
= =
= - = - -∑ ∑ (5.2.32)
式(5.2.29)与式(5.2.32)就是AR模型参数矩估计公式.可以发现它们是比较简单易行的.
下面用AR(1)与AR(2)模型的矩估计来说明.
AR(1)参数估计:
1
^(1)
^ ^(1)
^(0)
γ
α ρ
γ
= = (5.2.33a)
22
1
^ ^ ^ ^ ^(0) (1) (0)[1 (1) ]zα γ γ α γ ρ= - = - (5.2.33b)
AR(2)参数估计:
12
^^(1)[1 (2)]
^
^[1 (1) ]
ρ ρ
α
ρ
-
=
-
(5.2.34a)
2
22
^ ^[ (2) (1) ]
^
^[1 (1) ]
ρ ρ
α
ρ
-
=
-
(5.2.34b)
2
1 2
^ ^ ^ ^ ^(0)[1 (1) (2)]zα γ α ρ α ρ= - - (5.2.34c)
其中的^( )kγ与^( )kρp(k)由式(5.2.25)与(5.2.26)给出.
前面已指出AR模型的矩估计与最小二乘估计是很接近的.虽然Yule-Walker方程的导出
过程中并未引入"最小误差"的条件.但若以最小均方误差条件来推导参数估计关系,也得
到相同的结果.若以相同的样本序列(设N个样值)来作以上的两种估计,比如作出AR(1)或
AR(2)估计,就可发现其结果相差甚微,只是求和限有点差别.当N→ ∞,这个差别可略.
23
5.2.3 序列最小二乘法序列最小二乘法序列最小二乘法序列最小二乘法
上两节介绍了线性预测参数的两种估计方法.从原理上说.这两种方法是简单的而且有
效的,但从计算速度与实时处理这两方面来说,并不理想.这是因为在估值过程中需要计算
如方程式(5.2.12)与(5.2.28)中的逆矩阵,而逆矩阵的计算往往是复杂费时的;这两种方法都是
属于成批型的算法,只有在取得足够多(N足够大)的样值时,估值才有效,这对数据的实时
处理要求很不利,且要求设备有较大的存储容量.为克服上述困难,现在来介绍运用自然递
批的方法.本节讨论的是称作序列最小二乘(SLS)算法.
1. 最小二乘估计的递推结构最小二乘估计的递推结构最小二乘估计的递推结构最小二乘估计的递推结构
我们从最小二乘估计的结果出发来推导自然递推关系.为了避开矩阵的求逆运算.先对
式(5.2.12)作一些改造.为此,将式(5.4.11)的左右两边分别改写成:
1,
1
1,1 1,2 1,1
2,1 2,2 2,2 ,
1
,1 ,2 ,
,
1
1 1,1
1 2,1
1 ,1
n
t i ti
i
t t t ntn
t t t nt t i tiT
it
t p t p t p ntnn
t p i ti
i
t t
t t
t t p
x x
x x xx
x x xx x x
x x xx
x x
x x
x x
x x
-
=
- - -
- - --
=
- - -
-
=
-
-
-
= =
=
∑
∑
∑
x x
2 1,2 1,
2 2,2 2,
2 ,2 ,
t t tn t n
t t tn t n
t t p tn t p n
x x x x
x x x x
x x x x
- -
- -
- -
+ + +
若令
1, 2, ,[ , , , ]T
i t i t i t p ix x x- - -= x (5.4.35)
则有,
1
n
T
t ti i
i
x
=
=∑x x x (5.2.36)
以及
24
2
1, 1, 2, 1, ,
1 1 1
2
2, 1, 2, 2, ,
1 1 1
2
, 1, , 2, ,
1 1 1
n n n
t i t i t i t i t p i
i i i
n n n
t i t i t i t i t p iT
i i it
n n n
t p i t i t p i t i t p i
i i i
x x x x x
x x x x x
x x x x x
- - - - -
= = =
- - - - -
= = =
- - - - -
= = =
=
∑ ∑ ∑
∑ ∑ ∑
∑ ∑ ∑
x x
这是一个对称矩阵,故在式(5.2.16)的推证中1( )T
t
-x x的转置并不作倒序排列.若应用式(5.2.35)定义
的向量ix,则上式可简化为:
def
1
1
n
T T
t i i n
i
-
=
=∑==x x x x P (5.2.37)
于是式(5.2.11)与式(5.2.13)分别可写成:
1
1
^
n
n ti i
i
x-
=
=∑ααααP x (5.2.38)
与
1
^
n
n ti i
i
x
=
=∑ααααP x (5.2.39)
为了得到自然迭代关系,将式(5.2.38)左边分成两项,即,
1
1
1
^
n
n ti i tn n
i
x x
-
-
=
= +∑ααααP x x (5.2.40)
用1n-去代替n,式(5.2.38)仍应成立,即,
1 1
1
1 1
^
n n
T
i i n- ti i
i i
x
- -
= =
=
∑ ∑ααααx x x (5.2.41)
这里用1
^
n-αααα来表示1n-组观测值得到的参数估计,以区别于n组观测值得到的估计^
nαααα.下标n
也是为了区别才加上去的.应该强调的是,^
nαααα与前面公式中的^αααα是完全一样.
把式(5.2.41)代入式(5.2.40),便有:
1
1
1
1
^ ^
n
T
n n i i n- tn n
i
x
-
-
=
= +
∑α αα αα αα αP x x x (5.2.42)
若对式(5.2.42)左边同时加减一项1
^T
n n n-ααααx x,并应用1
n
-P的定义式(5.2.37),便得递推关系:
25
1 1
1 1
^ ^^( )T
n n n n- n tn n n-x- -= + -α α αα α αα α αα α αP P x x (5.2.43)
或者:
1 1
^ ^^( )T
n n- n n tn n n-x= + -α α αα α αα α αα α αP x x (5.2.44)
从上面的递推关系可以看到,估值^
nαααα可以由1
^
n-αααα,也就是它前面的估值项(在减少一个数
据集的基础上)以及一个修正项来计算.这个修正项是使用前面的估值引起的预测误差
1
^T
tn n n-x-ααααx的一个函数.应该指出,尽管这是个递推形式,但它仍然属于最小二乘估值算法,
因为上面的一切推导都是由最小二乘估值结果出发的.同时应注意式(5.2.44)远不理想,因为
从1
n
-P的定义来计算nP.仍然有逆矩阵运算问题.下面来解决这个问题.
2. 逆矩阵定理逆矩阵定理逆矩阵定理逆矩阵定理
所谓逆矩阵定理,就是不用计算逆矩阵而得到矩阵求逆的结果的一个递推算法.
定理定理定理定理 由式(5.2.37)定义的1
n
-P,它的逆只需由1n-P及nx(由式(5.2.35)所定义)按下述关系
递推得出,即:
1 1
1
11
T
n n n n
n nT
n n n
- -
-
-
= -
+
P x x P
P P
x P x
(5.2.45)
上式的证明这里不再讲述.
由式(5.2.44)与(5.2.45)构成了一组递推条件,并称此为序列最小二乘估计算法.该算法的
关键是用由式(5.2.35)定义的向量ix作为已知条件代原来已知矩阵x.而ix实际上是分别由各
组观测值组成的,从而构成一组一组递推估计的流程,即由第一组(1i=)观测值来作αααα的第一
—次估计1
^αααα;然后由1
^αααα与第二组(2i=)观测值作第二次估计2
^αααα;….那么流程的初值即作第
一次递推估计时的0
^αααα与0P应取什么值呢 下面分别加以说明.
3. 初始值初始值初始值初始值0
^αααα与与与与0P
初始值0
^αααα与0P是递推算法所必需的已知条件.然而在我们讨论的问题中,它们却是先验
未知的,所以还需对它们进行估值.
当对^αααα缺乏先验知识时,初始参数的估值0
^αααα可以设置为零向量:
0
^=0αααα (5.2.46)
26
关于0P的初始化值可以设置为一对角矩阵:
0β=P I 1β>> (5.2.47)
这里I为单位矩阵.为什么β要比1大很多 这是这样考虑的:从定义式(5.2.37)的求和限来看,
n最小为1,也就是说,1
0
-P应为0矩阵才行.但从迭代关系式(5.2.45)来看,对0P的任何实际
可行的估值都必须是有限的.因此,设置1
0ε-=P I和1ε,所有的( )k
kα都应为零,即偏相关函数是"截尾"的.还可以进一步证
明偏相关函数的截尾性也是平稳时间序列为AR模型的充分条件,即它是AR过程所持有的
标志.
利用AR过程的截尾特性,可以作出结论:"最佳"的阶数可被确定为这样的k值,这种
k值之后的{ }pπ的非常接近于0,此时p k=.
3. 与与与与Durbin递推程序同时进行判决递推程序同时进行判决递推程序同时进行判决递推程序同时进行判决
在Durbin递推算法公式(5.2.49)和(5.2.51)中,我们定义:
( ) ( )
1
(0) ( )
k
k k
j
j
E r r jα
=
= -∑ (5.2.74)
若信号属于AR(p)过程,是可以证明:
( )limp
N
N
E E
→∞
=
而NE定义为均方差(N为样值序值的样值个数)
2
1
0 1
1pN
N k i k i
k i
E x x
N
α
-
-
= =
= -
∑ ∑ (5.2.75)
从式(5.2.48),有1k=时,
(1)
1 1
(1)
(0)
K
γ
α
γ
= =
由式(5.2.75)定义,可得:
(1) (1)
1(0) (1)Eγ α γ= -
2k=时,有:
35
( 2) 2 (1)
2[1 ]E K E= -
此结果表时,由式(5.2.74)定义的( )kE也是具有递推特性.因此,算法式(5.2.51)可改写成:
1
( 1)
1( )
( 1)
( ) ( )
k
k
j
jk
k kk
k k j
K
E
γ α γ
α
-
-
=
-
- -
= =
∑
(5.2.76a)
( ) ( 1) ( 1)k k k
k j k k jKα α α- -
-= - (5.2.76b)
( ) 2 ( 1)[1 ]k k
kE K E-= - (5.2.76c)
式中,1, 2, , ; 1, 2, , 1k p j k= = - .
递推初值为:
(0 )(0)Eγ= (5.2.77a)
(1)
1 1
(1)
(0)
K
E
γ
α= = (5.2.78b)
(1) 2 ( 0)
1(1 )E K E= - (5.2.78c)
这是Durbin递推算法的另一种表示.我们的目的是借助于递推过程中的( )kE进行最佳阶数的
判定.判定准则有多种,这里只介绍用得较多的三种(它们之间并没有原则的差别).
(1) FPE准则
最终预测误差准则FPE(Final Prediction Error Criterion)是把使下式为最小值的k作为AR
模型的最佳阶数:
( )1
FPE( )
1
kN k
k E
N k
+ +
=
- -
(5.2.78)
般情况下,( )kE是k的单调递减函数,但系数
1
1
N k
N k
+ +
- -
是k的递增函数. 于是会在k的某个
值,FPE(k)出现最小值,这个k值就被定为AR模型的最佳阶数.
〔2) AIC准则
信息理论准则AIC (1nformtion Theoretic Criterion)是把使下式为最小值的k,用为AR模型
的最佳阶数:
( )2
AIC( ) ln[ ]kk
k E
N
= + (5.2.79)
式(5.2.79)与式(5.2.78)实质上是一样的.因为当N很大时,FPE的对数收敛于AIC,即有:
36
lim[ln FPE( )] AIC( )
N
k k
→∞
= (5.2.80)
这个极限关系是很易验证的.
(3) CAT准则
自回归传递函数准则CAT(Autoregressive Transfer Function Criterion)是把使下式为最小值
的k作为AR模型的阶数:
( ) ( )
1
1
CAT( )
k
j k
j
N k N k
k
N NE NE=
- -
=
∑ (5.2.81)
它相当于在估计滤波器的滤波误差的均方值上给出最佳阶数.
应该说,AR模型的准确阶数的判定,对于谱分析来说,比线性预测应用更为重要.在线
性预测中,阶数的不准确,引起的误差并不会大大,只要保持自协方差函数相差不太大;即
可.这就要特别注意在计算样本的自协方差函数,偏相关函数以及模型识别与参数化估计时,
—定要设法将均值项去掉.去均值项的常用办法是从ix中减去它的样本均值:
1
1
^
N
i
i
x
N
μ
=
=∑
即用^
t ty xμ= -作tx参数估计.但如果tx 的期望值( )tE X本身为零,再减去^μ就会降低估计
精度,可见.在对数据序列1 2, , ,Nx x x 作时序分析时,判断( )tE X是否为0是很有必要的.
5.2.7 窄带平稳随机过程窄带平稳随机过程窄带平稳随机过程窄带平稳随机过程的的的的MATLAB实现实现实现实现
这里窄带平稳随机过程定义为信号功率谱密度是一个带通型且带宽远小于中心频率的随
机过程,通常窄带平稳随机过程( )X t可以展开如下形式:
2( ) ( )cos 2 ( )sin 2 Re[ ( ) ]cj f t
c c s c lX t X t f t X t f t X t eππ π= - = (5.2.82)
其中,( ) ( ) ( )l c sX t X t jX t= +是等效基带信号,( ), ( )c sX t X t分别是平稳随机,其功率谱密度是
基带型.一个典型的窄带平稳随机过程是高斯白噪声经过带通系统的输出,可以证明,窄带
平稳随机过程具有如下性质:
(1) 如果( )X t平稳,则( ), ( )c sX t X t也平稳;
(2) 如果[ ( ) 0E X t=,则[ ( )] [ ( )] 0c sE X t E X t= =;
37
(3) 如果( )X t是高斯平稳过程,则( ), ( )c sX t X t也是高斯平稳过程;
(4) 如果( )X t是高斯平稳过程用均值为0,则( ), ( )c sX t X t相互正交,且功率与( )X t相同.
例例例例5.2.2 高斯白噪声是一个功率谱密度为常数,且各时刻满足高斯分布的随机过程,设
高斯白噪声的双边功率谱密度为0/ 2N,现将高斯白噪声经过一个中心频率为cf,带宽为B的
带通滤波器,得到的输出信号即为窄带随机过程.
(1) 用Matlab产生高斯平稳窄带随机过程,010, 1Hz, 1fc B N= = =;
(2) 画出窄高斯随机过程的等效基带信号;
(3) 求出窄高斯随机过程的功率,等效基带信号的实部功率,虚部功率.
解解解解::::窄高斯随机过程的Matlab的程序为:zppw.m
*********************************************************
clear all
close all
N0=1; % 单边功率谱密度
fc=10; % 中心频率
B=1; % 带宽
dt=0.01;
T=100;
t=0:dt:T-dt;
% 产生功率功率为N0*B的高斯白噪声
P=N0*B;
st=sqrt(P)*randn(1,length(t));
% 将上述白噪声经过窄带带通系统
[f,sf]=T2F(t,st); % 高斯信号频谱
figure(1)
plot(f,abs(sf)) % 高斯信号的幅频特性
38
[tt,gt]=bpf(f,sf,fc-B/2,fc+B/2); % 高斯信号经过带通系统
glt=hilbert(real(gt)); % 窄带信号的解析信号,调用hilbert函数
glt=glt.*exp(-j*2*pi*fc*tt); % 得到解析信号
[ff,glf]=T2F(tt,glt);
figure(2)
plot(ff,abs(glf));
xlabel('频率 (Hz)')
ylabel('窄带高斯过程样本的幅频率特性')
figure(3)
subplot(411)
plot(tt,real(gt));
title('窄带高斯过程样本')
subplot(412)
plot(tt, real(glt).*cos(2*pi*fc*tt)-imag(glt).*sin(2*pi*fc*tt))
title('由等效基带重构的窄带高斯过程样本')
subplot(413)
plot(tt, real(glt))
title('窄带高斯过程样本的同相分量')
subplot(414)
plot(tt,imag(glt));
xlabel('时间 t(秒)')
title('窄带高斯过程样本的正交分量')
% 求窄带高斯信号功率;注:由于样本的功率近似等于随机过程的功率,因此可能出现一些偏差
P_gt=sum(real(gt).^2)/T;
P_glt_real=sum(real(glt).^2)/T;
39
P_glt_imag=sum(imag(glt).^2)/T;
% 验证窄带高斯过程的同相分量,正交分量的正交性
a=real(glt)*(imag(glt)')/T;
********************************************************
用到的子函数bpf.m如下.
*******************************************************
function [t,st]=bpf(f,sf,B1,B2)
df=f(2)-f(1);
T=1/df;
hf=zeros(1,length(f));
bf=[floor(B1/df):floor(B2/df)];
bf1=floor(length(f)/2)+bf;
bf2=floor(length(f)/2)-bf;
hf(bf1)=1/sqrt(2*(B2-B1));
hf(bf2)=1/sqrt(2*(B2-B1));
yf=hf.*sf.*exp(-j*2*pi*f*0.1*T);
[t,st]=F2T(f,yf);
*******************************************************
用到的子函数F2T.m如下.
*******************************************************
function [t,st]=F2T(f,sf)
df=f(2)-f(1);
Fmx=(f(end)-f(1)+df);
dt=1/Fmx;
N=length(sf);
T=dt*N;
40
t=0:dt:T-dt;
sff=fftshift(sf);
st=Fmx*ifft(sff);
*******************************************************
程序运行结果如图5.2.2所示.
0102030405060708090100
-0.5
0
0.5
窄窄窄窄窄窄窄窄
0102030405060708090100
-0.5
0
0.5
由由由由窄由由由窄窄窄窄窄窄窄窄
0102030405060708090100
-0.5
0
0.5
窄窄窄窄窄窄窄窄由窄窄窄窄
0102030405060708090100
-0.5
0
0.5
时时 t(秒)
窄窄窄窄窄窄窄窄由窄窄窄窄
图图图图5.2.2 窄带高斯过程窄带高斯过程窄带高斯过程窄带高斯过程
5.3 通用的预测编通用的预测编通用的预测编通用的预测编码码码码(量化编码方法量化编码方法量化编码方法量化编码方法)
音频类信源(包括语音,音乐,广播等音频信号在内)的特点是发出连续信息,因此其编
码方法与字符信源有所不同.首先要在时间上进行抽样离散化过程,然后在数值(幅度)范围
内进行量化,量化后的信号就是多进制的数字化信息;最后才进行数字信号编码,其输出的
二进制代码(符号)在信道内传输.简单地说,对发出连续信息的音频信源要经过抽样,量化
和编码三个过程,才能进行信息传输.在这些信号处理过程中都会有信息受到损失,而这些
信息的损失又是无法弥补的,因此,对于音频信源来说,信源编码是属于限失真编码,其编
码效率受到信息率失真函数( )R D的限制,无法实现像字符信源那样的最佳编码效果.
抽样与量化,我们在第二章已经讨论过了,限失真编码原理在第三章讨论过了,这里我
们根据上述基础与原理,来讨论量化的编码方法.
41
5.3.1 脉冲编码调制脉冲编码调制脉冲编码调制脉冲编码调制
前面已经指出,模拟信息源输出的模拟信号经抽样和量化后得到的输出电平序列
{ ( )}q sm kT,可以利用M进制PAM直接进行传输,也可以将每一个量化电平用编码方式传输.
所谓编码是把量化后的信号变换成代码,其相反的过程称为译码.
将模拟信号的抽样量化值变换成代码,称为脉冲编码调制(PCM).图5.3.1和表5.3.1给
出了脉冲编码调制的一个实例.假设模拟信号m(t)的最大值| ( ) | 4m t时,出"1"码;反之出"0"码.由于在
13折线法中用了7位二进代码来代表段落和段内码,所以对一个输入信号的抽样值需要进行
7次比较.每次所需的标准电流wI均由本地译码电路提供.
本地译码电路包括记忆电路,7/11变换电路和恒流源.记亿电路用来寄存二进代码,
因除第一次比较外,其余各次比较都要依据前几次比较的结果来确定标准电流wI值.因此,
7位码组中的前6位状态均应由记忆电路寄存下来.7/11变换电路就是前面非均匀量化中谈
到的数字压缩器.因为采用非均匀量化的7位非线性编码等效于11位线性码,而比较器只能
编7位码,反馈到本地译码电路的全部码也只有7位.因为恒流源有11基本权值电流支路,
需要11个控制脉冲来控制,所以必须经过变换,把7位码变成11位码,其实质就是完成非
线性和线性之间的变换.恒流源用来产生各种标准电流值.为了获得各种标准电流wI,在恒
流源中有数个基本权值电流支路.基本的权值电流个数与量化级数有关,如上例中,128个
量化级需要编7位码,它要求11个基本的权值电流支路,每个支路均有一个控制开关.每次
该哪几个开关接通组成比较用的标准电流wI,由前面的比较结果经变换后得到的控制信号来
控制.
附带指出,保持电路的作用是保持输入信号的抽样值在整个比较过程中具有一定的幅度.
由于逐次比较型编码器编7位码(极性码除外)需要将sI与wI比较7次,在整个比较过程中部
应保持输入信号的幅度不变,故需要采用保持电路.下面我们通过一个例子来说明编码过程.
例例例例5.3.1 设输入信号抽样值为+1270个量化单位(这里量化单位是指以输入信号归一化的
1/2048为单位),采用逐次比较型编码将它按照13折线A律特性编成8位码.
设码组的8位码分别用1 2 3 4 5 6 7 8c c c c c c c c表示.编码过程如下:
(1) 确定极性码1c
因输入信号抽样值为正,故极性码11c=
(2) 确定段落码2 3 4c c c
13折线中,正半部分的8个段落以1/2048为单位的每个段落的起点电平如表5.3.5所
示.由于段落码中的2c是用来表示输入信号抽样值处于8个段落的前四段还是后四段的,故
输入比较器的标准电流应选择为128wI=个量化单位.现在输入信号抽样值1270sI=个量化单
47
位,大于标准电流,故第一次比较结果为s wI I>J,所以21c=.它表示输入信号抽样值处于8
个段落中的后四段(5~8段).
表表表表5.3.5 段落起点电平段落起点电平段落起点电平段落起点电平
3c用来进一步确定它属于5~6段还是7~8段.因此,标准电流应选择为512wI=个量化
单位.第二次比较结果为s wI I>,故31c=.它表示输入信号属于7~8段.
同理,确定4c的标准电流应为1024wI=个量化单位.第三次比较结果为s wI I>,故41c=.
由以上三次比较得段落码为"111",输入信号抽样值1270sI=个量化单位应属于第八段.
(3) 确定段内码5 6 7 8c c c c
由编码原理已经知道,段内码是在已经确定输入信号所处段落的基础上,用来表示输入
信号处于该段落的哪一量化级的.5 6 7 8c c c c的取值与量化级之间的关系见表5.3.4.上面已经确
定输入信号处于第八段,该段中的16个量化级之间的间隔均为64个量化单位,故确定5c的
标准电流应选为:
8 ( 1024 8 64 1536wI= + × = + × =段落起点电平 量化级间隔) 个量化单位
第四次比较结果为s wI I,故81c=.它说明输入信号处于第八段中的3量化级.
经上述七次比较,编出的8位码为11110011.它表示输入抽样位处于第八段3量化级,
其量化电乎为1216个量化单位,故量化误差等于54个量化单位.顺便指出,除极性码外的
7位非线性码组1110011,相对应的11位线性码组为10011000000.
上面较详细地讨论了逐次比较型编码器的原理,下面再来讨论PCM信号的译码原理.常
用译码器大致可分为三种类型:电阻网络型,级联型,级联一网络混合型等.这里仅就图5.3.4
所示的电阻网络型译码器加以讨论.
图图图图5.3.4 电阻网络型译码器电阻网络型译码器电阻网络型译码器电阻网络型译码器
电阻网络型译码器与逐次比较型编码器中的本地译码器基本相同.从原理上说,两者都
是用来译码,但编码器中的译码,只译出信号的幅度,不译出极性;而收端的译码器在译出
信号幅度值的同时,还要恢复出信号的极性.电阻网络型译码器各部分电路的作用简述如下.
记忆电路用来将接收的串行码变为并行码,故又称为"串/并变换"电路.7/11变换
电路用来将表示信号幅度的7位非线性码转变为11位线性码.极性控制电路用来恢复译码后
的脉冲极性.寄存读出电路把寄存的信号在一定时刻并行输出到恒流源中的译码逻辑电路上
去.使产生所需要的各种逻辑控制脉冲.这些逻辑控制脉冲加到恒流源的控制开关上,从而
驱动权值电流电路产生译码输出.
由以上电阻网络则译码器各部分电路的作用.不难理解这种译码器的工作原理.其译码
过程就是根据所收到的码组(极性码除外)产生相应的控制脉冲去控制恒流源的标推电流支
路,从而输出—个与发送端原抽样值接近的脉冲.该脉冲的极性受极性控制电路控制.
2. PCM系统的抗噪声性能系统的抗噪声性能系统的抗噪声性能系统的抗噪声性能
上面我们较详细地讨论了脉冲编码调制的原理.现在,将要分析图5.3.2所示的PCM系
统的抗噪声性能.由该图可以看出,接收端低通滤波器的输出为:
49
0
^( ) ( ) ( ) ( )q em t m t n t n t= + +
式中:0( )m t —— 输出信号成分;
( )qn t —— 由量化噪声引起的输出噪声;
( )en t —— 由信道加性噪声引起的输出噪声.
为了衡量PCM系统的抗噪声性能,通常将系统输出端总的信噪比定义为:
2
0 0
2 2
0
[ ( )]
[ ( )] [ ( )]q e
S E m t
N E n t E n t
=
+
(5.3.1)
式中 E —— 求统计平均.可见,分析PCM系统的抗噪声性能时,需要考虑量化噪声和
信道加性噪声的影响.不过由于量化噪声和信道加性噪声的来源不同,而且它们互不依赖,
故可以先讨论它们单独有在时的系统性能,然后再分析系统总的抗噪声性能.我们先分析仅
考虑量化噪声的系统性能.假设发送端采用理想冲激抽样,则PCM系统输出端平均信号量化
噪声功率经为:
2
20 0
2
[ ( )]
[ ( )]q q
S E m t
M
N E n t
= = (5.3.2)
式中,qN为量化噪声功率,M量化级数.
对二进制编码,式(5.3.2)又可写成:
202N
q
S
N
= (5.3.3)
式中,N—— 二进制代码位数.由上式可见,PCM系统输出端平均信号量化噪声功率比将
仅依赖于每一个编码组的位数N.上述比值将随N按指数增加.大家知道,对于一个频带限
制在Hf的信号,按照抽样定理,此时要求每秒钟最少传输的抽样脉冲数等于2Hf;若PCM
系统的编码位数为N,则要求系统每秒传输2HNf.为此,这时的系统总带宽B至少等于HNf.
故式(5.3.3)还可写成:
2
02H
B
f
q
S
N
= (5.3.4)
由此可见,PCM系统输出端的信号量化噪声功率比还与系统带宽B成指数关系.
下面我们来分析信道加性噪声对PCM系统性能的影响.由于信道中加性噪声对PCM信
50
号的干扰,将造成接收端判决器判决错误,二进制"1"码可能误判为"0"码,而"0"码可
能误判为'1"码,其错误概率将取决于信号的类型和接收机输入端的平均信号功率比.因为
PCM信号中每一码组代表着一定的量化抽样值,所以其中只要发生误码,接收端恢复的抽样
值就会与发端原抽样值不同.通常我们只需要考虑仅有一位错码的码组错误,而多于一个错
码的码组错误可以不予考虑.
仅考虑信道加噪声时PCM系统的输出信噪比为:
01
4e e
S
N P
= (5.3.5)
式中,eP为误码率,由上式可见,由误码引起的信噪比与误码率成反比.
前面已经指出,传输模拟信号的PCM系统的性能用接收端输出的平均信器噪功率比来度
量.将式(5.3.2),(5.3.3)和(5.3.5)代入式(5.3.11)得:
22 2
0 0
2 2 2 2
0
[ ( )]2
[ ( )] [ ( )] 1 4 2 1 4 2
N
N N
q e e e
S E m tM
N E n t E n t P P
= = =
+ + +
(5.3.6)
在接收端输入大信噪比条件下,即当24 2 1N
eP<>时,
2
0
2
0
2 1
4 2 4
N
N
e e
S
N P P
≈ = (5.3.8)
在用基带传输的PCM中继系统中,通常使误码率降到610-是容易实现的,此时,可按式(5.3.8)
来估计PCM系统的性能.
5.3.2 增量调制增量调制增量调制增量调制
增量调制( M)是在PCM方式的基础上发展起来的另一种模拟信号数字传输的方法. M
可以看成PCM的一个待例,因为它们都是用二进制代码形式去表示模拟信号的方式.但是,
在PCM中,信号的代码表示模拟信号的抽样值,而且,为了减小量化噪声,一般需要较长的
代码及较复杂的编译码设备.而 M是将模拟信号变换成仅由一位二进制码组成的数字信号
序列,并且在接收端也只需要用一个线性网络,便可复制出原模拟信号.因而, M有它自
51
己的特点,而且其编译码设备通常要比PCM的简单.
1. 增量调制原理增量调制原理增量调制原理增量调制原理
我们知道,一位二进制码只能代表两种状态,当然就不可能去表示抽样值的大小.可是,
用一位码却可以表示相邻抽样值的相对大小,而相邻抽样值的相对变化将能同样反映模拟信
号的变化规律.因此,由一位二进制码去表示模拟信号的可能性是存在的.为了确信这一点,
让我们通过下而的例子来说明.设一个频带有限的模拟信号如图5.3.5中的( )m t所示.现在把
横轴t分成许多相等的时间段t .此时可以看出,如果t 很小,则( )m t在间隔为t 的时刻上
得到的相邻值的差别(差值)也将很小.因此,如果把代表( )m t幅度的纵轴也分成许多相等的
小区间σ,那么,一个模拟信号( )m t就可用如图5.3.5所示的阶梯波形'( )m t来逼近.显然,
只要时间间隔t 和台阶σ都很小,则( )m t和'( )m t将会相当地接近.由于阶梯波形相邻间隔
上的幅度差不是σ+就是σ-.,因此,倘着用二进制码的"1"代表'( )m t在给定时刻上升一个
台阶σ,用"0"表示'( )m t下降一个台阶σ,则'( )m t就被一个二进码的序列所表征(见图5.3.5
中横轴下面的二进码序列).于是,该序列也相当于表征了( )m t.
图图图图5.3.5 增时调制波形示意图增时调制波形示意图增时调制波形示意图增时调制波形示意图
在讨论怎样得到发端的阶梯波形及由此波形又如何确定二进码序列之前,我们先讨论一
下在接收端怎样由二进码序列恢复出阶梯波形的问题.即 M信号的译码问题.不难看出,
接收端只要每收到一个"1"码就使输出上升一个值σ,每收到一个"0"码就使输出下降一
个σ值,连续收到"1"码(或"0"码)就使输出一直上升(成下降),这样就可以近似地复制出
阶梯波形'( )m t.这种功能的译码器可由一个积分器来完成,如图5.3.6所示.积分器遇到"1"
码(即有E+脉冲电压),就以固定斜率上升一个E ,并让E 等于σ,遇"0"码所表示的E-,
就以同样的斜率下降一个E .图5.3.6(b)表示了积分器的输入和输出波形.由此看到,积分
器的输出波形并不是阶梯波形,而是一个斜变波形.但因Eσ =,故在所有抽样时刻it上斜
52
变波形与阶梯波形有完全相同的值.因而,斜变波形同样与原来的模拟信号相似.由于积分
器实现起来容易,且能符合译码要求,故常被采用.最简单的积分器就是图5.3.6(c)所示的
RC积分器,其RC的乘积应远大于一个二进码的脉冲宽度.积分器输出虽已接近原来模拟信
号,但往往还包含有不必要的高次谐波分量,故需再经低通滤波器平滑,这样,便可得到十
分接近模拟信号的输出信号.
图图图图5.3.6 积分器译码示意积分器译码示意积分器译码示意积分器译码示意
现在回过来讨论, M的编码原理.一个简单的 M编码器组成如图5.3.7所示,它由相
减器,判决器,本地译码器及抽样脉冲产生器(脉冲源)组成.本地译码器与接收端的译码器
完全相同.判决器将在抽样脉冲到来时刻对输入信号的变化作出判决,并输出脉冲.这个编
码器的工作过程如下:将模拟信号( )m t与本地译码器输出的斜变波形'( )m t进行比较,为了获
得这个比较结果,先进行相减,然后在抽样脉冲作用下将相减结果进行极性判决.如果在给
定抽样时刻it有:
( ) | '( ) | 0
i it t t tm t m t
- -= =- >
则判决器输出"1"码;如有:
( ) | '( ) | 0
i it t t tm t m t
- -= =- <
则发"0"码.这里,it-是it时刻的前一瞬间,即相当于在阶梯波形跃变点的前一瞬间.于是,
编码器将输出一个如图5.3.5所示的二进码序列.
53
图图图图5.3.7 M编码器编码器编码器编码器
从上述讨论可以看出, M信号是按台阶σ来量化的(增,减一个σ值),因而同样存在
量化噪声问题. M系统中的量化噪声有两种形式;一种称为过载量化噪声,另一种称为一
般量化噪声,如图5.3.8所示.过载量化噪声(有时简称过载噪声)发生在模拟信号斜率陡变时,
由于台阶σ是固定的,而且每秒内台阶数也是确定的,因此,阶梯电压波形就跟不上信号的
变化,形成了很大失真的阶梯电压波形,这样的失真称为过载现象,也称过载噪声,如图5.3.8(b)
所示;如果无过载噪声发生,则模拟信号与阶梯波形之间的误差就是一般的量化噪声,如图
5.3.8(a)所示.图中的( ) ( ) '( )n t m t m t= -,可统称其为量化噪声.
图图图图5.3.8 两种形式的量化噪声两种形式的量化噪声两种形式的量化噪声两种形式的量化噪声....(a) 一般量化噪声一般量化噪声一般量化噪声一般量化噪声;;;;(b) 过载量化噪声过载量化噪声过载量化噪声过载量化噪声
设抽样时间间隔为t (抽样频率1/sf t= ),则一个台阶上的最大斜率K为:
sK f
t
σ
σ= =
这也就是译码器的最大跟踪斜率.当信号实际斜率超过这个最大跟踪斜率时,则将造成过载
噪声.因此,为了不发生过载现象,则必须使sf和σ的乘积达到一定的数值,以使信号实际
斜率不超过这个数慎.这个数值通常可以增大sf或σ来达到.
54
对于一般量化噪声,由图5.3.8(a)不难看出,σ小大则这个量化噪声大,σ小则噪声小.
用大的σ虽然能减小过载噪声,但却增大了一般量化噪声.因此,σ值应适当选取.
不过,我们看到, M系统的抽样频率必须选得足够高,因为这样,既能减小过载噪声,
又能降低一般量化噪声,从而使 M 系统的量化噪声减小到给定的容许数值.一般, M 系
统中的抽样频率要比PCM系统的抽样频率高得多(通常要高两倍以上).
2. M 系统中的量化噪声系统中的量化噪声系统中的量化噪声系统中的量化噪声
图5.3.9(a)是 M系统的组成方框图,它包括增量调制器,信道,检测器,积分器(译码
器)相低通滤波器等.为了简明起见,我们将图(a)中0( ), '( ), ( ), ( )qm t m t p t e t的波形画于图(b),并
设系统输出的信号和量化噪声分别用( )m t和( )qn t表示.如果信道的加性噪声足够小,以致不
造成误码,那么,用于接收的检测器将输出一个与0( )p t完全相同的信号,也即0 0' ( ) ( )p t p t=,
此时,系统的输出信号0( )m t与( )m t将有最好的近似(因为,量化噪声仍然存在);如果信道噪
声造成了误码,则在系统的输出噪声中不仅存在量化噪声,而且还存在由误码引起的噪声.
图图图图5.3.9 M系统的成及有关点的波形系统的成及有关点的波形系统的成及有关点的波形系统的成及有关点的波形
55
我们分析存在量化噪声时的系统性能.此时认为信道加性噪声很小,不造成误码.那么,
接收端检测器输出的0' ( )p t就是0( )p t,而解调积分器输出端的信号便是'( )m t(见图5.3.928).
容易看出,在这个积分器输出端的误差波形正是量化误差波形( )qe t.因此,如果求得( )qe t的
平均功率,则系统的输出量化噪声功率也就可以确定.我们还观察到,只要 M 系统不发生
过载现象(过载现象在设计时是需要克服的),那么( )qe t总是不大于σ±的.我们假设随时间随
机变化的( )qe t在区间( , )σ σ- +上均匀分布,于是( )qe t的一维概率密度( )qf e可表示为:
1
( )
2qf e
σ
= eσ σ- ≤ ≤ +
因而( )qe t的平均功率可表示成:
2
2 2 21
[ ( )] ( )
2 2q qE e t e f e de e de
σ σ
σ σ
σ
σ
+ +
- -
= = =∫ ∫ (5.3.9)
但要注意,上述的量化噪声功率并不是系统最终输出的量化噪声功率.这是因为,由图5.3.8
看出,( )qe t的最小周期大致是抽样频率sf的倒数,而且大于1/sf的任意周期都可能出现.因
此,从频谱的角度看,( )qe t的频谱将从很低频开始一直延伸到频率sf甚至更高.设( )qe t的功
率谱密度为( )eP f),则可以近似认为,
2
( )
3e
s
P f
f
σ
=, 0sf f< <
这就是说,( )qe t的平均功率被认为均匀地分布在频率范围(0,sf)之内.这样,具有功率谱密
度为( )eP f的噪声,通过低通滤波器(截止频率为mf)之后的量化噪声功率为:
2
( )
2
m
q e m
s
f
N P f f
f
σ
= =
(5.3.10)
由此可见, M系统输出的量化噪声功率与量化台阶σ及比值( / )m sf f有关,而与输入信号号
的幅度无关.当然,这后一条性质是在未过载的前提下才成立的.不发生过载现象,这实际
上是对输入信号的一个限制.
在临界条件最大的信噪比为:
2 2 3 3
0
2 2 2 2 2
3
0.04
8 8
s s s
q k k m k m
S f f f
N f f f f f
σ
π π
= = ≈i (5.3.11)
56
由此可见,最大信噪比0( / )qS S与抽样频率sf的三次方成正比.而与信号频率kf的二次
方成反比.因此,对于 M系统而言,提高抽样频率将能明显地提高信号与量化噪声的功率
比.
5.3.3 增量增量增量增量(差分差分差分差分)脉冲编码脉冲编码脉冲编码脉冲编码调制调制调制调制(DPCM)系统系统系统系统
前面对PCM和 M系统性能的分析可以看出,在不考虑信道误码率的情况下, M的
性能通常比PCM差.这主要是因为 M系统不管误差信号如何变化,传输的增量σ是固定
不变的.如果使增量的数值随误差信号( )tε的变化量化成M个电平之一,然后再进行编码,
这样,系统的性能就会得到改善.在这样的系统中,由于对传输的增量还要经过脉冲编码调
制,因面称它为增量脉冲编码调制或差分脉冲编码调制(DPCM).
图5.3.10示出了一增量脉冲编码调制系统的组成方框图.图中( )tε是输入信号( )m t和量
化信号'( )m t相减后的误差信号.该误差信号经抽样,量化与编码后,得到一二进制脉冲序列,
即DPCM信号.它一方面送往信道传输,另一方面经反馈支路的译码器和积分器,以便提供
所需的量化信号'( )m t.
图图图图5.3.10 增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图
差分脉冲编码调制信号在接收端首先经译码器恢复成量化的误差信号,再经积分器和低
通滤波器便可恢复出原发送信号.
下面我们以两位编码(即2N=)来说明差分脉冲编码调制的编码过程.由于2N=,故误
差信号的量化电平数4M=(因2NM=).假设4个量化电平分别为3 , , , 3v v v v+ + - - ,它
们分别由二进制脉冲(双极性), , ,++ +- -+ - -表示.误差信号( )tε的抽样,量化和编码过程如图
57
5.3.11所示.图(a)表示对误差信号的抽样与量化,符号"0"表示( )tε的实际抽样值,"."表
示量化后的值.由图可见,当误差信号的取值在0 ~ 2v 范围时,则被量化为v+ ,即当
0 ( ) 2t vε≤ < 时,量化为v+ .同理,当2 ( ) 0tε- ≤ <时,量化为v- .2 ( )v tε ≤时,量化
为3v+ .( ) 2t vε=0 & abs(x(i))<32
67
out(i,2)=0; out(i,3)=0; out(i,4)=0; step=2; st=0;
elseif 32<=abs(x(i)) & abs(x(i))<64
out(i,2)=0; out(i,3)=0; out(i,4)=1; step=2; st=32;
elseif 64<=abs(x(i)) & abs(x(i))<128
out(i,2)=0; out(i,3)=1; out(i,4)=0; step=4; st=64;
elseif 128<=abs(x(i)) & abs(x(i))<256
out(i,2)=0; out(i,3)=1; out(i,4)=1; step=8; st=128;
elseif 256<=abs(x(i)) & abs(x(i))<512
out(i,2)=1; out(i,3)=0; out(i,4)=0; step=16; st=256;
elseif 512<=abs(x(i)) & abs(x(i))<1024
out(i,2)=1; out(i,3)=0; out(i,4)=1; step=32; st=512;
elseif 1024<=abs(x(i)) & abs(x(i))<2048
out(i,2)=1; out(i,3)=1; out(i,4)=0; step=64; st=1024;
elseif 2048<=abs(x(i)) & abs(x(i))=4096)
out(i,2:8)=[1 1 1 1 1 1 1];
else
tmp=floor((abs(x(i))-st)/step);
t=dec2bin(tmp,4)-48; % 函数dec2bin输出的是ASCII字符串,48对应0
out(i,5:8)=t(1:4);
end
end
out=reshape(out',1,8*n);
***********************************
其中解码函数为pcm_decode.m.
*******************************
68
function [out]=pcm_decode(in,v)
n=length(in);
in=reshape(in',8,n/8)';
slot(1)=0; slot(2)=32; slot(3)=64; slot(4)=128;
slot(5)=256; slot(6)=512; slot(7)=1024; slot(8)=2048;
step(1)=2; step(2)=2; step(3)=4; step(4)=8;
step(5)=16; step(6)=32; step(7)=64; step(8)=128;
for i=1:n/8
ss=2*in(i,1)-1;
tmp=in(i,2)*4+in(i,3)*2+in(i,4)+1;
st=slot(tmp);
dt=(in(i,5)*8+in(i,6)*4+in(i,7)*2+in(i,8))*step(tmp)+0.5*step(tmp);
out(i)=ss*(st+dt)/4096*v;
end
*******************************
运行结果如图5.3.16和图5.3.17所示.
012345678910
-1
-0.5
0
0.5
1
sample sequence
012345678910
-1
-0.5
0
0.5
1
pcm decode sequence
图图图图5.3.16 A律律律律PCM编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图
69
-60-50-40-30-20-100
0
10
20
30
40
50
60
输输窄输原原输输由dB值
窄量原量量
图图图图5.3.17 A律律律律PCM量化信噪比的动态范围量化信噪比的动态范围量化信噪比的动态范围量化信噪比的动态范围
利用上述程序,将时间t的间隔变为1/ 4096t =,则A律PCM编码后正弦信号与输入正
弦信号的波形对比如图5.3.18所示.
00.10.20.30.40.50.60.70.8
-1
-0.5
0
0.5
1
sample sequence
00.10.20.30.40.50.60.70.8
-1
-0.5
0
0.5
1
pcm decode sequence
图图图图5.3.18 A律律律律PCM编码后波形与输入波形的对编码后波形与输入波形的对编码后波形与输入波形的对编码后波形与输入波形的对比比比比示意图示意图示意图示意图
70
5.3.4 MATLAB提供的信源编提供的信源编提供的信源编提供的信源编/译码函数译码函数译码函数译码函数
在MATLAB通信工具箱中,提供了两种信源编/译码的方法,即标量量化和预测量化.
下面介始一下其使用方法.
1. 标量量化标量量化标量量化标量量化
标量量化就是给每一个落入某一特定范围的输入信号分配一个单独值的过程,并且落入
不同范围内的信号分配的值也各不相同.
(1) 信源编码中的μ律或A压扩计算函数compand.
格式一:out=compand(in,param,V)
格式二:out=compand(in,param,V,method)
功能:格式一实现μ律压扩;格式二实现μ律或A压扩.
参数含义:param对格式一为μ律;对格式二为μ或A值.V为峰值.压扩方法由methode
参指定,有:'mu/compressor', μ律压缩;'mu/expander',μ律扩展;'A/compressor',A律
压缩;'A/expander', A律扩展.
(2) 产生量化索引和量化输出函数quantiz
格式一:indx=quantize(sig, partition)
功能:根据判断向量partition,对输入信号sig产生量化索引indx,indx的长度与sig矢
量的长度相同.向量partition则是由若干个边界判断点且各边界点的大小严格按升序排列组
成的实矢量.若partition的矢量长度小于1N-,则索引向量indx中的每个元素的大小为
[0, 1]N-范围内的一个整数.量化方法如下:
若信号sig小于或等于partition(1),则输出0;若信号sig大于partition(1),而小于等于
partition(i+1),则输出1;若信号sig大于partition(N-1),则输出N-1.
格式二:[indx, quant]=quantiz(sig, partition, codebook)
功能:根据码本codebook,产生量化索引indx和信号的量化值quant.codebook存放indx
每个partition的量化值,对应indx=i-1的值在codebook(i).若partition的矢量长度为1N-,
则codebook长度为N.
格式三:[indx, quant, distor]=quantiz(sig, partition, codebook)
功能:产生量化索引indx,信号的量化值quant及量化误差disor.
(3) 采用训练序列和Lloyd算法优化标量算法的函数lloyds
71
格式:[partition, codebook]=Lloyds(training_set, ini_codebook)
功能:用训练集矢量training_set优化标量量化参八partition和码本codebook.Ini_codebook
是码本codebook的初始值.码本长度大于等于2,输出码本的长度与初始码本长度相同.输
出量化比参数partition的长度较码本长度小1.当ini_codebook为整数时,该函数以其作为码
本长茺.当处理后相对误差小于710-时,停止处理.
程序示列程序示列程序示列程序示列::::
用训练序列Lloyd算法,对一个正弦信号数据进行标量量化.
程序prog5301.m如下.
**********************************************
clear all
clear close
N=2^3;
t=[0:100]*pi/20; % 以3比特传输信道
u=cos(t);
[p,c]=lloyds(u,N); % 生成分界点矢量和编码手册
[index,quant,distor]=quantiz(u,p,c); % 量化信号
plot(t,u,t,quant,'*')
**********************************************
运行上述程序将显示出正弦信号的量化效果图.如图5.3.19所示.
0246810121416
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
图图图图5.3.19 标量量化前标量量化前标量量化前标量量化前,,,,后信号比较后信号比较后信号比较后信号比较
72
2. 预测量化预测量化预测量化预测量化
如果知道一些发送信号的先验信息,就可以利用这些信息,根据过去发送信号来估计下
一个将要发送的信号值,这样的过程即称之为预测量化.
差分脉冲编码调制(DPCM)就是采取了预测量化技术.MATLAB通信工具箱提供了
dpcmenco( ),dpcmdeco( )和decmopt( )函数.
(1) 差分脉冲调制编码函数decmenco
格式一:indx=dpemenco(sig, codebook, partition, predictor)
功能:返回DPCM编码的编码索引indx.其中参数sig输入信号,predictor为预测器传
递函数,其形成为1[0, , , ]mt t .预测误差的量化参数由partition和predictor指定.
格式二:[indx,quant]=dpemenco(sig, codebook, partition, predictor)
功能:除产生DPCM编码的编码索引indx外,还产生量化值quant.输入参数codebook,
partition, predictor可以由函数decmopt函数估计.当预测器为一阶传递函数时,为DPCM码
调制.
(2) 信源编码中的DPCM解码函数decmdeco
格式一:sig=dpemdeco(indx, codebook, predictor)
功能:根据DPCM信号索引indx进行解码.predictor为指定的预测器,codebook为码本.
格式二:[sig,quant]=dpemdeco(indx, codebook, predictor)
功能:根据DPCM信号索引indx进行解码,同时输出量化的预测误差quant.参数
codebook,predictor可以用dpcmoot函数估计.通常m阶预测器传递函数的形式为1[0, , , ]mt t .
(3) 用训练数据优化差分脉冲调制参数的函数dpcmopt
格式一:predictor=dpcmopt(training_set, ord)
功能:对给定训练集的预测器进行估计,训练集及其顺序由training_set和ord指定,预
测器由predictor输出.
格式一:[predictor, codebook, partition]=dpcmopt(training_set, ord, ini_codebook)
功能:输出预测器predictor,优选码本codebook,预测误差partition.输入变量ini_codebook
可以是码本矢量的初值或其长度.
程序示例程序示例程序示例程序示例::::
用训练数据优化DPCM方法,对一个余弦信号数据进行标量量化.
MATLAB程序prog5302.m如下.
73
***********************************************
clear all
clear close
N=2^3;
t=[0:100]*pi/20; % 以3比特传输信道
u=cos(t);
[predictor,codebook,partition]=dpcmopt(u,1,N); % 优化的预测传递函数
[index,quant]=dpcmenco(u,codebook,partition,predictor); % 使用DPCM编码
[sig,quant]=dpcmdeco(index,codebook,predictor); % 使用DPCM解码
plot(t,u,t,quant,'*')
***********************************************
程序运行结果如图5.3.20所示
0246810121416
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
图图图图5.3.20 DPCM预测量化误差预测量化误差预测量化误差预测量化误差
5.4 语音的线性预测编语音的线性预测编语音的线性预测编语音的线性预测编码码码码
线性预测分析方法让最有效的语音分析技术之一,其特点是既能极为精确地估计语音参
数,又有比较快的计算速度.因此,应用线性预测方法对语音编码已属高效率编码技术,它
不同于常规的数字编码.在上面讨论的各种量化编码方法,诸如PCM(脉码调制),M (增量
74
调制),DPCM(差值脉码调制)等都是常规的编码方法.因为这些方法部都是通用的波形数字
编码方法,并未考虑语音本身的特点,所以数据压缩率不可能很高.
线性预测的理论与方法在语音信号处理中的应用,主要表现在声码器技术领域中.除了
对语音信号的线性模型参数进行估计外,还能对所有的语音基本参数,诸如基音周期,共振
峰,语音谱,声强等进行估计.本节的日的不是讨论语音信号处理的所有问题,而是对语音
的线性预测编码(常用缩写字母LPC表示),主要是对语音模型参数进行估计.因此,首先要
介绍语音的数字模型;然后应用上节讨论的线性预测参数估计的方法来讨论语音模型的参数
估计问题,最后介绍一下基本语音参数的线性预测估计及相应的声码器原理.
5.4.1 语音信号的数学模型语音信号的数学模型语音信号的数学模型语音信号的数学模型
语音的生成机构大致可分为三部分.第一是声源声源声源声源.大家都知道人的声带是一种声源,它
会发生自激振动,产生常称为元音的语音,有时就叫作声音.另外还有两种声源:由声道变
窄(通常由舌头来控制)形成湍流叹声的称为摩擦声源.它产生的声音就叫摩擦音;由闭合的
声道急开而形成的脉冲波产生湍流噪声的叫爆破声源.第二是共鸣机构共鸣机构共鸣机构共鸣机构.它由鼻腔,口腔与
舌头组成,有时也叫声道声道声道声道.第三是放射机构放射机构放射机构放射机构.它由嘴唇或鼻孔发出声音并向空间传播出去.
对于这样的一种人类发声机能,可用多种的模型来模拟.我们这里特别关心的是所谓数字模数字模数字模数字模
型型型型,即能用数字处理方法来实现上述发声机能的一种模型.这种模型用准周期的脉冲源模拟
声带的振动,用随机嗓声(一种采用伪随机序列)模拟摩擦声,用可变参数的数学滤波器来模
拟声道谐振特性与放射持性.发声器官的数字模型如图5.4.1所示,图中开关是作浊音与清音
转换控制用的,乘法器是完成音强(或叫增益)控制功能的.数字滤波器的频率特性随受控参
数的变化而改变.如果所有的控制信号都由真实的语言信导分析所得,那么该滤波器(也叫声
道滤波器)的输出便接近于原始语言信号的信号序列,可以恢复出语音.这也就是声码器的基
本原理.
图图图图5.4.1 发声器官的数字模型发声器官的数字模型发声器官的数字模型发声器官的数字模型
75
现在用n来表示离散时间变量,( )s n表示n时刻语音信号的样值,并假设语音信号样值
序列是符号AR模型的,则有:
1
( ) ( ) ( )
p
i
i
s n a n i w n
=
= - +∑ (5.4.1)
其中,iα为模型参数,( )w n为白噪声序列.这样,语音的线性预测编码与AR模型参数估计
的关系可陈述如下:
若对输入的语音样值序列作出模型参数与残差估计,则n时刻的样值估计为:
1
^^( ) ( )
p
i
i
s n s n iα
=
= -∑ (5.4.2)
其中,^{ }iα为模型参数估计值的集合.估计误差或称残差信号便为:
1
^^( ) ( ) ( ) ( ) ( )
p
i
i
e n s n s n s n s n iα
=
= - = - -∑ (5.4.3)
若将( )e n作为图5.4.1声道滤波器的输入,滤波器能输出语言信号( )s n的近似值,则滤波器的
传递函数必须与n时刻的的输入,输出信号相适配.为此,对式(5.4.3)两边作Z变换,得:
1
^( ) ( ) ( )
p
i
i
i
E Z S Z Z S Zα-
=
= -∑ (5.4.4)
于是传递函为:
1
1
( )
^( ) 1
( )
p
i
i
i
S Z
H Z Z
E Z
α
-
-
=
= = -
∑ (5.4.5)
可知,声道滤波器是—个全极点型数字滤波器.该滤波器的参数{ | 1, 2, , }ii pα= 由语音分析
所得到的参数估计集合^{ }iα确定.于是对语音信导的编码转变为对参数集合^{ }iα及残差e的编
码问题.在接收端可由图5.4.1所示的数学模型恢复出语音信号.因此,线性预测编码的基本
问题就是很据输入的语音信号样值作出相应的AR模型参数估计.
5.4.2 由短时相关函数估由短时相关函数估由短时相关函数估由短时相关函数估计模型参数计模型参数计模型参数计模型参数
在5.2节里已作过介绍的AR模型参数的估计l方法,当然可以应用于语音信号A R模型
参数估计.不论是最小二乘法还是矩估计方法,估计的主要依据都是信息样值的自协方函数
76
或自相关函数.从自协方差函数或自相关函数的定义来看,应该从语音样值的整个集合来计
算概率平均值(参看式(5.2.3)与(5.2.4)).实际上这是做不到的.我们只能在时间轴上作短时平
均计算,以求得相应的短时间内发出的某个音的最佳预测参效,也就是动态最佳参数.如果
求整个长时间的平均,那么,对长时间的预测误差来说,系统可达到动态最佳状态.但对某
一短时间段信号来说,未必是最好的.也就是说,没有适应动态变化的系统,就达不到传输
动态信息的目的.在语音信号处理中,短时平均值的计算,通常采用两种方法,即协方差法
与自相关法.本书主要讨论两种方法.
1. 协协协协方差法方差法方差法方差法
在整个语音取样序列中,从n时刻开始取一段N个样值,记作( ), 0,1, , 1ns m m N= - .由
于对i时刻取样的样值常记为( )s i,因此该N个样值又可记作( 0), ( 1), , ( 1)s n s n s n N+ + + - .
两种记法之间的关系为:
( ) ( )ns m s n m= + 0,1, , 1m N= - (5.4.6)
于是,最二乘法方程式(5.2.10)改写为:
11
0 1 0
( ) ( ) ( 1) ( )
pN N
n n n n j
m j m
s m s m k s m s m kα
--
= = =
- = - -
∑ ∑∑ 1, 2, ,k p= (5.4.7)
它们的对应关系为:式(5.2.10)中的, ,i n t分别相当于式(5.4.7)中的, ,m N n;tix相当于( )ns m;
,t k ix-相当于( ) ( ) ( )n k ns m s n k m s m k-= - + = -.若令,
1
2
0
( ) ( , )
N
n n
m
s m k k k
-
=
- =∑ 1k p≤ < (5.4.8)
1
0
( ) ( ) ( , )
N
n n n
m
s m j s m k k j
-
=
- - =∑ 1k p≤ <, 0j p≤ ≤ (5.4.9)
则式(5.4.7)又可改写成:
1
2
(1,1) (1, 2) (1, ) (1,0)
(2,1) (2, 2) (2, ) (2,0)
( ,1) ( , 2) ( , ) ( ,0)
n n n n
n n n n
pn n n n
p
p
p p p p p
α
α
α
=
(5.4.10)
或简化为:
1
( , ) ( ,0)
p
j n n
j
k j kα
=
=∑ 1, 2, ,k p= (5.4.11)
77
式中的( , )nk j 就是信号的协方差函数.如果求出这些协方差函数值,那么最佳线性参数的估
值问题就是解方程组(5.4.11)的问题了.它们可由递推解法解出.协方差的计算是由计算机程
序来完成的.为了正确地编制程序,需对式(5.4.8)与(5.4.9)有确切的理解.为此举例说明如下:
设输入信号取样序列为(0), (1), (2),s s s ,计算5, 8, 4n N p= = =时( , )nk j 的值.下面我
们来写出式(5.4.10)左边矩阵的第一行与右边列矢量的元素,它们分别为:
以及
图5.4.2给出计算这8个协方差的各样值间的关系.对其他, ,n N p值的( , )nk j 的数值也
可作类似的计算.
图图图图5.4.2 求协方差的样值差求协方差的样值差求协方差的样值差求协方差的样值差,,,,5, 8, 4n N p= = =
符号"^"表示该两时刻样值相乘;"
"表示所有值相加.
78
2. 自相关法自相关法自相关法自相关法
另一种求短时平均的方法是用窗函数( )w m截取信号序列( )( 0,1, 2, )ns m m= 中N长的一
段,并在此范围内求平均,即:
( ) ( ) ( )ns m s n m w m= + (5.4.12)
其中,若:
1 0 1
( )
0
m N
w m
≤ ≤ -
=
其他
.则意味着除( ), ( 1), , ( 1)s n s n s n N+ + - 这N个样值外,
其余样值均为0.于是协方差函数表达式(5.4.8)与(5.4.9)的具体形式也有所不同.它们的主要
变化在于求和限要改变,比如当用0m=时.式(5.4.9)中将出现( )s n k-与( )s n j-,显然它们
均为0.为了清楚起见,作变量代换,令l m k= -,并考虑到( 1)s n N+ -右侧信号为0.故要
求1l k j N+ - ≤ -才行.于是式(5.4.9)变成:
1 ( )
0
( ) ( ) ( , )
N k j
n n n
l
s l k j s l k j
- - -
=
+ - =∑ 1, 2, , ; 0,1, ,k p j p= = (5.4.14)
式(5.4.14)与式(5.2.25)自协方差函数的表达式有相同的涵义,只是系数1/N及符号定义不同.
自然也可以求出自相关函数( )kρ的表达式.由于被窗口函数截取的样本在求平均过程中不改
变,故有( , ) ( , )n nk j j k =的特点.在语音信号处理中,通用的记号为:
( , ) ( , ) ( ) ( )n n n nk j j k R k j R τ= = - = (5.4.15)
式中:k jτ= -.这样,式(5.4.14)可改写成:
1
0
( ) ( ) ( )
N
n n n
l
R s l s l
τ
τ τ
- -
=
= +∑ (5.4.16)
与此相对应的AR模型参数估计方法是矩估计方法,要求解的方程组便是Yule-Walker方程.
应用式(5.4.16)的记号,方程式(5.2.23)变为:
1
(| |) ( )
p
k n n
j
R k j R kα
=
- =∑ 1, 2, ,k p= (5.4.17)
展开矩阵形式为:
1
2
(0) (1) ( 1) (1)
(1) (0) ( 2) (2)
( 1) ( 2) (0) ( )
n n n n
n n n n
pn n n n
R R R p R
R R R p R
R p R p R R p
α
α
α
-
- =
- -
(5.4.18)
79
左边的R矩阵为Toepliz矩阵,故式(5.4.18)的求解可以应用Durbin迭代程序或格型滤波算法.
如前所述的.
为了使在计算R矩阵中对各样本的应用关系更清晰,仍以8, 5, 4N n p= = =为例为说明.
由式(5.4.16)可写出,
图5.4.3是这个计算公式的示意图,很明显,它比图5.2.2简单,计算量也少得多,但它
的估计精度度不如协方差方法.还应该指出的是,由式(5.4.13)表示的窗函数,由于边缘的突
变会引起边缘样本系数估计的误差增大,有时候把它叫"边缘效应".为避免这个现象,往往
需对样本作预加重处理和采用和采用边缘缓变的窗函数.在线性预测声码器(LPC)中,需要考
虑上述的技术措施来改善合成语音的质量.这里不作详细讨论.
图图图图5.4.3 求自相关函数的样值间关系求自相关函数的样值间关系求自相关函数的样值间关系求自相关函数的样值间关系,,,,8, 5, 4N n p= = =
80
5.4.3 基音周期检测基音周期检测基音周期检测基音周期检测
利用线性预测方法提取话音参数而组成的声码器称为线性预测声码器.由于话音信号为
非平稳的随机过程,因此用时相关系数来过行语音参数的估计.除了模型参数的估计外,由
图5.4.3的模型可知,还需估计出基音周期及激励信号.此外,为了传输连续变化的语音信号,
先要将语音信号分帧,一般以10~30 ms数据作为一帧.在接收端两逐帧地进行合成并连接起
来组成连续话音输山.图5.4.4是线性预测声码器的原理框图.利用这样的声码器来传输话音
信号便可达到压缩数据率的目的.
图图图图5.4.4 线性预测声码器原理框图线性预测声码器原理框图线性预测声码器原理框图线性预测声码器原理框图....(a) 语音分析器语音分析器语音分析器语音分析器;;;;(b) 语音合成器语音合成器语音合成器语音合成器
语音分析器中参数估计方法在前几节中已作了详细的讨论,现在需要讨论基音周期检测
的方法.所谓基音周期检测或提取基音,就是要确定信号周期的长度.从理论上说这可以有
许多方法.但是由于语音信号的随机性,音调变化范围很宽,且语音还有嗓音或称浊音(Voiced)
与非嗓音或称清音(Unvoiced)之分,前者含有基音而后者却无基音,这就使基音提取变得复杂
且易出错.若基音提取有错误,会直接导致合成语音质量的下降,因此提出了各种改进方法.
这里只介绍1972年由Markel提出的称作简单逆滤波器跟踪(SIFT)算法,这是基于LPC分析
的一种方法,可以作为线性预测方法应用的一个例子.
图5.4.5是SIFT算法的原理框图.由于基音周期信息是低频信息,因此输人信号( )s n
81
首先经过截止频卒为900 Hz的低通滤波器.经滤波后的信号频带变为0~1kHz,故可以将原
来10kHz的取样频率经5次分频,降列2kHz(每5个样值去掉4个),得到离散信号( )x n.然
后从信号( )x n提取基音周期,其主要过程为逆滤波与求自相关函数峰值位置.
图图图图5.4.5 基音检测的基音检测的基音检测的基音检测的SIFT算法框图算法框图算法框图算法框图
所谓逆滤波,是对声道滤波的逆过程,即以经过低通滤波的语声信号( )x n作为输入,式
(5.4.6)的分母作为传递函数的一种滤波处理.该传递函数可表示为:
1
( ) 1
p
i
i
i
A Z Zα-
=
= -∑ (5.4.19)
其中滤波器参数iα由线性预测,即由式(5.4.10)或式(5.4.18)中的iα来估计.
信号( )x n经过逆滤波得到( )y n,实质上它就是预测器的预测误差.由线性预测理论可知,
( )y n有接近白噪声的属性,它的谱是接近平坦的,各样值间基本上互不相关.这就表明由式
(4.3.15)与(5.4.16)所定义的自相关函数一般很小,只有当0τ=或信号具有周期性(以P为周期,
, 2 ,P P P= ± ± )时,其自相关函数达到最大值.也就是说,可不用考虑信号的起始时间,只
要确定自相关函数中第一个最大值的τ值就可以代表其周期估计.图5.4.6表示这个分析过程
的几个典型波形.其中图(a)为一段被分析的输入波形,(b)为逆滤波器输出信号的时域波形;
(c)为逆滤波器输出信号的归一化自相关函数.由图可见,当8τ= ms时有明显的峰值,它便
可大致地代表该信号段的基音周期.
为了获得进行周期值的更高分辨率,可以对最大峰值所处的一定范围内的自相关函数进
行内插.当自相关函数的归一化峰值电平低于给定门限值时,则判定该论语音段为清音.
82
图图图图5.4.6 SIFT算法的典型信号波形算法的典型信号波形算法的典型信号波形算法的典型信号波形
SIFT算法利用线性预测分析提供了一个平坦谱的信号( )y n,以便于进行基音检测.谱的
平坦化越好(预测很确),对基音检测效果也起好,然而,对于高基音频率的声源(如儿童语音)
的谱不易平坦化.这是由于在0~900 Hz范围内缺乏一个以上的基音谐波(对由电话线输入的
信号尤其如此).对这类说话者及传输条件.就要采用别的基音检测方法了.
5.4.4 共振峰检测共振峰检测共振峰检测共振峰检测
在语音信号的数字模型一节中所叙述的语音生成机构的第二部分 —— 声道.一般可简
化成一个非均匀截面的声道.当声音从声源出发顺着声管传播时,它的频谱形状会被声管的
选择性所改变.这一效应正如同管乐器中所发生的谐振现象一样.声道管的谐振频率被称作
共振峰频率,共振峰频率与声道的形状和大小有关,每种形状都有一套共振峰频率作为其特
征,依次叫作第一共振峰频中,第二共振峰频率,……,并分别记为1 2, ,F F .人们借助
于舌,口腔等器官来改变声道形状,从而改变共振峰频率特征,使声音变成了可辨别的语言.
无疑,除了基音周期之外,共振峰频率是语音信号的另一个重要参数.图5.4.7显示了两个英
语元音的声波形图,除了有明显的基音周期外,其波形的结构有较大的差别.其中图(a)为元
音[i: ]的声波形,在一个周期内可看出有一个低频的衰减振荡,其上还叠加一个较强的高频率
83
振荡.可见该元音的第一共振峰频率较低,第二,第三共振峰频率较高.(b)为元音[u: ]的声
波形,可看出它的第一,第二共振峰频率很低,所以它的高频能量很小.
图图图图5.4.7 元音元音元音元音[i: ]与与与与[u: ]的声波形图的声波形图的声波形图的声波形图
以上讨论可知,若以共振峰参数来代表语音信息,同时配以基音周期参数,也可以很好
地合成语音信号.这样的语音分析与合成设备被称为共振峰声码器,图5.4.8是它的原理框图.
比较图5.4.4与图5.4.8可以发现,除了用共振峰分析与合成器代替预测参数分析与控制器之
外,其他都是类似的.我们的讨论侧重于如何应用线性预测方法来确定共振峰参数.确定共
振峰的方法很多,且各有优缺点.如求根法,即求出由式(5.4.19)所表示的传递函数的多项式
的根,也就是由式(5.4.5)所表示的声道滤波器的极点.
图图图图5.4.8 共振峰声码器原理框图共振峰声码器原理框图共振峰声码器原理框图共振峰声码器原理框图
84
对于大多数语来说,全极点模型很好地代表了声道的响应.但声学理论指出,鼻音和摩
擦音还要有谐振与反谐振(极点零点),这给数学处理带来了麻烦.但数学上很易证明当| | 1a<
时,有:
1
1
0
1n n
n
aZ a Z
-∞
- -
=
- =
∑ (5.4.20)
因而可以用多个极点去逼近一个零点,这就可以用全极点模型来讨论声道的谐振动特性.因
为式(5.4.4)分母的系数都是实数,所以分母多项式的根将是实根或复共轭根.将LPC数估计iα
代入式(5.4.5)的分母,得到多项式( )A Z为:
1
( ) 1
p
k
k
k
A Z Zα-
=
= -∑ (5.4.21)
若{ | 1, 2, , }iZ i p= 是( )A Z的一组根,则式(5.4.21)又可表示成:
1
1
( ) (1 )
p
k
k
A Z Z Z-
=
= -∏ (5.4.22)
由于Z平面上点与S平面上点的对应关系为:
kS T
kZ e= (5.4.23)
其中2k k kS j Fσ π= - +是与Z平面中的iZ相对应的S平面上的根,一般可表示成示成复共轭对
形式,即:
*, 2k k k kS S j Fσ π= - ± (5.4.24)
其中,"*"表示复共轭.kS和*
kS为声道的复谐振频率的一般表达式,各个中心频率为2kFπ,
它们可以直接由式(5.4.21)的根求出.共振峰的求根法的基本原理也基于此.
将式(5.4.24)代入式(5.4.23)并展开,得:
( 2 )*, cos(2 ) sin(2 )
Re Im
k k kj F T T
k kk k
k k
Z Z e e F T j F T
Z j Z
σ π σπ π- ± -= = ±
= ±
(5.4.25)
因此,如果已经解出LPC参数1 2, , ,pα α α ,那么通过解多项式(5.4.21)的根,便可以找出Z平
面上的p个根{ | 1, 2, , }iZ i p= .于是可求现相应的共振峰频率为:
85
Im1
arctan
2 Re
k
k
k
Z
F
T Zπ
=
(5.4.26)
我们定义频率响应大于中心频率的峰值响应的0.707倍的频率间隔为谐振带宽谐振带宽谐振带宽谐振带宽,则存中
心频率为2kFπ时,声道的谐振带宽近似为2kσ.在Z平面中.从原点到极点的距离就决定了
共振峰的带宽,即
| |kT
kZ eσ-= (5.4.27)
如果以Hz为单位计算,则带宽为:
21
ln | |
2
k
k kB Z
T
σ
π π
= = (Hz) (5.4.28)
由此可见,用线性预测方法作共振峰分析,经过对多项式求根,可较精确地确定共振峰
的中心频率与带宽,也可较容易地找出极点与共振峰的对应关系.这就为构成低比特率的共
振峰声码器提供了条件. '
5.5 图像的预测图像的预测图像的预测图像的预测编码编码编码编码
在图像压缩编码方法中,预测编码仍是主要手段之一,它常常与差值脉码调制(DPCM)
同义.在图像信源中.除了静态图像外,很大一部分为活动图像,其中电视网像又是最为常
见的一种活动图像.活动图像属于三维数据的信号源了.也就是说,活动图像是在时间轴上
一幅一随地更换出现的.比如电视图像,为了应用人眼的视觉残留效应而使之产生连续感,
则采用每秒25帧的更换速率,这里用'帧"代替通常叫"幅"的术语,由于帧速很高,这就
使各帧图象之间相差甚微,或者说存在着帧间相关性,于是,出现了帧内预测编码与帧间预
测编码或两者联合使用的编码方式.从而提高了图形数据的压缩倍数.
由于线性预测的基本理论及语言信号的线性预测编码都已作详细讨论,本节只介绍图像
预测编码的一些特殊点.首先介绍图像的数字化表示,然后介绍帧内DPCM乃帧间DPCM的
基本思想.
5.5.1 图像的数学化图像的数学化图像的数学化图像的数学化
所谓图像,就是一种视觉信息.它是用一个所谓二维光强度函数( , )f x y来表示的.这
里,x y表示空间坐标,而该点的函数值( , )f x y则表示试点的"光强度光强度光强度光强度",或通常称之为亮度或
86
灰度.因为光是能量的一种形式,故( , )f x y必为非零值,并且是在有限值城内的连续函数,
即:
0 ( , )f x y< < ∞ (5.5.1)
数字化图像则是存二维守间中对)t强度面故/(x,y)取样与量化的结果,以八x Y坍
短阵表示为:
(0,0) (0,1) (0, 1)
(1,0) (1,1) (1, 1)
( , )
( 1,0) ( 1,1) ( 1, 1)
f f f N
f f f N
f x y
f N f N f N N
-
- ≈
- - - -
(5.5.2)
式中左边( , )f x y是连续量,而右边矩阵使是通常所谓的数字图像数字图像数字图像数字图像.矩阵的每一个元素称作像
元或像素,常用缩写字母"Pel (Picture E1ements)"来表示.
现在,我们需要关心的两个参数是矩阵的大小N与每个像素的量化层数G.为了数字处
理的方便,一般都使它们由2的整数幂来表示,即有:
2nN= (5.5.3)
及
2mG= (5.5.3)
式中,,n m均为整数.m实际上就是每个像素的量化比特数.因此,一帧N N×的数字图像
的总比特数便为:
B N N m= × × (5.5.5)
例如,灰度为64级的128×128数字图像就有98304 bit数据量,需要12288 B = 12 KB的存储
单元.
那么,N与m应该多大才能使式(5.5.2)中的近似程度符合晴晰度要求呢 显然N与m越
大,清晰度就越好.但从处理,存储与传输的角度考虑,不希望N与m太大,以致于数据量
不会太大.根据国际标准CCIR 601建议的数字电视标准数据,高度信号的取样频率为13.5
MHz,两个色度信号取样频率各为6.75MHz,8m= bit,于是,数字电视的比特率高达
216Mb/s,高清晰度电视的比特率更达1.188Gb/s.处理一般图像所用的最小系统,也应能显
示256×256个像素与64灰皮级,即N=256与m=6才行.即使如此,一帧图像的数据量仍
然十分大,故数字图像的数据压缩要求比语声信号的更为迫切.
87
5.5.2 帧内预测编码帧内预测编码帧内预测编码帧内预测编码
帧内预测编码利用的是一帧图像中各像素之间的相关性.正如式(5.5.2)所表示的数字图
像那样,它由一行行像素组成,一行中各像素之间又存在相关性,行与行之间也存在相关性.
如果只利用一行内像素间的相关性进行预测,则叫一维预测编码.如果同时利用行内与行间
相关性进行预测,则为二维预测编码.
1. 一维预测一维预测一维预测一维预测
一维预测的基本原理与线性预测的基本原理是相同的.但在图像预测编码中,往往不去
计算相关系数.而是由实验观察确定.这不仅是为计算方便,还因为按均方误差准则作出的
最佳预测,并不反映人眼的主观感觉.一维预测在传真图像,即二维图像编码中应用较多.
图图图图5.5.1 一维预测过程一维预测过程一维预测过程一维预测过程....(a) 扫描行扫描行扫描行扫描行;;;;(b) 预测器预测器预测器预测器;;;;(c) 逆预测器逆预测器逆预测器逆预测器....
给扫描行上的像素编上如图5.5.1(a)所示的号码.预测算法就是用已经扫描过的像素序列
1 2, ,i ix x- - 作为参考像素,按某种数学关系1 2( , , )i if x x- - 对现时像素ix的值作出预测,得到预
测值^
ix.然后将^
ix与ix作比较,若一致记作"0",不一致则记作"1".在二值图像酌情况下,
可以用模2加来完成这个操作.在多灰度值图像时,失得到ix与^
ix的差值,然后对差值作出
88
量化编码,这就是已介绍过的DPCM过程.这个关系的数字表达式为:
1 2
^( , , )
^( )
^( )
i i i
i i i
i i i
x f x x
y x x
y x x
- -=
= ⊕
= -
二值图像
多值图像
(5.5.6)
符号⊕表示模2加.
上面已指出,( )fi表达式由实验观察确定,一般仍为线性多项式,即,
1
^
l
i k i k
k
x a x-
=
=∑ (5.5.7)
其中系数ia就是根据实验确定的参数.若1l=,实际上就是前一像素预测.在二值传真图像
情况,由于扫描结果总是连续"0"与连续"1"相间,因此,误差图像代表了图像轮廓,如图5.5.2
所示.
图图图图5.5.2 前像素预测的误差图样前像素预测的误差图样前像素预测的误差图样前像素预测的误差图样....(a) 两相邻扫描行两相邻扫描行两相邻扫描行两相邻扫描行;;;;(b) 前像素预测误差图样前像素预测误差图样前像素预测误差图样前像素预测误差图样
2. 二维预测二维预测二维预测二维预测
虽然图像编码的二维DPCM可以由一些经验给出预测系数.但严格地说,图像编码的二
维DPCM应根据不同的图像所求得的自相关函数来求出最佳预测系数,然后设计出与图像相
匹配的预测编码系统.但这样做计算量与复杂性都将增加,人们经过大量的实验企图寻找一
种通用的算法,常见的有两种表达式比较接近实际:
一种模型是:
1 2| | | |( , )c cR eα βα β- -= (5.5.8)
其中,α和β分别代表水平与垂直方向的像素间距离像素间距离像素间距离像素间距离;1 2,c c则可根据不同的图像进行调整选
89
取的系数.
另一种模型是:
2 / 33 / 2 3 / 2( , ) ( )Rα β α γβ = + (5.5.9)
其中,α和β仍分别代表水平与垂直方向的像素间距离,γ则代表行间距离与水平方向像素
之间的距离之比值.
关于用作预测的像素选择,一般来说,1 2 3 4, , ,s s s s与0s相关性最强,因此用3到 4个预测
像素已有较高的准确性.再用更多的像素来预测,所得到的准确性改善就不很显著了.为了
使上述讨论有一个数量的概念,下面给山一组协方差函数与预测系数,预测误差的实验数据
(取自J.B. O'Neal于l966年发表的结果).预测像素分布结构如图5.5.3所示,该数据是对一
帧侧面头像取样的结果.
图图图图5.5.3 二维预测的邻近像素排列二维预测的邻近像素排列二维预测的邻近像素排列二维预测的邻近像素排列
协方差:
1 01 02 013
04 05 06 01
1, 0.934 0.960 0.905
0.919 0.832 0.814 0.829
R R R R
R R R R
= = = =
= = = =
下标"0i"表示表示0s与is之间的相关系数.
预测系数与预测误差如图5.5.1所示.
图图图图5.5.1 预测系数与预测误差预测系数与预测误差预测系数与预测误差预测系数与预测误差
90
还应该指山的是,以上讨论都是对单色图像或黑白电视图像进行的.至于彩色电视图像,
对于目前流行的体制,不管是PAL制还是NTSC制,色度信号都是由色副载波加到亮度信号
上去的.这就使全电视信号(或称复合电视信号)的样值之间的相关性,由于色度皮信号与亮
度信号的相位不一致而发生了变化.这就要重新设计它们的预测结构与预测系数,最直接的
办法就是将亮度信号与色度信号分开来处理,叫分量预测编码.由于它已用于专门技术,故
本课程不作详细讨论.但它所应用的基本原理还是本章叙述过的内容.
5.5.3 帧内预测编码帧内预测编码帧内预测编码帧内预测编码
帧间预测编码是对活动图像,特别是电视图像实施的一种压缩编码方法.它是利用活动
图像相邻帧间的相关性,即时间相关性来达到压缩数据的目的.显然,为了比较与处理,需
要能存储二帧以上图像的帧存储器,技术比较复杂,成本也较高,但效果较好.将电视信号
压缩至16Mbit/s的数码率,仍可保持较好的图像质量,因而有很大的吸引力.本节将介绍
帧间预测编码的基本思想及一些基本实现方法.
1. 帧间帧间帧间帧间预测的基本思想预测的基本思想预测的基本思想预测的基本思想
以电视信号为例来说明帧间预测的基本思想.现行的电视信号是由奇数场(由奇数行组成)
与偶数场(由偶数行组成)构成一帧.因此.帧间预测也包含有场间预测的含义.图5.5.4画出
了相邻三场及相邻像素的示意图,每个像素所标的坐标依次为(列,行,场), 当前场中待
预测的像素坐标为(0,0,0),坐标正方向如图所示.
图图图图5.5.4 相邻三场的相邻像素相邻三场的相邻像素相邻三场的相邻像素相邻三场的相邻像素
91
设I为所有用于预测的像素坐标的集合,三维矢量0=(0,0,0)为待预测像素坐标,则
0I.
对某帧图像来说,任何一个像素sn都可能是待预测像素,设{ | }s-∈n ii I是与sn 相邻的已
传输像素的集合,则sn的预测值可表示为:
^( )s a i s-
∈
=∑n n n i
i I
(5.5.10)
这里,i n都应为三维矢量.
类似于一维预测的讨论,关键是要选择( )ani使均方误差值2^( )E s s-n n最小,即设计出在
均方误差测量条件下的最佳预测器.设图像信号为广义平稳的随机过程,则预测器系数可由
以下方程的解得到:
( ) ( ) ( )s sa i R R
∈
- = ∈∑n
i I
i k k k I (5.5.11)
由此组系数组成的预测器,得到均方误差值为:
( ) ( ) ( )s sR a i R
∈
-∑0n
i I
i (5.5.12)
因为2( )s sRσ=0是信号的方差,故( ) ( )sa i R
∈
∑n
i I
i就表示预测使传输信号方差减小的量,即代表
了预测增益的大小.在极端情况下,若预测误差为0,则信号就不需传输,这表示一切都在"预
料之中",显然这对信息传输来说是无意义的.
式(5.5.10)与(5.5.11)包含了帧内预测(只用帧内像素预测),帧间预测(只用不同场或帧的像
素来预测)以及帧内帧间混合预测(同时用三维像素来预测).Haskell等人对不同预测器作了实
验研究,如表5.5.2所示.
表表表表5.5.2 不同的预测器的比较不同的预测器的比较不同的预测器的比较不同的预测器的比较
注:坐标参看图5.5.4
92
实验表明,一般来说,帧(场)间预测可得10 dB左右的增益,在运动不大的场面,帧间预
测性能较好,而运动较大的场面,还是帧内预测为好.在混合预测时,以场和豫素差的预测
为佳.
作为帧间预测的特例,对活动很慢的图像,可以用前一帧像素来代替(不是预测,即不计
算差值)当前帧,因而可作隔帧传输.对未传输的帧,可用简单重复的方法补上,或作插值补
上.由于人眼视觉特性,对静止或慢动图像要求较高的空间分辨率,而对时间分辨率要求饺
低.故上述做法带来的误差往往不易被觉察.
2. 运动补偿预测运动补偿预测运动补偿预测运动补偿预测
上面谈到人眼对画像的静止部分要求较高的空间分辨率与较低的时间分辨率,而对运动
部分则恰好相反.于是可以利用人眼的视觉特性,将图像分成静止部分与运动部分分别处理.
对静止部分只用前一帧数据的重复,对运动部分则设法测定其位移方向与位移量.以位移量
确定预测像素与预测误差,并对只作重复的静止部分作出补偿,构成当前帧的完整图像,故
称之为运动补偿预测.
运动补偿预测技术由以下几方面组成:
(1) 分割图像为静止的与运动的两部分;
(2) 估计运动部分的位移;
(3) 用位移估计进行运动补偿预测;
(4) 预测信息编码.
运动补偿预测的原理框图如图5.5.5所示,输入信号送至预测器,同时也加到分割器和位
移估计器.分割信息还送至位移估计器,二者又同时用控制预测器,并作编码传送.
图图图图5.5.5 运动补偿预测编码器原理运动补偿预测编码器原理运动补偿预测编码器原理运动补偿预测编码器原理
93
由上所述,可见图像分割是整个工作的基础,也是较为困难的一步.这里介始两种简单
的方法.一种称为位移估计法.即把图像分成若干矩形子块,适当选择块的大小,并确定子
块是动的或是不动的,估计出运动于块的位移,以便进行预测,另一种叫递归估计法,它是
对每个像素都进行递归估计,作出最佳逼近,确定位移.
(1) 位移估计匹配算法
匹配算法是在当前帧,比如说第k帧,取一矩形像素块(M N×个像素),与前一帧(第1k-
帧)相同位置的下一个M N×像素块作比较.若灰度完全一样,或变化最小,则该块为不动块.
否则要在前一帧(第1k-帧)图像中扩大范围再作比较.比如扩大至( 2 ) ( 2 )m mM d N d+ × +,其
中md表示在帧间隔时间内可能运动的最大位移,这要根据图物体运动速率来估算.若,取
10md=,则扩展区如图5.5.6所示.
图图图图5.5.6 扩展的匹配区域扩展的匹配区域扩展的匹配区域扩展的匹配区域
匹配程度的计算,可以用均方误差最小或绝对误差最小等准则来作比较.理论上应该在
匹配区内的所有可能取的M N×子块都与第k帧内的M N×子块亮度值作比较,作出差别最
小的判决.
(2) 递归估计算法
近年来对像素递归估计算法研究较多,出现了许多算法,比如Newton Raphson算法,
Bergmann算法,Robbins梯度算法等.这里不再较详细讨论.
应用帧间预测技术可获得高的数据压缩倍数,对窄带传输系统特别有吸引力,因而在会
议电视及可视电话等数字传输中很早就得到了的应用.20世纪未,由于微电子技术和计算机
94
技术的高速发展,使运动估计和运动补偿技术有了快速实现的可能,出现了许多实时,有效
的算法,从而使帧间预测技术在数字电视与高清晰度电视信号的数据压缩编码中得到了广泛
的应用.
第第第第五五五五章章章章 预测预测预测预测编码编码编码编码
对信号样值进行预测的概念在第2章讨论差值量化时已经提出.在那里,预测器只是一
个能输出预测信号的功能块,并指出对差值信号进行量化编码的比特率要低于对原信号的量
化编码的比特率.
本章将讨论预测器的工作原理,其中主要讨论线性预测的原理,算法以及它在语音,图
像信号编码中的应用.
由信号的取样值构成的信号序列中样值间存在相关性.利用相关性一个成功可追溯到
Wiener在1949年发表的论文"平稳时间序列的外推,内插,平滑及其在工程中的应用".由
于Wiener的开拓性工作,使滤波理论发展到一个新的阶段,即用统计学的观点来描述通信与
控制问题.整个滤波理论应包括对时间序列的滤波,预测与平滑三个方面.这些工作都是以
信号样值之间的相关性为依据的.如果在n时刻从噪声污染中复现出信号( )x n的近似值为
^( )x n,这便是滤波问题;如果在n时刻输出^( )x n k+,即对n + k时刻信号值的估计,则这是
预测问题;如果在n时刻输出^( )x n k-,即对n - k时刻信号的修正,则称作平滑问题.这里,
我们关心只是预测问题,但必须搞清楚时间序列分析的一些基本问题,诸如时间序列的概率
模型,时域估计等.
应用线性预测来进行数据压缩编码的一个种实现就是图像与语音信号的DPCM编码方
法.这个方法的侧重点在于寻找一个预测器,该预测器能够减小差值信号的方差,从而也能
够减小量化误差或降低数码率.在语音的线性预测编码中则另辟蹊径,用准周期脉冲与随机
噪声去激励一个线性时变系统所产生的输出作为语音主号.不同的语音所对应的线性时变系
统的参数也不同.应用线性预测的方法对这些时变的参数进行估值与编码来代替对语音波形
的直接编码,可以达到数据压缩的目的.这种能够用准周期脉冲与随机噪声合成语音信号的
时变参数系统,叫作线性预测声码器LPC(Linear Prediction Vocoder).目前用2400 b/s的速率
传输LPC参数已可获得较满意的话音质量,比较PCM编码所需64000 b/s的数据率,其压缩
比约为27.因此,这是一个很有效的数据压缩方法,在5.4节后面我们会更详细的讨论.5.5
节讨论图像线性预测编码的方法,主要讨论帧内与帧间的DPCM概念.
在讨论线性预测方法之前,先学习与介始随机过程的线性模型与线性预测的基本原理.
这是很重要的理论基础.上面已提到,它不仅应用于线性预测,而且也应用于整滤波理论与
2
现代谱分析研究.
5.1 时间序列的概率模型时间序列的概率模型时间序列的概率模型时间序列的概率模型
本节讲述时间序列的几种概率模型,它们总称为随机模型.根据信息率的观点,载有信
息的信号总是含有随机(也叫不确定)因素的.事实上,如果尚未听到声音就知道对方想说什
么,则从信息论意义上廛,这个语音信号没有意义,可以完全取消.但从其他方面考虑,比
如欣赏音乐,还是有意义的.
从数学上说,随机过程可定义为随机变量的一个集合{ ( ), }X t t T∈,其中T为定义随机过
程的时间点的集合.如果T是连续的(即,t-∞ <
= =
-
- <
∑
∑
(5.1.15)
由式(5.1.15)可见,自相关系数ρ是k的偶函数,且在k q=处被"截尾".这是MA过程的一
个特点.对于01β=的MA(1)过程,自相关系统为:
1
2
1
1 0
( ) 1
1
0
k
k k
β
ρ
β
=
= = ±
+
其他
(5.1.16)
3. MA过程的可逆性过程的可逆性过程的可逆性过程的可逆性
7
为了说明可逆性问题,先考察下面两个1阶的MA过程.
过程A:
1t t tX Z Zβ-= +
过程B:
1
1
t t tX Z Z
β-= +
由式(5.1.16)很容易看出,这两个不同的过程具有完全一样的自相关系数.这就说明,我们不
能从给定的自相关系数惟一地识另一个MA过程.
那么,{ }iβ满足什么条件下能使过程可逆呢 对此,让我们作进一步考察.对过程A和
B做如下改写:
过程A:
2
1 1 2 1 2( )t t t t t t t t tX Z Z Z Z Z Z Z Zβ β β β β- - - - -= - = - - = = - + -
过程B:
1 12
1 1
t t t tX Z Z Z
β β- -= - + -
如果| | 1β<,则A的级数收敛而B的级数不收敛,故称过程A是可逆的,而过程B不可逆的.
也就是说,条件| | 1β<使自相关系统能惟一对应一个MA过程.
现在来讨论一般情况下的可逆性条件.为了讨论方便,引入"后移算子"B,其定义为:
j
t t jB X X-= (5.1.17)
其中j为整数.于是,式(5.1.11)可写成,
0 1( ) ( )q
t q t tX B B Z B Zβ β β θ= + + + = (5.1.18)
仿照上面的改写方法,式(5.1.18)又可写成:
1( )t tZ B Xθ-= (5.1.19)
其中( )Bθ为B的q次多项式.若将B看作复变量而不再作算子,且设多项式,
0 1( )q
qB B Bθ β β β= + + + (5.1.20)
的根为1 2, , ,qb b b ,则多项式( )Bθ可以分解因式,并有:
8
1
1
1
1
( )
( )
q
i
q
ii
i
i
k
B
b B
b B
θ-
=
=
= =
-
-
∑
∏
(5.1.21)
其中ik是由1 2, , ,qb b b 确定的常数.
将式(5.1.21)代入式(5.1.19),并将1/( )ib B-展成,
2
2
1
1
i i i
B B
b b b
+ + +
得到
2
1
1
q
i
tt
ii i i
kB B
Z X
b b b=
= + + +
∑ (5.1.22)
现在应该恢复算子B的功能,可见式(5.1.22)中每个i值都对应1个无穷级数,q个无穷级数之
和代表了前面模型A与B的情况,它的可逆条件对应于无穷级数
2
2
1
i i
B B
b b
+ + +
的收敛条件,
即希望
1
1
ib
<或| |ib.于是有结论:
一个q阶的MA过程的可逆条件为:该过程的后移算子多项式( )Bθ的根,全部在B的平
均单位圆之外.
4.1.5 自回归过程自回归过程自回归过程自回归过程(AR模型模型模型模型)
假设{ }tZ是一个均值为0,方差为2
zσ的纯随机过程.如果有一个过程{ }tX满足:
1 1t t q t q tX X X Zα α- -= + + + (5.1.23)
则称该过程为p阶回归过程,记作( )AR p或称AR(Auto-Regressive Process)模型.这个模型很
象多元回归模型,但不是对于独立变量,而是对于tX的过去值的回归,故称之为自回归自回归自回归自回归.AR
模式在线性预测编码中起主角作用的一个线性模式,因此对它的讨论要深入一点.
1. AR(1)模型模型模型模型(Markov过程过程过程过程)
1阶自回归过程简写为AR(1)过程(模型),也叫Markov过程,即每个当前信号tX只与前
9
一时刻的信号值1tX-有关,其表达式为:
1t t tX X Zα-= + (5.1.24)
应用后移算子B(定义见式(5.1.17)),上式改写为:
(1 )t tB X Zα- =
因而有:
2 2
2
1 2
0
(1 )
1
t
tt
s
t t t t s
s
Z
X B B Z
B
Z Z Z Z
α α
α
α α α
∞
- - -
=
= = + + +
-
= + + + =∑
(5.1.25)
这就是MA(∞)的表达式.也就是说,1阶自回归过程可用无限多阶MA过程来表示.
随机变量tX的均值(1阶矩)为:
2[ ] [ ](1 )
1
z
t tE X E Z
μ
α α
α
= + + + =
-
| | 1α< (5.1.26)
在对Z作零均值假设条件下,可见[ ] 0tE X=.
2阶矩量有方差与自协方差,现分别计算如下:
2def
2 4 2
2
0 0
var[ ] var[ ](1 ) | | 1
1
cov[ , ] [ ]
z
t tx
s j
t t k t t k t s t k j
s j
X Z
X X E X X E Z Z
α
α α α α
α
α α
∞ ∞
+ + - + -
= =
= + + + = <
-
= =
==
∑ ∑
(5.1.27)
当j s k= +时,该两项相乘的均值为tZ的方差2
zσ;当j s k≠ +时,各项乘积的均值为0(因不
相关),故,
2
22
2
0
cov[ , ] ( )
1
k
s k s kz
t t k zz
s
X X k
α σ
γ σ α α α σ
α
∞
+
+
=
= = = =
-
∑ | | 1α< (5.1.28)
由此可见,1阶AR过程(Markov过程)的1阶矩与2阶矩都为常数,故该过程是2阶平
稳的,其条件为| | 1α
10
部分.由图可见,当α值减小时,自相关系数衰减加剧,当α为负值时,( )kρ变成交变状态,
即出负相关.负相关意味着两个随机变量中一个取较大值时另一个却取较小值.
从图中还可以看到AR模型的所谓"拖尾效应",即当k p≥(这里1p=)时,( )kρ并不为
0.虽然AR(1)模型的tX只与1tX-有关,但由于相关的传递性,1tX-又与2tX-有关,……,这
与MA模型的自相关系数的"截尾"特性相对应,构成了这两种模型的不同特点.
图图图图5.1.1 AR(1)过程的自相关函数例子过程的自相关函数例子过程的自相关函数例子过程的自相关函数例子
2. AR(p)过程过程过程过程
AR(p)过程的表达式见式(5.1.23),应用算子B仍可改写为:
1(1 )p
p t tB B X Zα α- - = (5.1.30)
我们所关心的问题,仍然是过程的平稳及过程中各样值之间的相关性.相关性问题将在5.4
节详细讨论,这里只讨论AR(p)过程的平稳条件.为此,令
1( ) (1 )p
pB B B α α= - - (5.1.31)
于是,式(5.1.30)可改写为:
( )t tB X Z = (5.1.32)
或者写成tX的表达式:
11
1( )t tX B Z -= (5.1.33)
所谓过程平稳的条件,就是要求表达式(5.1.33)有是限的,并且不是同时t的函数.由于式(5.1.31)
是一个以iα为系数的B的p次多项式,我们设该多项式有p个异实根1 2, , ,pb b b ,并进行类
似于式(5.1.21)与式(5.1.22)的讨论,即对( )B 作因式分解:
1
( ) ( )
p
i
i
B b B
=
= -∏ (5.1.34)
然后写出1( )B -为:
1
1
( )
p
i
ii
k
B
b B
-
=
=
-
∑ (5.1.35)
于是,式(5.1.33)可展开为:
2
2
1
1
p
i
tt
ii i i
kB B
X Z
b b b=
= + + +
∑ (5.1.36)
由此可见,过程tX平稳的条件,也就是上式中括号内无穷级数收敛的条件.显然,该条件就
是
1
1
ib
=的1 2,α α取值范围为:
2 1 2 1 2( ) 1,( ) 1, 1 1α α α α α+ < - < - < 时,用上述方法来确定AR过程的平稳域,会遇到计算上的困难,要寻
求更简单的判别方法.下面介绍常用一个方法.
3. Schur-Cohn准则准则准则准则
这里介绍的是一个比较简单的平稳性判别准则.由于从数据压缩应用出发来确定的模型
系数必须使系统能稳定工作,因此,我们主要关心的就是AR(p)模型的系数iα与系统稳定性
关系.
AR(p)模型所对应系数稳定的条件的充分必要条件就系统的特征多项式
1
1 1 0( )p p
p pF a a a aω ω ω ω-
-= + + + (5.1.38)
的根在复平面的单位圆内.这个结论可从( )B 的根的条件来推证.这里应该注意式(5.1.38)
与( )B 的区别,即i p iaα-=.
Schur-Cohn准则的关键是构造所谓的S-C行列式,它的定义为:
13
01 1
1 02
1 0
0 1 1
10 1
10
0
0
0
0
p p p k
p p k
kp
k
pk
p pk
p k p
a a a a
a a a a
a a a
a a a a
a a a a
a a a
- - +
- +
-
-
--
- +
=
,1, 2, ,k p= (5.1.39)
其中,ia为ia的共轭复数,且令01 = -.
S-C准则准则准则准则:( )Fω为稳定系统的特征多项式的充要条件是:它的S-C行行形式满足:
0
0k
k
k
为奇数
为偶数
(5.1.40)
显然,满足上述条件,意味着( ) 0Fω=的所有根都在单位圆内;否则,至少有一个根不
在单位圆内,但不能判定它们在圆外还是在圆上.
我们不去证明准则的正确性,只是举一个用它对AR(2)过程所对应的系统的稳定性判断
的例子,来说明它的应用.
AR(2)模型所对系统的特征多项式为:
2
2 1 0( )F a a aω ω ω= + +
S-C行列式为:
0 2
1
2 0
a a
a a
=
因为式(5.1.37)中的ia为实数,故这里i ia a=.然后用0
2 1 1 2, , 1a a aα α= - = - =代入上式,可得:
2
1 21α = -
由条件(5.1.40)应有10 <.因此,要求21 1α- +,不等式两边消去2
2(1 )α+后,
有:
2 2
2 1(1 )α α- >
得到稳定性条件为:2 11α α± 乘等式的两边,再各自求均值,便得到自协方差函数的关系式:
1 2( ) ( 1) ( 2) ( )pk k k k pγ α γ α γ α γ= - + - + + - (5.2.21)
相应的自相关函数的表达式为:
1 2( ) ( 1) ( 2) ( )pk k k k pρ α ρ α ρ α ρ= - + - + + - 1, 2 ,k p= (5.2.22)
这便是由p个方程组成的线性方程组,称之为Yule-Walker方程.这是一个在时间序列分析中
非常有用的方程,它反映了线性模型参数{ }iα与自相关系数{ }kρ之间的关系.写成矢量与矩
阵形式为:
=ρ αρ αρ αρ αR (5.2.23)
式中,
( (1), (2), , ( ))Tpρ ρ ρ= ρρρρ (5.2.24a)
1 2( , , , )T
pα α α= αααα (5.2.24b)
1 (1) ( 1)
(1) 1 ( 2)
( 1) ( 2) 1
p
p
p p
ρ ρ
ρ ρ
ρ ρ
-
- =
- -
R (5.2.24c)
可以看以,R是一个Toeplitz矩阵.
21
2. AR模型参数的估计模型参数的估计模型参数的估计模型参数的估计
若取得了信号的N个样值1 2, ,Nx x x ,它们是随机序列tX的一段样本值,则称N为样本
长度.为根据该段样本数据来估计tX,先要对自协方差( )kγ与自相关函数( )kρ作出估计.
我们采用:
1
1
^ ^( ) ( )
N k
l l k
l
k k x x
N
γ γ
-
+
=
= - =∑ 0,1, 2, , 1k N= - (5.2.25)
及
^( )
^ ^( ) ( )
^(0)
k
k k
γ
ρ ρ
γ
= - = 0,1, 2, , 1k N= - (5.2.26)
也有采用如下定义的:
* *
1
1
^ ^( ) ( )
N k
l l k
l
k k x x
N k
γ γ
-
+
=
= - =
-
∑ 0,1, 2, , 1k N= - (5.2.27)
从式(5.2.25)和(5.2.27)可以发现,
*^ ^( ) ( )
N k
k k
N
γ γ
-
=
因此,可见它们是渐近相等的,即当N→ ∞时有*^ ^( ) ( )k kγ γ=.
有了样本的自相关系数,便可由式(5.2.23)解出:
1^^^-=α ρα ρα ρα ρR (5.2.28)
展开得:
1
1
2
^^1 (1) ( 1) (1)
^^(1) 1 ( 2) (2)
^^( 1) ( 2) 1 ( )p
p
p
p p p
αρ ρ ρ
αρ ρ ρ
αρ ρ ρ
--
- =
- -
(5.2.29c)
常常称^αααα为αααα的Yule-Walker的估计.
为了得到残差能量的估计,将式(5.2.20)移项,两边平方后再取均值,得:
2 2
1 , 1
[ ] (0) 2 ( ) ( )
p p
z t j i j
j i j
E Z j j iσ γ α γ α α γ
= =
= = - + -∑ ∑
为了化简上式,可先把(5.2.23)按式(5.2.24)展开,并对第j个方程两边乘以jα;然后将方程组
的左边与右边分别加起来,这样就可得到:
22
1 , 1
( ) ( )
p p
j i j
j i j
j j iα γ α α γ
= =
= -∑ ∑ (5.2.30)
将上式代入2
zσ表达式,得:
2
1 , 1
(0) ( ) (0) ( )
p p
z j i j
j i j
j j iσ γ α γ γ α α γ
= =
= - = - -∑ ∑ (5.2.31)
将估计值^( )j iγ-与^
iα代入(5.2.31)便得到2
zσ的估计值2^
zσ:
2
1 , 1
^ ^ ^ ^ ^ ^ ^^(0) ( ) (0) ( )
p p
z j i j
j i j
j j iσ γ α γ γ α α γ
= =
= - = - -∑ ∑ (5.2.32)
式(5.2.29)与式(5.2.32)就是AR模型参数矩估计公式.可以发现它们是比较简单易行的.
下面用AR(1)与AR(2)模型的矩估计来说明.
AR(1)参数估计:
1
^(1)
^ ^(1)
^(0)
γ
α ρ
γ
= = (5.2.33a)
22
1
^ ^ ^ ^ ^(0) (1) (0)[1 (1) ]zα γ γ α γ ρ= - = - (5.2.33b)
AR(2)参数估计:
12
^^(1)[1 (2)]
^
^[1 (1) ]
ρ ρ
α
ρ
-
=
-
(5.2.34a)
2
22
^ ^[ (2) (1) ]
^
^[1 (1) ]
ρ ρ
α
ρ
-
=
-
(5.2.34b)
2
1 2
^ ^ ^ ^ ^(0)[1 (1) (2)]zα γ α ρ α ρ= - - (5.2.34c)
其中的^( )kγ与^( )kρp(k)由式(5.2.25)与(5.2.26)给出.
前面已指出AR模型的矩估计与最小二乘估计是很接近的.虽然Yule-Walker方程的导出
过程中并未引入"最小误差"的条件.但若以最小均方误差条件来推导参数估计关系,也得
到相同的结果.若以相同的样本序列(设N个样值)来作以上的两种估计,比如作出AR(1)或
AR(2)估计,就可发现其结果相差甚微,只是求和限有点差别.当N→ ∞,这个差别可略.
23
5.2.3 序列最小二乘法序列最小二乘法序列最小二乘法序列最小二乘法
上两节介绍了线性预测参数的两种估计方法.从原理上说.这两种方法是简单的而且有
效的,但从计算速度与实时处理这两方面来说,并不理想.这是因为在估值过程中需要计算
如方程式(5.2.12)与(5.2.28)中的逆矩阵,而逆矩阵的计算往往是复杂费时的;这两种方法都是
属于成批型的算法,只有在取得足够多(N足够大)的样值时,估值才有效,这对数据的实时
处理要求很不利,且要求设备有较大的存储容量.为克服上述困难,现在来介绍运用自然递
批的方法.本节讨论的是称作序列最小二乘(SLS)算法.
1. 最小二乘估计的递推结构最小二乘估计的递推结构最小二乘估计的递推结构最小二乘估计的递推结构
我们从最小二乘估计的结果出发来推导自然递推关系.为了避开矩阵的求逆运算.先对
式(5.2.12)作一些改造.为此,将式(5.4.11)的左右两边分别改写成:
1,
1
1,1 1,2 1,1
2,1 2,2 2,2 ,
1
,1 ,2 ,
,
1
1 1,1
1 2,1
1 ,1
n
t i ti
i
t t t ntn
t t t nt t i tiT
it
t p t p t p ntnn
t p i ti
i
t t
t t
t t p
x x
x x xx
x x xx x x
x x xx
x x
x x
x x
x x
-
=
- - -
- - --
=
- - -
-
=
-
-
-
= =
=
∑
∑
∑
x x
2 1,2 1,
2 2,2 2,
2 ,2 ,
t t tn t n
t t tn t n
t t p tn t p n
x x x x
x x x x
x x x x
- -
- -
- -
+ + +
若令
1, 2, ,[ , , , ]T
i t i t i t p ix x x- - -= x (5.4.35)
则有,
1
n
T
t ti i
i
x
=
=∑x x x (5.2.36)
以及
24
2
1, 1, 2, 1, ,
1 1 1
2
2, 1, 2, 2, ,
1 1 1
2
, 1, , 2, ,
1 1 1
n n n
t i t i t i t i t p i
i i i
n n n
t i t i t i t i t p iT
i i it
n n n
t p i t i t p i t i t p i
i i i
x x x x x
x x x x x
x x x x x
- - - - -
= = =
- - - - -
= = =
- - - - -
= = =
=
∑ ∑ ∑
∑ ∑ ∑
∑ ∑ ∑
x x
这是一个对称矩阵,故在式(5.2.16)的推证中1( )T
t
-x x的转置并不作倒序排列.若应用式(5.2.35)定义
的向量ix,则上式可简化为:
def
1
1
n
T T
t i i n
i
-
=
=∑==x x x x P (5.2.37)
于是式(5.2.11)与式(5.2.13)分别可写成:
1
1
^
n
n ti i
i
x-
=
=∑ααααP x (5.2.38)
与
1
^
n
n ti i
i
x
=
=∑ααααP x (5.2.39)
为了得到自然迭代关系,将式(5.2.38)左边分成两项,即,
1
1
1
^
n
n ti i tn n
i
x x
-
-
=
= +∑ααααP x x (5.2.40)
用1n-去代替n,式(5.2.38)仍应成立,即,
1 1
1
1 1
^
n n
T
i i n- ti i
i i
x
- -
= =
=
∑ ∑ααααx x x (5.2.41)
这里用1
^
n-αααα来表示1n-组观测值得到的参数估计,以区别于n组观测值得到的估计^
nαααα.下标n
也是为了区别才加上去的.应该强调的是,^
nαααα与前面公式中的^αααα是完全一样.
把式(5.2.41)代入式(5.2.40),便有:
1
1
1
1
^ ^
n
T
n n i i n- tn n
i
x
-
-
=
= +
∑α αα αα αα αP x x x (5.2.42)
若对式(5.2.42)左边同时加减一项1
^T
n n n-ααααx x,并应用1
n
-P的定义式(5.2.37),便得递推关系:
25
1 1
1 1
^ ^^( )T
n n n n- n tn n n-x- -= + -α α αα α αα α αα α αP P x x (5.2.43)
或者:
1 1
^ ^^( )T
n n- n n tn n n-x= + -α α αα α αα α αα α αP x x (5.2.44)
从上面的递推关系可以看到,估值^
nαααα可以由1
^
n-αααα,也就是它前面的估值项(在减少一个数
据集的基础上)以及一个修正项来计算.这个修正项是使用前面的估值引起的预测误差
1
^T
tn n n-x-ααααx的一个函数.应该指出,尽管这是个递推形式,但它仍然属于最小二乘估值算法,
因为上面的一切推导都是由最小二乘估值结果出发的.同时应注意式(5.2.44)远不理想,因为
从1
n
-P的定义来计算nP.仍然有逆矩阵运算问题.下面来解决这个问题.
2. 逆矩阵定理逆矩阵定理逆矩阵定理逆矩阵定理
所谓逆矩阵定理,就是不用计算逆矩阵而得到矩阵求逆的结果的一个递推算法.
定理定理定理定理 由式(5.2.37)定义的1
n
-P,它的逆只需由1n-P及nx(由式(5.2.35)所定义)按下述关系
递推得出,即:
1 1
1
11
T
n n n n
n nT
n n n
- -
-
-
= -
+
P x x P
P P
x P x
(5.2.45)
上式的证明这里不再讲述.
由式(5.2.44)与(5.2.45)构成了一组递推条件,并称此为序列最小二乘估计算法.该算法的
关键是用由式(5.2.35)定义的向量ix作为已知条件代原来已知矩阵x.而ix实际上是分别由各
组观测值组成的,从而构成一组一组递推估计的流程,即由第一组(1i=)观测值来作αααα的第一
—次估计1
^αααα;然后由1
^αααα与第二组(2i=)观测值作第二次估计2
^αααα;….那么流程的初值即作第
一次递推估计时的0
^αααα与0P应取什么值呢 下面分别加以说明.
3. 初始值初始值初始值初始值0
^αααα与与与与0P
初始值0
^αααα与0P是递推算法所必需的已知条件.然而在我们讨论的问题中,它们却是先验
未知的,所以还需对它们进行估值.
当对^αααα缺乏先验知识时,初始参数的估值0
^αααα可以设置为零向量:
0
^=0αααα (5.2.46)
26
关于0P的初始化值可以设置为一对角矩阵:
0β=P I 1β>> (5.2.47)
这里I为单位矩阵.为什么β要比1大很多 这是这样考虑的:从定义式(5.2.37)的求和限来看,
n最小为1,也就是说,1
0
-P应为0矩阵才行.但从迭代关系式(5.2.45)来看,对0P的任何实际
可行的估值都必须是有限的.因此,设置1
0ε-=P I和1ε,所有的( )k
kα都应为零,即偏相关函数是"截尾"的.还可以进一步证
明偏相关函数的截尾性也是平稳时间序列为AR模型的充分条件,即它是AR过程所持有的
标志.
利用AR过程的截尾特性,可以作出结论:"最佳"的阶数可被确定为这样的k值,这种
k值之后的{ }pπ的非常接近于0,此时p k=.
3. 与与与与Durbin递推程序同时进行判决递推程序同时进行判决递推程序同时进行判决递推程序同时进行判决
在Durbin递推算法公式(5.2.49)和(5.2.51)中,我们定义:
( ) ( )
1
(0) ( )
k
k k
j
j
E r r jα
=
= -∑ (5.2.74)
若信号属于AR(p)过程,是可以证明:
( )limp
N
N
E E
→∞
=
而NE定义为均方差(N为样值序值的样值个数)
2
1
0 1
1pN
N k i k i
k i
E x x
N
α
-
-
= =
= -
∑ ∑ (5.2.75)
从式(5.2.48),有1k=时,
(1)
1 1
(1)
(0)
K
γ
α
γ
= =
由式(5.2.75)定义,可得:
(1) (1)
1(0) (1)Eγ α γ= -
2k=时,有:
35
( 2) 2 (1)
2[1 ]E K E= -
此结果表时,由式(5.2.74)定义的( )kE也是具有递推特性.因此,算法式(5.2.51)可改写成:
1
( 1)
1( )
( 1)
( ) ( )
k
k
j
jk
k kk
k k j
K
E
γ α γ
α
-
-
=
-
- -
= =
∑
(5.2.76a)
( ) ( 1) ( 1)k k k
k j k k jKα α α- -
-= - (5.2.76b)
( ) 2 ( 1)[1 ]k k
kE K E-= - (5.2.76c)
式中,1, 2, , ; 1, 2, , 1k p j k= = - .
递推初值为:
(0 )(0)Eγ= (5.2.77a)
(1)
1 1
(1)
(0)
K
E
γ
α= = (5.2.78b)
(1) 2 ( 0)
1(1 )E K E= - (5.2.78c)
这是Durbin递推算法的另一种表示.我们的目的是借助于递推过程中的( )kE进行最佳阶数的
判定.判定准则有多种,这里只介绍用得较多的三种(它们之间并没有原则的差别).
(1) FPE准则
最终预测误差准则FPE(Final Prediction Error Criterion)是把使下式为最小值的k作为AR
模型的最佳阶数:
( )1
FPE( )
1
kN k
k E
N k
+ +
=
- -
(5.2.78)
般情况下,( )kE是k的单调递减函数,但系数
1
1
N k
N k
+ +
- -
是k的递增函数. 于是会在k的某个
值,FPE(k)出现最小值,这个k值就被定为AR模型的最佳阶数.
〔2) AIC准则
信息理论准则AIC (1nformtion Theoretic Criterion)是把使下式为最小值的k,用为AR模型
的最佳阶数:
( )2
AIC( ) ln[ ]kk
k E
N
= + (5.2.79)
式(5.2.79)与式(5.2.78)实质上是一样的.因为当N很大时,FPE的对数收敛于AIC,即有:
36
lim[ln FPE( )] AIC( )
N
k k
→∞
= (5.2.80)
这个极限关系是很易验证的.
(3) CAT准则
自回归传递函数准则CAT(Autoregressive Transfer Function Criterion)是把使下式为最小值
的k作为AR模型的阶数:
( ) ( )
1
1
CAT( )
k
j k
j
N k N k
k
N NE NE=
- -
=
∑ (5.2.81)
它相当于在估计滤波器的滤波误差的均方值上给出最佳阶数.
应该说,AR模型的准确阶数的判定,对于谱分析来说,比线性预测应用更为重要.在线
性预测中,阶数的不准确,引起的误差并不会大大,只要保持自协方差函数相差不太大;即
可.这就要特别注意在计算样本的自协方差函数,偏相关函数以及模型识别与参数化估计时,
—定要设法将均值项去掉.去均值项的常用办法是从ix中减去它的样本均值:
1
1
^
N
i
i
x
N
μ
=
=∑
即用^
t ty xμ= -作tx参数估计.但如果tx 的期望值( )tE X本身为零,再减去^μ就会降低估计
精度,可见.在对数据序列1 2, , ,Nx x x 作时序分析时,判断( )tE X是否为0是很有必要的.
5.2.7 窄带平稳随机过程窄带平稳随机过程窄带平稳随机过程窄带平稳随机过程的的的的MATLAB实现实现实现实现
这里窄带平稳随机过程定义为信号功率谱密度是一个带通型且带宽远小于中心频率的随
机过程,通常窄带平稳随机过程( )X t可以展开如下形式:
2( ) ( )cos 2 ( )sin 2 Re[ ( ) ]cj f t
c c s c lX t X t f t X t f t X t eππ π= - = (5.2.82)
其中,( ) ( ) ( )l c sX t X t jX t= +是等效基带信号,( ), ( )c sX t X t分别是平稳随机,其功率谱密度是
基带型.一个典型的窄带平稳随机过程是高斯白噪声经过带通系统的输出,可以证明,窄带
平稳随机过程具有如下性质:
(1) 如果( )X t平稳,则( ), ( )c sX t X t也平稳;
(2) 如果[ ( ) 0E X t=,则[ ( )] [ ( )] 0c sE X t E X t= =;
37
(3) 如果( )X t是高斯平稳过程,则( ), ( )c sX t X t也是高斯平稳过程;
(4) 如果( )X t是高斯平稳过程用均值为0,则( ), ( )c sX t X t相互正交,且功率与( )X t相同.
例例例例5.2.2 高斯白噪声是一个功率谱密度为常数,且各时刻满足高斯分布的随机过程,设
高斯白噪声的双边功率谱密度为0/ 2N,现将高斯白噪声经过一个中心频率为cf,带宽为B的
带通滤波器,得到的输出信号即为窄带随机过程.
(1) 用Matlab产生高斯平稳窄带随机过程,010, 1Hz, 1fc B N= = =;
(2) 画出窄高斯随机过程的等效基带信号;
(3) 求出窄高斯随机过程的功率,等效基带信号的实部功率,虚部功率.
解解解解::::窄高斯随机过程的Matlab的程序为:zppw.m
*********************************************************
clear all
close all
N0=1; % 单边功率谱密度
fc=10; % 中心频率
B=1; % 带宽
dt=0.01;
T=100;
t=0:dt:T-dt;
% 产生功率功率为N0*B的高斯白噪声
P=N0*B;
st=sqrt(P)*randn(1,length(t));
% 将上述白噪声经过窄带带通系统
[f,sf]=T2F(t,st); % 高斯信号频谱
figure(1)
plot(f,abs(sf)) % 高斯信号的幅频特性
38
[tt,gt]=bpf(f,sf,fc-B/2,fc+B/2); % 高斯信号经过带通系统
glt=hilbert(real(gt)); % 窄带信号的解析信号,调用hilbert函数
glt=glt.*exp(-j*2*pi*fc*tt); % 得到解析信号
[ff,glf]=T2F(tt,glt);
figure(2)
plot(ff,abs(glf));
xlabel('频率 (Hz)')
ylabel('窄带高斯过程样本的幅频率特性')
figure(3)
subplot(411)
plot(tt,real(gt));
title('窄带高斯过程样本')
subplot(412)
plot(tt, real(glt).*cos(2*pi*fc*tt)-imag(glt).*sin(2*pi*fc*tt))
title('由等效基带重构的窄带高斯过程样本')
subplot(413)
plot(tt, real(glt))
title('窄带高斯过程样本的同相分量')
subplot(414)
plot(tt,imag(glt));
xlabel('时间 t(秒)')
title('窄带高斯过程样本的正交分量')
% 求窄带高斯信号功率;注:由于样本的功率近似等于随机过程的功率,因此可能出现一些偏差
P_gt=sum(real(gt).^2)/T;
P_glt_real=sum(real(glt).^2)/T;
39
P_glt_imag=sum(imag(glt).^2)/T;
% 验证窄带高斯过程的同相分量,正交分量的正交性
a=real(glt)*(imag(glt)')/T;
********************************************************
用到的子函数bpf.m如下.
*******************************************************
function [t,st]=bpf(f,sf,B1,B2)
df=f(2)-f(1);
T=1/df;
hf=zeros(1,length(f));
bf=[floor(B1/df):floor(B2/df)];
bf1=floor(length(f)/2)+bf;
bf2=floor(length(f)/2)-bf;
hf(bf1)=1/sqrt(2*(B2-B1));
hf(bf2)=1/sqrt(2*(B2-B1));
yf=hf.*sf.*exp(-j*2*pi*f*0.1*T);
[t,st]=F2T(f,yf);
*******************************************************
用到的子函数F2T.m如下.
*******************************************************
function [t,st]=F2T(f,sf)
df=f(2)-f(1);
Fmx=(f(end)-f(1)+df);
dt=1/Fmx;
N=length(sf);
T=dt*N;
40
t=0:dt:T-dt;
sff=fftshift(sf);
st=Fmx*ifft(sff);
*******************************************************
程序运行结果如图5.2.2所示.
0102030405060708090100
-0.5
0
0.5
窄窄窄窄窄窄窄窄
0102030405060708090100
-0.5
0
0.5
由由由由窄由由由窄窄窄窄窄窄窄窄
0102030405060708090100
-0.5
0
0.5
窄窄窄窄窄窄窄窄由窄窄窄窄
0102030405060708090100
-0.5
0
0.5
时时 t(秒)
窄窄窄窄窄窄窄窄由窄窄窄窄
图图图图5.2.2 窄带高斯过程窄带高斯过程窄带高斯过程窄带高斯过程
5.3 通用的预测编通用的预测编通用的预测编通用的预测编码码码码(量化编码方法量化编码方法量化编码方法量化编码方法)
音频类信源(包括语音,音乐,广播等音频信号在内)的特点是发出连续信息,因此其编
码方法与字符信源有所不同.首先要在时间上进行抽样离散化过程,然后在数值(幅度)范围
内进行量化,量化后的信号就是多进制的数字化信息;最后才进行数字信号编码,其输出的
二进制代码(符号)在信道内传输.简单地说,对发出连续信息的音频信源要经过抽样,量化
和编码三个过程,才能进行信息传输.在这些信号处理过程中都会有信息受到损失,而这些
信息的损失又是无法弥补的,因此,对于音频信源来说,信源编码是属于限失真编码,其编
码效率受到信息率失真函数( )R D的限制,无法实现像字符信源那样的最佳编码效果.
抽样与量化,我们在第二章已经讨论过了,限失真编码原理在第三章讨论过了,这里我
们根据上述基础与原理,来讨论量化的编码方法.
41
5.3.1 脉冲编码调制脉冲编码调制脉冲编码调制脉冲编码调制
前面已经指出,模拟信息源输出的模拟信号经抽样和量化后得到的输出电平序列
{ ( )}q sm kT,可以利用M进制PAM直接进行传输,也可以将每一个量化电平用编码方式传输.
所谓编码是把量化后的信号变换成代码,其相反的过程称为译码.
将模拟信号的抽样量化值变换成代码,称为脉冲编码调制(PCM).图5.3.1和表5.3.1给
出了脉冲编码调制的一个实例.假设模拟信号m(t)的最大值| ( ) | 4m t时,出"1"码;反之出"0"码.由于在
13折线法中用了7位二进代码来代表段落和段内码,所以对一个输入信号的抽样值需要进行
7次比较.每次所需的标准电流wI均由本地译码电路提供.
本地译码电路包括记忆电路,7/11变换电路和恒流源.记亿电路用来寄存二进代码,
因除第一次比较外,其余各次比较都要依据前几次比较的结果来确定标准电流wI值.因此,
7位码组中的前6位状态均应由记忆电路寄存下来.7/11变换电路就是前面非均匀量化中谈
到的数字压缩器.因为采用非均匀量化的7位非线性编码等效于11位线性码,而比较器只能
编7位码,反馈到本地译码电路的全部码也只有7位.因为恒流源有11基本权值电流支路,
需要11个控制脉冲来控制,所以必须经过变换,把7位码变成11位码,其实质就是完成非
线性和线性之间的变换.恒流源用来产生各种标准电流值.为了获得各种标准电流wI,在恒
流源中有数个基本权值电流支路.基本的权值电流个数与量化级数有关,如上例中,128个
量化级需要编7位码,它要求11个基本的权值电流支路,每个支路均有一个控制开关.每次
该哪几个开关接通组成比较用的标准电流wI,由前面的比较结果经变换后得到的控制信号来
控制.
附带指出,保持电路的作用是保持输入信号的抽样值在整个比较过程中具有一定的幅度.
由于逐次比较型编码器编7位码(极性码除外)需要将sI与wI比较7次,在整个比较过程中部
应保持输入信号的幅度不变,故需要采用保持电路.下面我们通过一个例子来说明编码过程.
例例例例5.3.1 设输入信号抽样值为+1270个量化单位(这里量化单位是指以输入信号归一化的
1/2048为单位),采用逐次比较型编码将它按照13折线A律特性编成8位码.
设码组的8位码分别用1 2 3 4 5 6 7 8c c c c c c c c表示.编码过程如下:
(1) 确定极性码1c
因输入信号抽样值为正,故极性码11c=
(2) 确定段落码2 3 4c c c
13折线中,正半部分的8个段落以1/2048为单位的每个段落的起点电平如表5.3.5所
示.由于段落码中的2c是用来表示输入信号抽样值处于8个段落的前四段还是后四段的,故
输入比较器的标准电流应选择为128wI=个量化单位.现在输入信号抽样值1270sI=个量化单
47
位,大于标准电流,故第一次比较结果为s wI I>J,所以21c=.它表示输入信号抽样值处于8
个段落中的后四段(5~8段).
表表表表5.3.5 段落起点电平段落起点电平段落起点电平段落起点电平
3c用来进一步确定它属于5~6段还是7~8段.因此,标准电流应选择为512wI=个量化
单位.第二次比较结果为s wI I>,故31c=.它表示输入信号属于7~8段.
同理,确定4c的标准电流应为1024wI=个量化单位.第三次比较结果为s wI I>,故41c=.
由以上三次比较得段落码为"111",输入信号抽样值1270sI=个量化单位应属于第八段.
(3) 确定段内码5 6 7 8c c c c
由编码原理已经知道,段内码是在已经确定输入信号所处段落的基础上,用来表示输入
信号处于该段落的哪一量化级的.5 6 7 8c c c c的取值与量化级之间的关系见表5.3.4.上面已经确
定输入信号处于第八段,该段中的16个量化级之间的间隔均为64个量化单位,故确定5c的
标准电流应选为:
8 ( 1024 8 64 1536wI= + × = + × =段落起点电平 量化级间隔) 个量化单位
第四次比较结果为s wI I,故81c=.它说明输入信号处于第八段中的3量化级.
经上述七次比较,编出的8位码为11110011.它表示输入抽样位处于第八段3量化级,
其量化电乎为1216个量化单位,故量化误差等于54个量化单位.顺便指出,除极性码外的
7位非线性码组1110011,相对应的11位线性码组为10011000000.
上面较详细地讨论了逐次比较型编码器的原理,下面再来讨论PCM信号的译码原理.常
用译码器大致可分为三种类型:电阻网络型,级联型,级联一网络混合型等.这里仅就图5.3.4
所示的电阻网络型译码器加以讨论.
图图图图5.3.4 电阻网络型译码器电阻网络型译码器电阻网络型译码器电阻网络型译码器
电阻网络型译码器与逐次比较型编码器中的本地译码器基本相同.从原理上说,两者都
是用来译码,但编码器中的译码,只译出信号的幅度,不译出极性;而收端的译码器在译出
信号幅度值的同时,还要恢复出信号的极性.电阻网络型译码器各部分电路的作用简述如下.
记忆电路用来将接收的串行码变为并行码,故又称为"串/并变换"电路.7/11变换
电路用来将表示信号幅度的7位非线性码转变为11位线性码.极性控制电路用来恢复译码后
的脉冲极性.寄存读出电路把寄存的信号在一定时刻并行输出到恒流源中的译码逻辑电路上
去.使产生所需要的各种逻辑控制脉冲.这些逻辑控制脉冲加到恒流源的控制开关上,从而
驱动权值电流电路产生译码输出.
由以上电阻网络则译码器各部分电路的作用.不难理解这种译码器的工作原理.其译码
过程就是根据所收到的码组(极性码除外)产生相应的控制脉冲去控制恒流源的标推电流支
路,从而输出—个与发送端原抽样值接近的脉冲.该脉冲的极性受极性控制电路控制.
2. PCM系统的抗噪声性能系统的抗噪声性能系统的抗噪声性能系统的抗噪声性能
上面我们较详细地讨论了脉冲编码调制的原理.现在,将要分析图5.3.2所示的PCM系
统的抗噪声性能.由该图可以看出,接收端低通滤波器的输出为:
49
0
^( ) ( ) ( ) ( )q em t m t n t n t= + +
式中:0( )m t —— 输出信号成分;
( )qn t —— 由量化噪声引起的输出噪声;
( )en t —— 由信道加性噪声引起的输出噪声.
为了衡量PCM系统的抗噪声性能,通常将系统输出端总的信噪比定义为:
2
0 0
2 2
0
[ ( )]
[ ( )] [ ( )]q e
S E m t
N E n t E n t
=
+
(5.3.1)
式中 E —— 求统计平均.可见,分析PCM系统的抗噪声性能时,需要考虑量化噪声和
信道加性噪声的影响.不过由于量化噪声和信道加性噪声的来源不同,而且它们互不依赖,
故可以先讨论它们单独有在时的系统性能,然后再分析系统总的抗噪声性能.我们先分析仅
考虑量化噪声的系统性能.假设发送端采用理想冲激抽样,则PCM系统输出端平均信号量化
噪声功率经为:
2
20 0
2
[ ( )]
[ ( )]q q
S E m t
M
N E n t
= = (5.3.2)
式中,qN为量化噪声功率,M量化级数.
对二进制编码,式(5.3.2)又可写成:
202N
q
S
N
= (5.3.3)
式中,N—— 二进制代码位数.由上式可见,PCM系统输出端平均信号量化噪声功率比将
仅依赖于每一个编码组的位数N.上述比值将随N按指数增加.大家知道,对于一个频带限
制在Hf的信号,按照抽样定理,此时要求每秒钟最少传输的抽样脉冲数等于2Hf;若PCM
系统的编码位数为N,则要求系统每秒传输2HNf.为此,这时的系统总带宽B至少等于HNf.
故式(5.3.3)还可写成:
2
02H
B
f
q
S
N
= (5.3.4)
由此可见,PCM系统输出端的信号量化噪声功率比还与系统带宽B成指数关系.
下面我们来分析信道加性噪声对PCM系统性能的影响.由于信道中加性噪声对PCM信
50
号的干扰,将造成接收端判决器判决错误,二进制"1"码可能误判为"0"码,而"0"码可
能误判为'1"码,其错误概率将取决于信号的类型和接收机输入端的平均信号功率比.因为
PCM信号中每一码组代表着一定的量化抽样值,所以其中只要发生误码,接收端恢复的抽样
值就会与发端原抽样值不同.通常我们只需要考虑仅有一位错码的码组错误,而多于一个错
码的码组错误可以不予考虑.
仅考虑信道加噪声时PCM系统的输出信噪比为:
01
4e e
S
N P
= (5.3.5)
式中,eP为误码率,由上式可见,由误码引起的信噪比与误码率成反比.
前面已经指出,传输模拟信号的PCM系统的性能用接收端输出的平均信器噪功率比来度
量.将式(5.3.2),(5.3.3)和(5.3.5)代入式(5.3.11)得:
22 2
0 0
2 2 2 2
0
[ ( )]2
[ ( )] [ ( )] 1 4 2 1 4 2
N
N N
q e e e
S E m tM
N E n t E n t P P
= = =
+ + +
(5.3.6)
在接收端输入大信噪比条件下,即当24 2 1N
eP<>时,
2
0
2
0
2 1
4 2 4
N
N
e e
S
N P P
≈ = (5.3.8)
在用基带传输的PCM中继系统中,通常使误码率降到610-是容易实现的,此时,可按式(5.3.8)
来估计PCM系统的性能.
5.3.2 增量调制增量调制增量调制增量调制
增量调制( M)是在PCM方式的基础上发展起来的另一种模拟信号数字传输的方法. M
可以看成PCM的一个待例,因为它们都是用二进制代码形式去表示模拟信号的方式.但是,
在PCM中,信号的代码表示模拟信号的抽样值,而且,为了减小量化噪声,一般需要较长的
代码及较复杂的编译码设备.而 M是将模拟信号变换成仅由一位二进制码组成的数字信号
序列,并且在接收端也只需要用一个线性网络,便可复制出原模拟信号.因而, M有它自
51
己的特点,而且其编译码设备通常要比PCM的简单.
1. 增量调制原理增量调制原理增量调制原理增量调制原理
我们知道,一位二进制码只能代表两种状态,当然就不可能去表示抽样值的大小.可是,
用一位码却可以表示相邻抽样值的相对大小,而相邻抽样值的相对变化将能同样反映模拟信
号的变化规律.因此,由一位二进制码去表示模拟信号的可能性是存在的.为了确信这一点,
让我们通过下而的例子来说明.设一个频带有限的模拟信号如图5.3.5中的( )m t所示.现在把
横轴t分成许多相等的时间段t .此时可以看出,如果t 很小,则( )m t在间隔为t 的时刻上
得到的相邻值的差别(差值)也将很小.因此,如果把代表( )m t幅度的纵轴也分成许多相等的
小区间σ,那么,一个模拟信号( )m t就可用如图5.3.5所示的阶梯波形'( )m t来逼近.显然,
只要时间间隔t 和台阶σ都很小,则( )m t和'( )m t将会相当地接近.由于阶梯波形相邻间隔
上的幅度差不是σ+就是σ-.,因此,倘着用二进制码的"1"代表'( )m t在给定时刻上升一个
台阶σ,用"0"表示'( )m t下降一个台阶σ,则'( )m t就被一个二进码的序列所表征(见图5.3.5
中横轴下面的二进码序列).于是,该序列也相当于表征了( )m t.
图图图图5.3.5 增时调制波形示意图增时调制波形示意图增时调制波形示意图增时调制波形示意图
在讨论怎样得到发端的阶梯波形及由此波形又如何确定二进码序列之前,我们先讨论一
下在接收端怎样由二进码序列恢复出阶梯波形的问题.即 M信号的译码问题.不难看出,
接收端只要每收到一个"1"码就使输出上升一个值σ,每收到一个"0"码就使输出下降一
个σ值,连续收到"1"码(或"0"码)就使输出一直上升(成下降),这样就可以近似地复制出
阶梯波形'( )m t.这种功能的译码器可由一个积分器来完成,如图5.3.6所示.积分器遇到"1"
码(即有E+脉冲电压),就以固定斜率上升一个E ,并让E 等于σ,遇"0"码所表示的E-,
就以同样的斜率下降一个E .图5.3.6(b)表示了积分器的输入和输出波形.由此看到,积分
器的输出波形并不是阶梯波形,而是一个斜变波形.但因Eσ =,故在所有抽样时刻it上斜
52
变波形与阶梯波形有完全相同的值.因而,斜变波形同样与原来的模拟信号相似.由于积分
器实现起来容易,且能符合译码要求,故常被采用.最简单的积分器就是图5.3.6(c)所示的
RC积分器,其RC的乘积应远大于一个二进码的脉冲宽度.积分器输出虽已接近原来模拟信
号,但往往还包含有不必要的高次谐波分量,故需再经低通滤波器平滑,这样,便可得到十
分接近模拟信号的输出信号.
图图图图5.3.6 积分器译码示意积分器译码示意积分器译码示意积分器译码示意
现在回过来讨论, M的编码原理.一个简单的 M编码器组成如图5.3.7所示,它由相
减器,判决器,本地译码器及抽样脉冲产生器(脉冲源)组成.本地译码器与接收端的译码器
完全相同.判决器将在抽样脉冲到来时刻对输入信号的变化作出判决,并输出脉冲.这个编
码器的工作过程如下:将模拟信号( )m t与本地译码器输出的斜变波形'( )m t进行比较,为了获
得这个比较结果,先进行相减,然后在抽样脉冲作用下将相减结果进行极性判决.如果在给
定抽样时刻it有:
( ) | '( ) | 0
i it t t tm t m t
- -= =- >
则判决器输出"1"码;如有:
( ) | '( ) | 0
i it t t tm t m t
- -= =- <
则发"0"码.这里,it-是it时刻的前一瞬间,即相当于在阶梯波形跃变点的前一瞬间.于是,
编码器将输出一个如图5.3.5所示的二进码序列.
53
图图图图5.3.7 M编码器编码器编码器编码器
从上述讨论可以看出, M信号是按台阶σ来量化的(增,减一个σ值),因而同样存在
量化噪声问题. M系统中的量化噪声有两种形式;一种称为过载量化噪声,另一种称为一
般量化噪声,如图5.3.8所示.过载量化噪声(有时简称过载噪声)发生在模拟信号斜率陡变时,
由于台阶σ是固定的,而且每秒内台阶数也是确定的,因此,阶梯电压波形就跟不上信号的
变化,形成了很大失真的阶梯电压波形,这样的失真称为过载现象,也称过载噪声,如图5.3.8(b)
所示;如果无过载噪声发生,则模拟信号与阶梯波形之间的误差就是一般的量化噪声,如图
5.3.8(a)所示.图中的( ) ( ) '( )n t m t m t= -,可统称其为量化噪声.
图图图图5.3.8 两种形式的量化噪声两种形式的量化噪声两种形式的量化噪声两种形式的量化噪声....(a) 一般量化噪声一般量化噪声一般量化噪声一般量化噪声;;;;(b) 过载量化噪声过载量化噪声过载量化噪声过载量化噪声
设抽样时间间隔为t (抽样频率1/sf t= ),则一个台阶上的最大斜率K为:
sK f
t
σ
σ= =
这也就是译码器的最大跟踪斜率.当信号实际斜率超过这个最大跟踪斜率时,则将造成过载
噪声.因此,为了不发生过载现象,则必须使sf和σ的乘积达到一定的数值,以使信号实际
斜率不超过这个数慎.这个数值通常可以增大sf或σ来达到.
54
对于一般量化噪声,由图5.3.8(a)不难看出,σ小大则这个量化噪声大,σ小则噪声小.
用大的σ虽然能减小过载噪声,但却增大了一般量化噪声.因此,σ值应适当选取.
不过,我们看到, M系统的抽样频率必须选得足够高,因为这样,既能减小过载噪声,
又能降低一般量化噪声,从而使 M 系统的量化噪声减小到给定的容许数值.一般, M 系
统中的抽样频率要比PCM系统的抽样频率高得多(通常要高两倍以上).
2. M 系统中的量化噪声系统中的量化噪声系统中的量化噪声系统中的量化噪声
图5.3.9(a)是 M系统的组成方框图,它包括增量调制器,信道,检测器,积分器(译码
器)相低通滤波器等.为了简明起见,我们将图(a)中0( ), '( ), ( ), ( )qm t m t p t e t的波形画于图(b),并
设系统输出的信号和量化噪声分别用( )m t和( )qn t表示.如果信道的加性噪声足够小,以致不
造成误码,那么,用于接收的检测器将输出一个与0( )p t完全相同的信号,也即0 0' ( ) ( )p t p t=,
此时,系统的输出信号0( )m t与( )m t将有最好的近似(因为,量化噪声仍然存在);如果信道噪
声造成了误码,则在系统的输出噪声中不仅存在量化噪声,而且还存在由误码引起的噪声.
图图图图5.3.9 M系统的成及有关点的波形系统的成及有关点的波形系统的成及有关点的波形系统的成及有关点的波形
55
我们分析存在量化噪声时的系统性能.此时认为信道加性噪声很小,不造成误码.那么,
接收端检测器输出的0' ( )p t就是0( )p t,而解调积分器输出端的信号便是'( )m t(见图5.3.928).
容易看出,在这个积分器输出端的误差波形正是量化误差波形( )qe t.因此,如果求得( )qe t的
平均功率,则系统的输出量化噪声功率也就可以确定.我们还观察到,只要 M 系统不发生
过载现象(过载现象在设计时是需要克服的),那么( )qe t总是不大于σ±的.我们假设随时间随
机变化的( )qe t在区间( , )σ σ- +上均匀分布,于是( )qe t的一维概率密度( )qf e可表示为:
1
( )
2qf e
σ
= eσ σ- ≤ ≤ +
因而( )qe t的平均功率可表示成:
2
2 2 21
[ ( )] ( )
2 2q qE e t e f e de e de
σ σ
σ σ
σ
σ
+ +
- -
= = =∫ ∫ (5.3.9)
但要注意,上述的量化噪声功率并不是系统最终输出的量化噪声功率.这是因为,由图5.3.8
看出,( )qe t的最小周期大致是抽样频率sf的倒数,而且大于1/sf的任意周期都可能出现.因
此,从频谱的角度看,( )qe t的频谱将从很低频开始一直延伸到频率sf甚至更高.设( )qe t的功
率谱密度为( )eP f),则可以近似认为,
2
( )
3e
s
P f
f
σ
=, 0sf f< <
这就是说,( )qe t的平均功率被认为均匀地分布在频率范围(0,sf)之内.这样,具有功率谱密
度为( )eP f的噪声,通过低通滤波器(截止频率为mf)之后的量化噪声功率为:
2
( )
2
m
q e m
s
f
N P f f
f
σ
= =
(5.3.10)
由此可见, M系统输出的量化噪声功率与量化台阶σ及比值( / )m sf f有关,而与输入信号号
的幅度无关.当然,这后一条性质是在未过载的前提下才成立的.不发生过载现象,这实际
上是对输入信号的一个限制.
在临界条件最大的信噪比为:
2 2 3 3
0
2 2 2 2 2
3
0.04
8 8
s s s
q k k m k m
S f f f
N f f f f f
σ
π π
= = ≈i (5.3.11)
56
由此可见,最大信噪比0( / )qS S与抽样频率sf的三次方成正比.而与信号频率kf的二次
方成反比.因此,对于 M系统而言,提高抽样频率将能明显地提高信号与量化噪声的功率
比.
5.3.3 增量增量增量增量(差分差分差分差分)脉冲编码脉冲编码脉冲编码脉冲编码调制调制调制调制(DPCM)系统系统系统系统
前面对PCM和 M系统性能的分析可以看出,在不考虑信道误码率的情况下, M的
性能通常比PCM差.这主要是因为 M系统不管误差信号如何变化,传输的增量σ是固定
不变的.如果使增量的数值随误差信号( )tε的变化量化成M个电平之一,然后再进行编码,
这样,系统的性能就会得到改善.在这样的系统中,由于对传输的增量还要经过脉冲编码调
制,因面称它为增量脉冲编码调制或差分脉冲编码调制(DPCM).
图5.3.10示出了一增量脉冲编码调制系统的组成方框图.图中( )tε是输入信号( )m t和量
化信号'( )m t相减后的误差信号.该误差信号经抽样,量化与编码后,得到一二进制脉冲序列,
即DPCM信号.它一方面送往信道传输,另一方面经反馈支路的译码器和积分器,以便提供
所需的量化信号'( )m t.
图图图图5.3.10 增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图增量脉冲编码调制系统组成方框图
差分脉冲编码调制信号在接收端首先经译码器恢复成量化的误差信号,再经积分器和低
通滤波器便可恢复出原发送信号.
下面我们以两位编码(即2N=)来说明差分脉冲编码调制的编码过程.由于2N=,故误
差信号的量化电平数4M=(因2NM=).假设4个量化电平分别为3 , , , 3v v v v+ + - - ,它
们分别由二进制脉冲(双极性), , ,++ +- -+ - -表示.误差信号( )tε的抽样,量化和编码过程如图
57
5.3.11所示.图(a)表示对误差信号的抽样与量化,符号"0"表示( )tε的实际抽样值,"."表
示量化后的值.由图可见,当误差信号的取值在0 ~ 2v 范围时,则被量化为v+ ,即当
0 ( ) 2t vε≤ < 时,量化为v+ .同理,当2 ( ) 0tε- ≤ <时,量化为v- .2 ( )v tε ≤时,量化
为3v+ .( ) 2t vε=0 & abs(x(i))<32
67
out(i,2)=0; out(i,3)=0; out(i,4)=0; step=2; st=0;
elseif 32<=abs(x(i)) & abs(x(i))<64
out(i,2)=0; out(i,3)=0; out(i,4)=1; step=2; st=32;
elseif 64<=abs(x(i)) & abs(x(i))<128
out(i,2)=0; out(i,3)=1; out(i,4)=0; step=4; st=64;
elseif 128<=abs(x(i)) & abs(x(i))<256
out(i,2)=0; out(i,3)=1; out(i,4)=1; step=8; st=128;
elseif 256<=abs(x(i)) & abs(x(i))<512
out(i,2)=1; out(i,3)=0; out(i,4)=0; step=16; st=256;
elseif 512<=abs(x(i)) & abs(x(i))<1024
out(i,2)=1; out(i,3)=0; out(i,4)=1; step=32; st=512;
elseif 1024<=abs(x(i)) & abs(x(i))<2048
out(i,2)=1; out(i,3)=1; out(i,4)=0; step=64; st=1024;
elseif 2048<=abs(x(i)) & abs(x(i))=4096)
out(i,2:8)=[1 1 1 1 1 1 1];
else
tmp=floor((abs(x(i))-st)/step);
t=dec2bin(tmp,4)-48; % 函数dec2bin输出的是ASCII字符串,48对应0
out(i,5:8)=t(1:4);
end
end
out=reshape(out',1,8*n);
***********************************
其中解码函数为pcm_decode.m.
*******************************
68
function [out]=pcm_decode(in,v)
n=length(in);
in=reshape(in',8,n/8)';
slot(1)=0; slot(2)=32; slot(3)=64; slot(4)=128;
slot(5)=256; slot(6)=512; slot(7)=1024; slot(8)=2048;
step(1)=2; step(2)=2; step(3)=4; step(4)=8;
step(5)=16; step(6)=32; step(7)=64; step(8)=128;
for i=1:n/8
ss=2*in(i,1)-1;
tmp=in(i,2)*4+in(i,3)*2+in(i,4)+1;
st=slot(tmp);
dt=(in(i,5)*8+in(i,6)*4+in(i,7)*2+in(i,8))*step(tmp)+0.5*step(tmp);
out(i)=ss*(st+dt)/4096*v;
end
*******************************
运行结果如图5.3.16和图5.3.17所示.
012345678910
-1
-0.5
0
0.5
1
sample sequence
012345678910
-1
-0.5
0
0.5
1
pcm decode sequence
图图图图5.3.16 A律律律律PCM编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图编码后波形与输入波形的对比示意图
69
-60-50-40-30-20-100
0
10
20
30
40
50
60
输输窄输原原输输由dB值
窄量原量量
图图图图5.3.17 A律律律律PCM量化信噪比的动态范围量化信噪比的动态范围量化信噪比的动态范围量化信噪比的动态范围
利用上述程序,将时间t的间隔变为1/ 4096t =,则A律PCM编码后正弦信号与输入正
弦信号的波形对比如图5.3.18所示.
00.10.20.30.40.50.60.70.8
-1
-0.5
0
0.5
1
sample sequence
00.10.20.30.40.50.60.70.8
-1
-0.5
0
0.5
1
pcm decode sequence
图图图图5.3.18 A律律律律PCM编码后波形与输入波形的对编码后波形与输入波形的对编码后波形与输入波形的对编码后波形与输入波形的对比比比比示意图示意图示意图示意图
70
5.3.4 MATLAB提供的信源编提供的信源编提供的信源编提供的信源编/译码函数译码函数译码函数译码函数
在MATLAB通信工具箱中,提供了两种信源编/译码的方法,即标量量化和预测量化.
下面介始一下其使用方法.
1. 标量量化标量量化标量量化标量量化
标量量化就是给每一个落入某一特定范围的输入信号分配一个单独值的过程,并且落入
不同范围内的信号分配的值也各不相同.
(1) 信源编码中的μ律或A压扩计算函数compand.
格式一:out=compand(in,param,V)
格式二:out=compand(in,param,V,method)
功能:格式一实现μ律压扩;格式二实现μ律或A压扩.
参数含义:param对格式一为μ律;对格式二为μ或A值.V为峰值.压扩方法由methode
参指定,有:'mu/compressor', μ律压缩;'mu/expander',μ律扩展;'A/compressor',A律
压缩;'A/expander', A律扩展.
(2) 产生量化索引和量化输出函数quantiz
格式一:indx=quantize(sig, partition)
功能:根据判断向量partition,对输入信号sig产生量化索引indx,indx的长度与sig矢
量的长度相同.向量partition则是由若干个边界判断点且各边界点的大小严格按升序排列组
成的实矢量.若partition的矢量长度小于1N-,则索引向量indx中的每个元素的大小为
[0, 1]N-范围内的一个整数.量化方法如下:
若信号sig小于或等于partition(1),则输出0;若信号sig大于partition(1),而小于等于
partition(i+1),则输出1;若信号sig大于partition(N-1),则输出N-1.
格式二:[indx, quant]=quantiz(sig, partition, codebook)
功能:根据码本codebook,产生量化索引indx和信号的量化值quant.codebook存放indx
每个partition的量化值,对应indx=i-1的值在codebook(i).若partition的矢量长度为1N-,
则codebook长度为N.
格式三:[indx, quant, distor]=quantiz(sig, partition, codebook)
功能:产生量化索引indx,信号的量化值quant及量化误差disor.
(3) 采用训练序列和Lloyd算法优化标量算法的函数lloyds
71
格式:[partition, codebook]=Lloyds(training_set, ini_codebook)
功能:用训练集矢量training_set优化标量量化参八partition和码本codebook.Ini_codebook
是码本codebook的初始值.码本长度大于等于2,输出码本的长度与初始码本长度相同.输
出量化比参数partition的长度较码本长度小1.当ini_codebook为整数时,该函数以其作为码
本长茺.当处理后相对误差小于710-时,停止处理.
程序示列程序示列程序示列程序示列::::
用训练序列Lloyd算法,对一个正弦信号数据进行标量量化.
程序prog5301.m如下.
**********************************************
clear all
clear close
N=2^3;
t=[0:100]*pi/20; % 以3比特传输信道
u=cos(t);
[p,c]=lloyds(u,N); % 生成分界点矢量和编码手册
[index,quant,distor]=quantiz(u,p,c); % 量化信号
plot(t,u,t,quant,'*')
**********************************************
运行上述程序将显示出正弦信号的量化效果图.如图5.3.19所示.
0246810121416
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
图图图图5.3.19 标量量化前标量量化前标量量化前标量量化前,,,,后信号比较后信号比较后信号比较后信号比较
72
2. 预测量化预测量化预测量化预测量化
如果知道一些发送信号的先验信息,就可以利用这些信息,根据过去发送信号来估计下
一个将要发送的信号值,这样的过程即称之为预测量化.
差分脉冲编码调制(DPCM)就是采取了预测量化技术.MATLAB通信工具箱提供了
dpcmenco( ),dpcmdeco( )和decmopt( )函数.
(1) 差分脉冲调制编码函数decmenco
格式一:indx=dpemenco(sig, codebook, partition, predictor)
功能:返回DPCM编码的编码索引indx.其中参数sig输入信号,predictor为预测器传
递函数,其形成为1[0, , , ]mt t .预测误差的量化参数由partition和predictor指定.
格式二:[indx,quant]=dpemenco(sig, codebook, partition, predictor)
功能:除产生DPCM编码的编码索引indx外,还产生量化值quant.输入参数codebook,
partition, predictor可以由函数decmopt函数估计.当预测器为一阶传递函数时,为DPCM码
调制.
(2) 信源编码中的DPCM解码函数decmdeco
格式一:sig=dpemdeco(indx, codebook, predictor)
功能:根据DPCM信号索引indx进行解码.predictor为指定的预测器,codebook为码本.
格式二:[sig,quant]=dpemdeco(indx, codebook, predictor)
功能:根据DPCM信号索引indx进行解码,同时输出量化的预测误差quant.参数
codebook,predictor可以用dpcmoot函数估计.通常m阶预测器传递函数的形式为1[0, , , ]mt t .
(3) 用训练数据优化差分脉冲调制参数的函数dpcmopt
格式一:predictor=dpcmopt(training_set, ord)
功能:对给定训练集的预测器进行估计,训练集及其顺序由training_set和ord指定,预
测器由predictor输出.
格式一:[predictor, codebook, partition]=dpcmopt(training_set, ord, ini_codebook)
功能:输出预测器predictor,优选码本codebook,预测误差partition.输入变量ini_codebook
可以是码本矢量的初值或其长度.
程序示例程序示例程序示例程序示例::::
用训练数据优化DPCM方法,对一个余弦信号数据进行标量量化.
MATLAB程序prog5302.m如下.
73
***********************************************
clear all
clear close
N=2^3;
t=[0:100]*pi/20; % 以3比特传输信道
u=cos(t);
[predictor,codebook,partition]=dpcmopt(u,1,N); % 优化的预测传递函数
[index,quant]=dpcmenco(u,codebook,partition,predictor); % 使用DPCM编码
[sig,quant]=dpcmdeco(index,codebook,predictor); % 使用DPCM解码
plot(t,u,t,quant,'*')
***********************************************
程序运行结果如图5.3.20所示
0246810121416
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
图图图图5.3.20 DPCM预测量化误差预测量化误差预测量化误差预测量化误差
5.4 语音的线性预测编语音的线性预测编语音的线性预测编语音的线性预测编码码码码
线性预测分析方法让最有效的语音分析技术之一,其特点是既能极为精确地估计语音参
数,又有比较快的计算速度.因此,应用线性预测方法对语音编码已属高效率编码技术,它
不同于常规的数字编码.在上面讨论的各种量化编码方法,诸如PCM(脉码调制),M (增量
74
调制),DPCM(差值脉码调制)等都是常规的编码方法.因为这些方法部都是通用的波形数字
编码方法,并未考虑语音本身的特点,所以数据压缩率不可能很高.
线性预测的理论与方法在语音信号处理中的应用,主要表现在声码器技术领域中.除了
对语音信号的线性模型参数进行估计外,还能对所有的语音基本参数,诸如基音周期,共振
峰,语音谱,声强等进行估计.本节的日的不是讨论语音信号处理的所有问题,而是对语音
的线性预测编码(常用缩写字母LPC表示),主要是对语音模型参数进行估计.因此,首先要
介绍语音的数字模型;然后应用上节讨论的线性预测参数估计的方法来讨论语音模型的参数
估计问题,最后介绍一下基本语音参数的线性预测估计及相应的声码器原理.
5.4.1 语音信号的数学模型语音信号的数学模型语音信号的数学模型语音信号的数学模型
语音的生成机构大致可分为三部分.第一是声源声源声源声源.大家都知道人的声带是一种声源,它
会发生自激振动,产生常称为元音的语音,有时就叫作声音.另外还有两种声源:由声道变
窄(通常由舌头来控制)形成湍流叹声的称为摩擦声源.它产生的声音就叫摩擦音;由闭合的
声道急开而形成的脉冲波产生湍流噪声的叫爆破声源.第二是共鸣机构共鸣机构共鸣机构共鸣机构.它由鼻腔,口腔与
舌头组成,有时也叫声道声道声道声道.第三是放射机构放射机构放射机构放射机构.它由嘴唇或鼻孔发出声音并向空间传播出去.
对于这样的一种人类发声机能,可用多种的模型来模拟.我们这里特别关心的是所谓数字模数字模数字模数字模
型型型型,即能用数字处理方法来实现上述发声机能的一种模型.这种模型用准周期的脉冲源模拟
声带的振动,用随机嗓声(一种采用伪随机序列)模拟摩擦声,用可变参数的数学滤波器来模
拟声道谐振特性与放射持性.发声器官的数字模型如图5.4.1所示,图中开关是作浊音与清音
转换控制用的,乘法器是完成音强(或叫增益)控制功能的.数字滤波器的频率特性随受控参
数的变化而改变.如果所有的控制信号都由真实的语言信导分析所得,那么该滤波器(也叫声
道滤波器)的输出便接近于原始语言信号的信号序列,可以恢复出语音.这也就是声码器的基
本原理.
图图图图5.4.1 发声器官的数字模型发声器官的数字模型发声器官的数字模型发声器官的数字模型
75
现在用n来表示离散时间变量,( )s n表示n时刻语音信号的样值,并假设语音信号样值
序列是符号AR模型的,则有:
1
( ) ( ) ( )
p
i
i
s n a n i w n
=
= - +∑ (5.4.1)
其中,iα为模型参数,( )w n为白噪声序列.这样,语音的线性预测编码与AR模型参数估计
的关系可陈述如下:
若对输入的语音样值序列作出模型参数与残差估计,则n时刻的样值估计为:
1
^^( ) ( )
p
i
i
s n s n iα
=
= -∑ (5.4.2)
其中,^{ }iα为模型参数估计值的集合.估计误差或称残差信号便为:
1
^^( ) ( ) ( ) ( ) ( )
p
i
i
e n s n s n s n s n iα
=
= - = - -∑ (5.4.3)
若将( )e n作为图5.4.1声道滤波器的输入,滤波器能输出语言信号( )s n的近似值,则滤波器的
传递函数必须与n时刻的的输入,输出信号相适配.为此,对式(5.4.3)两边作Z变换,得:
1
^( ) ( ) ( )
p
i
i
i
E Z S Z Z S Zα-
=
= -∑ (5.4.4)
于是传递函为:
1
1
( )
^( ) 1
( )
p
i
i
i
S Z
H Z Z
E Z
α
-
-
=
= = -
∑ (5.4.5)
可知,声道滤波器是—个全极点型数字滤波器.该滤波器的参数{ | 1, 2, , }ii pα= 由语音分析
所得到的参数估计集合^{ }iα确定.于是对语音信导的编码转变为对参数集合^{ }iα及残差e的编
码问题.在接收端可由图5.4.1所示的数学模型恢复出语音信号.因此,线性预测编码的基本
问题就是很据输入的语音信号样值作出相应的AR模型参数估计.
5.4.2 由短时相关函数估由短时相关函数估由短时相关函数估由短时相关函数估计模型参数计模型参数计模型参数计模型参数
在5.2节里已作过介绍的AR模型参数的估计l方法,当然可以应用于语音信号A R模型
参数估计.不论是最小二乘法还是矩估计方法,估计的主要依据都是信息样值的自协方函数
76
或自相关函数.从自协方差函数或自相关函数的定义来看,应该从语音样值的整个集合来计
算概率平均值(参看式(5.2.3)与(5.2.4)).实际上这是做不到的.我们只能在时间轴上作短时平
均计算,以求得相应的短时间内发出的某个音的最佳预测参效,也就是动态最佳参数.如果
求整个长时间的平均,那么,对长时间的预测误差来说,系统可达到动态最佳状态.但对某
一短时间段信号来说,未必是最好的.也就是说,没有适应动态变化的系统,就达不到传输
动态信息的目的.在语音信号处理中,短时平均值的计算,通常采用两种方法,即协方差法
与自相关法.本书主要讨论两种方法.
1. 协协协协方差法方差法方差法方差法
在整个语音取样序列中,从n时刻开始取一段N个样值,记作( ), 0,1, , 1ns m m N= - .由
于对i时刻取样的样值常记为( )s i,因此该N个样值又可记作( 0), ( 1), , ( 1)s n s n s n N+ + + - .
两种记法之间的关系为:
( ) ( )ns m s n m= + 0,1, , 1m N= - (5.4.6)
于是,最二乘法方程式(5.2.10)改写为:
11
0 1 0
( ) ( ) ( 1) ( )
pN N
n n n n j
m j m
s m s m k s m s m kα
--
= = =
- = - -
∑ ∑∑ 1, 2, ,k p= (5.4.7)
它们的对应关系为:式(5.2.10)中的, ,i n t分别相当于式(5.4.7)中的, ,m N n;tix相当于( )ns m;
,t k ix-相当于( ) ( ) ( )n k ns m s n k m s m k-= - + = -.若令,
1
2
0
( ) ( , )
N
n n
m
s m k k k
-
=
- =∑ 1k p≤ < (5.4.8)
1
0
( ) ( ) ( , )
N
n n n
m
s m j s m k k j
-
=
- - =∑ 1k p≤ <, 0j p≤ ≤ (5.4.9)
则式(5.4.7)又可改写成:
1
2
(1,1) (1, 2) (1, ) (1,0)
(2,1) (2, 2) (2, ) (2,0)
( ,1) ( , 2) ( , ) ( ,0)
n n n n
n n n n
pn n n n
p
p
p p p p p
α
α
α
=
(5.4.10)
或简化为:
1
( , ) ( ,0)
p
j n n
j
k j kα
=
=∑ 1, 2, ,k p= (5.4.11)
77
式中的( , )nk j 就是信号的协方差函数.如果求出这些协方差函数值,那么最佳线性参数的估
值问题就是解方程组(5.4.11)的问题了.它们可由递推解法解出.协方差的计算是由计算机程
序来完成的.为了正确地编制程序,需对式(5.4.8)与(5.4.9)有确切的理解.为此举例说明如下:
设输入信号取样序列为(0), (1), (2),s s s ,计算5, 8, 4n N p= = =时( , )nk j 的值.下面我
们来写出式(5.4.10)左边矩阵的第一行与右边列矢量的元素,它们分别为:
以及
图5.4.2给出计算这8个协方差的各样值间的关系.对其他, ,n N p值的( , )nk j 的数值也
可作类似的计算.
图图图图5.4.2 求协方差的样值差求协方差的样值差求协方差的样值差求协方差的样值差,,,,5, 8, 4n N p= = =
符号"^"表示该两时刻样值相乘;"
"表示所有值相加.
78
2. 自相关法自相关法自相关法自相关法
另一种求短时平均的方法是用窗函数( )w m截取信号序列( )( 0,1, 2, )ns m m= 中N长的一
段,并在此范围内求平均,即:
( ) ( ) ( )ns m s n m w m= + (5.4.12)
其中,若:
1 0 1
( )
0
m N
w m
≤ ≤ -
=
其他
.则意味着除( ), ( 1), , ( 1)s n s n s n N+ + - 这N个样值外,
其余样值均为0.于是协方差函数表达式(5.4.8)与(5.4.9)的具体形式也有所不同.它们的主要
变化在于求和限要改变,比如当用0m=时.式(5.4.9)中将出现( )s n k-与( )s n j-,显然它们
均为0.为了清楚起见,作变量代换,令l m k= -,并考虑到( 1)s n N+ -右侧信号为0.故要
求1l k j N+ - ≤ -才行.于是式(5.4.9)变成:
1 ( )
0
( ) ( ) ( , )
N k j
n n n
l
s l k j s l k j
- - -
=
+ - =∑ 1, 2, , ; 0,1, ,k p j p= = (5.4.14)
式(5.4.14)与式(5.2.25)自协方差函数的表达式有相同的涵义,只是系数1/N及符号定义不同.
自然也可以求出自相关函数( )kρ的表达式.由于被窗口函数截取的样本在求平均过程中不改
变,故有( , ) ( , )n nk j j k =的特点.在语音信号处理中,通用的记号为:
( , ) ( , ) ( ) ( )n n n nk j j k R k j R τ= = - = (5.4.15)
式中:k jτ= -.这样,式(5.4.14)可改写成:
1
0
( ) ( ) ( )
N
n n n
l
R s l s l
τ
τ τ
- -
=
= +∑ (5.4.16)
与此相对应的AR模型参数估计方法是矩估计方法,要求解的方程组便是Yule-Walker方程.
应用式(5.4.16)的记号,方程式(5.2.23)变为:
1
(| |) ( )
p
k n n
j
R k j R kα
=
- =∑ 1, 2, ,k p= (5.4.17)
展开矩阵形式为:
1
2
(0) (1) ( 1) (1)
(1) (0) ( 2) (2)
( 1) ( 2) (0) ( )
n n n n
n n n n
pn n n n
R R R p R
R R R p R
R p R p R R p
α
α
α
-
- =
- -
(5.4.18)
79
左边的R矩阵为Toepliz矩阵,故式(5.4.18)的求解可以应用Durbin迭代程序或格型滤波算法.
如前所述的.
为了使在计算R矩阵中对各样本的应用关系更清晰,仍以8, 5, 4N n p= = =为例为说明.
由式(5.4.16)可写出,
图5.4.3是这个计算公式的示意图,很明显,它比图5.2.2简单,计算量也少得多,但它
的估计精度度不如协方差方法.还应该指出的是,由式(5.4.13)表示的窗函数,由于边缘的突
变会引起边缘样本系数估计的误差增大,有时候把它叫"边缘效应".为避免这个现象,往往
需对样本作预加重处理和采用和采用边缘缓变的窗函数.在线性预测声码器(LPC)中,需要考
虑上述的技术措施来改善合成语音的质量.这里不作详细讨论.
图图图图5.4.3 求自相关函数的样值间关系求自相关函数的样值间关系求自相关函数的样值间关系求自相关函数的样值间关系,,,,8, 5, 4N n p= = =
80
5.4.3 基音周期检测基音周期检测基音周期检测基音周期检测
利用线性预测方法提取话音参数而组成的声码器称为线性预测声码器.由于话音信号为
非平稳的随机过程,因此用时相关系数来过行语音参数的估计.除了模型参数的估计外,由
图5.4.3的模型可知,还需估计出基音周期及激励信号.此外,为了传输连续变化的语音信号,
先要将语音信号分帧,一般以10~30 ms数据作为一帧.在接收端两逐帧地进行合成并连接起
来组成连续话音输山.图5.4.4是线性预测声码器的原理框图.利用这样的声码器来传输话音
信号便可达到压缩数据率的目的.
图图图图5.4.4 线性预测声码器原理框图线性预测声码器原理框图线性预测声码器原理框图线性预测声码器原理框图....(a) 语音分析器语音分析器语音分析器语音分析器;;;;(b) 语音合成器语音合成器语音合成器语音合成器
语音分析器中参数估计方法在前几节中已作了详细的讨论,现在需要讨论基音周期检测
的方法.所谓基音周期检测或提取基音,就是要确定信号周期的长度.从理论上说这可以有
许多方法.但是由于语音信号的随机性,音调变化范围很宽,且语音还有嗓音或称浊音(Voiced)
与非嗓音或称清音(Unvoiced)之分,前者含有基音而后者却无基音,这就使基音提取变得复杂
且易出错.若基音提取有错误,会直接导致合成语音质量的下降,因此提出了各种改进方法.
这里只介绍1972年由Markel提出的称作简单逆滤波器跟踪(SIFT)算法,这是基于LPC分析
的一种方法,可以作为线性预测方法应用的一个例子.
图5.4.5是SIFT算法的原理框图.由于基音周期信息是低频信息,因此输人信号( )s n
81
首先经过截止频卒为900 Hz的低通滤波器.经滤波后的信号频带变为0~1kHz,故可以将原
来10kHz的取样频率经5次分频,降列2kHz(每5个样值去掉4个),得到离散信号( )x n.然
后从信号( )x n提取基音周期,其主要过程为逆滤波与求自相关函数峰值位置.
图图图图5.4.5 基音检测的基音检测的基音检测的基音检测的SIFT算法框图算法框图算法框图算法框图
所谓逆滤波,是对声道滤波的逆过程,即以经过低通滤波的语声信号( )x n作为输入,式
(5.4.6)的分母作为传递函数的一种滤波处理.该传递函数可表示为:
1
( ) 1
p
i
i
i
A Z Zα-
=
= -∑ (5.4.19)
其中滤波器参数iα由线性预测,即由式(5.4.10)或式(5.4.18)中的iα来估计.
信号( )x n经过逆滤波得到( )y n,实质上它就是预测器的预测误差.由线性预测理论可知,
( )y n有接近白噪声的属性,它的谱是接近平坦的,各样值间基本上互不相关.这就表明由式
(4.3.15)与(5.4.16)所定义的自相关函数一般很小,只有当0τ=或信号具有周期性(以P为周期,
, 2 ,P P P= ± ± )时,其自相关函数达到最大值.也就是说,可不用考虑信号的起始时间,只
要确定自相关函数中第一个最大值的τ值就可以代表其周期估计.图5.4.6表示这个分析过程
的几个典型波形.其中图(a)为一段被分析的输入波形,(b)为逆滤波器输出信号的时域波形;
(c)为逆滤波器输出信号的归一化自相关函数.由图可见,当8τ= ms时有明显的峰值,它便
可大致地代表该信号段的基音周期.
为了获得进行周期值的更高分辨率,可以对最大峰值所处的一定范围内的自相关函数进
行内插.当自相关函数的归一化峰值电平低于给定门限值时,则判定该论语音段为清音.
82
图图图图5.4.6 SIFT算法的典型信号波形算法的典型信号波形算法的典型信号波形算法的典型信号波形
SIFT算法利用线性预测分析提供了一个平坦谱的信号( )y n,以便于进行基音检测.谱的
平坦化越好(预测很确),对基音检测效果也起好,然而,对于高基音频率的声源(如儿童语音)
的谱不易平坦化.这是由于在0~900 Hz范围内缺乏一个以上的基音谐波(对由电话线输入的
信号尤其如此).对这类说话者及传输条件.就要采用别的基音检测方法了.
5.4.4 共振峰检测共振峰检测共振峰检测共振峰检测
在语音信号的数字模型一节中所叙述的语音生成机构的第二部分 —— 声道.一般可简
化成一个非均匀截面的声道.当声音从声源出发顺着声管传播时,它的频谱形状会被声管的
选择性所改变.这一效应正如同管乐器中所发生的谐振现象一样.声道管的谐振频率被称作
共振峰频率,共振峰频率与声道的形状和大小有关,每种形状都有一套共振峰频率作为其特
征,依次叫作第一共振峰频中,第二共振峰频率,……,并分别记为1 2, ,F F .人们借助
于舌,口腔等器官来改变声道形状,从而改变共振峰频率特征,使声音变成了可辨别的语言.
无疑,除了基音周期之外,共振峰频率是语音信号的另一个重要参数.图5.4.7显示了两个英
语元音的声波形图,除了有明显的基音周期外,其波形的结构有较大的差别.其中图(a)为元
音[i: ]的声波形,在一个周期内可看出有一个低频的衰减振荡,其上还叠加一个较强的高频率
83
振荡.可见该元音的第一共振峰频率较低,第二,第三共振峰频率较高.(b)为元音[u: ]的声
波形,可看出它的第一,第二共振峰频率很低,所以它的高频能量很小.
图图图图5.4.7 元音元音元音元音[i: ]与与与与[u: ]的声波形图的声波形图的声波形图的声波形图
以上讨论可知,若以共振峰参数来代表语音信息,同时配以基音周期参数,也可以很好
地合成语音信号.这样的语音分析与合成设备被称为共振峰声码器,图5.4.8是它的原理框图.
比较图5.4.4与图5.4.8可以发现,除了用共振峰分析与合成器代替预测参数分析与控制器之
外,其他都是类似的.我们的讨论侧重于如何应用线性预测方法来确定共振峰参数.确定共
振峰的方法很多,且各有优缺点.如求根法,即求出由式(5.4.19)所表示的传递函数的多项式
的根,也就是由式(5.4.5)所表示的声道滤波器的极点.
图图图图5.4.8 共振峰声码器原理框图共振峰声码器原理框图共振峰声码器原理框图共振峰声码器原理框图
84
对于大多数语来说,全极点模型很好地代表了声道的响应.但声学理论指出,鼻音和摩
擦音还要有谐振与反谐振(极点零点),这给数学处理带来了麻烦.但数学上很易证明当| | 1a<
时,有:
1
1
0
1n n
n
aZ a Z
-∞
- -
=
- =
∑ (5.4.20)
因而可以用多个极点去逼近一个零点,这就可以用全极点模型来讨论声道的谐振动特性.因
为式(5.4.4)分母的系数都是实数,所以分母多项式的根将是实根或复共轭根.将LPC数估计iα
代入式(5.4.5)的分母,得到多项式( )A Z为:
1
( ) 1
p
k
k
k
A Z Zα-
=
= -∑ (5.4.21)
若{ | 1, 2, , }iZ i p= 是( )A Z的一组根,则式(5.4.21)又可表示成:
1
1
( ) (1 )
p
k
k
A Z Z Z-
=
= -∏ (5.4.22)
由于Z平面上点与S平面上点的对应关系为:
kS T
kZ e= (5.4.23)
其中2k k kS j Fσ π= - +是与Z平面中的iZ相对应的S平面上的根,一般可表示成示成复共轭对
形式,即:
*, 2k k k kS S j Fσ π= - ± (5.4.24)
其中,"*"表示复共轭.kS和*
kS为声道的复谐振频率的一般表达式,各个中心频率为2kFπ,
它们可以直接由式(5.4.21)的根求出.共振峰的求根法的基本原理也基于此.
将式(5.4.24)代入式(5.4.23)并展开,得:
( 2 )*, cos(2 ) sin(2 )
Re Im
k k kj F T T
k kk k
k k
Z Z e e F T j F T
Z j Z
σ π σπ π- ± -= = ±
= ±
(5.4.25)
因此,如果已经解出LPC参数1 2, , ,pα α α ,那么通过解多项式(5.4.21)的根,便可以找出Z平
面上的p个根{ | 1, 2, , }iZ i p= .于是可求现相应的共振峰频率为:
85
Im1
arctan
2 Re
k
k
k
Z
F
T Zπ
=
(5.4.26)
我们定义频率响应大于中心频率的峰值响应的0.707倍的频率间隔为谐振带宽谐振带宽谐振带宽谐振带宽,则存中
心频率为2kFπ时,声道的谐振带宽近似为2kσ.在Z平面中.从原点到极点的距离就决定了
共振峰的带宽,即
| |kT
kZ eσ-= (5.4.27)
如果以Hz为单位计算,则带宽为:
21
ln | |
2
k
k kB Z
T
σ
π π
= = (Hz) (5.4.28)
由此可见,用线性预测方法作共振峰分析,经过对多项式求根,可较精确地确定共振峰
的中心频率与带宽,也可较容易地找出极点与共振峰的对应关系.这就为构成低比特率的共
振峰声码器提供了条件. '
5.5 图像的预测图像的预测图像的预测图像的预测编码编码编码编码
在图像压缩编码方法中,预测编码仍是主要手段之一,它常常与差值脉码调制(DPCM)
同义.在图像信源中.除了静态图像外,很大一部分为活动图像,其中电视网像又是最为常
见的一种活动图像.活动图像属于三维数据的信号源了.也就是说,活动图像是在时间轴上
一幅一随地更换出现的.比如电视图像,为了应用人眼的视觉残留效应而使之产生连续感,
则采用每秒25帧的更换速率,这里用'帧"代替通常叫"幅"的术语,由于帧速很高,这就
使各帧图象之间相差甚微,或者说存在着帧间相关性,于是,出现了帧内预测编码与帧间预
测编码或两者联合使用的编码方式.从而提高了图形数据的压缩倍数.
由于线性预测的基本理论及语言信号的线性预测编码都已作详细讨论,本节只介绍图像
预测编码的一些特殊点.首先介绍图像的数字化表示,然后介绍帧内DPCM乃帧间DPCM的
基本思想.
5.5.1 图像的数学化图像的数学化图像的数学化图像的数学化
所谓图像,就是一种视觉信息.它是用一个所谓二维光强度函数( , )f x y来表示的.这
里,x y表示空间坐标,而该点的函数值( , )f x y则表示试点的"光强度光强度光强度光强度",或通常称之为亮度或
86
灰度.因为光是能量的一种形式,故( , )f x y必为非零值,并且是在有限值城内的连续函数,
即:
0 ( , )f x y< < ∞ (5.5.1)
数字化图像则是存二维守间中对)t强度面故/(x,y)取样与量化的结果,以八x Y坍
短阵表示为:
(0,0) (0,1) (0, 1)
(1,0) (1,1) (1, 1)
( , )
( 1,0) ( 1,1) ( 1, 1)
f f f N
f f f N
f x y
f N f N f N N
-
- ≈
- - - -
(5.5.2)
式中左边( , )f x y是连续量,而右边矩阵使是通常所谓的数字图像数字图像数字图像数字图像.矩阵的每一个元素称作像
元或像素,常用缩写字母"Pel (Picture E1ements)"来表示.
现在,我们需要关心的两个参数是矩阵的大小N与每个像素的量化层数G.为了数字处
理的方便,一般都使它们由2的整数幂来表示,即有:
2nN= (5.5.3)
及
2mG= (5.5.3)
式中,,n m均为整数.m实际上就是每个像素的量化比特数.因此,一帧N N×的数字图像
的总比特数便为:
B N N m= × × (5.5.5)
例如,灰度为64级的128×128数字图像就有98304 bit数据量,需要12288 B = 12 KB的存储
单元.
那么,N与m应该多大才能使式(5.5.2)中的近似程度符合晴晰度要求呢 显然N与m越
大,清晰度就越好.但从处理,存储与传输的角度考虑,不希望N与m太大,以致于数据量
不会太大.根据国际标准CCIR 601建议的数字电视标准数据,高度信号的取样频率为13.5
MHz,两个色度信号取样频率各为6.75MHz,8m= bit,于是,数字电视的比特率高达
216Mb/s,高清晰度电视的比特率更达1.188Gb/s.处理一般图像所用的最小系统,也应能显
示256×256个像素与64灰皮级,即N=256与m=6才行.即使如此,一帧图像的数据量仍
然十分大,故数字图像的数据压缩要求比语声信号的更为迫切.
87
5.5.2 帧内预测编码帧内预测编码帧内预测编码帧内预测编码
帧内预测编码利用的是一帧图像中各像素之间的相关性.正如式(5.5.2)所表示的数字图
像那样,它由一行行像素组成,一行中各像素之间又存在相关性,行与行之间也存在相关性.
如果只利用一行内像素间的相关性进行预测,则叫一维预测编码.如果同时利用行内与行间
相关性进行预测,则为二维预测编码.
1. 一维预测一维预测一维预测一维预测
一维预测的基本原理与线性预测的基本原理是相同的.但在图像预测编码中,往往不去
计算相关系数.而是由实验观察确定.这不仅是为计算方便,还因为按均方误差准则作出的
最佳预测,并不反映人眼的主观感觉.一维预测在传真图像,即二维图像编码中应用较多.
图图图图5.5.1 一维预测过程一维预测过程一维预测过程一维预测过程....(a) 扫描行扫描行扫描行扫描行;;;;(b) 预测器预测器预测器预测器;;;;(c) 逆预测器逆预测器逆预测器逆预测器....
给扫描行上的像素编上如图5.5.1(a)所示的号码.预测算法就是用已经扫描过的像素序列
1 2, ,i ix x- - 作为参考像素,按某种数学关系1 2( , , )i if x x- - 对现时像素ix的值作出预测,得到预
测值^
ix.然后将^
ix与ix作比较,若一致记作"0",不一致则记作"1".在二值图像酌情况下,
可以用模2加来完成这个操作.在多灰度值图像时,失得到ix与^
ix的差值,然后对差值作出
88
量化编码,这就是已介绍过的DPCM过程.这个关系的数字表达式为:
1 2
^( , , )
^( )
^( )
i i i
i i i
i i i
x f x x
y x x
y x x
- -=
= ⊕
= -
二值图像
多值图像
(5.5.6)
符号⊕表示模2加.
上面已指出,( )fi表达式由实验观察确定,一般仍为线性多项式,即,
1
^
l
i k i k
k
x a x-
=
=∑ (5.5.7)
其中系数ia就是根据实验确定的参数.若1l=,实际上就是前一像素预测.在二值传真图像
情况,由于扫描结果总是连续"0"与连续"1"相间,因此,误差图像代表了图像轮廓,如图5.5.2
所示.
图图图图5.5.2 前像素预测的误差图样前像素预测的误差图样前像素预测的误差图样前像素预测的误差图样....(a) 两相邻扫描行两相邻扫描行两相邻扫描行两相邻扫描行;;;;(b) 前像素预测误差图样前像素预测误差图样前像素预测误差图样前像素预测误差图样
2. 二维预测二维预测二维预测二维预测
虽然图像编码的二维DPCM可以由一些经验给出预测系数.但严格地说,图像编码的二
维DPCM应根据不同的图像所求得的自相关函数来求出最佳预测系数,然后设计出与图像相
匹配的预测编码系统.但这样做计算量与复杂性都将增加,人们经过大量的实验企图寻找一
种通用的算法,常见的有两种表达式比较接近实际:
一种模型是:
1 2| | | |( , )c cR eα βα β- -= (5.5.8)
其中,α和β分别代表水平与垂直方向的像素间距离像素间距离像素间距离像素间距离;1 2,c c则可根据不同的图像进行调整选
89
取的系数.
另一种模型是:
2 / 33 / 2 3 / 2( , ) ( )Rα β α γβ = + (5.5.9)
其中,α和β仍分别代表水平与垂直方向的像素间距离,γ则代表行间距离与水平方向像素
之间的距离之比值.
关于用作预测的像素选择,一般来说,1 2 3 4, , ,s s s s与0s相关性最强,因此用3到 4个预测
像素已有较高的准确性.再用更多的像素来预测,所得到的准确性改善就不很显著了.为了
使上述讨论有一个数量的概念,下面给山一组协方差函数与预测系数,预测误差的实验数据
(取自J.B. O'Neal于l966年发表的结果).预测像素分布结构如图5.5.3所示,该数据是对一
帧侧面头像取样的结果.
图图图图5.5.3 二维预测的邻近像素排列二维预测的邻近像素排列二维预测的邻近像素排列二维预测的邻近像素排列
协方差:
1 01 02 013
04 05 06 01
1, 0.934 0.960 0.905
0.919 0.832 0.814 0.829
R R R R
R R R R
= = = =
= = = =
下标"0i"表示表示0s与is之间的相关系数.
预测系数与预测误差如图5.5.1所示.
图图图图5.5.1 预测系数与预测误差预测系数与预测误差预测系数与预测误差预测系数与预测误差
90
还应该指山的是,以上讨论都是对单色图像或黑白电视图像进行的.至于彩色电视图像,
对于目前流行的体制,不管是PAL制还是NTSC制,色度信号都是由色副载波加到亮度信号
上去的.这就使全电视信号(或称复合电视信号)的样值之间的相关性,由于色度皮信号与亮
度信号的相位不一致而发生了变化.这就要重新设计它们的预测结构与预测系数,最直接的
办法就是将亮度信号与色度信号分开来处理,叫分量预测编码.由于它已用于专门技术,故
本课程不作详细讨论.但它所应用的基本原理还是本章叙述过的内容.
5.5.3 帧内预测编码帧内预测编码帧内预测编码帧内预测编码
帧间预测编码是对活动图像,特别是电视图像实施的一种压缩编码方法.它是利用活动
图像相邻帧间的相关性,即时间相关性来达到压缩数据的目的.显然,为了比较与处理,需
要能存储二帧以上图像的帧存储器,技术比较复杂,成本也较高,但效果较好.将电视信号
压缩至16Mbit/s的数码率,仍可保持较好的图像质量,因而有很大的吸引力.本节将介绍
帧间预测编码的基本思想及一些基本实现方法.
1. 帧间帧间帧间帧间预测的基本思想预测的基本思想预测的基本思想预测的基本思想
以电视信号为例来说明帧间预测的基本思想.现行的电视信号是由奇数场(由奇数行组成)
与偶数场(由偶数行组成)构成一帧.因此.帧间预测也包含有场间预测的含义.图5.5.4画出
了相邻三场及相邻像素的示意图,每个像素所标的坐标依次为(列,行,场), 当前场中待
预测的像素坐标为(0,0,0),坐标正方向如图所示.
图图图图5.5.4 相邻三场的相邻像素相邻三场的相邻像素相邻三场的相邻像素相邻三场的相邻像素
91
设I为所有用于预测的像素坐标的集合,三维矢量0=(0,0,0)为待预测像素坐标,则
0I.
对某帧图像来说,任何一个像素sn都可能是待预测像素,设{ | }s-∈n ii I是与sn 相邻的已
传输像素的集合,则sn的预测值可表示为:
^( )s a i s-
∈
=∑n n n i
i I
(5.5.10)
这里,i n都应为三维矢量.
类似于一维预测的讨论,关键是要选择( )ani使均方误差值2^( )E s s-n n最小,即设计出在
均方误差测量条件下的最佳预测器.设图像信号为广义平稳的随机过程,则预测器系数可由
以下方程的解得到:
( ) ( ) ( )s sa i R R
∈
- = ∈∑n
i I
i k k k I (5.5.11)
由此组系数组成的预测器,得到均方误差值为:
( ) ( ) ( )s sR a i R
∈
-∑0n
i I
i (5.5.12)
因为2( )s sRσ=0是信号的方差,故( ) ( )sa i R
∈
∑n
i I
i就表示预测使传输信号方差减小的量,即代表
了预测增益的大小.在极端情况下,若预测误差为0,则信号就不需传输,这表示一切都在"预
料之中",显然这对信息传输来说是无意义的.
式(5.5.10)与(5.5.11)包含了帧内预测(只用帧内像素预测),帧间预测(只用不同场或帧的像
素来预测)以及帧内帧间混合预测(同时用三维像素来预测).Haskell等人对不同预测器作了实
验研究,如表5.5.2所示.
表表表表5.5.2 不同的预测器的比较不同的预测器的比较不同的预测器的比较不同的预测器的比较
注:坐标参看图5.5.4
92
实验表明,一般来说,帧(场)间预测可得10 dB左右的增益,在运动不大的场面,帧间预
测性能较好,而运动较大的场面,还是帧内预测为好.在混合预测时,以场和豫素差的预测
为佳.
作为帧间预测的特例,对活动很慢的图像,可以用前一帧像素来代替(不是预测,即不计
算差值)当前帧,因而可作隔帧传输.对未传输的帧,可用简单重复的方法补上,或作插值补
上.由于人眼视觉特性,对静止或慢动图像要求较高的空间分辨率,而对时间分辨率要求饺
低.故上述做法带来的误差往往不易被觉察.
2. 运动补偿预测运动补偿预测运动补偿预测运动补偿预测
上面谈到人眼对画像的静止部分要求较高的空间分辨率与较低的时间分辨率,而对运动
部分则恰好相反.于是可以利用人眼的视觉特性,将图像分成静止部分与运动部分分别处理.
对静止部分只用前一帧数据的重复,对运动部分则设法测定其位移方向与位移量.以位移量
确定预测像素与预测误差,并对只作重复的静止部分作出补偿,构成当前帧的完整图像,故
称之为运动补偿预测.
运动补偿预测技术由以下几方面组成:
(1) 分割图像为静止的与运动的两部分;
(2) 估计运动部分的位移;
(3) 用位移估计进行运动补偿预测;
(4) 预测信息编码.
运动补偿预测的原理框图如图5.5.5所示,输入信号送至预测器,同时也加到分割器和位
移估计器.分割信息还送至位移估计器,二者又同时用控制预测器,并作编码传送.
图图图图5.5.5 运动补偿预测编码器原理运动补偿预测编码器原理运动补偿预测编码器原理运动补偿预测编码器原理
93
由上所述,可见图像分割是整个工作的基础,也是较为困难的一步.这里介始两种简单
的方法.一种称为位移估计法.即把图像分成若干矩形子块,适当选择块的大小,并确定子
块是动的或是不动的,估计出运动于块的位移,以便进行预测,另一种叫递归估计法,它是
对每个像素都进行递归估计,作出最佳逼近,确定位移.
(1) 位移估计匹配算法
匹配算法是在当前帧,比如说第k帧,取一矩形像素块(M N×个像素),与前一帧(第1k-
帧)相同位置的下一个M N×像素块作比较.若灰度完全一样,或变化最小,则该块为不动块.
否则要在前一帧(第1k-帧)图像中扩大范围再作比较.比如扩大至( 2 ) ( 2 )m mM d N d+ × +,其
中md表示在帧间隔时间内可能运动的最大位移,这要根据图物体运动速率来估算.若,取
10md=,则扩展区如图5.5.6所示.
图图图图5.5.6 扩展的匹配区域扩展的匹配区域扩展的匹配区域扩展的匹配区域
匹配程度的计算,可以用均方误差最小或绝对误差最小等准则来作比较.理论上应该在
匹配区内的所有可能取的M N×子块都与第k帧内的M N×子块亮度值作比较,作出差别最
小的判决.
(2) 递归估计算法
近年来对像素递归估计算法研究较多,出现了许多算法,比如Newton Raphson算法,
Bergmann算法,Robbins梯度算法等.这里不再较详细讨论.
应用帧间预测技术可获得高的数据压缩倍数,对窄带传输系统特别有吸引力,因而在会
议电视及可视电话等数字传输中很早就得到了的应用.20世纪未,由于微电子技术和计算机
94
技术的高速发展,使运动估计和运动补偿技术有了快速实现的可能,出现了许多实时,有效
的算法,从而使帧间预测技术在数字电视与高清晰度电视信号的数据压缩编码中得到了广泛
的应用.
·上一篇：*SoCG~A'AuJP'A
·下一篇：多媒体技术模拟试卷