|
第7章 有限脉冲响应数字滤波器的设计 |
|
|
第7章有限脉冲响应数字滤波器的设计本章主要内容线性相位FIR滤波器的条件和特点利用窗函数法设计FIR滤波器利用频域采样法设计FIR滤波器 IIR与FIR数字滤波器的比较FIR数字滤波器的Matlab仿真实现引言假设N为FIR数字滤波器的单位脉冲响应h(n)的长 度,那么其系统函数H(z)为:IIR数字滤波器设计过程中只考虑了幅频特性,没有考虑相位特性,所设计的滤波器相位特性一般是非线性的。 为了得到线性相位特性,则要采用全通网络进行相位校正,从而使得设计复杂,成本变高,也难以得到严格的线性相位。引言FIR数字滤波器很 容易得到严格的线性相位。FIR数字滤波器的单位脉冲响应是有限长的,因此总是稳定的。FIR滤波器的设计方法:窗函数法频率取样法等纹波 逼近法7.1线性相位FIR数字滤波器的条件及特点1.线性相位FIR数字滤波器单位脉冲响应h(n)的长度为N,则其频率响应函数 为其中,为幅度特性,为相位特性。而线性相位FIR滤波器是指是w的线性函数。其中:第一类线性相 位的相位特性函数是ω的严格线性函数:第二类线性相位FIRDF的相位特性函数如下式:7.1线性相位FIR数字滤波器2. 线性相位FIR滤波器的时域约束条件1)第一类线性相位对h(n)的约束条件第一类线性相位FIRDF的相位函数,由(7.1.1)和( 7.1.2)得由式(7.1.5)得到7.1线性相位FIR数字滤波器将式(7.1.6)中两式相除得:即移项并用三角公式化 简得到要求和h(n)满足如下条件:7.1线性相位FIR数字滤波器2)第二类线性相位对h(n)的约束条件第二类线性相位 FIRDF的相位函数由式(7.1.1)和式(7.1.2),有:经过同样的推导过程可得到:要求和h(n)满足如下条件7.1线性 相位FIR数字滤波器3.线性相位FIR滤波器幅度特性Hg(ω)的特点对幅度特性Hg(ω)的约束情况1:h(n)=h(N-n- 1),N为奇数时Hg(ω)的特点:Hg(ω)关于0,π,2π三点偶对称,可以实现各种(低通、高通、带通、带阻)滤波器。情况2: h(n)=h(N-n-1),N为偶数时Hg(ω)的特点:Hg(ω)关于0,2π偶对称,而关于π奇对称,因此不能实现高通和带阻滤波 器。7.1线性相位FIR数字滤波器情况3:h(n)=-h(N-n-1),N为奇数时Hg(ω)的特点:Hg(ω)关于0,π,2 π三点奇对称,只能实现带通滤波器。情况4:h(n)=-h(N-n-1),N为偶数时Hg(ω)的特点:Hg(ω)关于0,2π奇对 称,而关于π偶对称,不能实现低通和带阻滤波器。线性相位FIRDF的时域和频域特性总结表7.1.1(a)线性相位FIRDF的时域 和频域特性一览表线性相位FIRDF的时域和频域特性总结表7.1.1(b)线性相位FIRDF的时域和频域特性一览表线性相位FIR DF的零点分布特点4.线性相位FIRDF的零点分布特点将代入上式,得到:由(7.1.1 4)式可以看出:如果是的零点,其倒数也必然是其零点;又因为h(n)是实序列,H(z)的零点必定共轭成 对,因此和也是其零点。这样线性相位FIR滤波器的零点确定其中一个,另外三个也确定了,当然也有特殊情况(共轭为本身 or倒数为本身的情况,如图所示)7.2利用窗函数法设计FIR滤波器设计思想寻找一个FIR滤波器,使其频率响应H(ejω)逼近理 想FIR滤波器的频率响应Hd(ejω).其单位脉冲响应用hd(n)表示。一般情况下,Hd(ejω)在边界频率处有不连续点,因此hd (n)是无限长的,且是非因果的序列。7.2利用窗函数法设计FIR滤波器理想低通滤波器的频率响应其对应的单位脉冲响应可以看到 hd(n)是一个无限长,非因果的序列。7.2利用窗函数法设计FIR滤波器将hd(n)截取长度为N的一段,构成h(n)为了保证设计 的滤波器具有第一类线性相位,必须满足截取的一段关于n=(N-1)/2偶对称可以将h(n)看作是hd(n)与矩形窗相乘 窗函数设计法的频域解释时域加窗频域卷积窗函数时域表示频域表示幅度相位FIR滤波器的幅频特性FIR滤波器的幅频特性理想低通与矩形窗频 谱函数卷积过程加窗对Hd(ω)的影响(1)在理想特性不连续点ωc附近形成过渡带。过滤带的宽度近似等于WR(θ)主瓣宽度,Δω= 4π/N。(2)通带内增加了波动,最大的峰值在ωc-2π/N处。阻带内产生了余振,最大的负峰在ωc+2π/N处。通带与阻带 中波动的情况与窗函数的幅度谱有关。WR(θ)波动愈快(加大时),通带与阻带内波动愈快,WR(θ)旁瓣的大小直接影响波动的大小。 这些影响是对hd(n)加矩形窗引起的,称之为吉布斯效应。减小吉布斯效应的方法增加矩形窗口的宽度N不能减少吉布斯效应的影响。N 的改变只能改变ω坐标的比例和的绝对大小,不能改变主瓣和旁瓣幅度相对值。加大N并不是减少吉布斯效应的有效方法。 寻找合适的窗函数形状,使其谱函数的主瓣包含更多的能量,相应旁瓣幅度就变小了;旁瓣的减少可使通带与阻带波动减少,从而加大阻带的衰减。 但这样总是以加宽过渡带为代价的。几种常用的窗函数矩形窗三角(Bartlett)窗汉宁(Hanning)窗哈明(Hamming) 窗布莱克曼(Blackman)窗凯塞-贝塞尔(Kaiser-Basel-)窗几种常用的窗函数矩形窗图7.2.3矩形窗的四种 波形三角(Bartlett)窗三角(Bartlett)窗其频谱函数为:其幅度函数为:汉宁(Hanning)窗汉宁(Hanning )窗哈明(Hamming)窗哈明(Hamming)窗布莱克曼(Blackman)窗布莱克曼(Blackman)窗凯塞-贝塞 尔(Kaiser-Basel)窗凯塞-贝塞尔(Kaiser-Basel)窗常用窗函数的波形常用窗函数的频谱6种窗函数的基本参数 窗函数类型旁瓣峰值?n(dB)过渡带宽度?B阻带最小衰减?s(dB)近似值精确值矩形窗-134?/N1.8?/N-21三角窗-25 8?/N6.1?/N-25汉宁窗-318?/N6.2?/N-44哈明窗-418?/N6.6?/N-53布莱克曼窗-5712?/N1 1?/N-74凯塞窗(?=7.865)-57?10?/N-80表中过渡带宽和阻带最小衰减是用对应的窗函数设计的FIR数字滤波器的 频率相响应指标。除了以上6种窗函数外,比较有名的还有Chebyshev窗、Gaussian窗。7.2.3用窗函数法设计FIR滤波器 的步骤1.根据允许的过渡带宽及阻带衰减,选定窗函数和N值。2.给出希望设计的滤波器的频率响应函数Hd(ejω)。3.计算hd(n )如果Hd(ejω)不能用简单函数表示,可以用求和代替积分。用窗函数法设计FIR滤波器的步骤4.将hd(n)与窗函数相乘得F IR数字滤波器的冲激响应h(n)5.计算FIR数字滤波器的频率响应,并验证是否达到所要求的指标例窗函数法设计FIR滤波器例 7.2.1用窗函数法设计线性相位高通FIRDF,要求通带截止频率ωp=?/2rad,阻带截止频率ωs=?/4rad ,通带最大衰减?p=1dB,阻带最小衰减?s=40dB。解:(1)选择窗函数w(n),计算窗函数长度N。已知阻带最小衰 减?s=40dB,由表(7.2.2)可知汉宁窗和哈明窗均满足要求,这里选择汉宁窗。本例中过渡带宽度,汉宁窗的精确过渡带宽度 ,解不等式得到:N≥24.8,而对高通滤波器N要取奇数,所以此处N取25.所以有:例窗函数法设计FIR滤波器(2 )构造式中,(3)求出将τ=12代入,得:公式中第一项对应全通滤波器,第二项是低通,两者之差即为高通。例窗函数法设计FIR滤波 器(4)加窗窗函数的MATLAB设计函数简介实际设计的时候一般用Matlab工具箱函数,可调用工具箱函数fir1实现窗函数法设计步 骤的(2)~(4)的过程。1)fir1是用窗函数法设计线性相位FIRDF的工具箱函数,以实现线性相位FIRDF的标准窗函数法设计。 2)用fir2函数设计时,可以指定任意形状的,所以称之为任意形状幅度特性窗函数法设计函数,它实质是一种频率采样法与窗函数法的综 合设计函数7.3利用频域采样法设计FIR滤波器1.基本思想频域采样法是在频率域对理想滤波器Hd(ejω)采样,在采样点上设计的 滤波器H(ejω)和理想滤波器Hd(ejω)幅度值相等,然后根据频率域的采样值求得实际设计的滤波器的频率特性H(ejω)。对理想滤 波器的频率特性Hd(ejω)在[0,2?]范围内等间隔地取样N个点7.3利用频域采样法设计FIR滤波器根据插值公式此式就是直 接利用频率采样值Hd(k)形成滤波器的系统函数,适合频率采样结构2.设计线性相位滤波器时对Hd(k)的约束条件当h(n)是偶对称, h(n)=h(N-n-1),N为奇数。H(ω)是偶对称的对H(ejω)在[0,2?]范围内等间隔地取样N个点2.设计线性相位滤波器 时对Hd(k)的约束条件当h(n)是偶对称,h(n)=h(N-n-1),N为偶数。H(ω)是偶对称的对H(ejω)在[0,2?]范 围内等间隔地取样N个点3.逼近误差及其改进措施采用频域取样法设计的FIR数字滤波器在阻带内的衰减很小,在实际应用中往往达不到要求。 产生这种现象的原因是由于在通带边缘采样点的陡然变化而引起的起伏振荡。增加阻带衰减的方法是在通带和阻带的边界处增加一些过渡的采样点, 从而减小频带边缘的突变,也就减小了起伏振荡,增大了阻带最小衰减。3.逼近误差及其改进措施的内插表示形式:式中因此采样点处H(e jwk)与H(k)相等,逼近误差为03.逼近误差及其改进措施图7.3.2采样点数N不同时的逼近误差比较(N=15,75)图7 .3.1频率采样法设计过程中的波形(N=15)3.逼近误差及其改进措施将过渡带采样点的个数m与滤波器阻带最小衰减?s的经验数 据列于表7.3.1中,我们可以根据给定的阻带最小衰减?s选择过渡带采样点的个数m。表7.3.1过渡带采样点的个数m与滤波器阻带 最小衰减?s的经验数据m123?s44~54dB65~75dB85~95dB4.频率采样法设计步骤(1)根据阻带最小衰减 选择过渡带采样点的个数m。(2)确定过渡带宽度,估算频域采样点数(即滤波器长度)N。如果增加m个过渡带采样点,则过渡带 宽度近似变成。当N确定时,m越大,过渡带越宽。如果给定过渡带宽度,则要求 ,即滤波器长度N必须满足以下等式:(3)构造一个希望逼近的频率响应函数:设计标准型片段常数特性的FIR 数字滤波器时,一般构造幅度特性函数为相应的理想频响特性,且满足表7.1.1要求的对称性。4.频率采样法设计步骤( 4)按照(7.3.1)式进行频域采样:并加入过渡采样。过渡带采样值可以设置为经验值,或用累试法确定,也可以采用优化方法估算。(5 )对进行N点IDFT,得到第一类线性相位FIR数字滤波器的单位脉冲响应:(6)检验设计结果。如果阻带最小衰减未达到指 标要求,则要改变过渡带采样值,直到满足指标要求为止。如果滤波器边界频率未达到指标要求,则要微调的边界频率。7. 4利用等波纹最佳逼近法设计FIR数字滤波器等波纹最佳逼近法的特点:(1)是一种优化设计法,它克服了窗函数设计法和频率采样法的缺 点,使最大误差(即波纹的峰值)最小化,并在整个逼近频段上均匀分布。即幅频响应在通带和阻带都是等波纹的;(2)可以分别控制通带和阻 带波纹幅度。(3)在滤波器长度给定的条件下,使加权误差波纹幅度最小化。(4)与窗函数设计法和频率采样法比较,设计的滤波器性能价 格比最高。阶数相同时,这种设计法使滤波器的最大逼近误差最小,即通带最大衰减最小,阻带最小衰减最大;指标相同时,这种设计法使滤波器阶 数最低。比较例题7.2.3和例题7.2.1,设计指标要求相同的带阻滤波器,等波纹最佳逼近法需要28阶,而窗函数法需要79阶!7. 4.1等波纹最佳逼近法的基本思想Hd(ω)——表示希望逼近的幅度特性函数。Hg(ω)——表示实际设计的滤波器幅度特性函数 。——加权误差函数W(ω)——称为误差加权函数,用来控制不同频段(一般指通带和阻带)的逼近精度。等波纹最佳逼近基于切比雪夫逼 近,在通带和阻带以的最大值最小化为准则,采用Remez多重交换迭代算法求解滤波器系数h(n)。W(ω)取值越大的频段逼近精度越 高,开始设计时应根据逼近精度要求确定W(ω),在Remez多重交换迭代过程中W(ω)是确知函数。由于切比雪夫(Chebyshev )和雷米兹(Remez)对解决该问题做出了贡献,所以又称之为切比雪夫逼近法,或雷米兹逼近法。7.4.1等波纹最佳逼近法的基本 思想等波纹最佳逼近设计中,把数字频段分为“逼近(或研究)区域”——一般指通带和阻带“无关区域”——一般指过渡带设计过程中 只考虑对逼近区域的最佳逼近。注意:无关区宽度不能为零,即Hd(ω)不能是理想滤波特性。否则,设计必将失败。7.4.1等波 纹最佳逼近法的基本思想等波纹滤波器的技术指标及其描述参数:等波纹最佳逼近设计法的数学证明复杂,所以本节略去其复杂的数学推导,只介 绍其基本思想和实现线性相位FIR数字滤波器的等波纹最佳逼近设计的MATLAB信号处理工具箱函数remez和remezord。[M, fo,mo,w]=remezord(f,m,rip,Fs)hn=remez(M,fo,mo,w)求滤波器阶数N和 误差加权函数W(ω)时,要求给出?1和?2工程实际中设计指标通常以?p和?s给出两种描述参数之间的换算关系:7.4.1等波纹最 佳逼近法的基本思想期望逼近的通带为[0,?/4],阻带为[5?/16,?],过渡带已知由图可见:和N由滤波器 设计指标(即?p和?s,以及过渡带宽度)确定等波纹最佳逼近法设计FIR数字滤波器的过程:(1)根据给定的逼近指标估算滤波器阶数 N和误差加权函数;(2)采用remez算法得到滤波器单位脉冲响应h(n)。MATLAB工具箱函数remezo rd和remez就是完成以上2个设计步骤的有效函数。7.4.2remez和remezord函数介绍1.Remez:采用re mez算法实现线性相位FIR数字滤波器的等波纹最佳逼近设计。其调用格式为:hn=remez(M,f,m,w)返回单位脉冲响 应向量hn。remez函数的调用参数(M,f,m,w)一般通过调用remezord函数来计算。调用参数含义如下:M为 FIR数字滤波器阶数,hn长度N=M+1。f和m给出希望逼近的幅度特性:f为边界频率向量,0≤f≤1,要求f为单调增向量(即 f(k)幅度向量,m与f长度相等,m(k)表示频点f(k)的幅度响应值。如果用命令Plot(f,m)画出幅频响应曲线,则k为奇数时,频段 [f(k),f(k+1)]上的幅频响应就是期望逼近的幅频响应值,频段[f(k+1),f(k+2)]为无关区。简言之,Plot(f, m)命令画出的幅频响应曲线中,起始频段为第一段,奇数频段为逼近区,偶数频段为无关区。7.4.2remez和remezord函 数介绍例如,通带频段,阻带频段,则,f=[0,1/4,5/1 6,1];m=[1,1,0,0]plot(f,m)画出的幅度特性曲线如图所示。图中奇数段(第1、3段)的水平幅度为希望 逼近的幅度特性,偶数段(第2段)的下降斜线为无关部分,逼近时形成过渡带,并不考虑该频段的幅频响应形状。w为误差加权向量,其长度为 f的一半。w(i)表示对m中第i个逼近频段的误差加权值。图7.4.2(a)中w=[1,10]。缺省w时,默认w为全1(即每个逼近 频段的误差加权值相同)。希望逼近的幅度特性曲线remez函数还可以设计两种特殊滤波器:希尔伯特变换器和数字微分器。调用格式:hn= remez(M,f,m,w,''hilbert'')hn=remez(M,f,m,w,''defferentiator'' )7.4.2remez和remezord函数介绍2.Remezord:根据逼近指标估算FIR数字滤波器的最低阶数M、误差加权向 量w和归一化边界频率向量f。返回参数作为remez函数的调用参数。调用格式:[M,fo,mo,w]=remezord(f ,m,rip,Fs)hn=remez(M,fo,mo,w)参数说明:f与remez中类似,这里f可以是模拟频率(Hz )或归一化数字频率,但必须以0开始,以Fs/2(用归一化频率时对应1)结束,而且省略了0和Fs/2两个频点。Fs为采样频率,缺省时 默认Fs=2Hz。注意:这里f的长度(包括省略的0和Fs/2两个频点)是m的两倍,即m中的每个元素表示f给定的一个逼近频段上希望逼 近的幅度值。例如,对上面希望逼近的幅度特性曲线图,f=[1/4,5/16];m=[1,0]。7.4.2remez和r emezord函数介绍注意:①省略Fs时,f中必须为归一化频率。②有时估算的阶数M略小,使设计结果达不到指标要求,这时要取M+ 1或M+2(必须注意对滤波器长度N=M+1的奇偶性要求)。所以必须检验设计结果。③如果无关区(过渡带)太窄,或截止频率太接近零 频率和Fs/2时,设计结果可能不正确。rip表示f和m描述的各逼近频段允许的波纹幅度(幅频响应最大逼近误差),f的长度是rip的 两倍。例:对上图所示的低通滤波器,对高通滤波器,remezord函数的调用参数f,m,rip和Fs调用remez和reme zord函数设计线性相位FIR数字滤波器的关键问题:设计指标remezord函数的调用参数f,m,rip和Fs其中,Fs一般 是根据实际信号处理要求(按照采样定理)确定。下面给出由各种滤波器设计指标确定remezord调用参数f,m和rip的公式,编程 时直接套用即可。1、低通滤波器设计指标逼近通带:,通带最大衰减:逼近阻带:,阻带最小衰减:remezord调用参数: remezord函数的调用参数f,m,rip和Fs2、高通滤波器设计指标 逼近通带:,通带最大衰减:;逼近阻带:,阻带最小衰 减:。调用参数:3、带通滤波器设计指标 逼近通带:,通带最大衰减:,逼近阻带:,阻带最小衰减:,remezord调用参 数:,4、带阻滤波器设计指标 ,,通带最小衰减:逼近通带:,阻带最大衰减:逼近阻带:remezord调用参数:注意:工程 实际中常常给出对模拟信号的滤波指标要求,设计数字滤波器,对输入模拟信号采样后进行数字滤波。这时,调用参数f可以用模拟频率表示。但 是,调用remezord时一定要加入采样频率参数Fs。等波纹最佳逼近法设计举例例7.4.1利用等波纹最佳逼近法重新设计FIR带 阻通滤波器。指标与例7.2.3相同,即: 逼近通带:,通带最大衰减: 逼近阻带:,阻带最小衰减:解:调用remezord和 remez函数求解。本例设计程序为ep741.m。%ep741.m:例7.4.1用remez函数设计带阻滤波器f=[0.2, 0.35,0.65,0.8]; %省略了0和1m=[1,0,1];rp=1;rs=60;%由式(7.4.4)和式(7.4. 5)求通带和阻带波纹幅度dat1、dat2和rip:dat1=(10^(rp/20)-1)/(10^(rp/20)+1);da t2=10^(-rs/20);rip=[dat1,dat2,dat1];[M,fo,mo,w]=remezord(f ,m,rip);hn=remez(M,fo,mo,w);%以下绘图检验部分略去等波纹最佳逼近法设计举例程序运行结果:M =28。即h(n)的长度N=29。注意:例7.2.3中N=80。等波纹最佳逼近法设计举例例7.4.2利用等波纹最佳逼近法设计 FIR数字低通滤波器,实现对模拟信号的数字滤波处理。要求与例7.2.2相同。求出h(n),并画出损耗函数曲线和相频特性曲线。解: 将例7.2.2中所给指标重写如下:通带截止频率:fp=1500Hz;通带最大衰减:αp=1dB;阻带截止频率:fs =2500Hz;阻带最小衰减:αs=40dB;对模拟信号的采样频率:Fs=10kHz。调用remezord和remez 函数设计的程序ep742.m在下页给出等波纹最佳逼近法设计举例%ep742.m:例7.4.2用remez函数设计低通滤波器F s=10000;%对模拟信号采样的频率为1kHzf=[1500,2500]; %边界频率为模拟频率(Hz)m=[1,0]; rp=1;rs=40;dat1=(10^(rp/20)-1)/(10^(rp/20)+1);dat2=10^(-rs/20 );rip=[dat1,dat2];[M,fo,mo,w]=remezord(f,m,rip,Fs);M=M+1; %边界频率为模拟频率(Hz)时,调用参数必须加入采样频率Fshn=remez(M,fo,mo,w);%以下绘图检验部分省略等 波纹最佳逼近法设计举例运行结果:滤波器阶数M=15。例7.2.2中用窗函数设计的滤波器阶数为23。7.5IIR与FIR数字 滤波器的比较首先,从性能上来说,IIR滤波器系统函数的极点可以位于单位圆内的任何地方,因此可用较低的阶数获得好的选择性,但是这个高效率是以相位的非线性为代价的。选择性越好,则相位非线性越严重。相反,FIR滤波器却可以得到严格的线性相位特性。然而由于FIR滤波器系统函数的极点固定在原点,所以只能用较高的阶数达到高的选择性;对于同样的滤波器幅频响应指标,FIR滤波器所要求的阶数可以比IIR滤波器高5~10倍,成本较高,运算量大,信号延时也较大。对相同的选择性和相同的线性相位要求来说,则IIR滤波器必须加全通网络进行相位校正,同样要大大增加滤波器的阶数和复杂性。7.5IIR与FIR数字滤波器的比较从结构上,IIR滤波器必须采用递归结构,极点位置必须在单位圆内,否则系统将不稳定。相反,FIR滤波器主要采用非递归结构,不论在理论上还是在实际的有限精度运算中都不存在稳定性问题。此外,FIR滤波器可以采用快速傅里叶变换算法实现,在相同阶数的条件下,运算速度可以大大提高。7.5IIR与FIR数字滤波器的比较从设计工具看,IIR滤波器可以借助于模拟滤波器的成果,因此一般都提供有效的封闭形式的设计公式进行准确计算,计算工作量比较小,对计算工具的要求不高。FIR滤波器设计则一般没有封闭形式的设计公式。另外,也应看到,IIR滤波器虽然设计简单,但主要用于设计具有片段常数特性的滤波器,如低通、高通、带通及带阻滤波器等,往往脱离不了模拟滤波器的局限性,而FIR滤波器则要灵活得多。 |
|
|
|
|
|
|
|
|
|
|