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

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

《MATLAB设计数字带通FIR滤波器的几种窗函数的比较要点.docx》由会员分享，可在线阅读，更多相关《MATLAB设计数字带通FIR滤波器的几种窗函数的比较要点.docx（16页珍藏版）》请在冰豆网上搜索。

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

课题介绍

通过平时所学、书本内容和网络上了解的知识，我们小组讨论决定对ＭATLAB结合窗函数设计一个数字带通ＦIR滤波器工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较，以MATLＡB软件为工具探讨出一个失真小高效的设计方法。

数字滤波器是一种用来过滤时间离散信号的数字系统,通过对抽样数据进行数学处理来达到频域滤波的目的。

与IＩR滤波器相比,ＦＩR的实现是非递归的，总是稳定的;更重要的是，FIR滤波器在满足幅频响应要求的同时,可以获得严格的线性相位特性。

根据ＦIR滤波器的原理，提出了FIR滤波器的窗函数设计法,给出了在MＡTLＡB环境下，用窗函数法设计FIＲ滤波器的过程和设计实例。

通过利用不同的窗函数方法设计FＩＲ滤波器,对所设计的滤波器进行分析比较，得出各种方法设计的滤波器的优缺点及其不同的使用场合,从而可以在设计滤波器时能够正确的选择FIR数字滤波器的窗函数的选取及设计方法。

数字滤波技术

　　数字滤波,就是通过一定的计算或判断程序减少干扰在有用信号中的比重，所以故实质上是一种程序滤波。

与此对应的就是模拟滤波,模拟滤波主要无源绿波（直接用电阻、电容、电感等不外接电源的元件组成的）与有源滤波（如运算放大器等需要外接电源组成的）,其目的是将信号中的噪音和干扰滤去或者将希望得到的频率信号滤出为我所用。

数字滤波的出现克服了模拟滤波的很多不足，具有以下优点：

1．是用程序实现的，不需要增加硬设备，所以可靠性高,稳定性好。

２.可以对频率很低的信号实现滤波,克服了模拟滤波的缺陷。

３.可以根据信号的不同,采用不同的滤波方法或参数,具有灵活、方便、功能强的特点。

本文主要对FIＲ滤波器加以介绍。

FIR滤波器

FIＲ滤波器是数字信号处理系统中最基本的元件,它可以在保证任意幅频特性的同时具有严格的线性相频特性,同时其单位冲激响应是有限的,没有输入到输出的反馈,是稳定的系统。

ＦIR滤波器具有以下主要优点：

1.ＦIR滤波器具有准确的线性相位;

2．FIＲ滤波器永远稳定；

３.FIR滤波器设计方法一般是线性的；

4.ＦIR滤波器在硬件上具有更高的运行效率；

5．FIR滤波器启动传输时间只需要有限时间。

FＩR滤波器的主要缺点有:

１.FIR滤波器为达到同样的性能要求需要比IIR滤波器高得多的阶数；

2.相应的FIＲ滤波器的时延比同等性能的IIＲ滤波器高很多。

FIＲ滤波器的硬件实现主要有数字集成芯片，ＤSP芯片FIR滤波器，可编程FIR滤波器,后两者的实际方法主要通过MATＬAB软件进行设计，其设计方法多样,形式灵活，能够满足各种要求，并且不受数字集成芯片规格的限制。

FIR滤波器的设计方法中最常用的是窗函数法。

　滤波器设计要求

利用MATLＡB仿真软件系统结合窗函数法设计一个数字带通FＩR滤波器。

要求：

分别使用矩形窗、三角形窗、汉明窗、布莱克曼窗、凯泽窗进行设计,并输出滤波器的频率特性。

然后进行比较得出优缺点。

参数：

采样频率为20kＨz;通带边缘频率：

fp1=4．5kHz,fp2=6.5kHz；通带峰值起伏:

αp<1dB;阻带边缘频率：

fｓ1=3k,ｆs２＝7.5ｋ;最小阻带衰减:

As>40dＢ

设计函数的选取

MＡTLAB信号处理工具箱提供了基于加窗的线性相位ＦIR滤波器设计函数fiｒ１和fｉr2,fｉr１函数的调用格式为:

b=fir1（n,Wn）

b=fｉｒ1（n，Wn,'ftｙpｅ'）

b=fiｒ1（ｎ,Ｗn,wｉndow）

b=fir1（n,Ｗn,'ftｙｐｅ＇,winｄow）

b=fir1（…..,'normａｌｉzatiｏn'）

函数参数说明如下：

1.n表示滤波器的阶数

2.'ftｙpe＇表示所设计滤波器的类型：

3.'high＇表示高通滤波器

4．＇stｏp'表示带阻滤波器

5．'DC-1'表示多通带滤波器，第一频带为通带

６.'DC-0'表示多通带滤波器,第一频带为阻带;默认时为低通或带通滤波器；

7．＇wｉndｏw'为窗函数，是长度为n+1的列向量，默认时函数自动取Hammiｎg窗。

该函数实现加窗的线性相位FＩR滤波器设计,可以设计标准低通、带通、高通和带阻滤波器

窗函数构造

MＡＴLAＢ工具箱已经提供了各种窗函数的构造函数，因而窗函数的构造十分方便,下面给出几种常用窗函数的构造方法：

1．矩形窗:

利用w=boxcar（n）的形式得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,它的元素由窗函数的值组成。

‘w=boxcaｒ（ｎ）’等价于‘w=ｏnes（1,n）’．

2．三角窗:

利用ｗ＝trianｇ（n）的形式得到窗函数，其中n为窗函数的长度,而返回值ｗ为一个n阶的向量，它的元素由窗函数的值组成。

w=ｔriang（Ｎ－2）等价于barｔletｔ（N）。

３.汉宁窗:

利用w=hanning（n）得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量，包含了窗函数的n个系数。

4．布拉克曼窗:

利用w＝blackｍan（n）得到窗函数,其中ｎ为窗函数的长度,而返回值ｗ为一个n阶的向量,包含了窗函数的n个系数。

它的主瓣宽度是矩形窗主瓣宽度的３倍,为12*pｉ/N,但是它的最大旁瓣值比主瓣值低５7dB。

５．凯塞窗：

利用ｗ=kａiser（n,betａ）的形式得到窗函数。

窗函数设计条件

在使用窗函数设计滤波器时要满足以下两个条件：

1．窗谱主瓣尽可能地窄,从而可以获得较陡峭的过渡带;

2．尽量减少窗谱的最大旁瓣的相对幅度,及尽可能是能量集中于主瓣，减少峰肩和波纹，进而增加阻带的衰减。

　窗函数设计条件

实际工程常用的窗函数有五种,即矩形窗、三角窗、汉宁窗、海明窗和凯塞窗。

这些窗函数之间的性能比较如下表所示。

5种窗函数性能比较

窗类型

旁瓣峰值

主瓣峰值

最小阻带衰减

矩形窗

13dＢ

4π/M

21ｄB

三角窗

25dB

8π/M

2５dB

汉宁窗

31dＢ

8π／M

4４dB

海明窗

41dB

8π/M

5３dB

凯塞窗

57dB

１２π/M

74dB

常用窗函数绘图比较:

在MＡTLＡB中运行以下代码：

代码2-1：

n=50;

x=1：

50;

juxing=boｘｃar（ｎ）;　ﻩ　%构造矩形窗

sanjｉaｏ=trianｇ（ｎ）;ﻩ％构造三角窗

hanmiｎｇ＝hamminｇ（ｎ）；　ﻩ%构造汉宁窗

bｕｌaiｋeｍaｎ＝blackｍａn（n）;　%构造布莱克曼窗

kａiｚｅ＝kaｉser（n）;　　ﻩ　％构造凯塞窗

plot（ｘ,juxing,'b．'，x,ｓａnjiao,'ｇx＇，ｘ,hanmiｎg,'r+',x,ｂulaikemａn,＇cｄ'，x,kａize,＇k＊＇）;

ｌeｇend（'矩形窗',＇三角窗',＇汉明窗'，'布莱克曼窗'，＇凯塞窗'）;

运行结果如图2-1所示:

图2-15种窗函数绘图比较

设计步骤

实际利用窗函数法进行FＩR滤波器设计时，依据所给的技术指标一般需要经过以下几个步骤进行设计：

1．给定理想的频率响应函数Hd（ejw）及技术指标；

2．求出理想的单位抽样响应hd（n）;

3．根据阻带衰减选择窗函数w（n）

4．根据过渡带宽度确定N值；

5．求出所设计的FＩR滤波器的单位抽样响应;

6．计算频率响应，验算指标是否满足要求。

设计方法

ﻩMATLＡB作为一款优秀的数值计算软件，本身就内置了丰富的函数,其中便有用于通信仿真的一系列函数,并且MＡＴLAB中还集成了通信设计的工具箱，不管是内置的函数,还是通信工具箱，均有专用于滤波器设计的工具，常用的主要有用函数法设计和用通信工具箱设计，下面分别予以介绍。

依据设计步骤编写

ﻩ此种方法不依赖MATLＡB中的滤波器设计函数，而是依据FＩＲ滤波器的设计步骤自己求解理想滤波器的冲击响应，然后用窗函数对冲击响应进行截短,从而得到FＩR滤波器，由于低通滤波器设计较为简单，因而可以通过两个低通的理想冲击响应函数相减得到理想带通的冲击响应，再通过窗函数对其截短，从而得到实际可行的ＦIR滤波器,这里以用汉明窗进行截短实现带通ＦIＲ滤波器，关键代码如下:

完整代码在后面

代码2-2：

wｐ1=０.3＊ｐi；　wｓ１=0．45*pi;

ｗp２=0.6５＊pi;　ws2=0.75*pi;

tr_width=ws1-wｐ1；ﻩ%求过渡带宽度

M=ceｉl（6.6*ｐｉ/tｒ_width）+1;ﻩﻩﻩ%求得所需窗函数的长度

n=[0:

１:

M-1];

wc1=（ws1＋wp1）／２;　wc2=（ws2＋ｗp２）/2;ﻩﻩ%求截止频率

hd=ideaｌ_lp（wｃ2,Ｍ）-ideaｌ_lｐ（ｗｃ1,M）;ﻩ%求得理想带通的冲击响应

w_ham=（hamming（Ｍ））';ﻩﻩﻩ%得到长度为M的汉宁窗

h=ｈｄ.*ｗ_ham;ﻩﻩﻩﻩﻩ%利用窗函数截短

程序运行结果：

图2-2ﻩ汉宁窗带通滤波器

利用MAＴLAB自带函数设计

利用原理说明介绍的fiｒ１函数进行设计，这种设计方法只需要给出滤波器的阶数，截止频率，窗函数等参数，MATＬAＢ即可自行完成设计，并可通过freqｚ函数查看滤波器的幅频响应和相频响应,已验证滤波器是否满足设计要求，下面给出利用ｆir１函数设计的不同窗函数的数字带通滤波器。

1.利用矩形窗进行设计

代码2-2:

ｆｓ=２000０;　　　　　　　%设定采样频率

fp1=450０;ｆp２＝６500;　　　　　　　ﻩ％第一截止频率

fｓ１=3000；fs２＝７500；　　　　　%第二截止频率

As=4０;　　　　　　　　　　　　　　　%最小阻带衰减

Ws1＝（fp１＋fｓ１）/fs;Ws２=（fp2+fs2）／ｆs;　　%截止频率归一化处理

w=（fp1-fs1）／fｓ;　　　　　　　　ﻩ%求归一化过渡带

M＝ceil（（Aｓ-7.95）/（14．3６＊w））　　　　ﻩ%计算所需滤波器的阶数

ｊｕxｉｎg=boxcar（M+1）;　　　　　ﻩ%生成长度为M+1的矩形窗

bｏxb=ｆir１（M,［Ｗs1,Ws2]，juｘing）；ﻩ％生成矩形窗设计的fir滤波器

ｆrｅｑz（boxb,１,fs,fs）;　　　　　　ﻩﻩ％绘制幅频和相频响应曲线

运行结果:

图2-3矩形窗fｉr滤波器幅频和相频响应曲线

从幅频响应上看,通带基本无波纹,阻带中波纹较大,因而阻带较不理想，相频响应曲线在通带内为直线,效果较好，信号失真小。

2.利用三角窗进行设计

利用三角窗进行设计时,原理与矩形窗基本相同，只不过生成窗函数时采用triａｎg（）函数生成三角窗,程序运行结果如下：

图２-4三角窗设计的fir滤波器幅频和相频响应曲线

3.利用汉宁窗进行设计

利用汉宁窗进行设计时,原理与矩形窗基本相同，只不过生成窗函数时采用hamming（）函数生成三角窗，程序运行结果如下:

图2-５汉宁窗设计的fｉr滤波器幅频和相频响应曲线

4.利用布莱克曼窗进行设计

利用布莱克曼窗进行设计时，原理与矩形窗基本相同,只不过生成窗函数时采用blackｍaｎ（）函数生成三角窗，程序运行结果如下：

图２-6布莱克曼窗设计的fｉr滤波器幅频和相频响应曲线

5.利用凯塞窗进行设计

利用凯塞窗进行设计时，滤波器的参数可以用以下函数求得:

调用函数[n,ｗn,bta,ftｙpｅ]＝kaisｅrｏrd（f,a,dev,fs）

参数:

1）ｆ为对应的归一化频率

2）a为由f指定的各个频带上的幅值向量，一般只有0和1,和f长度关系为（2*a的长度）—2=（f的长度）

３）ｄevs用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差,长度与a相等。

4）fs为信号的采样频率。

利用该函数修改代码得到凯塞窗设计ｆｉｒ滤波器的代码如下：

代码２-３

fs=20０00；　　　　　　　%设定采样频率

fp1=4500;fｐ２=６5００;　　　　ﻩ％第一截止频率

ｆｓ1=300０;fs２=75００;　　％第二截止频率

[ｎ,wn,btａ,ｆｔype]=kａｉseroｒd（[fs1,ｆp1,fｐ2,fs2],..．

　[０,1,0],［0.010.10870.0１］，fs）　ﻩ%求滤波器参数

b=fir1（n,wn,ｆtype，kaiser（n+1，bta））;　ﻩﻩ%生成fiｒ滤波器

freqz（b，1，fs，ｆｓ）;　　　　%绘制幅频和相频响应曲线

程序运行结果：

图2－7凯塞窗设计的ｆｉｒ滤波器幅频和相频响应曲线

利用MATＬＡＢ工具箱设计

ＭAＴLAB信号处理工具箱中已集成了用于滤波器设计和分析的工具:

ＦDATｏol,利用它可以实现数字滤波器的可视化设计与分析，操作简单方便，在MATＬＡB命令行中输入fｄatool命令即可打开滤波器设计工具

利用FＤAＴoｏl设计的凯塞窗fir滤波器幅频响应和相频响应如下:

图2－9FＤＡTool设计的凯塞窗ｆiｒ滤波器幅频响应

图2-1０FDATｏol设计的凯塞窗ｆir滤波器相频响应

滤波器滤波性能测试

由MATLAＢ模拟生成含有不同频率的数字信号，然后利用设计的滤波器对数字信号进行滤波,为方便观察，模拟生成的信号只含有包含在阻带的两个频率（2００0Hz，800０Ｈz）一个包含于通带的频率（5０00Hz）,测试代码如下:

代码2－4:

fｓ=200０0；

ｔ=０:

1/ｆs:

2;　　　　　

ｘ＝ｓｉn（２*pi*2000*ｔ）＋sin（2*pｉ*50０0*t）...

　＋ｓin（2*pi*８000＊t）;　　　　　%生成混合信号

xo=filter（ｂ,2,xｎ）；　　　　　　　％用滤波器对信号进行滤波

figuｒｅ;

nｎ＝50０0：

5100;　　　　　　　　%取一段信号

ｓubploｔ（２１1）;

tt=nn/ｆs;

pｌot（tt,x（ｎn））;ﻩﻩﻩﻩ%绘制原始信号

ａxｉs（[0.25，0.255,-4,4］）;

ｙlabｅｌ（'原始信号'）；xlabel（＇时间'）;

subｐlot（21２）；

pｌｏt（tｔ，xo（nn））;ﻩﻩﻩ　%绘制滤波后的信号

ａxｉs（[0.25,0.2５5，－0.5,0.5］）；

yｌabel（'滤波后的信号'）;xlａｂel（＇时间'）;

运行结果：

图3-１ﻩ滤波器滤波性能测试波形

从原始型号和滤波后的信号对比可以看出,在用设计的滤波器进行滤波后信号基本成单一频率的正弦波,滤波结果令人满意。

　滤波器时延测量

滤波器的时延特性也是滤波器的性能指标之一，为了观测所设计滤波器的时延，可以利用一个还有一单一冲击的数字信号序列通过滤波器观测滤波器的输出波形，相对于滤波性能的测量只是在原始信号产生方法上有所不同,这里仅给出运行结果：

图3-2ﻩ滤波器时延性能测试波形

从图上看,滤波器仍然有一定的时延，这也是fir滤波器的缺点之一,并且时延与滤波器滤波性能相矛盾,滤波器的时延随滤波器阶数的升高而增加,实际设计时要综合考虑两方面的因素。

　滤波器稳定性测量

FIR滤波器的一个突出优点便是它的稳定性，从信号与系统的理论可知，当一个系统的Z域的传递函数的极点都在单位圆内时系统是稳定的，FIＲ滤波器的传递函数的分母为1,即FIR滤波器一定是稳定的,ＭAＴLＡB中也提供了专门用于绘制零极点图的函数:

ｚplane（）,通过该函数可以轻松绘制系统的零极点图,以凯泽窗设计的FIR滤波器为例,绘制其零极点图如下:

图3-3凯塞窗FIＲ滤波器零极点图

从该零极点图可以看出,ＦIR设计的滤波器是属于无极点的系统,因而系统一定是稳定的。

从上面对滤波器三方面的测试可以看出FIＲ滤波器的特点，第一个便是其相位曲线，在通带内,相频曲线一定为直线，二是ＦIR滤波器由于阶数比较高,因而有一定的时延，三是ＦIR滤波器一定是稳定的。

　　　　收获及总结

当我们小组讨论决定对ＭATLAB结合窗函数设计一个数字带通FIＲ滤波器工程中涉及的几个窗函数方法优劣好坏用图形形象的进行比较的时候,我们小组认为这还是还是很有难度的，因为我们虽然做了实验但是对MATLAB中有关滤波器的设计知识了解还是比较少的,但正是这种有难度的题目才能提高我们的能力，我们小组决定这个探究课题后,从网上获得了很多有关MATLAB设计滤波器的资料,加上我们对书本的翻阅，我们的设计思路渐渐明朗，经过我们小组的协调合作，我们最终还是弄明白了FIＲ滤波器的工作原理，通过查阅网络、书本等相关资料，我们明白了FＩR滤波器设计相关的一些函数的使用方法，并了解了其中的优劣，在设计过程中,我们小组三人通力合作，分工明确,体会到了小组合作的乐趣，不但学到了知识,而且加深了感情,这是难能可贵的，也感谢王老师这匠心独运的作业方式，让我们受益匪浅。

　　　　　信息０8－1班徐小龙

汉宁窗带通滤波器设计源代码:

ｗｐ１＝0.3＊pi;ws1=0.45*pｉ;

wp2=0.65*ｐi;ws2=０.7５*pi;

ｔr_ｗidｔh＝ｗs1-wｐ1;ﻩﻩﻩ％求过渡带宽度

M=ceil（６.6*pi/tr＿wiｄth）+1;ﻩﻩﻩﻩ％求得所需窗函数的长度

n=［０:

1:

M-1];

wｃ1=（wｓ1+wp１）/2;　ｗｃ2＝（ws２+ｗp２）/２;ﻩ％求截止频率

hd＝ｉdeal_lp（wｃ2,Ｍ）－ideal＿lp（wc1,M）;ﻩﻩ%求得理想带通的冲击响应

ｗ_ｈaｍ=（hamｍing（M））＇；ﻩﻩﻩ%得到长度为M的汉明窗

h＝hｄ.＊ｗ_ham；ﻩﻩﻩﻩ%利用窗函数截短

％绘图部分

suｂｐｌｏt（1，1,１）

subplot（２，2，1）;stem（n,hｄ，＇.'）;tｉtｌｅ（'理想冲击响应＇）

axｉs（[0　Ｍ-1　-0．３0.4］）;xlabel（＇ｎ＇）;ylaｂel（'hｄ（n）＇）

subplot（2,2，２）;sｔｅｍ（ｎ,w_hａm,'.'）;titｌｅ（＇汉明窗'）

aｘｉs（［0M-１0　1.１]）;　xlaｂeｌ（'n'）;yｌabel（'ｗ（n）＇）

subｐlｏt（2,2，3）；ｓｔeｍ（n，h，'.'）；ｔitle（'实际冲击响应'）

ａxｉs（[0M-１　-0.30.4］）；xｌabeｌ（'ｎ'）;　ylabeｌ（'h（ｎ）'）

ｓｕbpｌｏt（2,２,4）;ploｔ（ｗ/pi，db）;ｔｉtｌe（＇幅频响应'）;

ａxis（[01　-１00　10］）;xlabel（'f＇）;ylａbel（'dB'）

理想冲击响应求解函数：

fuｎctionhｄ＝ideａl＿lｐ（wｃ,Ｍ）;

alpha＝（M-1）/2;

n　=[0:

１：

（M-1）];

m＝n－alpｈa+　eps;ﻩﻩ%加入一个无穷小量，避免除零

hd＝sin（wｃ*m）./　（pｉ*ｍ）;ﻩﻩﻩﻩ%求得理想冲击响应