窗函数设计FIR数字滤波器论文

1、摘 要在数字信号处理中，数字滤波器占有极其重要的地位。数字滤波器是一个能够完成待定任务的离散时间系统，它可以利用有限精度算法来实现。现代数字滤波器可以用软件或设计专用的数字处理硬件两种方式来实现。用软件实现数字滤波器的优点是，随着滤波器参数的改变，很容易改变滤波取的性能。有限长单位冲激响应（FIR）数字滤波器因其严格的线性相位而在数字通信、语音信号处理及自适应信号处理等领域得到了广泛的应用。本文讲述了有限长单位冲激响应（FIR）数字滤波器的原理，对FIR数字滤波器的窗函数设计方法进行了研究。窗函数法在FIR数字滤波器的设计中有着广泛的应用。传统的数字滤波器设计方法繁琐而且不直观，本文利用MAT

2、LAB具有强大的科学计算和图形显示这一优点，与窗函数法设计理论相结合共同设计FIR数字滤波器，不但使设计结果更加直观，而且提高了滤波器的设计精度，从而更好的达到预期效果。关键字：数字信号处理；FIR滤波器；窗函数。Abstract In digital signal processing, digital filter of the utmost importance. Digital filter is able to complete the pending tasks for discrete time systems, it can be achieved using the fini

3、te precision algorithm. Modern digital filter can be designed with software or dedicated digital signal processing hardware to achieve in two ways. Digital filter with the advantages of software, with the change of filter parameters, it is easy to change the filter to take the performance. Finite Im

4、pulse Response (FIR) digital filters because of the strict linear phase in digital communications, speech signal processing and adaptive signal processing and other fields has been widely used. This article describes the units of finite impulse response (FIR) digital filter theory, on the FIR digita

5、l filter design window function was studied. Window function in the design of FIR digital filter has a wide application. Conventional digital filter design method cumbersome and not intuitive, this paper MATLAB is a powerful scientific computing and graphical display of this advantage, and the windo

6、w method combines design theory to design FIR digital filter, the design results not only more intuitive, and improve the accuracy of the filter design to better achieve the desired results. Key words: Digital signal processing; FIR filter; window function.目 录引言 11 绪论1 数字信号处理的概述1数字信号处理的任务1数字信号处理系统21

7、.2 滤波器的的概念2滤波器的定义2滤波器的分类21.3 MATLAB简介22 数字滤波器的概述32.1 数字滤波器的分类与结构32.2.1 数字滤波器的分类3滤波器的基本结构42.2 数字滤波器的原理5数字滤波器的工作原理5滤波器的性能指标53 FIR数字滤波器设计63.1 窗函数设计法7设计方法7加窗处理产生的影响73.2 窗函数设计FIR数字滤波器的MATLAB实现8主程序框图及相关子程序89111315174 结论24谢 辞26参考文献27引言 随着信息时代和数字世界的到来，数字信号处理已成为当今时代极其重要的一门学科和技术领域。数字信号处理在通信、语音、图像、自动控制、雷达、军事、航

8、空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理应用中，数字滤波器十分重要并已获得广泛应用。数字滤波器是一个离散时间系统（按预定的算法，将输入离散时间信号转换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1/2抽样频率点呈镜像对称。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。由于滤波器的概念比较抽象，加上其数值计算又比较繁琐，所以需借助好的计算机软件来进行辅助设计，是

9、滤波器研究领域的一个发展趋势。这样的软件有很多种，其中最具代表性的就是MATLAB语言及其工具箱，其将一个优秀软件包的易用性与可靠性，通用性和专业性，以及一般目的的应用和高深的专业应用近乎完美地集成在一起.并凭借其功能强大，技术先进和应用之深广，使其逐渐成为国际性的计算标准，为世界各地超过20万名科学家和工程师所采用。利用MATLAB软件，可以使滤波器的设计和实现变得事半功倍，因此，本文借助MATLAB软件来进行应用研究。本文通过用窗函数法对数字滤波器的设计过程，充分说明在MATLAB环境中，滤波器的设计己变得异常简单和高效率，而且结果更为稳定、可靠。相信随着MATLAB软件在信息技术中的不断

10、推广和使用，将会使科技工程人员从各类琐碎的计算工作中解脱出来，有时间思考和研究更重要的问题，这必将为科技的进步起到巨大的推动作用。1 绪论1.1 数字信号处理的概述数字信号处理就是通过计算机或专用处理设备，用数值计算等数字的方式对数字序列或信号进行分析、滤波估计、检测及识别等加工处理，将信号变换成符合需要的某种形式。.1数字信号处理的任务信号是传递信息的载体，在信号处理中它通常是时间或者频率的函数。信号处理理论完成信号的滤波与检测、参数提取和估计、频谱分析等工作，使信号更便于使用，更容易识别。数字信号处理针对数字或符号表示的函数序列，用数值计算的方法，完成对数字信号的检测、滤波、参数估计等处理

11、，以得到用户需要的信号输出形式，如：通过信号滤波提取和增强信号的有用分量，削弱噪声分量。.2数字信号处理系统一个完整的模拟信号数字化处理系统如图所示。数字信号处理的研究对象要求输入和输出都是数字序列，因此信号处理系统并不要求包括框图的所有部分，如本身输入为数字序列的系统，则不需要抗混叠滤波器和A/D；最终只需要信号的数字参数时，不需要最后的D/A。Xa(n)Xa(n)X(n)y(t)ya(t)D/A转换器抗混叠滤波器A/D转换器数字信号处理器 数字信号处理系统框图1.2 滤波器的的概念滤波是信号处理中一种最基本而又重要的技术，利用滤波，可以从复杂的信号中提取所需的信号，抑制不需要的信号。所谓滤

12、波器，就是对已知激励，可以在时间域或频域产生规定响应的网络，要使滤波器能够提取有用信号，要求滤波器对信号与噪声有不同的增益，对有用信号尽量无失真放大，而对噪声尽量衰减。滤波器的定义所谓滤波器，就是对已知激励，可以在时间域或频域产生规定响应的网络，要使滤波器能够提取有用信号，要求滤波器对信号与噪声有不同的增益，对有用信号尽量无失真放大，而对噪声尽量衰减。滤波器的分类滤波器一般分为两类：模拟滤波器和数字滤波器；从功能上可分为低通滤波器（LPF）；高通滤波器（HPF）；带通滤波器（BPF）；带阻滤波器（BSF）；从实现方法上可分为IIR滤波器和FIR滤波器。本文主要讨论FIR数字滤波器。FIR滤波器

13、，即有限长单位冲激响应，有以下几个特点:(1) 系统的单位冲激响应h(n)在有限个n值处不为零； (2) 系统函数H(z)在|z|>0处收敛，极点全部z=0处因果系统）； (3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。1.3 MATLAB简介MATLAB 软件是由美国 Math Work公司于1984年推出，经过不断的完善与发展,如今已成为覆盖多个学科的国际公认的最优秀的数值计算仿真软件。MATLAB 具有强大的数值计算能力，许多复杂的计算问题只需短短几行代码就能在 MATLAB中实现。作为一个跨平台软件，MATLAB已推出

14、 Unix Windows Linux和 Mac等十多个操作系统的版本，大大方便了在不同操作系统平台下的研究工作。目前基于 Windows 系统的MATLAB7.0 它秉承以往版本的优点并且人机界面友好，非常容易使用。MATLAB 软件具有很强的开放性和适应性，在保持内核不变的情况下MATLAB 软件可以针对不同的学科推出相应的工具箱，目前已推出图像处理工具箱，信号处理工具箱，小波工具箱，神经网络工具箱以及通信工具箱等多个学科的专业工具箱，极大方便了不同学科的研究工作。国内有来越多的科研和技术人员认识到 MATLAB 的强大作用，并在不同领域内使用 MATLAB 来快速实现。充分发挥 MATL

15、AB语言的编程效率高、计算与图形出来能力非常强的特点，完成各种计算与图形绘制工作，可使模拟滤波器的设计变得简单快捷，为后继设计打下一个良好的基础。2 数字滤波器的概述2.1 数字滤波器的分类与结构数字滤波器是近几十年迅速发展起来的一门新技术和新学科，它在理论和实践中都有巨大的意义。数字滤波器理论在网络与滤波理论中占有越来越重要的地位。数字滤波器是数字信号处理中的基础部分，与模拟滤波器相比具有的优点有：精度高、可靠性高、灵活性高、便于大规模集成及多维过滤等。目前，在数字信号处理系统中，数字滤波器已经应用在地震法勘探、地震预报、雷达、声纳、通信、语言的分解与合成、空间技术、高空摄影、遥感、生物医疗

16、以及电视等很多方面。.1 数字滤波器的分类按时域特性，数字滤波器课分为无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。IIR滤波器的系统函数为： （21）该滤波器的脉冲响应h(n)在n=0到区间有值。FIR滤波器的系统函数为： （22）该滤波器的脉冲响应在n=0,1,2,3,4.,N-1的有限个点上有值。按频域来分数字滤波器可分为低通、高通、带通和带阻四种基本类型。本文重点讲解FIR滤波器。FIR滤波器的基本结构FIR滤波器的单位抽样响应为有限长度，一般采用非递归形式实现。通常的FIR数字滤波器有横截性和级联型两种。FIR滤波器实现的基本结构有：(1)FIR滤波器的横截型结构表示系

17、统输入输出关系的差分方程可写作： （23）  直接由差分方程得出的实现结构如图2.1所示：图2.1 横截型(直接型卷积型)若h(n)呈现对称特性，即此FIR滤波器具有线性相位，则可以简化加横截型结构，下面分情况讨论：图2.2 N为奇数时线形相位FIR滤波器实现结构 图2.3 N为偶数时线性相位FIR滤波器实现结构(2)FIR滤波器的级联型结构将H(z)分解成实系数二阶因子的乘积形式： （24） 这时FIR滤波器可用二阶节的级联结构来实现，每个二阶节用横截型结构实现。如图所示：图2.4 FIR滤波器的级联结构这种结构的每一节控制一对零点，因而在需要控制传输零点时可以

18、采用这种结构。2.2 数字滤波器的原理滤波器的工作原理在数字滤波中，我们主要讨论离散时间序列。如图2.5所示。设输入序列为，离散或数字滤波器对单位抽样序列的响应为。因在时域离散信号和系统中所起的作用相当于单位冲激函数在时域连续信号和系统中所起的作用。 数字滤波器,H(z)图2.5 数字滤波器原理数字滤波器的序列将是这两个序列的离散卷积，即 （25）同样，两个序列卷积的z变换等于个自z变换的乘积，即 （26）用代入上式，其中T为抽样周期，则得到 （27）式中和 分别为数字滤波器输入序列和输出序列的频谱，而为单位抽样序列响应的频谱。由此可见，输入序列的频谱经过滤波后，变为 ，按照的特点和我们处理信

19、号的目的，选取适当的使的滤波后的符合我们的要求。2.3.2滤波器的性能指标我们在进行滤波器设计时，需要确定其性能指标。一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波器特性为例，频率响应有通带、过渡带及阻带三个范围。在通带内：  1- AP| H(ej)| 1 |c在阻带中： |H(ej)| Ast st |c其中c 为通带截止频率, st为阻带截止频率，Ap为通带误差, Ast为阻带误差。 与模拟滤波器类似，数字滤波器按频率特性划分为低通、高通、带通、带阻等类型，由于数字滤波器的频率响应是周期性的，周期为2。各种理想数字滤波器的幅度频率响应如图所示：

20、图低通滤波器频率响应幅度特性的容限图 图2.7 各种理想数字滤波器的幅度频率响应3 FIR数字滤波器设计FIR滤波器的优点是可方便地实现线性相位。FIR滤波器单位冲激响应h(n)的特点：其单位冲激响应h(n)是有限长（1n N-1）,系统函数为： （31）在有限Z平面有(N-1)个零点，而它的(N-1)个极点均位于原点z=0处。FIR滤波器线性相位的特点：如果FIR滤波器的单位抽样响应h(n)为实数，而且满足以下任一条件：偶对称h(n)h(N-1-n)奇对称h(n)-h(N-1-n)其对称中心在n(N-1)/2处，则滤波器具有准确的线性相位。 窗函数设计法窗函数设计法的基本思想是用有限长单位冲

21、激响应数字滤波器去逼近希望的滤波器。假设希望逼近的滤波器的频率响应函数为，其单位冲激响应用表示。窗函数设计法就是截取为有线长的一断因果序列，并用合适的窗函数进行加权作为FIRDF的单位冲击响应h(n)。截取的长度和加权窗函数的类型都直接影响逼近精度。设计方法一般是先给定所要求的理想滤波器频率响应，由导出，我们知道理想滤波器的冲击响应是无限长的非因果序列，而我们要设计的是是有限长的FIR滤波器，所以要用有限长序列来逼近无限长序列。设计步骤如下：（1）根据要求求出： （32）（2）用有限长度的窗函数w(n)来截取，并进行加权处理，得到： （33）（3）将h(n)移位,使其成为因果序列：h1(n)=

22、 h(n-) （34）3.1.2加窗处理产生的影响加窗以后对理想低通滤波器的频率响应将产生什么样的影响呢？根据在时域是相乘关系,在频域则是卷积关系： （35） 其中, 为矩形窗谱, 是FIR滤波器频率响应。通过频域卷积过程看的幅度函数H()的起伏现象,可知,加窗处理后，对理想矩形的频率响应产生以下几点影响：(1)使理想频率特性不连续点处边沿加宽，形成一个过渡带，其宽度等于窗的频率响应的主瓣宽度。(2)在截止频率的两边的地方即过渡带的两边，出现最大的肩峰值，肩峰的两侧形成起伏振荡，其振荡幅度取决于旁瓣的相对幅度，而振荡的多少，则取决于旁瓣的多少。 (3)改变N，只能改变窗谱的主瓣宽度，改变的坐标

23、比例以及改变的绝对值大小，但不能改变主瓣与旁瓣的相对比例(此比例由窗函数的形状决定)。(4)对窗函数的要求a、窗谱主瓣尽可能窄，以获取较陡的过渡带；b、尽量减小窗谱的最大旁瓣的相对幅度；即能量集中于主瓣，使肩峰和波纹减小，增大阻带的衰减。3.1.3各种窗函数（1）矩形窗（rectangle window）在MATLAB中，实现矩形窗的函数为boxcar和rectwin，其调用格式为：w=boxcar(N)w=rectwin(N)（2）三角窗（triangular window）在MATLAB中，实现三角窗的函数为triang，其调用格式为：w=triang (N)（3）汉宁窗（hanning

24、window）在MATLAB中，实现汉宁窗的函数为hann，其调用格式为：w=hann(N)w=hann(N,sflag)（4）哈明窗（hamming window）在MATLAB中，实现哈明窗的函数为hamming，其调用格式为：w=hamming(N)w=hamming(N,sflag)（5）布莱克曼窗（blackman window）在MATLAB中，实现布莱克曼窗的函数为blackman，其调用格式为：w=blackman(N)w=blackman(N,sflag)几种窗函数的基本参数见下表表3.1。表3.1 各种窗函数的基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗

25、-134/N12三角形窗-258/N25汉宁窗-318/N44哈明窗-418/N53不莱克曼窗-5712/N743.2 窗函数设计FIR数字滤波器的MATLAB实现3.2.1主程序框图及相关子程序显示数字滤波器频域中的图形的子程序如下：functiondb,mag,pha,w=freqz_m2 (b,a)H,w=freqz (b,a,1000,'whole');H=(H(1:1:501)'w=(w(1:1:501)'mag=abs(H);db=20*log10 (mag+eps)/max(mag);pha=angle (H);开始读入技术指标计算hd(n)调用窗

26、函数子程序求w(n)调用子程序（函数）计算H(n)调用绘图子程序（函数）绘制H(n)幅度相位曲线结束图3.1 主程序框图计算h(n)= hd(n) w(n)主程序框图如图3.1所示。计算理想低通滤波器单位冲激响应的子函数：function hd=ideal_lp(Wc,N);alpha=(N-1)/2;n=0:(N-1);m=n-alpha+eps;hd=sin(Wc*m)./(pi*m);下面，利用窗函数和MATLAB子程序进行FIR数字低通滤波器的设计。其主程序如下所示：clear all;Wp=input('请输入频率Wp=');Ws =input('请输入频率W

27、s=');Ar =input('请输入阻带衰减Ar=');tr_width = Ws - Wp ;%过渡带宽度N=ceil (6.6*pi/tr_width)+1; %滤波器长度n=0:1:N-1;Wc=(Ws + Wp)/2;%理想低通滤波器的截止频率hd = ideal_lp (Wc ,N); %理想低通滤波器的单位冲激响应if (21<=Ar&&Ar<25); w_ham=( boxcar(N)' string='矩形窗','N=',num2str(N);else if (25<=Ar&am

28、p;&Ar<44); w_ham =(triang(N)' string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<53); w_ham =(hanning(N)' string='汉宁窗','N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' string='哈明窗','N=',num2str(

29、N);else if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end ;endh=hd.*w_ham;%截取得到实际的单位冲激响应db,mag,pha,w=freqz_m2 (h,1);%计算实际滤波器的幅度响应delta_w =2*pi/1000;Ap = -(min(db(1:1:Wp/delta_w+1)%实际通带波纹As = -round (max (db(Ws/delta_w+1:1:501)%实际阻带波纹subplot (3

30、,2,1); stem (n,hd); title ('理想脉冲响应hd(n)'); xlabel('n'); ylabel('hd(n)');subplot (3,2,2); stem (n,w_ham); ); xlabel('n'); ylabel('w(n)'); text(,1.3,string); subplot (3,2,3); stem(n,h); title ('实际脉冲响应h(n)'); xlabel('n'); ylabel('h(n)');su

31、bplot (3,2,4); plot(w,pha);title ('相频特性')axis(0,3.15,-4,4); xlabel('频率(pi)'); ylabel('相位()');subplot (3,2,5); plot (w/pi,db); title ('幅度特性(dB)');axis (0,1,-100,10); xlabel('频率(pi)'); ylabel('(分贝数)');subplot (3,2,6); plot (w,mag); title ('频率特性')

32、;axis (0,3.15, 0,); xlabel('频率(rad) ');ylabel('幅值');运行结果及分析：输入频率Wp=0.2*pi；Ws=0.4*pi； Ar=60后运行结果为：滤波器的长度为N=34，实际阻带衰减为52dB，实际的通带波动为0.0477，根据Ar=60程序自动选择了哈明窗。其相应的实验图形如下图3.2所示：图3.2 数字低通滤波器图3.2.3数字高通滤波器的窗函数设计数字高通滤波器有型和之分，其原理一样，在这我们只做型数字高通滤波器。计算型理想高通滤波器单位冲激响应的子函数：function hd = ideal_hp1 ( W

33、c ,N)alpha= (N-1)/2;n=0:1:N-1;m=n-alpha+eps;%hd = sin ( pi * m ) - sin ( Wc * m ) / (pi*m);hd=sin(pi*m)-sin(Wc*m)./(pi*m);其主程序如下：clear all;Wp =input('请输入频率Wp=');Ws=input('请输入频率Ws=');Ar =input('请输入阻带衰减Ar=');tr_width = Wp - Ws; %过渡带宽度N=ceil (6.2*pi/tr_width) ;%滤波器长度n=0:1:N-1;Wc

34、=(Ws + Wp)/2;%理想滤波器的截止频率hd = ideal_hp1 (Wc ,N);%理想滤波器的单位冲激响应if (21<=Ar&&Ar<25)；w_ham=( boxcar(N)' string='矩形窗','N=',num2str(N);else if (25<=Ar&&Ar<44)；w_ham =(triang(N)'string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<

35、;53); 分析果w_ham =(hanning(N)' string='汉宁窗','N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' string='哈明窗','N=',num2str(N);else if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end;en

36、d;h=hd .*w_ham; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2 (h,1); %计算实际滤波器的幅度响应 delta_w =2*pi/1000;Ap = -(min(db(Wp/delta_w+1:1:501);%实际通带波纹As = -round (max (db(1:1:Ws/delta_w+1);%实际阻带波纹subplot (3,2,1); stem (n,hd); title ('理想脉冲响应hd(n)');); xlabel('n'); ylabel('hd(n)');subplot (3,2,

37、2); stem (n,w_ham); );xlabel('n'); ylabel('w(n)'); text(,1.3,string);subplot (3,2,3); stem(n,h); title ('实际脉冲响应h(n)');); xlabel('n'); ylabel('h(n)');subplot (3,2,4); plot(w,pha); title ('相频特性');axis(0,3.15,-4,4); xlabel('频率(pi)'); ylabel('相

38、位()');subplot (3,2,5); plot (w/pi,db); title ('幅度特性(dB)');axis (0,1,-100,10); xlabel('频率(pi)'); ylabel('(分贝数)');subplot (3,2,6); plot (w,mag); title ('频率特性');axis (0,3.15, 0,); xlabel('频率(rad) '); ylabel('幅值');运行及结果分析：输入频率Wp=0.6*pi；Ws=0.4*pi； Ar=75

39、后运行结果为：滤波器的长度为N=31，实际的阻带衰减为26dB，实际的通带波动为0.4418，根据Ar=75程序自动选择了布莱克曼窗。其相应的实验图形如下图3.3所示：图3.3 数字高通滤波器图数字带通滤波器分型和型之分，这里我们做型，其理想高通滤波器单位冲激响应的子函数为：function hd = ideal_bp1 ( Wcl,Wch,N)alpha= (N-1)/2;n=0:1:N-1;m=n-alpha+eps;%hd=sin (Wch*m)-sin (Wcl*m)/ (pi*m);hd=sin(Wch*m)-sin(Wcl*m)./(pi*m);其主程序如下：clear all;W

40、ph =input('请输入频率Wph=');Wpl=input('请输入频率Wpl=');Wsl =input('请输入频率Wsl=');Wsh=input('请输入频率Wsh=');Ar =input('请输入阻带衰减Ar=');tr_width = min(Wpl - Wsl),(Wsh -Wph);%过渡带宽度N=ceil (6.1*pi/tr_width);%滤波器长度n=0:1:N-1;Wcl=(Wsl + Wpl)/2;%理想滤波器的截止频率Wch=(Wsh + Wph)/2;hd = ideal_b

41、p2 (Wcl ,Wch ,N);%理想滤波器的单位冲激响应if (21<=Ar&&Ar<25); w_ham=( boxcar(N)' string='矩形窗','N=',num2str(N);else if (25<=Ar&&Ar<44); w_ham =(triang(N)' string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<53); w_ham =(hanning(N)

42、9; string='汉宁窗','N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' string='哈明窗','N=',num2str(N);if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end; end;h=hd.*w_ham;%截取得到实际的单位冲激响应db,mag

43、,pha,w=freqz_m2 (h,1);%计算实际滤波器的幅度响应delta_w =2*pi/1000;Ap = -(min(db( Wpl/delta_w+1:1:Wph/delta_w+1);%实际通带波纹As = -round (max (db(Wsh/delta_w+1:1:501);%实际阻带波纹subplot (3,2,1); stem (n,hd); title ('理想脉冲响应hd(n)');); xlabel('n'); ylabel('hd(n)'); text(,1.3,string);subplot (3,2,2);

44、stem (n,w_ham); ); xlabel('n'); ylabel('w(n)'); text(,1.3,string);subplot (3,2,3); stem(n,h); title ('实际脉冲响应h(n)'); xlabel('n'); ylabel('h(n)');subplot (3,2,4); plot(w,pha); title ('相频特性');axis(0,3.15,-4,4); xlabel('频率(rad)'); ylabel('相位(ra

45、d)');subplot (3,2,5); plot (w/pi,db); title ('幅度特性(dB)');axis (0,1,-150,10); xlabel('频率(单位:pi)');ylabel('(分贝数)'); subplot (3,2,6)plot (w,mag); title ('频率特性')axis (0,3.15, 0,); xlabel('频率(rad) '); ylabel('幅值');运行及结果分析：输入频率Wph=0.6*pi，Wpl=0.4*pi，Wsh=0

46、.8*pi，Wsl=0.2*pi， Ar=23后运行结果为：滤波器的长度为N=31，实际阻带衰减27dB，实际的通带波动为1.4064，根据Ar=23，程序自动选择了矩形窗。其相应的实验图形如下图3.4所示：图3.4 带通滤波器图 计算理想数字带阻滤波器的单位冲激响应的子函数为：function hd = ideal_bs ( Wcl,Wch,N)alpha= (N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=sin(Wcl*m)+sin(pi*m)-sin( Wch*m) ./(pi*m);其主程序如下：clear all;Wph =input('请输入频率Wph

47、=');Wpl=input('请输入频率Wpl=');Wsl =input('请输入频率Wsl=');Wsh=input('请输入频率Wsh=');Ar =input('请输入阻带衰减Ar=');tr_width = min(Wsl - Wpl),( Wph - Wsh);%过渡带宽度N=ceil (6.2*pi/tr_width);%滤波器长度n=0:1:N-1;Wcl=(Wsl + Wpl)/2;%理想滤波器的截止频率Wch=(Wsh + Wph)/2;hd = ideal_bs (Wcl ,Wch ,N);%理想滤波

48、器的单位冲激响应if (21<=Ar&&Ar<25); w_ham=( boxcar(N)' string='矩形窗','N=',num2str(N);else if (25<=Ar&&Ar<44); w_ham =(triang(N)' string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<53); w_ham =(hanning(N)' string='汉宁窗'

49、;,'N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' string='哈明窗','N=',num2str(N);else if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end; end;h=hd.*w_ham;%截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2 (

50、h,1);%计算实际滤波器的幅度响应delta_w =2*pi/1000;Ap = -(min(db( 1:1:Wpl/delta_w+1);实际通带波纹As = -round (max (db(Wsl/delta_w+1:1: Wsh/delta_w+1);%实际阻带波纹subplot (3,2,1); stem (n,hd); title ('理想脉冲响应hd(n)');); xlabel('n'); ylabel('hd(n)');subplot (3,2,2); stem (n,w_ham); );xlabel('n');

51、 ylabel('w(n)'); text(,1.3,string);subplot (3,2,3); stem(n,h); title ('实际脉冲响应h(n)');); xlabel('n'); ylabel('h(n)');subplot (3,2,4); plot(w,pha); title ('相频特性');axis(0,3.15,-4,4); xlabel('频率(rad)'); ylabel('相位()');subplot (3,2,5); plot (w/pi,db)

52、; title ('幅度特性(dB)');axis (0,1,-80,10); xlabel('频率(pi)'); ylabel('(分贝数)');subplot (3,2,6); plot (w,mag); title ('频率特性');axis (0,3.15, 0,); xlabel('频率(rad) '); ylabel('幅值');运行结果及分析：输入频率Wph=0.8*pi，Wpl=0.2*pi，Wsh=0.6*pi，Wsl=0.4*pi， Ar=23后运行结果为：滤波器的长度为N=31

53、，实际阻带衰减27dB，实际的通带波动为1.4064，根据Ar=23，程序自动选择了矩形窗。其相应的实验图形如下图3.5所示：图3.5 带阻滤波器图3.2.6程序拓展及四种滤波器比较 利用循环语句将高通、低通、带通、带阻四种FIR数字滤波器合成为一个可以连续循环执行的程序。其主程序如下：clear all;i=0; b=1; while(b);temp=menu('选择窗函数类型','1=低通','2=高通','3=带通','4=带阻');menu1=1,2,3,4;MM=menu1(temp);if MM=1Wp

54、=input('请输入Wp=');Ws =input('请输入Ws='); Ar =input('请输入阻带衰减Ar=');tr_width = Ws - Wp ;%过渡带宽度N=ceil (6.6*pi/tr_width)+1;%滤波器长度n=0:1:(N-1);Wc=(Ws + Wp)/2;%理想滤波器的截止频率hd = ideal_lp (Wc ,N);%理想滤波器的单位冲激响应if (21<=Ar&&Ar<25); w_ham=( boxcar(N)' string='矩形窗','

55、;N=',num2str(N);else if (25<=Ar&&Ar<44); w_ham =(triang(N)' string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<53); w_ham =(hanning(N)' string='汉宁窗','N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' s

56、tring='哈明窗','N=',num2str(N);else if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end;end;h=hd.*w_ham;%截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2 (h,1);%计算实际滤波器的幅度响应delta_w =2*pi/1000;Ap = -(min(db(1:1:Wp/delta_w+1);%实际通带波纹As = -round(max

57、(db(Ws/delta_w+1:1:501);%实际阻带波纹endif MM=2Wp =input('请输入Wp=');Ws=input('请输入Ws=');Ar =input('请输入阻带衰减Ar=');tr_width = Wp - Ws; %过渡带宽度N=ceil (6.2*pi/tr_width); %滤波器长度n=0:1:(N-1);Wc=(Ws + Wp)/2; %理想滤波器的截止频率hd = ideal_hp1 (Wc ,N); %理想滤波器的单位冲激响应if (21<=Ar&&Ar<25); w_ha

58、m=( boxcar(N)' string='矩形窗','N=',num2str(N);else if (25<=Ar&&Ar<44); w_ham =(triang(N)' string='三角窗 ','N=',num2str(N);else if (44<=Ar&&Ar<53);w_ham =(hanning(N)' string='汉宁窗','N=',num2str(N);else if (53<=Ar&&Ar<74); w_ham =(hamming(N)' string='哈明窗','N=',num2str(N);else if (Ar>=74); w_ham =(blackman(N)' string='布莱克曼窗','N=',num2str(N);end; end; end; end; end;h=hd .*w_ham; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2 (h,1); %计算实际滤波器的幅度响应delta_w =2*pi