MATLAB在数字滤波器设计教学中的应用

1、MATLAB在数字滤波器设计教学中的应用自从 1965 年库利和图基发表了快速傅里叶变换算法以来， 数字信号处理这一学科蓬勃发展， 逐渐成为一门涉及许多学科而 广泛引用于许多领域的新兴学科。 数字滤波器具有精度高、 稳定 性好、灵活性强、便于集成、高性能等优点 1-2 ，在数字信号 处理中发挥着十分重要的作用。MATLAB件是美国MathWorks公司推出的一套高性能科学 计算软件，它具有强大的数值分析、矩阵运算、信号处理和图形 显示的功能。 信号处理工具箱中包含许多由信号处理领域的权威 专家编写的函数，可以直接调用，使编程变得很简单。特别在数 字滤波器的设计中，MATAL软件能够很方便的完成

2、数值计算以 及图形的绘制。下面我们具体说明 MATLAB件分别在IIR ( Infinite Impulse Response )和 FIR(Finite Impulse Response)数字滤波器设计中的应用。1 IIR 数字滤波器的设计IIR 数字滤波器的系统函数包含零点和单位圆内的极点，因 此可以用较低的阶数设计出频率选择性较高的滤波器， 所用的存 储单元少，计算量小、效率高。但 IIR 滤波器一般不具有线性相 位。设计 IIR 数字滤波器的思路是先把数字技术指标转换成模拟 技术指标， 然后设计模拟 IIR 滤波器， 最后映射成一个等效的数 字滤波器。我们用冲击响应不变法和双线性变换法

3、分别设计巴特沃斯滤波器和切比雪夫滤波器为例说明MATLAB件的运用方法。设计IIR数字带通滤波器，给定指标为（1）200Hz600Hz衰减20dB,（ 3）抽样频率 fs=2kHz。先将模拟信号经过AD转换得到数字信号，那么所设计数字 滤波器的的指标为：其中 Wc1=200 Hz Wc2=400 Hz为两个通带截止频率； Wst1=100 Hz，Wst2=600 Hz 为两个阻带截止频率，fs=2000 Hz为采样频率。 wc1， wc2， wst1 ， wst2 分别为转换得到的数字滤 波器的两个通带截止频率和两个阻带截止频率。通常的设计方法是将数字滤波器的技术指标转换成模拟滤 波器的技术指

4、标， 再将模拟带通滤波器的技术指标转换成模拟低 通滤波器的技术指标， 然后设计模拟低通滤波器， 最后将设计的 模拟低通滤波器转换成数字带通滤波器。 设计方法非常复杂， 然 而运用MATLAB件来设计非常方便，函数 buttord 和cheblord 直接根据给定的技术指标返回巴特沃斯滤波器和切比雪夫滤波 器的阶数和3dB的截止频率，再用butter和cheb1ord函数返回 所设计滤波器的系统函数。 用冲击响应不变法设计， 先设计出模 拟滤波器， 再用 impinvar 函数将模拟滤波器转换为数字滤波器； 用双线性变换法设计可以用该函数直接设计数字滤波器。 具体程 序代码如下：fs=2000

5、； %采样频率。Wc=2*pi*200 2*pi*400；Wst=2*pi*100 2*pi*600； %通带截止频率和阻带截止频率。Rp=2 Rst=20 ； %通带最大衰减和阻带最小衰减。wc=Wc/fs； wst=Wst/fs ； %所设计的数字滤波器的截止频率。N，Wn=buttord （Wc，Wst，Rp，Rst，s ）； %用巴特沃 斯滤波器设计，冲击响应不变法， s 表示设计模拟滤波器， N 为阶数，Wn为3dB的截止频率。B， A=butter （N， Wn， s ）； %返回模拟滤波器的系统函 数，B为分子多项式系数，A为分母多项式系数。b ， a=impinvar （B，

6、A， fs ）； %函数 impinvar 用冲激响 应不变法将模拟滤波器转化为数字滤波器， 返回 b 和 a 分别为数 字滤波器的分子和分母多项式系数。h1 ， w1=freqz （b， a， 100）； %函数 freqz 算出冲激响应 不变法设计巴特沃斯数字滤波器的频率响应。N， Wn=buttord （wc/pi ， wst/pi ， Rp， Rst）； %用巴特沃 斯滤波器直接设计数字滤波器， 得到滤波器的阶数和3dB截止频 率。b ， a=butter （N， Wn）； %返回数字滤波器的系数。h2 ， w2=freqz （b， a， 100）； %函数 freqz 算出双线性变

7、换法设计巴特沃斯数字滤波器的频率响应。N， Wn=cheb1ord（Wc， Wst， Rp， Rst， s ）； %用切比雪 夫滤波器设计模拟滤波器，得到滤波器的阶数和3dB截止频率。B，A=cheby1 (N，Rp，Wn，s )；%得到模拟滤波器的系 统函数。b ， a=impinvar ( B， A，fs )； %函数 impinvar 用冲激响 应不变法将模拟滤波器转化为数字滤波器h3 ，w3=freqz (b，a，100)； %函数 freqz 算出冲击响应 不变法设计切比雪夫数字滤波器的频率响应。N，Wn=cheb1ord(wc/pi ，wst/pi ，Rp，Rst)； %用切比 雪

8、夫滤波器直接设计数字滤波器， 得到滤波器的阶数和 3dB 截止 频率。b , a=cheby1 (N, Rp, Wn ;俺回数字滤波器的系数。h4， w4=freqz (b， a， 100)； %函数 freqz 算出双线性变 换法设计切比雪夫数字滤波器的频率响应。x=wc/pi , wst/pi ；y=-Rp , -Rp, -Rst , -Rst ； %技术指标plot(w1/pi , 20*log10(abs(h1), + , w2/pi, 20*log10(abs (h2), A , w3/pi , 20*log10 (abs (h3), v, w4/pi , 20*log10 (abs

9、(h4), o , x, y, * )； %画图 grid ； xlabel (w/pi )； ylabel (20Log(|H| ) )； axis (0 , 1, -50, 10)；legend ( 巴特沃斯,冲击响应不变法 , 巴特沃斯,双线 性变换法 , 切比雪夫,冲击响应不变法 , 切比雪夫,双线性 变换法 )；图 1 IIR 滤波器的幅频响应曲线图 1 为程序运行后的滤波器的幅频响应曲线。 可以看出所设 计的4种数字滤波器通带截止频率为0.2 n和0.4 n , w0.6 n为阻带，阻带的衰减均大于 20dB,满足了数字滤波器设计的技术 要求。巴特沃斯滤波器在通带没有出现波纹， 而

10、切贝雪夫滤波器 在通带出现波纹； 双线性变换法实现模拟滤波器到数字滤波器的 映射不出现混叠现象， 而冲击响应不变法实现该转换高频部分有 混叠现象，因此，在高频阻带，冲击响应不变法设计的滤波器的 幅频响应曲线会往上翘，所设计的滤波器的特征都符合要求。2 FIR 数字滤波器的设计FIR 滤波器的系统函数没有极点，保证了 FIR 滤波器是稳定 的。FIR滤波器很容易具有线性相位。设计 FIR滤波器通常有窗 函数法和频率采样法。窗函数法设计 FIR 滤波器的基本原理就是在时域逼近理想 滤波器的单位冲击响应 3 ，从而使所得到的频率响应与所要求 的理想频率响应尽可能接近。 这种时域逼近时通过在时域加窗的

11、 方法得到的， 通过给定的技术指标选择窗函数的类型和长度。 例 如利用窗函数法设计一个线性相位低通 FIR 数字滤波器， 通带截 止频率wp=0.4 n，阻带截止频率wst=0.6 n，阻带最小衰减50dB。 常用的窗函数有矩形窗、三角窗、汉宁窗、海明窗、布莱克曼窗 和凯泽窗，每个窗函数加窗后的滤波器的阻带最小衰减不一样， 根据阻带最小衰减50dB我们选择海明窗，然后根据海明窗过渡 带的宽度为6.6 n /N ,结合需要设计的过渡带宽度 wst-wp=0.2 n 计算出海明窗的长度 N。截止频率取通带截止频率和阻带截止频 率中心，通过 fir1 函数得到 FIR 滤波器的系统函数，程序代码 如

12、下：wp=0.4*pi ；wst=0.6*pi ；wc=（ wp+wst）/2 ； %wc为截止频率。N=ceil （ 6.6*pi/（ wst-wp）; %函数 ceil （ x）计算不小于x的最小整数，N为窗函数的长度。b=fir1 （ N-1， wc/pi ， hamming（ N）; %函数 fir1 由滤波 器的阶数， 截止频率以及窗函数的类型得到 FIR 滤波器的系统函 数。h ， w=freqz （ b， 1， 500）; %函数 freqz 由系统函数的系 数返回频率响应函数 h（ w）。plot （ w/pi ， 20*log10 （ abs（ h）， - ）; ylabel

13、 （ 20log （|H| ） ）; xlabel （ w/pi ）; grid ; %画出滤波器的幅频响 应对数图。axis （0 1 -100 10）;图 2 窗函数法设计 FIR 数字滤波器的幅频响应曲线程序运行后，结果如图2所示。滤波器的长度 N=33,阻带 截止频率0.6 n的衰减为50dB,达到了滤波器的设计指标。频率采样法对期望的频率响应在02n之间等间隔地采样，通过采样点的值还原出系统的频率响应。 虽然频率采样点与理想的频率响应精确的符合，但是无法控制采样点之间如何插 值。抽样点之间的理想频率特性变化越陡峭， 则内插值与理想值 之间的误差就越大。 在过渡带两边会产生肩峰， 通带

14、和阻带也会 产生波纹。 这可以通过增加过渡带抽样点进行优化， 例如设计一 个低通FIR数字滤波器，wc=0.5 n,抽样点数为N=33,要求滤 波器具有线性相位。（ 1）增加一点过渡带进行优化，过渡带抽 样点值为0.5 ;（2）增加抽样点数至N=65,加两点过渡带优化，过渡带抽样点值为 0.5886， 0.1065 。程序代码如下：N=33；absH=ones（ 1， 9）， zeros （ 1， 16）， ones（ 1， 8） ； % 理想低通滤波器的频率抽样点的幅度。n=0 ：N-1 ； angH=-2*pi/N*n*（N-1） /2； %理想低通滤波器的频率抽样点的相角。H=absH.

15、*exp（j*angH ）； %理想低通滤波器的频率抽样点。 b=ifft（ H， N）； %所设计滤波器的单位抽样响应。H1， w=freqz （ b， 1， 500）； %加一点过渡带优化 |H（ 9） |=0.5N=33；absH=ones（ 1， 9）， zeros （ 1， 16）， ones（ 1， 8） ； % 理想低通滤波器的频率抽样点的幅度。absH （10） =0.5； absH （ 25） =0.5；n=0 ：N-1 ；angH=-2*pi/N*n* （N-1）/2 ；%理想低通滤波器的频率抽样 点的相角。H=absH.*exp（j*angH ）； %理想低通滤波器的频率

16、抽样点。 b=ifft （ H， N）； %所设计滤波器的单位抽样响应。H2， w=freqz （ b，1， 500）；%增加抽样点数至N=65,加两点过渡带优化N=65； absH=ones（1，17），zeros（1，32），ones（1，16） ；% 理想低通滤波器的频率抽样点的幅度。absH（18）=0.5886；absH（19）=0.1065；absH（48）=0.1065；absH（49） =0.5886 ；n=0：N-1 ；angH=-2*pi/N*n* （N-1）/2 ； %理想低通滤波器的频率抽样 点的相角。H=absH.*exp（j*angH ）； %理想低通滤波器的频率抽

17、样点。 b=ifft （ H， N）； %所设计滤波器的单位抽样响应。 H3，w=freqz （b，1，500）；plot （ w/pi ，20*log10 （ abs（H1），+ ，w/pi ，20*log10 （abs（H2）， o ，w/pi ，20*log10 （ abs（ H3）， v ）； ylabel （ 20log|H| ）； xlabel （w/pi ）； grid ；lege nd （N=33 , N=33 优化,N=65 优化）;axis (0 1 -100 10 )；图 3 频率采样法设计 FIR 数字滤波器的幅频响应曲线 程序运行后，结果如图 3 所示。可以看出，增加一点过渡带 抽样后，过渡带的宽度由原来的 2n /33，增大为4 n /33，阻带 最小衰减由-17dB增大为-30dB。可见增加过渡带抽样点，阻带 衰减增大，但