1、中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 1课题介绍通过平时所学、书本内容和网络上了解的知识，我们小组讨论决定对MATLAB 结合窗函数设计一个数字带通 FIR 滤波器工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较，以 MATLAB 软件为工具探讨出一个失真小高效的设计方法。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。与 IIR 滤波器相比，FIR的实现是非递归的，总是稳定的；更重要的是，FIR 滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。根据 FIR 滤波器的原理，提出了FIR 滤波器的窗函数设计
2、法，给出了在 MATLAB 环境下，用窗函数法设计 FIR滤波器的过程和设计实例。通过利用不同的窗函数方法设计 FIR 滤波器，对所设计的滤波器进行分析比较，得出各种方法设计的滤波器的优缺点及其不同的使用场合，从而可以在设计滤波器时能够正确的选择 FIR 数字滤波器的窗函数的选取及设计方法。数字滤波技术数字滤波，就是通过一定的计算或判断程序减少干扰在有用信号中的比重，所以故实质上是一种程序滤波。与此对应的就是模拟滤波，模拟滤波主要无源绿波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算放大器等需要外接电源组成的） ，其目的是将信号中的噪音和干扰滤去或者将希望得到的频率信号滤
3、出为我所用。数字滤波的出现克服了模拟滤波的很多不足，具有以下优点：1是用程序实现的，不需要增加硬设备，所以可靠性高，稳定性好。2可以对频率很低的信号实现滤波，克服了模拟滤波的缺陷。3可以根据信号的不同，采用不同的滤波方法或参数，具有灵活、方便、功能强的特点。本文主要对 FIR 滤波器加以介绍。中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 2FIR 滤波器FIR 滤波器是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位冲激响应是有限的，没有输入到输出的反馈，是稳定的系统。FIR 滤波器具有以下主要优点：1.FIR 滤波器具有准确的线
4、性相位；2.FIR 滤波器永远稳定；3.FIR 滤波器设计方法一般是线性的；4.FIR 滤波器在硬件上具有更高的运行效率；5.FIR 滤波器启动传输时间只需要有限时间。FIR 滤波器的主要缺点有：1.FIR 滤波器为达到同样的性能要求需要比 IIR 滤波器高得多的阶数；2.相应的 FIR 滤波器的时延比同等性能的 IIR 滤波器高很多。FIR 滤波器的硬件实现主要有数字集成芯片，DSP 芯片 FIR 滤波器，可编程FIR 滤波器，后两者的实际方法主要通过 MATLAB 软件进行设计，其设计方法多样，形式灵活，能够满足各种要求，并且不受数字集成芯片规格的限制。FIR 滤波器的设计方法中最常用的是
5、窗函数法。滤波器设计要求利用 MATLAB 仿真软件系统结合窗函数法设计一个数字带通 FIR 滤波器。要求：分别使用矩形窗、三角形窗、汉明窗、布莱克曼窗、凯泽窗进行设计，并输出滤波器的频率特性。然后进行比较得出优缺点。参数：采样频率为 20kHz；通带边缘频率：fp1=4.5kHz,fp2=6.5kHz;通带峰值起伏：p40dB中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 3设计函数的选取MATLAB 信号处理工具箱提供了基于加窗的线性相位 FIR 滤波器设计函数fir1 和 fir2，fir1 函数的调用格式为：b=fir1(n,Wn)b=fir1(n,Wn,ftype)b=
6、fir1(n,Wn,window)b=fir1(n,Wn,ftype,window)b=fir1(,normalization)函数参数说明如下：1n 表示滤波器的阶数2ftype表示所设计滤波器的类型：3high表示高通滤波器4stop表示带阻滤波器5DC-1 表示多通带滤波器，第一频带为通带6DC-0 表示多通带滤波器，第一频带为阻带；默认时为低通或带通滤波器；7window为窗函数，是长度为 n+1 的列向量，默认时函数自动取 Hamming窗。该函数实现加窗的线性相位 FIR 滤波器设计，可以设计标准低通、带通、高通和带阻滤波器窗函数构造MATLAB 工具箱已经提供了各种窗函数的构造函
7、数，因而窗函数的构造十分方便，下面给出几种常用窗函数的构造方法：1矩形窗：利用 w=boxcar(n)的形式得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，它的元素由窗函数的值组成。 w=boxcar(n)等价于w=ones(1,n).中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 42三角窗：利用 w=triang(n)的形式得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，它的元素由窗函数的值组成。 w=triang(N-2)等价于 bartlett(N)。3汉宁窗：利用 w=hanning(n)得到窗函数，其中 n 为窗函数
8、的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。4布拉克曼窗：利用 w=blackman(n)得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。它的主瓣宽度是矩形窗主瓣宽度的 3 倍，为 12*pi/N，但是它的最大旁瓣值比主瓣值低 57dB。5凯塞窗：利用 w=kaiser(n,beta)的形式得到窗函数。窗函数设计条件在使用窗函数设计滤波器时要满足以下两个条件：1 窗谱主瓣尽可能地窄，从而可以获得较陡峭的过渡带；2尽量减少窗谱的最大旁瓣的相对幅度，及尽可能是能量集中于主瓣，减少峰肩和波纹，进而增加阻带的衰减。窗函数
9、设计条件实际工程常用的窗函数有五种，即矩形窗、三角窗、汉宁窗、海明窗和凯塞窗。这些窗函数之间的性能比较如下表所示。5 种窗函数性能比较窗类型 旁瓣峰值 主瓣峰值 最小阻带衰减矩形窗 13dB 4/M 21dB三角窗 25dB 8/M 25dB汉宁窗 31dB 8/M 44dB海明窗 41dB 8/M 53dB凯塞窗 57dB 12/M 74dB常用窗函数绘图比较：中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 5在 MATLAB 中运行以下代码：代码 2-1：n=50;x=1:50;juxing=boxcar(n); %构造矩形窗sanjiao=triang(n); %构造三角窗
10、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-1 5 种窗函数绘图比较设计步骤实际利用窗函数法进行 FIR 滤波器设计时，依据所给的技术指标一般需要经过以下几个步骤进行设计：1 给定理想的频率响应函数 Hd(ejw)及技术指标；2 求出理
11、想的单位抽样响应 hd(n);中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 63 根据阻带衰减选择窗函数 w(n)4 根据过渡带宽度确定 N 值;5 求出所设计的 FIR 滤波器的单位抽样响应；6 计算频率响应，验算指标是否满足要求。设计方法MATLAB 作为一款优秀的数值计算软件，本身就内置了丰富的函数，其中便有用于通信仿真的一系列函数，并且 MATLAB 中还集成了通信设计的工具箱，不管是内置的函数，还是通信工具箱，均有专用于滤波器设计的工具，常用的主要有用函数法设计和用通信工具箱设计，下面分别予以介绍。依据设计步骤编写此种方法不依赖 MATLAB 中的滤波器设计函数，而是
12、依据 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-
13、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 函数查看滤波器的
14、幅频响应和相频响应，已验证滤波器是否满足设计要求，下面给出利用 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 的矩形窗b
15、oxb=fir1(M,Ws1,Ws2,juxing); %生成矩形窗设计的 fir 滤波器中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 8freqz(boxb,1,fs,fs); %绘制幅频和相频响应曲线运行结果：图 2-3 矩形窗 fir 滤波器幅频和相频响应曲线从幅频响应上看，通带基本无波纹，阻带中波纹较大，因而阻带较不理想，相频响应曲线在通带内为直线，效果较好，信号失真小。 2. 利用三角窗进行设计利用三角窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用 triang（）函数生成三角窗，程序运行结果如下：图 2-4 三角窗设计的 fir 滤波器幅频和相频响应曲线
16、3. 利用汉宁窗进行设计利用汉宁窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用 hamming（）函数生成三角窗，程序运行结果如下：中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 9图 2-5 汉宁窗设计的 fir 滤波器幅频和相频响应曲线4. 利用布莱克曼窗进行设计利用布莱克曼窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用 blackman（）函数生成三角窗，程序运行结果如下：图 2-6 布莱克曼窗设计的 fir 滤波器幅频和相频响应曲线5. 利用凯塞窗进行设计利用凯塞窗进行设计时，滤波器的参数可以用以下函数求得：调用函数n,wn,bta,ftype=k
17、aiserord（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; %第
18、二截止频率n,wn,bta,ftype=kaiserord(fs1,fp1,fp2,fs2,.0,1,0,0.01 0.1087 0.01,fs) %求滤波器参数b=fir1(n,wn,ftype,kaiser(n+1,bta); %生成 fir 滤波器freqz(b,1,fs,fs); %绘制幅频和相频响应曲线程序运行结果：图 2-7 凯塞窗设计的 fir 滤波器幅频和相频响应曲线利用 MATLAB 工具箱设计MATLAB 信号处理工具箱中已集成了用于滤波器设计和分析的工具：FDATool，利用它可以实现数字滤波器的可视化设计与分析，操作简单方便，在 MATLAB 命令行中输入 fdatoo
19、l 命令即可打开滤波器设计工具利用 FDATool 设计的凯塞窗 fir 滤波器幅频响应和相频响应如下：中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 11图 2-9 FDATool 设计的凯塞窗 fir 滤波器幅频响应图 2-10 FDATool 设计的凯塞窗 fir 滤波器相频响应滤波器滤波性能测试由 MATLAB 模拟生成含有不同频率的数字信号，然后利用设计的滤波器对数字信号进行滤波，为方便观察，模拟生成的信号只含有包含在阻带的两个频率（2000Hz， 8000Hz）一个包含于通带的频率（ 5000Hz） ，测试代码如下：代码 2-4：fs=20000;t=0:1/fs:2
20、; x=sin(2*pi*2000*t)+sin(2*pi*5000*t).+sin(2*pi*8000*t); %生成混合信号xo=filter(b,2,xn); %用滤波器对信号进行滤波figure;nn=5000:5100; %取一段信号中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 12subplot(211);tt=nn/fs;plot(tt,x(nn); %绘制原始信号axis(0.25,0.255,-4,4);ylabel(原始信号 );xlabel(时间);subplot(212);plot(tt,xo(nn); %绘制滤波后的信号axis(0.25,0.255,
21、-0.5,0.5);ylabel(滤波后的信号 );xlabel(时间);运行结果：图 3-1 滤波器滤波性能测试波形从原始型号和滤波后的信号对比可以看出，在用设计的滤波器进行滤波后信号基本成单一频率的正弦波，滤波结果令人满意。滤波器时延测量滤波器的时延特性也是滤波器的性能指标之一，为了观测所设计滤波器的时延，可以利用一个还有一单一冲击的数字信号序列通过滤波器观测滤波器的输出波形，相对于滤波性能的测量只是在原始信号产生方法上有所不同，这里仅给出运行结果：中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 13图 3-2 滤波器时延性能测试波形从图上看，滤波器仍然有一定的时延，这也是
22、fir 滤波器的缺点之一，并且时延与滤波器滤波性能相矛盾，滤波器的时延随滤波器阶数的升高而增加，实际设计时要综合考虑两方面的因素。滤波器稳定性测量FIR 滤波器的一个突出优点便是它的稳定性，从信号与系统的理论可知，当一个系统的 Z 域的传递函数的极点都在单位圆内时系统是稳定的， FIR 滤波器的传递函数的分母为 1，即 FIR 滤波器一定是稳定的，MATLAB 中也提供了专门用于绘制零极点图的函数：zplane()，通过该函数可以轻松绘制系统的零极点图，以凯泽窗设计的 FIR 滤波器为例，绘制其零极点图如下：图 3-3 凯塞窗 FIR 滤波器零极点图从该零极点图可以看出，FIR 设计的滤波器是
23、属于无极点的系统，因而系统一定是稳定的。从上面对滤波器三方面的测试可以看出 FIR 滤波器的特点，第一个便是其相位曲线，在通带内，相频曲线一定为直线，二是 FIR 滤波器由于阶数比较高，中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 14因而有一定的时延，三是 FIR 滤波器一定是稳定的。收获及总结当我们小组讨论决定对 MATLAB 结合窗函数设计一个数字带通 FIR 滤波器工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较的时候，我们小组认为这还是还是很有难度的，因为我们虽然做了实验但是对 MATLAB 中有关滤波器的设计知识了解还是比较少的，但正是这种有难度的题目才能提高
24、我们的能力，我们小组决定这个探究课题后，从网上获得了很多有关 MATLAB设计滤波器的资料，加上我们对书本的翻阅，我们的设计思路渐渐明朗，经过我们小组的协调合作，我们最终还是弄明白了 FIR 滤波器的工作原理，通过查阅网络、书本等相关资料，我们明白了 FIR 滤波器设计相关的一些函数的使用方法，并了解了其中的优劣，在设计过程中，我们小组三人通力合作，分工明确，体会到了小组合作的乐趣，不但学到了知识，而且加深了感情，这是难能可贵的，也感谢王老师这匠心独运的作业方式，让我们受益匪浅。信息 08-1 班 徐小龙汉宁窗带通滤波器设计源代码：wp1=0.3*pi; ws1=0.45*pi;wp2=0.6
25、5*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; %求截止频率中国矿业大学信息与电气工程学院数字与信号处理探究学习 徐小龙 15hd=ideal_lp(wc2,M)-ideal_lp(wc1,M); %求得理想带通的冲击响应w_ham=(hamming(M); %得到长度为 M 的汉明窗h=hd .* w_ham; %利用窗函数截短%绘图部分subplot(1,1,1)subplot(2,2,1)
26、; stem(n,hd,.); title(理想冲击响应)axis(0 M-1 -0.3 0.4); xlabel(n); ylabel(hd(n)subplot(2,2,2); stem(n,w_ham,.);title(汉明窗)axis(0 M-1 0 1.1); xlabel(n); ylabel(w(n)subplot(2,2,3); stem(n,h,.);title(实际冲击响应)axis(0 M-1 -0.3 0.4); xlabel(n); ylabel(h(n)subplot(2,2,4); plot(w/pi,db);title(幅频响应);axis(0 1 -100 10); xlabel(f); ylabel(dB)理想冲击响应求解函数：function hd = ideal_lp(wc,M);alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha + eps; %加入一个无穷小量，避免除零hd = sin(wc*m) ./ (pi*m); %求得理想冲击响应
