数字信号处理上机实验 滤波器.docx
《数字信号处理上机实验 滤波器.docx》由会员分享，可在线阅读，更多相关《数字信号处理上机实验 滤波器.docx（36页珍藏版）》请在冰豆网上搜索。
数字信号处理上机实验滤波器
数字信号处理上机练习
7.12
（1）双线性变换法
程序：
%双线性变换法设计切比雪夫1型带通滤波器
%w计算式：
w=（f/ADfs）*2*pi单位：
rad
clc;clear;closeall;
ADfs=2000;%ADC的采样频率，单位：
Hz
ws1=（100/ADfs）*2*pi;%下阻带截止频率，单位：
rad
wp1=（200/ADfs）*2*pi;%通带下截止频率，单位：
rad
wp2=（400/ADfs）*2*pi;%通带上截止频率，单位：
rad
ws2=（600/ADfs）*2*pi;%上阻带截止频率，单位：
rad
display（'数字滤波器的边界频率单位：
rad'）;ws1,wp1,wp2,ws2
Rp=2;%通带最大衰减，单位：
dB
As=20;%阻带最小衰减，单位：
dB
T=1;%模拟滤波器数字化时的采样间隔，单位：
s
Ws1=（2/T）*tan（ws1/2）;Wp1=（2/T）*tan（wp1/2）;
Wp2=（2/T）*tan（wp2/2）;Ws2=（2/T）*tan（ws2/2）;%转换为模拟滤波器指标
Ws=[Ws1,Ws2];Wp=[Wp1,Wp2];
display（'转换后滤波器的边界频率单位：
rad/s'）;Ws1,Wp1,Wp2,Ws2
[n,wn]=cheb1ord（Wp,Ws,Rp,As,'s'）;%计算切比雪夫1型滤波器阶数和截止频率
display（'滤波器阶数'）;n
display（'2dB截止频率'）;wn
[b,a]=cheby1（n,Rp,wn,'s'）;%计算模拟带通滤波器的传输函数Ha（s）
display（'模拟原型系统函数Ha（s）分子多项式系数：
'）;b
display（'模拟原型系统函数Ha（s）分母多项式系数'）;a
[Ha,Wa]=freqs（b,a,512）;%计算并画出模拟原型滤波器的幅度频谱图
subplot（211）;
plot（Wa/pi,20*log10（abs（Ha））,'LineWidth',2）;
axis（[0,max（Wa/pi）,min（20*log10（abs（Ha）））,5]）;
xlabel（'角频率rad/s（×π）'）;
ylabel（'对数幅度dB'）;
title（'模拟原型带通滤波器的幅度谱（对数幅度）'）;
grid;
[bz,az]=bilinear（b,a,1/T）;%用双线性变换法转换为数字滤波器
display（'数字滤波器的系统函数H（z）分子多项式系数：
'）;bz
display（'数字滤波器的系统函数H（z）分母多项式系数'）;az
[Hz,Wz]=freqz（bz,az,512,1/T）;%计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位
subplot（212）;
plot（Wz*2*T,20*log10（abs（Hz））,'LineWidth',2）;
axis（[0,1,-100,5]）;
xlabel（'数字域频率（×π）'）;
ylabel（'对数幅度dB'）;
title（'双线性变换法设计的数字滤波器的幅度谱（对数幅度）'）;
grid;
结果：
数字滤波器的边界频率单位：
rad
ws1=
0.3142
wp1=
0.6283
wp2=
1.2566
ws2=
1.8850
转换后滤波器的边界频率单位：
rad/s
Ws1=
0.3168
Wp1=
0.6498
Wp2=
1.4531
Ws2=
2.7528
滤波器阶数
n=
2
2dB截止频率
wn=
0.64981.4531
模拟原型系统函数Ha（s）分子多项式系数：
b=
000.421800
模拟原型系统函数Ha（s）分母多项式系数
a=
1.00000.64572.41960.60970.8916
数字滤波器的系统函数H（z）分子多项式系数：
bz=
0.05120.0000-0.10240.00000.0512
数字滤波器的系统函数H（z）分母多项式系数
az=
1.0000-2.07332.4881-1.59440.6125
（2）冲激响应不变法
程序：
%冲激响应不变法设计切比雪夫1型带通滤波器
%w计算式：
w=（f/ADfs）*2*pi单位：
rad
clc;clear;closeall;
ADfs=2000;%ADC的采样频率，单位：
Hz
ws1=（100/ADfs）*2*pi;%下阻带截止频率，单位：
rad
wp1=（200/ADfs）*2*pi;%通带下截止频率，单位：
rad
wp2=（400/ADfs）*2*pi;%通带上截止频率，单位：
rad
ws2=（600/ADfs）*2*pi;%上阻带截止频率，单位：
rad
display（'数字滤波器的边界频率单位：
rad'）;ws1,wp1,wp2,ws2
Rp=2;%通带最大衰减，单位：
dB
Rs=20;%阻带最小衰减，单位：
dB
T=1;%模拟滤波器数字化时的采样间隔，单位：
s
Ws1=ws1/T;Wp1=wp1/T;
Wp2=wp2/T;Ws2=ws2/T;%转换为模拟滤波器指标
display（'转换后滤波器的边界频率单位：
rad/s'）;Ws1,Wp1,Wp2,Ws2
Ws=[Ws1,Ws2];Wp=[Wp1,Wp2];
[n,wn]=cheb1ord（Wp,Ws,Rp,Rs,'s'）;%计算切比雪夫1型滤波器阶数和截止频率
display（'滤波器阶数'）;n
display（'2dBfrequency'）;wn
[b,a]=cheby1（n,Rp,wn,'s'）;%计算模拟带通滤波器的传输函数Ha（s）
display（'模拟原型系统函数Ha（s）分子多项式系数：
'）;b
display（'模拟原型系统函数Ha（s）分母多项式系数'）;a
[Ha,Wa]=freqs（b,a,512）;%计算并画出模拟原型滤波器的幅度频谱图
subplot（211）;
plot（Wa/pi,20*log10（abs（Ha））,'LineWidth',2）;
axis（[0,1.5,-80,5]）;
xlabel（'角频率rad/s（×π）'）;
ylabel（'对数幅度dB'）;
title（'模拟原型带通滤波器的幅度谱（对数幅度）'）;
grid;
[bz,az]=impinvar（b,a,1/T）;%用冲激响应不变法转换为数字滤波器
display（'数字滤波器的系统函数H（z）分子多项式系数：
'）;bz
display（'数字滤波器的系统函数H（z）分母多项式系数'）;az
[Hz,Wz]=freqz（bz,az,512,1/T）;%计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位
subplot（212）;
plot（Wz*2*T,20*log10（abs（Hz））,'LineWidth',2）;
xlabel（'数字域频率（×π）'）;
ylabel（'对数幅度dB'）;
title（'冲激响应不变法设计的数字滤波器的幅度谱（对数幅度）'）;
axis（[0,1,-80,5]）
grid;
结果：
数字滤波器的边界频率单位：
rad
ws1=
0.3142
wp1=
0.6283
wp2=
1.2566
ws2=
1.8850
转换后滤波器的边界频率单位：
rad/s
Ws1=
0.3142
Wp1=
0.6283
Wp2=
1.2566
Ws2=
1.8850
滤波器阶数
n=
3
2dBfrequency
wn=
0.62831.2566
模拟原型系统函数Ha（s）分子多项式系数：
b=
0000.0811000
模拟原型系统函数Ha（s）分母多项式系数
a=
1.00000.46362.77220.81322.18890.28900.4922
数字滤波器的系统函数H（z）分子多项式系数：
bz=
0.00000.0272-0.05810.01090.0437-0.02370
数字滤波器的系统函数H（z）分母多项式系数
az=
1.0000-3.30306.0060-6.74635.1356-2.40930.6290
7.14
程序：
%双线性变换法设计切比雪夫1型带通滤波器
%w计算式：
w=（f/ADfs）*2*pi单位：
rad
clc;clear;closeall;
ADfs=10000;%ADC的采样频率，单位：
Hz
wp1=（500/ADfs）*2*pi;%下通带截止频率，单位：
rad
ws1=（1000/ADfs）*2*pi;%阻带下截止频率，单位：
rad
ws2=（2000/ADfs）*2*pi;%阻带上截止频率，单位：
rad
wp2=（3000/ADfs）*2*pi;%上通带截止频率，单位：
rad
display（'数字滤波器的边界频率单位：
rad'）;wp1,ws1,ws2,wp2
Rp=3;%通带最大衰减，单位：
dB
As=30;%阻带最小衰减，单位：
dB
T=1;%模拟滤波器数字化时的采样间隔，单位：
s
Ws1=（2/T）*tan（ws1/2）;Wp1=（2/T）*tan（wp1/2）;
Wp2=（2/T）*tan（wp2/2）;Ws2=（2/T）*tan（ws2/2）;%转换为模拟滤波器指标
display（'转换后滤波器的边界频率单位：
rad/s'）;Wp1,Ws1,Ws2,Wp2
Ws=[Ws1,Ws2];Wp=[Wp1,Wp2];
[n,wn]=cheb1ord（Wp,Ws,Rp,As,'s'）;%计算切比雪夫1型滤波器阶数和截止频率
display（'滤波器阶数'）;n
display（'3dB截止频率'）;wn
[b,a]=cheby1（n,Rp,wn,'stop','s'）;%计算模拟带阻滤波器的传输函数Ha（s）
display（'模拟原型系统函数Ha（s）分子多项式系数：
'）;b
display（'模拟原型系统函数Ha（s）分母多项式系数'）;a
[Ha,Wa]=freqs（b,a,512）;%计算并画出模拟原型滤波器的幅度频谱图
subplot（211）;
plot（Wa/pi,20*log10（abs（Ha））,'LineWidth',2）;
axis（[0,3,-120,5]）;
xlabel（'角频率rad/s（×π）'）;
ylabel（'对数幅度dB'）;
title（'模拟原型带阻滤波器的幅度谱（对数幅度）'）;
grid;
[bz,az]=bilinear（b,a,1/T）;%用双线性变换法转换为数字滤波器
display（'数字滤波器的系统函H（z）分子多项式系数：
'）;bz
display（'数字滤波器的系统函H（z）分母多项式系数'）;az
[Hz,Wz]=freqz（bz,az,512,1/T）;%计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位
subplot（212）;
plot（Wz*2*T,20*log10（abs（Hz））,'LineWidth',2）;
xlabel（'数字域频率（×π）'）;
ylabel（'对数幅度dB'）;
title（'双线性变换法设计的数字滤波器的幅度谱（对数幅度）'）;
axis（[0,1,-120,5]）
grid;
结果：
数字滤波器的边界频率单位：
rad
wp1=
0.3142
ws1=
0.6283
ws2=
1.2566
wp2=
1.8850
转换后滤波器的边界频率单位：
rad/s
Wp1=
0.3168
Ws1=
0.6498
Ws2=
1.4531
Wp2=
2.7528
滤波器阶数
n=
3
3dB截止频率
wn=
0.34302.7527
模拟原型系统函数Ha（s）分子多项式系数：
b=
1.000002.832602.674500.8418
模拟原型系统函数Ha（s）分母多项式系数
a=
1.00008.927016.671672.694415.74137.95840.8418
数字滤波器的系统函H（z）分子多项式系数：
bz=
0.0946-0.35080.7174-0.88020.7174-0.35080.0946
数字滤波器的系统函H（z）分母多项式系数
az=
1.0000-1.46010.3179-0.35070.68850.2289-0.3824
7.15
程序：
%双线性变换法设计切比雪夫2型带阻滤波器
%w计算式：
w=（f/ADfs）*2*pi单位：
rad
clc;clear;closeall;
ADfs=100000;%ADC的采样频率，单位：
Hz
wp1=（5000/ADfs）*2*pi;%下通带截止频率，单位：
rad
ws1=（10000/ADfs）*2*pi;%阻带下截止频率，单位：
rad
ws2=（20000/ADfs）*2*pi;%阻带上截止频率，单位：
rad
wp2=（30000/ADfs）*2*pi;%上通带截止频率，单位：
rad
display（'数字滤波器的边界频率单位：
rad'）;wp1,ws1,ws2,wp2
Rp=3;%通带最大衰减，单位：
dB
As=30;%阻带最小衰减，单位：
dB
T=1;%模拟滤波器数字化时的采样间隔，单位：
s
Ws1=（2/T）*tan（ws1/2）;Wp1=（2/T）*tan（wp1/2）;
Wp2=（2/T）*tan（wp2/2）;Ws2=（2/T）*tan（ws2/2）;%转换为模拟滤波器指标
display（'转换后模拟原型滤波器的边界频率单位：
rad/s'）;Wp1,Ws1,Ws2,Wp2
Ws=[Ws1,Ws2];Wp=[Wp1,Wp2];
[n,wn]=cheb2ord（Wp,Ws,Rp,As,'s'）;%计算切比雪夫2型滤波器阶数和截止频率
display（'滤波器阶数'）;n
display（'3dB截止频率'）;wn
[b,a]=cheby2（n,As,wn,'stop','s'）;%计算模拟带阻滤波器的传输函数Ha（s）
display（'模拟原型滤波器系统函数Ha（s）分子多项式系数：
'）;b
display（'模拟原型滤波器系统函数Ha（s）分母多项式系数'）;a
[Ha,Wa]=freqs（b,a,512）;%计算并画出模拟原型滤波器的幅度频谱图
subplot（211）;
plot（Wa/pi,20*log10（abs（Ha））,'LineWidth',2）;
axis（[0,1.5,min（20*log10（abs（Ha）））,5]）;
xlabel（'角频率rad/s（×π）'）;
ylabel（'对数幅度dB'）;
title（'模拟原型带阻滤波器的幅度谱（对数幅度）'）;
grid;
[bz,az]=bilinear（b,a,1/T）;%用双线性变换法转换为数字滤波器
display（'数字滤波器系统函数H（z）分子多项式系数：
'）;bz
display（'数字滤波器系统函数H（z）分母多项式系数'）;az
[Hz,Wz]=freqz（bz,az,512,1/T）;%计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位
subplot（212）;
plot（Wz*2*T,20*log10（abs（Hz））,'LineWidth',2）;
xlabel（'数字域频率（×π）'）;
ylabel（'对数幅度dB'）;
title（'双线性变换法设计的数字滤波器的幅度谱（对数幅度）'）;
grid;
结果：
数字滤波器的边界频率单位：
rad
wp1=
0.3142
ws1=
0.6283
ws2=
1.2566
wp2=
1.8850
转换后模拟原型滤波器的边界频率单位：
rad/s
Wp1=
0.3168
Ws1=
0.6498
Ws2=
1.4531
Wp2=
2.7528
滤波器阶数
n=
3
3dB截止频率
wn=
0.55721.6946
模拟原型滤波器系统函数Ha（s）分子多项式系数：
b=
1.0000-0.00003.8028-0.00003.5906-0.00000.8418
模拟原型滤波器系统函数Ha（s）分母多项式系数
a=
1.00004.245812.816119.644412.10093.78510.8418
数字滤波器系统函数H（z）分子多项式系数：
bz=
0.2263-0.76251.4500-1.74061.4500-0.76250.2263
数字滤波器系统函数H（z）分母多项式系数
az=
1.0000-1.94771.5590-1.02850.7650-0.28940.0286
7.17
（1）椭圆函数滤波器
程序：
%设计椭圆函数带通数字滤波器
%w计算式：
w/pi=2f/ADfs
clc;clear;closeall;
ADfs=25000;%ADC的采样频率，单位：
Hz
ws1=（3500/ADfs）*2;%下阻带截止频率，单位：
*pirad
wp1=（5000/ADfs）*2;%通带下截止频率，单位：
*pirad
wp2=（7000/ADfs）*2;%通带上截止频率，单位：
*pirad
ws2=（8500/ADfs）*2;%上阻带截止频率，单位：
*pirad
Rp=0.5;%通带最大衰减，单位：
dB
As=45;%阻带最小衰减，单位：
dB
ws=[ws1,ws2];wp=[wp1,wp2];
[n,wc]=ellipord（wp,ws,Rp,As）;%计算椭圆函数带通数字滤波器阶数和截止频率
display（'数字滤波器边界频率单位：
*pirad'）;ws1,wp1,wp2,ws2
display（'滤波器阶数'）;n
display（'0.5dB截止频率单位：
*pirad'）;wc
[b,a]=ellip（n,Rp,As,wc）;%计算椭圆函数带通数字滤波器的传输函数H（z）
display（'H（z）分子多项式系数：
'）;b
display（'H（z）分母多项式系数'）;a
hn=impz（b,a）;%计算单位冲激响应
[H,w]=freqz（b,a）;
%画出系统的幅频特性曲线
subplot（3,1,1）;
plot（12.5*w/pi,20*log10（abs（H））,'LineWidth',2）;%
xlabel（'f（khz）'）;
ylabel（'dB'）;
title（'椭圆函数带通数字滤波器的幅度响应（db）'）;
axis（[0,12.5,-70,3]）;
grid;
%画出系统的相频特性曲线
subplot（3,1,2）;
plot（w/pi,angle（H）,'LineWidth',2）;%
xlabel（'w/π'）;
ylabel（'rad'）;
title（'椭圆函数带通数字滤波器的相位响应'）;
grid;
%画出单位冲激响应
subplot（3,1,3）;
stem（hn,'.','Linewidth',2）;
title（'系统单位冲激响应'）;
xlabel（'n'）;
axis（[0,120,-0.2,0.2]）;
结果：
数字滤波器边界频率单位：
*pirad
ws1=
0.2800
wp1=
0.4000
wp2=
0.5600
ws2=
0.6800
滤波器阶数
n=
4
0.5dB截止频率单位：
*pirad
wc=
0.40000.5600
H（z）分子多项式系数：
b=
0.0113-0.00360.0110-0.00450.0196-0.00450.0110-0.00360.0113
H（z）分母多项式系数
a=
1.0000-0.46493.2535-1.14434.1557-0.98892.4381-0.29740.5525
（2）切比雪夫Ⅰ型滤波器
程序：
%设计切比雪夫1型带通数字滤波器
%w计算式：
w/pi=2f/ADfs
clc;clear;closeall;
ADfs=25000;%ADC的采样频率，单位：
Hz
ws1=（3500/ADfs）*2;%下阻带截止频率，单位：
*pirad
wp1=（5000/ADfs）*2;%通带下截止频率，单位：
*pirad
wp2=（7000/ADfs）*2;%通带上截止频率，单位：
*pirad
ws2=（8500/ADfs）*2;%上阻带截止频率，单位：
*pirad
Rp=0.5;%通带最大衰减，单位：
dB
As=45;%阻带最小衰减，单位：
dB
ws=[ws1,ws2];wp=[wp1,wp2];
[n,wc]=cheb1ord（wp,ws,Rp,As）;%计算切比雪夫1型带通数字滤波器阶数和截止频率
display（'数字滤波器边界频率单位：
*pirad'）;ws1,wp1,wp2,ws2
display（'滤波器阶数'）;n
display（'0.5dB截止频率单位：
*pirad'）;wc
[b,a]=cheby1（n,Rp,wc）;%计算切比雪夫1型带通数字滤波器的传输函数H（z）
display（'H（z）分子多项式系数：
'）;b
display（'H（z）分母多项式系数'）;a
hn=impz（b,a）;%计算单位冲激响应
[H,w]=freqz（b,a）;
%画出系统的幅频特性曲线
subplot（3,1,1）;
plot（12.5*w/pi,20*log10（abs（H））,'LineWidth',2）;%
xlabel（'f（khz）'）;
ylabel（'dB'）;
title（'切比雪夫1型带通数字滤波器的幅度响应（db）'）;
axis（[0,12.5,-70,3]）;
grid;
%画出系统的相频特性曲线
subplot（3,1,2）;
plot（w/pi,angle（H）,'LineWidth',2）;%
xlabel（'w/π'）;
ylabel（'rad'）;
title（'切比雪夫1型带通数字滤波器的相位响应'）;
grid;
%画出单位冲激响应
subplot（3,1,3）;
stem（hn,'.','Linewidth',2）;
title（'系统单位冲激响应'）;
xlabel（'n'）;
axis（[0,120,-0.2
升级会员 