课程设计
题 目 IIR数字滤波器对语音信号的处理
学生姓名 宋小庆 学号 **********
所在院(系) 物电学院
专业班级 电子信息科学与技术081班
指导教师 刘东
完成地点 实验楼513教室
2011 年 12 月 15 日
IIR数字滤波器对语音信号的处理
作者：***
指导老师：***
一 绪论
数字滤波是数字信号处理的重要内容，是由乘法器、加法器和单位延时器组成的一种运算过程，其功能是对输入离散信号进行运算处理，以达到改变信号频谱的目的。数字滤波器根据频域特性可分为低、高通、带通和带阻四个基本类型；根据时域特性可分为无限长单位冲激响应IIR(Infinite Impulse Response)滤波器和有限长单位冲激响应FIR(Finite Impulse Response)滤波器。
数字滤波在通信、图像编码、语音编码、雷达等许多领域中有着十分广泛的应用。目前，数字信号滤波器的设计图像处理、数据压缩等方面的应用取得了令人瞩目的进展和成就。鉴于此，数字滤波器的设计就显得尤为重要。
MATLAB是美国MathWorks公司推出的一套用于工程计算的可视化高性能语言与软件环境。MATLAB为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。MATLAB推出的工具箱使各个领域的研究人员可以直观方便地进行科学研究、工程应用，其中的信号处理(signal processing)、图像处理(image processing),小波(wavelet)等工具箱为数字滤波研究的蓬勃发展提供了有力的工具。
二 数字滤波器
2.1什么是数字滤波器
滤波器是指用来对输入信号进行滤波的硬件和软件。所谓数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器和模拟滤波器相比，因为信号的形式和实现滤波的方法不同，数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。
一般用两种方法来实现数字滤波器:一是采用通用计算机，把滤波器所要完成的运算编成程序通过计算机来执行，也就是采用计算机软件来实现;二是设计专用的数字处理硬件。
MATLAB的信号处理工具箱是专门应用于信号处理领域的专用工具箱，它的两个基本组成就是滤波器的设计与实现部分以及谱分析部分。工具箱提供了丰富而简便的设计，使原来繁琐的程序设计简化成函数的调用。只要以正确的指标参数调用相应的滤波器设计程序或工具箱函数，便可以得到正确的设计结果，使用非常方便。
2.2数字滤波器的分类
数字滤波器从功能上分类:可分为低通滤波器、高通滤波器、带通滤波器、带阻滤波器。
从滤波器的网络结构或者从单位脉冲响应分类:可分为IIR滤波器(即无限长单位冲激响应滤波器)和FIR滤波器(即有限长单位冲激响应滤波器)。它们的函数分别为:
第一个公式中的H (z)称为N阶IIR滤波器函数，第二个公式中的H (z)称为(N-1)阶FIR滤波器函数。
2.3数字滤波器的设计要求
滤波器的指标常常在频域给出。数字滤波器的频响特性函数一般为复函数，所以通常表示为：
其中，||称为幅频特性函数，Φ(w)称为相频特性函数。幅频特性表示信号通过该滤波器后各频率成分的衰减情况，而相频特性反映各频率通过滤波器后在时间上的延时情况。一般IIR数字滤波器，通常只用幅频响应函数||来描述设计指标，相频特性一般不作要求。
IIR滤波器指标参数如下图所示。图中，ωp和ωs分别为通带边界频率和阻带边界频率；δ1和δ2分别为通带波纹和阻带波纹；允许的衰减一般用dB数表示，通带内所允许的最大衰减（dB）和阻带内允许的最小衰减（dB）分别为αp和αs表示：
一般要求：
当0≤∣ω︳≤ωp时，-20lg︳︳H（ejω）︳≤αp;
当ωs≤︳ω︳≤π时，αs≤-20lg︳H（ejω）︳。
2.4数字滤波器设计方法概述
IIR数字滤波器设计流程图如下：
步骤流程图
IIR滤波器设计方法有两类，经常用到的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计思路是:先设计模拟滤波器得到传输函数Ha(s)，然后将Ha(s)按某种方法转换成数字滤波器的系统函数H (Z)。这一类方法是基于模拟滤波器的设计方法相对比较成熟，它不仅有完整的设计公式，也有完整的图表供查阅。更可以直接调用MATLAB中的对应的函数进行设计。另一类是直接在频域或者时域中进行设计的，设计时必须用计算机作辅助设计，直接调用MATLAB中的一些程序或者函数可以很方便地设计出所需要的滤波器。
三 IIR滤波器的设计
3.1典型的IIR数字滤波器的设计
模拟滤波器的理论和设计方法已发展得相当成熟，且有一些典型的模拟滤波器供我们选择，如巴特沃斯(Butterworth )滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer )滤波器、贝塞尔(Bessel )滤波器等，这些典型的滤波器各有特点。
用MATLAB进行典型的数字滤波器的设计，一般步骤如下:
（1）将设计指标归一化处理。如果采用双线性变换法，还需进行预畸变。
（2）根据归一化频率，确定最小阶数N和频率参数Wn。可供选用的阶数择函数有:buttord，cheblord，cheb2ord，ellipord等。
（3）运用最小阶数N设计模拟低通滤波器原型。模拟低通滤波器的创建函数有:buttap,cheblap, cheb2ap,ellipap和besselap，这些函数输出的是零极点式形式，还要用zp2tf函数转换成分子分母多项式形式。如果想根据最小阶数直接设计模拟低通滤波器原型，可用butter,chebyl,cheby2,ellip,bessel等函数，只是注意要将函数中的Wn设为1。
（4）根据第2步的频率参数Wn，模拟低通滤波原型转换模拟低通、高通、带通、带阻滤波器，可用函数分别是：lp21p,lp2hp,lp2bp,lp2bs。
（5）运用脉冲响应不变法或双线性变法把模滤波器转数字滤波器，调用的函数是impinvar和bilinear。脉冲响应不变法适用于采样频率大于4倍截止频率的锐截止低通带通滤波器，而双线性变换法适合于相位特性要求不高的各型滤波器。
（6）根据输出的分子分母系数，用tf函数生成H（z）的表达式，再用freqz函数验证设计结果。
设计巴特沃思数字低通滤波器和椭圆数字低通滤波器，要求通带边界频率fp=2.1kHZ，通带最大衰减Rp=0.5dB；阻带边界频率fs=8kHZ，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。
用脉冲响应不变法设计的巴特沃思数字低通滤波器的M程序如下：
fp=2100;
fs=8000;
Fs=20000;
Rp=0.5;
Rs=30;
T=1/Fs; %设计指标
W1p=fp/Fs*2;W1s=fs/Fs*2;%求归一化频率
[N,Wn]=buttord(W1p,W1s,Rp,Rs,'s');
%确定butterworth的最小介数N和频率参数Wn
[z,p,k]=buttap(N); %设计模拟低通原型的零极点增益参数
[bp,ap]=zp2tf(z,p,k); %将零极点增益转换成分子分母参数
[bs,as]=lp2lp(bp,ap,Wn*pi*Fs);%将低通原型转换为模拟低通
[bz,az]=impinvar(bs,as,Fs); %用脉冲响应不变法进行模数变换
sys=tf(bz,az,T); %给出传输函数H(Z)
[H,W]=freqz(bz,az,512,Fs); %生成频率响应参数
subplot(2,1,1);
plot(W,20*log10(abs(H))); %绘制幅频响应
grid on; %加坐标网格
xlabel('频率/Hz');
ylabel('振幅/dB');
subplot(2,1,2);
plot(W,abs(H)); grid on;
xlabel('频率/Hz');
ylabel('振幅/H');
运行后的波形如下：
运行结果：
N =4
bz = 0.0000 0.0999 0.1914 0.0252
az= 1.0000 -1.4336 1.0984 -0.4115 0.0627
用双线性变换法设计的椭圆数字低通滤波器的M程序如下：
fs=20000;
wp=2*pi*2100/fs;
ws=2*pi*8000/fs;
Rp=0.5;
Rs=30;
Ts=1/fs;
Wp=2/Ts*tan(wp/2);Ws=2/Ts*tan(ws/2); %按频率转换公式进行转换
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数
[z,p,k]=ellipap(N,Rp,Rs);%设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转化
[bz,az]=bilinear(b,a,fs); %运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,512,fs);
subplot(2,1,1);
plot(f,20*log10(abs(H)));
title('N=2 频率响应');
grid on;
xlabel('频率/Hz');
ylabel('振幅/dB');
subplot(2,1,2);
plot(f,abs(H)); grid on;
xlabel('频率/Hz');
ylabel('振幅/H');
运行结果：
N=2
bz= 0.1213 0.1662 0.1213
az= 1.0000 -0.98 0.4218
3.2完全滤波器设计
除了典型设计以外，MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，直接调用就可以设计滤波器，这为设计通用滤波器提供了方便。
设计Butterworth滤波器用函数butter()，可以设计低通、高通、带通和带阻的数字和模拟滤波器，其特性是通带内的幅度响应最大限度的平滑，但损失了截止频率处的下降斜度。
设计Chebyshev I型滤波器用函数chebyl()。可以设计低通、高通、带通和带阻的数字和模拟Chebyshev I型滤波器，其通带内为等波纹，阻带内为单调。Chebyshev I型滤波器的下降斜度比II型大，但其代价目是通带内波纹较大。
设计Chebyshev II型滤波器用函数cheby2()。可以设计低通、高通、带通和带阻的数字和模拟Chebyshev II型滤波器，其通带内为单调，阻带内等波纹。Chebyshev II型滤波器的下降斜度比I型小，但其阻带内波纹较大。
设计椭圆滤波器用函数ellip()，与chebyl, cheby2类似，可以设计低通、高通、带通和带阻的数字和模拟滤波器。与Butterworth和chebyshev滤波器相比，ellip函数可以得到下降斜度更大的滤波器，得通带和阻带均为等波纹。一般情况下，椭圆滤波器能以最低的阶实现指定的性能指标。
在使用各类滤波器函数时应当注意以下重点:
A、阶数和固有频率的选择:[N,Wn]=buttord(Wp,Ws,Rp,Rs)可得到符合要求性质的滤波器的最小阶数N以及数字Butterworth滤波器的固有频率Wn(即3dB )。设计的要求是在通带内的衰减不超过Rp,在阻带内的衰减不小于Rs，通带和阻带有截止频率分别是Wp, Ws,它们是归一化的频率，范围是[0, 1]，对应π弧度。
B、有关滤波器设计当中的频率归一化问题:信号处理工具箱中经常使用的频率是Nyquist频率，它被定义为采样频率的一半，在滤波器的阶数选择和设计中的截止频率均使用Nyquist频率进行归一化处理。例如对于一个采样频率为1000 Hz的系统，400Hz的归一化即为400/500=0.8。归一化频率的范围在[0, 1]之间。如果要将归一化频率转换为角频率，则将归一化频率乘以π;如果要将归一化频率转换为Hz,则将归一化频率乘以采样频率的一半。
C、设计一个N阶的低通Butterworth滤波器使用函数[B,A]=butter(N, Wn)，返回滤波器系数矩阵[B,A]。其中固有频率Wn必须是归一化频率。它的最大值是采样频率的一半。Fs缺省时默认为2Hz。如果Wn=[Wl,W2]是一个两元素的向量，则函数将设计出一个2N阶的带通滤波器，通带为[W1,W2]。
设计Chebyshev I型和Chebyshev II型数字低通滤波器，要求通带边界频率fp=2.1kHZ，通带最大衰减Rp=0.5dB；阻带边界频率fs=8kHZ，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。
Chebyshev I型的M程序如下：
Fs=20000; %抽样频率20KHz
Flp=2100;
Fls=8000;
Wp=2*Flp/Fs; %归一化的通带截止频率
Ws=2*Fls/Fs; %归一化的阻带截止频率
Rp=0.5; %通带最大衰减（单位：dB）
Rs=30; %阻带最小衰减（单位：dB）
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率
[b,a]=cheby1(N,Rp,Wn); %返回H(z)的分子分母系数
[hw,w]=freqz(b,a);
subplot(2,1,1);
plot(w/pi,20*log10(abs(hw)));grid on;
xlabel('ω/π');ylabel('幅度（dB）')
title('切比雪夫I型幅频响应');
subplot(2,1,2);plot(w/pi,abs(hw));
grid on;
xlabel('ω/π');ylabel('幅度（H）');
运行后的波形如下：
运行结果：
N=2
b= 0.1007 0.2014 0.1007
a= 1.0000 -0.9872 0.4140
Chebyshev II型M程序如下：
Fs=20000; %抽样频率20KHz
Flp=2100;
Fls=8000;
Wp=2*Flp/Fs; %归一化的通带截止频率
Ws=2*Fls/Fs; %归一化的阻带截止频率
Rp=0.5; %通带最大衰减（单位：dB）
Rs=30; %阻带最小衰减（单位：dB）
[N,Wn]=cheb2ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率
[b,a]=cheby2(N,Rs,Wn); %返回H(z)的分子分母系数
[hw,w]=freqz(b,a);
subplot(2,1,1);
plot(w/pi,20*log10(abs(hw)));grid on;
xlabel('ω/π');ylabel('幅度（dB）')
title('切比雪夫II型幅频响应');
subplot(2,1,2);plot(w/pi,abs(hw));
grid on;
xlabel('ω/π');ylabel('幅度（H）');
运行后的波形：
运行结果：
N=2
b= 0.2357 0.4241 0.2357
a= 1.0000 -0.2996 0.1950
3.3 结果分析
从频率响应图中可以看出:巴特沃斯滤波器具有单调下降的幅频特性，通带内平滑;切比雪I型滤波器的幅频特性在通带内有波动，阻带内单调;chebyshev II型滤波器的幅频特性在阻带内有波动，通带内单调;椭圆滤波器的选择性相对前三种是最好的，下降斜度比较大，通带和阻带内均为等波纹，同样的性能指标，椭圆滤波器可以最低的阶数来实现。这样根据不同的要求可以选用不同类型的滤波器。
四 对语音信号的滤波
利用函数wavread对语音信号的采集：
i=1;
[x,fs,bits]=wavread('1.wav');
%x：语音数据；fs：采样频率；bits：采样点数
sound(x,fs,bits); %话音回放
N=length(x);
n=0:N-1;
figure(i);
subplot(2,1,1);
plot(n,x); %画出原始语音信号的波形
xlabel('n');
ylabel('x(n)');
title('原始语音信号');
subplot(2,1,2);
[H,f]=freqz(x,1,512,fs);
plot(f,20*log10(abs(H))); %画出原始语音信号的频谱
xlabel('n');
ylabel('x(n)');
title('原始语音信号的频谱');
运行结果：
fs=44100
bits=16
加入上面设计的椭圆滤波器对其进行滤波：
fs=20000;
wp=2*pi*2100/fs;
ws=2*pi*8000/fs;
Rp=0.5;
Rs=30;
Ts=1/fs;
Wp=2/Ts*tan(wp/2);Ws=2/Ts*tan(ws/2); %按频率转换公式进行转换
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数
[z,p,k]=ellipap(N,Rp,Rs);%设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转化
[bz,az]=bilinear(b,a,fs); %运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,512,fs);
plot(f,20*log10(abs(H)));
title('N=2 频率响应');
grid on;
xlabel('频率/Hz');
ylabel('振幅/dB');
%-------滤波前后波形比较-----%
i =i+1;
figure(i);
y=filter(bz,az,x);
subplot(2,2,1);
plot(n,x); %画出原语音信号的波形
xlabel('n');
ylabel('x(n)');
title('原始语音信号');
subplot(2,2,3);
plot(n,y); %IIRLowPass后的波形
title('IIRLowPass后的波形');
%---------滤波前后语音信号频谱-----%
subplot(2,2,2);
[H,f]=freqz(x,1,512,fs);
plot(f,20*log10(abs(H))); %滤波前语音信号频谱
xlabel('频率/Hz');
ylabel('振幅/dB');
title('IIR低通滤波前语音信号频谱');
grid on;
subplot(2,2,4);
[H,f]=freqz(y,1,512,fs);
plot(f,20*log10(abs(H))); %滤波后语音信号频谱
xlabel('频率/Hz');
ylabel('振幅/dB');
title('IIR低通滤波后语音信号频谱');
grid on;
sound(y,fs,bits);
在matlab上运行后，发现语音信号在经过上面设计的椭圆数字低通滤波器后，声音发生了明显的改变。观察波形，因为导入的语音信号波形的频率分布绝大多数分布在低频区，故没有被滤波器过滤掉多少，所以语音信号波形变化不大。
但频谱发生了明显得改变。
同理，只要把之前设计的其他三个滤波器的程序替换掉如上的椭圆滤波器的程序，就可得到语音信号在不同滤波器滤波后的波形变化情况。
心得体会
通过2个星期紧张的课程设计，不用说，收获肯定是不少的。虽然之前已经多次接触过了matlab这款软件，单通过此次运用matlab来设计IIR滤波器，让我对matlab在数字信号处理中的应用有了更深的了解，真是受益匪浅。
本次课程设计中，有机地结合了理论与实践，既考察了我们对理论知识的掌握情况，还反映出我们实际动手能力和编程能力，更主要的是它激起我们创新思维，提高了自己分析问题和解决问题的能力，这就是它的吸引力所在，我会为了能编写出程序好几个小时都在研究如何编程，这在无形中以及提高了我各方面的能力。
时间虽然短暂，但这次课程设计所给的绝不会是短暂的效果，无论是在知识上，还是在思想上都给我烙下了深刻的印象，我想，这次对于这次课程设计我获得的已经满足了，谢谢老师在这段时间的指导。
参考文献
[1]赵树杰等，数字信号处理，西安电子科技大学出版社，1997.10
[2] 陈怀琛等，MATLAB及在电子信息课程中的应用，电子工业出版社出版，2002.4
[3]李丽 王振领，MATLAB工程计算及应用，人民邮电出版社，2001.9
[4] 陈怀琛等译，《数字信号处理及其MATLAB实现》，电子工业出版社；
/下载本文
