数字信号处理曹成茂无限长单位脉冲响应iir滤波器设计方法

1、第四章 无限长单位脉冲响应（IIR）滤波器设计方法4.1滤波器的基本原理 NiiMiiinybinxany00系统函数为： NiiiMiiizbzazH001bi0,i0,该系统为IIR系统。当所有bi=0, i0，该系统为FIR系统4.1.1滤波器的分类模拟滤波器数字滤波器rcjHAtjHjH,lg201lg201lg20,1lg2011lg20lg2022max11minmax4.1.2 可实现滤波器的特性通带波动最小阻带衰减At1maxjH式中，通带内2maxjH阻带内设计方法： 先设计一个合适的模拟滤波器，然后变换成满足预定指标的数字滤波器。 由于模拟的网络综合理论已经发展得很成熟模拟

2、滤波器有简单而严格的设计公式，设计起来方便、准确、可将这些理论推广到数字域，作为设计数字滤波器的工具。最优化设计方法分两步： a) 确定一种最优准则，如最小均方误差准则，即使设计出的实际频率响应的幅度特性 （与所要求的理想频率响应 的均方误差最小，| )(|jeH| )(|jdeHmin)()(21MijdjiieHeHiaibiaib此外还有其他多种误差最小准则， b) 在此最佳准则下，求滤波的系数 和 通过不断地迭代运算，改变 、 ，直到 满足要求为止。 滤波器的技术指标给定后，需要设计一个传输函数Ha(s)，希望其幅度平方函数满足给定的指标，一般滤波器的单位冲激响应为实数，因此2()(

3、)()()()aaasjaaHjHs HsHjHj s平面与z 平面间的映射TereererezjsTTjTjjsTez s平面的虚轴映射z平面的单位圆（r=1） ，s平面的左半平面映射z平面单位圆内（r1）。由 确定 的方法如下：（1）由 得到象限对称的s平面函数；（2）将 因式分解，得到各零极点，将左半平面极点归于 。 轴上的零点或者极点都为偶次，应取一半（应为共轭对）作为 的零点或极点。 （3）按照 与 的低频或高频特性的对比就可以确定出增益常数。 （4）由求出的零点，极点及增益常数，则可完全确定系统函数2)()()(jHsHsHajsaa2)( jHa)(sHa2)( jHa)(sHa

4、)(sHa)( jHa)(sHaj)(sHa例 根据以下幅度平方函数 确定系统函数2)( jHa)(sHa)36)(49()25(16)(22222jHa解：)36)(49()25(16)()()(22222sssjHsHsHajsaa)6)(7()25()(Hj5s67s20asssks，零点为，其极点为4)()(000kjHsHasa可以确定Step1:Step2:Step3:42131004)6)(7()25(4)(H222asssssss4.2 模拟滤波器设计方法 为了方便学习数字滤波器，先讨论几种常用的模拟低通滤波器设计方法，高通、带通 、带阻等模拟滤波器可利用变量变换方法，由低通滤

5、波器变换得到。 模拟滤波器的设计就是根据一组设计规范设计模拟系统函数Ha(s)，使其逼近某个理想滤波器特性。 因果系统中式中ha(t)为系统的冲激响应，是实函数。 不难看出 0)()(dtethjHtjaa0sincos)()(dttjtthjHaa)()(jHjHaa定义振幅平方函数 式中 Ha(s)模拟滤波器 系统函数 Ha(j)滤波器的频率响应 |Ha(j)|滤波器的幅频响应又 S=j，2=-S2 A(2)=A(-S2)|S=j) 1 ()()()()()()()()()(222jsaaaaaaasHsHjHjHAjHjHjHA 问题：由A(-S2)Ha(S) 对于给定的A(-S2)，先

6、在S复平面上标出A(-S2)的极点和零点，由(1)式知，A(-S2)的极点和零点总是“成对出现”，且对称于S平面的实轴和虚轴，选用A(-S2)的对称极、零点的任一半作为Ha(s)的极、零点，则可得到Ha(s)。 为了保证Ha(s)的稳定性，应选用 A(-S2)在S左半平面的极点作为Ha(s)的极点，零点可选用任一半。NcajjjHA22211)()()(2AN为滤波器阶数， 如图1其幅度平方函数：特点：具有通带内最大平坦的振幅特性，且随f ，幅频特 性 单调 。三种模拟低通滤波器的设计：4.2.1巴特沃兹滤波器巴特沃兹滤波器 (Butterworth 滤波器) (巴特沃兹逼近)图1 巴特沃兹滤

7、波器 振幅平方函数 通带: 使信号通过的频带 阻带：抑制噪声通过的频带 过渡带：通带到阻带间过渡的频率范围 c ：通带边界频率。 过渡带为零， 阻带|H(j)|=0 通带内幅度|H(j)|=const.， H(j)的相位是线性的。 理想滤波器 图1中，N增加，通带和阻带的近似性越好，过渡带越陡。 在过渡带内，阶次为的巴特沃兹滤波器的幅度响应趋于斜率为-6NdB/倍频程的渐近线。 通带内，分母/c1， ( /c)2N 1, 增加， A(2)快速减小。 =c, ， ,幅度衰减 ， 相当于3dB衰减点。 21)0()(2AAc2121)(2A 振幅平方函数的极点： 令分母为零，得 可见，Butter

8、worth滤波器 的振幅平方函数有2N个极点，它们均匀对称地分布在|S|=c的圆周上。 例：为N=3阶BF振幅平方函数的极点分布，如图2。NcaajSSHSH2)(11)()()()1(21cNPjS图2 三阶A(-S2)的极点分布 考虑到系统的稳定性，知DF的系统函数是由S平面左半部分的极点（SP3，SP4，SP5）组成的，它们分别为： 系统函数为： 令 ，得归一化的三阶BF： 反归一化，则有3254323,jcpcpjcpeSSeS)()()(5433pppcaSSSSSSsH1c1221)(23SSSsHa1)/(2)/(2)/(1)(23cccassssHMATLAB设计模拟设计模拟B

9、utterworth filterz,p,k=buttap(N)确定确定N阶归一化的阶归一化的Butterworth filter 的零点、极点和增益的零点、极点和增益(gain)( )(1)(2)()kH sspspsp Nnum,den=butter(N,Wc,s)确定阶数为确定阶数为N，3-db截频为截频为Wc(radian/s)的的Butterworth filter分子和分母多项式。分子和分母多项式。s 表示模拟域。表示模拟域。N,Wc=buttord(Wp,Ws,Ap,As,s)确定模拟确定模拟Butterworth filter的阶数的阶数N和和3-db截频截频Wc。Wc是由阻带

10、参数确定的。是由阻带参数确定的。 例：设计满足下列条件的模拟例：设计满足下列条件的模拟Butterworth低通滤波器低通滤波器 fp=1kHz, fs=5kHz, Ap=1dB, As=40dBWp=2*pi*1000;Ws=2*pi*5000;Ap=1;As=40;N,Wc=buttord(Wp,Ws,Ap,As,s);num,den = butter(N,Wc,s);omega = 0: 200: 12000*pi;h = freqs(num,den,omega);gain = 20*log10(abs(h);plot (omega/(2*pi),gain);xlabel(Frequen

11、cy in Hz); ylabel(Gain in dB);0100020003000400050006000-50-40-30-20-100Frequency in HzGain in dBN=4 Ap= 0.1098 dB As= 40.0000 dB4.2.2切比雪夫（Chebyshev）滤波器 (切比雪夫多项式逼近) 特点：误差值在规定的频段上等幅变化。 巴特沃兹滤波器在通带内幅度特性是单调下降的，如果阶次一定，则在靠近截止频率 处，幅度下降很多，或者说，为了使通内的衰减足够小，需要的阶次（N）很高，为了克服这一缺点，采用切比雪夫多项式逼近所希望的 。 切比雪夫滤波器的 在通带范围内是

12、等幅起伏的，所以同样的通带衰减，其阶数较巴特沃兹滤波器要小。可根据需要对通带内允许的衰减量（波动范围）提出要求，如要求波动范围小于1db。 c2()aHj2()aHj 振幅平方函数为)(11)()(2222cNaVjHAc1)coshcosh(1)coscos()(11xxNxxNxVN)(, 11)(,1xVxxxVxNN时有效通带截止频率与通带波纹有关的参量， 大 ，波纹大。 0 1 VN（x）N阶切比雪夫多项式，定义为 如图,通带内 变化范围1 c ，随/c ， 0 (迅速趋于零)当 =0时， N为偶数， ，min ， N为奇数， ， max，)2(cos11)0arccos(cos11

13、)(22220NNjHa22011)(jHa1)(02jHa1)2(cos2N0)2(cos2N2111c2a)( jH2a)( jH 切比雪夫滤波器的振幅平方特性 2minmax111lg20)()(lg20jHjHaa)1lg(1021101.02)(dB221)(,AjHar时 给定通带波纹值分贝数 后，可求 。有关参数的确定: a、通带截止频率c ，预先给定 b、通带波纹为 c、阶数N由阻带的边界条件确定。（ 、A事先给定） 222111AVcrNr221)(,AjHar时)/cosh(/1cosh()coshcosh()(,12crNarAarNxNarxVx得时MATLAB设计模拟

14、设计模拟type I Chebyshev filterz,p,k=cheb1ap(N,Ap);确定确定N阶归一化的阶归一化的Chebyshev filter 的零点、极点和增益的零点、极点和增益(gain)num,den=cheby1(N,Ap,Wc,s)确定阶数为确定阶数为N，通带截频为，通带截频为Wc(radian/s)的的Chebyshev filter。s 表示模拟域表示模拟域N,Wc=cheb1ord(Wp,Ws,Ap,As,s)确定模拟确定模拟Chebyshev filter的阶数的阶数N。Wc=Wp(rad/s)例：设计满足下列条件的模拟例：设计满足下列条件的模拟CB I型低通滤

15、波器型低通滤波器 fp=1KHz, fs=5kHz, Ap=1dB, As=40dBWp=2*pi*1000;Ws=2*pi*5000;Ap=1;As=40;N,Wc=cheb1ord(Wp,Ws,Ap,As,s);num,den = cheby1(N,Ap,Wc,s);omega=Wp Ws; h = freqs(num,den,omega); fprintf(Ap= %.4fn,-20*log10(abs(h(1); fprintf(As= %.4fn,-20*log10(abs(h(2);omega = 0: 200: 12000*pi;h = freqs(num,den,omega);

16、gain = 20*log10(abs(h);plot (omega/(2*pi),gain);xlabel(Frequency in Hz); ylabel(Gain in dB);Ap= 1.0000 As= 47.846711723311102187. 11088. 4102099. 6102187. 1)( ssssH0100020003000400050006000-60-50-40-30-20-100Frequency in HzGain in dB4.2.3椭圆滤波器（考尔滤波器） 特点：幅值响应在通带和阻带内都是等波纹的，对于给定的阶数和给定的波纹要求，椭圆滤波器能获得较其它滤

17、波器更窄的过渡带宽，就这点而言，椭圆滤波器是最优的。 其振幅平方函数为RN（，L）雅可比椭圆函数L表示波纹性质的参量),(11)()(2222LRjHANa N=5， 的特性曲线 可见，在归一化通带内（-11）， 在（0，1）间振荡，而超过L后 ， 在 间振荡。这一特点使滤波器同时在通带和阻带具有任意衰减量。 ),(25LR2L),(2LRN 下图为典型的椭圆滤波器振幅平方函数 椭圆滤波器的振幅平方函数 图中和A的定义 同切比雪夫滤波器当c、r、和A确定后，阶次N的确定方法为：)1()()1()(1/212121kKkKkKkKNAkkrc确定参量确定参数2/1222/1210)1 ()1 (

18、)(tktdtkK式中 为第一类完全椭圆积分 MATLAB设计椭圆低通滤波器椭圆低通滤波器N,Wc=ellipord(Wp,Ws,Ap,As,s)确定椭圆滤波器的阶数N。Wc=Wp。num,den=ellip(N,Ap,As,Wc,s)确定阶数为N，通带衰减为Ap dB，阻带衰减为As dB的椭圆滤波器的分子和分母多项式。Wc是椭圆滤波器的通带截频。 上面讨论了三种最常用的模拟低通滤波器的特性和设计方法，设计时按照指标要求，合理选用。 一般，相同指标下，椭圆滤波器阶次最低，切比雪夫次之，巴特沃兹最高，参数的灵敏度则恰恰相反。 4.2.4 模拟高通、带通及带阻滤波器的设计 模拟高通带通带阻滤波器

19、可以通过简单的频率变换来实现。设计过程为： 将模拟滤波器的参数指标通过频率变换转化为原型模拟低通滤波器的参数指标。 设计满足指标要求的原型模拟低通滤波器。 通过频率变换将原型模拟低通滤波器的系统函数变换为其他类型（高通、带通和带阻）的模拟滤波器MATLAB实现模拟滤波器的频率变换 低通到高通变换B，A=lp2hp(num, den, w0)给定归一化模拟低通滤波器分子和分母多项式的系数num和den，通过频率变换求得高通模拟滤波器，w0为高通滤波器的通带边界角频率。 低通到带通变换B，A=lp2bp(num, den, w0, Bw)给定归一化模拟低通滤波器分子和分母多项式的系数num和den

20、，通过频率变换求得带通模拟滤波器，w0为带通滤波器的中心角频率，w为带通滤波器的通带带宽（rad/s）。 低通到带通变换B，A=lp2bs(num, den, w0, Bw)给定归一化模拟低通滤波器分子和分母多项式的系数num和den，通过频率变换求得带阻模拟滤波器，w0为阻带的中心角频率，w为阻带带宽（rad/s）。 4.3 根据模拟滤波器设计IIR滤波器利用模拟滤波器设计数字滤波器，就是从已知的模拟滤波器系统函数Ha(s)设计数字滤波器系统函数H(z)，这归根到底是一个由S平面到Z平面的变换，这种映射变换应遵循两个基本原则：1）H（z）的频响要能模仿Ha(s)的频响，即S平面的虚轴应映射到Z平面的单位圆 上。2）H