数字信号处理实验七.docx
《数字信号处理实验七.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验七.docx（21页珍藏版）》请在冰豆网上搜索。
数字信号处理实验七
实验七数字滤波器设计
实验室名称：
信息学院2204实验时间：
2015年11月26日
姓名：
蒋逸恒学号：
20131120038专业：
通信工程指导教师：
陶大鹏
成绩
教师签名：
年月日
一、实验目的
1、学会使用matlab设计符合要求的数字滤波器的参数。
2、学会使用基本的matlab命令，把获得的数字滤波器参数带到传输函数中写出传输函数，之后利用之前的实验知识，画出滤波器级联框图，实现滤波器。
2、实验内容
Q7.1用MATLAB确定一个数字无限冲击响应低通滤波器所有四种类型的最低阶数。
指标如下：
40kHz的抽样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。
评论你的结果。
Q7.2用MATLAB确定一个数字无限冲击响应高通滤波器所有四种类型的最低阶数。
指标如下：
3500Hz的抽样率，1050Hz的通带边界频率，600Hz的阻带边界频率，1dB的通带波纹，50dB的最小阻带衰减。
评论你的结果。
Q7.5通过运行程序P7.1来设计巴特沃兹带阻滤波器。
写出所产生的传输函数的准确表达式。
滤波器的指标是什么？
你的设计符合指标吗？
使用MATLAB，计算并绘制滤波器的未畸变的相位响应及群延迟响应。
Q7.6修改程序P7.1来设计符合习题Q7.1所给指标的切比雪夫1型低通滤波器。
写出所产生的传输函数的准确表达式。
你的设计符合指标吗？
使用MATLAB，计算并绘制滤波器的未畸变的相位响应及群延迟响应。
Q7.20使用函数fir1，设计一个线性相位有限冲击响应低通滤波器，使其满足习题Q7.13给出的指标，并画出其增益和相位响应。
使用习题Q7.13中用凯泽公式估计出的阶数。
用表格形式显示滤波器系数。
你的设计满足指标吗？
若不满足，调整滤波器阶数直到设计满足指标。
满足指标的滤波器阶数是多少？
Q7.23用凯泽窗设计一个有限冲击响应低通滤波器。
滤波器的指标是：
wp=0.31,ws=0.41,As=50dB。
注意，函数kaiser需要参数β及阶数N的值，他们必须先用式（7.36）和式（7.37）分别算出。
你的设计满足指标吗？
Q7.25用fir2设计一个95阶有限冲击响应滤波器，它具有三个不同的常熟幅度级：
在频率范围0到0.25中为0.4，在频率范围0.3到0.45中为1.0，在频率范围0.5到1.0中为0.8。
画出所设计的滤波器的幅度响应。
你的设计满足指标吗？
Q7.27用remez设计具有如下指标的有限冲击响应带通滤波器：
通带边界为1.8kHz和3.0kHz，阻带边界为1.5kHz和4.2kHz，通带波纹δp=0.1，阻带波纹δs=0.02，抽样率为12kHz。
用kaiserord估计滤波器的阶数。
你的设计是一个最优有限冲击响应滤波器吗？
你的设计满足指标吗？
若不满足，增加滤波器阶数在满足指标方面有用吗？
指标由一个较低阶数的滤波器来满足而不是由kaiserord得到的来满足吗？
在不等过渡带的情形下，用设计的滤波器可能在较大的过度带宽中以增益响应表现不满意的行为。
改进该行为的一种方法是：
通过移动阻带边界减少过渡带宽，直到使设计在过渡带中以平滑的下降来满足指标。
在通带边界保持固定的情况下，尝试这种方法并确定新的指标，它在过渡带中提供平滑的下降。
三、实验器材及软件
1.微型计算机1台
2.MATLAB7.0软件
4、实验原理
IIR数字滤波器采用递归型结构，即结构上带有反馈环路。
IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。
由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。
相对的，FIR滤波器是直接型结构，没有递归结构，也就没有反馈环路，FIR滤波器运算结构通常也是由由延时、乘以系数和相加等基本运算组成，可以组合成各种结构形式的滤波器组来实现滤波。
IIR滤波器的单位脉冲响应为无限长，网络中有反馈回路。
FIR（FiniteImpulseResponse）滤波器的单位脉冲响应是有限长的，一般网络中没有反馈回路。
IIR滤波器的系统函数一般是一个有理分式，分母多项式决定滤波器的反馈网络。
FIR滤波器的系统函数一般是一个分母为1的多项式，所以没有反馈网络。
IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。
在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
FIR的设计在matlab中其实也同样很简单，其主要思想是通过先设计出无限长的IIR滤波器再通过截尾的方式获得一个FIR滤波器，但是简单的截尾会出现吉布斯现象，所以此时还需要对截尾之后的结果加窗；来减小吉布斯震荡现象。
IIR数字滤波器幅频特性精度很高，不是线性相位的，可以应用于对相位信息不敏感的音频信号上；FIR数字滤波器的幅频特性精度较之于IIR数字滤波器低，但是线性相位，就是不同频率分量的信号经过FIR滤波器后他们的时间差不变，这是很好的性质。
FIR数字滤波器是有限的单位响应也有利于对数字信号的处理，便于编程。
五、实验步骤
1、进行本实验，首先必须熟悉matlab的运用，所以第一步是学会使用matlab。
2、学习相关基础知识，根据《数字信号处理》课程的学习理解实验内容和目的。
3、在充分熟悉基础知识的情况下进行实验，利用matlab完成各种简单的波形产生和观察，理解各种波形产生的原理和方法。
4、从产生的图形中学习新的知识，掌握实验的目的，充分学习数字滤波器的设计和运用。
5、最后需要思考各种滤波器的联系和区别，以及如何根据要求设计出最符合要求且最容易实现的滤波器
6、小结本实验，通过小结发现问题并解决问题。
6、实验记录（数据、图表、波形、程序等）
7.1
输入：
[N1,wn1]=buttord（0.2,0.4,0.5,40）
[N2,wn2]=cheb1ord（0.2,0.4,0.5,40）
[N3,wn3]=cheb2ord（0.2,0.4,0.5,40）
[N4,wn4]=ellipord（0.2,0.4,0.5,40）
输出：
N1=8
wn1=0.2469
N2=5
wn2=0.2000
N3=5
wn3=0.4000
N4=4
wn4=0.2000
7.5
输入：
Ws=[0.40.6];
Wp=[0.20.8];
Rp=0.4;Rs=50;
[N1,Wn1]=buttord（Wp,Ws,Rp,Rs）;
[num,den]=butter（N1,Wn1,'stop'）;
disp（'分子系数：
'）;disp（num）;
disp（'分母系数：
'）;disp（den）;
[g,w]=gain（num,den）;
plot（w/pi,g）;grid
axis（[01-605]）;
xlabel（'\omega/\pi'）;ylabel（'增益，dB'）;
title（'巴特沃斯带阻滤波器的增益响应'）;
输出：
分子系数：
0.04930.00000.24650.00000.49300.00000.4930
0.00000.24650.00000.0493
分母系数：
1.00000.0000-0.08500.00000.63600.0000-0.0288
0.00000.05610.0000-0.0008
相位响应和群延迟响应：
w=0:
pi/255:
pi;
Ws=[0.40.6];Wp=[0.20.8];Rp=0.4;Rs=50;
[N1,Wn1]=buttord（Wp,Ws,Rp,Rs）;
[num,den]=butter（N1,Wn1,'stop'）;
h=freqz（num,den,w）;
subplot（1,2,1）;
plot（w/pi,unwrap（angle（h）））;
gridon;
k=[diff（angle（h））./diff（w）0];
subplot（1,2,2）;
plot（w/pi,k）;
axis（[0,1,3,10]）;
gridon;
7.6
输入：
Ws=0.4;Wp=0.2;Rp=0.5;Rs=40;
[N1,Wn1]=buttord（Wp,Ws,Rp,Rs）;
[num,den]=butter（N1,Wn1）;
disp（'分子系数：
'）;disp（num）;
disp（'分母系数：
'）;disp（den）;
[g,w]=gain（num,den）;
plot（w/pi,g）;grid
axis（[01-605]）;
xlabel（'\omega/\pi'）;
ylabel（'增益，dB'）;
title（'切比雪夫1型低通滤波器的增益响应'）;
输出：
分子系数：
0.00040.00200.00400.00400.00200.0004
分母系数：
1.0000-3.82696.2742-5.44642.4915-0.4797
相位响应和群延迟响应：
Ws=0.4;Wp=0.2;Rp=0.5;Rs=40;
%估计滤波器阶数
[N1,Wn1]=cheb1ord（Wp,Ws,Rp,Rs）;
%设计滤波器
[num,den]=cheby1（N1,Rp,Wn1）;
h=freqz（num,den,w）;
subplot（1,2,1）;
plot（w/pi,unwrap（angle（h）））;
gridon;
k=[diff（angle（h））./diff（w）0];
subplot（1,2,2）;
plot（w/pi,k）;
axis（[0,1,-5,20]）;
gridon;
7.20输入：
w=0:
pi/255:
pi;
Fs=2500;Fp=2000;
Rp=0.005;Rs=0.005;
Ft=10000;
N=kaiord（Fp,Fs,Rp,Rs,Ft）;
b=fir1（N,2*Fp/Ft）;
[g,w]=gain（b,1）;
subplot（2,1,1）;plot（w/pi,g）;grid;
axis（[01-1255]）
xlabel（'\omega/\pi'）;
ylabel（'增益，dB'）;
title（'FIR低通滤波器的增益响应'）;
h=freqz（b,1,w）;
subplot（2,1,2）;
plot（w/pi,unwrap（angle（h）））;gridon;
xlabel（'\omega/\pi'）;ylabel（'相位'）;
title（'FIR低通滤波器的相频响应'）;
输出：
b=
-0.00070.00070.0014-0.0000-0.0023-0.00190.0025
0.0052-0.0000-0.0083-0.00640.00790.0157-0.0000
-0.0233-0.01760.02150.0431-0.0000-0.0706-0.0600
0.09190.30130.39980.30130.0919-0.0600-0.0706
-0.00000.04310.0215-0.0176-0.0233-0.00000.0157
0.0079-0.0064-0.0083-0.00000.00520.0025-0.0019
-0.0023-0.00000.00140.0007-0.0007
7.23
输入：
w2=0:
pi/255:
pi;
Wp=0.31;Ws=0.41;
Wn=Wp+（Ws-Wp）/2;
As=50;Rs=10^（-As/20）;
Rp=Rs;
ifAs>21
N=round（（As-7.95）*2/（14.36*（abs（Wp-Ws）））+1）;
else
N=round（0.9222*2/abs（Wp-Ws）+1）;
end
ifAs>50
B=0.1102*（As-8.7）;
elseifAs>=21
B=0.5842*（As-21）^0.4+0.07886*（As-21）;
else
B=0;
end
s=kaiser（N+1,B）;
h=fir1（N,Wn,s）;
disp（'系数：
'）;
disp（h）;
[g,w]=gain（h,1）;
subplot（2,1,1）;
plot（w/pi,g）;
grid;
xlabel（'\omega/\pi'）;
ylabel（'增益，dB'）;
title（'增益响应'）;
Hz=freqz（h,1,w2）;
subplot（2,1,2）;
plot（w2/pi,unwrap（angle（Hz）））;
grid;
xlabel（'\omega/\pi'）;ylabel（'以弧度为单位的相位'）;
title（'解缠绕的相位响应'）;
输出：
系数：
0.00030.00080.0003-0.0011-0.00170.00000.0026
0.0027-0.0010-0.0049-0.00350.00330.00800.0034
-0.0074-0.0119-0.00180.01400.0161-0.0027-0.0241
-0.02010.01270.04060.0236-0.0354-0.0754-0.0258
0.12140.28710.35970.28710.1214-0.0258-0.0754
-0.03540.02360.04060.0127-0.0201-0.0241-0.0027
0.01610.0140-0.0018-0.0119-0.00740.00340.0080
0.0033-0.0035-0.0049-0.00100.00270.00260.0000
-0.0017-0.00110.00030.00080.0003
7.25
输入：
w2=0:
pi/255:
pi;
a=[0.40.41.01.00.80.8];
w=[00.250.30.450.51.0];
b=fir2（95,w,a）;
disp（'系数：
'）;
disp（b）;
[g,w]=gain（b,1）;
subplot（2,1,1）;
plot（w/pi,g）;
grid;
xlabel（'\omega/\pi'）;
ylabel（'增益，dB'）;
title（'增益响应'）;
h=freqz（b,1,w2）;
subplot（2,1,2）;
plot（w2/pi,unwrap（angle（h）））;grid;
xlabel（'\omega/\pi'）;ylabel（'以弧度为单位的相位'）;
title（'解缠绕的相位响应'）;
输出：
-0.00000.00000.00000.0000-0.0000-0.00000.0000
0.0000-0.0000-0.0000-0.0001-0.00010.00010.0003
0.0003-0.0001-0.0005-0.0005-0.0004-0.00000.0008
0.00160.0010-0.0013-0.0030-0.00190.00080.0025
0.00290.00280.0004-0.0053-0.0089-0.00340.0076
0.01210.0060-0.0024-0.0078-0.0136-0.0161-0.0003
0.03160.04440.0084-0.0537-0.0885-0.08100.7307
-0.0810-0.0885-0.05370.00840.04440.0316-0.0003
-0.0161-0.0136-0.0078-0.00240.00600.01210.0076
-0.0034-0.0089-0.00530.00040.00280.00290.0025
0.0008-0.0019-0.0030-0.00130.00100.00160.0008
-0.0000-0.0004-0.0005-0.0005-0.00010.00030.0003
0.0001-0.0001-0.0001-0.0000-0.00000.00000.0000
-0.0000-0.00000.00000.00000.0000-0.0000
7.27kaiserord获得阶数计算出的结果（73阶）
输入：
w2=0:
pi/255:
pi;
Fs1=1500;
Fp1=1800;
Fp2=3000;Fs2=4200;
Ft=12000;
ws1=2*Fs1/Ft;
wp1=2*Fp1/Ft;
wp2=2*Fp2/Ft;
ws2=2*Fs2/Ft;
Dp=0.1;Ds=0.02;
F=[Fs1Fp1Fp2Fs2];比kaiserord获得的阶数更低的
A=[010];阶数（20阶）得到的输出结果：
DEV=[DsDpDs];
[N,Wn,B,Ftap]=kaiserord（F,A,DEV,Ft）;
F2=[0ws1wp1wp2ws21];
A2=[001100];
h=remez（N,F2,A2）;
disp（'系数：
'）;
disp（h）;
[g,w]=gain（h,1）;
subplot（2,1,1）;
plot（w/pi,g）;
grid;
axis（[01-8060]）;
xlabel（'\omega/\pi'）;移动阻带边界（将4.2KHz移动到3.4KHz）来减少
ylabel（'增益，dB'）;过渡带宽后kaiserord获得阶数后计算的结果（73阶）
title（'增益响应'）;
Hz=freqz（h,1,w2）;
subplot（2,1,2）;
plot（w2/pi,unwrap（angle（Hz）））;
grid;
xlabel（'\omega/\pi'）;
ylabel（'以弧度为单位的相位'）;
title（'解缠绕的相位响应'）;
输出：
不使用kaiserord来获得阶数（73阶），而是用比kaiserord获得的阶数更低的一个阶数（20阶）得到的输出结果（图像记录是右边从上至下的第二个）：
系数：
0.0058-0.1330-0.02100.0709-0.00150.01460.1225
-0.0810-0.29530.04670.37890.0467-0.2953-0.0810
0.12250.0146-0.00150.0709-0.0210-0.13300.0058
移动阻带边界（将4.2KHz移动到3.4KHz）来减少过渡带宽，获得的输出，此时，阶数仍然是73阶（图像记录是右边从上至下的第三个）：
系数：
-0.00470.00660.00760.0007-0.0033-0.0002-0.0015
-0.00640.00030.01200.0069-0.0071-0.00610.0002
-0.0068-0.00730.01400.0219-0.0034-0.0194-0.0047
-0.0011-0.01450.00560.04270.0205-0.0368-0.0337
0.0046-0.0102-0.01760.07050.1129-0.0493-0.2139
-0.07750.20490.2049-0.0775-0.2139-0.04930.1129
0.0705-0.0176-0.01020.0046-0.0337-0.03680.0205
0.04270.0056-0.0145-0.0011-0.0047-0.0194-0.0034
0.02190.0140-0.0073-0.00680.0002-0.0061-0.0071
0.00690.01200.0003-0.0064-0.0015-0.0002-0.0033
0.00070.00760.0066-0.0047
7、实验思考题及解答
7.1归一化频率
，这里角频率取0.2；归一化频率
，这里角频率取0.4；理想通带波纹Rp=0.5dB；最小阻带衰减Rs=40dB。
巴特沃斯低通滤波器最低阶数N=8，截止频率Wn=0.2469
切比雪夫1型低通滤波器最低阶数N=5，截止频率Wn=0.2
切比雪夫2型低通滤波器最低阶数N=5，截止频率Wn=0.4
椭圆滤波器低通滤波器最低阶数N=4，截止频率Wn=0.2
从以上结果中观察到椭圆滤波器的阶数最低，阶数为4，截止频率也最低，为0.2，并且符合要求；此时的切比雪夫1型低通滤波器的截止频率也是最小的0.2；而巴特沃斯滤波器的阶数最高，实现最困难；截止频率最大的是切比雪夫2型，为0.4。
7.2归一化频率
，这里角频率取0.6；归一化频率
，这里角频率取0.34；理想通带波纹Rp=1dB；最小阻带衰减Rs=50dB。
巴特沃斯低通滤波器最低阶数N=8，截止频率Wn=0.5615
切比雪夫1型低通滤波器最低阶数N=5，截止频率Wn=0.6
切比雪夫2型低通滤波器最低阶数N=5，截止频率Wn=0.34
椭圆滤波器低通滤波器最低阶数N=4，截止频率Wn=0.6
从以上结果中观察到椭圆滤波器的阶数最低，阶数为4，并且符合要求；此时的切比雪夫2型低通滤波器的截止频率是最小的0.34；而巴特沃斯滤波器的阶数最高，实现最困难；截止频率最大的是切比雪夫1型和椭圆滤波器，为0.6。
7.5传输函数的表达式：
滤波器参数是：
通带归一化截至频率Wp1=0.2π，Wp2=0.8π，阻带归一化频率Ws1=0.4π，Ws2=0.6π，通带波纹Rp=0.4dB，最
升级会员 