第五章 列文森递归算法
§5.1 引言
第四章中已看到,很多种信号建模问题都需要求解如下形式的线性方程组:
x p =R a b (5.1)
其中x R 是Toeplitz 矩阵。例如在Pade 逼近法中,估计由矢量p a 所表示的分母系数()p a k 就是求解象式(5.1)的一组Toeplitz 方程,其中x R 是一个非对称的Toeplitz 矩阵,其第一列是信号值
(),(1),...,(1)x q x q x q p ++-,第一行是(),(1),...,(1)x q x q x q p --+,另外矢量b 由信号值(1),(2),...,(
x q x q x q p +++组成,因此与矩阵x R 的元素值密切相关。ARMA 过程建模中的修正
Y ule-Walker 方程也是这一类的线性方程组。在用Prony 法或自相关法对确定性信号进行全极点建模,以及用Y ule-Walker 法对随机过程的全极点建模中,都要遇到Toeplitz 方程。但与Pade 逼近时不同,这时的x R 是自相关值(0),(1),...,(1x x x r r r p -的哈密顿Toeplitz 矩阵。另外,由于这时的
[](1),,(1)T
x x r r p =-
b ,式(5.1)右边矢量还是与Toeplitz 矩阵的元素值密切相关。在求解分
子系数的Shanks 法中,又会遇到一组哈密顿Toeplitz 方程,但与前面不同的是这时的矢量b 中不含有矩阵x R 的元素值。在第七章的FIR 维纳滤波器设计中,还要遇到Toeplitz 方程,与Shanks 法中一样,其x R 是哈密顿Toeplitz 的,但b 的取值与x R 的值无关。
由于在很多问题中求解Toeplitz 方程的重要性,本章介绍求解这类方程的有效算法。在导出这些算法的过程中,我们将发现这类方程的解的一系列有趣的性质,由此可进一步获得信号建模的其它方法。在§5.2节,我们首先导出Levinson-Durbin 递归算法,该算法可用于求解Prony 全极点正则方程和自相关正则方程。Levinson-Durbin 递归将导致几个重要的结果,包括格型滤波器结构、数字滤波器的Schur-Cohn 稳定性测试、Toeplitz 矩阵的Cholesky 分解,以及Toeplitz 矩阵的递归求逆。§5.3节介绍求解一般的哈密顿Toeplitz 方程的Levinson 递归方法,注意这时对矢量b 没有约束。该算法可用在Shanks 法中,也可用于求解第七章中的一般FIR 维纳滤波问题。§5.4节将导出分基Levinson 递归算法,它比Levinson-Durbin 递归的效率要稍高些,并引入了奇异预测器多项式和线谱对的思想,它们在语音处理应用中很有用。
§5.2 Levinson-Durbin 递归
1947年Levinson 给出了求解一般的线性对称T oeplitz 方程组x =R a b 的递归算法。在一个有关维纳线性预测问题的评注论文中,Levinson 称该算法是一个“数学上平凡的过程”[16],但是,这一递归算法导致了若干重要的发现,其中包括格型滤波器结构,它在语音处理、谱估计、数字滤波器实现中获得了广泛的应用。后来在1961年,Durbin 针对方程右边是单位矢量的特例改进了Levinson 递归算法[7],本节我们就介绍该算法,称为Levinson-Durbin 递归。另外,将给出该递归的一些性质,说明它如何导致格型滤波器结构,并证明由自相关法所导出的全极点模型是稳定的。
5.2.1 递归式的推导
用Prony 法或自相关法进行全极点建模时,要求解正则方程,对p 阶模型,方程为:
1
()()()0p
x p
x l r k a
l r k l =+
-=∑; 1,2,...,k p = (5.2)
而建模误差为:
1
(0)()()p
p x p
x l r a
l r l ε==+
∑ (5.3)
将上面两个方程组合成矩阵,应有:
***
***1(0)(1)(2)()1(1)0(1)(0)(1)(1)(2)0(2)(1)(0)(2)0()()(1)(2)(0)x x x x p x
x x x p p
x x x x p x x x x r r r r p a r r r r p a r r r r p a p r p r p r p r ε??????
??????-????=??-????????????--??????
(5.4) 它有1p +个线性方程,也有1p +个未知量,即(1),(2),...,()p p p a a a p 和p ε。上式的等效矩阵形
式为:
p p ε=R a u
(5.5)
其中p R 是(1)(1)p p +?+的哈密顿Toeplitz 矩阵,[]1,0, 0
=u 是单位矢量。对实值数据,x R 将是对称Toeplitz 矩阵。
求解方程(5.5)的Levinson-Durbin 递归是一个按模型阶递归的算法,换言之,(1)j +阶全极点模型的系数1j +a 是由第j 阶模型的系数j a 而获得。因此,我们首先说明如何由j 阶正则方程的解导出
1j +阶正则方程的解。设()j a i 是如下的正则方程的解:
*
*
*
*
*
*1(0)
(1)(2)()(1)(1)(0)(1)(1)0(2)0(2)(1)(0)(2)()0()
(1)(2)
(0)x x x x j j x x x x j x
x x x j x x x x r r r r j a r r r r j a r r r r j a j r j r j r j r ε????????????-?
?
????=-???????
?????--??
????
(5.6)
方程的矩阵形式为:
1j j j ε=R a u
(5.7) 给定j a ,我们要导出如下的(1)j +阶正则方程的解:
111
1j j+j ε
++=R a u
(5.8)
推导过程如下。假定我们将j a 矢量补一个零,并用1j +R 左乘新矢量,则结果应为:
*
*
*
*
*
*
*
*
**
1(0)
(1)
(2)()(1)(1)(1)
(0)(1)(1)()(2)(2)(1)(0)(2)
(1)()()(1)(2)(0)(1)0(1)
()
(1)
(1)
(0)x x x x x j x
x x x x j x
x x x x j x x x x x x x x x x r r r r j r j a r r r r j r j a r r r r j r j a j r j r j r j r r r j r j r j r r ??+????-?
?
?--???
?
?--??+-?????
000j j εγ?????????=?????????
??
??? (5.9)
其中参数j γ为:
(5.10)
注意若0j
γ
=,则(5.9)式右边是一个(1)j +维单位矢量,且1=1, (1), ..., (), 0T
j j j a a j +????
a 是(j +1)阶正则方程(5.8)的解。但通常0j
γ
≠,1,(1),(2),...,(),0T
j j j a a a j ????
不是方程(5.8)的解。 推导Levinson-Durbin 递归的关键步骤是利用1j +R 的哈密顿Toeplitz 特性,使得可将式(5.9)写成如下的等效形式:
***
*
*
*
**
*
*
(0)(1)
(2)()(1)0()(1)(0)(1)(1)()(1)(2)(1)(0)(2)
(1)(1)()(1)(2)(0)(1)
1(1)()
(1)
(1)
(0)x x x x x j x x x x x j x x x x x j x
x x x x x x x x x r r r r j r j a j r r r r j r j a j r r r r j r j a r j r j r j r r r j r j r j r r +??????-??---????--??+-????
?
000j j γε??????????
=?????
???????
??
? (5.11)
取式(5.11)的复数共轭并将所得方程与式(5.9)进行加权组合,则对任意的(复)常数1j +Γ,有: ***111**01()(1)00(2)(1)00............00()()01j j j j j j j j j j j j j a j a a a j a j a j εγγε+++??????????????????????????????
??-+Γ=+Γ??????????????????????????????????????????????
R (5.12) 由于我们是要求解矢量1j+a ,使得它被
R 相乘后得到一个(1)j +维的单位矢量,显然,若取:
(5.13)
(1)j +(5.12)就变为:
111
1j j+j ε
++=R a u
其中
*
*1
1
*01()(1)(2)(1)()(1)01j j j j j+j j j a j a a a j a j a +??
??????????-=+Γ
??
??????????
????
??
a (5.14)
它就是(1)j +阶正则方程的解。另外,
(5.15)
是(1)j +阶建模的误差。(注:由于j ε是实数,式(5.13)中的复数共轭可以去掉。)如果我们定义
(0)1j a =,(1)0j a j +=,则被称为阶递归方程的式(5.14)可以简单表达为:
(5.16)
0j =阶模型的
解给出,即为:
0(0)1a =
0(0)x r ε= (5.17)
综上所述,可概括Levinson-Durbin 递归算法的各步骤如下。首先用零阶模型的解(5.17)对其初始化,然后对0,1,...,1j p =-,由第j 阶模型的解分三步求第(1)j +阶模型。第一步是用式(5.10)和(5.13)确定1j +Γ的值,1j +Γ也称为第(1)j +个反射系数;第二步是利用Levinson 阶更新方程由
()j a i 计算系数1()j a i +;最后一步是用式(5.15)更新误差1
j ε
+。该误差还有另外两个等效表达式,在
后面的讨论中要用到。第一个是:
1
2
2
111
1(0)1j j j j x i i r εε+++=????=-Γ=-Γ????
??
∏
(5.18)
第二个形式是来自式(5.3),表达为:
1
11
1
(0)()()j j x j x i r a
i r i ε+++==+
∑
(5.19)
完整的递归算法总结于表5.1。图5.1给出了其Matlab 程序。(注:不同的作者对Levinson-Durbin 递归中的变量定义有所不同,例如(5.12)中也可用1
j +-Γ或*
1j +Γ,结果就导致1j +Γ的符号不同。)
图5.1 列文森-杜宾递归的Matlab 程序
_____________________________________________________________________________________ 例5.2.1 求解自相关正则方程。
假定要用Levinson-Durbin 递归求解自相关正则方程,以确定一个信号的3阶全极点模型。设已知信号的自相关值为:
(0)1x r =,(1)0.5x r =,(2)0.5x r =,(3)0.25x r = 三阶全极点模型的正则方程为:
333(1)10.50.50.50.5
10.5(2)0.50.50.510.25(3)a a a ??????=-??????????????????
利用Levinson-Durbin 递归,我们得: 1. 一阶模型:
02
01110111(1)
(1)13, (0)1(0)24
11101 012x x x x r r r r γγεε=??Γ=-=-=-=-Γ=????
????????
=+Γ==Γ??????-????????
??
a 且 2. 二阶模型:
2
11122121
211
2(2)(1)(1), , 14
33
101-1/21/31/21/3011/3x x r a r γγεεε??=+=
Γ=-
=-
=-Γ=????????
--=-??????
??????-??????
a =且
3. 三阶模型:
2222
322323(3)(1)(2)(2)(1)1/12,
/1/8, 121/32x x x r a r a r γγεεε=++=-??Γ=-==-Γ=??
且 31
01
11/
31/33/81/31/33/880
11/8
???
?
??
---?
???
??=+
=??????---?????????
?
??
a
最后,取(0)8b =
=
,得全极点模型为:
31
2
3
8()33
1
1888
H z z
z
z
---=
-
-
+
获得了(3)x 的三阶全极点模型后,若3()H z 是()x n 的正确模型,即()x n 是3()H z 的逆Z 变换,
则让我们来确定()x n 的后面的自相关序列应是什么。这时,模型对3k >将保持不便,即对3k >,
0k Γ=,因此0k γ=。这样,由(5.10),并注意3j =,有:
3
3
1(4)()(4)0x x i r a
i r i =+
-=∑
求得(4)x r 为:
3
31
7(4)()(4)32
x x i r a i r i ==--=
∑
相继的()x r k 值可用类似的方式确定。例如,对4k ≥,
3
31
()()()x x i r k a i r k i ==--∑
它就是AR(3)过程的Yule-Walker 方程。
???????????????????????????????????????????????????????????????????????????????????????????????????????
由第一章§1.5.2节,若方差为2
w σ的白噪声激励一个一阶全极点滤波器:
1
1()1H z z
α-=
-
将输出一个AR(1)过程,且其自相关序列为:
2
2
()1k
w
x r k σα
α
=
-
在下面的例子中,将用Levision-Durbin 递归求有关该过程的反射系数k Γ和模型误差k ε的序列。
_____________________________________________________________________________________ 例5.2.2:AR(1)过程的反射系数。
设()x n 是一个AR(1)过程,其自相关为: 2
221,,, (1)
p
w
x σαααα??=
?
?-r 为求出有关该自相关矢量的反射系数,我们可以利用Levinson-Durbin 递归。首先初始化:
01=a ,2
02
(0)1w
x r σεα==
-
因此第一个反射系数为: 1
(1)(0)
x x r r αΓ=-=-
且 221011w εεσ??=-Γ=??
因此一阶模型为:
11α?
?
=-???
?a 对二阶模型,
11(2)(1)(1)0x x r a r γ=+=
且
20Γ=, 2
21w εεσ==
因此二阶模型为:
210α??=-??????
a
继续上述过程, 可导出, 对所有2
1,0,j j w j εσ>Γ==。具体说,假定第j 步的递归中0j Γ=,且
[]1
...
0T
j α=-a
则由于1
(1)()(1)(1)()
j
j
x j x x
x
i r j a i r j i
r j r j γα==++
-+=+-∑
,再看自相关序列的形式,有0j
γ
=,因此1
0j +Γ
=,且1j +a 是在j a 中加一个零而形成。
概括起来,对AR(1)过程,有:
[][]
2
222
2
0 (01)
...
0...1T
p T
p T
w p w w w αασσσσα=-=-??=??-??
a Γε
??????????????????????????????????????????????????????????????????????????????????????????????????????? 下面来看Levinson-Durbin 递归的计算复杂度,将它与求解p 阶自相关正则方程的高斯消去法进
行比较。用高斯消去法求解含p 个未知量的线性方程组约需
3
13
p 次乘法和除法,而Levinson-Durbin
递归在第j 步是需22j +次乘法,1次除法和21j +次加法,整个递归需p 步,总的乘法和除法次数是:
1
2
(23)2p j j p p -=+=+∑
总的加法次数是:
1
2
(21)p j j p -=+=
∑
因此,Levinson-Durbin 递归中的乘法和除法次数与2
p 成正比,而高斯消去法是与3
p 成正比。
该递归法优于高斯消去法的另一点是它所需的数据存储量较小。具体说,高斯消去法需2
p 的存储量,而Levinson-Durbin 递归只需21p +:1p +个用于自相关序列(0),...,()x x r r p ,p 个用于模型参数
(1),...,()p p a a p ,一个用于误差p ε。
尽管Levinson-Durbin 递归与高斯消去法相比,效率提高了,但应指出的是,求解正则方程可能只是整个建模过程计算量的一小部分。例如,对长度N 的信号,计算自相关值(0),...,()x x r r p 约需
(1)N p +次乘法和加法,因此,若N p >>,则求自相关序列的计算量将是建模算法所需计算量的
主要部分。
§5.2.2 格型滤波器
Levinson-Durbin 递归的一个副产品是数字滤波器的格型结构。它是数字滤波器的一个常规实现
方法,具有很多有趣和重要的特性,如模块结构、对参数量化效应不敏感、是保证滤波器稳定的简单方法等。本节将说明Levinson 阶递归方程如何可以被用于导出格型FIR 滤波器结构。第六章将介绍其它的格型滤波器结构,包括全极点格型和极零点格型滤波器,以及这些滤波器如何用于信号建模。 为推导格型滤波器的结构,应考察式(5.16)的Levinson 阶递归方程。但为了便于后续推导,首先定义倒序(reciprocal)矢量R
j a ,它是将矢量j a 的各元素倒序并取复数共轭而形成,即:
****
()1(1)(1)(2)(2) ......(1)(1)()1j j j R j j j j j j j a j a a j a a j a j a a j ????
????
-????-??=?=??????-????
????????
a a 或:
(5.20)
在Levinson 11()()(1)R
j j j j a i a i a i ++=+Γ-
(5.21)
用()j A z 表示()j a i 的Z 变换,()R
j A z 为倒序序列()R
j a i 的Z 变换,则由(5.20)得()j A z 和()R
j A z 的关系:
(5.22) 根据()j A z 和()R
j A z 重写式(5.21),将给出:
1
11()()()R
j j j j A z A z z A z -++=+Γ (5.23)
它是有关()j A z 的阶递归方程。
下一步推导1R
j a +和1()R
j A z +的阶递归方程。注意式(5.21),若两边取复数共轭,并将i 替换为
1j i -+,则有:
***
11(1)(1)()j j j j a j i a j i a i ++-+=-++Γ
(5.24)
再利用倒序矢量的定义式(5.20),可得所需的1()R
j a i +的更新方程:
*
11()(1)()R
R
j j j j a i a i a i ++=-+Γ
(5.25)
变换到Z 域,就得到有关1()R
j A z +的阶更新方程:
1
*
11()()()R
R
j j j j A z z A z A z -++=+Γ
(5.26)
总结起来,我们有一对耦合的差分方程:
1
1*
11()()(1)()(1)()
R
j j j j R
R j j j j a n a n a n a n a n a n ++++?=+Γ-??=-+Γ??
(5.27)
用于更新()j a n 和()R
j a n 。相应地有一对耦合方程更新系统函数()j A z 和()R
j A z , 写成矩阵形式应为:
1
11*1
11
()()1()()j j j R R j j
j A z A z z A z A z z
-++-++??Γ????
=??????Γ???????? (5.28)
这两种表示都描述了一个两端口的网络,如图5.2,它是FIR 格型滤波器的基本实现模块。级联p 个这样的模块,反射系数分别为12,,...,p ΓΓΓ,就形成p 阶的格型滤波器,如图5.3。注意输入()n δ和
输出()p a n 之间的系统函数是()p A z ,而输入和输出()R
p a n 之间的系统函数是()R
p A z 。
图5.2 FIR 格型滤波器模块
图5.3 p 阶FIR 格型滤波器
§5.2.3 Levinson-Durbin 递归算法的参数特性
本节介绍Levinson-Durbin 递归所产生的反射系数序列的一些重要性质,并给出自相关正则方程解的一个重要特性。具体说,我们将证明,当且仅当反射系数的幅度小于1,或等效地,当且仅当p R 是正定矩阵时,()p A z 的根在单位圆内。我们还将证明,由自相关方法所获得的全极点模型保证是稳定的。最后,我们将介绍一个自相关匹配性质,
它阐明了若在自相关法中取(0)b =且若()h n 是(0)()()
p b H z A z =
的逆Z 变换,则()h n 的自相关等于()x n 的自相关。
反射系数的第一个性质是:
该性质的获得是由于j ε是最小的平方误差,即:
{}
min
j j εξ= (5.29)
而
2
()
j n e n ξ∞
==
∑
其中()e n 为j 阶建模的误差。因此0j ε≥。而
2
1[1]j j j εε-=-Γ
由0j ε≥和10j ε-≥得(
)2
10j
-Γ≥,因而1j
Γ
≤。
必须牢记的是,性质1的合法性依赖于j ε的非负性,该非负性所隐含的假设是p R 中的自相关值()x r k 是根据式(4.121)计算的。如果p R 中所用的是任意的数字序列(0),(1),...,(1)x x x r r r p -,则0j ε≥和1j Γ≤并不一定为真。为说明这一点,设正则方程中的(0)1,(1)2x x r r ==,要求解一阶
的模型,这时:
1(1)2(0)
x x r r -Γ=
=-
其幅值大于1。但这个例子与性质1并不冲突,因为对任何合法的自相关序列,都应有(0)(1)x x r r ≥。后面我们还要揭示(性质7),反射系数j Γ的单位幅度约束以及序列j ε的非负性约束,都与矩阵p R 的正定性密切相关。 下一个性质说明了多项式()p A z 的根的位置与反射系数j Γ的幅度间的关系。
理(encirclement Principle)[3]。
闭路原理:给定z 的有理函数:
()
()()
B z P z A z =
设C 是z 平面上的一个简单闭合曲线,如图5.4,当路径C 沿逆时针方向环行一周时,将在()P z 平面产生一个闭合曲线,它沿逆时针方向环绕原点z p N N -次,其中z N 是C 中的零点数,p N 是C 中的极点数。(注:若z p N N -为负,则曲线()P z 按顺时针方向围绕原点z p N N -次。)
图5.4
Z 平面上含一个零点的简单闭合曲线和相应的P(Z )平面路径
闭路原理的几何解释如下。由于C 是一个起点、终点都在0z 的简单闭合曲线,因此()P z 也是一个闭合曲线,起点和终点都为0()P z 。若C 内含一个零点,则曲线绕环一周,零点到曲线上点的角度θ将增加2π。因此在()P z 平面所产生的曲线将沿逆时针方向环绕原点。类似地,若C 内有一个极点,则C 环绕一周将使θ减少2π,在()P z 平面将产生一个按顺时针方向环绕原点的曲线。当C 内有z N 个零点和p N 个极点时,角度的变化为2()z p N N π-,因此有闭路原理。虽然在闭路原理中没有说明()P z 的零点在曲线C 上的结论,但很显然,这时()P z 平面的曲线将要通过原点。 有了闭路原理,就可以用于推导性质2。具体说,对一阶模型
1
11()1A z z
-=+Γ
显然是当且仅当11Γ<时1()A z 是最小相位的。现在我们来证明,若()j A z 是最小相位的,则当且仅当11j +Γ<时1()j A z +是最小相位的。为完成证明,我们利用Levinson 阶更新方程的z 域表达式,即式(5.23),重写为:
1
11()()()R
j j j j A z A z z A z -++=+Γ
上式两边除以()
j
A z,则有:
11
1
()()
()1
()()
R
j j
j
j j
A z A z
P z z
A z A z
+-
+
==+Γ(5.30) 现在要证明若C是沿单位圆的一个闭合路径,则()
P z沿顺时针方向围绕原点的次数等于1
()
j
A z
+
在单位圆外的零点数。由于假定()
j
A z是最小相位的,则()
j
A z在单位圆内有j个零点,在
z=处有j个极点。现在假设
1
()
j
A z
+
在单位圆外有l个零点,圆内有1
j l
+-个零点,由于
1
()
j
A z
+
在0
z=处将有1
j+个极点,因此()
P z沿逆时针围绕原点的次数应为:
(1)(1)
z p
N N j l j l
-=+--+=-
即它顺时针围绕原点l次。现在注意C是单位圆,即j
z eω
=,因此,
1
111
()()
()()
R R j
j j
j j j
j
j j
C C
A z A e
z
A z A e
ω
ω
-
+++
Γ=Γ=Γ
再由式(5.30),显然()
P z是一个中心在1
z=,半径为
1
j+
Γ的圆(参看图5.5)。因而,若
1
1
j+
Γ<,则()
P z不包含且不通过原点。而相反若
1
1
j+
Γ>,则()
P z包含原点,
1
()
j
A z
+
的零点就不是都在单
位圆内。因此,若()
j
A z是最小相位的,则当且仅当
1
1
j+
Γ<时
1
()
j
A z
+
是最小相位的。
图5.5 曲线()
P z是中心在1
z=,半径为
1
j+
Γ的圆
在Toeplitz矩阵
p
R的正定性与正则方程解的最小相位特性之间也有一个本质的联系。尤其是下
面的性质说明了当且仅当
p
R是正定阵时,()
p
A z是最小相位的。(参看§5.2.7节)
性质3:若
p
a是Toeplitz正则方程
1
p p p
ε
=
R a u的解,则当且仅当
p
R是正定阵,即0
p
>
R时,
()
p
A z是最小相位的。
为证明该性质,设α是多项式()
p
A z的根。一般α是复数,我们可将()
p
A z因式分解为:
112(1)
121
1
()1()(1)(1)
p
k p
p p p
k
A z a k z z b z b z b z
α
------
-
=
=+=-++++
∑
我们将证明若0
p
>
R,则1
α<。将
p
a写成因式分解的矩阵形式,有:
1
21
1
110
(1)1
11
(2)
()
p
p
p
p
p
a b
a b b B
b
a p
αα
-
????
????
????
????
===
--
????
????
????
????
??
??
a(5.31) 因此有:
1
1
p p p p
Rε
α
??
==
-??
??
R a B u(5.32)
将上式左乘H B ,有:
110H
p p εα????=-????????
B R B
由于B 是满秩的,则若p R 是正定的,H
p B R B 也是正定的,因此:
1*100H
p s s s s ??
=
>????
B R B 这意味着
2
2
1s s >
(5.33)
又由于
01*
10110H
p p s s s s αεαα-??????
==??????--??????
B R B 因此*01s s α=,即
*
1
s s α=
(5.34)
由式(5.33),应有1α<,即()p A z 是最小相位的。相反,若()p A z 不是最小相位的,则()p A z 将有一个根α满足1α>,因此由(5.34)有10S S >,p R 就不是正定的。
性质1 表明求解自相关正则方程所产生的反射系数,其幅度小于1,性质2给出了若反射系数
幅度小于1,则多项式()p A z 的根在单位圆内。将两者相结合就得如下的性质。
()p A z 的根将会是什么情况。
该性质可以简单说明如下。设j Γ是反射系数序列,且对j p <有1j Γ<,而1p Γ=。由于对
所有j 有1j Γ≤,则由性质2,()p A z 的所有零点都在单位圆内或圆上。若k z 表示()p A z 的零点,则:
()1
1()1p
p k
k A z z
z
-==
-∏
上式两边p z -的系数应相等,因此有:
1
()p
p k k a p z ==
∏
(5.35)
但由于()p p a p =Γ,若1p Γ=,则有:
1
||1p
k
k z
==∏ (5.36)
显然,若有一个()p A z 的零点的幅值小于1,则至少必有一个零点的幅值大于1,否则所有根的乘积
就不会是1。而这与所有根都必须位于单位圆内或圆上的要求是矛盾的,因此性质5成立。
该性质可用于设计检测噪声中的正弦波的约束格型滤波器(第六章习题)。它也是谱估计中平稳随机过程Pisarenko 谐波分解的理论基础(参见第8章)。
下一个性质将给出()x n 的自相关序列与()x n 的全极点模型的单位采样响应()h n 的自相关序列间的关系。具体说,若选择()H z 的分子系数(0)b ,使得()x n 满足(0)(0)x h r r =的能量匹配约束,
响应,即:
1
(0)
()1()p
k
p
k b H z a
k z
-==
+
∑
其中系数()p a k 是如下正则方程的解:
1x p p ε=R a u
(5.37)
因此()h n 满足如下差分方程:
1
()()()(0)()p
p
k h n a
k h n k b n δ=+
-=∑
(5.38)
上式两端乘以*()h n l -并按n 求和,可得:
*
*
*
*
1
()()()()()(0)()()(0)()p
p
n k n n h n h
n l a
k h n k h n l b n h n l b h l δ∞
∞
∞
====-+
--=-=-∑∑∑∑ (5.39)
又由
*
*
()()()()h h n r l h n h
n l r l ∞
==
-=-∑
可用()h r l 重写式(5.39)为:
*
1
()()()(0)()p
h p
h k r l a
k r l k b h l =+
-=-∑ 由于()h n 是因果的,即0n <时()0h n =,而(0)(0)h b =,因此,对0l ≥,有
2
1
()()()(0)()p
h p h k r l a k r l k b l δ=+
-=∑
; 0l ≥ (5.40)
再利用()h r l 的共轭对称性,可将上式写成矩阵形式:
*
*
*
*
*
2**
1(0)
(1)
(2)()1(1)(1)(0)(1)(1)0(2)(0)0(2)(1)(0)(2)......0()()(1)
(2)
(0)h h h h p h h h h p
h
h h h p h h h h r r r r p a r r r r p a b r r r r p a p r p r p r p r ??????
????-??
?
???=??-???????
???????
--????
?
?
5.41) 或者写为:
2
1(0)h p b =R a u
(5.42)
因此,如下的两个Toeplitz 方程的解都是p a :
2
11
(0)h p x p p b ε==R a u R a u (5.43)
现在假设选择(0)b 满足能量匹配约束,即:
2
2
(0)()
()
(0)x h n n r x n h n r ∞
∞
===
=
=∑
∑
(5.44)
下面将证明对0k >有()()x h r k r k =。采用Levinson-Durbin 递归求解方程(5.43),应有:
1(1)(1)(1)(0)
(0)
x h x h r r a r r =-
=-
而(0)(0)x h r r =,因此(1)(1)x x r r =。现在假设对1,2,...,k j =,有()()x h r k r k =,则由(5.10)和(5.13),可得:
1111(1)()(1)(1)()(1)j
x
j j j x i j
h j j j h i r j a i r j i r j a i r j i εε+=+=?
+=-Γ-+-????+=-Γ-+-??
∑∑ (5.45)
由于对1,2,...,i j =,有()()x h r i r i =,因此(1)(1)x h r j r j +=+。最后,由(5.41)可知
2
1(0)
(0)()()p
h p
h k b r a
k r k ==+
∑
而()()h x r k r k =,因此有:
2
1
(0)
(0)()()p
x p
x p k b r a
k r k ε==+
=∑
(5.46)
从而(0)b =
§5.2.4 上行和下行递归
用Levinson-Durbin 递归求解自相关正则方程,就由信号的自相关序列()x r k 获得信号的全极点模型。除了模型参数之外,递归还产生一组反射系数12,,...,p ΓΓΓ以及最后的建模误差p ε。因此,递归可看作是由自相关序列到一组模型参数和一组反射系数的映射:
{}12(1),(2),...,(),(0)(0),(1),...,(),,...,,p p p x x x
p p
a a a p
b r r r p ε?????→?ΓΓΓ??L E V
在很多情况下,若能由滤波器系数()p a k 导出反射系数j Γ,或反过来推导,将会提供很多便利。本节我们介绍如何由一组参数递归地导出另一组。完成这种变换的递归就称为上行递归和下行递归。
式(5.16)给定的Levinson 阶递归方程就是由反射系数j Γ推导滤波器系数()p a k 的递归公式。具体说,由于
*
11()()(1)j j j j a i a i a j i ++=+Γ-+
(5.47)
因此滤波器系数1()j a i +很容易由()j a i 和1j +Γ而获得。递归开始时,取0(0)1a =;当系数()p a k 已被确定,递归结束时,取(0)p b ε=。由于(5.47)所定义的是给定1j +Γ时,j 阶滤波器模型参数如何被(向上地)更新为1j +阶滤波器,因此式(5.47)被称为上行递归。该递归可以图示为:
{}{}1
2,,,...,(1),(2),...,(),(0)p p p p p a a a p b εΓ
ΓΓ???→上行
表5.2总结了该递归。图5.6给出了其Matlab 程序。
图5.6 上行递归的Matlab 程序
______________________________________________________________________________________ 例5.2.3 上行递归
给定反射系数序列为:
[]1
2T
=ΓΓ2Γ
则由上行递归可得一阶和二阶全极点模型如下。由(5.14),一阶模型为:
1111110
(1)01a ????????==+Γ=Γ????????????????
a
类似地得二阶模型为:
***
2212112111222
001111(1)(1)(1)(2)0101a a a a ??????????
????????==+Γ=Γ+ΓΓ=Γ+ΓΓ??????????????????Γ????????????
a
因此,用1Γ和2Γ反射系数所表示的二阶全极点模型的一般形式为:
()*
1
2
21122()1A z z
z
--=+Γ+ΓΓ+Γ
给定了二阶模型2a ,注意若第三个反射系数3Γ再加到其序列中,则3a 很容易由2a 获得: ()()****
11232323222**33**3223222311233311011(1)(1)(2)(1)(2)(2)(2)(2)(1)(1)(3)01a a a a a a a a a a a ????????????Γ+ΓΓ+ΓΓ????+Γ????
==+Γ==??????????+ΓΓ+ΓΓ+ΓΓ??
????????Γ????????Γ????
a 因此,总可以不必再求解较低阶的滤波器而增加滤波器的阶数。
??????????????????????????????????????????????????????????????????????????????????????????????????????? 上行递归是反射系数序列p Γ到滤波器系数序列p a 的映射,即p p →a Γ。在很多情况下,该映
射可能要反过来,即1
:p p L
-→a Γ,是由模型参数p a 确定反射系数。完成这一变换的映射称为下
行或反向Levinson 递归。求反射系数是基于如下事实,即由于
()j j a j Γ=
(5.48)
因此反射系数可以通过反向地运行Levinson-Durbin 递归而计算。具体说,由给定的p a ,取
()p p a p Γ=,然后递归地求出每个较低阶的模型j a (1,2,...,1)j p p =--,并取()j
j a j Γ
=。具
体说明如下:
11111
2222
222
1
1
(1)(2)(2)(1)(1)(2)(2)
(1)(2)(1)p p p p p p
p p p p p p p p p a a a p a p a a a p a a a ---------=--Γ=-Γ=Γ=Γ=
Γ
a a a a a
为看出第j 阶模型是如何由第(1)j +阶模型导出的,我们假设系数1()p a i +已获得,式(5.16)的Levinson 阶递归方程是用()j a i 和(1)j a j i -+表示1()j a i +为:
*
11()()(1)j j j j a i a i a j i ++=+Γ-+
(5.49)
上式只给出了有关两个未知量()j a i 和*
(1)j a j i -+的一个方程。但若我们再利用Levinson 阶递归方程,用(1)j a j i -+和*
()j a i 表示1(1)j a j i +-+,则有:
*
11(1)(1)()j j j j a j i a j i a i ++-+=-++Γ
两边取复数共轭,有:
*
*
*
11(1)(1)()j j j j a j i a j i a i ++-+=-++Γ
(5.50)
它给出了有关两个未知量()j a i 和*
(1)j a j i -+的第二个方程。将式(5.49)和(5.50)联立成矩阵形式,有:
11
***11
()1
()(1)1(1)j j j j j j a i a i a j i a j i ++++Γ
??????
=??????-+Γ-+??????
(5.51)
若11j +Γ≠,则(5.51)中的矩阵可逆,()j a i 可唯一地求出。(注:若11j +Γ=,则(5.51)是奇异
(5.52)
这就是下行递归*11
*1112
1*
11()(1)(1)1
(2)(2)(1).........1()()(1)j j j j j j j j j j j a j a a a a a j a j a j a ++++++++??
????
??
????????
??-=-Γ??????
??-Γ???????
???
?????????
?
(5.53)
一旦获得了序列()j a i ,就可置()j j a j Γ=,再继续递归求出后面的较低阶多项式。下行递归可图示为:
{}{}12(1),(2),...,(),(0),,...,,p
p p p p a
a a p
b ε???→ΓΓΓ下行
表5.3总结了下行递归。
另一种也许是更巧妙的导出下行递归的方式是利用式(5.28),由1()j A z +求j 阶多项式()j A z 。
很容易看出,其解为:
(
)
1
1
11
*21
11
()()1
()()1
1j j j R R j j j j A z A z z z
A z A z z --+++++??-Γ
????
=????
??-Γ??????-Γ
?
? (5.54)
因此有:
(5.55)
它是式(5.52)的Z 域表示。完成下行递归的Matlab 程序如图5.7所示。
图5.7 下行递归的Matlab 程序
_____________________________________________________________________________________ 例5.2.4 下行递归。
假定我们要实现三阶FIR 滤波器:
1
2
3
()10.50.10.5H z z
z
z
---=+--
采用格型滤波器结构。利用下行递归,反射系数可由矢量[]31
0.5
0.1
0.5T
=--a 而获得,
然后用图 5.3所示的结构实现。下行递归首先取33(3)0.5a Γ==-,然后求二阶多项式的系数
[]2221,
(1),
(2)T
a a =a ,是在(5.53)中取2j =,并利用已知的3()a i ,即:
33232
3323
(1)(2)(1)1
(2)(1)(2)1a a a a a a ??????????
=-Γ????????-Γ??????????
由给定的3a 值和1j +Γ,可得:
22(1)0.50,10.61
0.5(2)0,10.50.210.25
a a ?-?????????=+=????
??????--?
????????? 从而22(2)0.2a Γ==。求出2a 后,再利用式(5.52)可求出1a 为:
[]122222
1(1)(1)(1)0.51a a a =
-Γ=-Γ
因此11(1)0.5a Γ==。最后得反射系数序列为:
[]0.5
0.2
0.5T
=-Γ
这样,一个具有给定系统函数的格型滤波器如图5.8所示。
图5.8
格型滤波器实现,123()10.50.10.5H z z z z ---=+--
??????????????????????????????????????????????????????????????????????????????????????????????????????? 下行递归的一个有趣的应用是数字滤波器的Schur-Cohn 稳定性测试[27]。该测试是基于性质2,即当且仅当反射系数的幅值小于1时多项式的根位于单位圆内。因此,若给定一个因果线性移不变滤波器的有理系统函数:
)()
()(z A z B z H =
可测试该滤波器的稳定性如下。首先对分母多项式()A z 的系数应用下行递归,得反射系数序列j Γ,然后,当且仅当所有反射系数的幅值小于1时,滤波器是稳定的。下面举一个例子。
_____________________________________________________________________________________ 例5.2.2 Schur-Cohn 稳定性测试。
试用Schur-Cohn 稳定性测试确定如下滤波器的稳定性:
1
23
1
()243H z z
z
z
---=
+-+
首先注意分母多项式的第一个系数不等于1,因此重写()H z 为:
1
2
3
0.5
()12 1.50.5H z z
z
z
---=
+-+
由33(3)0.5a Γ==,利用式(5.53)求解二阶多项式12222()1(1)(2)A z a z a z --=++为:
33232
3323
(1)(2)(1)2 1.511/31
40.5(2)(1)(2) 1.5210/331a a a a a a ???-??????????????
?=-Γ=-=????????????????---Γ??????????????????
由于22(2)10/3a Γ==-,即21Γ>,因此滤波器不稳定。
???????????????????????????????????????????????????????????????????????????????????????????????????????
§5.2.5 逆Levinson-Durbin 递归
前面已经介绍了如何由自相关序列导出反射系数和模型参数,还介绍了将反射系数序列映射到模型参数序列的递归,以及反向的递归。实际上,也可以由反射系数{}12,,...,p ΓΓΓ加上p ε,或由p 阶模型系数()p a k 加上(0)b ,递归地计算出自相关序列。完成这一变换的递归就称为逆Levinson 递归。为说明该递归如何进行,假设给定了反射系数序列(1,2,...,)j j p Γ=以及p 阶误差p ε。由于
(
)2
1
(0)1p
p x j
j r ε==-Γ∏
(5.56)
因此自相关序列的第一项为:
(5.57)
再加上零阶模型
(5.58)
就构成了该递归的初始化。
现在再假设自相关序列的前j 项已知,并已知j 阶滤波器系数()j a i ,我们将说明如何由此求序
列的下一项(1)x r j +。首先我们利用上行递归由1j +Γ和()j a i 求系数1()j a i +,然后置式(5.2)中的
1p j =+,1k j =+,就得如下的(1)
r j +表达式:
(5.59)
由于上式右边的求和只需要已知的自相关值,该表达式可用于确定(1)x r j +,并完成递归。逆Levinson-Durbin 递归可以图示为:
{}{}1
()
1
2,,...,,(0),(1),...,()p p x x x r r r p ε-Γ
ΓΓ????→LEV
算法总结于表5.4。
例5.2.6 逆Levinson-Durbin 递归
给定反射系数序列123
1/2Γ=Γ=Γ=,模型误差()3
323/4ε=,要求自相关序列[](0)(1)(2)(3)T
x x x x x r r r r =r 。首先是递归的初始化,即:
3
3
2
1
(0)2(1)
x i
i r ε==
=-Γ
∏
再求第一阶模型为:
11110.5????
==Γ???
?????a
从而,
1(1)(0) 1.0x x r r =-Γ=-
又利用上行递归更新模型参数得2a :
2121101011(1)(1)1/21/23/42011/201a a ??????????
=+Γ=+=??????????????????????????????
a
从而
2231(2)(1)(1)(2)(0)144
x x x r a r a r =--=
-=-
再利用上行递归由2a 得3a :
223322101011(1)(2)3/41/21(2)(1)1/23/47/82011/201a a a a ????????????????????=+Γ=+=???????????????
???????????????a
最后得(3)x r :
333171(3)(1)(2)(2)(1)(3)(0)1488
x x x x r a r a r a r =---=
+-=
要求的自相关序列为:
112,1,,48T
x ?
?=--???
?r
?????????????????????????????????????????????????????????????????????????????????????????????????????? 上面介绍的逆Levinson 递归给出了由反射系数和p ε求出自相关序列的步骤,若给定的不是反射
系数,而是滤波器的系数()p a k ,仍可以确定自相关序列,即先用下行递归求出反射系数,然后再利用上述的逆Levinson 递归。由反射系数j Γ或滤波器系数()p a k 求解()x r k 的Matlab 程序如图5.9所示。 现在来总结一下。将§5.2.4和§5.2.5结果结合到一起可看出,如下的三组参数可以互相等效: 1.自相关序列()x r k
2.全极点模型的()p a k 和(0)b
3.反射系数j Γ以及误差p ε
这种等效性如图5.10所示,它说明了如何由一组参数导出另一组。
图5.9
逆Levinson-Durbin 递归的Matlab 程序
图5.10 Levinson-Durbin 递归
§5.2.6 Schur 递归(暂省略)
§5.2.7 Cholesky 分解
我们已看到Levinson-Durbin 递归是求解自相关正则方程的有效算法。它也可用于完成哈密顿Toeplitz 自相关阵p R 的Cholesky 分解。哈密顿矩阵C 的Cholesky 分解(LDU)形式为:
H
=C LDL
(5.81) 其中L 是下三角阵,且对角元素为1,D 是对角阵。若D 的对角元素为非负(C 是半正定的),则D
可以分解为两个平方根矩阵的乘积形式:
1
1
22
=?D D
D
(5.82) 这样,C 的Cholesky 分解可以表达为一个上三角阵和一个下三角阵的乘积:
(
)(
)1
1
2
2
H
=C L D
D
L
(5.83)
利用自相关矩阵的Cholesky 分解,可以很容易地导出p R 的正定性、误差序列j ε的正实性以及反射系数幅值的单位约束这三者间的等效性。另外,我们还可以导出自相关阵的逆阵的闭合表达式,以及Toeplitz 矩阵求逆的递归算法。 为导出p R 的Cholesky 分解,考察一个(1)(1)p p +?+的上三角阵:
*
*
*
12*
*
2*
1
(1)(2)...()
01(1)...(1)0
1...(2) 0
1p p p p a a a p a a p a p ??
??
-??=-????????A
(5.84)
该矩阵是Levinson-Durbin 递归作用于自相关序列(0),...,()x x r r p 时所产生的矢量
01,,...,p a a a 所形成。注意p A 的第j 列是滤波器系数1R
j -a 再补足零。由于
R
j j j j ε=R a u
(5.85)
其中[]0,
0,,
1T
j = u 是长度为1j +的单位矢量,最后一个元素为1。这样就有
1
2
000
00p p p εεεε??*
????*
*=????**
*
?
?
R A (5.86)
它是一个下三角阵,对角线上为预测误差(*号表示那些非零的元素)。虽然式(5.85)的结果没有给出推导,但它可由p R 的哈密顿Toeplitz 特性而推知。例如,由于*
p p =JR J R (J 是幂等阵),则有
*
1()j j j j j ε==R a JR J a u
(5.87)
两边都左乘J ,并利用2
=J
I 和*
()R j j =Ja a ,再取复数共轭,(5.87)就变成了(5.85)。由于
两个下三角阵的乘积仍是下三角阵,因此若用H
p A 左乘p p R A ,则得到另一个下三角阵H
p p p A R A 。注意由于p A 的对角项为1,因此H
p p p A R A 的对角项与p p R A 相同,从而H
p p p A R A 的形式也是:
01
2
000
00H
p p p p εεεε??*
????*
*=????**
*
?
?
A R A (5.88)
但更重要的是注意H
p p p A R A 是哈密顿的,因此式(5.88)右边的矩阵也必是哈密顿的,显然其对
第五章离散时间信号的数字处理 Q5.1运行程序P5.1,产生连续时间序号及其抽样形式,并显示它们。 clf; t = 0:0.0005:1; f = 13; xa = cos(2*pi*f*t); subplot(2,1,1) plot(t,xa);grid xlabel('时间, msec');ylabel('振幅'); title('连续时间序号 x_{a}(t)'); axis([0 1 -1.2 1.2]) subplot(2,1,2); T = 0.1;n = 0:T:1; xs = cos(2*pi*f*n); k = 0:length(n)-1; stem(k,xs);grid; xlabel('时间 n');ylabel('振幅'); title('离散事件序号 x[n]'); axis([0 (length(n)-1) -1.2 1.2]) Q5.2 正弦信号的频率是多少赫兹?抽样周期是多少秒? 正弦信号的频率f=13Hz,抽样周期T=0.1s。 Q5.3 解释两个axis命令的效果。 给x,y轴标刻度。 Q5.4 以比在程序P5.1中列出的抽样周期低的两个抽样周期和高的两个抽样周期的四个其他值,运行程序P5.1.评论你的结果。 T=0.04s T=0.08s
T=0.15s T=0.3s 由上图可以发现:当取的T越小时,得到的图形越接近原图形。 Q5.5 通过将正弦信号的频率分别变为3HZ和7HZ,重做习题Q5.1。相应的等效离散时间信号与习题Q5.1中产生的离散时间信号之间有差别么?若没有,为什么没有? f=3Hz f=7Hz 由图可以看出,变换频率得到的两个图没有区别,因为他们的抽样周期一样。 Q5.6 运行程序P5.2,产生离散时间信号x[n]及其连续时间等效ya[t],并显示它们。 clf;T = 0.1;f = 13;n = (0:T:1)';
a=input('type in the first sequence ='); b=input('type in the second sequence ='); c=conv(a,b); M=length(c)-1; n=0:1:M; disp('output sequence =');disp(c) stem(n,c) xlabel('Time index n');ylabel('Amplitude'); type in the first sequence =[2 4 6 4 2 0 0 0] type in the second sequence =[3 -1 2 1] output sequence = Columns 1 through 9 6 10 18 16 18 12 8 2 0 Columns 10 through 11 0 0 ??? Undefined function or variable 'ylable'. Error in ==> E:\Matlab6p5FULL\bin\win32\Untitled.m On line 8 ==> xlabel('Time index n');ylabel('Amplitude'); type in the first sequence =[2 4 6 4 2 0 0 0] type in the second sequence =[3 -1 2 1] output sequence = Columns 1 through 9 6 10 18 16 18 12 8 2 0 Columns 10 through 11
DSP 原理与应用The Technology & Applications of DSPs 第五章: TMS320F28335片内外设 北京交通大学电气工程学院 夏明超郝瑞祥万庆祝 mchxia@https://www.doczj.com/doc/445223453.html, haorx@https://www.doczj.com/doc/445223453.html, qzhwan@https://www.doczj.com/doc/445223453.html, :TMS320F28335第五讲: TMS320F28335片内外设教学目标:
掌握TMS320F28335内核结构,例如A/D转换、串行通信接口、串行外设接口。 外设接重点: TMS320F28335A/DCS308335内部/C 的正确使用,串行通信接口应用。难点: TMS320F28335的ADC 寄存器操作和串行通信寄存器操作。教学内容分两部分 51§5.1:TMS320F28335内模拟/数字转换 §5.2 :TMS320F28335系列串行通信接口SCI 和Modbus 协议介绍DSP 原理与应用2
DSP 原理与应用3 ADC 有关引脚
§5.1 TMS320F28335 内模拟/数字转换§5.1 .1Features and functions of ADC module:◆core with built-in dual sample-and-hold◆Simultaneous sampling or sequential sampling modesp g q p g ◆Analog input: ◆Fast conversion time runs at ADC clock or Fast conversion time runs at , ADC clock, or 6.25 MSPS multiplexed inputs ◆, multiplexed inputs◆capability provides up to 16 " t i " i i l i E h i "autoconversions" in a single session. Each conversioncan be to select any 1 of 16 input channels.DSP 原理与应用4 Sequencer can be operated as two independent 8-state ◆Sequencer can be operated as two independent 8-state sequencers or as one large 16-state sequencer (i.e., two cascaded 8-state sequencers two cascaded 8state sequencers.◆(individually addressable to store conversion values store conversion values A/DC digital value:
第五章FIR 滤波器的设计附加题 1. 一FIR 数字滤波器的传输函数为 12341()[1242]30 H z z z z z ----=++++ 求()h n 、幅度()H ω、和相位()?ω。 2. 一个FIR 线性相位滤波器的h(n)是实数,且n<0和n>6时,h(n)=0。如果h(0)=1且系统函数在30.5j z e π=和z = 3处各有一个零点,求()H z 。 3. 一个FIR 线性相位滤波器的h(n)是实数,且n<0和n>6时,h(n)=0。如果h(0)=1且系统函数在30.5j z e π=和z = 3处各有一个零点,求()H z 。 4. 如图所示h 1(n)为N=8的偶对称序列,h 2(n)为其循环右移4位后的序列。 设H 1(k)=DFT [ h 1(n) ],H 2(k)=DFT [ h 2(n) ] (1)问 | H 1(k)| = | H 2(k)| 吗?1()k θ与2()k θ的关系是什么? (2)若h 1(n)、h 2(n)各构成一个低通滤波器,问它们是否是线性相位的?延时分别是多少? (3)两个滤波器的性能是否相同? 01234567 01234567 h 1(n) h 2(n) 5. 用矩形窗设计一线性相位FIR 低通数字滤波器 0()0j a c j d c e H e ωω ωωωωπ-?≤≤?=?<≤?? (1)求h d (n); (2)求h(n),并确定a 与N 的关系; (3)讨论N 取奇数和偶数对滤波器性能有什么影响。 6. 设计一个线性相位FIR 低通滤波器,给定抽样频率为42**1.5*10(/sec)s pi rad Ω=,通带截
第五章信号处理初步 测试工作的目的是获取反映被测对象的状态和特征的信息。但是有用的信号总是和各种噪声混杂在一起的,有时本身也不明显,难以直接识别和利用。只有分离信号与噪声,并经过必要的处理和分析、清除和修正系统误差之后,才能比较准确地提取测得信号中所含的有用信息。因此,信号处理的目的是: 1)分离信、噪,提高信噪比; 2)从信号中提取有用的特征信号; 3)修正测试系统的某些误差,如传感器的线性误差、温度影响等。 信号处理可用模拟信号处理系统和数字信号处理系统来实现。 模拟信号处理系统由一系列能实现模拟运算的电路,诸如模拟滤波器、乘法器、微分放大器等环节组成。其中大部分环节在前行课程和前面几章中已有讨论。模拟信号处理也作为数字信号处理的前奏,例如滤波、限幅、隔直、解调等预处理。数字处理之后也常需作模拟显示、记录等。 数字信号处理是用数字方法处理信号,它即可在通用计算机上借助程序来实现,也可以用专用信号处理机来完成。数字信号处理机具有稳定、灵活、快速、高效、应用范围广、设备体积小、重量轻等优点,在各行业中得到广泛的应用。 第一节数字信号处理的基本步骤 1.数字信号处理的基本步骤如图5-I所示。 信号的预处理是把信号变成适于数字处理的形式,以减轻数字处理的困难。 预处理包括: 1)电压幅值调理为便于采样,总是希望电压峰-峰值足够大,以便充分利用A/D换器的精确度。如12位的A/D转换器,其参考电压为 5V。由于2l2=4096,故其末位数的当量电压为2.5mV。若信号电平较低,转换后二进制数的高位都为0,仅在低位有值,转换后的信噪比将很差。若信号电平绝对值超过5V,则转换中又将发生溢出,这是不允的。所以进入A/D 转换的信号的电平应适当调整。 2)必要的滤波,以提高信噪比,并滤去信号中的高频噪声。
《机械工程测试技术》 第五章 数字信号处理初步 主讲:王建军 山东理工大学?机械工
第五章信号处理初步 ●测试的目的:获取被测对象的状态和特征的信息。但信号 总是与噪声混杂在一起。所以,有必要进行信号处理。●信号处理的目的: ?1)分离信、噪,提高信噪比。 ?2)从信号中提取有用的特征信息。 ?3)修正测试系统的某些误差,如:传感器的线性误差、温度影响。 ●信号分析:研究信号的构成和特征值。 ●信号处理:信号经过必要的变换以获取所需信息的过程。 ●信号处理分为两类:模拟信号处理和数字信号处理
模拟信号处理: ●实现模拟运算的电路,如模拟滤波器、乘法器、 微分放大器等。 ●模拟信号处理也可用于数字信号处理的前奏 (如滤波、限幅、隔直、解调)及后续处理 (如模拟显示、记录)。
数字信号处理: ●用数字方法处理信号,可采用通用计算机, 或专用的信号处理机实现。 ●数字信号处理技术目前正处于迅速的发展 阶段,如DSP芯片的开发与使用,势头很 好。
第一节数字信号处理的基本步骤 预处理A/D 转换数字信号处理器或计算机 A/D 转换 结果显示 预处理 x(t)y(t) 物理信号 x(t) 传感器 电信号信 号调理 电信号 A/D 转换数字信号数字信号 分析仪或计 算机 显示 物理信号 y(t) 传感器 电信号 信号调理 电信号 A/D 转换数字信号
?1、信号的预处理:把信号变成适于数字处理的形式,减轻数字处理的困难。 ●1)电压幅值调理,便于采样。例如:12位A/D 转换器,参考电压为±5V ,其末位数字的当量电压为2.5mV 。●2)必要的滤波,提高信噪比,虑去信号中的高频噪声。●3)隔离信号中的直流分量(如果所测信号不允许有直流分量)。 ●4)对调制信号进行预先解调。 预处理A/D 转换数字信号处理器或计算机 A/D 转换 结果显示 预处理 x(t) y(t)
第五章 数字滤波器 一、数字滤波器结构 填空题: 1.FIR 滤波器是否一定为线性相位系统?( )。 解:不一定 计算题: 2.设某FIR 数字滤波器的冲激响应,,3)6()1(,1)7()0(====h h h h 6)4()3(,5)5()2(====h h h h ,其他n 值时0)(=n h 。试求)(ωj e H 的幅频响应和 相频响应的表示式,并画出该滤波器流图的线性相位结构形式。 解: {}70,1,3,5,6,6,5,3,1)(≤≤=n n h ∑-=-=1 )()(N n n j j e n h e H ωω ??? ? ??++???? ??++???? ??++???? ??+=+++++++=---------------ωωωωωωωωωωωωωωωωωωω21 21272323272525272727 2 7 7654326533566531j j j j j j j j j j j j j j j j j j j e e e e e e e e e e e e e e e e e e e )(27 )(27c o s 225c o s 623c o s 102cos 12ωφω ωωωωωj j e H e =????? ???? ??+??? ??+??? ??+??? ??=- 所以)(ω j e H 的幅频响应为 ω ωωωωω27 27cos 225cos 623cos 102cos 12)(j e H -????????? ??+??? ??+??? ??+??? ??= )(ωj e H 的相频响应为 ωωφ2 7 )(-= 作图题: 3.有人设计了一只数字滤波器,得到其系统函数为: 2 11 2113699.00691.111455 .11428.26949.02971.114466.02871.0)(------+-+-++--=z z z z z z z H 2112570.09972.016303.08557.1---+--+z z z
第五章信号处理初步 一、知识要点及要求 (1)了解信号处理的目的和分类,及数字信号处理的基本步骤; (2)掌握模拟信号数字化出现的问题、原因和措施; (3)掌握信号的相关分析及其应用; (4)掌握信号的功率谱分析及其应用。 二、重点内容及难点 (一)信号处理 1、信号处理的目的 (1)分离信号和噪声,提高信噪比; (2)从信号中提取有用的特征信号; (3)修正测试系统的某些误差,如传感器的线性误差、温度影响等。 2、信号处理的分类 模拟信号处理:对模拟信号进行处理,由一系列能实现模拟运算的电路来实现。 数字信号处理:对数字信号进行处理,可以在通用计算机上借助程序来实现,或由专用数字信号处理机(DSP芯片)来实现。 (二)数字信号处理的基本步骤 1、 (1)电压幅值调整;(2)必要的滤波;(3)隔直;(4)解调。 2、A/D转换的作用:把模拟信号转换为数字信号,以便能用数字方法进行处理。 (1)采样:时间离散;(2)量化:幅值离散;(3)截断。 3、计算机或数字信号处理器的作用对数字化之后的信号进行处理。 (三)模拟信号的数字化 1、时域采样和混叠 时域采样,就是等时间间隔地取点。从数学处理上看,就是乘以采样函数,时域相乘相当于频域作卷积,就相当于频谱的周期延拓,即频谱的搬移。 在频域中,如果频谱的搬移距离过小,搬移后的频谱就会有一部分相互交叠,从而使新合成的频谱与原频谱不一致,无法准确地恢复原时域信号,这种现象称为混叠。 2、时域截断和泄漏 时域截断,就是取有限长的信号。从数学处理上看,就是乘以有限宽矩形窗函数。时域相乘相当于频域作卷积,就相当于频谱的周期延拓,即频谱的搬移。 在频域中,由于矩形窗函数的频谱是一个无限带宽的sinc函数,即使原模拟信号是有限带宽的,截断后也必然成为无限带宽的,这种信号的能量在频率轴分布扩展的现象称为泄漏。 3、频域采样和栅栏效应 频域采样,就是在频率轴上等间隔地取点,使频率离散化。从数学处理上看,就是乘以频率采样函数。频域相乘相当于时域作卷积,就相当于时域波形的周期延拓,即频域波形的