FIR带阻滤波器课程设计.doc

摘 要无限长脉冲数字滤波器的设计方法只考虑了幅度特性，没有考虑相位特性，所设的滤波器一般是某种确定的非线性相位特性。有限脉冲响应（fir）滤波器在保证了幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。本课题利用matlab软件实现。matlab是“矩阵实验室”（matrix laboratoy）的缩写，是一种科学计算软件，它使用方便，输入简捷，运算高效，内容丰富，因此利用matlab软件，通过一系列较为系统的函数法，根据已知的技术指标，就可以设计出满足要求的滤波器。关键词：matlab；窗函数；fir带阻数字滤波器；线性相位目 录1课题描述12 matlab简介23 窗函数设计法原理23.1线性相位53.2 基本窗函数53.2.1 矩形窗函数53.2.2 三角窗函数63.2.3汉宁窗函数73.2.4 哈明窗函数73.2.5 布莱克曼窗函数83.2.6 凯塞窗84 方案设计与分析104.1方案设计程序104.2分析125 总结与体会146 参考文献1 课程设计说明书（论文） 1课题描述现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因，使得具有线性相位的fir数字滤波器得到大力发展和广泛应用。在实际进行数字信号处理时，往往需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进行分析。这样，取用有限个数据，即将信号数据截断的过程，就等于将信号进行加窗函数操作。而这样操作以后，常常会发生频谱分量从其正常频谱扩展开来的现象，即所谓的“频谱泄漏”。当进行离散傅立叶变换时，时域中的截断是必需的，因此泄漏效应也是离散傅立叶变换所固有的，必须进行抑制。而要对频谱泄漏进行抑制，可以通过窗函数加权抑制dft的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。而在后面的fir滤波器的设计中，为获得有限长单位取样响应，需要用窗函数截断无限长单位取样响应序列。另外，在功率谱估计中也要遇到窗函数加权问题。由此可见，窗函数加权技术在数字信号处理中的重要地位。第 13 页 共 15页2 matlab简介matlab是“矩阵实验室”（matrix laboratoy）的缩写，是一种科学计算软件，主要适用于矩阵运算及控制和信息处理领域的分析设计，它使用方便，输入简捷，运算高效，内容丰富，因此很多专家在自己擅长的领域用它编写了许多专门的matlab工具包，由于matlab功能的不断扩展，所以是科学研究中最常用必不可少的工具。matlab由一系列工具组成。这些工具方便用户使用matlab的函数和文件，其中许多工具采用的是图形用户界面。包括matlab桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着matlab的商业化以及软件本身的不断升级，matlab的用户界面也越来越精致，更加接近windows的标准界面，人机交互性更强，操作更简单。而且新版本的matlab提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。matlab一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（m文件）后再一起运行。新版本的matlab语言是基于最为流行的c语言基础上的，因此语法特征与c语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是matlab能够深入到科学研究及工程计算各个领域的重要原因。3 窗函数设计法原理数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（iir）滤波器和有限长冲激响应（fir）滤波器。iir数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以iir滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。fir数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。因此设计fir滤波器的方法之一可以从时域出发，截取有限长的一段冲击响应作为h(z)的系数，冲击响应长度n就是系统函数h(z)的阶数。只要n足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的fir数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数： （3-1）其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计fir df用的h(n)，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用h(n)来近似。窗函数的基本思想：先选取一个理想滤波器（它的单位抽样响应是非因果、无限长的），再截取（或加窗）它的单位抽样响应得到线性相位因果fir滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。设x(n)是一个长序列，w(n)是长度为n的窗函数，用w(n)截断x(n)，得到n点序列xn(n)，即 xn(n) = x(n) w(n) （3-2）在频域上则有 （3-3）由此可见，窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形，而且也会影响到频域内的形状。matlab信号工具箱主要提供了以下几种窗函数，如表下所示：表1 matlab窗函数窗窗 函 数矩形窗boxcar三角窗barlett布莱克曼窗blackman哈明窗hamming汉宁窗hanning凯塞窗kaiser切比雪夫窗chebwin加矩形窗后的频谱和理想频谱可得到以下结论：加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。n越大，过渡带越窄、越陡；过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95，与n无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为 （3-4）这样两个肩峰点的相对衰耗分别是0.74db和-21db。其中（-0.0895）对应的点的值定义为阻带最小衰耗。以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波器的关键在于改进窗函数。窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为： 旁瓣峰值衰耗20lg(第一旁瓣峰值主瓣峰值) （3-5）为了改善滤波器的性能，需使窗函数谱满足：主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带；第一副瓣面积相对主瓣面积尽可能小，即能量尽可能集中在主瓣，外泄少，使设计出来的滤波器的肩峰和余振小。但上面两个条件是相互矛盾的，实际应用中，折衷处理，兼顾各项指标。3.1线性相位一个单一频率的正弦信号通过一个系统，假设它通过这个系统的时间需要t，则这个信号的输出相位落后原来信号wt的相位。从这边可以看出，一个正弦信号通过一个系统落后的相位等于它的w*t；反过来说，如果一个频率为w的正弦信号通过系统后，它的相位落后delta，则该信号被延迟了delta/w的时间。在实际系统中，一个输入信号可以分解为多个正弦信号的叠加，为了使得输出信号不会产生相位失真，必须要求它所包含的这些正弦信号通过系统的时间是一样的。因此每一个正弦信号的相位分别落后，w1*t，w2*t，w3*t。因此，落后的相位正比于频率w，如果超前，超前相位的大小也是正比于频率w。从系统的频率响应来看，就是要求它的相频特性是一条直线。在fir滤波器的设计中，为了得到线性相位的性质，通常利用实偶对称序列的相频特性为常数0和实奇对称序列为相频特性为常数90度的特点。因此得到的是对称序列，不是因果序列，是不可实现系统，为了称为物理可实现系统，需要将它向右移动半个周期，这就造成了相移特性随时间的变化，同时也是线性变化。单位脉冲响应h(n)（为实数）具有偶对称或奇对称性，则fir数字滤波器具有严格的线性相位特性。 数字滤波器中，iir数字滤波器方便简单，但它相位的非线性，要求采用全通网络进行相位校正，且稳定性难以保障。fir滤波器具有很好的线性相位特性，使得它越来越受到广泛的重视。3.2 基本窗函数数字信号处理领域中所用到的基本窗函数主要有：矩形窗函数、三角窗函数和汉宁窗函数，哈明窗函数布莱克窗函数，凯塞窗函数等。下面就对这些窗函数展开介绍。3.2.1 矩形窗函数矩形窗(rectangular window)函数的时域形式可以表示为： （3-2-1）它的频域特性为 （3-2-2）boxcar函数：生成矩形窗调用方式w = boxcar (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。从功能上讲，该函数又等价于w = ones(n,1)。3.2.2 三角窗函数三角窗是最简单的频谱函数为非负的一种窗函数。三角窗函数的时域形式可以表示为：当n为奇数时 （3-2-3）当n为偶数时 （3-2-4）bartlett函数：生成巴特利特窗调用方式w = bartlett(n)：(1) 输入参数n是窗函数的长度；(2) 输出参数w是由窗函数的值组成的n阶向量。(3) 三角窗也是两个矩形窗的卷积。三角窗函数的首尾两个数值通常是不为零的。当n是偶数时，三角窗的傅立叶变换总是非负数。3.2.3汉宁窗函数汉宁窗函数的时域形式可以表示为： （3-2-5）它的频域特性为： （3-2-6）其中，为矩形窗函数的幅度频率特性函数。汉宁窗函数的最大旁瓣值比主瓣值低31db，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8/n。hanning函数：生成汉宁窗调用方式(1) w = hanning(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。注意：此函数不返回是零点的窗函数的首尾两个元素。 (2) w = hanning(n,symmetric)：与上面相类似。(3) w = hanning(n,periodic)：此函数返回包括为零点的窗函数的首尾两个元素。3.2.4 哈明窗函数海明窗函数的时域形式可以表示为 （3-2-7）它的频域特性为 （3-2-8）其中，为矩形窗函数的幅度频率特性函数。海明窗函数的最大旁瓣值比主瓣值低41db，但它和汉宁窗函数的主瓣宽度是一样大的。hamming函数：生成海明窗调用方式(1) w = hamming(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。(2) w = hamming(n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。3.2.5 布莱克曼窗函数布莱克曼窗函数的时域形式可以表示为 (3-2-9)它的频域特性为 （3-2-10）其中，为矩形窗函数的幅度频率特性函数。布莱克曼窗函数的最大旁瓣值比主瓣值低57db，但是主瓣宽度是矩形窗函数的主瓣宽度的3倍，为12/n。blackman函数：生成海明窗调用方式(1) w = blackman (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。(2) w = blackman (n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。3.2.6 凯塞窗上面所讨论的几种窗函数，在获得旁瓣抑制的同时却增加了主瓣的宽度。而凯塞窗定义了一组可调的窗函数，它是由零阶贝塞尔函数构成的，其主瓣能量和旁瓣能量的比例是近乎最大的。而且，这种窗函数可以在主瓣宽度和旁瓣高度之间自由选择它们的比重，使用户的设计变得非常灵活。凯塞窗函数的时域形式可表示为 (3-2-11)其中，是第1类变形零阶贝塞尔函数，是窗函数的形状参数，由下式确定： （3-2-12）其中，为凯塞窗函数的主瓣值和旁瓣值之间的差值(db)。改变的取值，可以对主瓣宽度和旁瓣衰减进行自由选择。的值越大，窗函数频谱的旁瓣值就越小，而其主瓣宽度就越宽。kaiser函数：生成凯塞窗调用方式w = kaiser(n,beta)：输入参数n是窗函数的长度；输入参数beta用于控制旁瓣的高度；输出参数w是由窗函数的值组成的n阶向量。n一定时，beta越大，其频谱的旁瓣就越小，但主瓣宽度相应的增加；当beta一定时，n发生变化，其旁瓣高度不会发生变化。4 方案设计与分析用窗函数法设计一个fir带阻滤波器。指示如下：下通带截至频率 ；上通带截止频率 ；阻带下限频率 ；阻带上限频率 通带最大衰减 阻带最小衰减 =4.1方案设计程序6中窗函数的基本参数如下窗函数类型旁瓣峰值/db过渡带宽度阻带最小衰减/db近似值精确值矩形窗-134/n1.8/n-12三角窗-258/n6.1/n-25汉宁窗-318/n6.2/n-44哈明窗-418/n6.6/n-53布莱克曼窗-5712/n11/n-74凯塞窗-5710/n-80因为阻带最小衰减=,所以选择布莱克曼窗或凯塞窗都可以设计，先以布莱克窗进行设计。程序步骤如下：wlp=0.2*pi; wls=0.35*pi; wus=0.65*pi; wup=0.8*pi;wc=(wlp+wls)/2/pi,(wus+wup)/2/pi; b=wls-wlp; n=ceil(12*pi/b)-1; n=0:n-1; window=kaiser(n); h1,w=freqz(window,1) subplot(2,2,1)stem(window,.); xlabel(n); title( kaiser窗函数); subplot(2,2,2) plot(w/pi,20*log(abs(h1)/abs(h1(1); grid; xlabel(w/pi); ylabel(幅度(db); title( kaiser窗函数的频谱); hn = fir1(n-1,wc,stop); h2,w=freqz(hn,1,512); subplot(2,2,3) stem(n,hn,.); xlabel(n); ylabel(h(n); title( kaiser窗函数的单位脉冲响应); subplot(2,2,4) plot(w/pi,20*log(abs(h2)/abs(h2(1); grid; xlabel(w/pi); ylabel(幅度(db); title( kaiser带阻滤波器的幅度特性);运行图形如下：4.2分析由kaiser窗函数程序可知，n-1=80,所以n=81，通过图形“kaiser带阻滤波器的幅度特性”可知通带内最大衰减