脉冲响应不变法设计数字低通滤波器

1、燕山大学课程设计说明书题目脉冲响应不变法设计数字低通滤波器学院（系）：电气工程学院年级专业:10级精密仪器及机械2班学 号: 学生姓名:指导教师:刘永红教师职称:讲师电气工程学院课程设计任务书课程名称:数字信号处理课程设计13学号学生姓名1（专业）班级设计题目7、脉冲响应不变法设计数字低通滤波器设 计 技 术 参 数给定技术指标为：fp 100Hz，fs 300Hz， p 3dB， s 20dB， 采样频率Fs 1000 Hz。设 计 要 求设计Butterworth低通滤波器，用脉冲响应不变法转换成数字滤波器。参 考 资 料数字信号处理方面资料MATLAB方面资料周次前半周后半周应 完 成

2、内 容收集消化资料、学习 MATLAB软件， 进行相关参数计算编写仿真程序、调试指导教师签字基层教学单位主任签字基层教学单位：仪器科学与工程系指导教师:目录第一章 前言 第二章 数字信号处理基础知识部分2.1 巴特沃斯滤波器的设计2.2 数字滤波器IIR第三章 matlab 部分基础知识14143.1 MATLAB 介绍3.2 巴特沃斯 matlab 函数介绍15第四章 仿真部分及仿真图164.1 数字滤波器的设计步骤164.2 MATLAB 程序174.3 仿真结果18第五章 结论19第六章 参考文献20第一章前言数字信号处理课程设计是在学生完成数字信号处理和MATLAB 的结合后的基本实验

3、以后开设的。 本课程设计的目的是 为了让学生综合数字信号处理和 MATLAB 并实现一个较为完整 的小型滤波系统。滤波器是自动控制、信号处理和通信领域的重要组成部分 广泛地应用于各种系统中。 MATLAB 语言是一种简单、 高效的高 级语言,是一种内容丰富、功能强大的分析工具，其应用范围几乎 覆盖了所有的科学和工程计算领域。 MATLAB 中提供了丰富的用 于模拟滤波器设计的函数 ,通过编程可以很容易实现低通、高通、带通、带阻滤波器 ,并能画出滤波器的幅频特性曲线，大大简化了模拟滤波器设计。本文将通过利用 MATLAB 滤波器设计函数直 接实现巴特沃斯低通滤波器的设计，找到应用 MATLAB

4、来设计 低通滤波器的方法。介绍了用脉冲响应不变法设计巴特沃斯低通 滤波器的基本理论和设计思想，给出了基于 MATLAB 设计巴特 沃斯低通滤波器的具体步骤，文中还给出了利用 MATLAB 产生 的滤波器的幅频及相频特性曲线。详细介绍了在基于 MATLAB设计低通滤波器过程中常用到的工具和命令。第二章数字信号处理基础知识部分2.1巴特沃斯滤波器的设计1 .巴特沃斯滤波器的幅度平方函数及其特点巴特沃斯模拟滤波器幅度平方函数的形式是Ha j21 (j /j c)2N 11(/ c)2N(2-1-1 )中N为整数，是滤波器的阶次。=0 时，| Ha(j )=1 时；当时，1 Ha(j c)l=1，所以

5、称c为半功率点。此时通带内衰减20lgHa(j c) 3dB，所以 c又被称为3dB截止频率。巴特沃斯滤波器的振幅特性如图2.1.1所示。由图可见，当 时，随 加大，幅度迅速下降。下降的速度与阶次 N有关，N越大，幅度下降速度越快，过渡带越窄。Q2.幅度平方函数的极点分布及H a（S）的构成2将幅度平方函数Ha（j ）写成S的函数Ha(S)Ha( S)1 (S/J c)2N(2-1-2)此式表明幅度平方函数有2N个极点,极点sk用下式表示Z 八1/2N ZSk( 1) (J2N-1 (2-1-3 )c)2k 1eJ (E*eJ2J (丄竺2)ce 2 2N k=0,1,2，这2N个极点分布在S

6、平面半径为c的圆（称为巴特沃斯圆）上,角度间隔是n /N弧度。N=3时，极点间隔为n /3弧度或60度，如图2.1.2所示。图2.1.2三阶巴特沃斯滤波器极点分布极点对虚轴是对称的，且不会落在虚轴上。当 N是奇数时，实轴 上有极点；当N为偶数时，则实轴上没有极点。巴特沃斯滤波器 的N个极点为SkcJ冇）；k 0,1,，N 1（ 2-1-4）则Ha（S）的表达式即滤波器的系统函数为NHa（s）亍一（2-1-5）（S Sk）k 03. 频率归一化问题式（2-1-5 ）即为所求滤波器的系统函数，可看出Ha（s）与c有关，即使滤波器的幅度衰减特性相同，只要c不同，Ha（S）就不一样。为使设计统一，可将

7、所以的频率归一化。这里采用对3dB截止频率 c见归一化，归一化后的Ha（s）表示为1Ha（S）= （丄k 0 c（2-1-6），令频率，则S）c，称为归一化j，记为P称归一化复变量，这样归一化巴特沃斯的传输函数为a（P）1（2-1-7 ）式中 Pk 为（P Pk）k 0-/2k 1归一化极点，Pk S/c=e样，只要根据技术指标求出阶数 N，便可按式（2-1-8）求出N个归一化极点，再按照式（2-1-7）得到归一化的传输函数，k 0,1, N 1（ 2-1-8）。这a（p）。4.阶数N的确定阶数N的大小主要影响幅度特性下降速度，它应该由和s确定。令spp，,则N由下式表示:PJspN dlgs

8、p(2-1-9)关于3dB截止频率C，如果技术指标未给出，则按下式求出:p(100'1 p 1)1丽(2-1-10 )。5.设计步骤1）根据技术指标总结以上讨论，低通巴特沃斯滤波器的设计步骤如下:p, s, p和 s，用式（2-1-9，求出 No2）按照（2-1-8），求出归一化极点pk，代入（2-1-7 ），得到归 一化传输函数Ha（ P）。也可以直接查表.3）将Ha（ P）去归一化。将p=S/ C代入Ha（P），得到实际的滤波器传输函数Ha（s）。如果技术指标没有给出3dB截止频率c，可 以按照式（2-1-10 ,求出。2.2 数字滤波器IIR1.IIR 滤波器基础IIR 滤波器的

9、设计就是在给定了滤波器的技术指标后，确定滤波 器的阶数N和系数a,bi o在设计IIR滤波器时，常用的方法是利用模拟滤波器来设计数字滤波器，设计方法是首先将数字滤波 器技术指标转换为对应的模拟滤波器技术指标，然后设计满足技 术指标的模拟滤波器HaS，最后将设计出的模拟滤波器Ha(S)转 换为满足技术指标的数字滤波器 H z。将系统函数H(Z)从s平面转换到 z 平面的方法很多，主要有两种：一种是使数字滤波器的h n近似于模拟滤波器的ha t，可导出脉冲响应不变法；另一种是使数字滤波器的差分方程近似于滤波器的微分方程，由此可 导出双线性变换法。为保证转换后的 H z 稳定且满足技术要求， 对转换

10、关系提出以下 两点要求：(1)因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的。 模拟滤波器因果稳定要求其系统函数 Ha S的极点全部位于S平面的左半平面； 数字滤波器因果稳定则要求 H z 的极点全部在单 位圆内。因此，转换关系应是 S平面的左半平面映射到Z平面的 单位圆内。(2)数字滤波器的频率响应模仿模拟滤波器频率响应，S平面的虚轴映射Z平面单位圆。2.脉冲响应不变法脉冲响应不变法的核心是通过对连续函数ha(t)等间隔采样得到离散序列ha nT，使h(n) ha(nT)(其中T为采样间隔)，因此脉冲响应不变法是一种时域上的转换方法，转换步骤如下f 、拉氏逆变换.仕、等间隔采样 .Z.

11、 Z 、z变换./ 、Ha(S)ha(t)ha(nT)h(n)h(z)设模拟滤波器Ha(S)只有单阶极点，且分母多项式的阶次高于分子多项式的阶次，将Ha(S)用部分分式表示，则N AHa(S)(2-2-1)i 1 SSi式中，Si为Ha(S)的单阶极点。将Ha(S)进行拉氏逆变换得到ha(t)NAeSitu(t)(2-2-2)i式中，u(t)单阶阶跃函数。对ha(t)进行等间隔采样，采样间隔为T,得到Nh(n) ha(nT)Ae Tu(nt)i 1(2-2-3)对上式进行Z变换，得到数字滤波器的系统函数N AH(z) 彳鳥1i 11 e z(2-2-4)由这一转换过程看出，它对部分分式表达的模

12、拟系统函数更为方便，对任一极点S，Ha(s)到H(Z)得转换可直接用下式来完成丄1 AT 1(2-2-5)S si1 e z将式(2-2-1)的Ha S和式(2-2-4)的Hz加以比较，可以看出:1) S平面的单极点S Si变换到z平面上z eST处的单极点。2)Ha(S)与H(z)的部分分式的系数是相同的，都是 A 03)如果模拟滤波器是稳定的，所有极点Si位于s平面的左半平面,及极点的实部小于零，则变换后的数字滤波器的全部极点在单位 圆内，即模小于1，因此数字滤波器也是稳定的。3. H(ej )与Ha(j )的关系S平面和z平面 这里以采样下面分析从模拟滤波器转换到数字滤波器， 之间的映射

13、关系，从而找到这种转换方法的优缺点。 信号ha(t )作为桥梁，推导其映射关系。将ha(t )的米样信号用ha(t)表示,忆ha(t) (tnnT)对ha(t)进行拉式逆变换，得到H?a(s)f?(t)estdtha(t) (t nT)nst, e dtha(nT)e snTn式中，ha(nT)是ha(t)在采样点t=nT时的幅度值,它与序列的幅度相等，即h(n) ha(nT)，因此得到Ha(s)h(n)e snT h(n)z nnn上式说明采样信号的拉氏变换与相应的z esTH z esT (2-2-6)Z变换之间的映射关系为esT(2-2-7)已知模拟信号ha(t)的傅立叶变换Ha(j)和

14、其采样信号ha(t)的傅立叶变换Ha(j )之间的关系满足1Ha(j )- Ha(jjk s)(2-2-8)T k其中,s 2 /T，将s=j带入上式，得1Ha(s)- Ha(s jk s)(2-2-9)T k1H(z) Ha(s jk s)(2-2-10)T k上式表明将模拟信号ha(t)的拉氏变换在s平面上沿虚轴按照周期s 2 /T延拓后,再按上式的映射关系映射到 Z平面上，就得26到H(z)。 z esT可称为标准映射关系。c-jrlTR屮则得到:那么图2.2.1冲击响应不变法的映射关系jwre(2-2-11)0,r0,r上面关系式说明,s平面左半平面映射 z平面单位圆内，s平面的0,r

15、虚轴映射z平面的单位圆上，s平面右半平面映射z平面单位圆 外。这说明：第一，如果Ha(s)因果稳定，转换后得到H(z)仍因果稳定；第二，数字滤波器频率响应H(ejw)模仿模拟滤波器的频率响应Ha(j )，满足转换关系的两点要求。ej与H a (j )二者关系可由式(2-2-10)导出，考虑sejjk sHa(j9 kT-)(:-2-12)。这s为周就是说，数字滤波器的频率响应是模拟滤波器频率响应以期的周期延拓。因而正如采样定理所讨论的，只有当模拟滤波器 的频率响应是有限的，且带限于折叠频率以内时，即Ha0,-，才能使数字滤波器的频率响应在折叠T 2频率以内重现模拟滤波器的频率响应，而不产生混叠

16、失真，即H(ej )¥Ha(j)，4. 频率混叠现象实际上，任何一个模拟滤波器的频率响应都不是严格带限的，变 换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真。如果原模拟信号ha t的频带不是限于+/- n /T之间，则会在他的奇数倍附近产生频率混叠，从而映射到z平面上，在3 =+/-n附近产生频率混叠。脉冲响应不变法的频率混叠现象如下图所 示5. 脉冲响应不变法的优点脉冲响应不变法的优点是频率坐标变换是线性的，即一个优点是数字滤波器的单位脉冲响应完全模仿模拟滤波器的单 位脉冲响应，时域特性逼近好。第三章 matlab部分基础知识3.1 MATLAB 介绍是美国MathW

17、orks公司生产交互式大型软件，是一个可设计平台 MATLAB7.1 MATLAB的一个为科学和工程计算专门设计的以完成各种精确计算和数据处理的、可视化的、强大 的计算工具。 它集图示和精确计算于一身，在应用数学、物理、化工、机 电工程、医学、金融和其他需要进行复杂数值计算的领域得到了 广泛应用。它不仅是一个 在各类工程设计中便于使用的计算工 具，而且也是一个在数学、数值分析和工程 计算等课程教学中的 优秀的数学工具， 在世界各地的高和大型计算机上运行， 适 用于 Windows、 UNIX 等多种系统平台。 MATLAB 作为一种科学计 算的高级语言之所以受欢迎， 就是因为它有丰富的 函数资

18、源和 工具箱资源，编程人员可以根据自己的需要选择函数，而无需再 去编 写大量繁琐的程序代码，从而减轻了编程人员的工作负担， 被称为第四代编程语 言的 MATLAB 最大的特点就是简洁开放 的程序代码和直观实用的开发环境3.2 巴特沃斯 matlab 函数介绍1. buttordN,wc=buttord(wp , ws,a p,a s)用于计算巴特沃斯数字滤波 器的阶数N和3dB截止频率WC。调用参数wp，ws分别为数字 滤波器的通带、阻带截止频率的归一化值，要求：0Wwp< 1 , 0< ws< 1且对于低通滤波器 wpvws。1表示数字频率pi。2. butterb， a=

19、butter (N， wc，ftype')计算 N 阶巴特沃斯数字滤波器系统函数分子、 分母多项式的系数 向量 b、 a。调用参数 N 和 wc 分别为巴特沃斯数字滤波器的阶数和 3dB 截止 频率的归一化值(关于 pi 归一化)，一般是调用 N,wc=buttord(wp ， ws，a p，a s)格式计算N和wc。系数b、a是按照z-1的升幕排 列。3freqz范围内的 N 个H, freqz B, A,N,Fs 计算离散系统在 0频率等分点的频率响应的值。4. absAbs ( H)为取H得幅度。5. angleAngle ( H)为取H的相位角，单位为弧度，用pha显示计算的相

20、 位。6. axisaxis(xmin xmax ymin ymax) 中分别给出 x 轴和 y 轴的最大值、 最小值7. subplotsubplot(m,n,p) 将多个图画到一个平面上的工具。 其中， m 表示是图排成m行，n表示图排成n列，P是指你现在要把曲线画到figure 中哪个图上。8. plotplot(x,y)是绘制二维图形的最基本函数，以x元素为横坐标值，y元素为纵坐标值绘制曲线。9. Xlabel x 轴注解Ylabel y轴注解Title 标题注解Grid on 开启网格格式第四章仿真部分及仿真图4.1数字滤波器的设计步骤设计流程图4.2 MATLAB 程序Fs=100

21、0; fp=100; fs=300; Ap=3; As=20;Wp=fp/(Fs/2);Ws=fs/(Fs/2);N,Wc=buttord(Wp,Ws,Ap,As);b,a=butter(N,Wc);H,F=freqz(b,a,500,Fs);subplot(2,2,2) plot(F,20*log10(abs(H) xlabel(' 频率 (Hz)'); ylabel(' 幅度 (dB)') axis(0 500 -40 3);title(' 数字滤波器分频曲线 ') grid on subplot(2,2,1) plot(F,abs(H);xlabel(' 频率 (Hz)');ylabel(' 幅度 ');title(' 数字滤波器幅频曲线 ') grid on;subplot(2,2,3) pha=angle(H);plot(F,pha); xlabel(' 频率 (Hz)'); ylabel(' 相位 (rad)') grid ontitleC数字低通滤波器相频曲线')4.3仿真结果06S ncu01.2a、ifK II;I uIIIII：KI：IIIqIIII:II:;IIL' xi、IJ10C 2W3 CO4M 测1020