FIR滤波器的窗函数设计法及性能比较.docx
《FIR滤波器的窗函数设计法及性能比较.docx》由会员分享，可在线阅读，更多相关《FIR滤波器的窗函数设计法及性能比较.docx（27页珍藏版）》请在冰点文库上搜索。
FIR滤波器的窗函数设计法及性能比较
MATLAB课程设计报告
学院：
地球物理与石油资源学院
班级：
姓名：
学号：
班内编号：
指导教师：
完成日期：
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.应用。
常用的窗函数有
3．窗函数的选择标准
1.较低的旁瓣幅度，尤其是第一旁瓣；
2.旁瓣幅度要下降得快，以利于增加阻带衰减；
3.主瓣宽度要窄，这样滤波器过渡带较窄。
4.常用窗函数的参数
5.FIR滤波器的MATLAB实现方式
在MATLAB信号分析与处理工具箱中提供了大量FIR窗函数的设计函数，本次用到主要有以下几种：
hanning（N）hanning窗函数的调用
hamming（N）hamming窗函数的调用
blackman（N）blackman窗函数的调用
kaiser（n+1,beta）kaiser窗函数的调用
kaiserord计算kaiser窗函数的相关参数
freqz求取频率响应
filter对信号进行滤波的函数
6.实验具体步骤
本次实验分别通过调用hanning,hamming,Blackman,kaiser窗函数，给以相同的技术参数，来设计低通，带通，高通滤波器，用上述窗函数的选择标准来比较各种窗函数的优劣，并给以一个简谐波进行滤波处理，比较滤波前后的效果。
达到综合比较的效果。
二、源代码
1.利用hanninghammingblackmankaiser窗,设计一个低通FIR
functionlowpassfilter
clc;
clearall;
Fs=100;%采样频率
fp=20;%通带截止频率
fs=30;%阻带起始频率
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=N0+mod（N0+1,2）;
window1=hanning（N）;%使用hanning窗函数
window2=hamming（N）;%使用hamming窗函数
window3=blackman（N）;%使用blackman窗函数
[n,Wn,beta,ftype]=kaiserord（[2025],[10],[0.010.01],100）;
window4=kaiser（n+1,beta）;%使用kaiser窗函数
%设计加窗函数fir1
b1=fir1（N-1,wn,window1）;
b2=fir1（N-1,wn,window2）;
b3=fir1（N-1,wn,window3）;
b4=fir1（n,Wn/pi,window4,'noscale'）;
%求取频率响应
[H1,W1]=freqz（b1,1,512,2）;
[H2,W2]=freqz（b2,1,512,2）;
[H3,W3]=freqz（b3,1,512,2）;
[H4,W4]=freqz（b4,1,512,2）;
figure
（1）;
subplot（2,2,1）,plot（W1,20*log10（abs（H1）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'低通hanning窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'低通hamming窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,3）,plot（W3,20*log10（abs（H3）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'低通blackman窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,4）,plot（W4,20*log10（abs（H4）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'低通kaiser窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
T=1/Fs;
L=100;%信号长度
t=（0:
L-1）*T;%定义时间范围和步长
y=sin（2*pi*5*t）+5*sin（2*pi*15*t）+8*sin（2*pi*40*t）;%滤波前的图形
NFFT=2^nextpow2（L）;%Nextpowerof2fromlengthofy
Y=fft（y,NFFT）/L;%将时域信号变换到频域
f=Fs/2*linspace（0,1,NFFT/2+1）;%频域采样
figure
（2）;
plot（f,2*abs（Y（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;%滤波前频谱
title（'滤波前的频谱'）;
%滤波后频谱
%采用hanning窗滤波器
yy1=filter（b1,1,y）;%调用滤波函数
YY1=fft（yy1,NFFT）/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace（0,1,NFFT/2+1）;
figure（3）;
subplot（2,2,1）,plot（f1,2*abs（YY1（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'hanning窗的滤波效果'）;
%采用hammning窗滤波器
yy2=filter（b2,1,y）;
YY2=fft（yy2,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,2）,plot（f1,2*abs（YY2（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'hamming窗的滤波效果'）;
%采用blackman窗滤波器
yy3=filter（b3,1,y）;
YY3=fft（yy3,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,3）,plot（f1,2*abs（YY3（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'blackman窗的滤波效果'）;
%采用kaiser窗滤波器
yy4=filter（b4,1,y）;
YY4=fft（yy4,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,4）,plot（f1,2*abs（YY4（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'kaiser窗函数滤波效果'）;
%滤波前后的信号的时域对比
figure（4）;
plot（y）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'滤波前振幅特性'）;
figure（5）;
subplot（2,2,1）,plot（yy1）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'hanning窗函数滤波振幅特性'）;
subplot（2,2,2）,plot（yy2）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'hamming窗函数滤波振幅特性'）;
subplot（2,2,3）,plot（yy3）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'blackman窗函数滤波振幅特性'）;
subplot（2,2,4）,plot（yy4）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'kaiser窗函数滤波振幅特性'）;
%滤波前后的信号的相位对比
figure（6）;
plot（angle（Y））;xlabel（'时间/s'）;ylabel（'相位'）;title（'滤波前的相位特性'）;
figure（7）;
subplot（2,2,1）,plot（angle（YY1））;xlabel（'时间/s'）;ylabel（'相位'）;title（'hanning窗函数滤波相位特性'）;
subplot（2,2,2）,plot（angle（YY2））;xlabel（'时间/s'）;ylabel（'相位'）;title（'hamming窗函数滤波相位特性'）;
subplot（2,2,3）,plot（angle（YY3））;xlabel（'时间/s'）;ylabel（'相位'）;title（'blackman窗函数滤波相位特性'）;
subplot（2,2,4）,plot（angle（YY4））;xlabel（'时间/s'）;ylabel（'相位'）;title（'kaiser窗函数滤波相位特性'）;
2.设计一个hanninghammingblackmankaiser窗函数bandpass_FIR
%设计一个hanninghammingblackmankaiser窗函数bandpass_FIR
functionbandpassfilter
Fs=100;%采样频率
fp1=15;%通带下限截止频率
fp2=20;%通带上限截止频率
fs1=10;
fs2=25;
wp1=2*pi*fp1/Fs;%将通带下限截止频率转换为数字滤波器频率
wp2=2*pi*fp2/Fs;%将通带上限截止频率转换为数字滤波器频率
ws1=2*pi*fs1/Fs;%将通带下限截止频率转换为数字滤波器频率
ws2=2*pi*fs2/Fs;%将通带上限截止频率转换为数字滤波器频率
Bt=wp1-ws1;
N0=ceil（6.2*pi/Bt）;
N=N0+mod（N0+1,2）;
wn=[（wp1+ws1）/2/pi,（wp2+ws2）/2/pi];
window1=hanning（N）;%使用hanning窗函数
window2=hamming（N）;%使用hamming窗函数
window3=blackman（N）;%使用blackman窗函数
%设过渡带宽度为5Hz
[n,Wn,beta,ftype]=kaiserord（[10152025],[010],[0.010.010.01],100）;%求阶数n以及参数beta
window4=kaiser（n+1,beta）;%使用kaiser窗函数
%设计加窗函数fir1
b1=fir1（N-1,wn,window1）;
b2=fir1（N-1,wn,window2）;
b3=fir1（N-1,wn,window3）;
b4=fir1（n,Wn,window4,'noscale'）;
%求取频率响应
[H1,W1]=freqz（b1,1,512,2）;
[H2,W2]=freqz（b2,1,512,2）;
[H3,W3]=freqz（b3,1,512,2）;
[H4,W4]=freqz（b4,1,512,2）;
figure
（1）;
subplot（2,2,1）,plot（W1,20*log10（abs（H1）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'带通hanning窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'带通hamming窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,3）,plot（W3,20*log10（abs（H3）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'带通blackman窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,4）,plot（W4,20*log10（abs（H4）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'带通kaiser窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
T=1/Fs;
L=100;%信号长度
t=（0:
L-1）*T;%定义时间范围和步长
y=sin（2*pi*5*t）+5*sin（2*pi*15*t）+8*sin（2*pi*40*t）;%滤波前的图形
NFFT=2^nextpow2（L）;%Nextpowerof2fromlengthofy
Y=fft（y,NFFT）/L;%将时域信号变换到频域
f=Fs/2*linspace（0,1,NFFT/2+1）;%频域采样
figure
（2）;
plot（f,2*abs（Y（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;%滤波前频谱
title（'滤波前的频谱'）;
%滤波后频谱
%采用hanning窗滤波器
yy1=filter（b1,1,y）;%调用滤波函数
YY1=fft（yy1,NFFT）/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace（0,1,NFFT/2+1）;
figure（3）;
subplot（2,2,1）,plot（f1,2*abs（YY1（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'hanning窗的滤波效果'）;
%采用hammning窗滤波器
yy2=filter（b2,1,y）;
YY2=fft（yy2,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,2）,plot（f1,2*abs（YY2（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'hamming窗的滤波效果'）;
%采用blackman窗滤波器
yy3=filter（b3,1,y）;
YY3=fft（yy3,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,3）,plot（f1,2*abs（YY3（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'blackman窗的滤波效果'）;
%采用kaiser窗滤波器
yy4=filter（b4,1,y）;
YY4=fft（yy4,NFFT）/L;
f1=Fs/2*linspace（0,1,NFFT/2+1）;
subplot（2,2,4）,plot（f1,2*abs（YY4（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'kaiser窗函数滤波效果'）;
%滤波前后的信号的时域对比
figure（4）;
plot（y）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'滤波前振幅特性'）;
figure（5）;
subplot（2,2,1）,plot（yy1）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'hanning窗函数滤波振幅特性'）;
subplot（2,2,2）,plot（yy2）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'hamming窗函数滤波振幅特性'）;
subplot（2,2,3）,plot（yy3）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'blackman窗函数滤波振幅特性'）;
subplot（2,2,4）,plot（yy4）;xlabel（'时间/s'）;ylabel（'振幅'）;title（'kaiser窗函数滤波振幅特性'）;
%滤波前后的信号的相位对比
figure（6）;
plot（angle（Y））;xlabel（'时间/s'）;ylabel（'相位'）;title（'滤波前的相位特性'）;
figure（7）;
subplot（2,2,1）,plot（angle（YY1））;xlabel（'时间/s'）;ylabel（'相位'）;title（'hanning窗函数滤波相位特性'）;
subplot（2,2,2）,plot（angle（YY2））;xlabel（'时间/s'）;ylabel（'相位'）;title（'hamming窗函数滤波相位特性'）;
subplot（2,2,3）,plot（angle（YY3））;xlabel（'时间/s'）;ylabel（'相位'）;title（'blackman窗函数滤波相位特性'）;
subplot（2,2,4）,plot（angle（YY4））;xlabel（'时间/s'）;ylabel（'相位'）;title（'kaiser窗函数滤波相位特性'）;
3.分别设计hanninghammingblackmankaiser窗函数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）;
%调用窗函数
window1=hanning（N）;
window2=hamming（N）;
window3=blackman（N）;
[n,Wn,beta,ftype]=kaiserord（[35,40],[01],[0.010.01],100）;
window4=kaiser（n+1,beta）;
%设计加窗函数fir1
b1=fir1（N-1,wn,'high',window1）;
b2=fir1（N-1,wn,'high',window2）;
b3=fir1（N-1,wn,'high',window3）;
b4=fir1（n,Wn,'high',window4,'noscale'）;
%求取频率响应
[H1,W1]=freqz（b1,1,512,2）;
[H2,W2]=freqz（b2,1,512,2）;
[H3,W3]=freqz（b3,1,512,2）;
[H4,W4]=freqz（b4,1,512,2）;
figure
（1）;
subplot（2,2,1）,plot（W1,20*log10（abs（H1）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'高通hanning窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,2）,plot（W2,20*log10（abs（H2）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'高通hamming窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,3）,plot（W3,20*log10（abs（H3）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'高通blackman窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
subplot（2,2,4）,plot（W4,20*log10（abs（H4）））;%绘制频率响应图形
axis（[0,1,-100,100]）;
title（'高通kaiser窗的频率响应图形'）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
T=1/Fs;
L=100;%信号长度
t=（0:
L-1）*T;%定义时间范围和步长
y=sin（2*pi*5*t）+5*sin（2*pi*15*t）+8*sin（2*pi*40*t）;%滤波前的图形
NFFT=2^nextpow2（L）;%Nextpowerof2fromlengthofy
Y=fft（y,NFFT）/L;%将时域信号变换到频域
f=Fs/2*linspace（0,1,NFFT/2+1）;%频域采样
figure
（2）;
plot（f,2*abs（Y（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;%滤波前频谱
title（'滤波前的频谱'）;
%滤波后频谱
%采用hanning窗滤波器
yy1=filter（b1,1,y）;%调用滤波函数
YY1=fft（yy1,NFFT）/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace（0,1,NFFT/2+1）;
figure（3）;
subplot（2,2,1）,plot（f1,2*abs（YY1（1:
NFFT/2+1）））;xlabel（'frequency/Hz'）;ylabel（'Amuplitude'）;
title（'hanning窗的滤波效
升级会员 