离散时间滤波器的设计

8.1技术指标8.2由连续时间滤波器设计IIR滤波器8.2.1脉冲响应不变法8.2.2双线性变换法8.3窗函数法设计FIR滤波器8.3.1设计思想8.3.1布莱克曼窗族8.3.2凯泽窗族8.4FIR滤波器的等波纹最佳逼近设计,离散时间滤波器的设计,滤波器的设计与实现包括以下几个步骤：(1)根据需要制定滤波器的技术指标即频率响应，包括幅度和相位两方面的要求；(2)用因果稳定的离散时间系统逼近这些技术指标。对于FIR滤波器就是找到一个频率响应满足指标要求的因果有限长的单位脉冲响应；而对IIR滤波器则是找到一个频率响应满足指标要求的因果稳定的有理系统函数；(3)为实现该滤波器选择一有限精度的软件算法或硬件结构；(4)验证在选定的实现算法和精度情况下频响特性是否符合指标要求，如果产生了偏差则对步骤（3）的结果进行修正；(5)实现该滤波器。,8.1技术指标,等效的模拟系统,通带容差,阻带容差,通带截止频率,阻带截止频率,3dB/半功率截止频率,幅度响应容限图,若给出的是等效模拟系统的幅度响应,即横坐标是,给出的截止频率是,则需转换得到数字的截止频率:,对数增益:,3dB截止频率对应的对数增益:,在工程实际中，相对指标最受欢迎。,8.2由连续时间滤波器设计IIR滤波器,首先设计一个连续时间滤波器（称为原型连续时间滤波器），得到系统函数Hc(s)；再将系统函数Hc(s)转换成离散时间滤波器的系统函数H(z)。,Hc(s)到H(z)的转换过程实际上是s平面到z平面的映射过程，该映射需要满足两个要求：一是要求所得到的离散时间滤波器的频率响应保持连续时间滤波器频率响应的基本特征，即s平面的虚轴映射到z平面的单位圆，且离散时间滤波器与连续时间滤波器的频率响应的类型保持不变；二是要求将一个因果稳定的连续时间滤波器转换成一个因果稳定的离散时间滤波器，即将在s左半平面的极点全部转换成在z平面单位圆内的极点。,虚轴映射到单位圆,左半平面的极点映射到单位圆内,S平面映射到Z平面,8.2.0连续时间滤波器设计简介,区别：（1）波动不同（2）相同阶数，性能增加（3）设计复杂度增,3种常用的连续时间滤波器的幅度响应：,巴特沃斯滤模拟低通波器的设计公式：,幅度平方函数：,计算阶数和3分贝截止频率：,OR,计算极点:,写出系统函数：,向上取整,全极点型,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)gridon,设计模拟低通滤波器:,举例,输出:N=4Wc=8.1932e+003Bs=1.0e+015*00004.5063As=1.0e+015*0.00000.00000.00000.00144.5063,幅度响应,设计模拟低通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)gridon,举例,设计模拟低通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)gridon,举例,设计模拟高通滤波器:,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)gridon,举例,带通带阻滤波器的通带和阻带截止频率分别是2维,带阻用参数stop,8.2.1脉冲响应不变法1.设计思想：,使离散时间滤波器的单位脉冲响应为原型连续时间滤波器的单位脉冲响应的取样:,就可以根据原型连续时间滤波器的系统函数确定离散时间滤波器的系统函数:,先将连续时间滤波器的系统函数分解成部分分式：,2.系统函数转换方法,再写出离散时间滤波器的系统函数:,证明见课堂笔记,3.因果稳定性,极点的映射关系:,4频率响应的映射特点,如果连续时间滤波器是带限的，即则,混迭小的情况下，近似认为频响是线性关系。,离散时间滤波器与原型连续时间滤波器的频率响应的自变量之间的映射关系是：,虚轴映射到单位圆是一对多，频率轴每绕一圈频响累加一次，所以频响有混迭。,(1)离散时间滤波器的指标原型连续时间滤波器的指标,(2)设计,(3),注意：如果给出的等效模拟指标，则必须首先按照比例T线性转换成数字指标。,5.设计步骤,6.关于Td,若从离散时间指标开始，则Td的取值大小不影响混迭程度和离散时间系统。所以可以任意取值。为了简单起见，常令,优点:频率的线性映射;缺点：频率响应有混迭，应用只限于带限滤波器的设计，应用于高通或带阻滤波器设计时会产生严重的频响混叠。,7.特点,wp=0.2*pi;ws=0.4*piap=1;as=12Td=1;Wp=wp/Td;Ws=ws/TdN,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*)holdonH,w=freqz(Bz,Az);plot(w/pi,20*(log10(abs(H)axis(0.2,0.4,-20,0)grid,举例,输出：Bs=0000.5150As=1.00001.60311.28490.5150Bz=00.14530.08550Az=1.0000-1.47820.9106-0.2013,手算：,N,Wc=buttord(0.6*pi,0.5*pi,1,12,s)Bs,As=butter(N,Wc,high,s)Bz,Az=impinvar(Bs,As,1)H,W=freqs(Bs,As);plot(W/pi,20*(log10(abs(H),r)figureH,w=freqz(Bz,Az);plot(w/pi,20*(log10(abs(H),举例,思考本例的问题出在哪里。,8.2.2双线性变换法,1.设计思想,s平面的整个虚轴非线性地一对一地映射到z平面的单位圆上。,2.系统函数转换方法,如果0,则|z|1，即s左半平面的所有点都映射到z平面的单位圆内，包括极点，所以如果原型连续时间滤波器因果稳定则离散时间滤波器也一定因果稳定。,3.因果稳定性,|z|=1.即虚轴映射到单位圆上(一对一映射).,将带入上式得到,如果则,4.频率响应的映射特点,频率响应间的关系,需要采用预畸变确定原型连续时间滤波器的技术指标:,1)离散时间滤波器指标原型连续时间滤波器指标,2),3),5.设计步骤,注意：如果给出的等效模拟指标，则必须首先用比例T线性转换成数字指标。,6.关于Td,可任取，一般取1,优点：避免了脉冲响应不变法中的频率响应混叠问题；缺点:引入了频率失真，频率失真问题可以通过预畸变加以解决。只适用于滤波器具有近似理想的分段恒定幅度响应的情况，并且它无法实现线性的幅度或相位映射,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)holdonH,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=20N,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.450.55,0.40.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);gridon,Output:N=2wc=0.44100.5590B=0.02710-0.054100.0271A=1.000001.483800.5920,8.2节总结,1.设计步骤,2.冲击响应不变法：频率轴线性多对一映射,频响有混迭,不适用于高通等.双线性变换法：频率轴有畸变一对一映射,频响无混迭,不适用于微分器.,8.3窗函数法设计FIR滤波器,8.3.1设计思想,（1）先对具有广义线性相位的理想频率响应作傅里叶反变换，得到无限长相对M/2偶对称的单位脉冲响应,（2）然后对hdn加长度为M的相对M/2偶对称的窗进行截短，得到有限长的相对M/2偶对称的序列,hn就是一个逼近理想频率响应的因果FIR滤波器的单位脉冲响应，其频率响应是：,加窗的方式（包括窗形状和窗长）决定了滤波器的频率响应对理想频率响应的逼近程度。,窗函数的主瓣宽（与窗长成反比）越大，则滤波器的过渡带宽越宽；窗函数的旁瓣面积与主瓣面积之比越大，即旁瓣峰值与主瓣峰值之比（称最大旁瓣相对幅度：对某种窗形状，该值与窗长无关）越大，则滤波器的通带和阻带误差越大，即阻带衰减越小。因此，使滤波器性能好的窗函数应满足：主瓣宽尽量窄，且旁瓣相对幅度尽量小，即接近一个脉冲函数，以使滤波器的频率响应尽可能接近理想频率响应。,8.3.1布莱克曼窗族,布莱克曼窗族（除了三角形窗）,（a）矩形窗,（b）巴特利特窗,（c）汉宁窗,（d）海明窗,（e）布莱克曼窗,M=20,8.3.2凯泽窗族,M=20,可取任意正实数,修正的零阶贝塞尔函数,M=20,滤波器的阻带衰减和过度带宽与窗函数的形状及长度的关系：,过度带宽,同样的阻带衰减，采用凯泽窗比布窗过度带宽略小。,阻带衰减,1.求理想单位脉冲响应,窗函数法设计步骤,将4个指标转换成：理想截止频率，阻带衰减和过度带宽。,2.根据通带阻带误差确定窗形状:,（1）Blackman窗查表;（2）window窗计算:,3.根据过度带宽估计窗长(对高通和带阻滤波器M取偶数）,4.加窗,(1)Blackman窗查表:,(2)Kaiser床计算:,5.MATLAB验证频响,调整直到满足指标,举例,(2),(4),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);gridon,（5）MATLAB验证,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);gridon,h=0.00010.0023-0.00400.00040.0104-0.01700.00090.0358-0.05380.00140.1288-0.27270.3357-0.27270.12880.0014-0.05380.03580.0009-0.01700.01040.0004-0.00400.00230.0001,（6）调整后,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),如果M取奇数：M=25,8.4FIR滤波器的等波纹最佳逼近设计,阶数M给定情况下性能“最好的”FIR滤波器的设计。最优的准则：最大最小准则（使最大误差最小化）算法：Parks-McClellan算法或称Remez算法。,设计思想,以零相位型广义线性相位FIR滤波器为例。求单位脉冲响应,要求其加权逼近误差函