实验四 数字低通巴特沃斯滤波器的设计.doc

实验四 数字低通巴特沃斯滤波器的设计实验目的 1.掌握IIR数字滤波器的设计方法实验原理 IIR数字滤波器的设计，主要采用间接法，即：首先设计出低通模拟滤波器H(S)；进行频率变换，将其转换为高通、带通、带阻滤波器；再用脉冲响应不变法或双线性变换法从模拟滤波器转换为数字滤波器。对单极点的N阶H(S)用部分展开式：冲激响应不变法取H(S)的单阶极点Sk的指数函数作为H（Z）的极点双线性变换法是用代换H（S）中的S得到H(Z),双线性变换法可完全消除频率混叠失真但存在非线性频率失真，而冲激响应不变法存在混叠失真。在不同的设计阶段MATLAB的信号处理工具箱都给出了相应的滤波器设计函数，这些函数代表了不同类型的逼近函数的滤波器，常用的有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器类型。设计流程图下图所示为IIR滤波器设计过程及所用到的MATLAB工具箱中的函数。模拟原型低通滤波器设计（buttap，cheb1ap，cheb2ap，bess1lap，ellipap函数）求最小阶数N（buttord，cheb1ord，cheb2ord，ellipord函数）频率转换（低通转低通lp2lp，低通转高通lp2hp，低通转带通lp2bp， 低通转带阻函数lp2bs）模拟滤波器转数字滤波器（bilinear，impinvar函数）合成一步的设计函数（butter，cheb1，cheb2，ellip，besself）滤波器系数A,B实验内容1.设计一个低通巴特沃斯模拟滤波器，绘制滤波器的的幅频响应及零极点图。指标如下：通带截止频率：WP1000HZ, 通带最大衰减：RP=3dB阻带截止频率：Ws2000HZ, 阻带最小衰减：Rs=40 dB参考程序butter1.m2. 用冲激响应不变法和双线性变换法将一模拟低通滤波器转换为数字低通滤波器并图释H(S)和H(Z)，采样频率Fs1000Hz参考程序butter2.m实验报告1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明频率特性应包括幅频、相频特性2.理论计算模拟低通原型滤波器的阶数和极点分布，说明实验中所用的冲激响应不变法和双线性变换法的原理及其编程根据3.总结实验体会及实验中存在的问题思考：1.双线性变换法和冲激响应不变法比较有何优点？IIR设计函数：熟悉MATLAB巴特沃斯滤波器为例介绍设计函数的功能和用法，其它设计函数的用法类似。（1）求最小阶数N的函数buttord调用格式1：N, Wn = BUTTORD(Wp, Ws, Rp, Rs, s) 功能：求出巴特沃斯模拟滤波器的最小阶数N和频率参数Wn 说明：Wp, Ws,为通带截止频率和阻带截止频率如式所示；，Rp, Rs为通带最大衰减 Apass(dB)和阻带最小衰减Astop(dB)如式s对应模拟滤波器。（2）模拟低通滤波器设计函数buttap调用格式：Z,P,K = BUTTAP(N) 功能：按给定的阶数N设计出巴特沃斯模拟低通滤波器 说明：Z,P,K为返回的N阶模拟滤波器的零点、极点和增益系数。如： 式3.5.4其中k为增益常数，zj和pi为系统函数H(Z)的N个零点和N个极点。 (3)模拟数字变换函数双线性变换函数bilinear或脉冲响应不变法函数impinvar调用格式：NUMd,DENd = BILINEAR(NUM,DEN,Fs) 功能： 把模拟滤波器系数为NUM,DEN变为近似等价的数字滤波器系数为NUMd,DENd调用格式：BZ,AZ = IMPINVAR(B,A,Fs) 功能：设计出数字滤波器其单位抽样响应h(n)为模拟滤波器的单位冲激响应h(t)的采样值，采样频率为Fs 说明：B,A为模拟滤波器的系统函数H(S)的分子和分母多项式的系数，如下所示：B=b1,b2,bm,bm+1,A=a1,a2,an,an+1(4) 合为一步的数字滤波器设计函数butter调用格式1：B,A = BUTTER(N,Wn) 功能：设计N阶低通数字滤波器，系数矢量B，A按Z的降幂排列，当Wn为二元矢量时，Wn = W1 W2，返回一个2N阶带通滤波器，通带范围：W1 W W2，如式3.5.6调用格式2：B,A = BUTTER(N, Wn, high) 设计高通数字滤波器.调用格式3：B,A = BUTTER(N, Wn, stop) 设计带阻滤波器 Wn = W1,W2.调用格式4：Z,P,K = BUTTER(N,Wn) 返回低通数字滤波器的N阶零、极点矢量. FREQS 调用格式： H,W = FREQS(B,A,M) 功能：此函数可以求出系统频率响应的数值解，并列出此系统的幅频及相频响应曲线。参考程序：熟悉MATLAB1.butter1.m巴特沃兹滤波器的幅频响应图subplot(1,2,1)； 分两个窗口，幅频图在第一个窗口wp=1000;ws=2000;rp=3;rs=40; 设置指标N,wn=buttord(wp,ws,rp,rs,s) 计算巴特沃斯低通滤波器的阶数和3dB截止频率B,A=butter(N,wn,s); 代入N和Wn设计巴特沃斯模拟低通滤波器Z,P,K=buttap(N); 计算滤波器的零、极点 h,w=freqs(B,A,1024); 计算1024点模拟滤波器频率响应h，和对应的频率点w画频率响应幅度图plot(w,20*log10(abs(h)/abs(h(1)grid;xlabel(频率Hz);ylabel(幅度（dB）);给x轴和y轴加标注title(巴特沃斯幅频响应) 给图形加标题axis(0,3000, -40,3);line(0,2000,-3,-3);line(1000,1000,-40,3);%绘制巴特沃斯滤波器的极点图subplot(1,2,2) 在第二个窗口画极点图p=P;q=Z;x=max(abs(p,q);x=x+0.1;y=x;axis(-x,x,-y,y);axis(square)plot(-x,x,0,0);hold onplot(0,0,-y,y);hold onplot(real(p),imag(p),x)程序运行结果：N = 7wn = 1.0359e+003Z = P = -0.2225 + 0.9749i -0.2225 - 0.9749i -0.6235 + 0.7818i -0.6235 - 0.7818i -0.9010 + 0.4339i -0.9010 - 0.4339i -1.0000 K = 1.0000图例12.butter2.m设模拟低通滤波器的系统函数为：b=1;a=1,1000;w=0:1000*2*pi;模拟频率为2f，其中f取01000Hzh,w=freqs(b,a,w);计算模拟滤波器的频率响应subplot(2,2,1)plot(w/2/pi,abs(h)/abs(h(1);grid;画模拟滤波器幅频特性title(模拟频率响应);xlabel(f(Hz);ylabel(幅度); Fs=1000;bz,az=impinvar(b,a,Fs); 冲激响应不变法设计数字滤波器 bzl,azl=bilinear(b,a,Fs); 双线性变换法设计数字滤波器wz=0:pi/512:pi;hz1=freqz(bz,az,wz);hz2=freqz(bzl,azl,wz);subplot(2,2,2);plot(wz/pi,abs(hz1)/hz1(1);grid; 画出冲激响应不变法滤波器的幅频图，axis(0,1,0,1) 数字频率wz归一化为01title(冲激响应不变法数字频率响应) subplot(2,2,3);plot(wz/pi,abs(hz2)/hz2(1);grid; 画出双线性变换法滤波器的幅频图，axis(0,1,0,1) 数字频率wz归一化为01title(双线性变换法数字频率响应); 图例2 实验五 窗函数设计FIR滤波器(选做)实验目的 1.掌握利用窗函数法设计线性相位FIR数字滤波器的方法实验原理 用窗函数法设计FIR滤波器时，先根据给定的Wp和N通过傅立叶反变换求得相应的理想滤波器单位脉冲响应，由于是一个无限长序列，所以第二步要选择一个合适的窗函数w(n)来截取到合适的长度（阶数）以保证实现要求的阻带衰减，最后得到FIR滤波器单位冲激响应h(n)w(n)。即：其中实验内容1. 利用窗生成函数boxcar，hanning，bartlett设计FIR低通数字滤波器，指标如下：采样频率Fs1000Hz，截止频率ws200Hz，滤波器阶数N65，画出滤波器的幅度频率特性和窗的形状，参考程序chfir1.m2. 用基于窗函数的滤波器设计函数fir1，重新设计上述滤波器实验报告1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明2.写出理论计算低通滤波器的h(n)的方法3.总结实验体会及实验中存在的问题思考：1.比较矩形窗、汉宁窗、巴特里特窗的滤波器特性FIR设计函数：熟悉MATLAB（1）矩形窗生成函数hanning调用格式：W =HANNING(N) 功能：生成N点的汉宁窗 说明：N为矩形窗时域的长度，W为行矢量，其元素为对应矩形窗时域每一采样点的幅度值。（2）基于窗函数的滤波器设计函数fir1调用格式1：B = FIR1(N,Wn,ftype，window) 功能：以给定的窗函数（window例如：bartlett(N+1)），和滤波器类型（high、low、stop、band）设计出预期理想频率特性的滤波器。 说明：N为滤波器阶数，Wn为通带截止频率如式3.5.3的Wp，B为滤波器系数，即单位抽样响应h(n)b(n+1) n0,1,2,N；如果ftype和window缺省，该函数得到截止频率为Wn且满足线性相位条件的N阶FIR低通滤波器，window默认为hamming窗。FREQZ 功能：求H(Z)的系统频率响应的数值解，并画出此系统的幅频及相频响应曲线。 调用格式： H = FREQZ(B,A,W) 说明：此调用格式计算由向量w(范围在0 之间)指定的频率点上系统的频率响应，向量B和A分别为系统函数H(Z)的分子和分母多项式的系数。例如Z域的系统函数的形式为：则：B=b0,b1,bM-1,bM,A=a0,a1,aN-1,aN，要求A的首项不为0参考程序：熟悉MATLABchfir.mwp=200/1000*2*pi;将频率转换为数字角频率，单位为n=1:65;hd=sin(wp*(n-32)./(pi*(n-32);算出理想滤波器的冲激响应hd(32)=wp/pi; 避免除法溢出w1=boxcar(65); 生成矩形窗subplot(2,3,1) h1=hd.*rot90(w1); 给理想低通数字滤波器加矩形窗mag1,r1=freqz(h1); 求加窗后的频率特性plot(r1,20*log(abs(mag1);画出加矩形窗后的幅频特性图title(矩形窗) w2=triang(65); subplot(2,3,2)h2=hd.*rot90(w2); 给理想低通数字滤波器加巴特里特窗 mag2,r2=freqz(h2);plot(r2,20*log(abs(mag2); 画出加巴特里特窗后的幅频特性图