MATLAB窗函数法实现FIR的高通带通和低通滤波器的程序.docx
《MATLAB窗函数法实现FIR的高通带通和低通滤波器的程序.docx》由会员分享，可在线阅读，更多相关《MATLAB窗函数法实现FIR的高通带通和低通滤波器的程序.docx（33页珍藏版）》请在冰点文库上搜索。
MATLAB窗函数法实现FIR的高通带通和低通滤波器的程序
MATLAB课程设计报告
学院：
地球物理与石油资源学院
班级：
测井（基）11井1
姓名：
大牛啊啊啊
学号：
班内编号：
指导教师：
陈义群
完成日期：
2013年6月3日
题目FIR滤波器的窗函数设计法及性能比较
1.FIR滤波器简介
数字波波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。
根据其单位冲激响应函数的时域特性可分为两类：
无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。
与IIR滤波器相比，FIR滤波器的主要特点为：
a.线性相位；b.非递归运算。
2.FIR滤波器的设计
FIR滤波器的设计方法主要有三种：
a.窗函数设计法；b.频率抽样发：
c.最小平法抽样法:
这里我主要讨论在MATLAB环境下通过调用信号分析与处理工具箱的几类窗函数来设计滤波器并分析与比较其性能。
窗函数法设计FIR滤波器的一般步骤如下：
a.根据实际问题确定要设计的滤波器类型：
b.根据给定的技术指标，确定期望滤波器的理想频率特性；
c.求期望滤波器的单位脉冲响应：
d.求数字滤波器的单位脉冲响应；
e.应用。
常用的窗函数有
（1）Hanningwindow
（9.2.2）
（9.2.3）
w（n}=[0.5-0.5cos（-^，l）]RN（n）
N—1
（2）Hamnungwindow
w（n）=[0.54-0.46cos（-^-^-）]7?
v（/?
）
N—1
3.）Balckmanwindow
w（〃）=[0,42-0.5cos（^^）+0.08cos（i^^^）]Rv（〃）（9.2.4）
N—1N-1
（4）Kaiserwindow
式中/。
⑴是零阶8css”函数，可定义为
当x=0时与矩形窗一致；当,=5.4414H寸与海明窗结果相同;
当x=8.885时与布莱克曼窗结果札同。
4.常用窗函数的参数
常用窗函数的几个参照数据
窗类型
旁瓣峰值（dB）
主辩宽度
最小阻带衰减（dB）
矩形窗
-13
4%AV
-21
三角窗
-25
8%/丫
-25
汉宁窗
-31
脉/N
-44
哈明窗
-41
WN
-53
布莱克曼窗
-57
12，N
-74
5.FIR滤波器的MATLAB实现方式
在MATLAB信号分析与处理工具箱中提供了大量FIR窗函数的设计函数，本次用到主
计算kaiser窗函数的相关参数求取频率响应
对信号进行滤波的函数
本次实验分别通过调用hanning,hamming,Blackman,kaiser窗函数，给以相同
的技术参数，来设计低通，带通，高通滤波器，用上述窗函数的选择标准来比较各种窗函数的优劣，并给以一个简谐波进行滤波处理，比较滤波前后的效果。
达到综合比较的
效果。
二、源代码
L利用banninghammingblackmailkaiser窗,设计一个低通FIRfunctionlowpassfilterclc;
clearall;
Fs=100;%采样频率
fp=2O;%通带截止频率
fs=3O;%阻带起始频率
wp=2*pi*fp/Fs;%将模拟通带截止频率转换为数字滤波器频率ws=2*pi*fs/Fs;%将模拟阻带起始频率转换为数字波波器频率wn=（wp+ws）/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（[2025],[1OJJO.O1O.O1JJOO）;
window4=kaiser（n+l,beta）;%使用kaiser窗函数
%设计加窗函数fkl
b1=fir1（N-1,wn.window1）;
b2=fir1（N-1,wn,window2）;
b3=fir1（N-1avilwindow3）;
b4=fir1（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.2J）.plot（WL20*log10（abs（H1）））;%绘制频率响应图形axis（[0Jr100J00]）;
title。
低通banning窗的频率响应图形,）;
xlabelf频率（Hz）'）;
ylabdC幅值,）;
subplot（252,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形axis（[OJrlOOJOOJ）;
titled氐通hamming窗的频率响应图形'）；
xlabel（'频率（Hz））;
ylabdC幅值）;
subplot（223）.plot（W3,20*logl0（abs（H3）））;%绘制频率响应图形axis（[0Jr100,100]）;
title。
低通blackman窗的频率响应图形'）；
xlabelf频率（Hz）'）;
ylabel（,幅值）;
subplot（22,4）.plot（W4,20*log10（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）;%Nextpowerof2fromlengthofy
Y=fft（y,NFFT"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（b1,1,y）;%调用滤波函数
YYl=fft（yy1,NFFT）/L;%进行傅里叶变换，下同。
fl=Fs/2*linspace（0.1.NFFT/2+1）;
figured）;
subplot（2.2,1）,plot（f1,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,）;title（力amming窗的滤波效果）；
%采用blackman窗滤波器
yy3=filter（b3J,y）;
YY3=fft（yy3,NFFT）/L：
fl=Fs/2*linspace（0.1.NFFT/2+1）;
subplot（2,23hplot（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（0J.NFFT/2+1）;
subplot（2,24）,plot（f1,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）；xlabel（时间/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（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）,plot（angle（YY3））;xlabel（时间/sRylabel。
相位XtitleCblackman窗函数滤波相位特性）；
subplot（224）,plot（angle（YY4））;xlabelC时间/s）;ylabel（'相位）;titleCkaiser窗函数滤波相位特性
2.设-个lianninghammingblackmailkaiser窗函数bandpass_FIR
%设计•一个banninghammingblackmankaiser窗函数bandpass_FIR
functionbandpassfilter
Fs=100;%采样频率
fpl=15;%通带下限截止频率
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=[（wp1+ws1）/2/pi,（wp2+ws2）/2/pi];
windowl=hanning（N）;%使用banning窗函数
window2=haniming（N）;%使用hamming窗函数
window3=blackn】an（N）;%使用blackman窗函数
%设过渡带宽度为5Hz
[n,Wn,beta,ftype]=kaiserord（[10152025].[010]5[0.010.010.01]JOO）;%求阶数n以及参数beta
window4=kaiser（n+l,beta）;%使用kaiser窗函数
%设计加窗函数fi”
b1=fir1（N-1,wn.window1）;
b2=fir1（N-1,wn.window2）;
b3=fir1（N-1,wn,window3）;
b4=firl（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.2J）.plot（WL20*log10（abs（H1）））;%绘制频率响应图形
axis（[04<100,100]）;
title。
带通banning窗的频率响应图形'）；
xlabelC频率（Hz）'）;
ylabdC幅值"）;
subplot（252,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形
axis（[OJr100,100]）;
title。
带通hamming窗的频率响应图形'）；
xlabel（'频率（Hz））;
ylabelC幅值）;
subpk）t（223）.plot（W3,20*logl0（abs（H3）））;%绘制频率响应图形
axis（[0Jr100,100]）;
title。
带通blackman窗的频率响应图形
xlabelC频率（Hz）'）;
ylabelf幅值）;
subplot（2,2,4）,plot（W4,20*log10（abs（H4）））;%绘制频率响应图形
axis（[0Jr100,100|）;
title。
带通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）;%Nextpowerof2fromlengthofy
Y=fft（y,NFFT）/L;%将时域信号变换到频域
f=Fs/2*linspace（0J.NFFT/2+1）;%频域采样
figure
（2）;
plot（f,2*abs（Y（1:
NFFT/2+1）））;xlabel（,frequency/Hz,）;ylabel（,Amuplitude,）;%滤波前频谱title（'滤波前的频谱'）；
%滤波后频谱
%采用banning窗波波器
yyl=filter（b1,1,y）;%调用滤波函数
YYl=fft（yyLNFFT）/L;%进行傅里叶变换，下同。
fl=Fs/2*linspace（0.1.NFFT/2+1）;
figured）;
subplot（2,24）,plot（f1,2*abs（YYl（l:
NFFT/2+1）））;xIabel（,frequency/Hz,）;yIabel（,Amuplitude,）;
titleChanning窗的滤波效果'）；
%采用hammning窗滤波器
yy2=filter（b2J,y）;
YY2=fft（yy2.NFFT）/L;
fl=Fs/2*linspace（0,1.NFFT/2+1）;
subplot（2,22）,plot（f1,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.1.NFFT/2+1）;
subplot（223）,plot（f1,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,）;
xlabel（,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（'haniming窗函数滤波振幅特性）;
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（'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（湘位"itleCkaiser窗函数滤波相位特性'）;
3.分别设计hanninghammingblackmailkaiser窗函数highpass_FIR
functionhighpassfilter
clc;
clearall;
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）;
window3=blackman（N）;
hLWn,beta.ftype]=kaiserord（[35.40],[0l],[0.010.01]JOO）;
window4=kaiser（n+Lbeta）;
%设计加窗函数firl
b1=fir1（N-1,wn,'high\window1）;
b2=fir1（N-1,wn,'high\window2）;
b3=fir1（N-1,wn.'high\window3）;
b4=firl（n,Wn；high\window4.'noscale'）;
%求取频率响应
|HLWl]=freqz（bl,L512,2）;
[H2,W2]=freqz（b2,1,512,2）;
[H3,W3]=freqz（b3,l,5122）；
[H4,W4]=freqz（b4,1,512,2）;
figure（l）;
subplot（22.1）,plot（W1,20*log10（abs（H1）））;%绘制频率响应图形
axis（[0,1,-100.100]）;
title（'高通hanning窗的频率响应图形'）；
xlabel（'频率（Hz）'）;
ylab<幅值）;
subpk）t（2.2,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形axis（[0J<100J00]）;
title。
高通hamming窗的频率响应图形'）；
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*log10（abs（H4）））;%绘制频率响应图形axis（[0Jr100.100]）;
titleC高通kaiser窗的频率响应图形'）；
xlabelf频率（Hz）'）;
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）;%Nextpowerof2fromlengthofy
Y=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,）;ylabel（,Amuplitude,）;%滤波前频谱title（'滤波前的频谱）;
%滤波后频谱
%采用banning窗滤波器
yyl=filter（b1,1,y）;%调用滤波函数
YYl=fft（yy1,NFFT）/L;%进行傅里叶
升级会员 