第6章 有限脉冲响应数字滤波器的设计 IIR的设计 Specifications Desired IIR 脉冲响应不变法 阶跃响应不变法

Presentation on theme: "第6章 有限脉冲响应数字滤波器的设计 IIR的设计 Specifications Desired IIR 脉冲响应不变法 阶跃响应不变法"— Presentation transcript:

1 第6章 有限脉冲响应数字滤波器的设计 IIR的设计 Specifications Desired IIR 脉冲响应不变法 阶跃响应不变法
双线性变换法 Butterworth, Chebyshev, Ellipse LP to HP, BP, BS

2 回顾IIR 引入FIR 优点： 利用模拟滤波器成熟的理论和设计图表，保留了模拟滤波器优良的幅度特性 缺点：
只考虑了幅度 特性，没考虑相位特性，一般具有非线性相位特性 如果要得到线性相位特性，需要增加相位校正网络，使滤波器变得复杂，增加了成本，而且难以保证严格的线性相位特性 引入FIR 优点： 在满足幅度特性技术要求的同时，很容易实现严格的线性相位特性

3 FIR滤波器概述 系统函数 FIR滤波器设计方法 H(z)在z平面有N-1个零点，在原点处有一个N-1重极点，因此H(z)永远稳定。
不同于IIR。IIR借助模拟滤波器，FIR直接选择有限长度的h(n)，使得频响函数满足要求。 三种主要设计方法：窗函数法、频率采样法和Chebyshev等波纹法

4 6.1 线性相位FIR数字滤波器回顾 6.2 利用窗函数法设计FIR滤波器 6.3 利用频率采样法设计FIR滤波器 6.4 利用切比雪夫逼近法设计FIR滤波器 6.5 IIR和FIR数字滤波器的比较 6.6 FDATool

5 6.1 线性相位FIR数字滤波器回顾 什么是线性相位FIR？ 考虑长度为N的h(n)，系统函数为： 频率响应函数为： (6.1.1)
(6.1.2) Hg(ω)称为幅度特性，θ(ω)称为相位特性。注意，Hg(ω)不同于|H(ejω)|，Hg(ω)为ω的实函数，可能取负值，而|H(ejω)|总是正值。

6 严格地说， (7.1.4) 中θ(ω)不具有线性相位，但以上两种情况都满足群时延是一个常数，即
H(ejω)线性相位是指： θ(ω)是ω的线性函数，即 τ为常数 (6.1.3)   或θ(ω)满足下式： ,θ0是起始相位 (6.1.4) 严格地说， (7.1.4) 中θ(ω)不具有线性相位，但以上两种情况都满足群时延是一个常数，即 第一类线性相位 第二类线性相位

7 线性相位条件： 注意：充分条件 第一类线性相位： h(n)是实序列且对(N-1)/2偶对称，即 第二类线性相位：
(6.1.5) 第二类线性相位： h(n)是实序列且对(N-1)/2奇对称，即 (6.1.6)

8 7.2 利用窗函数法设计FIR滤波器 理想低通、高通、带通、带阻数字滤波器幅度特性

9 基本思路 理想滤波器为Hd(z) hd(n)是与其对应的单位脉冲响应。频率响应为 求和上下限为什么要取无穷？？？
对ω是连续的，需要无穷个这样的函数，才能逼近分段连续的理想滤波器的频率响应

10 基本思路 Unfortunately, hd(n)序列不可能无限长，我们只能选取其中的N个，组成实际的序列h(n)，比如：
wN(n)是长度为N的窗口。得到： 可以选择多种形状的窗口 注意： 通常选择构造长度为N的第一类线性相位FIR滤波器，即要求h(n)关于n=(N-1)/2偶对称。

11 如果选择矩形窗口RN(n) 图 理想低通的单位脉冲响应及矩形窗

12 窗函数法的一些性质 在时域上，截取了理想滤波器的部分单位单位响应序列，作为实际滤波器的单位脉冲序列
用有限长的h(n)代替无限长的hd(n) ，必然引起误差，即Gibbs效应，使得在频域上 过渡带加宽 通带、阻带波动 N越大，误差越小，但系统越复杂、成本越高 设计的基本目标 构造窗函数，减少Gibbs效应 满足线性相位特性

13 窗函数法导致的频域误差 h(n)对应的频率响应是什么样子？ 可以这样计算： 也可以这样计算：根据傅立叶变换的频域卷积定理，有
该用哪种方法呢？？？

14 矩形窗 wN(n) =RN(n) RN(ω)称为矩形窗的幅度函数

15 RN(ω)的特性 主瓣 旁瓣 RN(ω)在[-2π/N, 2π/N]区间的一段波形称为主瓣，其余较小的波动称为旁瓣。

16 设Hd(z)为希望逼近的线性相位低通滤波器，其频率响应为：
得到实际滤波器H(z)的频率响应为： 线性相位 Hd(ω)和RN(ω)的卷积 最终得到实际滤波器H(z)的幅度特性H(ω)为：

17 Hd(ω)和RN(ω)的卷积形成的波形：
ω=0，H(0)为(a)、(b)两波形的积分，即RN(ω)在[-ωc,ωc]之间积分。当ωc >>2π/N，近似于RN(ω)在[-π,π]之间积分，归一化该积分为1。 ω=ωc，当ωc >>2π/N，近似为RN(ω)一半波形的积分，近似值为1/2 ω=ωc-2π/N，主瓣在[-ωc,ωc]之内，最大负旁瓣在[-ωc,ωc]之外，Hd(ω)有一个最大正峰 ω=ωc+2π/N，主瓣在[-ωc,ωc]之外，最大负旁瓣在[-ωc,ωc]之内，Hd(ω)有一个最大负峰 最大正负峰距离4π/N

18 通过以上分析可知，对hd(n)加矩形窗处理后，H(ω)和原理想低通Hd(ω)差别有以下三点：
(1)在理想特性不连续点ω=ωc附近形成过渡带。过渡带的宽度，近似等于RN(ω)主瓣宽度，即4π/N。 (2)通带内增加了波动，最大的峰值在ωc-2π/N处。阻带内产生了余振，最大的负峰在ωc+2π/N处。 (3)通带和阻带中波纹的情况与窗函数的幅度函数有关， RN(ω)旁瓣幅度的大小直接影响了H(ω)波纹幅度的大小。 选择合适的窗函数会改变旁瓣幅度，因此可以改变H(ω)波纹幅度

19 增加矩形窗长度（即加大N）的效果 在主瓣附近，RN(ω)可近似为 N增大时： 主瓣幅度增加，旁瓣也增加，相对值不变 增加N并不能减小波纹幅度
主瓣、旁瓣距离变小，波动频率加快 增加N并不能减小波纹幅度

20 几种典型窗函数 四种波形图： 窗函数时域波形、幅度特性曲线 理想低通滤波器加窗后的单位脉冲响应、幅度特性曲线 定义几个参数：
旁瓣峰值αn——最大旁瓣的最大值相对于主瓣最大值的衰减值(dB) 过渡带宽度Bg——该窗函数设计的FIR数字滤波器的过渡带宽度 阻带最小衰减αs——该窗函数设计的FIR数字滤波器的阻带最小衰减 20

21 1. 矩形窗(Rectangle Window)
其频率响应为 幅度特性为 指标： αn=-13dB; Bg=4π/N; αs=-21dB 21

22 2. 三角形窗(Bartlett Window)
(6.2.8) 其频率响应为 (6.2.9) 幅度特性为 22

23 αn=-25dB; Bg=8π/N; αs=-25dB
指标： αn=-25dB; Bg=8π/N; αs=-25dB 23

24 当N>>1时，N-1≈N，幅度特性为
3. 汉宁(Hanning)窗——升余弦窗 其频率响应为 当N>>1时，N-1≈N，幅度特性为 旁瓣对消，能量集中在主瓣 24

25 图 汉宁窗的幅度特性 25

26 αn=-31dB; Bg=8π/N; αs=-44dB
指标： αn=-31dB; Bg=8π/N; αs=-44dB 26

27 能量更集中在主瓣，约占99.96%。是Matlab默认窗函数
4. 哈明(Hamming)窗——改进的升余弦窗 其频响函数WHm(ejω)为 其幅度函数WHmg(ω)为 当N>>1时，可近似表示为 能量更集中在主瓣，约占99.96%。是Matlab默认窗函数 27

28 αn=-41dB; Bg=8π/N; αs=-53dB
指标： αn=-41dB; Bg=8π/N; αs=-53dB 28

29 5. 布莱克曼(Blackman)窗 其频率响应函数为 其幅度函数为 旁瓣进一步抵消，但过渡带增宽 29

30 αn=-57dB; Bg=12π/N; αs=-74dB
指标： αn=-57dB; Bg=12π/N; αs=-74dB 30

31 常用的窗函数的时域波形 31

32 (a)矩形窗；(b)巴特利特窗(三角形窗)；(c)汉宁窗； (d)哈明窗；(e)布莱克曼窗
图 常用窗函数的幅度特性 (a)矩形窗；(b)巴特利特窗(三角形窗)；(c)汉宁窗； (d)哈明窗；(e)布莱克曼窗 32

33 图6.2.6 理想低通加窗后的幅度特性(N=51,ωc=0.5π) (a)矩形窗；(b)巴特利特窗(三角形窗)；(c)汉宁窗；
(d)哈明窗；(e)布莱克曼窗 33

34 6. 凯塞—贝塞尔窗(Kaiser-Basel Window)
前五种为参数固定窗函数，旁瓣幅度是固定的。凯塞—贝塞尔窗参数可调，是一种最优窗函数。 式中 I0(x)是零阶第一类修正贝塞尔函数，可用下面级数计算： 34

35 一般I0(x)取15~25项，便可以满足精度要求。α参数可以控制窗的形状。一般α加大，主瓣加宽，旁瓣幅度减小，典型数据为4<α<9。当α=5.44时，窗函数接近哈明窗。α=7.865时，窗函数接近布莱克曼窗。凯塞窗的幅度函数为 35

36 凯塞窗参数对滤波器的性能影响 36

37 六种窗函数的基本参数 37

38 用窗函数设计FIR滤波器的步骤 (1)根据技术要求确定待求滤波器的单位取样响应hd(n)。如果给出待求滤波器的频响为Hd(ejω)，那么单位取样响应用下式求出： 根据频率采样定理，hM(n)与hd(n)应满足如下关系： 38

39 例如，理想低通滤波器的单位取样响应hd(n)为：
(3) 计算滤波器的单位取样响应h(n)， h(n)=hd(n)w(n) (4)验算设计出的滤波器频率响应是否满足指标： 39

40 例6.2.1 用矩形窗、汉宁窗和布莱克曼窗设计FIR低通滤波器，设N=11,ωc=0.2πrad。
解 用理想低通作为逼近滤波器，有 40

41 用汉宁窗设计： 用布莱克曼窗设计： 41

42 图 例6.2.1的低通幅度特性 42

43 Using Matlab 例：用Hanning窗设计高通FIRDF，要求ωp=π/2，ωs=π/4，ap=1dB，as=40dB。
%HP_Window wp=.5; ws=.25; Bt=wp-ws; N0=ceil(6.2/Bt); N=N0+mod(N0+1,2); wc=(wp+ws)/2; hn=FIRl(N-1,wc,'high',hanning(N)); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)'); 查表可知，Hanning窗满足40dB要求。 43

44 N=25, wc=0.375pi; hn = 44

45 例：对模拟信号进行低通滤波处理。要求Ωp=1.5kHz, Ωs=2.5kHz, ap=1dB, as=40dB。采样频率Fs=10kHz。
用凯塞窗 %LP_Window wp=2*1.5/10; ws=2*2.5/10; as=40; Bt=ws-wp; alpha=0.5842*(as-21) *(as-21); M=ceil((as-8)/2.285/Bt); wc=(wp+ws)/2; hn=FIRl(M,wc,kaiser(M+1,alpha)); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)'); 45

46 N=24; hn = 46

47 例：设计带阻FIRDF，要求ωpl=.2π, ωsl=.35π, ωsu=.65π, ωpu=.8π， ap=1dB，as=60dB。
%BS_Window wp=[.2, .8]; ws=[.35, .65]; Bt=ws(1)-wp(1); M=ceil(12/Bt)-1; wp=[(ws(1)+wp(1))/2, (ws(2)+wp(2))/2]; hn=FIRl(M,wp,’stop’,blackman (M+1)); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)'); 47

48 N=80; 48

49 6.3 利用频率采样法设计FIR滤波器 设待设计的滤波器的传输函数用Hd(ejω)表示，对它在ω=0到2π之间等间隔采样N点，得到Hd(k)， 再对N点Hd(k)进行IDFT，得到h(n)， 49

50 式中，h(n)作为所设计的滤波器的单位脉冲响应，其系统函数H(z)为
50

51 FIR滤波器具有线性相位的条件是h(n)是实序列，且满足h(n)=h(N-n-1)，在此基础上我们已推导出其传输函数应满足的条件是：
1.用频率采样法设计线性相位滤波器的条件 FIR滤波器具有线性相位的条件是h(n)是实序列，且满足h(n)=h(N-n-1)，在此基础上我们已推导出其传输函数应满足的条件是： (6.3.4) (6.3.5) (6.3.6) (6.3.7) 51

52 将ω=ωk代入(7.3.4)~(7.3.7)式中，并写成k的函数：
在ω=0~2π之间等间隔采样N点， 将ω=ωk代入(7.3.4)~(7.3.7)式中，并写成k的函数： (6.3.8) (6.3.9) (6.3.10) (6.3.11) 52

53 设用理想低通作为希望设计的滤波器，截止频率为ωc，采样点数N，Hg(k)和θ(k)用下面公式计算：
注意： Hg(N) = Hg(0) (6.3.12) kc是通带内最后一个采样点序号 ωc 53

54 N=偶数时， (6.3.13) 54

55 如果待设计的滤波器为Hd(ejω),对应的单位脉冲响应为hd(n)，
2. 逼近误差及其改进措施 如果待设计的滤波器为Hd(ejω),对应的单位脉冲响应为hd(n)， 则由频率域采样定理知道，在频域0~2π之间等间隔采样N点，利用IDFT得到的h(n)应是hd(n)以N为周期，周期性延拓乘的主值区序列，即 由于理想滤波器对应的hd(n)序列无限长，因此有限的周期N必然带来时域混叠。N越大，混叠越小。 55

56 由采样定理表明，频率域等间隔采样H(k)，经过IDFT得到h(n)，其Z变换H(z)和H(k)的关系为
频率响应可以写成内插的形式： 56

57 在频响间断点附近区间内插几个过渡采样点，是不连续点变成过渡带，可增大阻带衰减，但会增加过渡带宽度
在采样点上，逼近误差为0 在采样点之间，有误差 在Hdg(ω)平滑的区域，误差较小 在Hdg(ω)突变的区域，误差较大 在频响间断点附近区间内插几个过渡采样点，是不连续点变成过渡带，可增大阻带衰减，但会增加过渡带宽度 阻带衰减和过渡带宽度是一对矛盾 57

58 过渡带采样点数m 1 2 3 as 44~54dB 65~75dB 85~95dB 图 理想低通滤波器增加过渡点 58

59 例6.3.1 利用频率采样法设计线性相位低通滤波器，要求截止频率ωc=π/3，as=40dB，过渡带宽度Bt≤π/16。
解： as=40dB时需要增加一个过渡带采样点，即m＝1。 计算N： 59

60 Using Matlab %LP_FreqSample T=input('T='); wp=1/3;Bt=1/16; m=1;
N=ceil((m+1)*2/Bt)+1; N=N+mod(N+1,2); F=[0,wp,wp+2/N,wp+4/N,1]; A=[1,1,T,0,0]; hn=fir2(N-1,F,A,boxcar(N)); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)'); 60

61 T=.38 as= dB 61

62 过渡带插值T不同，as不同 T=.5 as= dB T=.6 as= dB 62

63 窗函数和频率采样法的特点 优点： 缺点： 简单方便，易于实现 边界频率不易精确控制
窗函数法使得通带、阻带波纹幅度相等，频率采样法只能控制阻带波纹幅度，它们都不能分别控制通带和阻带的幅度 所设计的滤波器在阻带边界附近衰减最小，距离阻带边界越远，衰减越大。如果阻带边界的衰减刚好满足要求，那么阻带其它频段有很大富裕量，性价比较低 63

64 6.4 切比雪夫最佳逼近法 什么是最佳逼近？ 设希望设计的滤波器幅度特性为Hd(ω)，实际设计的滤波器幅度特性为Hg(ω)，其加权误差E(ω)用下式表示： E(ω)=W(ω)[Hd(ω)-Hg(ω)] 设计具有线性相位的FIR滤波器，使E(ω)最优 其单位脉冲响应序列h(n) 必须满足一定条件，例如：h(n)=h(n-N-1)且N为奇数 64

65 选择M+1个系数a(n),使加权误差E(ω)的最大值为最小
M=(N-1)/2 最佳逼近问题： 选择M+1个系数a(n),使加权误差E(ω)的最大值为最小 W(ω)可以在不同频段选择不同权重 65

66 该定理指出最佳一致逼近的充要条件是E(ω)在A上至少呈现M+2个“交错”(alternations)，使得
Alternation Theorem 该定理指出最佳一致逼近的充要条件是E(ω)在A上至少呈现M+2个“交错”(alternations)，使得 属于多项式逼近问题，证明已超出要求 66

67 对于低通滤波器 性能指标 67

68 Using Matlab 设计带阻滤波器，通带[0, 0.2pi]、[0.8pi, pi]，ap=1dB；阻带[0.35pi, 0.65pi]，as=60dB。 %BS_Equirip f=[0.2,0.35,0.65,0.8]; m=[1,0,1]; ap=1;as=60; delta1=(10^(ap/20)-1)/(10^(ap/20)+1); delta2=10^(-as/20); rip=[delta1,delta2,delta1]; [M,fo,mo,w]=remezord(f,m,rip); hn=remez(M,fo,mo,w); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)');

69 M=28, N=29, fo=[0, .2, .35, .65, .8, 1], mo=[1, 1, 0, 0, 1, 1], w=[1, ,1] 69

70 设计数字低通滤波器, Ωp=1.5kHz, ap=1dB; Ωs=2.5kHz, as=40dB; 采样频率10kHz
%LP_Equirip Fs=10000; f=[1500,2500]; m=[1,0]; ap=1;as=40; delta1=(10^(ap/20)-1)/(10^(ap/20)+1); delta2=10^(-as/20); rip=[delta1,delta2]; [M,fo,mo,w]=remezord(f,m,rip,Fs); M=M+1; hn=remez(M,fo,mo,w); figure(1);stem(0:(length(hn)-1),hn); fk=0:1/200:1; hk=freqz(hn,1,fk*pi); figure(2);plot(fk,20*log10(abs(hk))); xlabel('Frequency(pi)');ylabel('Magnitude(dB)');

71 M=15, N=16, w=[1, ] 71

72 6.5 IIR和FIR数字滤波器的比较 性能： IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的 对于同样的性能指标，FIR滤波器所要求的阶数一般比IIR滤波器高5-10倍，成本高，信号时延大，但是为线性相位 结构： IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定；另外，由于运算过程中的舍入处理，有限字长效应可能会引起寄生振荡 FIR滤波器采用非递归结构，在理论和实际中都不存在稳定性问题，运算误差引起的噪声较小；此外，FIR可以采用FFT算法实现，运算速度可以大大提高 72

73 IIR滤波器可以借助于模拟滤波器的成果，因此一般都有有效的封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。
设计工具： IIR滤波器可以借助于模拟滤波器的成果，因此一般都有有效的封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。 FIR滤波器计算通带、阻带衰减等无显式表达式，其边界频率不易精确控制；FIR的设计依赖计算工具 灵活性： IIR滤波器虽然设计简单，但主要用于低通、高通、带通和带阻等滤波器，脱离不了典型滤波器的频响特性 FIR滤波器则要灵活很多，易于适应某些特殊的应用 73

74 6.6 FDATool Matlab：FDATool 74

75

76