《信号分析与处理》-第 9 章 滤波器设计

9.1滤波器概述9.1.1滤波器的基本概念滤波器是指具有一定传输特性的信号处理装置。按照输入滤波器的信号是模拟的还是离散的,滤波器可以分为模拟滤波器和数字滤波器两大类。模拟滤波器用于处理输入的模拟信号,其作用是让信号中有用的频率成分通过系统,消除或抑制不必要的频率成分。模拟滤波器从功能上一般可分为四种,即低通(lowpass)、高通(highpass)、带通(bandpass)、带阻(bandstop)滤波器。图9-1所示为四种模拟滤波器的幅频特性图,图中所给的滤波器的幅频响应都是理想情况,在实际上是不可能实现的。在实际工作中,我们设计的滤波器都是在某些准则下对理想滤波器的近似,以保证滤波器在物理上是稳定的、可实现的。下一页返回9.1滤波器概述数字滤波器设计是数字信号处理的重要任务之一。从离散时间系统的滤波功能上同样可把系统分为低通、高通、带通、带阻滤波器,与模拟滤波器不同的是,数字滤波器以2π为周期,在主周期[0,2π]内其幅度特性关于ω=π呈偶对称,ω从0至π变化时,频率越来越高,ω从π至2π变化时,频率越来越低,其幅频特性如图9-2所示。根据系统的单位样值响应h[n]的长度分类,可把离散系统分为有限长单位样值响应滤波器(FiniteImpulseResponse,FIR)和无限长单位样值响应滤波器(InfiniteImpulseResponse,IIR)。上一页下一页返回9.1滤波器概述9.1.2滤波器的技术指标从图9-1和图9-2可见,四种类型的模拟滤波器与数字滤波器(数字角频率ω∈[0,π])的幅频特性图是对应类似的。下面以数字低通滤波器为例讲述其技术指标,模拟低通滤波器技术指标与数字低通滤波器技术指标是类似的,不再赘述。典型的数字滤波器的技术指标以容限图的形式给出,如图9-3所示。H(ejω)为数字低通滤波器的幅度响应,ωp

为通带截止频率,ωs为阻带截止频率,αp、αs分别为通带和阻带的容限。在频带[0,ωp]内的信号频率成分可以通过滤波器,该频带称为通带(Pass-band)。在频带[ωp,ωs]内的信号频率成分经过滤波器时受到抑制,该段频带称为过渡带。在频带[ωs,π]内的信号频率成分经过滤波器后有较大衰减,该段频带称为阻带(Stopband)。上一页下一页返回9.1滤波器概述在工程中设计滤波器时,通带和阻带内的允许衰减一般用分贝(dB)数表示,通带内允许的最大衰减用δp

表示,阻带内允许的最小衰减用δs

。δp

和δs分别定义为9.1.3滤波器通用设计步骤无论是模拟滤波器还是数字滤波器,其设计过程大致由如下几个步骤完成:(1)根据设计要求,确定滤波器的技术指标;上一页下一页返回9.1滤波器概述(2)用一个因果稳定的连续(或离散)时间线性时不变系统函数去逼近这组技术指标;(3)用有限精度算法来实现这个系统函数(包括滤波器结构选择、有限字长效应等);(4)通过电路系统设计、计算机软件或专用数字滤波器硬件来实际技术实现;(5)测试并验证滤波器性能是否达到规定的技术指标。上一页返回9.2模拟滤波器设计9.2.1Butterworth模拟滤波器的设计1.由幅度平方函数确定系统函数设模拟滤波器幅度响应用幅度平方函数Ha(jΩ)2来表示,即由于滤波器的单位样值响应ha(t)是实函数,因而H(jΩ)满足所以下一页返回9.2模拟滤波器设计2.Butterworth低通滤波器设计模拟滤波器的设计都是依据模拟原型滤波器进行的。所谓模拟原型滤波器指的是截止频率为1rad/s的滤波器。后面所讲的各类模拟滤波器和数字滤波器可通过这些低通原型滤波器变换得到。N阶Butterworth模拟低通滤波器的平方幅频响应函数定义为Butter-worth滤波器具有如下特征:(1)当Ω=0时,H(jΩ)=1,即在Ω=0处无衰减。上一页下一页返回9.2模拟滤波器设计(2)当Ω=Ωc时,H(jΩ)=2/2,幅度衰减3dB,此时Ωc称为3dB衰减点。无论阶数N为多少,滤波器的幅频特性都经过该3dB频率点,这就是3dB不变性,如图9-4所示。Butterworth模拟滤波器的设计步骤为:(1)根据设计要求,确定滤波器的技术指标[Ωp,δp]和[Ωs,δs];(2)根据技术指标[Ωp,δ

p]和[Ωs,δs],求出滤波器的阶数N和3dB截止频率Ωc;(3)根据阶数N查表9-1,求得原型函数Han(s);(4)将s=s/Ωc

代入Han(s),得到Butterworth低通滤波器的系统函数H(s)。上一页返回9.3IIR数字滤波器设计9.3.1冲激响应不变法设计IIR滤波器1.变换原理假设模拟滤波器的单位样值响应信号为ha(t),其时域波形如图9-5(a)所示,以Ts为采样周期对ha(t)进行理想采样。Ha(jΩ)是ha(t)的傅里叶变换,其幅度谱如图9-5(b)所示。则采样信号ha(nTs)=ha(t)t=nTs

的傅里叶变换为假设所需设计的数字滤波器的单位样值响应信号为h[n],若令h[n]=Tsha(nTs),对等式两边同时取傅里叶变换可得上一页下一页返回9.3IIR数字滤波器设计由于两种滤波器的单位样值响应ha(t)和h[n]除了发生采样和倍乘Ts外,波形形状基本未发生其他变化,因此该数字滤波器设计方法称为冲激响应不变法,其实现流程如下:冲激响应不变法实现了Ha(s)向H(z)的转换,即S平面向Z平面的映射,其映射关系如图9-6所示。2.数字化过程设模拟滤波器的系统函数Ha(s)只有N个单阶极点,且分母的阶次大于分子的阶次。Ha(s)的部分分式形式为上一页下一页返回9.3IIR数字滤波器设计该滤波器的单位样值响应为按冲激响应不变法的实现流程,可得数字滤波器的单位样值响应为上一页下一页返回9.3IIR数字滤波器设计3.混叠失真由图9--5可见,冲激响应不变法设计的数字滤波器频率响应H(ejω)是模拟滤波器频率响应Ha(jΩ)以2π为周期的周期延拓。时域采样定理表明,只有模拟滤波器的频率响应是限带的,即满足Ha(jΩ)=0,Ω>π/Ts才能使数字滤波器的频率响应不失真的重现模拟滤波器的频率响应,而不会发现频谱混叠失真。但是,任何一个实际的模拟滤波器频率响应都不是严格限带的,将不可避免的出现频谱混叠失真,图9--7说明了这一现象。图9-6表明冲激响应不变法实现S平面向Z平面的多值映射,S平面上每个2π/Ts的横条向Z平面整个平面映射一次,不严格的限带滤波器在映射过程中造成了频谱混叠,与图9-5描述的本质是一样的。上一页下一页返回9.3IIR数字滤波器设计4.冲激响应不变法的优缺点从以上讨论可以看出,使用冲激响应不变法来设计数字滤波器,具有如下几点优点:(1)数字滤波器的单位样值响应完全模仿模拟滤波器的单位样值响应,时域特性逼近较好;(2)两种滤波器的角频率满足ω=ΩTs的线性关系,因此一个线性相位的模拟滤波器可以映射成线性相位的数字滤波器。冲激响应不变法的缺点是只能设计低通和带通等限带的滤波器,而不能设计高通、带阻和全通滤波器。此外,由于模拟滤波器不是真正的带限频谱,所以S平面向Z平面的多值映射会导致频谱混叠等现象,从而导致设计出的滤波器阻带实际的衰减达不到设计的要求。上一页下一页返回9.3IIR数字滤波器设计9.3.2双线性法设计IIR滤波器冲激响应不变法是对模拟滤波器的单位样值响应采样所得到的数字滤波器,其缺点是无法设计非带限滤波器,且容易产生频率响应的混叠失真。为了避免冲激响应不变法的缺点,对于模拟滤波器的数字化方法,人们又提出了双线性变换方法。1.变换原理冲激响应不变法使S平面的jΩ轴每隔2π/Ts便映射到Z平面的单位圆上一周,这种多重映射可能会导致频域混叠现象。为了克服这一问题,需要找到由S平面到Z平面一对一的映射关系,这种关系应保证:(1)S平面的整个jΩ轴仅映射为Z平面单位圆上的一周;上一页下一页返回9.3IIR数字滤波器设计(2)若H(s)是稳定的,由H(s)映射得到的H(z)也应该是稳定的;(3)这种映射是可逆的,既能由H(s)得到H(z),也能由H(z)得到H(s);双线性变换方法满足上述3种条件,其基本思路是:首先将整个S平面压缩到S1

平面内的一条从-π/Ts

到π/Ts的带宽为2π/Ts的横带里,然后再将S1平面内的该横带映射到整个Z平面上去,这样就得到S平面与Z平面间的一一对应的单值关系,如图9-8所示。由此建立S平面与Z平面一一对应的单值关系,消除多值性,也就消除了频率混叠现象。上一页下一页返回9.3IIR数字滤波器设计2.双线性变换法的实现步骤利用双线性变换法设计IIR滤波器可由以下步骤完成:(1)确定数字滤波器指标要求[ωp,δp]和[ωs,δs]。(2)利用Ω=ctan(ω/2)将数字滤波器的技术指标([ωp,δp]和[ωs,δs])变换成模拟滤波器的技术指标([Ωp,δp]和[Ωs,δs]),即将特征频率ωp和ωs转换为Ωp和Ωs,通带最大衰减δp和阻带最小衰减δs不变。(3)根据技术指标[Ωp,δp]、[Ωs,δs]设计模拟滤波器,得到模拟滤波器的系统函数H(s)。(4)利用式(929)将S平面转换到Z平面,得到数字滤波器的系统函数H(z)。上一页下一页返回9.3IIR数字滤波器设计3.双线性变换法的优缺点双线性变换法的优点有:(1)避免了频率响应的混叠现象,但这是以牺牲频率变换的线性关系为代价的。(2)由于双线性变换法无须满足采样定理,因此任何类型的滤波器都可以设计。(3)双线性法利用常数c的选择,可以较为准确地控制滤波器截止频率的位置。双线性变换法的缺点在于设计的数字滤波器无法承袭模拟滤波器的线性相位特性,通常需要在数字滤波器之后通过相位校正网络来获得整体系统的线性相位。双线性变换法和冲激响应不变法的优缺点是相互补充的,在具体的设计过程中可以根据实际需要进行选择上一页下一页返回9.3IIR数字滤波器设计9.3.3频率变换法设计IIR滤波器频率变换法设计IIR滤波器主要包含频率变换和数字化两个过程,根据两个过程的先后次序不同,又分为两种实现方法,即模拟域频率变换法和数字域频率变换法。模拟域频率变换法是先频率变换后数字化,数字域频率变换法则是先数字化后频率变换,两种设计方法都属于原型滤波器设计方法,本节重点介绍模拟域频率变换法。1.模拟域频率变换法的实现步骤模拟域频率变换法的基本思想首先将所要设计的数字滤波器的性能指标按某种频率转换关系,转换成模拟低通滤波器的性能指标,然后利用模拟滤波器设计技术设计出模拟滤波器,最后采用数字化方法(如冲激响应不变法、双线性法等)把模拟滤波器数字化为数字滤波器。上一页下一页返回9.3IIR数字滤波器设计表9--2给出了在模拟域从原型滤波器到其他类型的滤波器的频率变换关系以及平面转换关系。其中Ω代表实际模拟滤波器的角频率,Ω代表归一化模拟原型低通滤波器角频率。s代表实际模拟滤波器系统函数所在的S平面,是jΩ的解析延拓,s代表归一化原型低通模拟滤波器系统函数所在的S平面,是jΩ的解析延拓。2.数字高通滤波器的频率变换法设计设归一化模拟原型低通滤波器Han(jΩ)与模拟高通滤波器HHP(jΩ)的幅度特性如图9-9所示。分析图9-9可得Ω与的对应关系上一页下一页返回9.3IIR数字滤波器设计3.数字带通滤波器的频率变换法设计设归一化模拟原型低通滤波器Han(jΩ)与模拟带通滤波器HBP(jΩ)的幅度特性如图9-10所示。Ωp2和Ωp1分别称为模拟带通滤波器通带的上、下限频率。4.数字带阻滤波器的频率变换法设计设归一化模拟原型低通滤波器Han(jΩ)与模拟带通滤波器HBS(jΩ)的幅度特性如图9-11所示。Ωp2和Ωp1分别称为模拟带阻滤波器通带的上、下限频率。上一页返回9.4FIR数字滤波器设计9.4.1FIR数字滤波器的线性相位特点假设FIR数字滤波器的单位样值响应为h[n],长度为N,则其系统函数为该式表明FIR滤波器在Z平面上有N-1个零点,在原点处有N-1个极点,系统是绝对稳定的。根据系统频率响应的定义,FIR滤波器的频率响应为下一页返回9.4FIR数字滤波器设计1.h(n)为实偶对称情况所谓偶对称是指FIR滤波器的单位样值响应h[n]满足此时称h[n]关于(N-1)/2偶对称。将式(943)代入式(944)可得令m=N-n-1,则有上一页下一页返回9.4FIR数字滤波器设计将式(948)和式(949)两端同时相加,可得因此,FIR滤波器的频率响应为幅度响应:相位响应:群延时:此时FIR滤波器具有线性相位,相位特性如图9-12(a)所示。上一页下一页返回9.4FIR数字滤波器设计2.h(n)为实奇对称情况所谓奇对称是指FIR滤波器的单位样值响应h[n]满足幅度响应:相位响应:群延时:此时FIR滤波器也具有线性相位,相位特性如图9-12(b)所示。上一页下一页返回9.4FIR数字滤波器设计9.4.2FIR数字滤波器的幅度函数特点和零点位置1.幅度函数根据FIR滤波器的单位样值响应h(n)是奇对称还是偶对称,以及h(n)的长度N是奇数还是偶数,可以把上述具有奇偶对称特性的线性相位FIR滤波器分为4种类型。(1)h(n)实偶对称,长度N为奇数(1型FIR滤波器)。FIR滤波器的幅度特性如式(948)所示,当N为奇数时,h(n)和式中的余弦项关于中心点(N-1)/2对称,把式(948)关于中心点合并处理后可得上一页下一页返回9.4FIR数字滤波器设计令m=(N-1)/2-n,则有分析式(959)可知,1型FIR滤波器具有如图9-14所示的时间序列、幅度特性和相位特性。上一页下一页返回9.4FIR数字滤波器设计(2)h(n)实偶对称,长度N为偶数(2型FIR滤波器)。若h(n)实偶对称,长度N为偶数,则使用同样的方法可得2型FIR滤波器的幅度函数为(3)h(n)实奇对称,长度N为奇数(3型FIR滤波器)。若h(n)实奇对称,长度N为奇数,则使用同样的方法可得3型FIR滤波器的幅度函数为上一页下一页返回9.4FIR数字滤波器设计分析式(9-61)可知,3型FIR滤波器具有如图9-15所示的时间序列、幅度特性和相位特性。(4)h(n)实奇对称,长度N为偶数(4型FIR滤波器)。若h(n)实奇对称,长度N为偶数,则使用同样的方法可得4型FIR滤波器的幅度函数为2.零点位置利用式(949)的类似推导过程,当FIR滤波器的单位样值响应是实奇对称或实偶对称时,滤波器的系统函数上一页下一页返回9.4FIR数字滤波器设计9.4.3窗函数法设计FIR滤波器1.设计原理窗函数设计法是从单位样值响应序列着手,使h[n]去逼近理想的单位冲激序列hd[n]。假设理想低通滤波器的频率响应为式中,Hd(ω)为理想滤波器的幅度函数,通常取Hd(ω)=1;ωc

为理想滤波器的截止频率,α为群延迟。则该理想滤波器的单位样值响应为上一页下一页返回9.4FIR数字滤波器设计由式(965)分析可知理想滤波器的单位样值响应无限长非因果序列,关于中心点n=α呈偶对称,波形如图9-17(a)所示。FIR滤波器的单位样值响应h[n]是有限长的,因此需要把无限长度的hd[n]进行截短处理,处理方法就是用一个有限长度的窗函数序列w[n]来截取hd[n]的一段作为h[n],即假设选取矩形窗函数,w(n)的时域波形图如917(b)所示,截断后所得h[n]的波形图如图9-17(c)所示。可见窗函数设计法的基本原理就是将满足滤波器技术要求的无限长单位样值响应加窗截断作为FIR滤波器的单位样值响应。上一页下一页返回9.4FIR数字滤波器设计2.常用窗函数及选取原则为了改善滤波器的性能,通常要求窗函数的幅度特性满足:(1)频率特性的主瓣宽度应尽量窄,且尽可能将能量集中在主瓣内;(2)最大旁瓣相对主瓣值尽可能小,以改善通带的平稳度和增大阻带的衰减;(3)窗函数的相位函数应与理想滤波器的相位函数相同。通常上述几点很难同时满足,在具体设计时需折中处理,下面给出常用的几种窗函数的时域表达式w(n)和频域表达式W(ejω),其幅度图如图9-18所示,基本参数详见表9-3。上一页下一页返回9.4FIR数字滤波器设计3.窗函数法设计步骤用窗函数设计FIR滤波器可分为如下几个步骤。(1)首先由给定滤波器的理想频率响应Hd(ejω),利用傅里叶逆变换,求理想滤波器的单位样值响应,即(2)根据过渡带以及阻带衰减的要求,参考表9--3的基本参数,选择窗函数的形式w[n],并估计窗口的大小N;(3)计算滤波器的单位样值响应h[n]=w[n]hd[n];(4)检验设计出的数字滤波器是否满足技术指标,若不满足所需性能指标,重复步骤(2)至(4),直到满足指标为止。上一页下一页返回9.4FIR数字滤波器设计9.4.4频率采样法设计FIR滤波器频率采样法是从频域出发,把给定的理想频率响应Hd(ejω)等间隔采样N个点得到再对Hd(k)作IDFT,得到为了保证所设计的FIR滤波器具有线性相位,采样频率采样法设计FIR滤波器时需满足如表9-4所列的约束条件。上一页下一页返回9.4FIR数字滤波器设计9.4.5FIR滤波器与IIR滤波器的比较FIR滤波器和IIR滤波器两种滤波器在其自身特点和设计方法上都是各有所长,在此对他们做简单比较,可作为实际设计滤波器时的参考。从滤波器设计来看,IIR滤波器的最大优点是可以利用模拟滤波器的设计技术,而模拟滤波器的设计有大量图表可查,设计技术成熟,所以IIR滤波器的设计比较简单,主要用于设计具有片断常数特性的数字滤波器,但是它在设计过程中往往脱离不了模拟滤波器的模式影响,此外IIR滤波器的不足之处是不容易实现线性相位。FIR滤波器在满足单位样值响应对称性的情况下可以实现严格的线性相位。上一页下一页返回9.4FIR数字滤波器设计从实现成本来看,IIR滤波器可以用较低的阶数实现较高的技术指标,设计成本较低。对于同样的滤波器设计指标,FIR滤波器要求的阶数通常要比IIR滤波器的阶数高5至10倍,因此FIR滤波器的设计成本较高。从滤波器结构上看,IIR滤波器一般采用递归结构,为使系统稳定,极点必须全部位于单位圆内,这种结构在有限字长处理是容易引起较大误差。FIR滤波器一般不采用递归结构,总是稳定的系统,因此误差很小。从运算角度看,FIR滤波器可以利用FFT算法实现,运算速度很快,运算量很小。IIR滤波器则不可能这样计算。总之,FIR滤波器和IIR滤波器各有所长,在实际运用时需多方面考虑。在对相位特性要求不高的场合可以偏向于选择IIR滤波器,对于要求线性相位的场合则最后选择FIR滤波器。上一页返回9.5本章相关的MATLAB命令及应用9.5.1IIR数字滤波器设计【例96】设计一个butterworth低通滤波器,通带截止频率为100Hz,阻带截止频率为200Hz,通带衰减rp小于2dB,阻带衰减rs大于15dB,取样频率为fs=1500Hz。(1)请用冲激响应不变法设计,并作出其频响特性图;(2)请用双线性法设计,并作出其频响特性图;1)冲激响应不变法程序clc;clear;closefp=100;fss=200;fs=1500;rp=2;rs=15;wp1=2*pi*fp/fs;%求出通带数字角频率wp1wss1=2*pi*fss/fs;%求出阻带数字角频率wss1下一页返回9.5本章相关的MATLAB命令及应用[N,wc]=buttord(2*pi*fp,2*pi*fss,rp,rs,'s')%求出最小阶数N和截止频率wc[z,p,k]=buttap(N);[ba,aa]=zp2tf(z,p,k)%求出归一化模拟滤波器的分子分母系数ba和aa[b,a]=lp2lp(ba,aa,wc)%求出截止频率为wc的模拟滤波器的分子分母系数b和a[bz,az]=impinvar(b,a,fs)[H,f]=freqz(bz,az);ma=20*log10(abs(H));mp=angle(H)*180/pi;subplot(2,1,1);plot(f/pi,ma)subplot(2,1,2);plot(f/pi,mp)上一页下一页返回9.5本章相关的MATLAB命令及应用2)双线性法程序clc;clear;closefp=100;fss=200;fs=1500;rp=2;rs=15;wp1=2*pi*fp/fs;%求出通带数字角频率wp1wss1=2*pi*fss/fs;%求出阻带数字角频率wss1op1=2*fs*tan(wp1/2);os1=2*fs*tan(wss1/2);%预畸变处理,求出模拟角频率op1和os1[N,wc]=buttord(op1,os1,rp,rs,'s')%求出最小阶数N和截止频率wc[z,p,k]=buttap(N);[ba,aa]=zp2tf(z,p,k)%求出归一化模拟滤波器的分子分母系数ba和aa上一页下一页返回9.5本章相关的MATLAB命令及应用[b,a]=lp2lp(ba,aa,wc)%求出截止频率为wc的模拟滤波器的分子分母系数b和a[bz,az]=bilinear(b,a,fs)[H,f]=freqz(bz,az);ma=20*log10(abs(H));mp=angle(H)*180/pi;subplot(2,1,1);plot(f/pi,ma)subplot(2,1,2);plot(f/pi,mp)9.5.2FIR数字滤波器设计【例97】用窗函数法设计一个线性相位的FIR低通滤波器,指标为:上下边带截止频率为w1=0.2π,w2=0.3π,阻带衰减不小于40dB,通带衰减不大于3dB。上一页下一页返回9.5本章相关的MATLAB命令及应用程序:clc;clear;closewp=0.2*pi;ws=0.3*pi;detaw=ws-wp;%由rs确定选用汉宁窗n=ceil(3.1*2*pi/detaw);wc=(wp+ws)/2;b=fir1(n-1,wc/pi,hanning(n))freqz(b,1,512)上一页下一页返回9.5本章相关的MATLAB命令及应用9.5.3F