MATLAB设计数字带通FIR滤波器的几种窗函数的比较.DOC.pdf
中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 1 课题介绍 通过平时所学、书本内容和网络上了解的知识，我们小组讨论决定对 MATLAB 结合窗函数设计一个数字带通FIR 滤波器工程中涉及的几个窗函数方 法优劣好坏用图形形象的进行比较，以MATLAB 软件为工具探讨出一个失真小 高效的设计方法。 数字滤波器是一种用来过滤时间离散信号的数字系统，通过对 抽样数据进行数学处理来达到频域滤波的目的。与IIR 滤波器相比， FIR 的实现 是非递归的，总是稳定的；更重要的是， FIR 滤波器在满足幅频响应要求的同时， 可以获得严格的线性相位特性。根据FIR 滤波器的原理，提出了FIR 滤波器的 窗函数设计法，给出了在MATLAB 环境下，用窗函数法设计FIR 滤波器的过程 和设计实例。 通过利用不同的窗函数方法设计FIR 滤波器，对所设计的滤波器进 行分析比较， 得出各种方法设计的滤波器的优缺点及其不同的使用场合，从而可 以在设计滤波器时能够正确的选择FIR 数字滤波器的窗函数的选取及设计方法。 数字滤波技术 数字滤波，就是通过一定的计算或判断程序减少干扰在有用信号中的比重， 所以故实质上是一种程序滤波。 与此对应的就是模拟滤波， 模拟滤波主要无源绿 波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算 放大器等需要外接电源组成的） ，其目的是将信号中的噪音和干扰滤去或者将希 望得到的频率信号滤出为我所用。数字滤波的出现克服了模拟滤波的很多不足， 具有以下优点： 1是用程序实现的，不需要增加硬设备，所以可靠性高，稳定性好。 2可以对频率很低的信号实现滤波，克服了模拟滤波的缺陷。 3可以根据信号的不同，采用不同的滤波方法或参数，具有灵活、方便、 功能强的特点。 本文主要对 FIR 滤波器加以介绍。 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 2 FIR 滤波器 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;通带峰值 起伏： p40dB 设计函数的选取 MATLAB 信号处理工具箱提供了基于加窗的线性相位FIR 滤波器设计函数 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 3 fir1和 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') 函数参数说明如下： 1n 表示滤波器的阶数 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 种窗函数性能比较 窗类型旁瓣峰值主瓣峰值最小阻带衰减 矩形窗13dB 4 /M 21dB 三角窗25dB 8 /M 25dB 汉宁窗31dB 8 /M 44dB 海明窗41dB 8 /M 53dB 凯塞窗57dB 12 /M 74dB 常用窗函数绘图比较： 在 MATLAB 中运行以下代码： 代码 2-1： n=50; x=1:50; 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 5 juxing=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-1 5种窗函数绘图比较 设计步骤 实际利用窗函数法进行FIR 滤波器设计时， 依据所给的技术指标一般需要经 过以下几个步骤进行设计： 1 给定理想的频率响应函数Hd(e jw )及技术指标； 2 求出理想的单位抽样响应hd(n); 3 根据阻带衰减选择窗函数w(n) 4 根据过渡带宽度确定N 值; 5 求出所设计的 FIR 滤波器的单位抽样响应； 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 6 6 计算频率响应，验算指标是否满足要求。 设计方法 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 的长度） 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 10 3）devs用于指定各个频带输出滤波器的频率响应与其期望幅值 之间的最大输出误差或偏差，长度与a 相等。 4）fs 为信号的采样频率。 利用该函数修改代码得到凯塞窗设计fir 滤波器的代码如下： 代码 2-3 fs=20000; %设定采样频率 fp1=4500;fp2=6500; %第一截止频率 fs1=3000;fs2=7500; %第二截止频率 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 命令行中输入 fdatool 命令即可打开滤波器设计工具 利用 FDATool 设计的凯塞窗 fir 滤波器幅频响应和相频响应如下： 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 11 图 2-9 FDATool 设计的凯塞窗 fir 滤波器幅频响应 图 2-10 FDATool 设计的凯塞窗 fir 滤波器相频响应 滤波器滤波性能测试 由 MATLAB 模拟生成含有不同频率的数字信号，然后利用设计的滤波器对 数字信号进行滤波， 为方便观察，模拟生成的信号只含有包含在阻带的两个频率 （2000Hz，8000Hz）一个包含于通带的频率（5000Hz） ，测试代码如下： 代码 2-4： fs=20000; t=0:1/fs:2; 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; %取一段信号 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 12 subplot(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,-0.5,0.5); ylabel('滤波后的信号 ');xlabel('时间'); 运行结果： 图 3-1 滤波器滤波性能测试波形 从原始型号和滤波后的信号对比可以看出，在用设计的滤波器进行滤波后信 号基本成单一频率的正弦波，滤波结果令人满意。 滤波器时延测量 滤波器的时延特性也是滤波器的性能指标之一，为了观测所设计滤波器的时 延， 可以利用一个还有一单一冲击的数字信号序列通过滤波器观测滤波器的输出 波形，相对于滤波性能的测量只是在原始信号产生方法上有所不同，这里仅给出 运行结果： 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 13 图 3-2 滤波器时延性能测试波形 从图上看，滤波器仍然有一定的时延，这也是fir 滤波器的缺点之一，并且 时延与滤波器滤波性能相矛盾， 滤波器的时延随滤波器阶数的升高而增加，实际 设计时要综合考虑两方面的因素。 滤波器稳定性测量 FIR 滤波器的一个突出优点便是它的稳定性，从信号与系统的理论可知，当 一个系统的 Z 域的传递函数的极点都在单位圆内时系统是稳定的，FIR 滤波器的 传递函数的分母为1，即 FIR 滤波器一定是稳定的， MATLAB 中也提供了专门 用于绘制零极点图的函数： zplane()，通过该函数可以轻松绘制系统的零极点图， 以凯泽窗设计的 FIR 滤波器为例，绘制其零极点图如下： 图 3-3 凯塞窗 FIR 滤波器零极点图 从该零极点图可以看出， FIR 设计的滤波器是属于无极点的系统，因而系统 一定是稳定的。 从上面对滤波器三方面的测试可以看出FIR 滤波器的特点， 第一个便是其相 位曲线，在通带内，相频曲线一定为直线，二是FIR 滤波器由于阶数比较高，因 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 14 而有一定的时延，三是FIR 滤波器一定是稳定的。 收获及总结 当我们小组讨论决定对MATLAB 结合窗函数设计一个数字带通FIR 滤波器 工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较的时候，我们小组 认为这还是还是很有难度的，因为我们虽然做了实验但是对MATLAB 中有关滤 波器的设计知识了解还是比较少的， 但正是这种有难度的题目才能提高我们的能 力，我们小组决定这个探究课题后，从网上获得了很多有关MATLAB 设计滤波 器的资料， 加上我们对书本的翻阅， 我们的设计思路渐渐明朗， 经过我们小组的 协调合作， 我们最终还是弄明白了FIR 滤波器的工作原理， 通过查阅网络、 书本 等相关资料， 我们明白了 FIR 滤波器设计相关的一些函数的使用方法，并了解了 其中的优劣，在设计过程中，我们小组三人通力合作，分工明确，体会到了小组 合作的乐趣，不但学到了知识，而且加深了感情，这是难能可贵的，也感谢王老 师这匠心独运的作业方式，让我们受益匪浅。 信息 08-1 班 徐小龙 汉宁窗带通滤波器设计源代码： 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; %求截止频率 中国矿业大学信息与电气工程学院数字与信号处理探究学习徐小龙 15 hd=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); 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); %求得理想冲击响应
