实验四IIR数字滤波器的设计

实验四IIR数字滤波器的设计及网络结构一、实验目的1．了解IIR数字滤波器的网络结构。2.掌握模拟滤波器、IIR数字滤波器的设计原理和步骤。3．学习编写数字滤波器的设计程序的方法。二、实验内容数字滤波器：是数字信号处理技术的重要内容。它的主要功能是对数字信号进行处理，保留数字信号中的有用成分，去除信号中的无用成分。1．数字滤波器的分类滤波器的种类很多，分类方法也不同。（1）按处理的信号划分：模拟滤波器、数字滤波器（2）按频域特性划分；低通、高通、带通、带阻。（3）按时域特性划分：FIR、IIR2．IIR数字滤波器的传递函数及特点数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行处理，改变输入信号的波形或频谱，达到保留信号中有用成分去除无用成分的目的。如果加上A/D、D/A转换，则可以用于处理模拟信号。设IIR滤波器的输入序列为x(n)，则IIR滤波器的输入序列x(n)与输出序列y(n)之间的关系可以用下面的方程式表示：01()()()MNijijynbxniaynj（5-1）其中，ja和ib是滤波器的系数，其中ja中至少有一个非零。与之相对应的差分方程为：10111....()()()1....MMNNbbzbzYzHZXzazaz（5-2）由传递函数可以发现无限长单位冲激响应滤波器有如下特点：（1）单位冲激响应h(n)是无限长的。（2）系统传递函数H(z)在有限z平面上有极点存在。（3）结构上存在着输出到输入的反馈，也就是结构上是递归型的。3．IIR滤波器的结构IIR滤波器包括直接型、级联型和并联型三种结构：①直接型：优点是简单、直观。但由于系数bm、ak与零、极点对应关系不明显，一个bm或ak的改变会影响H(z)所有零点或极点的分布，所以一方面，bm、ak对滤波器性能的控制关系不直接，调整困难；另一方面，零、极点分布对系数变化的灵敏度高，对有限字长效应敏感，易引起不稳定现象和较大误差。Matlab实现：filter()函数实现IIR数字滤波器直接形式。格式为：y=filter(b,a,x)b，a为差分方程输入、输出系数向量（或系统函数的分子、分母多项式，降幂），x为输入序列，y为输出序列。其中，传递函数（tf）形式NNMMzazazbzbbzH111101)(若则a=[1a1a2…aN]b=[b0b1b2…bM]②级联型：基于因式分解，将系统函数H(z)分解为因子乘积的形式。(5-3）级联型结构：Matlab实现：tf2zp()函数用于求系统函数的零、极点和增益常数，zp2sos()函数则根据tf2zp()函数结果求出各基本节系数。格式为：[z,p,K]=tf2zp(b,a)；sos=zp2sos(z,p,K)；b，a为差分方程输入、输出系数向量（系统函数的分子、分母多项式，降幂）。其中，零极点增益形式（zp）：101101)1()1()(NkkMiizzzzkzH若则零点向量z=[z1z2…zM-1];极点向量p=[z1,z2,…,zN-1]001201212111211()()11MmNNmmkkkNkkkkkkkbzzzHzKKHzzzazk为系统增益。二阶分式形式（sos）为：把H(z)划成二阶因式NkNkkkkkkkkzzzzzHzH112211022110)()(则其二阶因式为：NNNNNbbbbbaabbbaabbbsos2121022122212022111211101111③并联型：基于部分分式展开，将系统函数H(z)分解为部分分式和的形式。（5-4）并联型结构：Matlab实现：residue()函数可以实现并联型结构，有两种格式：[K,r,p]=residue(b,a);[b,a]=residue(b,a);其中，部分分式形式：)(1)1()1()(1)()1(1)1()(NMnzNMkkznpnrzprzH若则极点向量p=[p(1)p(2)…p(n)]其对应系数向量r=[r(1)r(2)…r(n)]余数多项式系数向量k=[k(1)k(2)…k(M-N+1)]【实例5-1】已知三阶IIR数字滤波器的系统函数001001001211121()()11MmNNmmkkkNkkkkkkkbzzHzKKHzzzaz32121613161132353)(zzzzzzH求：①直接形式的单位采样响应h(n)；②级联型结构的各基本节系数；③并联型结构的部分分式系数。解：MATLAB源程序为①b=[3,5/3,2/3];a=[1,1/6,1/3,-1/6];x=[1,zeros(1,50)];y=filter(b,a,x);n=0:50;plot(n,y);②b=[3,5/3,2/3,0];a=[1,1/6,1/3,-1/6];[z,p,K]=tf2zp(b,a);sos=zp2sos(z,p,K);③b=[3,5/3,2/3];a=[1,1/6,1/3,-1/6];[K,r,p]=residue(b,a);KK1=[K(1),K(2)];zz1=[z(1),z(2)];[b2,a2]=residue(KK1,zz1,0);5．IIR数字滤波器的具体设计（1）巴氏模拟原型滤波器的设计巴氏模拟低通滤波器幅度平方函数为NCjH211)(MATLAB工具箱函数buttap,buttord和butter是巴氏滤波器设计函数。1）[Z,P,K]=buttap(N)该格式用于计算N阶巴氏归一化（3dB截止频率Ωc=1）模拟低通原型滤波器系统函数的零、极点和增益。得到的系统函数为：)())(()())(()(2121NNappppppzpzpzpKpG如果要从计算得到的零、极点得到系统函数的分子和分母向量B和A，可以调用结构转换函数[B,A]=zp2tf(Z,P,K)。2）[N,wc]=buttord(wp,ws,Rp,Rs,’s’)该格式用于计算巴氏滤波器的阶数N和3dB截止频率wc。其中：N——滤波器的阶数。Wc——3dB截止频率的归一化值。wp、ws——通带、阻带边界频率的归一化值。要求：1ws0,1wp0，1表示数字频率π。Rp、Rs——通带最大、阻带最小衰减。s——可选项，直接设计模拟滤波器，此时wp和ws均为实际模拟角频率。3）[B,A]=butter(N,wc,’ftype’，’s’)计算N阶巴氏数字滤波器系统函数分子和分母多项式的系数向量B和A。【实例5-2】已知通带截止频率kHzfp5，通带最大衰减dBap5，阻带截止频率kHzfs12，阻带最小衰减dBas30，设计巴特沃斯型模拟低通滤波器。解：MATLAB源程序为Wp=2*pi*5000;Ws=2*pi*12000;Rp=5;As=30;[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');[B,A]=butter(N,Wc,'s');freqs(B,A);【实例5-3】设计阶数为3，5，10，15的巴氏模拟原型滤波器。并画出幅频响应曲线。解：MATLAB源程序为fori=1:4;switchicase1N=3;case2N=5;case3N=10;case4;N=15;end;[z,p,k]=buttap(N);[b,a]=zp2tf(z,p,k);[h,w]=freqs(b,a,n);Ah=abs(h);subplot(2,2,i),plot(w,Ah);axis([0201]);xlabel('w/wc');ylabel('|H(jw)|.^2');title(['filerN=',num2str(N)]);grid;end;【实例5-4】设通带、阻带截止频率fp=0.5kHz、fs=1.2kHz，通带、阻带最大衰减Rp=1dB,Rs=30dB，要求设计巴氏低通滤波器。解：MATLAB源程序为:OmegaP=2*pi*500;OmegaS=2*pi*1200;Rp=1;Rs=30;[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,Rs,'s');%确定阶数N[z0,p0,k0]=buttap(N);%确定传递函数z=z0*OmegaC;%去归一化k=k0*OmegaC^N;p=p0*OmegaC;bs=k*real(poly(z));as=real(poly(p));freqs(bs,as);（2）IIR数字滤波器的设计模拟滤波器Ha(s)转换成数字滤波器H(z)应满足要求：(1)因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的；(2)数字滤波器的频率响应模仿模拟滤波器的频响，s平面的虚轴映射z平面的单位圆，相应的频率之间成线性关系。脉冲响应不变法和双线性变换法都满足如上要求。①脉冲响应不变法用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t),让h(n)正好等于ha(t)的采样值，即h(n)=ha(nT)，其中T为采样间隔。②双线性变换法s平面与z平面之间满足以下映射关系：1111zzss平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。双线性变换时一种非线性变换)2/(tg，这种非线性引起的幅频特性畸变可通过预畸而得到校正。以低通数字滤波器为例，将设计步骤归纳如下：(1)确定数字滤波器的性能指标：通带临界频率fp、阻带临界频率fs；通带内的最大衰减Rp；阻带内的最小衰减Rs；(2)确定相应的数字角频率，ωp=2πfp；ωs=2πfs；(3)计算经过预畸的相应模拟低通原型的频率，)2/(tg；(4)根据Ωp和Ωs计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数Ha(s)；(5)用上面的双线性变换公式代入Ha(s)，求出所设计的传递函数H(z)；(6)分析滤波器特性，检查其指标是否满足要求。1)冲激响应不变法MATLAB实现：采用impinvar函数[bz,az]=impinvar(b,a,fs,tol)bz,az――转换的分子分母向量fs――采样频率，默认值为1Hztol――误差容限，表示转换后的离散系统函数是否有重复的极点。2）双线性变换法MATLAB实现采用bilinear函数：[bz,az]=bilinear(bs,as,fs)――把模拟滤波器的传递函数转换为数字的其中：bs,as——模拟滤波器的分子、分母向量bz,az——数字滤波器的分子、分母向量【实例5-5】基于Butterworth模拟滤波器原型，使用双线性转换设计数字滤波器，其中参数指标为：通带截止频率：0.2p通带波动值：dBRp1阻带截止频率：0.3s阻带波动值：dBAs15解：首先确定滤波器的阶数N，同时根据sp和确定c＝0.5。接着使用bilinear进行双线性转换，最后绘制在频域上的各种图像，其源代码如下：％数字滤波器指标wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;%将数字滤波器指标反转变化为模拟滤波器的参数T=1;fs=1/T;omegap=(2/T)*tan(wp/T);omegas=(2/T)*tan(ws/T);%butterworth原型模拟滤波器的设计[N,Wc]=buttord(omegap,omegas,Rp,As,'s');[B,A]=butter(N,Wc,'s');%双线性变换[b,a]=bilinear(B,A,fs);%频域图像的绘制freqz(b,a);程序运行后，产生4阶的butterworth数字滤波器，频率响应如图5-5所示的波形。图5-5Butterworth