MATLAB设计数字带通FIR滤波器的几种窗函数的比较

中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙1课题介绍通过平时所学、书本内容和网络上了解的知识，我们小组讨论决定对MATLAB结合窗函数设计一个数字带通FIR滤波器工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较，以MATLAB软件为工具探讨出一个失真小高效的设计方法。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。与IIR滤波器相比，FIR的实现是非递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。根据FIR滤波器的原理，提出了FIR滤波器的窗函数设计法，给出了在MATLAB环境下，用窗函数法设计FIR滤波器的过程和设计实例。通过利用不同的窗函数方法设计FIR滤波器，对所设计的滤波器进行分析比较，得出各种方法设计的滤波器的优缺点及其不同的使用场合，从而可以在设计滤波器时能够正确的选择FIR数字滤波器的窗函数的选取及设计方法。数字滤波技术数字滤波，就是通过一定的计算或判断程序减少干扰在有用信号中的比重，所以故实质上是一种程序滤波。与此对应的就是模拟滤波，模拟滤波主要无源绿波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算放大器等需要外接电源组成的），其目的是将信号中的噪音和干扰滤去或者将希望得到的频率信号滤出为我所用。数字滤波的出现克服了模拟滤波的很多不足，具有以下优点：1．是用程序实现的，不需要增加硬设备，所以可靠性高，稳定性好。2．可以对频率很低的信号实现滤波，克服了模拟滤波的缺陷。3．可以根据信号的不同，采用不同的滤波方法或参数，具有灵活、方便、功能强的特点。本文主要对FIR滤波器加以介绍。中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙2FIR滤波器FIR滤波器是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位冲激响应是有限的，没有输入到输出的反馈，是稳定的系统。FIR滤波器具有以下主要优点：1.FIR滤波器具有准确的线性相位；2.FIR滤波器永远稳定；3.FIR滤波器设计方法一般是线性的；4.FIR滤波器在硬件上具有更高的运行效率；5.FIR滤波器启动传输时间只需要有限时间。FIR滤波器的主要缺点有：1.FIR滤波器为达到同样的性能要求需要比IIR滤波器高得多的阶数；2.相应的FIR滤波器的时延比同等性能的IIR滤波器高很多。FIR滤波器的硬件实现主要有数字集成芯片，DSP芯片FIR滤波器，可编程FIR滤波器，后两者的实际方法主要通过MATLAB软件进行设计，其设计方法多样，形式灵活，能够满足各种要求，并且不受数字集成芯片规格的限制。FIR滤波器的设计方法中最常用的是窗函数法。滤波器设计要求利用MATLAB仿真软件系统结合窗函数法设计一个数字带通FIR滤波器。要求：分别使用矩形窗、三角形窗、汉明窗、布莱克曼窗、凯泽窗进行设计，并输出滤波器的频率特性。然后进行比较得出优缺点。参数：采样频率为20kHz；通带边缘频率：fp1=4.5kHz,fp2=6.5kHz;通带峰值起伏：αp1dB;阻带边缘频率：fs1=3k,fs2=7.5k;最小阻带衰减：As40dB设计函数的选取MATLAB信号处理工具箱提供了基于加窗的线性相位FIR滤波器设计函数中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙3fir1和fir2，fir1函数的调用格式为：b=fir1(n,Wn)b=fir1(n,Wn,'ftype')b=fir1(n,Wn,window)b=fir1(n,Wn,'ftype',window)b=fir1(…..,'normalization')函数参数说明如下：1．n表示滤波器的阶数2．'ftype'表示所设计滤波器的类型：3．'high'表示高通滤波器4．'stop'表示带阻滤波器5．'DC-1'表示多通带滤波器，第一频带为通带6．'DC-0'表示多通带滤波器，第一频带为阻带；默认时为低通或带通滤波器；7．'window'为窗函数，是长度为n+1的列向量，默认时函数自动取Hamming窗。该函数实现加窗的线性相位FIR滤波器设计，可以设计标准低通、带通、高通和带阻滤波器窗函数构造MATLAB工具箱已经提供了各种窗函数的构造函数，因而窗函数的构造十分方便，下面给出几种常用窗函数的构造方法：1．矩形窗：利用w=boxcar(n)的形式得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，它的元素由窗函数的值组成。‘w=boxcar(n)’等价于‘w=ones(1,n)’.2．三角窗：利用w=triang(n)的形式得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，它的元素由窗函数的值组成。w=triang(N-2)等价于bartlett(N)。3．汉宁窗：利用w=hanning(n)得到窗函数，其中n为窗函数的长度，而返中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙4回值w为一个n阶的向量，包含了窗函数的n个系数。4．布拉克曼窗：利用w=blackman(n)得到窗函数，其中n为窗函数的长度，而返回值w为一个n阶的向量，包含了窗函数的n个系数。它的主瓣宽度是矩形窗主瓣宽度的3倍，为12*pi/N，但是它的最大旁瓣值比主瓣值低57dB。5．凯塞窗：利用w=kaiser(n,beta)的形式得到窗函数。窗函数设计条件在使用窗函数设计滤波器时要满足以下两个条件：1．窗谱主瓣尽可能地窄，从而可以获得较陡峭的过渡带；2．尽量减少窗谱的最大旁瓣的相对幅度，及尽可能是能量集中于主瓣，减少峰肩和波纹，进而增加阻带的衰减。窗函数设计条件实际工程常用的窗函数有五种，即矩形窗、三角窗、汉宁窗、海明窗和凯塞窗。这些窗函数之间的性能比较如下表所示。5种窗函数性能比较窗类型旁瓣峰值主瓣峰值最小阻带衰减矩形窗13dB4π/M21dB三角窗25dB8π/M25dB汉宁窗31dB8π/M44dB海明窗41dB8π/M53dB凯塞窗57dB12π/M74dB常用窗函数绘图比较：在MATLAB中运行以下代码：代码2-1：n=50;x=1:50;中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙5juxing=boxcar(n);%构造矩形窗sanjiao=triang(n);%构造三角窗hanming=hamming(n);%构造汉宁窗bulaikeman=blackman(n);%构造布莱克曼窗kaize=kaiser(n);%构造凯塞窗plot(x,juxing,'b.',x,sanjiao,'gx',x,hanming,'r+',x,bulaikeman,'cd',x,kaize,'k*');legend('矩形窗','三角窗','汉明窗','布莱克曼窗','凯塞窗');运行结果如图2-1所示：图2-15种窗函数绘图比较设计步骤实际利用窗函数法进行FIR滤波器设计时，依据所给的技术指标一般需要经过以下几个步骤进行设计：1．给定理想的频率响应函数Hd(ejw)及技术指标；2．求出理想的单位抽样响应hd(n);3．根据阻带衰减选择窗函数w(n)4．根据过渡带宽度确定N值;5．求出所设计的FIR滤波器的单位抽样响应；中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙66．计算频率响应，验算指标是否满足要求。设计方法MATLAB作为一款优秀的数值计算软件，本身就内置了丰富的函数，其中便有用于通信仿真的一系列函数，并且MATLAB中还集成了通信设计的工具箱，不管是内置的函数，还是通信工具箱，均有专用于滤波器设计的工具，常用的主要有用函数法设计和用通信工具箱设计，下面分别予以介绍。依据设计步骤编写此种方法不依赖MATLAB中的滤波器设计函数，而是依据FIR滤波器的设计步骤自己求解理想滤波器的冲击响应，然后用窗函数对冲击响应进行截短，从而得到FIR滤波器，由于低通滤波器设计较为简单，因而可以通过两个低通的理想冲击响应函数相减得到理想带通的冲击响应，再通过窗函数对其截短，从而得到实际可行的FIR滤波器，这里以用汉明窗进行截短实现带通FIR滤波器，关键代码如下：完整代码在后面代码2-2：wp1=0.3*pi;ws1=0.45*pi;wp2=0.65*pi;ws2=0.75*pi;tr_width=ws1-wp1;%求过渡带宽度M=ceil(6.6*pi/tr_width)+1;%求得所需窗函数的长度n=[0:1:M-1];wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;%求截止频率hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);%求得理想带通的冲击响应w_ham=(hamming(M))';%得到长度为M的汉宁窗h=hd.*w_ham;%利用窗函数截短程序运行结果：中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙7图2-2汉宁窗带通滤波器利用MATLAB自带函数设计利用原理说明介绍的fir1函数进行设计，这种设计方法只需要给出滤波器的阶数，截止频率，窗函数等参数，MATLAB即可自行完成设计，并可通过freqz函数查看滤波器的幅频响应和相频响应，已验证滤波器是否满足设计要求，下面给出利用fir1函数设计的不同窗函数的数字带通滤波器。1.利用矩形窗进行设计代码2-2：fs=20000;%设定采样频率fp1=4500;fp2=6500;%第一截止频率fs1=3000;fs2=7500;%第二截止频率As=40;%最小阻带衰减Ws1=(fp1+fs1)/fs;Ws2=(fp2+fs2)/fs;%截止频率归一化处理w=(fp1-fs1)/fs;%求归一化过渡带M=ceil((As-7.95)/(14.36*w))%计算所需滤波器的阶数juxing=boxcar(M+1);%生成长度为M+1的矩形窗boxb=fir1(M,[Ws1,Ws2],juxing);%生成矩形窗设计的fir滤波器freqz(boxb,1,fs,fs);%绘制幅频和相频响应曲线中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙8运行结果：图2-3矩形窗fir滤波器幅频和相频响应曲线从幅频响应上看，通带基本无波纹，阻带中波纹较大，因而阻带较不理想，相频响应曲线在通带内为直线，效果较好，信号失真小。2.利用三角窗进行设计利用三角窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用triang（）函数生成三角窗，程序运行结果如下：图2-4三角窗设计的fir滤波器幅频和相频响应曲线3.利用汉宁窗进行设计利用汉宁窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用hamming（）函数生成三角窗，程序运行结果如下：中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙9图2-5汉宁窗设计的fir滤波器幅频和相频响应曲线4.利用布莱克曼窗进行设计利用布莱克曼窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用blackman（）函数生成三角窗，程序运行结果如下：图2-6布莱克曼窗设计的fir滤波器幅频和相频响应曲线5.利用凯塞窗进行设计利用凯塞窗进行设计时，滤波器的参数可以用以下函数求得：调用函数[n,wn,bta,ftype]=kaiserord（f，a，dev，fs）参数：1）f为对应的归一化频率2）a为由f指定的各个频带上的幅值向量，一般只有0和1，和f长度关系为（2*a的长度）—2=（f的长度）中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙103）devs用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等。4）fs为信号的采样频率。利用该函数修改代码得到凯塞窗设计fir滤波器的代码如下：代码2-3fs=20000;%设定采样频率fp1=4500;fp2=6500;%第一截止频率fs1=3000;fs2=7500;%第二截止频率[n,wn,bta,ftype]=kaiserord([fs1,fp1,fp2,fs2],...[0