基于MATLAB的数字滤波器设计

本科 毕业设计 机械毕业设计下载球磨机的毕业设计下载关于网络爬虫的毕业设计下载关于网络爬虫的毕业设计下载河南城建学院毕业设计论文下载 (论文) 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 目: 基于MATLAB的数字滤波器设计 姓 名: 学 号： 02452102 专 业: 电子信息工程 指导教师: 职 称: 讲 师 2006年 月 日 摘 要 本文分别研究了在MATLAB环境下IIR数字滤波器和FIR数字滤波器的设计方法及实现方法，并进行图形用户界面设计，以显示本文所介绍滤波器的特性。 在无限脉冲响应（IIR）数字滤波器设计中，先进行模拟滤波器的设计，然后进行模拟-​数字滤波器转换，即采用脉冲响应不变法及双线性Z变化法设计数字滤波器，最后进行滤波器的频带转换。在有限脉冲响应（FIR）数字滤波器设计中，讨论了FIR线性相位滤波器的特点和用窗函数法设计FIR滤波器两个问题。两类滤波器整个设计过程都是按照理论分析、编程设计、具体实现的步骤进行的。 为方便分析者直观，形象，方便的分析滤波器的特性，创新的设计出了图形用户界面——滤波器分析系统。整个系统分为两个界面，其内容主要包含四部分：System(系统)、Analysis(分析)、Tool(工具)、Help(帮助)。 关键词：数字滤波器、MATLAB、无限脉冲响应、有限脉冲响应、图形用户界面 Abstract In this thesis, designs of the Infinite Impulse Response digital filter (IIR) and Finite Impulse Response digital filter (FIR) under MATLAB are studied. And the Graphical User Interfaces (GUI) to analysis the characteristics of filter is designed. The design of IIR filters can be achieved through three steps: firstly, the design of analog low-pass filter; secondly, it is analog-to-digital filter conversion; lastly, it is the conversion of filter frequency band. In design of FIR filters, two questions are discussed: the characteristics of FIR linear phase filter and reasoning of related formulas; the other is about the design of the FIR filters by means of window functions. The design of FIR and IIR follows the procedures of theoretical analysis, programming design and realization. We design the Graphical User Interfaces (GUI) of the digital filter analysis system makes sure of people can analysis the characteristics of the design directly and easily. The whole system divide into two graphical interfaces, it contains four parts: System, Analysis, Tool and Help. Keywords: Digital Filter, MATLAB, IIR, FIR, Graphical User Interfaces 目 录 第1章 绪论………………………………………………………………………………1 1.1 数字滤波技术………………………………………………………………………1 1.1.1 滤波器原理………………………………………………………………………2 1.1.2 数字滤波器设计方法概述………………………………………………………2 1.2 MATLAB软件简介……………………………………………………………………3 第2章 IIR滤波器设计及其MATLAB实现………………………………………………5 2.1 模拟滤波器设计及其MATLAB实现…………………………………………………5 2.1.1 巴特沃斯低通滤波器设计………………………………………………………5 2.1.2 切比雪夫低通滤波器设计………………………………………………………7 2.1.3 椭圆低通滤波器设计……………………………………………………………11 2.2 模拟—数字滤波器变换及其MATLAB实现………………………………………13 2.2.1 脉冲响应不变法…………………………………………………………………13 2.2.2 双线性Z变换法…………………………………………………………………17 2.3 频带变换及其MATLAB实现………………………………………………………19 2.3.1 模拟低通滤波器转换成数字高通滤波器………………………………………19 2.3.2 模拟低通滤波器转换成数字带通滤波器………………………………………21 2.3.3 模拟低通滤波器转换成数字带阻滤波器………………………………………23 2.4 小结…………………………………………………………………………………25 第3章 FIR滤波器设计及其MATLAB实现……………………………………………26 3.1 线性相位FIR数字滤波器的条件和特点…………………………………………26 3.1.1 线性相位条件……………………………………………………………………26 3.1.2 线性相位FIR滤波器幅度特性的特点…………………………………………27 3.1.3 线性相位FIR滤波器零点分布特点……………………………………………27 3.2 常用窗函数及其MATLAB实现……………………………………………………27 3.2.1 常用窗函数介绍…………………………………………………………………27 3.2.2 各种窗函数的实现与比较………………………………………………………29 3.3 基于窗函数的FIR数字滤波器设计及其MATLAB实现……………………………30 3.3.1 海明窗设计数字低通滤波器…………………………………………………31 3.3.2 汉宁窗设计数字高通滤波器…………………………………………………33 3.3.3 布拉克曼窗及三角窗设计数字带通滤波器…………………………………35 3.3.4 汉宁窗设计数字带阻滤波器………………………………………………37 3.4 小结…………………………………………………………………………………38 第4章 图形用户界面设计--滤波器分析系统………………………………………39 4.1 MATLAB中GUI设计技术特点……………………………………………………39 4.1.1 GUIDE简介………………………………………………………………………39 4.1.2 GUI设计规范及设计步骤………………………………………………………40 4.2 滤波器分析系统及其MATLAB实现………………………………………………40 4.2.1 滤波器分析系统的初始界面……………………………………………………41 4.2.2 滤波器分析系统的主界面………………………………………………………43 4.3 小结…………………………………………………………………………………48 第5章 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf ……………………………………………………………………………49 参考文献…………………………………………………………………………………50 致谢………………………………………………………………………………………51 第1章 绪论 1.1 数字滤波技术 数字滤波是数字信号处理的重要基础，数字信号处理主要是研究用数字或符号的序列来表示信号波形，并用数字的方式去处理这些序列，把它们改变成在某种意义上更希望的形式，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。数字滤波器在对信号的过滤、检测与参数估计等处理过程中，是使用最为广泛的一种线性系统。 数字滤波器（Digital Filter，简称为DF）是指完成信号滤波处理功能的、用有限精度算法实现的离散时间线性非时变系统。数字滤波器的数学运算通常有两种实现方式。一种是频域法，即利用FFT快速运算办法对输入信号进行离散傅立叶变换，分析其频谱，然后根据希望的频率特性进行滤波，再利用傅立叶反变换恢复出时域信号。这种方法具有较好的频域选择特性和灵活性，并且由于信号频率与所希望的频谱特性是简单的相乘关系，所以它比计算等价的时域卷积要快得多。另一种方法是时域法，这种方法是通过对离散抽样数据作差分数学运算来达到滤波目的的。 数字滤波器的输入是一组（由模拟信号取样和量化的）数字量，其输出是经过数字变换的另一组数字量。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。随着数字技术的发展，用数字技术实现滤波器的功能愈来愈受到人们的重视，并得到了广泛的应用。 数字信号处理学科的一项重大进展是关于数字滤波器设计方法的研究。60年代中期，开始形成关于数字滤波器的一套完整的正规理论。这一时期，提出了各种各样的数字滤波器结构，有的以运算误差最小为特点，有的则以运算速度高见长；出现了数字滤波器的各种逼近方法和实现方法，对递归和非递归两类滤波器作了全面的比较；统一了数字滤波器的基本概念和理论，对有限冲激响应（IIR）和无限冲激响应（FIR）的认识有了完整理论。70年代后，科学技术蓬勃发展，数字信号处理开始与大规模和超大规模集成电路技术、微处理技术等新工艺新技术结合起来，并引进计算机辅助设计方法，大大丰富了数字滤波器的分析与设计，各种新的数字信号处理系统，也都能用专用数字硬件实时加以实现。相信在未来，随着电子仪器与电子技术应用系统朝着数字化、小型化、自动化以及多功能化等方向发展，包括数字滤波器在内的数字信号处理技术会有以惊人的速度进行飞跃式发展。 1.1.1 滤波器原理 滤波器，顾名思义，其作用是对输入信号起到滤波的作用。对于图1-1所示的LSI系统，其时域输入输出关系为： (1-1) x(n) y(n) 图1-1 若y(n),x(n)的傅立叶变换存在，则输入输出的频域关系是： (1-2) 当输入信号x(n) 通过滤波器系统h(n) 后，其输出 y(n) 中不再含有|ω|>ωc的频率成分，仅使|ω|<ωc的信号成分通过。因此，滤波器的形状不同，其滤波后的信号结果也不一样。 若滤波器的输入、输出都是离散时间信号，那么该滤波器的单位冲激响应h(n)也必然是离散的，这种滤波器称为数字滤波器（DF，Digital Filter）。当用硬件实现一个DF时，所需元件是延迟器、乘法器和加法器；而利用MATLAB软件时，它仅需线性卷积程序便可实现。而模拟滤波器（AF，Analog Filter）只能用硬件实现。因此DF比AF容易实现，且更容易获得理想的滤波性能。数字滤波器还具有以下优点：精度和稳定性高；系统函数容易改变；灵活性高；不存在阻抗匹配问题；便于大规模集成；可实现多维滤波。 1.1.2 数字滤波器设计方法概述 数字滤波器从功能上分类：可以分为低通滤波器(LP，Low Pass)、高通滤波器(HP，High Pass)、带通滤波器(BP，Band Pass)、带阻滤波器(BS，Band Stop)。 从滤波器的网络结构或者从单位脉冲响应分类：如同模拟滤波器的性能可由g(t)和G(s)来表征一样，数字滤波器的性能完全取决于h(n)和H(z)。因此，数字滤波器可以按照单位取样响应（或称脉冲响应，冲激响应等）h(n)的性质分为两类：有限脉冲响应(Finite Impulse Response)数字滤波器，简称FIR数字滤波器，它的h(n)序列长度是有限的；无限脉冲响应(Infinite Impulse Response)数字滤波器，简称IIR数字滤波器，它的h(n)序列长度是无限的，即当 时，h(n)仍有效。 数字滤波器的设计与实现，通常按下述步骤进行： 1、根据不同用途提出数字滤波器的技术指标、性能要求。 2、设计一个稳定的、因果的数学模型H(z)来逼近所要求的技术指标，并用有限精度的运算实现所设计的系统。本文将重点介绍此步骤。 3、设计专用的数字硬件来实现这个数字模型，或者通过电子计算机运行程序软件予以实现。本文将侧重于软件实现。 IIR滤波器和FIR滤波器的设计方法很不相同： IIR滤波器设计方法有两类，经常用到的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计思路是：先设计模拟滤波器得到传输函数G(s)，然后将G(s)按某种方法转换为数字滤波器的系统函数H(z)。这一类方法是基于模拟滤波器的设计方法相对比较成熟，它不仅有完整的设计公式，也有完整的图标供查阅，更可以直接调用MATLAB中的对应的函数进行设计。另一种是直接在频域或者时域中进行设计，设计时必须使用计算机辅助，直接调用MATLAB中的程序或函数即可设计。 FIR滤波器不能采用由模拟滤波器设计进行转换的方法，而经常使用的是窗函数法和频率采样法。也可以借助计算机辅助设计软件采用切比雪夫等波纹逼近法进行设计。 1.2 MATLAB软件简介 MATLAB是英文Matrix Laboratory（矩阵实验室）的缩写，它是由美国MathWorks公司推出的用于数值计算和图形处理的数学计算环境。在MATLAB环境下，用户可以集成的进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。它优秀的数值计算能力和卓越的数据可视化能力使其在同类软件中脱颖而出。 MATLAB系统最初由Cleve Moler用FORTRAN语言设计，现在的MATLAB程序是由MathWorks公司用C语言开发的。它的第一版（DOS版本1.0）发行于1984年，经过多年的改进，版本不断升级，其所包含的工具箱功能也越来越丰富，应用越广泛。 MATLAB语言之所以能如此迅速的普及，显示出如此旺盛的生命力，是由于它有着不同于其他语言的特点，归纳如下： 1、简单易学：MATLAB不仅是一个开发软件，也是一门编程语言。其语法规则与结构化高级编程语言(如C语言等)大同小异，而且使用更为简便。用MATLAB编程写程序犹如在演算纸上排列出公式与求解问题，所以被称为演算纸式科学算法语言。 2、计算功能强大：MATLAB拥有庞大的数学、统计及工程函数，可使用户立刻实现所需的强大数学计算功能。由各领域的专家学者们开发的数值计算程序，使用了安全、成熟、可靠的算法，从而保证了最快的运算速度和可靠的结果。此外，MATLAB还有数十个工具箱，可解决应用中的大多数数学、工程问题。 3、先进的可视化工具：MATLAB提供功能强大的、交互式的二维和三维绘图功能，可使用户创建富有表现力的彩色图形。可视化工具包括曲面渲染、线框图、伪彩图、光源、图像显示、动画等。 4、开放性、可扩展性强：M文件是可见的MATLAB程序，所以用户可以查看源代码。开放的系统设计使用户能够检查算法的正确性，修改已存在的函数，或者加入自己的新部件。 5、特殊应用工具箱：MATLAB的工具箱加强了对工程及科学中特殊应用的支持。工具箱也和MATLAB一样是完全用户化的，可扩展性强。将某个或几个工具箱与MATLAB联合使用，可以得到一个功能强大的计算组合包，满足用户的特殊要求。 MATLAB数字信号处理工具箱和滤波器设计工具箱专门应用于信号处理领域。工具箱提供了丰富而简便的设计，使原来繁琐的程序设计简化成函数的调用。只要以正确的指标参数调用相应的滤波器设计程序或工具箱函数，便可以得到正确的设计结果，使用非常方便。 接下来的两章我们将分别对IIR数字滤波器和FIR数字滤波器的设计实现进行分析讨论。在具体的设计过程中，将充分发挥MATLAB软件在数值计算、图像处理中的优势，并充分利用其强大的滤波器设计工具箱(Filter Design Toolbox)及图像用户界面工具辅助设计分析。 第2章 IIR滤波器设计及其MATLAB实现 IIR数字滤波器具有无限持续时间脉冲响应，而所有的模拟滤波器一般都具有无限长的脉冲响应，因此它与模拟滤波器相匹配。于是，设计IIR数字滤波器可以先设计模拟低通滤波器G(s)，然后按一定规则将G(s)进行滤波器变换而得到数字域的IIR低通滤波器H(z)，最后通过频带变换转换成其它频带滤波器。本章按按此思路将IIR数字滤波器设计分为三部分：设计模拟滤波器；进行滤波器转换，得到数字滤波器；最后再简单分析数字滤波器的频带转换问题。下面，将以此进行分别讨论。 2.1 模拟滤波器设计及其MATLAB实现 IIR数字滤波器技术设计最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法，它不但有完整的设计公式，而且还有较为完整的图表可供查询，因此充分利用这些已有的资源将会给数字滤波器的设计带来很大方便。滤波器设计最重要的是寻找一个稳定、因果的系统函数去逼近滤波器的技术指标，因此模拟滤波器的设计十分重要。实际中，有三种广泛应用的滤波器，即巴特沃斯低通滤波器、切比雪夫低通滤波器，椭圆低通滤波器，下面，将分别予以介绍。 2.1.1 巴特沃斯低通滤波器设计 巴特沃斯滤波器拥有最平滑的频率响应，在截止频率以外，频率响应单调下降。在通带中是理想的单位响应，在阻带中响应为0。在截止频率处有3dB的衰减。巴特沃斯滤波器除了具有平滑单调递减的频率响应的优点之外，其过渡带的陡峭程度正比于滤波器的阶数。高阶巴特沃斯滤波器的频率响应近似于理想低通滤波器。 巴特沃斯(Butterworth)低通滤波器是将巴特沃斯函数作为滤波器的传递函数 ，它的平方幅度函数为： (2-1) 式中, 为滤波器频率， 为3dB截止频率，N表示滤波器的阶次。 MATLAB信号处理工具箱为低通模拟巴特沃斯滤波器的产生提供了函数buttap,其调用的格式为：[z,p,k]=buttap(N)，其中，z表示零点,p表示极点,k表示增益，N表示阶次。下面举例设计巴特沃斯低通模拟滤波器： %巴特沃斯低通模拟滤波器 clear all; n=0:0.01:2; for i=1:4 switch i case 1 N=2; case 2 N=5; case 3 N=10; case 4 N=20; end [z,p,k]=buttap(N); %函数buttap--设计巴特沃斯低通滤波器 [b,a]=zp2tf(z,p,k); %函数zp2tf--零极点增益模型转换为传递函数模型 [H,w]=freqs(b,a,n); %函数freqs--求解模拟滤波器频率响应 magH2=(abs(H)).^2; %函数abs--取模值函数 hold on %函数hold--控制是否保持当前图形 plot(w,magH2) %函数plot--画二维线性图 axis([0 2 0 1]); %函数axis--控制坐标轴比例和外观 end xlabel('w/wc'); ylabel('|H(jw)|^2'); title('巴特沃斯低通模拟滤波器'); grid on 由此得出的巴特沃斯低通模拟滤波器的平方幅度响应曲线如图所示： 图2-1 巴特沃斯滤波器的幅频特性(基于MATLAB实现) 由图2-1可清楚地分析出，巴特沃斯滤波器拥有平滑的频率响应，在截止频率以外，频率响应单调下降。其过渡带的陡峭程度正比于滤波器的阶数，高阶巴特沃斯滤波器的频率响应近似于理想低通滤波器。由图2-1所示，对N=2阶滤波器和N=20阶滤波器的平方幅度响应进行比较后，证明了高阶巴特沃斯滤波器有着更好的幅度特性，更接近理想低通滤波器。 2.1.2 切比雪夫低通滤波器设计 巴特沃斯滤波器的频率特性在通带和阻带内部都是随频率单调变化的，如果在通带的边缘能够满足指标，那么在通带的内部会超过设计的指标要求，造成滤波器的阶次N比较高。因此，有效地设计方法是将指标的精确度均匀的分布在整个通带内，或者均匀分布在整个阻带内，或者同时分布在两者之内，这样就可用阶数较低的系统来满足要求。这可通过选择具有等波纹特性的逼近函数来达到。 切比雪夫滤波器的振幅特性就具有这种等波纹特性。它有两种形式：振幅特性在通带内是等波纹的、在阻带内是单调的切比雪夫I型滤波器；振幅特性在通带内是单调的、在阻带内是等波纹的切比雪夫II型滤波器。采用何种形式的切比雪夫滤波器取决于实际用途。 切比雪夫I型(Chebyshev-I)滤波器平方幅度响应函数为： (2-2) 其中 是一个小于1的正数，它与通带波纹有关， 越大，波纹也越大； 是的切比雪夫多项式，它的定义为： (2-3) 而切比雪夫II型(Chebyshev-II)滤波器平方幅度响应函数为： (2-4) 其中 是一个小于1的正数，表示波纹变化情况； 为截止频率；N为滤波器的阶次，也是 的阶次。 N为奇数 N为偶数 图2-2 理想切比雪夫I型滤波器的幅频特性 N为奇数 N为偶数 图2-3 理想切比雪夫II型滤波器的幅频特性 图2-2和图2-3分别画出了理想时的切比雪夫I型与切比雪夫II型滤波器阶次N为奇数与偶数时的幅频特性。而通过MATLAB信号处理工具箱中的函数cheb1ap及cheb2ap，可以实现切比雪夫滤波器设计，其调用格式为： [z,p,k]=cheb1ap(N,Rp) [z,p,k]=cheb2ap(N,Rs) 其中，z表示零点,p表示极点,k表示增益，N表示阶次，Rp为通带波纹(dB), Rs为阻带波纹(dB)。 %切比雪夫I型低通模拟滤波器 clear all; n=0:0.01:2; for i=1:2 switch i case 1 N=7; case 2 N=8; end Rp=1; %通带纹波系数为1dB [z,p,k]=cheb1ap(N,Rp); %函数cheb1--设计切比雪夫I型低通滤波器 [b,a]=zp2tf(z,p,k); [H,w]=freqs(b,a,n); %函数freqs--求解模拟滤波器频率响应 magH2=(abs(H)).^2; %绘制图形 posplot=['12' num2str(i)]; %输出图形 subplot(posplot) %函数subplot--创建子图图区 plot(w,magH2) axis([0 2 0 1]); xlabel('w/wc'); ylabel('|H(jw)|^2'); title(['N=' num2str(N)]); end 图2-4 切比雪夫I型滤波器的幅频特性(基于MATLAB实现) %切比雪夫II型低通模拟滤波器 clear all; n=0:0.01:2; for i=1:2 switch i case 1 N=7; case 2 N=8; end Rs=10; %阻带文波系数为10dB [z,p,k]=cheb2ap(N,Rs); %函数cheb2---设计切比雪夫II型低通滤波器 [b,a]=zp2tf(z,p,k); [H,w]=freqs(b,a,n); magH2=(abs(H)).^2; %输出图形 posplot=['12' num2str(i)]; subplot(posplot) plot(w,magH2) axis([0 2 0 1.1]); 图2-5 切比雪夫II型滤波器的幅频特性(基于MATLAB实现) 用图2-4、图2-5 与图2-2、2-3分别作比较可清楚地分析得出，切比雪夫滤波器有效地将技术指标的精确度均匀分布，切比雪夫I型滤波器的振幅特性在通带内是等波纹的、在阻带内是单调的；切比雪夫II型滤波器的振幅特性在通带内是单调的、在阻带内是等波纹的。 2.1.3 椭圆低通滤波器设计 切比雪夫I型滤波器在通带内成等波纹振荡，在阻带内却仍是单调下降的，切比雪夫II型在阻带内是等波纹的，在通带内却是单调下降的。因此过渡带的特性有所提高，但是并不理想。它的主要原因在于两者的系统函数在截止频率附近没有有限个零点，其零点在无限远处。1931年，考尔提出了采样有限零点设计的滤波器，因为这种方法在确定零点位置时与椭圆函数的许多特性有关，所以称之为椭圆低通滤波器。它的平方幅度响应函数为： (2-5) 式中 是雅可比椭圆函数， 是与通带衰减有关的函数。滤波器阶次N等于通带和阻带内最大点和最小点的和。MATLAB信号处理工具箱为低通模拟椭圆滤波器的产生提供了函数ellipap,其调用的格式为： [z,p,k]= ellipap (N,Rp,Rs)， 其中，z表示零点,p表示极点,k表示增益，N表示阶次,Rp为通带波纹(dB), Rs为阻带波纹(dB)。下面举例设计椭圆低通模拟原型滤波器： %椭圆低通滤波器 clear all; n=0:0.01:2; for i=1:4 switch i case 1 N=2; case 2 N=3; case 3 N=4; case 4 N=5; end Rp=1; %Rp--通带纹波系数为1 dB Rs=10; %Rs--阻带纹波系数为10dB [z,p,k]=ellipap(N,Rp,Rs); %函数ellipap--设计椭圆低通滤波器 [b,a]=zp2tf(z,p,k); [H,w]=freqs(b,a,n); magH2=(abs(H)).^2; %输出图形 posplot=['22'num2str(i)]; subplot(posplot) plot(w,magH2) axis([0 2 0 1.1]); xlabel('w/wc'); ylabel('|H(jw)|^2'); title(['N='num2str(N)]); grid on end 图2-6 椭圆滤波器的幅频特性(基于MATLAB实现) 由图2-6可得出，椭圆滤波器在通带及阻带内均成等波纹振荡，并具有最窄的过渡带，其性质比较优越，但椭圆滤波器的转移函数不但有极点，而且在 轴上还有零点，所以椭圆滤波器的设计比较复杂。 2.2 模拟-数字滤波器变换及其MATLAB实现 在设计了模拟低通滤波器后，就可以把它们变成数字滤波器了。这些变换均是复值映射，许多文献对此都有研究，根据数字滤波器所保持的模拟滤波器的不同特性，研究出不同的变换技术。其中，最重要的有两种：脉冲响应不变法（保持脉冲响应不变，又叫冲激响应不变法）和双线性Z变换法（保持系统函数不变）。 2.2.1 脉冲响应不变法 脉冲响应不变法的设计原理是使数字滤波器的单位抽样响应序列h(n)，模仿模拟滤波器的脉冲响应g(t)。 设系统传递函数为G(s)的模拟滤波器的单位脉冲响应为g(t),并将脉冲响应g(t)进行等间隔采样，使得数字滤波器的单位抽样响应h(n)刚好等于g(t)的采样值，即： (2-6) 其中的 为采样周期。 G(s)是模拟滤波器的系统传递函数，又令H(z)是数字滤波器的系统传递函数。采样信号的拉式变换与相应的采样序列Z变换的映射关系为： (2-7) 所以系统函数G(s)和H(z)的关系为： (2-8) 式(2-8)的物理意义为首先将模拟滤波器的系统函数G(s)作周期的延拓，在经过式(2-7)的映射变换，映射到Z平面上，从而得到数字滤波器的系统函数H(z)。且模拟和数字频率满足下列关系：ω=ΩT。经过式(2-7)的映射，s平面的左半平面映射为Z平面的单位圆内，因此，一个因果的和稳定的模拟滤波器映射成因果的和稳定的数字滤波器。 经过以上分析，按照脉冲响应不变法，通过模拟滤波器的系统传递函数G(s)，可直接求得数字滤波器系统函数H(Z),其设计具体步骤归纳如下： (1)利用ω=ΩT（可由关系式 推出），将数字滤波器指标 ， 转换为模拟滤波器指标 ， (2)根据指标 ， 来设计模拟滤波器G(s) (3)利用部分分式展开法，把G(s)展成 (2-9) (4)最后把模拟极点 转换为数字极点 ，得到数字滤波器： (2-10) 根据上述理论，将举例在MATLAB环境下用函数实现脉冲响应不变法设计一数字低通滤波器。其函数为[b,a]=impinvar(c,d,T)，其中，b表示数字滤波器自变量为 的分子多项式，a表示数字滤波器自变量为 的分母多项式，c表示模拟滤波器自变量为s的分子多项式,d表示模拟滤波器自变量为s的分母多项式，T表示采样变换参数。 %设模拟低通巴特沃斯滤波器，通带纹波为Rp=1dB,通带上限角频率ωp=0.2п,阻带下限角频率ωs=0.3п，带阻最小衰减αs=15dB，根据该低通模拟滤波器，利用冲激响应不变法设计相应的数字低通滤波器 wp=0.2*pi; ws=0.3*pi; Rp=1; As=15; T=1; %采样周期 %性能指标 Rip=10^(-Rp/20); Atn=10^(-As/20); OmgP=wp*T; OmgS=ws*T; [N,OmgC]=buttord(OmgP,OmgS,Rp,As,'s'); %选取模拟滤波器的阶数 [cs,ds]=butter(N,OmgC,'s'); %设计出所需的模拟低通滤波器[b,a]=impinvar(cs,ds,T); %应用脉冲响应不变法进行转换 [db,mag,pha,grd,w]=freqz_m(b,a); %求得相对，绝对频响及相位、群延迟响应 %绘制各条曲线 subplot(2,2,1); plot(w/pi,mag); title('幅频特性'); xlabel('w(/pi)'); ylabel('|H(jw)|'); axis([0,1,0,1.1]); set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); set(gca,'YTickMode','manual','YTick',[0 Atn Rip 1]); grid subplot(2,2,2); plot(w/pi,db); title('幅频特性(dB)'); xlabel('w(/pi)'); ylabel('dB'); axis([0,1,-40,5]); set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); set(gca,'YTickMode','manual','YTick',[-40 -As -Rp 0]); grid subplot(2,2,3); plot(w/pi,pha/pi); title('相频特性'); xlabel('w(/pi)'); ylabel('pha(/pi)'); axis([0,1,-1,1]); set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); grid subplot(2,2,4); plot(w/pi,grd); title('群延迟'); xlabel('w(/pi)'); ylabel('Sample'); axis([0,1,0,12]); set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); grid %本设计用到的M文件函数 function[db,mag,pha,grd,w]=freqz_m(b,a) %滤波器幅值响应（绝对、相对）、相位响应及群延迟 %Usage:[db,mag,pha,grd,w]=freqz_m(b,a) %500点对应[0,pi] %db 相对幅值响应； mag 绝对幅值响应； pha 相位响应； grd 群延迟响应； %w 采样频率； b 系统函数H(z)的分子项 %a 系统函数H(z)的分母项 [H,w]=freqz(b,a,500); %500点的幅频响应 mag=abs(H); db=20*log10(mag/max(mag)); pha=angle(H); grd=grpdelay(b,a,w); 图2-7 脉冲响应不变法设计数字低通滤波器的特性(基于MATLAB实现) 总结以上，脉冲响应不变法的优点是频率坐标变换是线性的，即ω=ΩT，如不考虑频率混叠现象，用这种方法设计数字滤波器会很好的重现原模拟滤波器的频率响应。另外一个优点是数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域逼近好。但其也具有很大的缺点，若抽样频率不高或其它原因将产生混叠失真，不能重现原模拟滤波器频率响应。所以，脉冲响应不变法适合低通、带通滤波器设计，不适合高通、带阻滤波器的设计。 2.2.2 双线性Z变换法 利用脉冲响应不变法设计数字滤波器时，由于ω=ΩT的频率关系是根据 推导的，所以是 轴每隔2π/T便映射到单位圆上一周，引起了频域混叠。为克服这一现象，人们找到了另一种映射关系： (2-11) 此关系称为双线性Z变换法。 双线性Z变换法的基本思路是：首先将整个s平面压缩到s1平面的一条带宽为2π/T（丛-π/T到π/T）的横带里，然后通过标准的变化关系 将横带变换成整个Z平面上去，这样就得到s平面与Z平面间的一一对应的单值关系，整个过程如图2-8所示： jΩ jΩ1 jIm(Z) π/T 0 б 0 б 0 1 б -π/T s平面 s1平面 Z平面 图2-8 双线性Z变换法的映射关系 由式(2-11)得 (2-12) 及 (2-13) (2-14) 式(2-11)及式(2-12)给出了s和z之间的映射关系，而式(2-13) 和式(2-14)给出了Ω和ω之间的映射关系，但这是一种非映射关系，双线性Z变换法正是利用了正切函数的非线性特点，把整个jΩ轴压缩到了单位圆的一周上。 在MATLAB中，双线性Z变换可通过bilinear函数实现,其调用格式为：[Bz,Az]=bilinear(B,A,Fs)，其中B,A为模拟滤波器的传递函数G(s)的分子分母多项式系数分量，而Bz,Az为数字滤波器的传递函数H(Z)的分子分母多项式的系数分量。 %使用双线性Z变换设计一低通数字滤波器，fp=100HZ,fs=300HZ，αp=3dB,αs=20dB,抽样频率为Fs=1000HZ fp=100;fs=300;Fs=1000; rp=3;rs=20; wp=2*pi*fp/Fs; ws=2*pi*fs/Fs; Fs=Fs/Fs %使Fs为一 wap=tan(wp/2); was=tan(ws/2); [n,wn]=buttord(wap,was,rp,rs,'s') [z,p,k]=buttap(n); [bp,ap]=zp2tf(z,p,k) [bs,as]=lp2lp(bp,ap,wap) [bz,az]=bilinear(bs,as,Fs/2) [H,w]=freqz(bz,az,256,Fs*1000); 图2-9 双线性Z变换设计数字低通滤波器的特性(基于MATLAB实现) 2.3 频带变换及其MATLAB实现 除了低通数字滤波器之外，实际中还常常需要高通、带通、及带阻数字滤波器。设计这三种滤波器的方法有很多，例如基于模拟滤波器转换法的IIR数字滤波器设计，基于直接数字域法的IIR数字滤波器设计等等。下面将具体介绍基于滤波器转换法的IIR数字滤波器设计—直接由模拟低通滤波器转换成数字高通、带通或带阻滤波器。 2.3.1 模拟低通滤波器转换成数字高通滤波器 若已知模拟低通滤波器的系统传递函数为G(s)，则模拟低通滤波器s平面到数字高通滤波器的Z平面的变换公式为： (2-15) 频率变换公式为： (2-16) 从而得到数字高通滤波器的传递函数表达式为： (2-17) 接下来，基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字高通滤波器，要求通带截止频率为通带截至频率为0.6π,通带内衰减不大于1dB,阻带起始频率为0.4π，阻带内衰减不小于15dB,采样周期为Ts=1。 clear all; Wp=0.6*pi; Ws=0.4*pi; Ap=1; As=15; [N,wn]=buttord(Wp/pi,Ws/pi,Ap,As) %计算巴特沃斯滤波器阶次和截至频率 [b,a]=butter(N,wn,'high'); %频率变换法设计巴特沃斯高通滤波器 [b0,B,A]=dir2cas(b,a) %数字高通滤波器级联型 [db,mag,pha,grd,w]=freqz_m(b,a) %数字滤波器响应 subplot(211); plot(w/pi,mag); title('高通数字滤波器幅频响|H(ejOmega)|') subplot(212); plot(w/pi,db); title('高通数字滤波器幅频响(db)') 其MATLAB运行结果为： b0 =0.0751 B = 1.0000 -2.0000 1.0000 1.0000 -2.0000 1.0000 A = 1.0000 0.1562 0.4488 1.0000 0.1124 0.0425 所以由巴特沃斯模拟滤波器设计的数字高通滤波器的传递函数表达式为： 图2-10 基于巴特沃斯模拟滤波器设计的数字高通滤波器的幅频响应曲线 2.3.2 模拟低通滤波器转换成数字带通滤波器 若已知模拟低通滤波器的系统传递函数为G(s)，则模拟低通滤波器s平面到数字带通滤波器的Z平面的变换公式为： (2-18) 频率变换公式为： (2-19) 其中： (2-20) 式中 为下截止频率， 为上截止频率。 从而可得数字带通滤波器的传递函数表达式为： (2-21) 接下来，基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字带通滤波器，要求带通上下截止频率为0.4π,0.3π,带通内衰减不大于3dB，阻带上下起始频率为0.5π，0.2π，阻带内衰减不小于18dB。 clear all; Wp=[0.3*pi,0.4*pi]; Ws=[0.2*pi,0.5*pi]; Ap=3; As=18; [N,wn]=buttord(Wp/pi,Ws/pi,Ap,As); %计算巴特沃斯滤波器阶次和截至频率 [b,a]=butter(N,wn,'bandpass'); %频率变换法设计巴特沃斯带通滤波器 [b0,B,A]=dir2cas(b,a) %数字带通滤波器级联型 [db,mag,pha,grd,w]=freqz_m(b,a); %数字滤波器响应 subplot(211); plot(w/pi,mag); title('带通数字滤波器幅频响|H(ejOmega)|') subplot(212); plot(w/pi,db); title('带通数字滤波器幅频响(db)') 其MATLAB运行结果为： b0 =0.0213 B = 1.0000 -2.0000 1.0000 1.0000 2.0000 1.0000 A = 1.0000 -0.6249 0.7840 1.0000 -1.0053 0.8061 所以由巴特沃斯模拟滤波器设计的数字带通滤波器的传递函数表达式为： 图2-11 基于巴特沃斯模拟滤波器设计的数字带通滤波器的幅频响应曲线 2.3.3 模拟低通滤波器转换成数字带阻滤波器 若已知模拟低通滤波器的系统传递函数为G(s)，则模拟低通滤波器s平面到数字带阻滤波器的Z平面的变换公式为： (2-22) 频率变换公式为： (2-23) 其中： (2-24) 式中 为下截止频率， 为上截止频率。 从而可得数字带阻滤波器的传递函数表达式为： (2-25) 接下来，基于MATLAB利用巴特沃斯模拟滤波器举例，设计一数字带通滤波器，要求通带上下截止频率为0.8π,0.2π,通带内衰减不大于1dB，阻带上下起始频率为0.7π，0.4π，阻带内衰减不小于30dB。 clear all; Wp=[0.4*pi,0.7*pi]; Ws=[0.2*pi,0.8*pi]; Ap=1; As=30; [N,wn]=buttord(Wp/pi,Ws/pi,Ap,As); %计算巴特沃斯滤波器阶次和截至频率 [b,a]=butter(N,wn,'stop'); %频率变换法设计巴特沃斯带阻滤波器 [b0,B,A]=dir2cas(b,a) %数字带阻滤波器级联型 [db,mag,pha,grd,w]=freqz_m(b,a); %数字滤波器响应 subplot(211); plot(w/pi,mag); title('数字带阻滤波器幅频响|H(ejOmega)|') subplot(212); plot(w/pi,db); title('数字带阻滤波器幅频响(db)') 其MATLAB运行结果为： b0 =0.1054 B =1.0000 0.3563 1.0021 1.0000 0.3548 0.9968 1.0000 0.3527 1.0053 1.0000 0.3496 0.9948 1.0000 0.3475 1.0032 1.0000 0.3460 0.9979 A =1.0000 1.1450 0.8087 1.0000 0.8469 0.5087 1.0000 0.4608 0.2943 1.0000 -0.0222 0.2631 1.0000 -0.3959 0.4508 1.0000 -0.6617 0.7750 图2-12 基于巴特沃斯模拟滤波器设计的数字带阻滤波器的幅频响应曲线 2.4 小结 这一章主要是用MATLAB语言进行IIR滤波器的设计和实现。IIR滤波器的设计步骤分为三步，即模拟低通滤波器设计，模拟-数字滤波器变换，滤波器的频带变换。 在模拟低通滤波器的设计中，主要讨论了三种设计方法；在模拟-数字滤波器变换中，讨论了两种变换方法，即脉冲响应不变法和双线性Z变换法；在频带变换的实现中，主要以巴特沃斯滤波器为例并结合MATLAB信号处理工具箱中提供的几个相关函数来进行分析设计。 整个设计过程都是在理论分析的基础上，用MATLAB语言来进行编程设计，并最终通过具体滤波器指标来加以实现的。 第3章 FIR滤波器设计及其MATLAB实现 IIR数字滤波器的设计方法是利用模拟滤波器成熟的理论及设计图表进行的，因而保留了一些典型模拟滤波器优良的幅度特性，但设计中只考虑到了幅度特性，没考虑到相位特性，所设计的滤波器相位特性一般是非线性的。为得到线性相位特性，必须增加相位校正网络，使滤波器设计变得复杂。而FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数H(z)为 H(z)是 的(N-1)次多项式，它在z平面上有(N-1)个零点，原点z=0是(N-1)阶重极点。因此，H(z)永远稳定。稳定和线性相位特性是FIR滤波器突出的特点。 FIR滤波器的设计方法与IIR数字滤波器的设计方法有很大的不同。FIR滤波器的设计任务是选择有限长度的h(n)，使传递函数 满足技术要求。 3.1 线性相位FIR数字滤波器的条件和特点 3.1.1 线性相位条件 对于长度为N的h(n),传输函数为 (3-1) (3-2) 式中， 称为幅度特性， 称为相位特性。 的线性相位是指 是 的线性函数，即 ， 为常数 (3-3) 如果 满足下式： ， 是起始相位 (3-4) 以上两种情况都满足群延迟是一个常数，即 一般称满足(3-3)为第一类线性相位；满足(3-4)为第二类线性相位。第一类线性相位特性是h(n)是实序列且对(N-1)/2偶对称，即h(n)=h(N-n-1)。第二类线性相位特性是h(n)是实序列且对(N-1)/2奇对称，即h(n)=-h(N-n-1)。 3.1.2 线性相位FIR滤波器幅度特性 的特点 由于h(n)的长度N取奇数还是偶数，对 的特性有影响。因此，对于两类线性相位，下面将分四种情况讨论其幅度特性特点。 (1) h(n)=h(N-n-1)，N=奇数，其幅度特性的特点是对ω=0，π，2π是偶对称的。 (2) h(n)=h(N-n-1)，N=偶数，其幅度特性的特点是对ω=π奇对称的，且在ω=π处有一零点，使 ，这样，对于高通和带阻不适合采用这种情况。 (3) h(n)=-h(N-n-1)，N=奇数，其幅度特性的特点在ω=0，π，2π处为零，即在z= 1处是零点，且 对ω=0，π，2π呈奇对称形式。 (4) h(n)=-h(N-n-1)，N=偶数，其幅度特性 在ω=0，2π处为零，即在z=1处有一个零点，且对ω=0，2π奇对称，对ω=π呈偶对称。 3.1.3 线性相位FIR滤波器零点分布特点 第一类和第二类线性相位的系统函数综合起来表示为： (3-5) 上式表明，如 是H(z)的零点，其倒数 也必然是其零点；又因为h(n)是实序列，H(z)的零点必定共轭成对，因此 和 也是其零点。这样，线性相位FIR滤波器零点分布特点是零点必须是互为倒数的共轭对，确定其中一个，另外三个零点也就确定了。 3.2 常用窗函数及MATLAB实现 窗函数在设计FIR数字滤波器中有很重要的作用，正确的选择窗函数可以提高所设计的数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶数。因此必须对各种窗函数有相应的了解。 3.2.1 常用窗函数介绍 1、矩形窗(Rectangular window) 这是一种最简单的窗函数，但从阻带衰减的角度来看，其性能最差。它的定义为: (3-6) 其频率响应为： (3-7) 主瓣宽度为4π/N，第一副瓣比主瓣低13dB。 在MATLAB中，实现矩形窗的函数为rectwin，其调用格式为w=rectwin (N)，其中N是窗函数的长度，返回值w是一个N阶的向量，它由窗函数的值组成。 2、三角形窗(Triangular window) (3-8) 其频率响应为： (3-9) 其主瓣宽度为8π/N，第一副瓣比主瓣低26dB。 在MATLAB中，实现矩形窗的函数为triang，其调用格式为w=triang(N)。 3、汉宁窗(Hanning window) 汉宁窗函数又称升余弦函数，其时域表达式为： (3-10) 其频率响应为： (3-11) 其主瓣宽度为8π/N，能量更集中在主瓣之中。 在MATLAB中，实现矩形窗的函数为hann，其调用格式为w=hann(N)。 4、海明窗(Hamming window) 海明窗函数是一种改进的升余弦函数，其时域表达式为： (3-12) 其幅度频率特性为： (3-13) 其主瓣宽度为8π/N，能量更集中在主瓣之中，主瓣的能量约占99.96%，第一旁瓣的峰值比主瓣小40dB。 在MATLAB中，实现矩形窗的函数为hamming，其调用格式为w=hamming(N)。 5、布拉克曼(Blackman window) (3-14) 布拉克曼的主瓣宽度是矩形窗主瓣宽度的三倍，为12π/N。它的最大旁瓣比主瓣值低57dB。 在MATLAB中，实现矩形窗的函数为blackman，其调用格式为w=blackman(N)。 3.2.2 各种窗函数的实现与比较 在FIR数字滤波器的设计中，要寻找适当的窗函数来进行处理，调整窗口长度N可以有效地控制过渡带宽度；调整窗函数形状可使带内波动减少以及加大阻带衰减。如果找到的窗函数形状，使其谱函数的主瓣包含更多的能量，则相应的旁瓣幅度就能减小，而旁瓣的减小可使通带阻带波动减小，加大阻带衰减，但这样总是以加宽过渡带为代价的。所以要根据实际需要和技术要求来选取窗函数进行设计。 在MATLAB中，各个窗函数的调用函数均已给出(见3.2.1节)且用法比较简单，所以为了更好的比较、区别、分析及选择各个窗函数，特对各函数的幅度特性进行比较。 图3-1给出了以上五种窗函数的时域幅度波形，图3-2给出了当N=64时五种窗函数的频域幅度谱。可以看出，随着旁瓣的减小，主瓣宽度相应增加了。 图3-1 常用窗函数时域幅度曲线比较 (a)矩形窗 (b)三角窗 (c)汉宁窗 (d)