MATLAB窗函数法实现FIR的高通-带通和低通滤波器的程序

1、MATLAB课程设计报告学 院：地球物理与石油资源学院班级：测井（基）11井1姓名： 大牛啊啊啊学 号：班内编号：指导教师：陈义群完成日期： 2013年6月3日题目 FIR滤波器的窗函数设计法及性能比较1. FIR滤波器简介数字波波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来 达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR） 滤波器和有限冲激响应（FIR）滤波器。与IIR滤波器相比，FIR滤波器的主要特点为： a.线性相位；b.非递归运算。2. FIR滤波器的设计FIR滤波器的设计方法主要有三种：a.窗函数设计法；b.频率抽样发：c

2、.最小平法抽样法: 这里我主要讨论在MATLAB环境下通过调用信号分析与处理工具箱的几类窗函数来设计滤 波器并分析与比较其性能。窗函数法设计FIR滤波器的一般步骤如下： a.根据实际问题确定要设计的滤波器类型：b.根据给定的技术指标，确定期望滤波器的理想频率特性；c.求期望滤波器的单位脉冲响应：d.求数字滤波器的单位脉冲响应；e.应用。常用的窗函数有(1)Hanning window(9.2.2)(9.2.3)w(n = 0.5-0.5cos(-，l )RN(n)N 1(2)Hamnung windoww(n) = 0.54-0.46cos(-)7?v(/?)N 13. ) Balck man

3、 windoww()= 0,42-0.5cos() + 0.08cos(i)Rv()(9.2.4)N 1N-1(4)Kaiser window乂'()=/。（(9.2.5)式中/。是零阶8css”函数，可定义为Z0(x) = l + £ /H-l(x/2)阳当x = 0时与矩形窗一致；当, = 5.4414H寸与海明窗结果相同;当x = 8.885时与布莱克曼窗结果札同。4. 常用窗函数的参数常用窗函数的几个参照数据窗类型旁瓣峰值（dB）主辩宽度最小阻带衰减（dB）矩形窗-134%AV-21三角窗-258%/丫-25汉宁窗-31脉/N-44哈明窗-41WN-53布莱克曼窗-5

4、712，N-745. FIR滤波器的MATLAB实现方式在MATLAB信号分析与处理工具箱中提供了大量FIR窗函数的设计函数，本次用到主要有以下几种： harming (N) hamming(N) blackman(N) kaiser(n+lz beta)h a rm i n g窗函数的调ffl hamming窗函数的调用 blac kma n窗函数的调ffl kaiser窗函数的调用kaiserordf reqz filter6.实验具体步骤计算kaiser窗函数的相关参数 求取频率响应对信号进行滤波的函数本次实验分别通过调用hanning , hamming , Blackman, kai

5、ser窗函数，给以相同的技术参数，来设计低通，带通，高通滤波器，用上述窗函数的选择标准来比较各种窗 函数的优劣，并给以一个简谐波进行滤波处理，比较滤波前后的效果。达到综合比较的效果。二、源代码L利用 banning hamming blackmail kaiser 窗,设计一个低通 FIR function lowpassfilter clc;clear all;Fs=100;%采样频率fp=2O;%通带截止频率fs=3O;%阻带起始频率wp=2*pi*fp/Fs;%将模拟通带截止频率转换为数字滤波器频率 ws=2*pi*fs/Fs;%将模拟阻带起始频率转换为数字波波器频率 wn=（wp+ws

6、）/2/pi; %标准化的截止频率响应Bt=ws-wp;N0=ceil(6.2*pi/Bt);% 滤波器长度N=NO+mod(NO+ L2);windowl=hanning(N);%使用 banning 窗函数window2=hamming(N);%使用 hamming 窗函数window3=blackman(N);%使用 blackman 窗函数 n,Wn,beta,ftype=kaiserord(20 25,1 OJJO.O1 O.O1JJOO);window4=kaiser(n+l,beta);%使用 kaiser 窗函数%设计加窗函数fklb 1 =fir 1 (N-1 ,wn. wi

7、ndow 1);b2=fir 1 (N-1, wn, window2);b3=fir 1 (N-1 avil window3);b4=fir 1 (nAVn/pi,window4 ,'noscale');%求取频率响应|HLWl=freqz(bl,L512,2);H2.W2=freqz(b2,1,512,2);H3,W3=freqz(b3,1,512,2);|H4.W4=freqz(b4,L512,2);figure(l);subplot(2.2 J ).plot( W L20*log 10(abs(H 1);% 绘制频率响应图形 axis(0Jr100J00);title。低

8、通banning窗的频率响应图形,);xlabelf频率(Hz)');ylabdC 幅值,);subplot(252,2),plot(W2,20*log 10(abs(H2);% 绘制频率响应图形 axis(OJrlOOJOOJ);titled氐通hamming窗的频率响应图形')；xlabel('频率(Hz);ylabdC 幅值);subplot(223).plot(W3,20*logl0(abs(H3);% 绘制频率响应图形 axis(0Jr 100,100);title。低通blackman窗的频率响应图形')；xlabelf频率(Hz)');yl

9、abel(,幅值);subplot(22,4).plot(W4,20*log 10(abs(H4);% 绘制频率响应图形 axis(0Jr100J00J);title。低通kaiser窗的频率响应图形')；xlabel(1频率(Hz);ylabelC 幅值);T=l/Fs;L=100:%信号长度t=(O:L-l)*T;%定义时间范围和步长y=sin(2*pi*5*t)+5*sin(2*pi*15*l)+8*sin(2*pi*40*t);% 滤波前的图形NFFT = 2Anextpow2(L); % Next power of 2 from length of yY = fft(y,NF

10、FT"L;%将时域信号变换到频域 f = Fs/2*linspace(0.1,NFH72+l);% 频域采样 figure(2);plot(f,2*abs(Y( 1 :NFFT/2+l );xIabel(,frequency/Hz,);ylabeI(,AnnipIitudc,);%滤波前频谱 title ('滤波前的频谱);%滤波后频谱%采用banning窗波波器 yyl=filter(b 1,1 ,y);% 调用滤波函数YYl=fft(yy 1 ,NFFT)/L;%进行傅里叶变换，下同。 f l=Fs/2*linspace(0.1 .NFFT/2+1);figured);s

11、ubplot(2.2,1 ),plot(f 1,2*abs(YY1 (1 :NFFT/2+1) ;xIabel(,frequency/Hz,);ylabel(,Amuplitude,); titleChanning窗的滤波效果')；%采用hammning窗滤波器 yy2=filter(b2J,y);YY2=fft(yy2,NFFT)/L:fl =Fs/2*linspace(0,1 .NFFT/2+1);subplot(2,22hplot(fl,2*abs(YY2(l:NFFT/2+l) ;xlabel(,frequency/Hz,);ylabeI(,Amuplitude,); titl

12、e(力amming窗的滤波效果)；%采用blackman窗滤波器yy3=filter(b3J,y);YY3=fft(yy3,NFFT)/L： f l=Fs/2*linspace(0.1 .NFFT/2+1);subplot(2,23h plot(fL2*abs(YY3(l:NFFT/2+l) ;xlabel(Trequency/Hz');ylabcl('Amuplitude'); titlefblackman窗的滤波效果');%采用kaiser窗滤波器 yy4=filter(b4J,y);YY4=fft(yy4,NFFT)/L;fl=Fs/2*linspace(

13、0J .NFFT/2+1);subplot(2,24),plot(f 1,2*abs(YY4( 1 :NFFT/2+1) ;xlabel(Trequency/Hz');ylabel('Amuplitude'); xlabel(,frequency/Hz,):ylabel(,Amuplitude,);titleCkaiser窗函数滤波效果')；%滤波前后的信号的时域对比figure(4);plot(y);xlabel(时间/s);ylabclC 振幅);title('滤波前振幅特性)；figure(5);subplot(22D,plot(yyl)；xlab

14、el(时间/s');ylabel(振幅”itleChanning 窗函数滤波振幅特性);subpk)t(222),plot(yy2);xlabel(时间/s);ylabel(，振幅');title('hamming 窗函数滤波振幅特性);subplot(2,2,3),plot(yy3);xlabel(时间/s);ylabel(振幅，):title(blackman 窗函数滤波振幅特性)； subplot(224),plot(yy4);xlabel(时间/s);ylabelf振幅);titlefkaiser 窗函数滤波振幅特性,);%滤波前后的信号的相位对比figure(

15、6);ploi(angle(Y);xlabel('时间/s);ylabelC 相位以titlef 滤波前的相位特性);figure(7);subplot(252,l),plot(angle(YYl);xlabel('时间/s');ylabcl('相位'):title('hanning 窗函数滤波相位特 性)；subplot(222),plot(angle(YY2);xlabel('时间/s');ylabel('相位');title('hamming 窗函数滤波相位特性)；subpk)t(2,2,3),plo

16、t(angle(YY3);xlabel(时间/sRylabel。相位XtitleCblackman 窗函数滤波相位特 性)；subplot(224),plot(angle(YY4);xlabelC时间/s);ylabel('相位);titleCkaiser 窗函数滤波相位特性2 .设-个 lianning hamming blackmail kaiser 窗函数 bandpass_FIR%设计一个 banning hamming blackman kaiser 窗函数 bandpass_FIRfunction bandpassfilterFs=100;%采样频率fpl = 15;%通带

17、下限截止频率fp2=2O;%通带上限截止频率fsl=10;fs2=25;wpl=2*pi*fpl/Fs;%将通带下限截止频率转换为数字滤波器频率wp2=2*pi*fp2/Fs;%将通带上限截止频率转换为数字滤波器频率wsl=2*pi*fsl/Fs;%将通带下限截止频率转换为数字滤波器频率ws2=2*pi*fs2/Fs;%将通带上限截止频率转换为数字滤波器频率Bt=wpl-wsl;N0=ceil(6.2*pi/Bt):N=N0+mod(N0+1,2)；wn=(wp 1 +ws 1 )/2/pi,(wp2+ws2)/2/pi;windowl=hanning(N);%使用 banning 窗函数wi

18、ndow2=haniming(N);%使用 hamming 窗函数window3=blackn】an(N);%使用 blackman 窗函数%设过渡带宽度为5Hzn,Wn,beta,ftype=kaiserord( 10 15 20 25.0 1 050.01 0.01 0.01JOO);%求阶数n 以及参数 betawindow4=kaiser(n+l,beta);%使用 kaiser 窗函数%设计加窗函数fi”b 1 =fir 1 (N-1 ,wn. window 1);b2=fir 1 (N-1 ,wn. window2);b3=fir 1 (N-1 ,wn, window3);b4=f

19、irl(n.Wn.window4；noscale,);%求取频率响应HLWl=freqz(bLL512,2);H2.W2=freqz(b2,1,512,2);H3,W3=freqz(b34,512,2);H4.W4=freqz(b4,1,512,2);figured);subplot(2.2 J ).plot( W L20*log 10(abs(H 1);% 绘制频率响应图形axis(04<100,100);title。带通banning窗的频率响应图形')；xlabelC频率(Hz)');ylabdC 幅值");subplot(252,2),plot(W2,2

20、0*log 10(abs(H2);% 绘制频率响应图形axis(OJr 100,100);title。带通hamming窗的频率响应图形')；xlabel('频率(Hz);ylabelC 幅值);subpk)t(223).plot(W3,20*logl0(abs(H3);% 绘制频率响应图形axis(0Jr 100,100);title。带通blackman窗的频率响应图形xlabelC频率(Hz)');ylabelf 幅值);subplot(2,2,4),plot(W4,20*log 10(abs(H4);% 绘制频率响应图形axis(0Jr100,100|);tit

21、le。带通kaiser窗的频率响应图形')；xlabel(濒率(Hz)');ylabelC 幅值);T=l/Fs;L=100:%信号长度t=(0: L- 定义时间范围和步长y=sin(2*pi*5*t)+5*sin(2*pi*15 %)+8*sin(2*pi*40*t);% 滤波前的图形NFFT = 2Anextpow2(L); % Next power of 2 from length of yY = fft(y,NFFT)/L;%将时域信号变换到频域f = Fs/2*linspace(0J.NFFT/2+1);% 频域采样figure(2);plot(f,2*abs(Y(

22、1 :NFFT/2+1 );xlabel(,frequency/Hz,);ylabel(,Amuplitude,);%滤波前频谱 title ('滤波前的频谱')；%滤波后频谱%采用banning窗波波器yyl=filter(b 1,1 ,y);% 调用滤波函数YYl=fft(yyLNFFT)/L;%进行傅里叶变换，下同。fl =Fs/2*linspace(0.1 .NFFT/2+1);figured);subplot(2,24),plot(f 1,2*abs(YYl(l:NFFT/2+1) ;xIabel(,frequency/Hz,);yIabel(,Amuplitude,

23、);titleChanning窗的滤波效果')；%采用hammning窗滤波器yy2=filter(b2J,y);YY2=fft(yy2.NFFT)/L;fl =Fs/2*linspace(0,1 .NFFT/2+1);subplot(2,22),plot(f 1,2*abs(YY2( 1 :NFFT/2+1) ;xlabel(Trequency/Hz,):ylabel(,Amuplituder);title(力amming窗的滤波效果)；%采用blackman窗滤波器yy3=filter(b3,l,y);YY3=fft(yy3,NFFT)/L：fl =Fs/2*linspace(0.

24、1 .NFFT/2+1);subplot(223), plot(f 1,2*abs(YY3( 1:NFFT/2+1) ;xlabel(*frequency/Hz,);ylabcl(,Amuplitude,);titlefblackman窗的滤波效果);%采用kaiser窗滤波器yy4=filter(b4J,y);YY4=fft(yy4,NFFT)/L;fl=Fs/2*linspace(0J .NFFT/2+1);subplot(2.24),plot(fl,2*abs(YY4( 1 :NFFT/2+1) ;xlabel(,frequency/Hz,):ylabel(,Amuplitude,);x

25、label(,frequency/Hz,):ylabel(,Amuplitude,);title('kaiser窗函数滤波效果')；%波波前后的信号的时域对比figure(4);plot(y);xlabel(时间/sRylabclC 振幅RtitleC 滤波前振幅特性)；figure(5);subplot(2,2.1 ),plot(yyl);xlabel(时间/s);ylabelC振幅);titleChanning 窗函数滤波振幅特性);subplot(222),plot(yy2)：xlabel(时间/s);ylabel('振幅');title('han

26、iming 窗函数滤波振幅特性);subplot(223),plol(yy3)：xlabel(时间/sRylabelf振幅);titleCblackman 窗函数滤波振幅特性)； subplot(224),plot(yy4)：xlabelC时间/sRylabelC振幅);titleCkaiser 窗函数滤波振幅特性)； %滤波前后的信号的相位对比figure(6);plot(angle(Y);xlabel(时间/s);ylabelC 相位滤波前的相位特性)；figure(7);subpk)t(2,2),plot(angle(YYl);xlabelC时间/s);ylabelC相位);title(

27、'hanning 窗函数滤波相位特 性)；subplot(222),plot(ai】gle(YY2);xlabel('时间/s');ylabelC相位);title('hamming 窗函数滤波相位特 性')；subpk)t(2,2,3),plot(angle(YY3);xlabel('时间/s');ylabel('相位');title('blackn】an 窗函数滤波相位特 性)；subplot(2,2,4),pk)t(angle(YY4);xlabelC时间/s);ylabel(湘位"itleCkai

28、ser 窗函数滤波相位特性 ');3 .分别设计 hanning hamming blackmail kaiser 窗函数 highpass_FIRfunction highpassfilterclc;clear all;Fs=100;%采样频率fs=35;%高通阻带模拟截止频率fp=40;%高通通带模拟起始频率ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;wn=(wp+ws)/2/pi;Bt=wp-ws;N0=ceil(55*pi/Bt);N=N0+mod(N0+1,2);%调用窗函数windowl =hanning(N);window2=hamming(N);window

29、3=blackman(N);hLWn,beta.ftype=kaiserord(35.40,0 l,0.01 0.01JOO);window4=kaiser(n+ Lbeta);%设计加窗函数firlb 1 =fir 1 (N-1 ,wn,'highwindow 1);b2=fir 1 (N-1 ,wn,'highwindow2);b3=fir 1 (N-1 ,wn.'highwindow3);b4=firl(n,Wn；highwindow4 .'noscale');%求取频率响应|HLWl=freqz(bl,L512,2);H2,W2=freqz(b2

30、,1,512,2);H3,W3=freqz(b3,l,5122)；H4,W4=freqz(b4,1,512,2);figure(l);subplot(22.1),plot(W 1,20*log 10(abs(H 1);% 绘制频率响应图形axis(0,1,-100.100);title('高通hanning窗的频率响应图形')；xlabel('频率(Hz)');ylab< 幅值);subpk)t(2.2,2),plot( W2,20*log 10(abs(H2);% 绘制频率响应图形 axis(0J<100J00);title。高通hamming窗的

31、频率响应图形')；xlabel('频率(Hz)');ylabelC 幅值");subplot(223),plot(W3,20*k)gl0(abs(H3);% 绘制频率响应图形 axis(0Jr100.1001);title，高通blacknian窗的频率响应图形)；xlabcl('频率(Hz)');ylabd(1幅值);subpk)t(224),plot( W4.20*log 10(abs(H4);% 绘制频率响应图形 axis(0Jr100.100);titleC高通kaiser窗的频率响应图形')；xlabelf频率(Hz)'

32、;);ylabelC 幅值);T=l/Fs;L=100:%信号长度t=(O:L-l)*T;%定义时间范围和步长y=sin(2*pi*5*t)+5*sin(2*pi*15 %)+8*sin(2*pi*40*t);% 滤波前的图形NFFT = 2Anextpow2(L); % Next power of 2 from length of yY = fft(y,NFFT"L：%将时域信号变换到频域f = Fs/2*linspace(0.1,NFH72+l);% 频域采样 figure(2);plot(f,2*abs(Y( 1 :NFFT/2+1 );xlabel(,frequency/Hz

33、,);ylabel(,Amuplitude,);%滤波前频谱 title ('滤波前的频谱);%滤波后频谱%采用banning窗滤波器yyl=filter(b 1,1 ,y);% 调用滤波函数YYl=fft(yy 1 ,NFFT)/L;%进行傅里叶变换，下同。f l=Fs/2*linspace(0,1 .NFFT/2+1);figure(3);subplot(2.2r 1 ),plot(fL2*abs(YY1 (1 :NFFT/2+1) ;xlabel(,frequency/Hz,):ylabel(,Amuplituder);titleChanning窗的滤波效果')；%采用h

34、ammning窗滤波器yy2=filter(b2J,y);YY2=fft(yy2,NFFT)/L;fl =Fs/2*linspace(0,1 .NFFT/2+1);subplot(222),plot(fl,2*abs(YY2( 1 :NFFT/2+1) ;xlabel(*frequency/Hz,);ylabeI(,Amuplitude,); title(,hamming窗的滤波效果)；%采用blackman窗滤波器yy3=filter(b3,hy);YY3=fft(yy3,NFFT)/L:fl =Fs/2*linspace(0,1 .NFFT/2+1);subplot(2,23), plot

35、(f 1,2*abs(YY3( 1 :NFFT/2+1) ;xlabcl(,frequency/Hz,);ylabcI('Amuplitude,);titlefblackman窗的滤波效果%采用kaiser窗滤波器yy4=filter(b4,Ly);YY4=fft(yy4,NFFT)/L;f 1 =Fs/2*linspace(0,1 .NFFT/2+1);subplot(2.24),plot(fL2*abs(YY4( 1 :NFFT/2+1) ;xlabel(,frequency/Hz,);ylabel(,Amuplituder);xlabeK'frequency/Hzylab

36、elCAmuplitude*);titleCkaiser窗函数滤波效果')；%滤波前后的信号的时域对比figure(4);plot(y);xlabel(时间/s);ylabel(,振幅滤波前振幅特性)figure(5);subplot(221),plot(yyl)；xabel(时间/s);ylabel(振幅。;titleChanning 窗函数滤波振幅特性')； subplot(222),pk)t(yy2);xlabel(时间/s);ylabel(振幅)titleChamming 窗函数滤波振幅特性)； subplot(2,2,3),plot(yy3);xlabel(时间/s，

37、);ylabelf振幅);title(blackman 窗函数滤波振幅特性)； subplot(224),plot(yy4)：xlabel(时间/s');ylabel('振幅');title('kaiser 窗函数滤波振幅特性)； %滤波前后的信号的相位对比figure(6);plot(angle(Y);xlabelC时间/s');ylabel('相位);title。滤波前的相位特性figure(7);subplot(2,2.1 ),plot(angle(YYl);xlabel('时间/s');ylabcl('相位'

38、;);title('hanning 窗函数滤波相位特 性)；subpk)t(2,2,2),plot(angle(YY2);xlabel(时间/sRylabel。相位XtitleChamming 窗函数滤波相位特 性)；subplot(223),pk)t(angle(YY3);xlabelC时间/sr);ylabel(，相位RtitleCblackman 窗函数滤波相位特 性)；subplot(224),plot(angle(YY4);xlabelC时间/s);ylabe如相位RtitleCkaiser 窗函数滤波相位特性三、运行结果1 .给定的简谐信号：y = sin(2r) + 5

39、sin(2?) + 8 sin(2Z)其 W； =5z,于2 =15Hz,于3 =40Hz图一：输入简谐信号滤波前的频谱图二：输入简谐信号滤波前的振幅图三：输入简谐信号滤波前的相位涌波前的相位特性4低通blackman®的频率响应图形频率（Hz）低通kaiser®的频率响应图形 100频率（Hz）020406080100120140时间/S2 .低通滤波器的设计低通滤波器的技术指标：采样频率Fs= 100Hz：通带截止频率fp=20Hz；阻带起始频率fs=30HzHanning Hamming Blackman Kaiser采用相同的技术指标。以下即是四个窗函数的频响图及

40、对简谐信号滤波后的效果图。图四不同低通窗函数低通滤波器的归一化频响图500-50hanning窗的源波效果6.204060frequency/Hzblackman窗的滤波效果 64 28pn】_-dnE<02040frequency/Hz04 2 9pn1-dnE4hamming®的源波效果 600204060frequency/Hzkaiser©函数漉波效果15O8pn=-dnE40600204060frequency/Hz图六不同窗函数低通滤波器对信号的滤波后时间域效果banning©函数漉波振幅特性10rhamming窗函数漉波振幅特性 10-101c

41、050100时间/Sblackman窗函数源波振幅特性10»-10 050100时间/skaise僮函数漉波振幅特性2-101050100时间/s-1 050100时间/s图七不同窗函数低通滤波器对信号的滤波后相位变化banning®函数漉波相位特性 hamming窗因数漉波相位特性4«- 4 cc050100150时间/Sblackman窗函数源波相位特性4,- 4cc050100150时间/Skaiser窗函数漉波相位特性4- 4 cc050100150时间/S-4 050100150时间Zs由以上滤波后频率，相位，振幅变化以观察到：hanning窗hammi

42、ng窗及blackman窗 的滤波效果基本相当，但三者相比：hamming窗的过渡带衰减最快，blackman窗旁箫幅度 最小。而kaiser窗只有5Hz信号，15Hz信号被截断，与设计要求有出入。图十不同窗函数带通滤波器对信号的滤波后时间域振幅效果3.带通滤波器的设计采样频率Fs=100Hz：阻带截止频率1： fsl = 10Hz：通带起始频率1： fpl = 15：通带截止频率2： fp2=20；阻带截止频率 2： fs2=25Hz：Hanning Hamming Blackman Kaiser采用相同的技术指标。以下即是四个窗函数的频响图及 对简谐信号滤波后的效果图。图八不同带通窗函数低

43、通滤波器的归一化频响图帝通banning®的频率响应图形50带通hamming®的频率响应图形 1000带通blackman窗的频枣晌应图形-100 0频率（Hz）帝通kaiser®的频率响应图形 100500.5频室(Hz)图九不同窗函数带通滤波器对信号的滤波后频率域效果9Pm=dnE<<Dpn=ldnEVbanning©的滤波效果 0.80.60.40.200204060frequency/Hzblackman©的滤波效果0.80.60.40.200204060frequency/Hzhamming面的滤波效果 0.8 0.6 0.4 0.200204060frequency/Hzkaiser函数滤波效果0.80.6 .0.4 ,I；0.2 / j.pl J V-_-0204060frequency/Hzbanning®函数漉波振幅特性hamming窗函数漉波振幅特性blackman©函数滤波振幅特性kaiser窗函数漉波振幅特性图十一不同窗函数