北京邮电大学电信工程学院多媒体通信中心门爱东电信工程学院多媒体通信中心
门爱东教授
menad@https://www.wendangku.net/doc/a713295627.html,
数字信号处理
Digital Signal Processing
第5 章FIR 数字滤波器设计和实现
2
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
主题概述
1 -绪论
2 -离散时间信号和离散时间系统
3-离散傅里叶变换及其快速计算方法4-IIR 数字滤波器设计和实现
5–FIR 数字滤波器设计和实现
5.1) 概述
5.2) 线性相位FIR DF 约束条件和频率响应5.3) 窗函数法5.4) 频率取样法
5.5) FIR 数字滤波器的优化设计5.6) FIR 数字滤波器的实现结构5.7) 附录
5.8)本章小结
6 –数字信号处理中的有限字长效应
3北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.1 概述:IIR 和FIR 比较
IIR 与FIR 性能特性比较
IIR 数字滤波器:
幅频特性较好;但相频特性较差; 有稳定性问题;
FIR 数字滤波器:
可以严格线性相位,又可任意幅度特性 因果稳定系统 可用FFT 计算
但阶次比IIR 滤波器要高得多
4
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.1 概述:IIR 和FIR 比较
IIR 与FIR 设计方法比较
IIR DF :
无限冲激响应,H(Z) 是z -1的有理分式,借助于模拟滤波器设计方法,阶数低(同样性能要求)。其优异的幅频特性是以非线性相位为代价的。
缺点:只能设计特定类型的滤波器,不能逼近任意的频响。
FIR DF :
有限冲激响应,系统函数H(Z) 是z -1的多项式,采用直接逼近要求的频率响应。设计灵活性强
缺点:①设计方法复杂;②延迟大;③阶数高。
(运算量比较大,因而在实现上需要比较多的运算单元和存储单元)
FIR DF 的技术要求:
通带频率ωp ,阻带频率ωs 及最大衰减αp ,最小衰减αs 很重要的一条是保证H(z) 具有线性相位。
5
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.1 概述:FIR DF 设计方法
FIR 数字滤波器
设计FIR 滤波器的任务:
给定要求的频率特性,按一定的最佳逼近准则,选定h(n) 及阶数N 。
三种设计方法:
n 窗函数加权法o 频率采样法
p FIR DF 的CAD --切比雪夫等波纹逼近法
6
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.1 概述:FIR DF 零极点
FIR 滤波器的I/O 关系:
1
0N r y(n)h(r)x(n r)
?==?∑0121
(), ,,,...,=?h n n N FIR 滤波器的系统传递函数:
121
1
011N N N r
N r h()z h()z .....h(N )H(z)h(r)z
z ?????=++?==
∑?在Z 平面上有N-1 个零点;在原点处有一个(N-1)阶极点,永远稳定。
FIR 系统定义:一个数字滤波器DF 的输出y(n),如果仅取决于
有限个过去的输入和现在的输入x(n), x(n-1),. ......, x(n-N+1),则称之为FIR DF 。
FIR 滤波器的单位冲激响应:
7北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
FIR DF 的频率响应为:
FIR 滤波器的最重要特点是能实现线性相位。具有线性相移特性的FIR 滤波器是FIR 滤波器中应用最广泛的一种。
1
j ()
r 0
()()H ()e N j j n n H e h n e ω
ωθωω??==
=∑H r (ω):振幅响应,它是一个取值可正可负的实函数。θ(ω) = arg [H(e jw )] 为数字滤波器的相位函数。
5.1 概述:FIR DF 频率响应
8
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
信号通过线性滤波器时,其幅度和相位都将发生改变,滤波器增益|H(ω)|和相位θ(ω)都会随频率的变化而改变。
如:输入正弦信号Acos(n ω0)
则:输出为|H(ω0)| Acos(n ω0+θ),其中相移θ=θ(ω0)
输出频率和输入频率相同,但幅度和相位都发生了变化 输出信号比输入信号滞后的样点数n (位移)可由下式求得:
设:n ω0+θ=0
000
()n -
-θωθ
ωω==-滤波器在数字频率ω0 处的相位延迟(位移)
由于相位延迟n 的不同,最终产生了相位失真。
确保不产生相位失真的办法:使不同频率的信号通过滤波器时有
相同的延迟。
5.1 概述:相位失真
9
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
对不同的频率有恒定的相移,会产生相位失真.
如:方波y(t) 可以用无数正弦波的叠加来得到
4
1111
y(t)[sin(t)sin(3t)sin(5t)sin(7t)sin(9t)3579
π
=
Ω+Ω+Ω+Ω+Ω+
若每个正弦波相移π/2 弧度:
确保所有频率具有相同相位延迟的简单方法:
随着频率的变化而改变相位,使滤波器具有线性相位特性,即使所有频率的相位延迟保持恒定,这种方法可通过使系统的相位函数θ(ω)为频率ω的线性函数来实现。
5.1 概述:相位失真
p 4
1111
y (t)=[cos(t)cos(3t)cos(5t)cos(7t)cos(9t)]
3579
πΩ+Ω+Ω+Ω+Ω+ 可见相移之后正弦波之和已不再是方波。10
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
主题概述
1 -绪论
2 -离散时间信号和离散时间系统
3-离散傅里叶变换及其快速计算方法4-IIR 数字滤波器设计和实现
5–FIR 数字滤波器设计和实现
5.1) 概述
5.2) 线性相位FIR DF 约束条件和频率响应5.3) 窗函数法5.4) 频率取样法
5.5) FIR 数字滤波器的优化设计5.6) FIR 数字滤波器的实现结构5.7)本章小结5.8) 附录
6 –数字信号处理中的有限字长效应
11北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2 线性相移FIR DF 约束条件和频率响应
三个内容:
n 约束条件
恒延时滤波
偶对称:恒相延时和恒群延时同时成立 奇对称:仅恒群延时成立
o 频率响应
Type I :h(n) 偶对称、N 为奇数 Type II :h(n) 偶对称、N 为偶数 Type III :h(n) 奇对称、N 为奇数 Type IV :h(n) 奇对称、N 为偶数
p FIR DF 零极点分布
12
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
相延时:()()p θωτωω
=?
群延时:()()g d d θωτωω
=?
5.2.1 线性相移FIR DF 约束条件:恒延时滤波
恒延时滤波
滤波器的延时有相延时和群延时两种
1
()()
()()|()|()N j jn j j j r n H e h n e H e e H e ωωωβωθωω??===±=∑令
j ()=
arg [H(e )]
ωθω恒延时滤波器:τp (ω) 或τg (ω) 是不随ω变化的常量,这时滤波器具有线性相位特性。
13北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
()θωτω
=?(负号是因为系统必有时延)
由于FIR 滤波器的传递函数为:
10
10
()()()[cos sin ]
N j j n
n N n H e h n e
h n n j n ω
ωωω??=?===?∑∑w
θ(w)
10
1
0()sin ()arg[()]arctan ()cos N j n N n h n n H e h n n ωωθωτωω?=?=????==?=?????????
∑∑故:
5.2.1 线性相移FIR DF 约束条件:恒延时
恒相延时和恒群延时同时成立
要使τp 、τg 都不随ω变化,θ(ω) 必须是一条过原点直线
14
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
于是:
1
01
()sin sin tan()cos ()cos N n N n h n n
h n n
ωτω
τωτω
ω?=?===
∑∑1
10
()sin cos ()cos sin N N n n h n n h n n
τωωτωω??===∑∑1
()in()N n h n s n τωω?=?=∑][][]{}[][][][]{}h(0)sin()h(N-1)sin h(n)sin -n h(N-0
{ h(1)sin h(0)sin()h(N-1 h (1 1 -( h N-1)(N-2)s )s )s in in -(-n)si (n N (N-1-n -2) N-1}in 0
))++=+ττωτωωω+
τωω+τωτωωτωω+τ++=ωωτωωω
--+--5.2.1 线性相移FIR DF 约束条件:恒延时
15
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
可以证明,当
1
12
()[()] (0n N-1)
N h n h N n ?τ==??≤≤且 1
2
()()p g N ?τω=τω=τ=
5.2.1 线性相移FIR DF 约束条件:恒延时
上式成立,此时
恒相延时和恒群延时同时成立的线性相位滤波器的必要条件是:
不管N 为偶数,还是N 为奇数,系统冲激响应h(n) 都关于中心点(N-1)/2 偶对称。当N 为奇数时对称中心轴位于整数样点上; 当N 为偶数时对称中心轴位于非整数样点上。
h(n) 为偶对称,N 为偶数
1
2
N ?7n
h(n)h(n) 为偶对称,N 为奇数
12
N ?6
n
h(n)16
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
02
()π
θω=θ?τω=
?τω于是有:
[]1
1
()sin cos()tan cot 2sin()
()cos N n N n h n n
h n n
ωπτωτωτωτωω?=?=???
?==
=
????
∑∑110
j ()r ()H ()e ()()[cos sin ]N N j j n n n H e h n e h n n j n ωθωωωωω???=====?∑∑5.2.1 线性相移FIR DF 约束条件:恒群延时
只要求恒群延时成立
若只要求群延时τg (ω) 为一常数,则相移特性为不过原点的直线。
ω
θ(ω)
2
π10
1
0()sin ()arg[()]arctan 2()cos N j n N
n h n n H e h n n ωωπθωτωω?=?=??????==?=???????
∑∑故
17
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 11
()cos cos ()sin sin N N n n h n n h n n ωτωωτω
??===?∑∑1
()cos()0
N n h n n τωω?=?=∑可以证明,当
1
12
()[()] (0n N-1)N h n h N n ?τ=
=???≤≤且 1
2
()()p g N ?τω=τω=τ=
上式成立,此时
故
5.2.1 线性相移FIR DF 约束条件:恒群延时
18
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
FIR 滤波器单独满足恒定群延时的必要条件为:
冲激响应h(n) 对中心点(N-1)/2 成奇对称。此时,无论N 为奇数或偶数,滤波器的相频特性均为线性,并包含有π/2 的固定相移:
因此,信号通过此类滤波器时不仅产生(N-1)/2 个取样点的延迟,还将
产生90o 的相移,通常这类滤波器又被称为90o 移相器,并具有很好的应用价值。N-1
()2
2
πθωω
?=1111222N N N h h N h ?????????
=???=??????????
???当N 为奇数时
,故102N h ???
=?
???
12
N ?7h(n) 为奇对称,N 为偶数n
h(n)0
12
N ?6h(n) 为奇对称,N 为奇数
n
h(n)
5.2.1 线性相移FIR DF 约束条件:恒群延时
19
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
00,()1 2()(1)
N h n h N n θθωτω
τ==???
??
=??
=????相时延和群时延同时成立 奇对称:θ(ω) 对所有的频率成分都有一个90°相移。
因此,有四种类型的FIR DF :
5.2.1 线性相移FIR DF 约束条件
线性相位约束条件
对于任意给定的值N ,当FIR 滤波器的h(n) 相对其中心点(N-1)/2 是对称时,不管是偶对称还是奇对称,此时滤波器的相移特性是线性的,且群延时都是τ= (N-1)/2 。
偶对称:θ(ω) 为过原点的,斜率为-τ的一条直线
0,()221 2()(1) N h n h N n ππθθωτωτ?
==???
??
=
??
=??????
仅群时延同时成立I N II N III N IV N ????
???类型 : h(n)偶对称,为奇数类型: h(n)偶对称,为偶数类型: h(n)奇对称,为奇数类型: h(n)奇对称,为偶数
20
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type I
h(n) 偶对称,N 为奇数(恒相时延、恒群时延〕
此时,由于h(n) 序列的长度为奇数,因此滤波器的频率响应函数可进行以下拆分(前后对称部分、中心点):
1
11
1
1
21
2
00
12
1()()()()2N N N N j w
jw
jnw
j N n nw jnw n n N H e h n e
h n e h n e h e ?????????=?+==???==
+
+????∑∑
∑
h(n) 为偶对称,N 为奇数
1
2
N ?6
n
h(n)对上式的第二和式作变量替换(n=N-1-m) 后得到:
1
1
111
2
2
(1)20
1()()(1)2N N N j w
jw
jnw
j N w jnw n n N H e h n e
h N n e e h e ?????????==???=
+
??+??
??
∑
∑
由对称条件1()()
h n h N n =??则H(e j ω) 表示为:
21
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type I
1
11
2
2
01
11
22
1
11
22
(1)11
22
1()()21()21()12c s 2o 2
N N j
j n N N j
n N N jn j N jn N N j j j j w jn n n N H e h n h e
N e
h h n N e
h h e
e
e
e e n e e N n ω
ω
ω
ωωωωω
ωω
ω???????????=????=????=?????=
+??????
?????????=+??????????++?????????????????????=+??????∑
∑∑?
??
?????
令1
'2
N n n ?=
?则上式为
1
122112()2
'10
11()2(')cos '22()cos ()
N N j j N N j j r n n N N H e e h h n n e
a n n e H ωω
ωθωωωω=?????=?????????=+?????????
??
==∑∑22
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
由此可以看出其线性相位特性。由于cos(n ω) 对于ω=0、π、2π都是偶对称,所以振度响应H r (ω) 对ω=0、π、2π也是偶对称。
5.2.2 线性相移FIR DF 频率响应:Type I
其中
1()()21()2()2N a n h n 0N a n h n n 0??
==???
??=?≠??
振幅响应:
12
r H ()()cos N n a n n ωω
?==
∑相频响应:
N-1
()--2
θωτωω==0
2
46
8
00.20.4
0.60.8
1h(n)
2
46
8
00.20.40.60.81a(n)
1
2
123
frequency Unit:pi M a g n i t u d e
Magnitude Response
1
2
-30-20-100frequency Unit:pi
P h a s e
Phase Response
N=9
Hr (w)
23
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
h(n) 偶对称,N 为偶数(恒相时延、恒群时延〕
由于h(n) 序列的长度为偶数,因此滤波器的频率响应函数可拆分成如下两部分(前后对称部分,中心点处无值):
5.2.2 线性相移FIR DF 频率响应:Type II
h(n) 为偶对称,N 为偶数
12
N ?7n
h(n)对上式的第二和式作变量替换(n=N-1-m) 后得到:
由对称条件
1()()
h n h N n =??则H(e j ω) 表示为:
11
1
2
2
()()()()N n N N N j jn jn jn n n H e h n e
h n e
h n e
ω
ω
ω
ω
??????=====+
∑∑∑1122
10
1()()()()N N j jn j N jn n n H e h n e h N n e e ωωωω
?????===+??∑∑24
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type II
令'2
N
n n
=?,则上式为:
12
1
122
11
22
()()()()cos[(
)]jn j N w j n N
n n N j N j
H e h n N e
n e e n e h ω
ωωωω?=?????=???=??
?+?=∑∑1
1
2
1
22
1
21
2221
2()
()()cos[()]
()cos[()]()
N
N j
w jw n N
N j
w j n r
N H e e
h n n w e
b n n w e H θωω??=??==??=?=∑∑其中
222
12
()(
),,.....,/N
b n h n n N =?=(注意n 从1 开始,即b(0)=0,或没有定义)
25
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type II
与所设计的b(n) 或h(n) 无关,恒为0。这种类型(即h(n) 偶对称,N 为偶数)不能用于高通或带阻滤波器。
2)由于cos[(n-1/2)ω] 对于ω=π是奇对称,所以,H r (w) 对ω=π也是奇对称;以ω=0、2π为偶对称。
振幅响应:
N 2
r n 1
1H ()b(n) cos (n-)2ωω=?
?=??
??∑相频响应:
N-1
()--2
θωτωω==02
46
0.2
0.40.60.81
h(n)
2
46
8
00.20.40.60.8
1b(n)
1
2
00.511.52
frequency Unit:pi M a g n t u d e
Magnitude Response
1
2
-30
-20-100frequency Unit:pi
P h a s e
Phase Response
N=8n 从1开始
Hr (w)
注意:
1) 在ω= π处,有:
2
110
2()()cos N
r n H b n n ππ=??
??=?=?????
???∑26
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
1
1
1
1112
2
2
(1)(1)00
1
1
111
12
2
2
2
1()2
2()()(1)()[]
1
1
()2sin[()]2()sin[(
)]2
2
12()sin 2N N N j jn j N jn jn j N jn n n n N N N N j
j n n N j H e h n e h N n e e h n e e e N N e
h n j n e h n n n j N e
h n ωωωωωωωωωπωωω????????????===????????==??=
+
??=
???=?=???=?∑
∑
∑
∑
∑
1120N n ω??=??
????
?????????
∑ h(n) 奇对称,N 为奇数(恒群时延〕
h(n) 长度为奇数,拆分成前后两部分:
5.2.2 线性相移FIR DF 频率响应:Type III
对上式的第二和式作变量替换,并利用对称条件h(n)=-h(N-1-n),得:
1
111210
2
()()()()N N N j jn jn jn N n n n H e h n e
h n e
h n e ω
ω
ω
ω
????????===
==
+
∑∑∑
1
2
N ?6h(n) 为奇对称,N 为奇数
n
h(n)27
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
Hr (w)
5.2.2 线性相移FIR DF 频率响应:Type III
1
'2
N n n ?=
?,则上式为:1
1
1
22
2
[
]()
()()sin ()
N N j
j j r
n H e e
c n n e H πωωθωωω??
=?==∑其中1212212
(),,.....,()/N c n h n n N ???
=?????=?令振幅响应:
1
2
1()()sin
N r n H c n n ωω?==
∑()相频响应:
2N-1
()-2
πθωω
=
02
468-0.4
-0.200.20.40.6h(n)
2
46
8
0.2
0.4
0.6
0.81c(n)
1
2
00.511.5
frequency Unit:pi M a g n i t u d e
Magnitude Response
1
2
-30-20-10010
frequency Unit:pi
P h a s e
Phase Response 0.5pi
n 从1开始
28
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
与c(n) 或h(n) 的值无关,因此,这种类型的滤波器不适用于低通、带阻或高通滤波器设计,而且,这说明jH r (w) 是纯虚数,对于逼近理想数字希尔伯特变换和微分器,它是很有用的。理想的希尔伯特变换是一个全通滤波器,它对输入信号产生90 度的相移,它频繁用于通信系统中的调制。微分器广泛用于模拟和数字系统中对信号求导。
2)由于sin(n ω) 对于ω=0、π、2 π都是奇对称,所以,H r (w) 以ω=0、π、2π为奇对称。
注意:
1) 在ω=0 和π处,有:
121
()()sin N j r n H e c n n ω
ω?==
=∑5.2.2 线性相移FIR DF 频率响应:Type III
29北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type IV
h(n) 奇对称,N 为偶数(恒群时延〕
1
1222
1
2()()()()sin()()n N
N j w j j r H e e
d n n
e H πω
θωωω??==?=∑0
1
2
N ?7h(n) 为奇对称,N 为偶数
n
h(n)其中
212322(),,,,...,
N N d n h n n ??
=?=????
2
11()d()sin -2N r n H n n ωω=??
??=????
????
∑1
2
2
N πθωω
??()=30
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.2.2 线性相移FIR DF 频率响应:Type IV
2
46
8
-0.4
-0.200.20.40.6h(n)
2
46
8
00.20.40.60.81d(n)
1
2
00.5
11.5frequency Unit:pi M a g n i t u d e
Magnitude Response
1
2
-30-20-10010frequency Unit:pi
P h a s e
Phase Response 0.5pi
Hr (w)
与d(n) 或h(n) 的取值无关,因此传输函数H(z) 在z = 1 处为零点。显然,这种类型不能用于实现低通滤波器。又有,所以这类滤波器适用于设计希尔伯特变换和微分器。
2)由于sin[(n-1/2)ω] 在ω=π处偶对称,在0、2π是奇对称,所以,H r (w) 以ω=π偶对称,0、2π为奇对称。
注意:
1) 在ω=0 处,有:
2
1
1
2()()sin()N
j r n H e d n n ωω==?=∑
31
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
一般形式:
()()()
j j r H e e H ωθωω=偶对称:1
2
()N ?θω=?
ω奇对称:
122
()N π?θω=
?ω(两个恒时延条件)
(一个恒时延条件)
( H(ω) 为ω的实函数)
5.2.2 线性相移FIR DF 频率响应:小结
32
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
一般的FIR DF 的零、极点:
1
1
111
()
()()()()()N N n
N N n N n n f z H z h n z
z
h n z z ????????=====
∑∑n 在z=0处,有一个(N-1)阶的极点,故滤波器稳定;o 其零点要求f(z)=0,根据代数上的理论,它为N-1阶多项式,应有N-1 个根,所以有N-1 个零点。如果h(n) 为实数值,其根肯定是共轭对称的。
5.2.3 线性相移FIR DF 零极点分布
33
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
1()()
h n h N n =±??1
,......,0?=N n 110
1()()()N N n n
n n H z h n z h N n z ????====±??∑∑令:m=N-1-n
)
()()()(1
)
1(1
)
1(1
0)
1(????=???=???±=±=±=∑∑z H z
z
m h z
z
m h z H N N m m
N N m m N 于是:
)
()(1)1(???±=z H z z H N 5.2.3 线性相移FIR DF 零极点分布
线性相移FIR DF 的零极点:
9如果z i 是H(z) 的零点,即H(z i ) = 0则H(z -1) =0,即z i -1亦为H(z) 的零点。
34
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
上面提到Z i 肯定是共轭的,故Z i * 亦必为其零点
于是零点有:
11
,*,,*
i i i i z z
z z 1-1
Z a1
Z a2
1/b b
①
②
③
④
5.2.3 线性相移FIR DF 零极点分布
总结:
1) 一般情况,i j i i
e r z ω=,有四个零点:
i
j i i e r z ω=i
j i i e r z ω?=*i j i i e r z ω???=11i
j i i e
r z ω11
*)(??=2) r=1,单位圆上的零点:
i
i j i j i e z e z ωω- ==(共轭对)
3) 位于实轴上的实数:b, 1/b (实轴上的倒数对)。4) z i =±1:单零点
()1
1**i i i i
z z z z ??===35北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
主题概述
1 -绪论
2 -离散时间信号和离散时间系统
3-离散傅里叶变换及其快速计算方法4-IIR 数字滤波器设计和实现
5–FIR 数字滤波器设计和实现
5.1) 概述
5.2) 线性相位FIR DF 约束条件和频率响应5.3) 窗函数法5.4) 频率取样法
5.5) FIR 数字滤波器的优化设计5.6) FIR 数字滤波器的实现结构5.7)本章小结5.8) 附录
6 –数字信号处理中的有限字长效应
36
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
思路:
理想数字滤波器
求出的FIR 数字滤波器
()()j jn d d
n H e h n e
∞
ω?ω
=?∞
=
∑1
()()N j jn n H e h n e ωω
??==∑()()d h n h n ????→窗函数
截短
要求:线性相位
尽可能降低逼近误差
5.3 FIR DF 窗口法(傅里叶级数法)
h d (n) 无限长,且非因果
h(n) 有限长,且因果
37
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
设所要求的DF 的频率响应是H d (e jw ),需要注意:它可能是低通、高通、带通和带阻FIR DF ,没有特指某种类型的数字滤波器。
不管是何种FIR DF ,它的频率响应是频域中的周期函数,周期为2π,所以它可以展开为傅氏级数形式:
()()j jn d d n H e h n e ∞
ω?ω
=?∞
=
∑
12()()j jn d d
h n H e e d πωω
?π
=ωπ∫5.3.1 窗口法:基本原理
式中h d (n) 是傅里叶系数,也是单位取样响应序列。由傅里叶级数理论可得:
38
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
因此,所要求的DF 的系统函数便可求得:
显然,H d (z) 是非因果的,且h d (n) 的持续时间为-∞~+∞,物理上不可实现。
我们可以采用逼近H d (e jw ) 的方法
n 首先把h d (n) 先截短为有限项,把h d (n) 截为2M+1项,得:
()()n
d d n H z h n z ∞
?=?∞
=
∑
5.3.1 窗口法:基本原理
1()()M
n
d n M
H z h n z ?=?=
∑
()()n
d d n H z h n z
∞
?=?∞
=
∑
39
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT o 然后把截短后的h d (n) 右移,使之变成因果性的序列。
令H(z) 等于H 1(z) 乘以z -M 得:
p 令h(n)= h d (n-M), n=0, 1, 2, ......, 2M ,则
210
()(()(()))M
M
M n M n
d n n d M
h H z z H z h n M n z z ??+?=?=?==
=∑
∑20
()()M
n
n H z h n z
?==∑20
()()M
j jn n H e h n e
ω
?ω
==∑频率响应z=e j ω
5.3.1 窗口法:基本原理
显然
H(z) 是物理可实现的
其冲激响应h(n) 的持续时间也是有限的
选择h d (n) = ±h d (N-1-n),保证H(z) 具有线性相位。
40
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
对h d (n) 的截短必然产生误差,即以|H(e jw )| 近似|H d (e jw )| 。定义逼近误差为均方误差:
2212|()()|j j d H e H e d π
ωω?πε=
?ωπ
∫而H d (e jw ) 可以展开为:
1
012()()cos()sin()j jn d d n n n n n a H e h n e a n b n ωωωω∞
∞
∞
=?∞
==?=
=
++∑∑∑式中:
020(); ()()[()()]
d n d d n d d a h a h n h n b j h n h n ==+?=??5.3.2 窗口法:性能分析
|H(e jw )| 对|H d (e jw )| 的逼近
41
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 因为|H(e jw )| 是对h d (n) 截短而产生的,假定:
011
2()cos()sin()
M
M
j n n n n A H e A n B n ω
===+ω+ω∑∑即当|n|>M 时,A n = 0, B n =0。
所以把上述两式代入逼近误差中,利用三角函数的正交性可得:
()()()()
2
2
2
2
2
01
1
1
2
M
M
n n n n n
n n n n M a A a A b B a
b ∞
===+?ε
=+?+?+
+∑∑∑由于上式中每一项都是正的,所以,只有当
0012,,,,,......, n n n n A a A a B b n M ====时
2min
ε=最小。
5.3.2 窗口法:性能分析
42
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
说明:
①当用|H(e jw )|≈|H d (e jw )| 时,要使ε2 =min ,|H(e jw )| 的单位取样响应h(n) 必须是|H d (e jw )|的傅里叶系数h d (n)。②有限项傅氏级数是在最小平方意义上对原信号的最佳逼近其逼近误差为:
22||
()
d n M h n ε>=
∑
M 越大,ε2愈小(因为h d (n) 值愈小)。
5.3.2 窗口法:性能分析
43
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
将h d (n) 截短:
0(),||(),d
h n n M
h n else
≤?=??相当于将h d (n ) 与一窗函数W R (n) 相乘,即
()()()
d R h n h n w n =5.3.2 窗口法:Gibbs 效应
其中
10,||(),R n M
w n else
≤?=?
?在一定意义上来看,窗函数决定了我们能够“看到”多少个原来的冲激响应,“窗”这个用词的含义也就在此。
44
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
窗函数的频谱:
5.3.2 窗口法:Gibbs 效应
jM -jM -j M
j jn jn R R -j n n -M
2M+12M+1-j j -j 2
2
2
-j
j -j 2
22e -e e W (e )w (n)e e 1-e (2M 1)N e e -e
sin sin 22sin sin
e
e -e 22
ωωω
∞
ω?ω?ωω
=?∞
=ωωω
ωωω=
==
??
+ωω????==
ωω??????
∑
∑=
此矩形窗谱为一钟形偶函数,在ω=+2π/N 之间为其主瓣,主瓣宽度Δω= 4π/N ,在主瓣两侧有无数幅度逐渐减小的旁瓣, 见图所示。
22sin(/)()sin(/)
j R N W e ωωω=
ω
2π/N -2π/N
主瓣
第1个旁瓣
第2个旁瓣
45
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
1212()()()()()[]j j j d R j j d
R H e H e W e H e W e d ω
ωωπ
θω?θ?π??=???π=
θπ∫ 截短,根据时域相乘映射为频域卷积,得:
5.3.2 窗口法:Gibbs 效应
为便于分析,我们假定|H d (e jw )| 是理想低通滤波器LPF 。
式中积分等于θ由–ωc 到ωc 区间内
W N [e j(w-θ)] 下的面积,随着
ω变化,不同正负、不同大小的旁瓣移入和移出积分区间,使得此面积发生变化,也即|H(e jw )| 的大小产生波动。
-w c
w c
θ
()[]
j R W e ωθ?46
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
θ
2N
π2N π?
W R (θ)
-ωc
ωc
ω
H d (θ)
0.5
0.5
0.0895
0.0895
0.0468
0.0468
卷积
?
5.3.2 窗口法:Gibbs 效应
-w c
w c
θ
π-π
H d (θ)
47北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
现在分析几个特殊频率点的滤波器性能:n ω= 0 时:由于一般情况下都满足ωc >> 2π/ N ,因此,H(0) 的值近似等于窗谱函数W R (e jw ) 与θ轴围出的整个面积。
5.3.2 窗口法:Gibbs 效应
θ
2N
π2N
π?
W R (θ)
①
-w c
0w c
θ
π-π
H d (θ)
①
c
c
j01
R
1H (e )W ()d 2ω?ω=θθ
π∫o ω=ωc 时:
此时窗谱主瓣一半在积分区间内一半在区间外,因此,窗谱曲线围出的面积,近似为ω=0 时所围面积的一半,即。
011122
()
()()?=?≈∫c
c
c j j R c H e H e
w d ωωωωθθπ111()H (0)2
c H ω≈-w c
w c θ
H d (θ)
w =w c
W R (w-θ)
②
48
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.2 窗口法:Gibbs 效应
p ω=ωc -2π/N 时,正肩峰
此时窗谱主瓣全部处于积分区间内,而其中一个最大负瓣刚好移出积分区间,这时得到最大值,形成正肩峰。之后,随着ω值的不断增大,H(e jw ) 的值迅速减小,此时进入滤波器过渡带。
q ω=ωc + 2π/N 时,负肩峰
此时窗谱主瓣刚好全部移出积分区间,而其中一个最大负瓣仍全部处于区间内,因此得到最小值,形成负肩峰。之后,随着ω值的继续增大,H(e jw ) 的值振荡并不断减小,形成滤波器阻带波动。
112121089502()().()c
c
c R
c H w
d H Max N
N
ωωππ
ωωθθπ?
?
=?
?==∫112008950().()c H H Min N
π
ω+
=?=-w c
w c
θ
H d (θ)
W R (w-θ)
2c w w N
π=?
③
-w c
w c
θ
H d (θ)
W R (w-θ)
2c w w N
π=+
④
49
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
n 理想滤波器的不连续点演化为过渡带
o 通带与阻带内出现起伏
p Gibbs 现象
过渡带:正负肩峰之间的频带。其宽度等于窗口频谱的主瓣宽度。
对于矩形窗W R (e jw ), 此宽度为4π/N 。
肩峰及波动:这是由窗函数的旁瓣引起的。
旁瓣越多,波动越快、越多。相对值越大,波动越厉害,肩峰越强。肩峰和波动与所选窗函数的形状有关,要改善阻带的衰减特性只能通过改变窗函数的形状。
?在对h d (n) 截短时,由于矩形窗函数的频谱有较大的旁瓣,这些旁瓣在与H d (e jw ) 卷积时产生了吉布斯现象。
?长度N 的改变只能改变w 坐标的比例及窗函数W R (e jw ) 的绝对大小,但不能改变肩峰和波动的相对大小(因为不能改变窗函数主瓣和旁瓣的相对比例,波动是由旁瓣引起的),即增加N ,只能使通、阻带内振荡加快,过渡带减小,但相对振荡幅度却不减小。
加窗处理对理想矩形频率响应的影响:注意:过渡带宽度与窗的宽度N 有关,随之增减而变化。
阻带最小衰减(与旁瓣的相对幅度有关)只由窗函数决定,与N 无关。
5.3.2 窗口法:Gibbs 效应
50
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
Gibbs 现象;
5.3.2 窗口法:Gibbs 效应
51
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
设计FIR DF 时,窗函数不仅可以影响过渡带宽度,还能影响肩峰和波动的大小,因此,选择窗函数应使其频谱:
n 主瓣宽度尽量小,以使过渡带尽量陡。
o 旁瓣相对于主瓣越小越好,这样可使肩峰和波动减小,即能量尽可能集中于主瓣内。
对于窗函数,这两个要求是相互矛盾的,要根据需要进行折衷的选择,5.3.3 窗口法:常用窗函数
w
020lg|W(w)/W(0)|
B
3dB
A
D (dB/Oct)
为了定量地比较各种窗函数的性能,给出三个频域指标:
n 3db 带宽B ,单位为(最大可
能的频率分辨力)
o 最大旁瓣峰值A(dB),A 越小,由旁瓣引起的谱失真越小
p 旁瓣谱峰渐进衰减速度D (dB/oct ) 一个好的窗口,应该有最小的B 、A 及最大的D 。
2N
π
ωΔ=52
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
以下介绍的窗函数
均为偶对称函数,都具有线性相位特性。
设窗的宽度为N ,窗函数的对称中心点在(N-1)/2处。因此,均为因果函数。
矩形窗
最简单的窗函数,从阻带衰减的角度看,其性能最差。
它的频率响应函数为:
5.3.3 窗口法:基本窗函数_矩形窗
101
0,(), n N w n else
≤≤??=?
?1
2
22sin ()sin N j N W e
ωωωω????????=??????
22sin ()sin r N W ωωω????
??=
??????
振幅响应
53
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
为了对过渡带和阻带衰减进行精确分析,对窗振幅响应进行连续积分(或累积振幅响应),即
矩形窗函数w(n) 以及它的振幅响应、累积振幅响应如下图所示。
5.3.3 窗口法:基本窗函数_矩形窗
1121222sin ()()sin w wc w wc j r N H e Wr d d , N ωππλλλλλππ++????????==>>??????∫∫ 0
22450
1
Rectangular Window : N=45n
w (n )
-1
010
45
Amplitude Response
frequency in pi units W r
-1
01
4013
0Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1 1
50
21
Accumulated Amplitude Response
frequency in pi units D e c i b e l s
Width=(1.8)*pi/N
性能指标
3dB 带宽B=0.89Δω 最大旁瓣峰值
A= -13dB
旁瓣谱峰渐进衰减速度
D=-6dB/oct 在Matlab 中,实现矩形窗的函数为w = boxcar(n)。
54
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
振幅响应
在ω= ω1处具有第一个零点:
因而主瓣的宽度为2Δω,所以过渡带宽也近似为2Δω。
大约在w = 3π/N 处,出现第一个旁瓣(即主旁瓣),其幅度为:
将它与主瓣振幅N 比较,则最大旁瓣峰值A(dB) 为A= -13db 。
累积振幅响应
第一个旁瓣为21dB ,这个21dB 的阻带衰减与窗长度N 无关。 根据最小阻带衰减,可以精确地计算出过渡带宽为: 它大约是近似带宽的一半。1122
N
N
ωπ
πωωΔ==
=或33221332sin()() ,
sin()N Wr M N N π
πππ=≈≥时
18.s p N
π
ωω?=
5.3.3 窗口法:基本窗函数_矩形窗
55北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
三角窗(或巴特利特Bartlett 窗)
由于矩形窗从0 到1 (或1 到0 )有一个突变的过渡带,这造成了吉布斯现象。Bartlett 提出了一种逐渐过渡的三角窗形式,它是两个矩形窗的卷积。
B=1.28Δω, A=-27dB, D= -12dB/oct ,近似过渡带宽8π/N ,精确过渡带宽6.1π/N ,最小阻带衰减25dB 。与矩形窗来比较,阻带衰减性能有所改善,但代价是过渡带的加宽。
5.3.3 窗口法:基本窗函数_三角窗
2012
1
2,,,...,()(),,...,n
N n N w n N w N n n N ?=??=?
??=???2
12
242()sin ()sin N j N W e N ωωωω????
????????=??
????
???
????
?56
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
2245
1
Triangular Window : N=45n
w (n )
-1
1
022
Amplitude Response
frequency in pi units
W r
-1
01
6027
Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
6025
Accumulated Amplitude Response
frequency in pi units
D e c i b e l s
Width=(6.1)*pi/N
在Matlab 中,函数bartlett(n)和triang(n)用来计算相似的三角窗,但它
们有两个重要的区别:bartlett 函数返回的序列两端总是0,因此,对于奇数n ,语句bartlett(n+2) 的中间部分等于triang(n);对于偶数n ,bartlett 仍然是两个矩形序列的卷积,但n 为偶数时的三角窗没有标准定义。
5.3.3 窗口法:基本窗函数_三角窗
57
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 余弦窗
B=1.2Δω, A=-23dB, D=-12dB/oct 。近似过渡带宽8π/N ,精确过渡带宽6.5π/N ,最小阻带衰减34dB 。
5.3.3 窗口法:基本窗函数_余弦窗
0121
1()sin ,,,...,n w n n N N π??
==??????
10122()cos ,...,,,...,
n N N w n n N π??
==???????
2
1
2
11()N
j W e
u u N N ωππωωω???????=?++??????
???
?????222sin ()sin j
N u e ωωωω??????=??????
或
其中
频率响应58
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.3 窗口法:基本窗函数_余弦窗
22450
1
Triangular W indow : N=45n
w (n )
-1
1
027.9994
Amplitude Response
frequency in pi units
W r
-1
01
60
23
0Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
6034
Accumulated Amplitude Response
frequency in pi units
D e c i b e l s
Width=(6.5)*pi/N
59
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
升余弦窗函数
汉宁窗、汉明窗、布莱克曼窗都是升余弦窗的特例。它们都是频率为0~ 2π/(N-1) 和4π/(N-1) 的余弦序列的组合。升余弦窗的频率特性比矩形窗有很大改善。
其中A 、B 、C 为常数。
n 当A = 0.5,B=0.5,C=0 时,为汉宁(Hanning)窗。Matlab 中,w = hanning(n)
o 当A = 0.54,B=0.46,C=0 时,为汉明(Hamming) 窗。Matlab 中,w = hamming(n)
p 当A = 0.42,B=0.5,C=0.08 时,为布莱克曼窗。Matlab 中,w = blackman(n)
5.3.3 窗口法:升余弦窗函数
2()cos()cos()
w n A B n C n =?+60
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
-22
022
01
Window Function
n
w (n )
Rectangular
Triangle -.-Hanning -o-
Hamming -x-
Blackman -*-5.3.3 窗口法:升余弦窗函数
61
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
-1
1
Amplitude Response
frequency in pi units
W r
W(ω)0.5u(ω)
0.25u(ω-2π/(N-1))
0.25u(ω+2π/(N-1))
Hanning 窗(升余弦窗)
5.3.3 窗口法:升余弦窗函数_汉宁窗
2205050121
11()sin ..cos ,,,...,n n w n n N N N ππ????
==?=???????????
220502511().().W u u u N N ππωωωω?????
?=+?++?????????
?????B=1.44Δω, A=-32db,
D=-18db/oct ,近似过渡带宽8π/N ,精确过渡带宽 6.2π/N ,最小阻带衰减44dB 。与矩形窗来比,最小阻带衰减性能明显提高,但过渡带也明显增大。
62
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.3 窗口法:升余弦窗函数_汉宁窗
22450
1
Hanning Window : N=45n
w (n )
-1
1
023
Amplitude Response
frequency in pi units
W r
-1
01
6032
Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
6044
Accumulated Amplitude Response
frequency in pi units
D e c i b e l s
Width=(6.2)*pi/N
63
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
Hamming 窗(改进的升余弦窗)
5.3.3 窗口法:升余弦窗函数_汉明窗
20540460121
1()..cos ,,,...,n w n n N N π??
=?=??????2205402302311().()..W u u u N N ππωωωω???
?=+?++?????????
?B=1.3Δω, A=-43dB, D=-6dB/oct ,近似过渡带宽
8π/N ,精确过渡带宽6.6π/N ,最小阻带衰减53dB 。通过这一系数调整,使能量的99.963% 都集中在了窗谱的主瓣内。
64
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.3 窗口法:升余弦窗函数_汉明窗
-22
0220
1
Hamming Window : N=45n
w (n )
22
45
023.84
Amplitude Response
frequency in pi units
W r
-1
01
6043
Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
7053
Accumulated Amplitude Response
frequency in pi units
D e c i b e l s
Width=(6.6)*pi/N
65北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
Blackman 窗(二阶升余弦窗)
5.3.3 窗口法:升余弦窗函数_布莱克曼窗
22042050082012111()..cos .cos ,,,...,n w n n n N N N ππ????
=?+=???
????????
22042025114400411().()..W u u u N N u u N N ππωωωωππωω?????
?=+?++?????????
??????????
?+?++?????????
?????B=1.68Δω, A=-58db, D=-18db/oct ,近似过渡带宽12π/N ,
精确过渡带宽11π/N ,最小阻带衰减74dB 。通过增加余弦的二次谐波分量,能够进一步抑制旁瓣,但主瓣宽度却比矩形窗谱的主瓣宽度大三倍。
66
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.3 窗口法:升余弦窗函数_布莱克曼窗
-22
0220
1
Blackman W indow : N=45n
w (n )
22
45
018.48Amplitude Response
frequency in pi units
W r
-101
58
Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
74
0Accumulated Amplitude Response frequency in pi units
D e c i b e l s
Width=(11)*pi/N
比较以上窗函数,可以看到,矩形窗函数具有最窄的主瓣B,但也有最大的旁瓣峰值A 和最慢的衰减速度D。
汉宁窗主瓣稍宽,但有着较小的旁瓣和较大的衰减速度,因而被认为是较好的窗口。
67北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.3 窗口法:凯瑟窗函数
凯瑟(Kaiser )窗
上面讨论的几种窗函数以牺牲主瓣宽度,换取旁瓣抑制; Kaiser 窗全面反映了这种主瓣和旁瓣衰减之间的互换关系;
定义了一组可调的由零阶贝塞尔Bessel 函数构成的窗函数;
通过调整参数β可以在主瓣宽度和旁瓣衰减之间自由选择它们的比重。从而实现以同一种窗类型来满足不同窗性能需求的目的。
Kaiser 窗函数由J.F. Kaiser 提出,由下式给出:
其中I 0是修正过的零阶贝塞尔Bessel 函数 β是用来调整窗形状的参数,β依赖于参数N 。
()
2001211(),I w n 0n N 1I n N ββ?
??
????
??
?
?=≤≤?????????68
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
对于相同的N ,Kaiser 窗可以提供不同的过渡带宽,这是其他窗函数做不到的。
通过调整参数β,就可以方便地完成对过渡带宽度和阻带衰减的调整。参数β越高,其频谱的旁瓣越小,但主瓣宽度也相应增加。
x
I 0(x)
I 0(β)
1
β
零阶贝塞尔函数
Kaiser 窗函数依参数β而变化
Matlab 中,函数w = kaiser(n, beta) 实现Kaiser 窗
5.3.3 窗口法:凯瑟窗函数
69
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
下面图固定β,当窗的长度变化时,相应的旁瓣的高度保持不变。
5.3.3 窗口法:凯瑟窗函数
70
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
β= 5.658,则过渡带宽等于7.8pi/N ,最小阻带衰减为60dB ,如下图所示:
5.3.3 窗口法:凯瑟窗函数
-22
0220
1
Kaiser Window : N=45
n
w (n )
22
45
022.6383Amplitude Response
frequency in pi units
W r
-101
42
0Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
60
0Accumulated Amplitude Response frequency in pi units
D e c i b e l s
Width=(7.8)*pi/N
71北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
凯瑟窗的计算
由于Bessel 函数的复杂性,这种窗的设计公式很难推导,为此,Kaiser 提出了经验公式。
给定ωp 、ωs 、R p 和A s ,参数β定义如下:
对于过渡带宽△ω= ωs -ωp (rad/s),滤波器阶数为
需要强调的是:阶数为N 的滤波器大致能满足要求,但最后的结果还必须要演算以便证明这一点。
在Matlab 中,函数w = kaiser(n, beta) 实现Kaiser 窗。
5.3.3 窗口法:凯瑟窗函数
040110287058422100788621..(.),.().(),s s s s s s
A A 50A A 21A 50
0, A 2β?>=?+?≤≤<1??
???795
1
2285..s A N ω
?=
+Δ795
1
1436..s A N f
?=
+Δ72
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
切比雪夫(Chebyshev )窗
在给定旁瓣高度下,Chebyshev 窗的主瓣宽度最小,具有等波动性,也就是说,其所有的旁瓣都具有相等的高度。 Matlab 函数w = chebwin(n,r)
以窗长度n 和旁瓣高度rdB 为参数计算切比雪夫窗。
Chebyshev 仅对奇数长度的窗有定义,若n 为偶数,函数w = chebwin(n,r) 先将它加1,然后设计长为n+1 的切比雪夫窗。 其傅里叶变换的旁瓣幅度低于主瓣r dB 。
其它窗函数
Papoulis 窗、Parzen 窗、Poisson 窗、Cauchy 窗、Gaussian
窗、Bartlett-Hann 、Blackman-Harris 、Nuttall's Blackman-Harris 、Bohman 、Flat Top window 、Hann 、Parzen (de la Valle-Poussin)、Tapered cosine 等。
Matlab 窗设计和分析工具(WinTool)具有GUI 界面,可以用来设计和分析窗函数,其用法:
>> wintool
5.3.3 窗口法:切比雪夫窗函数
73
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 5.3.3 窗口法:切比雪夫窗函数
-22
0220
1
Hanning Window : N=45n
w (n )
-1
1
025.9905
Amplitude Response
frequency in pi units
W r
-1
01
6040 0
Amplitude Response in dB
frequency in pi units
D e c i b e l s
-1
1
7055
Accumulated Amplitude Response
frequency in pi units
D e c i b e l s
74
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
-80
10π/N
10π/N
-57
凯塞(β=7.865)
-7411π/N 12π/N -57布莱克曼窗(二阶升余弦窗)-536.6π/N 8π/N -41汉明窗
(改进升余弦窗)-446.2π/N 8π/N -31汉宁窗(升余弦窗)-211.8π/N 4π/N -13矩形窗加窗后滤波器阻带最小衰减(dB )
加窗后滤波器过渡带宽(△ω)窗函数主瓣宽度旁瓣峰值衰减(dB )窗函数5.3.3 窗口法:常用窗函数的性能指标
75
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
FIR DF 窗口法设计步骤
性能要求?H d (e j ω)
n 把H d (e j ω) 展成傅里叶级数,得到h d (n);o 把h d (n) 自然截短到所需的长度N=2M+1;p 将截短后h d (n) 右移M 个取样间隔,得h(n);q 将h(n) 乘以合适的窗函数,即得所需的滤波器的冲激响应,这时窗函数以n = M 对称(当然窗函数也可直接加在h d (n) 上,这时窗函数以原点为对称);r 利用h(n),既可用硬件构成滤波器的系统函数H(z),也可直接用计算机软件实现滤波。
5.3.4 窗口法:设计步骤
76
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
数字低通滤波器的设计
下面通过一个例题来阐述低通滤波器所涉及的一些问题。
例5.1一个理想低通数字滤波器的频率响应如图所示,为:
假定w c =0.25π,分别取N =11、21、31,观察加窗后对滤波器幅频特性的影响。
5.3.4 窗口法:设计步骤_ 数字低通
100,||(),c
j d c w H e w ωωωπ
≤≤?=?
<≤?|()|
j d H e ω0
c ωc ω?π
?π077
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
解:由于H d (e j ω)是一个实周期函数,把它展成为付氏级数:
()()j jn d d n H e h n e ω
ω
∞
?=?∞
=
∑()
()111221
2()()sin c c c c c c w w jw jn jn d d w w jnw jnw c h n H e e d e d nw e e j n n ω
ωωωππ
ππ
???=
?=?=
?=∫∫5.3.4 窗口法:设计步骤_ 数字低通
式中
将h d (n) 截短为N=2M+1,并将截短后的h d (n) 移位,得:
sin()()()()c d n M w h n h n M n M π
?=?=
?78
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
然后乘以窗函数w R (n),得到最后得h(n)。对于w c = 0.25π,由上式得:
025sin[().]
()()()d n M h n h n M n M ππ
?×=?=
?(当w c =π,就得到一个全通滤波器)
当N=11 时,M= 5,求得
h(0) = h(10) = -0.045,h(1) = h(9) = 0,h(2) = h(8) = 0.075,h(3) = h(7) = 0.1592,
h(4) = h(7) = 0.2251,h(5) = 0.25。
当N=11 时,乘以汉明窗:
20540460121
1()..cos ,,,...,n w n n N N π??
=?=??????
79
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
-5
05
-0.1
00.10.20.3hd(n)
510
-0.1
00.10.20.3h(n)
510
0.20.40.60.81w(n)
510
-0.1
00.1
0.20.3
h(n)
80
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
o 当N 取不同值时,H(e jw ) 都不同程度上近似于H d (e jw )。
o N 过小时,通带过窄,且阻带内波纹较大
o 当N 增加时,通频带接近于0.25pi,阻带内波纹减小,但在通带内出现了波纹,随着N 的继续增加,这些波纹并不能消失。
o 由图中还可以看到,使用汉明窗后,通带内的振荡基本消失,阻带内的纹波也大大减小,从这一点上来说,滤波器的性能得到了改善,但是,这是以过渡带的加宽为代价的。
矩形窗汉明窗
5.3.4 窗口法:设计步骤_ 数字低通
81
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
适用范围
对于能用解析式表达,且傅里叶级数的系数容易求解的滤波器:
此时,窗函数法是设计FIR DF 较为方便的一个方法;
如果h d (n) 不易求,则使用该方法较为困难。
窗函数
用窗口法设计FIR DF ,一个重要问题是选用窗函数w(n) 及决定截短的长度N 。
窗函数的选择:阻带衰减指标
满足设计给定的阻带衰减和其它滤波器性能要求; 能量尽量集中于主瓣内; 个人的经验及喜好有关。
5.3.4 窗口法:设计步骤_ 数字低通
1
2()()j jn d d
h n H
e e
d π
ωω
πω
π
?
=
∫sin() ()() c
c n w n n h n w n τττπ
τπ
??≠???=?
?=??82
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
窗函数长度N 的选择:过渡带宽指标
采用试验方法,即逐渐增大M ,直至H(e jw ) 在通带和阻带内部达
到指标要求。
若对|H d (e jw )|的过渡带提出了具体要求,因为FIR DF 的过渡带
等于窗函数的主瓣宽度,那么通过查表计算N :
例如5.2?6,5.8?6。一般选择N 为奇数。
注意:根据所要设计线性相位FIR DF 类型来决定最终N 取奇数还是偶数。
截止频率w c 的确定:
截止频率w c 对应于明确的0.5 增益点,而不再标志某个增益点。 对于非理想滤波器,其截止频率w c 不采用通带边缘频率w p 或阻
带边缘频率w s ,而使用过渡带的中点(即通带边缘和阻带边缘之间的中点。因此,窗函数法不能精确地确定其通带和阻带的边缘频率):
5.3.4 窗口法:设计步骤_ 数字低通
N ??=???>????
相应的窗函数的主瓣宽度或精确过渡带上取整
要求的滤波器过渡带222
s p p s c p p ωωωωωωω?+=+
=+=
过渡带宽度
83
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
频率
增益
所要求的通带边缘频率ωp
所要求的阻带边缘频率ωs
设计中所用的截止频率ωc
过渡带宽度
0.5
84
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
例5.2
根据下列指标设计低通滤波器
通带边缘频率f p =2kHz 阻带边缘频率f stop =3kHz 阻带衰减40dB 取样频率f s =10kHz
解:
(1) 求对应的理想数字频率:
过渡带宽= 3kHz –2kHz = 1kHz 。转换为数字频率过渡带:截止频率:
数字截止频率:2120210
.c s f f πωπ
π×Δ===kHz
f f f stop
p c 5.22
=+=
205.c
c
s
f f ωπ
π==5.3.4 窗口法:设计步骤_ 数字低通
频率
|H(e jw )|
所要求的通带边缘频率ωp
所要求的阻带边缘频率ωs
设计中所用的截止频率ωc
过渡带宽度
-6dB 0dB -40dB
2kHz
3kHz
85
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
(2) 设理想线性相位滤波器为:由此可得脉冲响应:
??
?≤=其它
,0||,e )(c
-j e H j d ωωωτω[]05-j j n sin (-)1sin(-).() e e d 2(-)(-)c d n w n h n n n π
ωτωπ
ττπ
ωπ
τπτπ
=
=
=
∫
-205051()..cos ,n w n n 0,1,2,N 1
N π??
=?=??????
)
()51(]
5.0)51sin[()(n w n n n h ×?×?=
π
π62623102...N ππ
ωπ=
==Δ(4) 由过渡带宽确定窗口长度:
(3) 由阻带衰减确定窗函数:因为阻带衰减40dB ,通过查表知道,可以
选择Hanning 窗:
5.3.4 窗口法:设计步骤_ 数字低通
τ=(N-1)/2=15
则此滤波器的脉冲响应为:
86
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
h d (n) =
-0.0212 0.0000 0.0245 -0.0000 -0.0289 0.0000 0.0354 -0.0000 -0.0455 0.0000 0.0637 -0.0000 -0.1061 0.0000 0.3183 0.5 0.3183 0.0000 -0.1061 -0.0000 0.0637 0.0000 -0.0455 -0.0000 0.0354 0.0000 -0.0289 -0.0000 0.0245 0.0000 -0.0212w(n) =
0 0.0109 0.0432 0.0955 0.1654 0.2500 0.3455 0.4477 0.5523 0.6545 0.7500 0.8346 0.9045 0.9568 0.9891 1.0000 0.9891 0.9568 0.9045 0.8346 0.7500 0.6545 0.5523 0.44770.3455 0.2500 0.1654 0.0955 0.0432 0.01090
5.3.4 窗口法:设计步骤_ 数字低通
计算得:
87
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
h(n) =
0 0.0000 0.0011 -0.0000 -0.0048 0.0000 0.0122 -0.0000 -0.0251 0.0000 0.0477 -0.0000-0.0960 0.0000 0.3148 0.5000 0.3148 0.0000 -0.0960 -0.0000 0.0477 0.0000 -0.0251 -0.0000 0.0122 0.0000 -0.0048 -0.0000 0.0011 0.00000
由h(n) 可以得到H(z)、H(e jw ) 和差分方程,其中差分方程为:y(n) = 0.0011x(n-2) -0.0048x(n-4) +0.0122x(n-6)
-0.0251x(n-8) + 0.0477x(n-10) -0.0960x(n-12) + 0.3148x(n-14) + 0.5x(n-15) + 0.3148x(n-16) -0.0960x(n-18) + 0.0477x(n-20) -0.0251x(n-22) + 0.0122x(n-24) -0.0048x(n-26) + 0.0011x(n-28)
5.3.4 窗口法:设计步骤_ 数字低通
88
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
89
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
90
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
例5.3一段乐曲中夹杂着高频噪声,严重影响收听质量。
下图给出了噪声污染后的信号及其频谱。系统的取样频率为16kHz 。设计一个滤波器来提高声音质量。
5.3.4 窗口法:设计步骤_ 数字低通
1
2
34
5
6
x 10
4
-1-0.500.51x(n)
1
2
3
45
6
7
8
02004006008001000X(k)
kHz
05-03-noisy.wav
91
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 此低通滤波器的参数总结如下:
取样频率:16kHz 通带边缘频率:2kHz 过渡带宽:545Hz 截止频率:数字截止频率:
解:从图中频谱可以看出,噪声从2kHz 的频率点开始占据主导地
位,因为这个频率很低,消除噪声的同时也会损失一部分音乐信息。
可用通带边缘频率为2kHz 的低通滤波器对此段音乐进行滤波。由于没有特殊的阻带衰减要求,任何具有合理特性的窗函数均可,在此选择Hamming 窗。
本例中,对过渡带宽也没有特殊要求。为了得到合理的陡峭滚降,选择N = 101,则过渡带宽为:
Hz
545216000
1016.6=×ππ200054522272/.c f kHz
=+=2272
22028416..c c s
f f ωπππ
===5.3.4 窗口法:设计步骤_ 数字低通
92
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
窗函数为Hamming 窗
1
N ,0,1,2,3,n N n n w ?=??
?
?????= 12cos 46.054.0)(π)
()50(]
284.0)50sin[()(n w n n n h ×?×?=
π
π-1
-0.5
0.5
1 1.5
-1
-0.8-0.6-0.4-0.200.20.40.60.81100
Real Part
I m a g n a r y P a r t
Pole-Zero Plot
其设计过程同前例。
得到的脉冲响应、幅频特性和零极点图如图所示。
窗口长度:N=101
则滤波器的脉冲响应为:
5.3.4 窗口法:设计步骤_ 数字低通
93北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
20
40
60
80
100
-0.1
00.10.20.3h(n)
1
2
3
45
6
7
8
-100
-80-60-40-200
kHz
M a g n i t u d e
Magnitude Response 5.3.4 窗口法:设计步骤_ 数字低通
94
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_ 数字低通
1
2
34
5
6
x 10
4
-1-0.500.51x(n)
1
2
3
45
6
7
8
0500100015002000X(k)
kHz
05-03-filtered.wav
伴有噪声的音乐经过上述滤波器滤除噪声后的信号及其频谱如上图所示。滤波后的信号几乎没有噪声,但歌曲听起来有点压抑,这是因为歌曲中高频分量也随噪声一起被滤除的缘故。
可以想象,如果噪声存在于所有的频率分量上,则不可能在不严重降低信号质量的情况下滤除噪声。
95
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
以上设计的是数字低通滤波器,若希望设计数字高通、带通和带阻滤波器,只需要改变付氏级数系数中积分的上、下限即可。 数字高通滤波器
5.3.4 窗口法:设计步骤_数字高通
|()|
j d H e ω
ω
c
ωc ω?ππ
?0令其幅频特性为:
0()j M j c d c
e
||H e 0||ωω
ωωπ
ωω??≤≤=?
≤(频域为e -jM ω,表示时域右移M 位)
96
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_数字高通
则
()()
1122121
2()()()()()()()()cos()cos(cos()co [sin()()sin()sin()sin()s())]
c c
c
c j n M j n M
d j n M j n M j n M j n M c c c c h n
e d e
d e e e e j n M j n M j n M j n M j n M n M n M j n M n M n M ωπ
π
ωωω
ωωππωω
π
π
πωπ
ππωωωππ??
????????=
+
??
=
?+??
??=????+?++????????∫∫求得
s sin[()]sin in[()]([()]
()()(sin[()))](()
)c d ap lp c n M n M h n n M h n n h M n M n M n M n πωπππ
ωπ
????=
?==?????
97
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_数字高通从这个结果可以看出:一个高通滤波器相当于用一个
全通滤波器(即ωc =π)减去一个低通滤波器。
传输函数:脉冲响应:()()()hp ap lp H z H z H z =?()()()
hp ap lp h n h n h n =?H ap (z)
H lp (z)
X(z)
Y(z)
98
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
数字带通滤波器
5.3.4 窗口法:设计步骤_数字带通
令其幅频特性为:
0() j M j l h
d e
||H e others
ωω
ωωω??≤≤=?
?(频域为e -jM ω,表示时域右移M 位)
|()|
j d H e ωω
l ωl ω?π
π?0h ω?h ω则
1122()()()l
h
h
l
j n M j n M d h n e d e d ωωωωω
ω
ωω
π
π
????
=
+
∫∫99
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_数字带通
求得
sin[()]sin[()]
()()h l d n M n M h n n M ωωπ
???=
?从这个结果可以看出:一个带通滤波器相当于两个低通滤波器相减,其中一个截止频率为ωh ,另一个为ωl 。
传输函数:脉冲函数:
()()()bp lph lpl H z H z H z =?()()()
bp lph lpl h n h n h n =?H lph (z)
H lpl (z)
X(z)
Y(z)
100
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_数字带通
或者一个带通滤波器相当于一个低通滤波器和一个高通滤波器相乘,即先经过一个LP DF ,再经过一个HP DF 。
传输函数:脉冲函数:
()()()
bp lp hp H z H z H z =()()()
bp lp hp h n h n h n =?(频域相乘,时域卷积)
H lp (z)H hp (z)
X(z)
Y(z)
101
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 数字带阻滤波器
5.3.4 窗口法:设计步骤_数字带阻
令其幅频特性为:
0,||,
||() j M
j l h h
d e || H e others
ωωωωωωππωω??≤≤≤?≤≤?=?
?(频域为e -jM ω,表示时域右移M 位)
则
111222()()()()h
l l
h
w w j n M w j n M w j n M w d w w h n e dw e dw e dw
ππ
π
π
π
?????
?=
+
+
∫∫
∫
|()|
j d H e ω
ω
l ωl ω?π
π?0h ω?h ω102
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤_数字带阻
求得
sin[()]sin[()]sin[()]
()()l h d n M w n M n M w h n n M ππ
?+???=
?从这个结果可以看出:一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器,低通滤波器的截止频率为ωl ,高通在ωh 。
传输函数:脉冲函数:
()()()bp lp hp H z H z H z =+()()()
bp lp hp h n h n h n =+H lp (z)
H hp (z)
+
X(z)
Y(z)
103
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT 例5.4给取样频率为22kHz 的系统设计一个FIR 带通滤
波器,中心频率为4kHz ,通带边缘在 3.5kHz 和4.5kHz ,过渡带宽为500Hz ,阻带衰减50dB 。5.3.4 窗口法:设计步骤
|()|
j d H e ω()
f kHz 3.5
h
f 4.5
l f 4
解:
过渡带宽:500Hz ,转换为数字频率为:
20520045522..c s
f f πωπ
π
×Δ===截止频率:
0523502532505245025475./..../...l lc h hc f f kHz f f kHz
=?=?==+=+=数字截止频率:
203045.l
l s
f f ωπ
π==204318.h
h s
f f ωπ
π==104
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤
脉冲响应:
sin[()]sin[()]
()()h l d n M n M h n n M ωωπ
???=
?窗函数:因为阻带衰减50dB ,可以选择Hamming 窗,即
20540461()..cos ,n w n n 0,1,2,N 1
N π??
=?=??????
并且窗口长度为:
6600455.146.N ππ??
=?
???
=,一般N 取奇数,因此N =147,M=73
则此滤波器的脉冲响应为:
7304318730304573sin[().]sin[().]
()()
()n n h n w n n πππ
?×??×=
×?105
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤
计算得:
h(n) = -0.0006 -0.0000 0.0007 0.0006 -0.0002 -0.0007 -0.0003
0.0003 0.0004 0.0001 -0.0000 0.0002 -0.0001 -0.0008-0.0008 0.0007 0.0019 0.0008 -0.0017 -0.0024 -0.00010.0023 0.0018 -0.0006 -0.0014 -0.0004 -0.0001 -0.0011-0.0005 0.0027 0.0039 -0.0007 -0.0063 -0.0048 0.0035
…………………………………………………………………….
-0.0097 0.0930 0.0771 -0.0210 -0.0763 -0.0368 0.02990.0439 0.0086 -0.0184 -0.0124 0.0004 -0.0019 -0.00580.0053 0.0175 0.0084 -0.0145 -0.0212 -0.0022 0.01780.0149 -0.0036 -0.0129 -0.0057 0.0037 0.0039 0.00030.0011 0.0033 -0.0003 -0.0064 -0.0054 0.0031 0.00830.0035 -0.0048 -0.0063 -0.0007 0.0039 0.0027 -0.0005-0.0011 -0.0001 -0.0004 -0.0014 -0.0006 0.0018 0.0023-0.0001 -0.0024 -0.0017 0.0008 0.0019 0.0007 -0.0008-0.0008 -0.0001 0.0002 -0.0000 0.0001 0.0004 0.0003-0.0003 -0.0007 -0.0002 0.0006 0.0007 -0.0000 -0.0006
106
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤
02040
6080100
120140
-0.1
-0.0500.050.10.15h(n)
2 3.544.5
681012
-100
-80-60-40-200
kHz
M a g n i t u d e
Magnitude Response
107
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
5.3.4 窗口法:设计步骤
-1
-0.5
00.5
1
-1
-0.8-0.6-0.4-0.200.20.40.60.81146
Real Part
I m a g i n a r y P a r t
Pole-Zero Plot
108
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
FIR DF Matlab 设计函数
b=fir1(n,wn,options),单带FIR 滤波器 b = fir2(n,f,m,options),多带FIR 滤波器
两者可设计设计低通、高通、带通、带阻和通用多带FIR 滤波器
Fir1 具有以下多种形式:
b = fir1(n,Wn)
b = fir1(n,Wn,'ftype')b = fir1(n,Wn,window)
b = fir1(n,Wn,'ftype',window)b = fir1(...,'noscale')
参数
向量b 是n 阶FIR 滤波器的系数
截止频率Wn 是从0 到1 之间的数,1 对应着奈氏频率。 对于高通滤波器,ftype 为‘high’
5.3.5 窗口法:Matlab 实现
109北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
对于带通或带阻滤波器,Wn 为包含通频带边带频率的一个二元素向量,‘ftype’为‘stop’
参量‘window’表示所采用的窗函数类型。window 的长度必须为n+1,若window 缺省,则fir1 使用汉明窗。
注意因为奇数阶的II 型滤波器(h(n) 为偶对称,长度N 为偶数)在高频段的频率响应为零,所以fir1 函数在高通和带阻情况下不设计II 型滤波器,因此,如果n 为奇数时,fir1 将阶次加1 并返回I 型滤波器。
函数b = fir2(n,f,m)
也可设计加窗FIR 滤波器,但它针对任意形状的分段(piece-wise )线性频率响应。
向量f 由从0 到1 的频率点组成,其中1 表示奈氏频率,第一个点必须是0,最后一个点必须是1,频率点必须是递增的。 m 是对应于频率点f 处的期望的频率幅值响应。 f 和m 的长度必须相等。
5.3.5 窗口法:Matlab 实现
110
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
Matlab 频率响应函数
在Matlab 中提供了一个freqz 函数,利用这个函数开发一个新的函数freqz_m ,它给出了绝对的和相对的dB 值幅度响应、相位响应以及群延时响应。
5.3.5 窗口法:Matlab 实现
function [db,mag,pha,grd,w] = freqz_m(b,a);% Modified version of freqz subroutine % ------------------------------------% [db,mag,pha,grd,w] = freqz_m(b,a);
% db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians % grd = Group delay over 0 to pi radians
% w = 501 frequency samples between 0 to pi radians % b = numerator polynomial of H(z) (for FIR: b=h)% a = denominator polynomial of H(z) (for FIR: a=[1])%
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))';mag = abs(H);
db = 20*log10((mag+eps)/max(mag));pha = angle(H);
grd = grpdelay(b,a,w);
111
北京邮电大学电信工程学院多媒体通信中心门爱东D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
例5.5
根据下列技术指标,设计一个数字FIR 低通滤波器wp = 0.2π, Rp = 0.25dB ws = 0.3π,As = 50dB
采用汉明窗,确定脉冲响应,并给出所设计的滤波器的频率响应图。解:在设计中,没有使用通带波动值Rp=0.25dB ,但必须检查设计的实际波动,验证它是否确实在给定容限内。
5.3.5 窗口法:Matlab 实现
% Lowpass filter design -Hamming window wp = 0.2*pi; ws = 0.3*pi;tr_width = ws -wp
N = ceil(6.6*pi/tr_width) + 1n=[0:1:N];
wc = (ws+wp)/2h = fir1(N,wc/pi);
[db,mag,pha,grd,w] = freqz_m(h,[1]);delta_w = 2*pi/1000;
Rp = -(min(db(1:1:wp/delta_w+1))) % Passband Ripple
As = -round(max(db(ws/delta_w+1:1:501))) % Min Stopband attenuation
112
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
运行结果如下:
tr_width = 0.3142 (过渡带宽)
N = 67 (滤波器的阶数,长度为68)wc = 0.7854 (理想LPF 截止频率)Rp = 0.0364 (实际通带波动)As = 53 (最小阻带衰减)
5.3.5 窗口法:Matlab 实现
10
20
3040
50
60
-0.1
00.10.20.3
Actual Impulse Response
n
h (n )
00.20.3
1
50
0Magnitude Response in dB
frequency in pi units
D e c i b e l s
从结果看,67 阶Hamming 窗的FIR 数字滤波器的实际阻带衰减为53dB ,通带波动为0.0364dB ,显然满足上面所提的技术要求,其时域和频域响应曲线如下图所示。
113
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
例5.6利用例5.5 给出的设计技术指标,选择Kaiser 窗,
设计所需的低通滤波器。解:
5.3.5 窗口法:Matlab 实现
% Lowpass filter design -Kaiser window wp = 0.2*pi; ws = 0.3*pi; As = 50;tr_width = ws -wp;
N = ceil((As-7.95)/(14.36*tr_width/(2*pi))+1) + 1n=[0:1:N];
beta = 0.1102*(As-8.7)wc = (ws+wp)/2;
h = fir1(N,wc/pi,Kaiser(N+1,beta));[db,mag,pha,grd,w] = freqz_m(h,[1]);delta_w = 2*pi/1000;
As = -round(max(db(ws/delta_w+1:1:501))) % Min Stopband Attenuation
114
北京邮电大学电信工程学院多媒体通信中心门爱东
D igital S ignal P rocessing , Men Aidong, Multimedia Telecommunication Centre, BUPT
运行结果如下:
N = 61 (滤波器阶数)beta = 4.5513
As = 51 (实际的阻带衰减)
5.3.5 窗口法:Matlab 实现
10
203040
50
60
-0.1
00.10.20.3Actual Impulse Response
n
h (n )
00.20.3
1
50
0Magnitude Response in dB
frequency in pi units
D e c i b e l s