IIR高通、带通及带阻数字滤波器设计.doc

西南科技大学
课 程 设 计 报 告

课程名称： 数字通信课程设计
设计名称：IIR高通、带通和带阻数字滤波器设计
姓 名： 张 威 威
学 号: 20074864
班 级： 通 信 0701
指导教师： 龙 惠 民
起止日期： 2010.6.28-2010.7.9

西南科技大学信息工程学院制
课 程 设 计 任 务 书
学生班级： 通信0701 学生姓名： 张威威 学号： 20074864
设计名称： IIR高通、带通和带阻数字滤波器设计
起止日期： 2010.6.28—2010.7.9 指导教师： 龙惠民

设计要求：

用MATLAB设计以下的IIR数字滤波器：
1、巴特沃思数字高通滤波器设计：
抽样频率为10kHZ,,通带截止频率为2.5 kHZ，通带衰减不大于2dB，阻带上限截止频率1.5kHZ ,阻带衰减不小于15 dB
2、巴特沃思数字带通滤波器设计：
抽样频率为10kHZ,,通带范围是1.5 kHZ到2.5 kHZ，通带衰减不大于3dB，在1kHZ和4kHZ处衰减不小于20 dB
3、巴特沃思数字带阻滤波器设计：
抽样频率为10kHZ,,在-2 dB衰减处的边带频率是1.5 kHZ，4 kHZ， 在-13dB衰减处频率是2kHZ和3kHZ
4、分别绘制这三种数字滤波器的幅度响应曲线和相位响应曲线；
5、采用切比雪夫Ⅰ型滤波器为原型重新设计上述三种数字滤波器；
6、分别绘制这三种数字滤波器的幅度响应曲线和相位响应曲线；
7、对两种滤波器原型的设计结果进行比较。

课 程 设 计 学 生 日 志

时间
设计内容
2009.06.28
查阅资料，确定方案
2009.06.30
设计总体方案
2009.07.03
编写程序
2009.07.06
编写程序
2009.07.08
撰写报告
2009.07.09
撰写报告
2009.07.10
答辩

课 程 设 计 考 勤 表

周
星期一
星期二
星期三
星期四
星期五

课 程 设 计 评 语 表

指导教师评语：

成绩： 指导教师：

年 月 日

IIR高通、带通和带阻数字滤波器设计
一、 设计目的和意义
数字滤波是数字信号分析中最重要的组成部分之一，数字滤波与模拟滤波相比，具有精度和稳定性高、系统函数容易改变、灵活性高、不存在阻抗匹配问题、便于大规模集成、可实现多维滤波等优点。
  本次主要设计高通、带通和带阻数字滤波器，比较利用巴特沃兹与切比雪夫1型设计的优缺点。并在实际应用中比较利弊选择使用。
二、 设计原理
IIR数字滤波器又叫做无限冲击响应数字滤波器。一般的做法是：先将给定的数字滤波器的技术指标转化为对应的模拟滤波器的技术指标，然后按照此指示设计模拟滤波器，再按照一定的规则将模拟滤波器离散化，这样就得到了数字滤波器。
设计的基本思路是：将低通转换为高通、带通和带阻，低通滤波器的映射函数为：

1、 低通——高通（LP-HP）
基本思想是：上述LP变换中的Z代以-Z，则LP=>HP。

2、 低通——带通
基本思想：把带通的中心频率

由以上分析得变换关系：

3、 低通——带阻
把带阻的中心频率

由此得到变换关系：

三、 详细设计步骤
1、巴特沃思数字高通滤波器设计：
抽样频率为10kHZ,,通带截止频率为2.5 kHZ，通带衰减不大于2dB，阻带上限截止频率1.5kHZ ,阻带衰减不小于15 dB。
设计步骤：
（1）所给的技术指标： fp=2500; fs=1500;Fs=10000; Rp=2;Rs=15
（2）归一化处理：wp=fp*2/Fs; ws=fs*2/Fs;
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=buttord(wp,ws,Rp,Rs);
（4）设计巴特沃兹高通滤波器，调用系统函数：[B,A]=butter(N,Wn,'high');
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。
2、巴特沃思数字带通滤波器设计：
抽样频率为10kHZ,,通带范围是1.5 kHZ到2.5 kHZ，通带衰减不大于3dB，在1kHZ和4kHZ处衰减不小于20 dB
（1）所给的技术指标：fp1=1500;fp2=2500;fs1=1000;fs2=4000;Fs=10000;Rp=3;Rs=20;
（2）归一化处理：Wp1=fp1/(Fs/2);Wp2=fp2/(Fs/2);Ws1=fs1/(Fs/2);Ws2=fs2/(Fs/2);
Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=buttord(wp,ws,Rp,Rs);
（4）设计巴特沃兹带通滤波器，调用系统函数：[B,A]=butter(N,Wn);
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。
3、巴特沃思数字带阻滤波器设计：
抽样频率为10kHZ,,在-2 dB衰减处的边带频率是1.5 kHZ，4 kHZ， 在-13dB衰减处频率是2kHZ和3kHZ
（1）所给的技术指标：fp1=1500;fp2=4000;fs1=2000;fs2=3000;Fs=10000;Rp=-2;Rs=-13;
（2）归一化处理：Wp1=fp1/(Fs/2);Wp2=fp2/(Fs/2);Ws1=fs1/(Fs/2);Ws2=fs2/(Fs/2);
Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=buttord(wp,ws,Rp,Rs);
（4）设计巴特沃兹带阻滤波器，调用系统函数：[B,A]=butter(N,Wn,'stop');
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。
4、切比雪夫1型数字高通滤波器设计：
抽样频率为10kHZ,,通带截止频率为2.5 kHZ，通带衰减不大于2dB，阻带上限截止频率1.5kHZ ,阻带衰减不小于15 dB。
设计步骤：
（1）所给的技术指标： fp=2500; fs=1500;Fs=10000; Rp=2;Rs=15
（2）归一化处理：wp=fp*2/Fs; ws=fs*2/Fs;
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=cheb1ord(wp,ws,Rp,Rs);
（4）设计切比雪夫1型高通滤波器，调用系统函数：[B,A]=cheby1(N,Rp,Wn,'high');
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。
5、切比雪夫1型数字带通滤波器设计：
抽样频率为10kHZ,,通带范围是1.5 kHZ到2.5 kHZ，通带衰减不大于3dB，在1kHZ和4kHZ处衰减不小于20 dB
（1）所给的技术指标：fp1=1500;fp2=2500;fs1=1000;fs2=4000;Fs=10000;Rp=3;Rs=20;
（2）归一化处理：Wp1=fp1/(Fs/2);Wp2=fp2/(Fs/2);Ws1=fs1/(Fs/2);Ws2=fs2/(Fs/2);
Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
（4）设计切比雪夫1型带阻滤波器，调用系统函数：[B,A]=cheby1(N,Rp,Wn);
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。
6、巴特沃思数字带阻滤波器设计：
抽样频率为10kHZ,,在-2 dB衰减处的边带频率是1.5 kHZ，4 kHZ， 在-13dB衰减处频率是2kHZ和3kHZ
（1）所给的技术指标：fp1=1500;fp2=4000;fs1=2000;fs2=3000;Fs=10000;Rp=2;Rs=13;
（2）归一化处理：Wp1=fp1/(Fs/2);Wp2=fp2/(Fs/2);Ws1=fs1/(Fs/2);Ws2=fs2/(Fs/2);
Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];
（3）获得数字滤波器的最小阶次N和滤波器的固有Wn
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
（4）设计切比雪夫1型带阻滤波器，调用系统函数：[B,A]=cheby1(N,Rp,Wn,'stop');
（5）数字滤波器的频率响应：[h,f]=freqz(B,A,256,Fs);
（6）显示幅频响应与相频响应。

注:设计到的函数用法
（1）巴特沃兹阶次选择：buttord
[N,Wn]=buttord(Wp,Ws,Rp,Rs):返回符合要求的数字滤波器的最小阶次N和滤波器的固有频率Wn（3dB频率）。参数Wp为通带截止频率，Ws为阻带截止频率，Rp为通带允许的最大衰减，Rs为阻带达到的最小衰减。Wp和Ws为归一化频率，其值在0到1之间，1对应抽样频率的一半。Rp和Rs的单位为dB。对于低通和高通滤波器，Wp和Ws都是标量；对于带通和带阻滤波器，Wp和Ws为1×2的向量。
（2）巴特沃兹滤波器：butter 函数
[B,A]=butter(N,Wn，‘ftype’)：设计一个N阶的数字滤波器，返回值B和A分别为滤波器的转移函数的分子分母多项式系数向量，Wn为滤波器的固有频率，为归一化频率，其范围在0-1之间，1对应抽样频率的一半。参数‘ftype’指定设计滤波器的类型。当参数ftype为‘high’的时候，设计的事巴特沃兹高通滤波器；当它为‘stop’时候，设计的是巴特沃兹带通滤波器，此时Wn必须为二维向量，Wn=[w1,w2]。
（3）切比雪夫1型滤波器： cheb1ord函数
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs):返回数字滤波器的最小阶次N和滤波器的固有频率Wn（3dB频率）。参数Wp为通带截止频率，Ws为阻带截止频率，Rp为通带允许的最大衰减，Rs为阻带达到的最小衰减。Wp和Ws为归一化频率，其值在0到1之间，1对应抽样频率的一半。Rp和Rs的单位为dB。对于低通和高通滤波器，Wp和Ws都是标量；对于带通和带阻滤波器，Wp和Ws为1×2的向量。
（4）[B,A]=cheby1(N,R,Wn,’ftype’):返回值B和A分别为滤波器的转移函数的分子分母多项式的系数向量，Wn为滤波器的固有频率，为归一化的频率，其值在0到1之间。参数‘ftype’来制定滤波器的类型，当参数ftype为high的时候，设计的是切比雪夫1型高通数字滤波器；当它为stop时候，设计的是切比雪夫1型数字带通滤波器，此时Wn必须为二维向量，Wn=[w1,w2].
（5）数字滤波器的频率响应：freqz
[h，f]=freqz（b，a，n，fs）：返回分子分母系数向量分别为b和a的数字滤波器的频率响应h和频率点向量f。其中向量h和f的长度均为n，采样频率指定为fs，单位为Hz，频率向量f的范围是[0,fs/2]。
（6）求相角：angle函数
P=angle（h）：以弧度为单位，返回向量h的相角。如果向量h为复数向量，则返回值在-π到π之间。

四、 设计结果及分析
1、 巴特沃兹高通滤波器主要程序代码：

clear all;
%巴特沃兹高通滤波器
Fsamp = 10e3;
Rp = 2.5e3;
Rs = 1.5e3;
Ap = 2;
As = 15;
wp = 2*pi*Fp/Fsamp; %归一化
ws = 2*pi*Fs/Fsamp;
Wp = tan(wp/2);
Ws = tan(ws/2);

[N, WN] = buttord(Ws,Wp,Ap,As,'s'); %获得滤波器的阶数与固有频率
[B, A] = butter(N, WN,'s'); %低通原型
[BT,AT] = lp2hp(B,A,Ws);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);

figure(1);
subplot(321)
plot(W/pi,20*log10(abs(H)));grid;
title('巴特沃兹带通滤波器')
ylabel('gain');
axis([0 1 -80 5]);
subplot(322)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');
2、 巴特沃兹带通和带阻数字滤波器与此类似，故省略。
3、 切比雪夫1型带阻滤波器的主要程序代码：
Fsamp = 10e3;
Fp1 = 1.5e3;
Fp2 = 4e3;
Fs1 = 2e3;
Fs2 = 3e3;
Ap = 2;
As = 13;
wp1 = 2*pi*Fp1/Fsamp;
wp2 = 2*pi*Fp2/Fsamp;
ws1 = 2*pi*Fs1/Fsamp;
ws2 = 2*pi*Fs2/Fsamp;
Wp1 = tan(wp1/2);
Wp2 = tan(wp2/2);
Ws1 = tan(ws1/2);
Ws2 = tan(ws2/2);
Bw = Ws2 - Ws1;
Wp2 = Ws1*Ws2/Wp1; %调整截止频率参数。
W0 = (Ws1*Ws2)^0.5;
Ws = 1;
Wp = (Bw*Wp1)/((Ws1*Ws2)-Wp1^2); %归一化，求低通原型的参数。
[N, WN] = cheb1ord(Wp,Ws,Ap,As,'s');
[B, A] = cheby1(N,Ap,WN,'s'); %低通原型
[BT,AT] = lp2bs(B,A,W0,Bw);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
subplot(325)
plot(W/pi,20*log10(abs(H)));grid;
title('切比雪夫1型带阻滤波器')
xlabel('w/pi'),ylabel('gain');
axis([0 1 -80 5]);
subplot(326)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');
4、切比雪夫1型高通和带阻滤波器与带通的程序代码类似。
5、图（1）为巴特沃兹高通、带通和带阻数字滤波器的幅频响应和相频响应。
图（2）为切比雪夫1型高通、带通和带阻数字滤波器的幅频响应和相频响应。

图（1）

图（2）
结论：由以上两图可以看出，巴特沃兹滤波器在通带内有最大化的平坦幅度响应，转换带没有切比雪夫陡峭，没有震荡带，具有一定的过冲和震荡，但不及切比雪夫。切比雪夫通带内有纹波。在相同的技术指标下的具体分析如下：
（1）采用巴特沃兹设计的高通滤波器的通带比较平淡，没有纹波，其相位变化也比较平淡。而采用切比雪夫设计的高通滤波器由于其通带内有纹波，所以其相位会有比较大的相对变化，并且其转换带比较陡。
（2）采用巴特沃兹设计的带通滤波器的通带比较平坦，没有纹波，其相位变化也比较平坦。而采用切比雪夫设计的带通滤波器由于其通带内有纹波，所以其相位会有比较大的相对变化。
（3）采用巴特沃兹设计的带阻滤波器的阻带比较陡峭，其相位会有比较大的相对变化。而采用切比雪夫设计的带阻滤波器通带变化相对缓慢，故其相位变化相对比较下。
（4）当采用相同的技术指标设计高通滤波器时最好使用切比雪夫设计，其过渡带比较短;当设计带通滤波器时，巴特沃兹与切比雪夫没有太大差别；当设计带阻滤波器时，最好采用巴特沃兹设计，其过渡带比较小。

五、 体会
通过本次设计，我对巴特沃兹与切比雪夫的设计指标有了更加明确的认识，以及在什么情况下适合选用巴特沃兹，在什么情况下采用切比雪夫设计有了深刻了解。并在学习中加深了对各个函数的使用，也学习到了如何写好一篇报告。

六、 参考文献
[1] 刘会灯、朱飞.MATLAB编程基础与典型应用.人民邮电出版社，2008.7
[2] 程佩青.数字信号处理教程.清华大学出版社，2009.11
[4] 万永革.数字信号处理的MATLAB实现.科学出版社，2008.9
[5] 郭仕剑、王宝顺、贺志国.MATLAB数字信号处理.人民邮电出版社，2007.11

附：
%高通IIR滤波器设计
Fsamp = 10e3;
Rp = 2.5e3;
Rs = 1.5e3;
Ap = 2;
As = 15;
wp = 2*pi*Fp/Fsamp;
ws = 2*pi*Fs/Fsamp;
Wp = tan(wp/2);
Ws = tan(ws/2);
[N, WN] = buttord(Ws,Wp,Ap,As,'s');
[B, A] = butter(N, WN,'s'); %低通原型
[BT,AT] = lp2hp(B,A,Ws);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
figure(1);
subplot(321)
plot(W/pi,20*log10(abs(H)));grid;
title('巴特沃兹带通滤波器')
ylabel('gain');
axis([0 1 -80 5]);
subplot(322)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');

%带通IIR滤波器设计
Fsamp = 10e3;
Fp1 = 1.5e3;
Fp2 = 2.5e3;
Fs1 = 1e3;
Fs2 = 4e3;
Ap = 3;
As = 20;
wp1 = 2*pi*Fp1/Fsamp;
wp2 = 2*pi*Fp2/Fsamp;
ws1 = 2*pi*Fs1/Fsamp;
ws2 = 2*pi*Fs2/Fsamp;
Wp1 = tan(wp1/2);
Wp2 = tan(wp2/2);
Ws1 = tan(ws1/2);
Ws2 = tan(ws2/2);
Bw = Wp2 - Wp1;
Ws2 = Wp1*Wp2/Ws1; %调整截止频率参数。
W0 = (Wp1*Wp2)^0.5;
Wp = 1;
Ws = -((Wp1*Wp2)-Ws2^2)/(Bw*Ws2); %归一化，求低通原型的参数。
[N, WN] = buttord(Wp,Ws,Ap,As,'s');
[B, A] = butter(N, WN,'s'); %低通原型
[BT,AT] = lp2bp(B,A,W0,Bw);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
subplot(323)
plot(W/pi,20*log10(abs(H)));grid;
title('巴特沃兹带通滤波器')
ylabel('gain');
axis([0 1 -80 5]);
subplot(324)
plot(W/pi,angle(H));grid;
ylabel('phase');

%带阻IIR滤波器设计
Fsamp = 10e3;
Fp1 = 1.5e3;
Fp2 = 4e3;
Fs1 = 2e3;
Fs2 = 3e3;
Ap = 2;
As = 13;
wp1 = 2*pi*Fp1/Fsamp;
wp2 = 2*pi*Fp2/Fsamp;
ws1 = 2*pi*Fs1/Fsamp;
ws2 = 2*pi*Fs2/Fsamp;
Wp1 = tan(wp1/2);
Wp2 = tan(wp2/2);
Ws1 = tan(ws1/2);
Ws2 = tan(ws2/2);
Bw = Ws2 - Ws1;
Wp1 = Ws1*Ws2/Wp2; %调整截止频率参数。
W0 = (Ws1*Ws2)^0.5;
Ws = 1;
Wp = (Bw*Wp1)/((Ws1*Ws2)-Wp1^2); %归一化，求低通原型的参数。
[N, WN] = buttord(Wp,Ws,Ap,As,'s');
[B, A] = butter(N, WN,'s'); %低通原型
[BT,AT] = lp2bs(B,A,W0,Bw);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
subplot(325)

plot(W/pi,20*log10(abs(H)));grid;
title('巴特沃兹带阻滤波器')
xlabel('w/pi'),ylabel('gain');
axis([0 1 -80 5]);
subplot(326)
plot(W/pi,angle(H));grid;

xlabel('w/pi'),ylabel('phase');
%高通IIR滤波器设计
Fsamp = 10e3;
Fp = 2.5e3;
Fs = 1.5e3;
Ap = 2;
As = 15;
wp = 2*pi*Fp/Fsamp;
ws = 2*pi*Fs/Fsamp;
Wp = tan(wp/2);
Ws = tan(ws/2);
[N, WN] = cheb1ord(Ws,Wp,Ap,As,'s');
[B, A] = cheby1(N,Ap,WN,'s'); %低通原型
[BT,AT] = lp2hp(B,A,Ws);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
figure(2)
subplot(321)
plot(W/pi,20*log10(abs(H)));grid;
ylabel('gain');
title('切比雪夫1型高通滤波器')
axis([0 1 -80 5]);
subplot(322)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');

%带通IIR滤波器设计
Fsamp = 10e3;
Fp1 = 1.5e3;
Fp2 = 2.5e3;
Fs1 = 1e3;
Fs2 = 4e3;
Ap = 3;
As = 20;
wp1 = 2*pi*Fp1/Fsamp;
wp2 = 2*pi*Fp2/Fsamp;
ws1 = 2*pi*Fs1/Fsamp;
ws2 = 2*pi*Fs2/Fsamp;
Wp1 = tan(wp1/2);
Wp2 = tan(wp2/2);
Ws1 = tan(ws1/2);
Ws2 = tan(ws2/2);
Bw = Wp2 - Wp1;
Ws2 = Wp1*Wp2/Ws1; %调整截止频率参数。
W0 = (Wp1*Wp2)^0.5;
Wp = 1;
Ws = -((Wp1*Wp2)-Ws2^2)/(Bw*Ws2); %归一化，求低通原型的参数。
[N, WN] = cheb1ord(Wp,Ws,Ap,As,'s');
[B, A] = cheby1(N,Ap,WN,'s'); %低通原型
[BT,AT] = lp2bp(B,A,W0,Bw);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
subplot(323)
plot(W/pi,20*log10(abs(H)));grid;
ylabel('gain');
title('切比雪夫1型带通滤波器')
axis([0 1 -80 5]);
subplot(324)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');

%带阻IIR滤波器设计
Fsamp = 10e3;
Fp1 = 1.5e3;
Fp2 = 4e3;
Fs1 = 2e3;
Fs2 = 3e3;
Ap = 2;
As = 13;
wp1 = 2*pi*Fp1/Fsamp;
wp2 = 2*pi*Fp2/Fsamp;
ws1 = 2*pi*Fs1/Fsamp;
ws2 = 2*pi*Fs2/Fsamp;
Wp1 = tan(wp1/2);
Wp2 = tan(wp2/2);
Ws1 = tan(ws1/2);
Ws2 = tan(ws2/2);
Bw = Ws2 - Ws1;
Wp2 = Ws1*Ws2/Wp1; %调整截止频率参数。
W0 = (Ws1*Ws2)^0.5;
Ws = 1;
Wp = (Bw*Wp1)/((Ws1*Ws2)-Wp1^2); %归一化，求低通原型的参数。
[N, WN] = cheb1ord(Wp,Ws,Ap,As,'s');
[B, A] = cheby1(N,Ap,WN,'s'); %低通原型
[BT,AT] = lp2bs(B,A,W0,Bw);
[num,den] = bilinear(BT,AT,0.5); %双线性变换
[H,W] = freqz(num,den,512);
subplot(325)
plot(W/pi,20*log10(abs(H)));grid;
title('切比雪夫1型带阻滤波器')
xlabel('w/pi'),ylabel('gain');
axis([0 1 -80 5]);
subplot(326)
plot(W/pi,angle(H));grid;
xlabel('w/pi'),ylabel('phase');