数字信号处理实验三.docx
《数字信号处理实验三.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验三.docx（27页珍藏版）》请在冰点文库上搜索。
数字信号处理实验三
《数字信号处理》上机实验三实验报告
学号：
14020710014姓名：
黄新凯
IIR滤波器设计
（1）用matlab确定一个数字IIR低通滤波器所有四种类型的最低阶数。
指标如下：
40kHz的采样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
（1）程序：
fc=40;fp=4;fs=8;rp=0.5;rs=40;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'buttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'low','s'）;
w=0:
0.002:
5;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'b-'）
disp（'cheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'r-'）
disp（'cheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'ellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'g-'）
legend（'butter','cheby1','cheby2','ellip'）
gridon
xlabel（'w'）
ylabel（'h'）
运行结果：
buttord
n=
9
wc=
0.7533
cheb1ord
n=
5
wpo=
0.6283
cheb2ord
n=
5
wso=
1.2069
ellipord
n=
4
wc=
0.6283
（2）用matlab确定一个数字IIR高通滤波器所有四种类型的最低阶数。
指标如下：
3500Hz的采样率，1050Hz的通带边界频率，600Hz的阻带边界频率，1dB的通带波纹，50dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
程序：
fc=3500;fp=1050;fs=600;rp=1;rs=50;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'buttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'high','s'）;
w=0:
0.001:
6;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'b-'）
disp（'cheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'r-'）
disp（'cheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'ellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'g-'）
legend（'butter','cheby1','cheby2','ellip'）
axis（[0,5,-120,0]）
gridon
xlabel（'w'）
ylabel（'h'）
运行结果：
buttord
n=
12
wc=
1.7402
cheb1ord
n=
7
wpo=
1.8850
cheb2ord
n=
7
wso=
1.2049
ellipord
n=
5
wc=
1.8850
（3）用matlab确定一个数字IIR带通滤波器所有四种类型的最低阶数。
指标如下：
7kHz的采样率，1.4kHz和2.1kHz的通带边界频率，1.05kHz和2.45kHz的阻带边界频率，0.4dB的通带波纹，50dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
程序：
fc=7;fp=[1.4,2.1];fs=[1.05,2.45];rp=0.4;rs=50;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'buttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'s'）;
w=0:
0.002:
5;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'b-'）
disp（'cheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'r-'）
disp（'cheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'ellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'g-'）
legend（'butter','cheby1','cheby2','ellip'）
axis（[0,5,-120,20]）
gridon
xlabel（'w'）
ylabel（'h'）
buttord
运行结果：
n=
12
wc=
1.23051.9250
cheb1ord
n=
7
wpo=
1.25661.8850
cheb2ord
n=
7
wso=
1.10502.1437
ellipord
n=
5
wc=
1.25661.8850
（4）用matlab确定一个数字IIR带阻滤波器所有四种类型的最低阶数。
指标如下：
12kHz的采样率，2.1kHz和4.5kHz的通带边界频率，2.7kHz和3.9kHz的阻带边界频率，0.6dB的通带波纹，45dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
程序：
fc=12;fp=[2.1,4.5];fs=[2.7,3.9];rp=0.5;rs=40;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'buttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'stop','s'）;
w=0:
0.002:
5;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'b-'）
disp（'cheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'stop','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'r-'）
disp（'cheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'stop','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'ellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'stop','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'g-'）
legend（'butter','cheby1','cheby2','ellip',4）
gridon
axis（[0,5,-120,20]）
xlabel（'w'）
ylabel（'h'）
运行结果：
buttord
n=
10
wc=
1.27262.2684
cheb1ord
n=
6
wpo=
1.22522.3561
cheb2ord
n=
6
wso=
1.38452.0851
ellipord
n=
5
wc=
1.9962.3562
用到的函数：
butter，buttord，cheb2ord，chebl1，cheby2，
ellip，ellipord.
2.FIR滤波器设计
分别用矩形窗，Blackman窗，Hamming窗，Hanning窗和Bartlett窗设计截止频率为0.3Pi，窗长为M（M=11，41，81，121）的FIR低通滤波器。
在图中画出：
（1）理想低通滤波器的冲激响应；
程序：
wc=0.3*pi;
M=[11,41,81,121];
fori=1:
1:
4
tao=（M（i）-1）/2;
h_d1n=zeros（1,N（i））;
forn=1:
1:
N（i）
ifn==tao
h_d1n（n）=wc/pi;
else
h_d1n（n）=（sin（wc*（n-tao）））/（pi*（n-tao））;
end
end
subplot（2,2,i）
n1=1:
1:
M（i）;
stem（n1,h_d1n,'b'）;
title（{['a=',num2str（i）]}）;
end
运行结果：
（2）所加窗函数；
程序：
M=[11,41,81,121];
fori=1:
1:
4
n=[0:
1:
M（i）-1];
wn_boxcar=boxcar（M（i））;
wn_bartlett=bartlett（M（i））;
wn_hanning=hanning（M（i））;
wn_hamming=hamming（M（i））;
wn_blackman=blackman（M（i））;
figure（i）
subplot（3,2,1）,stem（n,wn_boxcar,'r.'）;title（['M=',num2str（i）,'矩形窗函数']）;
subplot（3,2,2）,stem（n,wn_bartlett,'r.'）;title（['M=',num2str（i）,'三角窗函数']）;
subplot（3,2,3）,stem（n,wn_hanning,'r.'）;title（['M=',num2str（i）,'汉宁窗函数']）;
subplot（3,2,4）,stem（n,wn_hamming,'r.'）;title（['M=',num2str（i）,'海明窗函数']）;
subplot（3,2,5）,stem（n,wn_blackman,'r.'）;title（['M=',num2str（i）,'布莱克曼窗函数']）;
end
运行结果：
（3）加窗后的滤波器冲激响应
（4）滤波器的幅频特性；
程序：
clearall;
clc;
wc=0.3;
N=[11,41,81,121];
k=1;
fori=1:
1:
4
M=N（i）-1;
h_boxcar=fir1（M,wc,boxcar（N（i）））;
h_bartlett=fir1（M,wc,bartlett（N（i）））;
h_hanning=fir1（M,wc,hanning（N（i）））;
h_hamming=fir1（M,wc,hamming（N（i）））;
h_blackman=fir1（M,wc,blackman（N（i）））;
n=[0:
1:
N（i）-1];
figure（k）
k=k+1;
subplot（3,2,1）,stem（n,h_boxcar,'.'）;title（['N=',num2str（N（i））,'矩形窗函数冲击响应']）;
subplot（3,2,2）,stem（n,h_bartlett,'.'）;title（['N=',num2str（N（i））,'三角窗函数冲击响应']）;
subplot（3,2,3）,stem（n,h_hanning,'.'）;title（['N=',num2str（N（i））,'汉宁窗函数冲击响应']）;
subplot（3,2,4）,stem（n,h_hamming,'.'）;title（['N=',num2str（N（i））,'海明窗函数冲击响应']）;
subplot（3,2,5）,stem（n,h_blackman,'.'）;title（['N=',num2str（N（i））,'布莱克曼窗函数冲击响应']）;
Hx_boxcar=freqz（h_boxcar）;dBHx_boxcar=20*log10（（abs（Hx_boxcar）/max（abs（Hx_boxcar））））;
Hx_bartlett=freqz（h_bartlett）;dBHx_bartlett=20*log10（（abs（Hx_bartlett）/max（abs（Hx_bartlett））））;
Hx_hanning=freqz（h_hanning）;dBHx_hanning=20*log10（（abs（Hx_hanning）/max（abs（Hx_hanning））））;
Hx_hamming=freqz（h_hamming）;dBHx_hamming=20*log10（（abs（Hx_hamming）/max（abs（Hx_hamming））））;
Hx_blackman=freqz（h_blackman）;dBHx_blackman=20*log10（（abs（Hx_blackman）/max（abs（Hx_blackman））））;
w=0:
pi/512:
511*pi/512;
figure（k）
k=k+1;
subplot（3,2,1）,plot（w/pi,dBHx_boxcar）;title（['N=',num2str（N（i））,'矩形窗幅度响应（dB）']）;xlabel（'\omega/\pi'）;ylabel（'20log|H（e^j^\omega）|（dB）'）;gridon;
subplot（3,2,2）,plot（w/pi,dBHx_bartlett）;title（['N=',num2str（N（i））,'三角窗幅度响应（dB）']）;xlabel（'\omega/\pi'）;ylabel（'20log|H（e^j^\omega）|（dB）'）;gridon;
subplot（3,2,3）,plot（w/pi,dBHx_hanning）;title（['N=',num2str（N（i））,'汉宁窗幅度响应（dB）']）;xlabel（'\omega/\pi'）;ylabel（'20log|H（e^j^\omega）|（dB）'）;gridon;
subplot（3,2,4）,plot（w/pi,dBHx_hamming）;title（['N=',num2str（N（i））,'海明窗幅度响应（dB）']）;xlabel（'\omega/\pi'）;ylabel（'20log|H（e^j^\omega）|（dB）'）;gridon;
subplot（3,2,5）,plot（w/pi,dBHx_blackman）;title（['N=',num2str（N（i））,'布莱克曼幅度响应（dB）']）;xlabel（'\omega/\pi'）;ylabel（'20log|H（e^j^\omega）|（dB）'）;gridon;
end
运行结果：
（5）根据结果比较不同长度对应的滤波器特性；
截断长度N增加时，主瓣宽度变窄，过渡带减小，起伏振荡变密。
（6）比较不同的窗对应的效果。
随着窗函数越来越平滑，主瓣宽度加宽，通带及阻带波纹
改善，且阻带最小衰减增大，但过渡带加宽，可见窗函数
设计中，通带、阻带波纹的改善是以加宽过渡带为代价的。
用到的函数：
fir1
3.IIR滤波器的实现结构
（1）用级联形式实现该系统，画出实现框图。
（2）用并联形式实现该系统，画出实现框图。
用到的函数：
tf2zp，residue，residuez
程序：
b=[3,8,12,7,2,-2];
a=[16,24,24,14,5,1];
[sos,G]=tf2sos（b,a）
[r,p,k]=residuez（b,a）
R1=[r
（1）,r
（2）];
P1=[p
（1）,p
（2）];
[b1,a1]=residuez（R1,P1,0）
R2=[r（4）,r（5）];P2=[p（4）,p（5）];
[b2,a2]=residuez（R2,P2,0）
运行结果：
sos=
1.0000-0.333301.00000.50000
1.00002.00002.00001.00000.50000.2500
1.00001.00001.00001.00000.50000.5000
G=
0.1875
r=
-0.4219+0.6201i
-0.4219-0.6201i
2.3438
0.3437-2.5079i
0.3437+2.5079i
p=
-0.2500+0.6614i
-0.2500-0.6614i
-0.5000
-0.2500+0.4330i
-0.2500-0.4330i
k=
-2
b1=
-0.8437-1.03120
a1=
1.00000.50000.5000
b2=
0.68752.34370
a2=
1.00000.50000.2500
升级会员 