实验三 数字滤波器设计 ( Filter Design) 01:32

数字滤波器设计 数字滤波器的分类 窗函数法设计FIR数字滤波器 IIR滤波器设计 01:32

数字滤波器分类 经典数字滤波器从功能上可以分为四种 四种理想滤波器的幅频响应 低通(LP)、高通(HP)、带通(BP)、带阻(BS)滤波器。 01:32

理想滤波器的单位脉冲响应 低通 高通 带通 带阻 很容易看出，以上的理想滤波器都是物理不可实现的，其根本的原因是通带和阻带之间频率响应有阶跃突变。为了物理上可实现，在一个频带到另一个频带之间设置过渡带。 01:32 4

若给出的是等效模拟系统的幅度响应,即横坐标是 , 给出的截止频率是 ,则需转换得到数字的截止频率: 通带容差 幅度响应容限图 3dB/半功率截止频率 阻带容差 阻带截止频率 通带截止频率 若给出的是等效模拟系统的幅度响应,即横坐标是 , 给出的截止频率是 ,则需转换得到数字的截止频率:

1、基于窗函数法的FIR数字滤波器设计 01:32

FIR数字滤波器的特点 长度为N的FIR数字滤波器 传递函数： 滤波器系数： 01:32

FIR滤波器的线性相位条件 当FIR系统的单位脉冲响应满足奇对称或者偶对称时，能保证我们设计出来的FIR滤波器是线性相位的 相位延迟： 群延迟： 相位延迟为常数 严格线性相位 群延迟为常数 广义线性相位 当FIR系统的单位脉冲响应满足奇对称或者偶对称时，能保证我们设计出来的FIR滤波器是线性相位的 第I类： h(n)偶对称，N为奇数；第II类： h(n)偶对称，N为偶数； 第III类：h(n)奇对称，N为奇数；第IV类：h(n)奇对称，N为偶数。 01:32

表 1 四种线性相位 FIR 滤波器的性质 类型 I II III IV 长度 N 奇 偶 h [ n ] 的对称性 偶对称 奇对称 |H ( w )| 关于 =0 =p 的周期 2 p 4 b 0.5 |H (0)| 任意 |H (p)| 可适用的 滤波器类型 LP,HP,BP BS 等 LP, BP 微分器 ,Hilbert 变换器 变换器， BP,HP ,BP 01:32

8.3 窗函数法设计FIR滤波器 8.3.1 设计思想 （1）先对具有广义线性相位的理想频率响应作傅里叶反变换，得到无限长相对M/2偶对称的单位脉冲响应 （2）然后对hd[n]加长度为M的相对M/2偶对称的窗进行截短，得到有限长的相对M/2偶对称的序列 h[n]就是一个逼近理想频率响应的因果FIR滤波器的单位脉冲 响应，其频率响应是： 加窗的方式（包括窗形状和窗长）决定了滤波器的频率响应对理想频率响应的逼近程度。

加方窗得到的实际低通滤波器 单位脉冲响应 幅频响应 由幅频响应图看出，时域上截短和移位使得实际低通滤波器的频率响应和理想低通滤波器的有两点不同 在截止频率 附近缓慢下降，即出现了“过渡带”； 通带和阻带都出现了“纹波”。

理想低通滤波器频谱和方窗频谱的卷积过程 * * * * 矩形窗频谱幅度函数： 加窗后低通滤波器频谱幅度函数的特点： ①改变了理想频响的边沿特性，形成过渡带，宽度取决于WR(ω)的主瓣宽度。（方窗： ） （决定于窗长和窗的形状） ②过渡带两旁产生肩峰和余振（带内、带外起伏），取决于WR(ω)的旁瓣，旁瓣多，余振多；旁瓣相对幅度大，肩峰强，与 N无关。 （通带和阻带波纹的高度，取决于窗形状） 正肩峰产生的主要原因是：与第一张图相比，第三张图中右边的旁瓣，在积分中不参加运算，因此积分中负的少了。 * * * *

8.3.1 布莱克曼窗族

布莱克曼窗族（除了三角形窗）

（a）矩形窗 （b）巴特利特窗 （d）海明窗 （c）汉宁窗 M=20 提问：哪种窗最好，没有绝对好的 （e）布莱克曼窗

8.3.2 凯泽窗族 修正的零阶贝塞尔函数 可取任意正实数 M=20

M=20

滤波器的阻带衰减和过度带宽与窗函数的形状及长度的关系： 同样的阻带衰减，采用凯泽窗比布窗过度带宽略小。

窗函数法设计步骤 将4个指标转换成：理想截止频率，阻带衰减和过度带宽。 1. 求理想单位脉冲响应

2.根据通带阻带误差确定窗形状 : （1）Blackman窗查表; （2） Kaiser窗计算: （1）需人做，（2）MATLAB做

3.根据过度带宽估计窗长 (对高通和带阻滤波器M取偶数） (1) Blackman窗查表: (2) Kaiser窗计算: 4. 加窗 5. MATLAB验证频响 调整 直到满足指标

举例 (1) (2) 选海明 (3) 只有兰线需手工 (4)

（5）MATLAB验证 h=fir1(26,0.62, 'high',hamming(27)) H=fft(h,512); %或者freqz(h,1) plot([0:511]/256,20*log10(abs(H))) axis([0.5,0.7,-50,0]); grid on

（6）调整后 h=fir1(24,0.665, 'high',hamming(25)) H=fft(h,512); plot([0:511]/256,20*log10(abs(H))) axis([0.5,0.7,-80,0]); grid on h= 0.0001 0.0023 -0.0040 0.0004 0.0104 -0.0170 0.0009 0.0358 -0.0538 0.0014 0.1288 -0.2727 0.3357 -0.2727 0.1288 0.0014 -0.0538 0.0358 0.0009 -0.0170 0.0104 0.0004 -0.0040 0.0023 0.0001

如果M取奇数：M=25 n=0:25; h=(sinc(n-12.5)-sinc((n-12.5)*0.665)*0.665).*hamming(26)'; H=fft(h,512); plot([0:511]/256, 20*log10(abs(H)))

fdatool

FIR数字滤波器的优点 (与IIR数字滤波器比较) 容易设计成线性相位 避免被处理的信号 产生相位失真，这一特点在宽频带信号处理、阵 列信号处理、数据传输等系统中 ，非常重要。 系统总是稳定的，无稳定性问题 因果性总是满足 任何一个非因果的有限长序列，总可以通过一定的延时，转变为因果序列。 可利用FFT实现滤波运算 无反馈运算，运算误差小 01:32 27

FIR数字滤波器的缺点(与IIR数字滤波器比较) 因为没有极点，要获得好的过渡带特性，需以较高的阶数为代价； 无法利用模拟滤波器的设计结果，一般无解 析设计公式，要借助计算机辅助设计程序完成。 01:32 28

IIR滤波器设计 01:32

8.1 引言 设计IIR滤波器的过程即为求系统函数H(z)的过程 可见，关键为确定H(z)的系数组{ar}、 {bk} 8.1 引言 设计IIR滤波器的过程即为求系统函数H(z)的过程 可见，关键为确定H(z)的系数组{ar}、 {bk} 直接设计法是从期望的IIR数字滤波器的设计指标出发，直接计算IIR数字滤波器的零极点{zr}、{pk}或H(z)的系数组{ar}、{bk}的值。 直接设计法通常用最优化方法计算，即给定目标频率响应后，基于某一误差准则，不断调整滤波器的系数，使目标频率响应和实际频率响应之间的误差最小。当误差不能进一步减小或者达到预先定义的迭代次数时，算法结束并给出最终滤波器系数。 或零、极点{zr} 、 {pk} 。 IIR滤波器设计方法 借助模拟滤波器的设计方法 直接设计法 01:32

借助模拟滤波器的设计方法 模拟滤波器的设计方法非常成熟，许多典型系统有成熟的公式、图表可以查阅，便于设计 为了应用模拟滤波器设计已有成果，需要从 模拟域变换到数字域 借助模拟滤波器的IIR数字滤波器设计法就是： 先设计一个合适的模拟滤波器，然后将它变换成数字滤波器 01:32

(1)通过频率变换，将期望得到的数字滤波器技术指标转换为模拟低通滤波器的技术指标； 借助模拟滤波器的设计数字滤波器的步骤 (1)通过频率变换，将期望得到的数字滤波器技术指标转换为模拟低通滤波器的技术指标； (2)根据转换后的技术指标设计模拟低通滤波器的系统函数 ； (3)利用从模拟域到数字域的变换式，将模拟低通滤波器的系统函数 转换成所要设计的数字滤波器的系统函数 01:32

8.2.0 连续时间滤波器设计简介 3种常用的连续时间滤波器的幅度响应： 区别： （1）波动不同 （2）相同阶数， 性能增加 （3）设计 复杂度增

巴特沃斯滤模拟低通波器的设计公式： 幅度平方函数： 计算阶数和3分贝截止频率： 向上取整 OR 计算极点: 写出系统函数： 全极点型

举例 设计模拟低通滤波器: [N,Wc]=buttord(2000*pi,4000*pi,1,15, 's')   设计模拟低通滤波器: 举例 [N,Wc]=buttord(2000*pi,4000*pi,1,15, 's') [Bs,As]=butter(N,Wc, 's') [H,W]=freqs(Bs,As); plot(W/2/pi,20*(log10(abs(H)))) axis([1000,2000,-16,0]) grid on 输出: N = 4 Wc = 8.1932e+003 Bs = 1.0e+015 * 0 0 0 0 4.5063 As = 1.0e+015 * 0.0000 0.0000 0.0000 0.0014 4.5063

幅度响应

举例 设计模拟低通chebyI型滤波器: [N,Wc]=cheb1ord(2000*pi,4000*pi,1,15, 's')   设计模拟低通chebyI型滤波器: [N,Wc]=cheb1ord(2000*pi,4000*pi,1,15, 's') [Bs,As]=cheby1(N,1,Wc, 's') [H,W]=freqs(Bs,As); plot(W/2/pi,20*(log10(abs(H)))) axis([0,4000,-30,0]) grid on

举例 设计模拟低通chebyII 型滤波器: [N,Wc]=cheb2ord(2000*pi,4000*pi,1,15, 's')   设计模拟低通chebyII 型滤波器: 举例 [N,Wc]=cheb2ord(2000*pi,4000*pi,1,15, 's') [Bs,As]=cheby2(N,15,Wc, 's') [H,W]=freqs(Bs,As); plot(W/2/pi,20*(log10(abs(H)))) axis([0,4000,-30,0]) grid on

举例 设计模拟高通滤波器: [N,Wc]=buttord(4000*pi,2000*pi,1,15, 's')   设计模拟高通滤波器: [N,Wc]=buttord(4000*pi,2000*pi,1,15, 's') [Bs,As]=butter(N,Wc, 'high', 's') [H,W]=freqs(Bs,As); plot(W/2/pi,20*(log10(abs(H)))) axis([0,4000,-16,0]) grid on 带通带阻滤波器的通带和阻带 截止频率分别是2维,带阻用参 数’stop’

从模拟滤波器到数字滤波器 01:32

S平面映射到Z平面 虚轴映射到单位圆 左半平面的极点映射到单位圆内

8.3 IIR数字滤波器设计的脉冲响应不变法 脉冲响应不变法变换原理 数字滤波器的单位脉冲响应 T—采样周期 模拟滤波器的单位脉冲响应 01:32

脉冲响应不变法的基本转换关系 时域采样等效于从S平面到Z的直接映射： 考虑 只有一个极点 ： 由 映射成Z平面上的极点 对应的数字滤波器： 考虑 只有一个极点 ： 由 映射成Z平面上的极点 对应的数字滤波器： 01:32

多个极点的情况 ， 01:32

极点： 系数相同： 稳定性不变：s 平面 z 平面 s 平面的单极点 z 平面上 处的单极点 虽然脉冲响应不变法能保证s 平面极点与z 平面极点有这种代数对应关系；但并不等于整个s 平面与z 平面有这种代数对应关系；特别是DF的零点位置与AF的零点位置没有这种代数对应关系。 01:32

当T 很小时，数字滤波器增益很大，易溢出，需修正 令： 则： 01:32

抽样率为 ，试用脉冲响应不变法设计一性能近似上述模拟滤波器的IIR数字滤波器。 例8.3.2：模拟滤波器的系统函数为 抽样率为 ，试用脉冲响应不变法设计一性能近似上述模拟滤波器的IIR数字滤波器。 解： 对应地 乘以因子 ，整理得 01:32

举例 wp=0.2*pi; ws=0.4*pi ap=1; as=12 Td=1; Wp=wp/Td; Ws=ws/Td [N,Wc]=buttord(Wp,Ws, ap , as, 's') [Bs,As]=butter(N,Wc, 's') [Bz,Az]=impinvar(Bs,As,1/Td) [H,W]=freqs(Bs,As); plot(W/pi,20*(log10(abs(H))), 'r*') hold on [H,w]=freqz(Bz,Az); plot(w/pi,20*(log10(abs(H)))) axis([0.2,0.4,-20,0]) grid

输出： Bs = 0 0 0 0.5150 As =1.0000 1.6031 1.2849 0.5150 Bz = 0 0.1453 0.0855 0 Az =1.0000 -1.4782 0.9106 -0.2013

手算：

脉冲响应不变法的缺点 从S平面到Z平面的直接映射： Z平面 S平面 映射关系不是一一对应的 01:32

时域 频域 模拟 离散化 数字 周期化 采样率越低，混叠越严重 混迭失真 01:32

h(n)完全模仿模拟滤波器的单位抽样响应 时域逼近良好 优缺点 优点： h(n)完全模仿模拟滤波器的单位抽样响应 时域逼近良好 频率映射保持线性关系： 线性相位模拟滤波器转变为线性相位数字滤波器 缺点： 频率响应混迭 只适用于限带的低通、带通滤波器 01:32

8.4 IIR数字滤波器设计的双线性变换法 脉冲响应不变法的缺点：频域混叠 造成频域混叠的原因：从s平面到z平面的直接映射 解决方法：压缩频率轴 Z平面 S平面 因为小omiga =大omiga乘以Ts ； 所以，这里频率轴的“ 大omiga” 的 Pi/T，对应的是小omiga的pi 映射关系不是一一对应的 01:32

双线性变换的角度关系 为避免混叠，将s平面上的 j 轴 的(  ,)先压缩到 再映射到z平面： W1 : W : ， Ts 为采样间隔，c为常数 W1 : W : 01:32

上式在s平面上解析延拓,令

S平面 S1平面 Z平面 一一对应 01:32

双线性变换 双线性变换： 双线性变换的频率变换关系： 上式表示的频率变换关系是非线性的，通常称为频率预畸 c: 变换常数 01:32

1.变换常数c 变换常数c的选择并不影响最终的设计结果H(z)，但是用不同的常数c会得到不同的模拟滤波器 ，具体原因在后面频率变换中解释。因此如果我们并不关心 ，常数c 可以取1。 ， 01:32

2.用双线性变换设计IIR数字滤波器的步骤 1） 将数字滤波器的频率指标{wk}转换为模拟滤波器的频率指标{Wk}； 2）设计模拟滤波器的Ha(s)； 3）用双线性变换将Ha(s)转换为H(z)。 01:32

7. 特点 优点：避免了脉冲响应不变法中的频率响应混叠问题； 缺点:引入了频率失真，频率失真问题可以通过预畸变加以解决。只适用于滤波器具有近似理想的分段恒定幅度响应的情况，并且它无法实现线性的幅度或相位映射 例如模拟微分器（幅度响应线性）不能通过双线性变换法 得到数字微分器。

举例 wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20; Td=1; Ws=2/Td*tan(ws/2); Wp=2/Td*tan(wp/2) [N,Wc]=buttord(Wp,Ws,ap,as,’s’) [Bs,As]=butter(N,Wc, ‘s’) [Bz,Az]=bilinear(Bs,As,1/Td) [H,W]=freqs(Bs,As); plot(W/pi,20*(log10(abs(H))),’Rx‘) hold on [H,w]=freqz(Bz,Az); plot(w/pi,20*(log10(abs(H)))) ylabel(‘虚线：模拟滤波器幅度[dB] 实线：数字滤波器幅度[dB]’) xlabel(‘虚线：模拟角频率[*π弧度/秒] 实线:数字角频率[*π弧度]’) axis([0.25,0.5,-20,-0.45]) grid

或者 wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20 [N,wc]=buttord(wp/pi,ws/pi,ap,as) [Bz,Az]=butter(N,wc) [H,w]=freqz(Bz,Az); plot(w/pi,20*(log10(abs(H))))

举例 [N,Wc]=buttord(0.6,0.5,1,12); [b,a]=butter(N,Wc, 'high'); [H,w]=freqz(b,a); plot(w/pi,20*(log10(abs(H)))) grid

举例 设计离散时间带通滤波器

[N,wc]=buttord([0.45 0.55],[0.4 0.6],3,10) %双线性变换法 [B,A]=butter(N,wc) [H,w]=freqz(B,A); plot(w/pi,20*(log10(abs(H)))) ylabel(‘20log|H(ejω)| [dB]’) xlabel(‘数字角频率[*π弧度]’) axis([0.4,0.6,-10,0]); grid on Output: N = 2 wc = 0.4410 0.5590 B = 0.0271 0 -0.0541 0 0.0271 A = 1.0000 0 1.4838 0 0.5920

8.2节 总结 1. 设计步骤 2.脉冲响应不变法： 频率轴线性多对一映射,频响有混迭,不适用于高通等. 双线性变换法： 频率轴有畸变一对一映射,频响无混迭,不适用于微分器.