基于MATLAB的FIR数字滤波器的设计与优化.pdf

I 摘 要 随着信息时代和数字世界的到来，数字信号处理已成为今一门极其重要的学科和技术领域。数字信号处理在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理应用中，数字滤波器十分重要并已获得广泛应用。数字滤波是数字信号处理的重要内容，数字滤波器可分为 IIR 和 FIR 两大类。对于 FIR 数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。本文根据 FIR 滤波器的特点，在 MATLAB 坏境下用窗函数设计 FIR 数字滤波器，并对信号进行分析，最后给出了 FIR 带通滤波器对信号滤波的效果。关键词：数字滤波器；FIR；MATLAB；窗函数 Abstract With the information age and the advent of the digital world,digital signal processing has become today a very important discipline and technical field.Digital signal processing in communications,voice,image,automatic control,radar,military,aerospace,medical and household appliances and many other fields has been widely used.In digital signal processing applications,the digital filter is very important and has been widely used.Digital filtering is an important part of digital signal processing,digital IIR and FIR filters can be divided into two broad categories.For the FIR digital filter design,can be given according to the frequency characteristics of the direct design,design method used in the text is a window function method.Based on the characteristics of FIR filters in MATLAB under the bad with the window function throughout the design FIR digital filters,and signal analysis,and finally gives a FIR band-pass filter the signal filtering effect.Key words：digital filter;FIR;MATLAB;window function II 目 录 第一章 绪 论 1 1.1 数字滤波器的研究背景与意义 1 1.2 数字滤波器的发展趋势 2 第二章 数字滤波器的设计原理 4 2.1 数字滤波器的基本结构 4 2.2 滤波器的性能指标 6 2.3 FIR 滤波器的基本结构 7 第三章 FIR 数字滤波器设计方法 9 3.1 窗函数法设计 FIR 滤波器 9 3.2 频率抽样法设计 FIR 滤波器 11 3.3 切比雪夫逼近法设计 FIR 滤波器 13 3.4 FIR 数字滤波器的最优化设计 15 3.4.1 等波纹切比雪夫逼近准则 15 3.4.2 Remez 算法 15 3.5 基于等波纹切比雪夫逼近准则的 FIR DF 的最优化设计步骤 16 第四章 FIR 数字滤波器的 MATLAB 仿真 17 4.1 典型窗函数及其调用格式 17 4.2 基于窗函数的 FIR 滤波器的 MATLAB 实现 17 4.2.1 理想低通滤波器实现 17 4.2.2 系统各响应函数 18 4.2.3 滤波器主函数 18 4.3 滤波器主程序解析 21 第五章 仿真实验结果及分析 23 5.1 原信号的时域波形和处理后时域波形 23 5.2 原信号频域波形与处理后信号的频域波形 23 结束语 25 参考文献 26 基于 MATLAB 的 FIR 数字滤波器的设计与优化 绪论 1 第一章 绪 论 模拟滤波器与数字滤波器的设计对工程，应用数学及计算机科学都是非常重要的。对设计人员来说，滤波器是控制，信号处理和通信领域的重要组成部分。近年来，不论是在科学技术研究，还是在产品的开发等方面，数字信号处理应用越来越广泛，并取得了丰硕的成果。数字滤波占有极其重要的地位。像处理、模式识别、谱分析等应用中的一个基本处理算法。数字滤波是语音和图在许多信号处理应用中用数字滤波器替代模拟滤波器具有许多优势。数字滤波器容易实现不同的幅度和相位频率特性指标，克服了与模拟滤波器器件性能相关的电压漂移、温度漂移和噪声问题。在数字信号处理中，数字滤波是其基本处理方法之一，占有极其重要的地位。数字信号发展过程中的另一个重大进展是数字滤波器按单位脉冲响应 h(n)的长度分类可分有限脉冲响应(FIR)滤波器和无限脉冲响应(IIR)滤波器。两者各有优缺点：IIR 滤波器能以较低的阶次获得相同的幅度滤波性能，但一般为非线性相位；FIR滤波器单位脉冲响应是有限长的，系统必定稳定，且可以做成严格的线性相位，故在图像处理、数据传输等需要信道具有线性相位特性的场合应用广泛。FIR 滤波器的设计方法有窗函数法、频率抽样法等，两种方法分别从时域和频域为出发点来进行设计。有限长单位冲激响应(FIR)数字滤波器，与传统的通过硬件电路实现的模拟滤波器相比有以下优点:(1)简化了硬件电路的设计，提高了硬件电路的集成度和可靠性。(2)对干扰信号的抑制能力有了明显提高，这对系统的控制精度和稳定性的提高起到了促进作用。(3)数字滤波器的参数调节比起模拟滤波器来更加方便、灵活。(4)数字滤波器可以实现数据的并行处理，提高了系统运行速度。1.1 数字滤波器的研究背景与意义 当今，数字信号处理(DSP：Digtal Signal Processing)技术正飞速发展，它不但自成一门学科，更是以不同形式影响和渗透到其他学科：它与国民经济息息相关，与国防建设紧密相连；它影响或改变着我们的生产、生活方式，因此受到人们普遍的关注。数字化、智能化和网络化是当代信息技术发展的大趋势，而数字化是智能化和基于 MATLAB 的 FIR 数字滤波器的设计与优化 绪论 2 网络化的基础，实际生活中遇到的信号多种多样，例如广播信号、电视信号、雷达信号、通信信号、导航信号、射电天文信号、生物医学信号、控制信号、气象信号、地震勘探信号、机械振动信号、遥感遥测信号，等等。上述这些信号大部分是模拟信号，也有小部分是数字信号。模拟信号是自变量的连续函数，自变量可以是一维的，也可以是二维或多维的。大多数情况下一维模拟信号的自变量是时间，经过时间上的离散化(采样)和幅度上的离散化(量化)，这类模拟信号便成为一维数字信号。因此，数字信号实际上是用数字序列表示的信号，语音信号经采样和量化后，得到的数字信号是一个一维离散时间序列；而图像信号经采样和量化后，得到的数字信号是一个二维离散空间序列。数字信号处理，就是用数值计算的方法对数字序列进行各种处理，把信号变换成符合需要的某种形式。例如，对数字信号经行滤波以限制他的频带或滤除噪音和干扰，或将他们与其他信号进行分离；对信号进行频谱分析或功率谱分析以了解信号的频谱组成，进而对信号进行识别；对信号进行某种变换，使之更适合于传输，存储和应用；对信号进行编码以达到数据压缩的目的。数字滤波技术是数字信号分析、处理技术的重要分支。无论是信号的获取、传输，还是信号的处理和交换都离不开滤波技术，它对信号安全可靠和有效灵活地传输是至关重要的。在所有的电子系统中，使用最多技术最复杂的是数字滤波器。数字滤波器的优劣直接决定产品的优劣。1.2 数字滤波器的发展趋势 在信号处理过程中，所处理的信号往往混有噪音，从接收到的信号中消除或减弱噪音是信号传输和处理中十分重要的问题。根据有用信号和噪音的不同特性，提取有用信号的过程称为滤波，实现滤波功能的系统称为滤波器。在近代电信设备和各类控制系统中，数字滤波器应用极为广泛，这里只列举部分应用最成功的领域。1.语音处理 语音处理是最早应用数字滤波器的领域之一，也是最早推动数字信号处理理论发展的领域之一。该领域主要包括 5 个方面的内容：第一，语音信号分析。即对语音信号的波形特征、统计特性、模型参数等进行分析计算；第二，语音合成。即利用专用数字硬件或在通用计算机上运行软件来产生语音；第三，语音识别。即用专用硬件或计算机识别人讲的话，或者识别说话的人；第四，语音增强。即从噪音或干扰中提取被掩盖的语音信号。第五，语音编码。主要用于语音数据压缩，目前已经建立了一系列语音编码的国际标准，大量用于通信和音频处理。近年来，这 5 个方面都取得了不少研究成果，并且，在市场上已出现了一些相关的软件和硬件产品。2.图像处理 基于 MATLAB 的 FIR 数字滤波器的设计与优化 绪论 3 数字滤波技术以成功地应用于静止图像和活动图像的恢复和增强、数据压缩、去噪音和干扰、图像识别以及层析 X 射线摄影，还成功地应用于雷达、声纳、超声波和红外信号的可见图像成像。3.生物医学信号处理 数字滤波器在医学中的应用日益广泛，如对脑电图和心电图的分析、层析 X 射线摄影的计算机辅助分析、胎儿心音的自适应检测等。4.音乐制作 数字滤波器为音乐领域开辟了一个新局面，在对音乐信号进行编辑、合成、以及在音乐中加入交混回响、合声等特殊效果特殊方面，数字滤波技术都显示出了强大的威力。数字滤波器还可用于作曲、录音和播放，或对旧录音带的音质进行恢复等。基于 MATLAB 的 FIR 数字滤波器的设计与优化 数字滤波器的设计原理 4 第二章 数字滤波器的设计原理 2.1 数字滤波器的基本结构 作为线性时不变系统的数字滤波器可以用系统函数来表示，而实现一个系统函数表达式所表示的系统可以用两种方法：一种方法是采用计算机软件实现；另一种方法是用加法器、乘法器、和延迟器等元件设计出专用的数字硬件系统，即硬件实现。不论软件实现还是硬件实现，在滤波器设计过程中，由同一系统函数可以构成很多不同的运算结构。对于无限精度的系数和变量，不同结构可能是等效的，与其输入和输出特性无关；但是在系数和变量精度是有限的情况下，不同运算结构的性能就有很大的差异。因此，有必要对离散时间系统的结构有一基本认识。数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。IIR 滤波器的特征是，具有无限持续时间冲激响应。这种滤波器一般需要用递归模型来实现，因而有时也称之为递归滤波器。FIR 滤波器的冲激响应只能延续一定时间，在工程实际中可以采用递归的方式实现，也可以采用非递归的方式实现。数字滤波器的设计方法有多种，如双线性变换法、窗函数设计法、插值逼近法和 Chebyshev 逼近法等等。随着 MATLAB 软件尤其是 MATLAB 的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。一个数字滤波器可以用系统函数表示为：(2.1)由这样的系统函数可以得到表示系统输入与输出关系的常系数线形差分程为：00()()()NMkkkky na y nkb x nk(2.2)可见数字滤波器的功能就是把输入序列x(n)通过一定的运算变换成输出序列y(n)。不同的运算处理方法决定了滤波器实现结构的不同。无限冲激响应滤波器的单位抽样响应h(n)是无限长的，其差分方程如(2.2)式所示，是递归式的，即结构上存在着输出信号到输入信号的反馈，其系统函数具有(2.1)式的形式，因此在z平面的有限区间(0)z 有极点存在。前面已经说明，对于一个给定的线形时不变系统的系统函数，有着各种不同的01()()()1MkkkNkkkb zY zH zX za z基于 MATLAB 的 FIR 数字滤波器的设计与优化 数字滤波器的设计原理 5 等效差分方程或网络结构。由于乘法是一种耗时运算，而每个延迟单元都要有一个存储寄存器，因此采用最少常熟乘法器和最少延迟支路的网络结构是通常的选择，以便提高运算速度和减少存储器。然而，当需要考虑有限寄存器长度的影响时，往往也采用并非最少乘法器和延迟单元的结构。IIR 滤波器实现的基本结构有：(1)IIR 滤波器的直接型结构；优点：延迟线减少一半，变为N 个，可节省寄存器或存储单元；缺点：其它缺点同直接I型。通常在实际中很少采用上述两种结构实现高阶系统，而是把高阶变成一系列不同组合的低阶系统(一、二阶)来实现。x(n)y(n)x(n-M)y(n-M)x(n-1)y(n-1)Z-1 Z-1Z-1Z-1Z-1Z-1Z-1Z-1b1b2bM-a1 -a2-a3b0 图 2.1 直接型 Z-1Z-1Z-1Z-1A01A11b11b21b1Lb2LA2LA1LA1 图 2.2 并联型(2)IIR 滤波器的级联型结构；基于 MATLAB 的 FIR 数字滤波器的设计与优化 数字滤波器的设计原理 6 优点：系统实现简单，只需一个二阶节系统通过改变输入系数即可完成；极点位置可单独调整；运算速度快(可并行进行)；各二阶网络的误差互不影响，总的误差小，对字长要求低。缺点：不能直接调整零点，因多个二阶节的零点并不是整个系统函数的零点，当需要准确的传输零点时，级联型最合适。(3)IIR 滤波器的并联型结构。优点：简化实现，用一个二阶节，通过变换系数就可实现整个系统；极、零点可单独控制、调整，调整1i、2i只单独调整了第i对零点，调整1i、2i则单独调整了第i对极点；各二阶节零、极点的搭配可互换位置，优化组合以减小运算误差；可流水线操作。缺点：二阶阶电平难控制，电平大易导致溢出，电平小则使信噪比减小。b1b2 -b1L -b2Lb0Z-1Z-1Z-1Z-1a1La2La21a11-b11-b21 图 2.3 级联型 2.2 滤波器的性能指标 我们在进行滤波器设计时，需要确定其性能指标。一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波器特性为例，频率响应有通带、过渡带及阻带三个范围。在通带内：1)(1jpeHA在阻带中：AeHj)(t 其中c为通带截止频率,st为阻带截止频率，AP为通带误差,AST为阻带误差。与模拟滤波器类似，数字滤波器按频率特性划分为低通、高通、带通、带阻、全通等类型，由于数字滤波器的频率响应是周期性的，周期为2。各种理想数字滤波器的幅度频率响应如图 2.4 所示：基于 MATLAB 的 FIR 数字滤波器的设计与优化 数字滤波器的设计原理 7 图 2.4 各种理想数字滤波器的幅度频率响应 2.3 FIR 滤波器的基本结构 FIR 滤波器的单位抽样响应为有限长度，一般采用非递归形式实现。通常的 FIR数字滤波器有横截性和级联型两种。FIR 滤波器实现的基本结构有：1.FIR 滤波器的横截型结构 表示系统输入输出关系的差分方程可写作：)()()(10mnxmhnyNm (2.3)直接由差分方程得出的实现结构如图 2.5 所示：x(n)Z-1h(1)h(N-2)h(N-1)h(0)Z-1Z-1Z-1y(n)图 2.5 横截型（直接型、卷积型）若 h(n)呈现对称特性，即此 FIR 滤波器具有线性相位，则可以简化加横截型结构，下面分情况讨论：基于 MATLAB 的 FIR 数字滤波器的设计与优化 数字滤波器的设计原理 8 x(n)y(n)h(N-1)/2Z-1Z-1Z-1Z-1Z-1Z-1Z-1h(0)h(1)h(2)h(N-3)/2 图 2.6 N 为奇数时 FIR 滤波器实现结构 Z-1Z-1Z-1Z-1Z-1Z-1Z-1h(0)h(1)h(2)h(N/2-1)图 2.7 N 为偶数时 FIR 滤波器实现结构 FIR 滤波器的级联型结构将 H(z)分解成实系数二阶因子的乘积形式：221121010)()(zbzbbznhzHkknkknnn (2.4)这时FIR滤波器可用二阶节的级联结构来实现，每个二阶节用横截型结构实现。如图所示，这种结构的每一节控制一对零点，因而在需要控制传输零点时可以采用这种结构。b01h(1)Z-1Z-1Z-1Z-1Z-1Z-1b02b0(N-1)b11b21b12b22b1(N-1)b2(N-1)图 2.8 FIR 滤波器的级联结构 基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 9 第三章 FIR 数字滤波器设计方法 设计 FIR 数字数字滤波器的方法通常有三种:窗函数法、频率抽样法、等波纹逼近法。下面我们分别讨论这三种设计方法。3.1 窗函数法设计 FIR 滤波器 工程中比较常用的窗函数11有：矩形窗函数、三角形(Bartlett)窗函数、汉宁(Harming)窗函数、海明(Hamming)窗函数、布莱克曼(Blackman)窗函数和凯塞(Kaiser)窗函数。这里我们主要介绍常用的矩形窗函数、汉宁(Harming)窗函数、。1、矩形窗（Rectangle Window）矩形窗：)()(nRnwN (3.1)窗 谱：121)2sin()2sin()(WjjwReNeW (3.2)幅度函数：)2sin()2sin()(NWR (3.3)用矩形窗设计低通数字滤波器的程序及运行结果如下：omegac=0.37;N=81;m=(N-1)/2;n=0:2*m+10;h=omegac/pi*sinc(omegac*(n-m)/pi);w=ones(1,N)zeros(1,length(n)-N);hd=h.*w;omega=-pi:2*pi/300:pi;Hd=freqz(hd,1,omega);plot(omega,abs(Hd);基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 10 图 3.1 矩形窗设计低通数字滤波器 2、汉宁窗（升余弦窗）其中：)()(25.0)(5.0)()12cos(1 5.0)(1212nRggnRnRNnnWNNnjNnjNN (3.4)利用傅氏变换的移位特性，汉宁窗频谱的幅度函数 W（）可用矩形窗的幅度函数表示为：三部分矩形窗频谱相加，使旁瓣互相抵消，能量集中在主瓣，旁瓣大大减小，主瓣宽度增加 1 倍。)12()12(25.0)(5.0)(NWNWWWRRR 图 3.2 滤波器频率响应 图 3.3 汉宁窗函数的脉冲响应 窗函数：(3.5)频率响应：2)1()(NjhamjhameWeW (3.6)当 N1 时，频率响应的幅度函数)(jhameW的主瓣宽度为8N,第一旁瓣比主瓣低 31dB。表 4.1 几种常用窗函数对比，窗函数的选择原则是：(1)具有较低的旁瓣幅度，尤其是第一旁瓣幅度；(2)旁瓣幅度下降速度要大，以利增加阻带衰减；(3)主瓣的宽度要窄，以获得较陡的过渡带。通常上述三点很难同时满足。当选用主瓣宽度较窄时，虽然得到较陡的过渡带，2()0.51cos()()1HnNnnRnN基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 11 但通带和阻带的波动明显增加：当选用最小的旁瓣幅度时，虽能得到匀滑的幅度响应和较小的阻带波动，但过渡带加宽。因此，实际选用的窗函数往往是它们的折衷。在保证主瓣宽度达到一定要求的条件下，适当牺牲主瓣宽度来换取旁瓣波动的减少。表 3.1 几种窗函数的比较 窗函数 旁瓣峰值幅度/dB 过渡带宽 阻带最小衰减/dB 矩形窗-13 4/N-12 三角形窗-25 8/N-25 汉宁窗-31 8/N-44 哈明窗-41 8/N-53 布莱克曼窗-57 12/N-74 凯塞窗-57 10/N-80 3.2 频率抽样法设计 FIR 滤波器 FIR 滤波器具有线性相位的条件是h(n)是实序列，且满足h(n)=h(N-n-1)，从而推出传输函数所满足的条件13：)()()(jgjdeHeH （3.7)2)1()(N （3.8)2()(ggHH，N=奇数 （3.9）1,.,2,1,0Nk，N=偶数 （3.10)在20之间等间隔采样 N 点 Nkk2 1,.,2,1,0Nk 如果待设计的滤波器为)(jdeH，对应的单位取样响应为)(hnd deeHnhnjjdxxd)(21)(（3.11）则由频率采样定理知道，在频域20之间等间隔采样N点，利用 IDFT 得到()dHn基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 12 应是()dHn以 N 为周期，周期性延拓乘以()NRn，即 deeHnhnjjdxxd)(21)(（3.12）由于时域混叠，因其所设计的h()n和h()dn有偏差。为此，如果所需的频率采样点数 N 加大。N 愈大，所设计的滤波器愈逼近待设计的滤波器jdH(e)若从频域分析，由采样定理表明，频率域等间隔采样H(k)，经过 IDFT 得到h(n)，其 Z 变换H(z)和H(k)的关系为：10121)()1()(NKNkjNZekHNzzH （3.13）频率采样法的 MATLAB 程序如下：M=40;%取滤波器的阶数为 40 al=(M-1)/2;%群时延 n=0:M-1;T2=0.59417456;T1=0.109021;Hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4);%采样值的幅值 k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);%采样值的相位 H=Hrs.*exp(j*angH);h=real(ifft(H,M);%长度为 M 的单位脉冲响应 stem(n,h);Title(滤波器的实际单位脉冲响应);freqz(h,1,512);Title(幅度响应和相位响应);基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 13 00.10.20.30.40.50.60.70.80.91-2000-1500-1000-5000500Normalized Frequency (rad/sample)Phase(degrees)00.10.20.30.40.50.60.70.80.91-150-100-50050Normalized Frequency (rad/sample)Magnitude(dB)幅 度 响 应 和 相 位 响 应 图 3.4 所设计的带通滤波器的幅度和相位响应 3.3 切比雪夫逼近法设计 FIR 滤波器 除了窗函数设计法（时间窗口法）和频率采样法，还可以用切比雪夫逼近法设计 FIR 滤波器，且切比雪夫逼近法是一种等波纹逼近法，在用切比雪夫逼近法设计 FIR 滤波器时，需要用到雷米兹（Remez）交替算法，并且需要遵从两个准则：均方误差最小准则和最大误差最小化准则。雷米兹（Remez）交替算法14：能很好的解决通带截止频率pw和阻带截止频率sw不能精确控制的问题。下图就是雷米兹（Remez）交替算法的流程图：基于切比雪夫一致逼近法设计 FIR 数字低通滤波器的程序如下：clear all;f=0 0.6 0.7 1;%给定频率轴分点 A=1 1 0 0;%给定在这些频率分点上理想的幅频响应 weigh=1 10;%给定在这些频率分点上的加权 b=remez(32,f,A,weigh);%设计出切比雪夫最佳一致逼近滤波 h,w=freqz(b,1,256,1);h=abs(h);h=20*log10(h);subplot(211)stem(b,.);grid;title (切比雪夫逼近滤波器的抽样值)subplot(212)基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 14 plot(w,h);grid;title (滤波器幅频特性(dB)(r+1)个极值频率的初始猜测值计算极值频率集上的最优q通过 r点内插得到 P(w)计算 E(w)并求出|E(w)|q 处得局部极大值 极值点是否多于（r+1）个极值点位置是否变化舍去=0，中较小的一个局部极值，保留（r+1）个最大极值最优逼近YYNN 图 3.5 雷米兹（Remez）交替算法的流程图 基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 15 05101520253035-0.200.20.40.60.8切 比 雪 夫 逼 近 滤 波 器 的 抽 样 值00.050.10.150.20.250.30.350.40.450.5-100-50050滤 波 器 幅 频 特 性(dB)图 3.6 切比雪夫逼近滤波器的抽样值和滤波器的幅频特性 3.4 FIR 数字滤波器的最优化设计 最优化设计方法15是指采用最优化准则来设计的方法。在 FIR DF 的最优化设计中,最优化准则有均方误差最小化准则和等波纹切比雪夫逼近(也称最大误差最小化)准则两种。实际设计中,只有采用窗函数法中的矩形窗才能满足前一种最优化准则,但由于吉布斯(Gibbs)效应的存在,使其根本不能满足设计的要求。为了满足设计的要求,可以采用其它的窗函数来消除吉布斯效应,但此时的设计已 经不能满足该最优化准则了。因此,要完成 FIR DF 的最优化设计,只能采用后一种优化准则来实现。3.4.1 等波纹切比雪夫逼近准则 在滤波器的设计中,通常情况下通带和阻带的误差要求是不一样的。等波纹切比雪夫逼近准则就是通过对通带和阻带使用不同的加权函数,实现在不同频段(通常指的是通带和阻带)的加权误差最大值相同,从而实现其最大误差在满足性能指标的条件下达到最小值。3.4.2 Remez 算法 Remez 算法是由 Parks 和 McClellan 等人在 1972 年推导出来的。它是将 FIR 基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器设计方法 16 DF 的五个参数(N,1，2,pw,sw,1，2,pw,sw)中的N,pw，sw和12固定,而视1 (或2)为变量的一种迭代方法。尽管 Herrmann 等人在 1971 年也推导出来了将N和1,2固定,而将pw,sw设为变量的另一种迭代方法,但由于前一种算法最灵活且最有效。因而成为最优化设计的主要方法。该方法求解过程为:（1）求解：首先在滤波器的通带和阻带内等间隔地取r 1个频率点(,1.)kkor 作为交错点的初始值,然后利用交错定理中)cos()()(10wnnaPrn的表达式及的解析式求解出满足下式的 值:rkPHWkd,.,1,0)1()()()(3.14)（2）求解()P w:利用求解出来的值和预先假定的r 1个频率点求出()iP(其中0,11ir)的值,然后根据拉格朗日插值公式求出()P的最终表达式。（3）求解()E:将求得的()P代入下式判定其是否满足不等式的要求。)()()()(PHWEd (3.15)若满足要求,则说明已经获得了最优解，若在某些频率点不满足要求,则需要将这些频率点作为新的极值点重新计算,经过反复的迭代直到在所有的频率点上都满足不等式的要求。这时的值就是最终所求的值。这样便获得了最佳逼近。3.5 基于等波纹切比雪夫逼近准则的 FIR DF 的最优化设计步骤（1）给出所需的频率响应)(dH,加权函数)(w以及滤波器的单位抽样响应h(n)的长度N。（2）由中给定的参数来形成所需的,)(dH,)(w,)(P的表达式。（3）根据 Remez 算法,求解逼近问题。（4）根据求得的)(p表达式,利用傅立叶逆变换计算出单位抽样响应)(nh的表达式即可获解。（5）基于 MATLAB 的 FIR DF 最优化设计的仿真实现利用数字信号处理工具箱中的remezord 和 remez 函数可以实现 FIR DF 的最优化设计。在此先介绍这两个函数功能:1.利用remezord 函数可以通过估算得到滤波器的近似阶数n,归一化频率带边界0f,频带内幅值0a及各个频带内的加权系数weights。输入参数f为频带边缘频率,a为各个频带所期望的幅度值,dev是各个频带允许的最大波动。2.利用 remez 函数可以得到最优化设计的 FIR DF 的)(h n系数,输入参数n 是滤波器的阶数,0f,0a,weights参数含义说明同。ftype是所设计的滤波器类型,它除了可以设计普通的滤波器外,它还可以设计数字希尔钞特变换器以及数字微分器。实际设计中,由于 remezord 函数可跑高估或低估滤波器的阶数 n,因此在得到滤波器的系数后,必须检查其阻带最小衰减是否满足设计要求。如果此时的技术指标不能满足设计要求，则必须提高滤波器的阶数到2,1nn等。基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器的 MATLAB 仿真 17 第四章 FIR 数字滤波器的 MATLAB 仿真 4.1 典型窗函数及其调用格式 1.在 MATLAB 中，实现三角窗的函数为 boxcar，调用格式为：w=boxcar(N)%数组 w 中返回 N 点矩形窗函数 2.在 MATLAB 中，实现三角窗的函数为 triang，调用格式为：w=triang(N)%数组 w 中返回 N 点三角窗函数 3.在 MATLAB 中，实现汉宁窗的函数为 hann，调用格式如下：w=hann(N)w=hann(N,sflag)Hann 函数中的参数 sflag 为采样方式，其值可取 symmetric(默认值)或 periodic。当 sflagsymmetric 时，为对称采样；当 sflagperiodic 时，为周期采样，此时 hann函数计算 N+1 个点的窗，但是仅返回前 N 个点。4.在 MATLAB 中，实现海明窗的函数为 hamming，调用格式分别如下：w=hamming(N)w=hamming(N,sflag)其中 sflag 的用法同上。5.在 MATLAB 中，实现布拉克曼窗的函数为 blackman，调用格式如下：w=blackman(N)w=blackman(N,sflag)6.在 MATLAB 中，实现切比雪夫窗的函数为 chebwin，调用格式为：w=chebwin(N,r)其中 r 表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低 rdB(其默认值为100dB)，且旁瓣是等纹波的。7.在 MATLAB 中，实现巴特里特窗的函数为 bartlett，调用格式为：w=bartlett(N)4.2 基于窗函数的 FIR 滤波器的 MATLAB 实现 4.2.1 理想低通滤波器实现 通过先产生一个理想的低通滤波器变换为一理想的带通滤波器,作为实际滤波基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器的 MATLAB 仿真 18 器的模拟对象.该理想低通滤波器由 ideal_lp.m 文件产生.该程序代码如下所示:ideal_lp.m 文件:function hd=ideal_lp(w,N);alpha=(N-1)/2;n=0:(N-1);m=n-alpha+eps;hd=sin(w*m)./(pi*m);4.2.2 系统各响应函数 该函数为求取系统的绝对幅度响应、相对的 db 值幅度响应、相位响应和群延时响应的函数,该函数由 Freqz_m.m 文件实现.该程序代码如下所示:Freqz_m.m 文件 function db,mag,pha,grd,w=freqz_m(b,a)%求取系统的绝对幅度响应、相对的 db 值幅度响应、相位响应和群延时响应的函数%db 为相对振幅(dB)%mag 为绝对振幅%pha 为相位响应%grd 为群延时%w 为频率样本点向量%b 为 Ha(z)分子多项式系数(对 FIR 而言，b=h)%a 为 Hz(z)分母多项式系数(对 FIR 而言，a=1)H,w=freqz(b,a,1000,whole);%freqz 显示数字滤波器频域中的图形 H=(H(1:501);w=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);4.2.3 滤波器主函数 该带通滤波器由 Daitong.m 文件实现.该程序的具体代码如下:Function daitong N1=1024;t=0:1:N1-1;fs=5000;基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器的 MATLAB 仿真 19 s=(sin(2*100*pi*t/fs)+sin(2*pi*1000*t/fs)+2*sin(2*pi*2000*t/fs)+randn(1,length(t/fs);%处理前的信号函数 Y=fft(s);%原信号的傅里叶变换 ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;As=60;tr_width=min(wp1-ws1),(ws2-wp2)M=ceil(11*pi/tr_width)+1 n=0:1:M-1;wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);%理想带通滤波器 w_bla=(blackman(M);%调用窗函数 h=hd.*w_bla;%加窗 db,mag,pha,grd,w=freqz_m(h,1);%求取系统的绝对幅度响应、相对的db 值幅度响应、相位响应和群延时响应的函数 delta_w=2*pi/1000;Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w)As=-round(max(db(ws2/delta_w+1:1:501)Ts=t(2)-t(1);Ws=2*pi/Ts;Wn=Ws/2;y1=filter(h,1,s);%求处理后信号时域波形 Y1=fft(y1);%求处理后信号的频谱 Ya=abs(Y(1:length(t)/2);Ya1=abs(Y1(1:length(t)/2);f=linspace(0,Wn/(2*pi),length(t)/2);figure(1);%打印理想脉冲响应波形和实际脉冲响应的波形 subplot(1,2,1);stem(n,hd);Title (理想脉冲响应)axis(0 M-1-0.4 0.5);基于 MATLAB 的 FIR 数字滤波器的设计与优化 FIR 数字滤波器的 MATLAB 仿真 20 ylabel(hd(n)subplot(1,2,2);stem(n,w_bla);title(Blackman Window)axis(0 M-1 0 1.1);ylabel(w(n)subplot(2,2,3);stem(n,h);Title (实际脉冲响应)axis(0 M-1-0.4 0.5);ylabel(h(n)figure(2);%打印滤波器频率响应和相位响应 freqz(h,1);grid;ylabel(Decibels)axis(0 1-150 10);figure(3);%打印原信号时域波形和处理后的时域波形 subplot(1,2,1)plot(t,s);title (原信号时域图);axis(0,300,-5,5);grid;ss=ifft(Y1(1:length(t);subplot(122)plot(t,ss);title (滤波后的时域图);axis(0,300,-5,5);grid;figure(4);%打印原信号频谱和处理后信号的频谱