数字信号处理实验五FIR滤波器设计汇总.docx
《数字信号处理实验五FIR滤波器设计汇总.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验五FIR滤波器设计汇总.docx（43页珍藏版）》请在冰豆网上搜索。
数字信号处理实验五FIR滤波器设计汇总
实验五FIR滤波器设计
一．实验内容
（1）认真复习FIR数字滤波器的基本概念，线性相位FIR滤波器的条件和特点、幅度函数特点、零点位置的基本特点与性质；窗函数设计法的基本概念与方法，各种窗函数的性能和设计步骤，线性相位FIR低通、高通、带通和带阻滤波器的设计方法，频率采样设计法的基本概念和线性相位的实现方法。
（2）掌握几种线性相位的特点，熟悉和掌握矩形窗、三角形窗、汉宁窗、海明窗、布莱克曼窗、凯塞窗设计IIR数字滤波器的方法，熟悉和掌握频率抽样设计法的线性相位的设计方法，并对各种线性相位的频率抽样法的设计给出调整和改进。
（3）熟悉利用MATLAB进行各类FIR数字滤波器的设计方法。
二．实验内容
a.设线性相位FIR滤波器单位抽样响应分别为
h（n）={-4,1,-1,-2,5,6,5,-2,-1,1,-4}
h（n）={-4,1,-1,-2,5,6,6,5,-2,-1,1,-4}
h（n）={-4,1,-1,-2,5,0,-5,2,1,-1,4}
h（n）={-4,1,-1,-2,5,6,-6,-5,2,1,-1,4}
分别求出滤波器的幅度频率响应H（ω），系统函数H（z）以及零极点分布，并绘制相应的波形和分布图。
在matlab中新建函数amplres，代码如下:
function[A,w,type,tao]=amplres（h）
N=length（h）;tao=（N-1）/2;
L=floor（（N-1）/2）;
n=1:
L+1;
w=[0:
500]*2*pi/500;
ifall（abs（h（n）-h（N-n+1））<1e-10）
A=2*h（n）*cos（（（N+1）/2-n）'*w）-mod（N,2）*h（L+1）;
type=2-mod（N,2）;
elseifall（abs（h（n）+h（N-n+1））<1e-10）&（h（L+1）*mod（N,2）==0）
A=2*h（n）*sin（（（N+1）/2-n）'*w）;
type=4-mod（N,2）;
elseerror（'错误，这不是线性相位滤波器!
'）
end
对第一个单位抽样响应，在matlab中新建函数a1，代码如下：
h1=[-4,1,-1,-2,5,6,5,-2,-1,1,-4];
M=length（h1）;n=0:
M-1;
[A,w,type,tao]=amplres（h1）;type
subplot（2,1,1）,stem（n,h1）;
title（'冲激响应h1'）;
ylabel（'h（n）'）;xlabel（'n'）;
subplot（2,1,2）,plot（w/pi,A）;
ylabel（'A'）;xlabel（'\pi'）;
title（'·幅频响应'）;
figure
rz=roots（h1）
fori=1:
8
r（i）=1/rz（i）;
end
r'
zplane（h1,1）;
title（'h1零极点图'）;
生成结果如下：
>>a1
type=
1
rz=
-0.9807+0.1956i
-0.9807-0.1956i
-0.5578+0.8300i
-0.5578-0.8300i
0.4052+1.2374i
0.4052-1.2374i
1.2169+0.0000i
0.8218+0.0000i
0.2390+0.7299i
0.2390-0.7299i
ans=
-0.9807+0.1956i
-0.9807-0.1956i
-0.5578+0.8300i
-0.5578-0.8300i
0.2390+0.7299i
0.2390-0.7299i
0.8218+0.0000i
1.2169+0.0000i
对第二个单位抽样响应，在matlab中新建函数a2，代码如下：
h2=[-4,1,-1,-2,5,6,6,5,-2,-1,1,-4];
M=length（h2）;n=0:
M-1;
[A,w,type,tao]=amplres（h2）;type
subplot（2,1,1）,stem（n,h2）;
title（'冲激响应h2'）;
ylabel（'h（n）'）;xlabel（'n'）;
subplot（2,1,2）,plot（w/pi,A）;
ylabel（'A'）;xlabel（'\pi'）;
title（'幅频响应'）;
figure
rz=roots（h2）
fori=1:
8
r（i）=1/rz（i）;
end
r'
zplane（h2,1）;
title（'·h2零极点图'）;
生成结果如下：
>>a2
type=
2
rz=
1.3120+0.0000i
-1.0000+0.0000i
-0.8868+0.4622i
-0.8868-0.4622i
0.4778+1.1851i
0.4778-1.1851i
-0.2957+0.9553i
-0.2957-0.9553i
0.7622+0.0000i
0.2926+0.7258i
0.2926-0.7258i
ans=
0.7622+0.0000i
-1.0000+0.0000i
-0.8868+0.4622i
-0.8868-0.4622i
0.2926+0.7258i
0.2926-0.7258i
-0.2957+0.9553i
-0.2957-0.9553i
对第三个单位抽样响应，在matlab中新建函数a3，代码如下：
h3=[-4,1,-1,-2,5,0,-5,2,1,-1,4];
M=length（h3）;n=0:
M-1;
[A,w,type,tao]=amplres（h3）;type
subplot（2,1,1）,stem（n,h3）;
title（'冲激响应h3'）;
ylabel（'h（n）'）;xlabel（'n'）;
subplot（2,1,2）,plot（w/pi,A）;
ylabel（'A'）;xlabel（'\pi'）;
title（'幅频响应¦'）;
figure
rz=roots（h3）
fori=1:
8
r（i）=1/rz（i）;
end
r'
zplane（h3,1）;
title（'h3零极点图'）;
生成结果如下：
>>a3
type=
3
rz=
-1.0000+0.0000i
-0.8732+0.4874i
-0.8732-0.4874i
0.1010+1.2041i
0.1010-1.2041i
1.0000+0.0000i
0.8280+0.5607i
0.8280-0.5607i
0.0692+0.8247i
0.0692-0.8247i
ans=
-1.0000+0.0000i
-0.8732+0.4874i
-0.8732-0.4874i
0.0692+0.8247i
0.0692-0.8247i
1.0000+0.0000i
0.8280+0.5607i
0.8280-0.5607i
对第四个单位抽样响应，在matlab中新建函数a4，代码如下：
h4=[-4,1,-1,-2,5,6,-6,-5,2,1,-1,4];
M=length（h4）;n=0:
M-1;
[A,w,type,tao]=amplres（h4）;type
subplot（2,1,1）,stem（n,h4）;
title（'冲激响应h4'）;
ylabel（'h（n）'）;xlabel（'n'）;
subplot（2,1,2）,plot（w/pi,A）;
ylabel（'A'）;xlabel（'\pi'）;
title（'·幅频响应¦'）;
figure
rz=roots（h4）
fori=1:
8
r（i）=1/rz（i）;
end
r'
zplane（h4,1）;
title（'h4零极点图'）;
生成结果如下：
>>a4
type=
4
rz=
0.2631+1.3394i
0.2631-1.3394i
-0.9505+0.3108i
-0.9505-0.3108i
-0.7309+0.6825i
-0.7309-0.6825i
0.9021+0.4315i
0.9021-0.4315i
1.0000+0.0000i
0.1412+0.7189i
0.1412-0.7189i
ans=
0.1412+0.7189i
0.1412-0.7189i
-0.9505+0.3108i
-0.9505-0.3108i
-0.7309+0.6825i
-0.7309-0.6825i
0.9021+0.4315i
0.9021-0.4315i
b.设计FIR数字低通滤波器，技术指标为：
ωp=0.2π，ωst=0.3π，δ1=0.25dB，δ2=50dB。
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
（3）选择凯塞窗函数设计该滤波器，并绘制相应的波形图。
在matlab中新建函数ideal_lp，代码如下：
functionhd=ideal_lp（wc,M）;
%IdealLowPassfiltercomputation
%--------------------------------
%[hd]=ideal_lp（wc,M）;
%hd=idealimpulseresponsebetween0toM-1
%wc=cutofffrequencyinradians
%M=lengthoftheidealfilter
%
alpha=（M-1）/2;
n=[0:
1:
（M-1）];
m=n-alpha+eps;%addsmallestnumbertoavoidividedbyzero
hd=sin（wc*m）./（pi*m）;
end
在matlab中新建函数freqz_m，代码如下：
function[db,mag,pha,grd,w]=freqz_m（b,a）;
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
1:
501））';w=（w（1:
1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
grd=grpdelay（b,a,w）;
在matlab中新建函数b1，代码如下：
%用Hamming窗函数设计FIR数字滤波器
wp=0.2*pi;ws=0.3*pi
N=61
n=[0:
1:
N-1]
wc=（ws+wp）/2;%理想低通滤波器
hd=ideal_lp（wc,N）;%理想低通的冲激响应
w_ham=（hamming（N））';
h=hd.*w_ham;%FIR滤波器冲激响应
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））%实际的通带衰减
As=-round（max（db（ws/delta_w+1:
1:
501）））%实际的最小阻带衰减
subplot（221）;stem（n,hd）;title（'理想冲激响应'）
axis（[0N-1-0.10.3]）;xlabel（'n'）;ylabel（'hd（n）'）
subplot（222）;stem（n,w_ham）;title（'hamming窗'）
axis（[0N-101.1]）;xlabel（'n'）;ylabel（'w（n）'）;
subplot（223）;stem（n,h）;title（'实际冲激响应'）
axis（[0N-1-0.10.3]）;xlabel（'n'）;ylabel（'h（n）'）
subplot（224）;plot（w/pi,db）;
axis（[00.8-1000]）;xlabel（'以PI为单位的频率'）;ylabel（'对数幅度/db'）;
实验结果如下：
在matlab中新建函数b2，代码如下：
%用Kaiser窗函数设计FIR数字滤波器
wp=0.2*pi;ws=0.3*pi;As=50
tr_width=ws-wp
N=ceil（（As-7.95）/（14.36*tr_width/（2*pi））+1）+1
n=[0:
1:
N-1]
beta=0.1102*（As-8.7）
wc=（wp+ws）/2%理想低通的截止频率
hd=ideal_lp（wc,N）
w_kai=（kaiser（N,beta））'
h=hd.*w_kai
[db,mag,pha,grd,w]=freqz_m（h,[1]）
delta_w=2*pi/1000
Rp=-（min（db（1:
1:
wp/delta_w+1）））%Ê实际的通带衰减
As=-round（max（db（ws/delta_w+1:
1:
501）））%实际的最小阻带衰减
subplot（211）;plot（w/pi,db）;title（'凯森窗幅度响应（dB）'）;grid
axis（[00.5-1000]）
ylabel（'对数幅度/db'）;xlabel（'以/pi为单位的频率'）
subplot（212）;plot（w/pi,pha）;title（'相位响应'）;grid
axis（[00.5-44]）
ylabel（'相位'）;xlabel（'以/pi为单位的频率'）
实验结果如下：
c.设计FIR数字带通滤波器，技术指标为：
下阻带边缘：
ωst1=0.2π，δs1=60dB，下通带边缘：
ωp1=0.35π，δp1=1dB；上通带边缘：
ωp2=0.65π，δp1=1dB，上阻带边缘：
ωst2=0.8π，δs2=60dB；
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
在matlab中新建函数c，代码如下：
ws1=0.2*pi;
wp1=0.35*pi;
ws2=0.8*pi;
wp2=0.65*pi;
Ap=60;
Rp=1;
tr_width=min（（wp1-ws1）,（ws2-wp2））;
M=ceil（11*pi/tr_width）;
n=[0:
1:
M-1];
wc1=（ws1+wp1）/2;
wc2=（wp2+ws2）/2;
hd=ideal_lp（wc2,M）-ideal_lp（wc1,M）;
w_bla=（blackman（M））';
h=hd.*w_bla;
[H,W]=freqz（h,1）;
subplot（2,2,1）;stem（n,hd）;title（'理想脉冲抽样'）;
subplot（2,2,2）;stem（n,w_bla）;title（'布莱克曼窗'）;
subplot（2,2,3）;stem（n,h）;title（'实际脉冲抽样'）;
subplot（2,2,4）;plot（W/pi,20*log10（abs（H）））;title（'幅度响应（db）'）;
生成结果如下：
d.设计FIR数字带通滤波器，技术指标为：
下阻带边缘：
ωst1=0.2π，δs1=60dB，下通带边缘：
ωp1=0.4π，δp1=1dB；上通带边缘：
ωp2=0.6π，δp1=1dB，上阻带边缘：
ωst2=0.8π，δs2=60dB；
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
在matlab中新建函数d，代码如下：
ws1=0.2*pi;
wp1=0.4*pi;
ws2=0.8*pi;
wp2=0.6*pi;
Ap=60;
Rp=1;
tr_width=min（（wp1-ws1）,（ws2-wp2））;
M=ceil（11*pi/tr_width）;
n=[0:
1:
M-1];
wc1=（ws1+wp1）/2;
wc2=（wp2+ws2）/2;
hd=ideal_lp（wc2,M）-ideal_lp（wc1,M）;
w_bla=（hamming（M））';
h=hd.*w_bla;
[H,W]=freqz（h,1）;
subplot（2,2,1）;stem（n,hd）;title（'理想脉冲抽样'）;
subplot（2,2,2）;stem（n,w_bla）;title（'海明窗'）;
subplot（2,2,3）;stem（n,h）;title（'实际脉冲抽样'）;
subplot（2,2,4）;plot（W/pi,20*log10（abs（H）））;title（'幅度响应（db）'）;
生成结果如下：
e.设计FIR数字带通滤波器，技术指标为：
下阻带边缘：
ωst1=0.2π，δs1=20dB，下通带边缘：
ωp1=0.4π，δp1=1dB；上通带边缘：
ωp2=0.6π，δp1=1dB，上阻带边缘：
ωst2=0.8π，δs2=20dB；
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
在matlab中新建函数e，代码如下：
ws1=0.2*pi;
wp1=0.4*pi;
ws2=0.8*pi;
wp2=0.6*pi;
Ap=20;
Rp=1;
tr_width=min（（wp1-ws1）,（ws2-wp2））;
M=ceil（11*pi/tr_width）;
n=[0:
1:
M-1];
wc1=（ws1+wp1）/2;
wc2=（wp2+ws2）/2;
hd=ideal_lp（wc2,M）-ideal_lp（wc1,M）;
w_bla=（boxcar（M））';
h=hd.*w_bla;
[H,W]=freqz（h,1）;
subplot（2,2,1）;stem（n,hd）;title（'理想脉冲抽样'）;
subplot（2,2,2）;stem（n,w_bla）;title（'矩形窗'）;
subplot（2,2,3）;stem（n,h）;title（'实际脉冲抽样'）;
subplot（2,2,4）;plot（W/pi,20*log10（abs（H）））;title（'幅度响应（db）'）;
生成结果如下：
f.设计FIR数字高通滤波器，技术指标为：
通带截止频率为ωp=15π/27，阻带截止频率为ωst=11π/27，通带最大衰减为δ1=2.5dB，阻带最小衰减为δ2=55dB。
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
系统一：
在matlab中新建函数f，代码如下：
As=55;
ws=11*pi/27;
wp=15*pi/27;
tr_width=wp-ws;%计算过渡带
M=ceil（（As-7.95）*2*pi/（14.36*tr_width）+1）+1;%按凯泽窗计算滤波器的长度
disp（['滤波器的长度',num2str（M）]）;
beta=0.1102*（As-8.7）;%计算凯泽窗的beta值
n=[0:
1:
M-1];
disp（['线性相位滤波器',num2str（beta）]）;
w_kai=（kaiser（M,beta））';%求凯泽窗函数
wc=（ws+wp）/2;
hd=ideal_lp（pi,M）-ideal_lp（wc,M）;%求理想脉冲响应
h=hd.*w_kai;
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（wp/delta_w+1:
1:
501）））;
disp（['实际通带波动为',num2str（Rp）]）;
As=-round（max（db（1:
1:
ws/delta_w+1）））;
disp（['最小阻带衰减为ª',num2str（As）]）;
subplot（2,2,1）;
stem（n,hd）;
title（'理想脉冲响应'）;
axis（[0M-1-0.40.8]）;
ylabel（'hd（n）'）;
subplot（2,2,2）;
stem（n,w_kai）;
title（'凯泽窗'）;
axis（[0M-101.1]）;
ylabel（'wd（n）'）;
subplot（2,2,3）;
stem（n,h）;
title（'实际脉冲响应'）;
axis（[0M-1-0.40.8]）;
xlabel（'n'）;ylabel（'h（n）'）;
subplot（2,2,4）;
plot（w/pi,db）;
title（'·幅度响应/dB'）;
axis（[01-10010]）;
grid;
xlabel（'以pi为单位的频率'）;
ylabel（'·分贝数/dB'）;
生成结果如下：
滤波器的长度47
线性相位滤波器5.1023
实际通带波动为0.025017
最小阻带衰减为?
55
滤波器的长度47
线性相位滤波器5.1023
实际通带波动为0.025017
最小阻带衰减为55
g.设计FIR数字高通滤波器，技术指标为：
通带截止频率为ωp=0.6π，阻带截止频率为ωst=0.4π，通带最大衰减为δ1=0.25dB，阻带最小衰减为δ2=40dB。
（1）通过技术指标，选择一种窗函数进行设计；
（2）求滤波器的单位抽样响应、频率响应，并绘制波形。
在matlab中新建函数c，代码如下：
function[db,mag,pha,w]=freqz_m2（b,a）;
%Modifiedversionoffreqzsubroutine
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';w=（w（1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
在matlab中新建函数g，代码如下：
Wp=0.6*pi;
Ws=0.4*pi;
tr_width=Wp-Ws;
M=ceil（6.2*pi/tr_width）;
n=0:
1:
M-1;
Wc=（Ws+Wp）/
升级会员 