课程设计（论文）任务及评语
院（系）教研室：通信教研室
学号学生姓名专业班级
课程设计
（论文）
题目
音频滤波中ButterWorth法IIR数字滤波器的设计
课程设计（论文）任务
1．音频信号的频率范围为30-20000Hz。设计一IIR滤波器，完成对带噪音频信号的滤波。噪声信号频率为200Hz；设采样频率为441000Hz；
2. 理解并掌握ButterWorth法设计IIR滤波器的工作原理。
3．实现ButterWorth法IIR滤波器的仿真设计；
4．用该滤波器完成对带噪音频信号的滤波；
5．调试源文件，观察、分析并打印设计的波形。
6．完成本次设计，填写设计指导书。
指
导
教
师
评
语
及
成
绩
平时成绩：论文成绩：指导教师签字：
答辩成绩：总成绩：年月日
目录
1 引言 (1)
1.1数字滤波器简介 (1)
1.2 MATLAB软件介绍 (1)
1.3巴特沃斯滤波器简介 (2)
2 设计要求与思路 (3)
2.1设计要求 (3)
2.2设计思路 (3)
3 基于MATLAB的数字滤波器设计 (6)
3.1 语音信号的采集与频谱分析 (6)
3.2 噪音信号的构建与分析 (7)
3.3 ButterWorth法IIR滤波器的设计 (8)
3.4 加噪信号滤波与频谱分析 (9)
4 总结 (11)
参考文献 (12)
附录 (13)
1 引言
1.1数字滤波器简介
数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。
数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，其差分方程为：
)()()(1
i n y i n x n y N
i i N
i i b a -+==∑∑==
系统函数为： ∑∑=---+=
N
k K
k M
r r
r Z
a Z
b z H 1
01)(
1.2 MATLAB 软件介绍
Matlab 是一个交互式的系统，它的基本运算单元是不需指定维数的矩阵，按照IEEE 的数值计算标准（能正确处理无穷数Inf(Infinity)、无定义数NaN(not-a-number)及其运算）进行计算。系统提供了大量的矩阵及其它运算函数，可以方便地进行一些很复杂的计算，而且运算效率极高。Matlab 命令和数学中的符号、公式非常接近，可读性强，容易掌握，还可利用它所提供的编程语言进行编程完成特定的工作。除基本部分外，Matlab 还根据各专门领域中的特殊需要提供了许多可选的工具箱，如应用于自动控制领域的Control System 工具箱和神经网络中Neural Network 工具箱等。
MATLAB 主要特点：
(1)语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB 程序书写形式自由，利用起丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于
库函数都由本领域的专家编写，用户不必担心函数的可靠性。可以说，用MATLAB 进行科技开发是站在专家的肩膀上。
(2)运算符丰富。由于MATLAB 是用C 语言编写的，MATLAB 提供了和C 语言几乎一样多的运算符，灵活使用MATLAB 的运算符将使程序变得极为简短。
(3)MATLAB 既具有结构化的控制语句（如for 循环，while 循环，break 语句和if 语句），又有面向对象编程的特性。
(4)程序限制不严格，程序设计自由度大。例如，在MATLAB 里，用户无需对矩阵预定义就可使用。
(5)程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。
1.3巴特沃斯滤波器简介
巴特沃斯逼近又称最平幅度逼近。巴特沃斯低通滤波器幅度平方函数定义为
式中N 为正整数，代表滤波器的阶次。
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在
阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为 每倍频12分贝、 三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的 形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。
N
c a j H 22
11)(???
?
??ΩΩ+=Ω
2 设计要求与思路
2.1设计要求
1．音频信号的频率范围为30-20000Hz 。设计一IIR 滤波器，完成对带噪音频信号的滤波。 噪声信号频率为200Hz ；设采样频率为441000Hz ；
2. 理解并掌握ButterWorth 法设计IIR 滤波器的工作原理。 3．实现ButterWorth 法IIR 滤波器的仿真设计； 4．用该滤波器完成对带噪音频信号的滤波； 5．调试源文件，观察、分析并打印设计的波形。 6．完成本次设计，填写设计指导书。 2.2设计思路
IIR 滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab 都提供相应的工具箱函数，使得IIR 数字滤波器设计变得非常简单。
从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S 域变换到Z 域，而得到所需类型的数字滤波器。
图2.1 先频率变换再离散
冲击响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响冲激相应不变法是从时域出发，要求数字滤波器的激响应h(n)对应于模拟滤波器ha(t)的等间隔抽样，h(n)=ha(nT) ，其中T 是抽样周期，因此时域逼近良好。
如果令Ha(s)是ha(t)的拉普拉斯变换，H(z)为h(n)的Z 变换，利用采样序列的Z 变换与模拟信号的拉普拉斯变换的关系得
模拟低通
滤波器原
模拟带阻滤波器
数字带阻
滤波器
模拟域 频率变换 冲激响应不变法 双线性变换法
则可看出，脉冲响应不变法将模拟滤波器的S 平面变换成数字滤波器的Z 平面，这个从s 到z 的变换z=esT 是从S 平面变换到Z 平面的标准变换关系式。
图2.2脉冲响应不变法的映射关系
由（1）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为
(2) 这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率
(3) 才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即
|ω|<π (4) 但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图2所示。这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。
??
? ??-=Ω-=∑∑∞-∞=∞-∞==k T j s X T jk s X T z X k a s k a e z sT
π21)(1)(j Ω3π / T
π / T
－3π / T －π / T o o σ－11jIm[z ]
Re[z ]
Z 平面
S 平面
?
?
? ??-=∑∞-∞=T k j H T e H k a j πωω
21)(0
)(=Ωj H a ?
?
?
??=T j H T e H a j ωω1)(2
||s T Ω=≥Ωπ
图2.3脉冲响应不变法中的频响混叠现象
对某一模拟滤波器的单位冲激响应ha(t)进行采样，采样频率为fs ，若使fs 增加，即令采样时间间隔（T=1/fs ）减小，则系统频率响应各周期延拓分量之间相距更远，因而可减小频率响应的混叠效应。
冲激相应不变法是从时域出发，要求数字滤波器的激响应h(n)对应于模拟滤波器ha(t)的等间隔抽样，h(n)=ha(nT) ，其中T 是抽样周期，因此时域逼近良好。
优点：（1）h(n)完全模仿模拟滤波器的单位抽样响应时域逼近良好。 （2）线性相位模拟滤波器转变为线性相位数字滤波器。
缺点：（1）对时域的采样会造成频域的“混叠效应”，故有可能使所设计数字滤波器的频率响应与原来模拟滤波器的频率响应相差很大。
（2）不能用来设计高通和带阻滤波器。只适用于限带的低通、带通滤波器的设计。
－3
π－2π
……
)
j (a ΩH Ω
o o
－π
2π
3π
π
ω＝Ω T
)
(e j ωH T
π
2T
πT πT
π2－
3 基于MATLAB的数字滤波器设计
3.1 语音信号的采集与频谱分析
随即录取一段.wav语音文件，以文件名shiyan存进D盘内。在MATLAB软件平台下，利用wavread函数对语音信号进行采样，并记住采样频率和采样数。
wavread 函数调用格式：
y=wavread('D:\shiyan.wav'); %读取file所规定的wav文件，返回采样值放在向量y 中。
[y,fs,nbits]=wavread('D:\shiyan.wav'); %采样值放在向量y中，其中fs表示采样频率（hz），nbits 表示采样位数。
Y=fft(y); %快速傅里叶变换，对语音信号进行频谱分析。
首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在MATLAB中利用FFT函数对信号进行快速傅里叶变换，得到信号的频谱特性。
程序：
[y,fs,nbits]=wavread('D:\shiyan.wav');
M=length(y);
y=wavread('D:\shiyan.wav');%把语音信号进行加载入Matlab 仿真软件平台中
sound(y,fs,nbits); %对加载的语音信号进行回放
Y=fft(y); %快速傅里叶变换
figure(1)
subplot(2 ,1 ,1),plot(y);title('原始信号波形');
subplot(2 ,1 ,2),plot(abs(Y));title('原始信号频谱');
其中，采样频率fs=441000hz,nbits=32Bit。
程序的结果如下图所示：
图3.1 语音信号的采集
3.2 噪音信号的构建与分析
利用MATLAB中的随机函数(randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其进行频谱分析。程序如下
N=length(y); %求出语音信号的长度
noise=0.04*rand(N,2); %噪声信号的函数
z=fft(noise); %快速傅里叶变换
subplot(2 ,1 ,1),plot(noise);title('噪声信号波形');
subplot(2 ,1 ,2),plot(abs(z));title('噪声信号频谱');
axis([0,1200,0,5]);
图3.2 噪声信号
构建噪音信号后，就可以添加到到语音信号中，构成实验用的待滤波信号：s=y+noise; %噪声信号的叠加
figure(3)
subplot(2,1,1);plot(s);title ('滤波前的时域波形');
S=fft(s);
subplot(2,1,2);plot(abs(S));title ('滤波前的频域波形');
图3.3 滤波前波形分析
3.3 ButterWorth法IIR滤波器的设计
程序：
fp1=30;
fp2=20000; %通带频率
fb1=190;
fb2=2100; %阻带频率
Fs = 441000; %抽样频率
Rp=3; %通带最大衰减
Rs=20; %阻带最小衰减
wp1=2*pi*fp1/Fs;
wp2=2*pi*fp2/Fs;
ws1=2*pi*fb1/Fs;
ws2=2*pi*fb2/Fs;
wp=[wp1 wp2],ws=[ws1 ws2];
[N,wc]=buttord(wp,ws,Rp,Rs); %巴特沃斯滤波器 [B,A]=butter(N,wc,'stop'); % 带阻滤波器 [h,w]=freqz(B,A); figure(4)
subplot(211),plot(abs(h)); grid;
title('滤波器的性能分析'); axis([0,20,0,2]); 程序结果如下图：
0.1
0.2
0.30.40.50.60.70.80.9
1
-400
-300-200-1000
Normalized Frequency (?π rad/sample)
P h a s e (d e g r e e s )
00.10.2
0.30.40.50.60.70.80.9
1
-100
-50
50Normalized Frequency (?π rad/sample)
M a g n i t u d e (d B )
图3.4 加噪语音信号的波形
3.4 加噪信号滤波与频谱分析 o=filter(B,A,s) figure(5)
subplot(2,1,1);plot(o);title('滤波后时域波形');
subplot(2,1,2);plot(s);title('滤波前时域波形');
图3.5 滤波前后对比分析
4 总结
通过近两周多系统学习与设计，我投入了最大的热情和精力，从构建框架到程序设计，再进行MATLAB仿真，这其中也遇到不少困难与问题，经过同学的帮助与老师不断的细心指导让我完成了这次课程设计。通过这次课程设计，让我知道了MATLAB的实用性与重要性，也增加了我的动手能力，让我在独立思考中发现问题，解决问题。通过MATLAB 与数字信号处理的系统结合，让我更好的领会学习中不能理解的问题，加深自己对课本知识的更好理解，让我充分理解了理论结合于实践的道理，不再死学，做到了最好的活学活用。这次课设，也让我知道了文献的查找，不再拘泥于见得百度搜索，增进了我的视角，加大了我的知识面。
通过这次课设，感谢老师的指导与督促，让我明白了什么叫学到了真东西，我会在以后的生活与学习中最好的体现这种精神。
参考文献
[1]程佩青.数字信号处理教程.北京：清华大学出版社出版，2001
[2]丁玉美.数字信号处理学习指导与题解.西安：电子工业出版社出版, 2007
[3]吴大正.MATLAB及在电子信息教程中的应用.北京：电子工业出版社,2005
[4]王家文,王皓等.MATLAB7.0编程基础.北京：机械工业出版社,2005
[5]柳春锋.用MATLAB语言实现IIR滤波器的设计[J]. 齐齐哈尔大学学报, 2001,(04) .
[6]张玉. 利用MATLAB设计巴特沃斯数字滤波器[J]. 本溪冶金高等专科学校学报, 2003,(03)
附录：
[y,fs,nbits]=wavread('D:\shiyan.wav');
M=length(y);
y=wavread('D:\shiyan.wav'); %把语音信号进行加载入Matlab 仿真软件平台中sound(y,fs,nbits); %对加载的语音信号进行回放
Y=fft(y); %快速傅里叶变换
figure(1)
subplot(2 ,1 ,1),plot(y);title('原始信号波形');
subplot(2 ,1 ,2),plot(abs(Y));title('原始信号频谱');
M=length(y); %求出语音信号的长度
f=200;t=2;
noise=sin(2*pi*f*t)*rand(M,2); %噪声信号的函数
z=fft(noise); %快速傅里叶变换
figure(2)
subplot(2 ,1 ,1),plot(noise);title('噪声信号波形');
subplot(2 ,1 ,2),plot(abs(z));title('噪声信号频谱');
axis([0,250000,0,100]);
s=y+noise; %噪声信号的叠加
figure(3)
subplot(2,1,1);plot(s);title ('滤波前的时域波形');
S=fft(s);
subplot(2,1,2);plot(abs(S));title ('滤波前的频域波形');
fp1=30;
fp2=20000; %通带频率
fb1=190;
fb2=2100; %阻带频率
Fs = 441000; %抽样频率
Rp=3; %通带最大衰减
Rs=20; %阻带最小衰减
wp1=2*pi*fp1/Fs;
wp2=2*pi*fp2/Fs;
ws1=2*pi*fb1/Fs;
ws2=2*pi*fb2/Fs;
wp=[wp1 wp2],ws=[ws1 ws2];
[N,wc]=buttord(wp,ws,Rp,Rs); %巴特沃斯滤波器[B,A]=butter(N,wc,'stop'); % 带阻滤波器h=freqz(B,A);
figure(4)
subplot(211),plot(abs(h));
grid;
title('滤波器的性能分析');
axis([0,20,0,2]);
o=filter(B,A,s)
figure(5)
subplot(2,1,1);plot(o);title('滤波后时域波形');
subplot(2,1,2);plot(s);title('滤波前时域波形');