用MATLAB设计FIR数字滤波器

实验八用MATLAB设计FIR数字滤波器（二）一、实验目的：1、加深对窗函数法设计FIR数字滤波器的基本原理的理解。2、学习用MATLAB语言的窗函数法编写设计FIR数字滤波器的程序。3、了解MATLAB语言有关窗函数法设计FIR数字滤波器的常用函数用法。二、实验原理：1、用窗函数法设计FIR数字滤波器FIR数字滤波器的系统函数为N-1-nn=0H(z)=h(n)z这个公式也可以看成是离散LSI系统的系统函数M-m-1-2-mmm=0012mN-1-2-k-k12kkk=1bzb+bz+bz++bzY(z)b(z)H(z)====X(z)a(z)1+az+az++az1+az分母a0为1，其余ak全都为0时的一个特例。由于极点全部集中在零点，稳定和线性相位特性是FIR滤波器的突出优点，因此在实际中广泛使用。FIR滤波器的设计任务是选择有限长度的h(n)，使传输函数H(ejω)满足技术要求。主要设计方法有窗函数法、频率采样法和切比雪夫等波纹逼近法等。本实验主要介绍窗函数法。用窗函数法设计FIR数字滤波器的基本步骤如下：（1）根据过渡带和阻带衰减指标选择窗函数的类型，估算滤波器的阶数N。（2）由数字滤波器的理想频率响应H(ejω)求出其单位脉冲响应hd(n)。可用自定义函数ideal_lp实现理想数字低通滤波器单位脉冲响应的求解。程序清单如下：functionhd=ideal_lp(wc,N)%点0到N-1之间的理想脉冲响应%wc=截止频率（弧度）%N=理想滤波器的长度tao=(N-1)/2;n=[0:(N-1)];m=n-tao+eps;%加一个小数以避免0作除数hd=sin(wc*m)./(pi*m);其它选频滤波器可以由低通频响特性合成。如一个通带在ωc1～ωc2之间的带通滤波器在给定N值的条件下，可以用下列程序实现：Hd=ideal_lp(wc2,N)-ideal_lp(wc1,N)（3）计算数字滤波器的单位冲激响应h(n)=w(n)hd(n)。（4）检查设计的滤波器是否满足技术指标。如果设计的滤波器不满足技术指标，则需要重新选择或调整窗函数的类型，估算滤波器的阶数N。再重复前面的四个步骤，直到满足指标。常用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、切比雪夫窗、布莱克曼窗、凯塞窗等，MATLAB均有相应的函数可以调用。另外，MATLAB信号处理工具箱还提供了fir1函数，可以用于窗函数法设计FIR滤波器。由于第一类线性相位滤波器（类型Ⅰ）能进行低通、高通、带通、带阻滤波器的设计，因此，本实验所有滤波器均采用第一类线性相位滤波器。三、实验内容：1、阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。2、选择合适的窗函数设计FIR数字低通滤波器，要求：ωp=0.2π，Rp=0.05dB；ωs=0.3π，As=40dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解：程序清单如下：wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;N0=ceil(6.6*pi/deltaw);N=N0+mod(N0+1,2)%为实现FIR类型1偶对称滤波器，应确保N为奇数windows=(hamming(N))';wc=(ws+wp)/2;%截止频率hd=ideal_lp(wc,N);b=hd.*windows;[db,mag,pha,grd,w]=freqz_m(b,1);%dbn=0:N-1;dw=2*pi/1000;Rp=-(min(db(1:wp/dw+1)))%检验通带波动As=-round(max(db(ws/dw+1:501)))%检验最小阻带衰减subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');xlabel('n');ylabel('h(n)');subplot(2,2,2);stem(n,windows);axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)');subplot(2,2,3);plot(w/pi,db);axis([0,1,-80,10]);title('幅度频率响应');xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})');set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]);set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);gridsubplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]);title('相位频率响应');xlabel('频率（单位：\pi）');ylabel('\phi(\omega)');set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]);set(gca,'YTickMode','manual','YTick',[-3.1416,0,3.1416,4]);grid程序运行结果：N=67Rp=0.0394As=52波形如下：0102030405060-0.0500.050.10.150.20.25实际脉冲响应nh(n)020406000.20.40.60.81窗函数特性nwd(n)00.20.31-50-20-30幅度频率响应频率（单位：）H(ej)00.20.31-3.141603.14164相位频率响应频率（单位：）()3、用凯塞窗设计一个FIR数字高通滤波器，要求：ωp=0.3π，Rp=0.1dB；ωs=0.2π，As=50dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解：程序清单如下：wp=0.3*pi;ws=0.2*pi;deltaw=wp-ws;N0=ceil(6.1*pi/deltaw);N=N0+mod(N0+1,2)%为实现FIR类型1偶对称滤波器，应确保N为奇数windows=(kaiser(N))';wc=(ws+wp)/2;hd=ideal_lp(pi,N)-ideal_lp(wc,N);b=hd.*windows;[db,mag,pha,grd,w]=freqz_m(b,1);n=0:N-1;dw=2*pi/1000;Rp=-(min(db(wp/dw+1:501)))%检验通带波动As=-round(max(db(1:ws/dw+1)))%检验最小阻带衰减subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');xlabel('n');ylabel('h(n)');subplot(2,2,2);stem(n,windows);axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)');subplot(2,2,3);plot(w/pi,db);axis([0,1,-40,2]);title('幅度频率响应');xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})');set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]);set(gca,'YTickMode','manual','YTick',[-20,-3,0]);gridsubplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]);title('相位频率响应');xlabel('频率（单位：\pi）');ylabel('\phi(\omega)');set(gca,'XTickMode','manual','XTick',[0,ws/pi,wp/pi,1]);set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid运行结果：N=61Rp=1.1756As=29波形如下：0102030405060-0.200.20.40.60.8实际脉冲响应nh(n)010203040506000.20.40.60.81窗函数特性nwd(n)00.20.31-20-30幅度频率响应频率（单位：）H(ej)00.20.31-3.141603.1416相位频率响应频率（单位：）()4、选择合适的窗函数设计一个FIR数字带通滤波器，要求：fp1=3.5kHz，fp2=6.5kHz，Rp=0.05dB；fs1=2.5kHz，fs2=7.5kHz，As=60dB。滤波器采样频率Fs=20kHz。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解：程序清单如下：fp1=3500;fp2=6500;fs1=2500;fs2=7500;Fs=20000;ws1=fs1/(Fs/2)*pi;ws2=fs2/(Fs/2)*pi;wp1=fp1/(Fs/2)*pi;wp2=fp2/(Fs/2)*pi;deltaw=wp1-ws1;N0=ceil(11*pi/deltaw);N=N0+mod(N0+1,2)%为实现FIR类型1偶对称滤波器，应确保N为奇数windows=blackman(N);wc1=(ws1+wp1)/2/pi;wc2=(ws2+wp2)/2/pi;b=fir1(N-1,[wc1,wc2],windows);[db,mag,pha,grd,w]=freqz_m(b,1);n=0:N-1;dw=2*pi/1000;Rp=-(min(db(wp1/dw+1:wp2/dw+1)))%检验通带波动ws0=[1:ws1/dw+1,ws2/dw+1:501];%建立阻带频率样点数组As=-round(max(db(ws0)))%检验最小阻带衰减subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');xlabel('n');ylabel('h(n)');subplot(2,2,2);stem(n,windows);axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)');subplot(2,2,3);plot(w/pi,db);axis([0,1,-150,10]);title('幅度频率响应');xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})');set(gca,'XTickMode','manual','XTick',[0,fs1,fp1,fp2,fs2,8000]);set(gca,'YTickMode','manual','YTick',[-150,-40,-3,0]);gridsubplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]);title('相位频率响应');xlabel('频率（单位：\pi）');ylabel('\phi(\omega)');set(gca,'XTickMode','manual','XTick',[0,fs1,fp1,fp2,fs2,500]);set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid运行结果：N=111Rp=0.0034As=74波形图如下：020406080100-0.200.20.4实际脉冲响应nh(n)02040608010000.20.40.60.81窗函数特性nwd(n)0-150-40-30幅度频率响应频率（单位：）H(ej)00.51-4-2024相位频率响应频率（单位：）()5、选择合适的窗函数设