基于Matlab的IIR数字滤波器的设计别人实验

1、基于Matlab的IIR数字滤波器的设计一 绪论数字滤波是数字信号处理的重要内容，是由乘法器、加法器和单位延时器组成的一种运算过程，其功能是对输入离散信号进行运算处理，以达到改变信号频谱的目的。数字滤波器根据频域特性可分为低、高通、带通和带阻四个基本类型；根据时域特性可分为无限长单位冲激响应IIR(Infinite Impulse Response)滤波器和有限长单位冲激响应FIR(Finite Impulse Response)滤波器。数字滤波在通信、图像编码、语音编码、雷达等许多领域中有着十分广泛的应用。目前，数字信号滤波器的设计图像处理、数据压缩等方面的应用取得了令人瞩目的进展和成就。鉴

2、于此，数字滤波器的设计就显得尤为重要。 MATLAB是美国MathWorks公司推出的一套用于工程计算的可视化高性能语言与软件环境。MATLAB为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。MATLAB推出的工具箱使各个领域的研究人员可以直观方便地进行科学研究、工程应用，其中的信号处理(signal processing)、图像处理(image processing),小波(wavelet)等工具箱为数字滤波研究的蓬勃发展提供了有力的工具。二 数字滤波器2.1什么是数字滤波器滤波器是指用来对输入信号进行滤波的

3、硬件和软件。所谓数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器和模拟滤波器相比，因为信号的形式和实现滤波的方法不同，数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。 一般用两种方法来实现数字滤波器:一是采用通用计算机，把滤波器所要完成的运算编成程序通过计算机来执行，也就是采用计算机软件来实现;二是设计专用的数字处理硬件。MATLAB的信号处理工具箱是专门应用于信号处理领域的专用工具箱，它的两个基本组成就是滤波器的设计与实现部分以及谱分析部分。工具箱提供了丰富而简便的设计，使原来

4、繁琐的程序设计简化成函数的调用。只要以正确的指标参数调用相应的滤波器设计程序或工具箱函数，便可以得到正确的设计结果，使用非常方便。2.2数字滤波器的分类 数字滤波器从功能上分类:可分为低通滤波器、高通滤波器、带通滤波器、带阻滤波器。 从滤波器的网络结构或者从单位脉冲响应分类:可分为IIR滤波器(即无限长单位冲激响应滤波器)和FIR滤波器(即有限长单位冲激响应滤波器)。它们的函数分别为:第一个公式中的H (z)称为N阶IIR滤波器函数，第二个公式中的H (z)称为(N-1)阶FIR滤波器函数。2.3数字滤波器的设计要求 滤波器的指标常常在频域给出。数字滤波器的频响特性函数一般为复函数，所以通常表

5、示为：其中，|称为幅频特性函数，(w)称为相频特性函数。幅频特性表示信号通过该滤波器后各频率成分的衰减情况，而相频特性反映各频率通过滤波器后在时间上的延时情况。一般IIR数字滤波器，通常只用幅频响应函数|来描述设计指标，相频特性一般不作要求。IIR滤波器指标参数如下图所示。图中，p和s分别为通带边界频率和阻带边界频率；1和2分别为通带波纹和阻带波纹；允许的衰减一般用dB数表示，通带内所允许的最大衰减（dB）和阻带内允许的最小衰减（dB）分别为p和s表示： 一般要求：低通滤波器的技术要求2.4数字滤波器设计方法概述IIR数字滤波器的设计步骤流程图如下：步骤流程图 IIR滤波器设计方法有两类，经常

6、用到的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计思路是:先设计模拟滤波器得到传输函数Ha(s)，然后将Ha(s)按某种方法转换成数字滤波器的系统函数H (Z)。这一类方法是基于模拟滤波器的设计方法相对比较成熟，它不仅有完整的设计公式，也有完整的图表供查阅。更可以直接调用MATLAB中的对应的函数进行设计。另一类是直接在频域或者时域中进行设计的，设计时必须用计算机作辅助设计，直接调用MATLAB中的一些程序或者函数可以很方便地设计出所需要的滤波器。三 IIR滤波器的设计3.1典型的IIR数字滤波器的设计模拟滤波器的理论和设计方法已发展得相当成熟，且有一些典型的模拟滤波器供我们选择，如

7、巴特沃斯(Butterworth )滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer )滤波器、贝塞尔(Bessel )滤波器等，这些典型的滤波器各有特点。用MATLAB进行典型的数字滤波器的设计，一般步骤如下:（1） 将设计指标归一化处理。如果采用双线性变换法，还需进行预畸变。（2） 根据归一化频率，确定最小阶数N和频率参数Wn。可供选用的阶数择函数有:buttord，cheblord，cheb2ord，ellipord等。（3） 运用最小阶数N设计模拟低通滤波器原型。模拟低通滤波器的创建函数有:buttap,cheblap, cheb2ap,ellipap和besselap，

8、这些函数输出的是零极点式形式，还要用zp2tf函数转换成分子分母多项式形式。如果想根据最小阶数直接设计模拟低通滤波器原型，可用butter,chebyl,cheby2,ellip,bessel等函数，只是注意要将函数中的Wn设为1。（4） 根据第2步的频率参数Wn，模拟低通滤波原型转换模拟低通、高通、带通、带阻滤波器，可用函数分别是：lp21p,lp2hp,lp2bp,lp2bs。（5） 运用脉冲响应不变法或双线性变法把模滤波器转数字滤波器，调用的函数是impinvar和bilinear。脉冲响应不变法适用于采样频率大于4倍截止频率的锐截止低通带通滤波器，而双线性变换法适合于相位特性要求不高的

9、各型滤波器。（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);%

10、确定butterworth的最小介数N和频率参数Wnz,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; %加坐标

11、网格xlabel(频率/Hz); ylabel(振幅/dB);subplot(2,1,2); plot(W,abs(H); grid on; xlabel(频率/Hz); ylabel(振幅/H);运行后的波形如下：运行结果：N =4bz = 0.0000 0.0999 0.1914 0.0252az= 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/

12、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(

13、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=2bz= 0.1213 0.1662 0.1213az= 1.0000 -0.9889 0.42183.2完全滤波器设计除了典型设计以外，MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，直接调用就可以设计滤波器，这为设计通用滤波器提供了方便。 设计Butterworth滤波器用函数butter()，可以设计低

14、通、高通、带通和带阻的数字和模拟滤波器，其特性是通带内的幅度响应最大限度的平滑，但损失了截止频率处的下降斜度。 设计Chebyshev I型滤波器用函数chebyl()。可以设计低通、高通、带通和带阻的数字和模拟Chebyshev I型滤波器，其通带内为等波纹，阻带内为单调。Chebyshev I型滤波器的下降斜度比II型大，但其代价目是通带内波纹较大。 设计Chebyshev II型滤波器用函数cheby2()。可以设计低通、高通、带通和带阻的数字和模拟Chebyshev II型滤波器，其通带内为单调，阻带内等波纹。Chebyshev II型滤波器的下降斜度比I型小，但其阻带内波纹较大。 设

15、计椭圆滤波器用函数ellip()，与chebyl, cheby2类似，可以设计低通、高通、带通和带阻的数字和模拟滤波器。与Butterworth和chebyshev滤波器相比，ellip函数可以得到下降斜度更大的滤波器，得通带和阻带均为等波纹。一般情况下，椭圆滤波器能以最低的阶实现指定的性能指标。在使用各类滤波器函数时应当注意以下重点:A、阶数和固有频率的选择:N,Wn=buttord(Wp,Ws,Rp,Rs)可得到符合要求性质的滤波器的最小阶数N以及数字Butterworth滤波器的固有频率Wn(即3dB )。设计的要求是在通带内的衰减不超过Rp,在阻带内的衰减不小于Rs，通带和阻带有截止频

16、率分别是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)，返回滤

17、波器系数矩阵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; %抽样频率20KHzFlp=2100; Fls=8000;Wp=2*Flp/Fs; %归一化的通带截止频率Ws=2*

18、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(/);ylab

19、el(幅度（H）);运行后的波形如下：运行结果：N=2b= 0.1007 0.2014 0.1007a= 1.0000 -0.9872 0.4140Chebyshev II型M程序如下：Fs=20000; %抽样频率20KHzFlp=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

20、,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=2b= 0.2357 0.4241 0.2357a= 1.0000 -0.2996 0.19503.3 结果分析从频率响应图中可以看出:巴特沃斯滤波器具有单调下降的幅频特性，通带内平滑;切比雪I型滤波器的幅频特性在通带内有波动，

21、阻带内单调;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); %画出

22、原始语音信号的波形 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=44100bits=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*ta

23、n(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); s