数字信号处理上机实验 滤波器.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