基于窗函数法的FIR带阻滤波器的软件实现.docx
《基于窗函数法的FIR带阻滤波器的软件实现.docx》由会员分享,可在线阅读,更多相关《基于窗函数法的FIR带阻滤波器的软件实现.docx(18页珍藏版)》请在冰豆网上搜索。
基于窗函数法的FIR带阻滤波器的软件实现
基于窗函数法的FIR带阻滤波器的软件实现
————————————————————————————————作者:
————————————————————————————————日期:
课程设计任务书
2010—2011学年第一学期
专业:
通信工程学号:
姓名:
课程设计名称:
数字信号处理课程设计
设计题目:
基于窗函数法的FIR带阻滤波器的软件实现
完成期限:
自2011年1月3日至2011年1月9日共1周
一.设计目的
1.巩固所学的理论知识。
2.提高综合运用所学理论知识独立分析和解决问题的能力。
3.更好地将理论与实践相结合。
4.掌握信号分析与处理的基本方法与实现。
5.熟练使用MATLAB语言进行编程实现。
二.设计内容
编写MATLAB程序实现FIR带阻滤波器的设计。
指标如下:
下通带截至频率
;上通带截止频率
;阻带下限频率
;阻带上限频率
;通带最大衰减
;阻带最小衰减
.
三.设计要求
1、根据指标要求选择合适的窗函数进行设计;
2、绘出
及其幅频响应并分析设计结果,验证所设计的滤波器是否满足要求。
四.设计条件
计算机、MATLAB语言环境
五、参考资料
[1]《数字信号处理》(第三版),丁玉美,高西全.西安电子科技大学出版社,2000。
[2]《MATLAB及在电子信息课程中的应用》,陈怀堔,吴大正,高西全.电子工业出版社,2006。
[3]《MATLAB7.0从入门到精通》,求是科技。
人民邮电出版社,2006。
[4]《数字信号处理(第三版)》学习指导,高西全,丁玉美。
西安科技大学出版社,2001。
指导教师(签字):
教研室主任(签字):
批准日期:
年月日
摘要
无限长脉冲数字滤波器的设计方法只考虑了幅度特性,没有考虑相位特性,所设的滤波器一般是某种确定的非线性相位特性.有限脉冲响应(FIR)滤波器在保证了幅度特性满足技术要求的同时,很容易做到有严格的线性相位特性。
本课题利用MATLAB软件实现。
MATLAB是“矩阵实验室”(MATrixLABoratoy)的缩写,是一种科学计算软件,它使用方便,输入简捷,运算高效,内容丰富,因此利用MATLAB软件,通过一系列较为系统的函数法,根据已知的技术指标,就可以设计出满足要求的滤波器。
关键词:
MATLAB;窗函数;FIR带阻数字滤波器;线性相位
1课题描述
现代图像、语声、数据通信对线性相位的要求是普遍的。
正是此原因,使得具有线性相位的FIR数字滤波器得到大力发展和广泛应用。
在实际进行数字信号处理时,往往需要把信号的观察时间限制在一定的时间间隔内,只需要选择一段时间信号对其进行分析。
这样,取用有限个数据,即将信号数据截断的过程,就等于将信号进行加窗函数操作。
而这样操作以后,常常会发生频谱分量从其正常频谱扩展开来的现象,即所谓的“频谱泄漏”。
当进行离散傅立叶变换时,时域中的截断是必需的,因此泄漏效应也是离散傅立叶变换所固有的,必须进行抑制。
而要对频谱泄漏进行抑制,可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣,或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。
而在后面的FIR滤波器的设计中,为获得有限长单位取样响应,需要用窗函数截断无限长单位取样响应序列。
另外,在功率谱估计中也要遇到窗函数加权问题。
由此可见,窗函数加权技术在数字信号处理中的重要地位。
2MATLAB简介
MATLAB是“矩阵实验室”(MATrixLABoratoy)的缩写,是一种科学计算软件,主要适用于矩阵运算及控制和信息处理领域的分析设计,它使用方便,输入简捷,运算高效,内容丰富,因此很多专家在自己擅长的领域用它编写了许多专门的MATLAB工具包,由于MATLAB功能的不断扩展,所以是科学研究中最常用必不可少的工具。
MATLAB由一系列工具组成。
这些工具方便用户使用MATLAB的函数和文件,其中许多工具采用的是图形用户界面。
包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。
随着MATLAB的商业化以及软件本身的不断升级,MATLAB的用户界面也越来越精致,更加接近Windows的标准界面,人机交互性更强,操作更简单。
而且新版本的MATLAB提供了完整的联机查询、帮助系统,极大的方便了用户的使用。
简单的编程环境提供了比较完备的调试系统,程序不必经过编译就可以直接运行,而且能够及时地报告出现的错误及进行出错原因分析。
MATLAB一个高级的矩阵/阵列语言,它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步,也可以先编写好一个较大的复杂的应用程序(M文件)后再一起运行.新版本的MATLAB语言是基于最为流行的C++语言基础上的,因此语法特征与C++语言极为相似,而且更加简单,更加符合科技人员对数学表达式的书写格式。
使之更利于非计算机专业的科技人员使用。
而且这种语言可移植性好、可拓展性极强,这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。
3窗函数设计法原理
数字滤波器可以理解为是一个计算程序或算法,将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列,并在转化过程中,使信号按预定的形式变化.数字滤波器有多种分类,根据数字滤波器冲激响应的时域特征,可将数字滤波器分为两种,即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。
IIR数字滤波器具有无限宽的冲激响应,与模拟滤波器相匹配。
所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。
FIR数字滤波器的单位脉冲响应是有限长序列。
它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题,设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。
因此设计FIR滤波器的方法之一可以从时域出发,截取有限长的一段冲击响应作为H(z)的系数,冲击响应长度N就是系统函数H(z)的阶数.只要N足够长,截取的方法合理,总能满足频域的要求。
一般这种时域设计、频域检验的方法要反复几个回合才能成功.要设计一个线性相位的FIR数字滤波器,首先要求理想频率响应
。
是w的周期函数,周期为
,可以展开成傅氏级数:
=
(3—1)
其中
是与理想频响对应的理想单位抽样响应序列。
但不能用来作为设计FIRDF用的h(n),因为
一般都是无限长、非因果的,物理上无法实现。
为了设计出频响类似于理想频响的滤波器,可以考虑用h(n)来近似
。
窗函数的基本思想:
先选取一个理想滤波器(它的单位抽样响应是非因果、无限长的),再截取(或加窗)它的单位抽样响应得到线性相位因果FIR滤波器。
这种方法的重点是选择一个合适的窗函数和理想滤波器.
设x(n)是一个长序列,w(n)是长度为N的窗函数,用w(n)截断x(n),得到N点序列xn(n),即
xn(n)=x(n)w(n)(3-2)
在频域上则有
(3—3)
由此可见,窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形,而且也会影响到频域内的形状。
MATLAB信号工具箱主要提供了以下几种窗函数,如下表所示:
表1MATLAB窗函数
窗
窗函数
矩形窗
Boxcar
三角窗
Barlett
布莱克曼窗
Blackman
哈明窗
Hamming
汉宁窗
Hanning
凯塞窗
Kaiser
切比雪夫窗
Chebwin
加矩形窗后的频谱和理想频谱可得到以下结论:
加窗使过渡带变宽,过渡带的带宽取决于窗谱的主瓣宽度。
矩形窗情况下的过渡带宽是
。
N越大,过渡带越窄、越陡;
过渡带两旁产生肩峰,肩峰的两侧形成起伏振荡。
肩峰幅度取决于窗谱主瓣和旁瓣面积之比.矩形窗情况下是8.95%,与N无关。
工程上习惯用相对衰耗来描述滤波器,相对衰耗定义为
(3—4)
这样两个肩峰点的相对衰耗分别是0。
74dB和-21dB。
其中(-0。
0895)对应的点的值定义为阻带最小衰耗.
以上的分析可见,滤波器的各种重要指标都是由窗函数决定,因此改进滤波器的关键在于改进窗函数。
窗函数谱的两个最重要的指标是:
主瓣宽度和旁瓣峰值衰耗。
旁瓣峰值衰耗定义为:
旁瓣峰值衰耗=20lg(第一旁瓣峰值/主瓣峰值)(3-5)
为了改善滤波器的性能,需使窗函数谱满足:
主瓣尽可能窄,以使设计出来的滤波器有较陡的过渡带;
第一副瓣面积相对主瓣面积尽可能小,即能量尽可能集中在主瓣,外泄少,使设计出来的滤波器的肩峰和余振小。
但上面两个条件是相互矛盾的,实际应用中,折衷处理,兼顾各项指标。
3.1线性相位
一个单一频率的正弦信号通过一个系统,假设它通过这个系统的时间需要t,则这个信号的输出相位落后原来信号wt的相位。
从这边可以看出,一个正弦信号通过一个系统落后的相位等于它的w*t;反过来说,如果一个频率为w的正弦信号通过系统后,它的相位落后delta,则该信号被延迟了delta/w的时间.在实际系统中,一个输入信号可以分解为多个正弦信号的叠加,为了使得输出信号不会产生相位失真,必须要求它所包含的这些正弦信号通过系统的时间是一样的.因此每一个正弦信号的相位分别落后,w1*t,w2*t,w3*t。
因此,落后的相位正比于频率w,如果超前,超前相位的大小也是正比于频率w.从系统的频率响应来看,就是要求它的相频特性是一条直线。
在FIR滤波器的设计中,为了得到线性相位的性质,通常利用实偶对称序列的相频特性为常数0和实奇对称序列为相频特性为常数90度的特点。
因此得到的是对称序列,不是因果序列,是不可实现系统,为了称为物理可实现系统,需要将它向右移动半个周期,这就造成了相移特性随时间的变化,同时也是线性变化。
单位脉冲响应h(n)(为实数)具有偶对称或奇对称性,则FIR数字滤波器具有严格的线性相位特性.
数字滤波器中,IIR数字滤波器方便简单,但它相位的非线性,要求采用全通网络进行相位校正,且稳定性难以保障。
FIR滤波器具有很好的线性相位特性,使得它越来越受到广泛的重视.
3。
2基本窗函数
数字信号处理领域中所用到的基本窗函数主要有:
矩形窗函数、三角窗函数和汉宁窗函数,哈明窗函数布莱克窗函数,凯塞窗函数等。
下面就对这些窗函数展开介绍。
3.2.1矩形窗函数
矩形窗(RectangularWindow)函数的时域形式可以表示为:
(3—2-1)
它的频域特性为
(3—2-2)
Boxcar函数:
生成矩形窗
调用方式w=boxcar(N):
输入参数N是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
从功能上讲,该函数又等价于w=ones(n,1)。
3.2.2三角窗函数
三角窗是最简单的频谱函数
为非负的一种窗函数。
三角窗函数的时域形式可以表示为:
当n为奇数时
(3—2—3)
当n为偶数时
(3—2-4)
Bartlett函数:
生成巴特利特窗
调用方式w=bartlett(n):
(1)输入参数n是窗函数的长度;
(2)输出参数w是由窗函数的值组成的n阶向量。
(3)三角窗也是两个矩形窗的卷积。
三角窗函数的首尾两个数值通常是不为零的.当n是偶数时,三角窗的傅立叶变换总是非负数。
3。
2。
3汉宁窗函数
汉宁窗函数的时域形式可以表示为:
(3-2—5)
它的频域特性为:
(3—2-6)
其中,
为矩形窗函数的幅度频率特性函数.
汉宁窗函数的最大旁瓣值比主瓣值低31dB,但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍,为8π/N。
hanning函数:
生成汉宁窗
调用方式
(1)w=hanning(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
注意:
此函数不返回是零点的窗函数的首尾两个元素。
(2)w=hanning(n,’symmetric'):
与上面相类似。
(3)w=hanning(n,’periodic'):
此函数返回包括为零点的窗函数的首尾两个元素。
3。
2.4哈明窗函数
海明窗函数的时域形式可以表示为
(3-2—7)
它的频域特性为
(3-2-8)
其中,
为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB,但它和汉宁窗函数的主瓣宽度是一样大的.
Hamming函数:
生成海明窗
调用方式
(1)w=hamming(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2)w=hamming(n,sflag):
参数sflag用来控制窗函数首尾的两个元素值;其取值为symmetric或periodic;默认值为symmetric。
3。
2。
5布莱克曼窗函数
布莱克曼窗函数的时域形式可以表示为
(3-2-9)
它的频域特性为
(3—2-10)
其中,
为矩形窗函数的幅度频率特性函数。
布莱克曼窗函数的最大旁瓣值比主瓣值低57dB,但是主瓣宽度是矩形窗函数的主瓣宽度的3倍,为12π/N.
Blackman函数:
生成海明窗
调用方式
(1)w=blackman(n):
输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2)w=blackman(n,sflag):
参数sflag用来控制窗函数首尾的两个元素值;其取值为symmetric或periodic;默认值为symmetric。
3.2。
6凯塞窗
上面所讨论的几种窗函数,在获得旁瓣抑制的同时却增加了主瓣的宽度。
而凯塞窗定义了一组可调的窗函数,它是由零阶贝塞尔函数构成的,其主瓣能量和旁瓣能量的比例是近乎最大的。
而且,这种窗函数可以在主瓣宽度和旁瓣高度之间自由选择它们的比重,使用户的设计变得非常灵活。
凯塞窗函数的时域形式可表示为
(3-2-11)
其中,
是第1类变形零阶贝塞尔函数,
是窗函数的形状参数,由下式确定:
(3-2—12)
其中,
为凯塞窗函数的主瓣值和旁瓣值之间的差值(dB)。
改变β的取值,可以对主瓣宽度和旁瓣衰减进行自由选择.β的值越大,窗函数频谱的旁瓣值就越小,而其主瓣宽度就越宽。
Kaiser函数:
生成凯塞窗
调用方式w=kaiser(n,beta):
输入参数n是窗函数的长度;输入参数beta用于控制旁瓣的高度;输出参数w是由窗函数的值组成的n阶向量。
n一定时,beta越大,其频谱的旁瓣就越小,但主瓣宽度相应的增加;当beta一定时,n发生变化,其旁瓣高度不会发生变化。
4方案设计与分析
用窗函数法设计一个FIR带阻滤波器。
指示如下:
下通带截至频率
;
上通带截止频率
;
阻带下限频率
;
阻带上限频率
通带最大衰减
阻带最小衰减
=
4。
1方案设计程序
6中窗函数的基本参数如下
窗函数类型
旁瓣峰值
/dB
过渡带宽度
阻带最小衰减
/dB
近似值
精确值
矩形窗
-13
4
/N
1。
8
/N
—12
三角窗
—25
8
/N
6。
1
/N
-25
汉宁窗
-31
8
/N
6。
2
/N
—44
哈明窗
—41
8
/N
6。
6
/N
—53
布莱克曼窗
-57
12
/N
11
/N
—74
凯塞窗
-57
10
/N
-80
因为阻带最小衰减
=
,所以选择布莱克曼窗或凯塞窗都可以设计,先以布莱克窗进行设计。
程序步骤如下:
wlp=0.2*pi;
wls=0。
35*pi;
wus=0。
65*pi;
wup=0。
8*pi;
wc=[(wlp+wls)/2/pi,(wus+wup)/2/pi];
B=wls—wlp;
N=ceil(12*pi/B)—1;
n=0:
N-1;
window=kaiser(N);
[h1,w]=freqz(window,1)
subplot(2,2,1)
stem(window,’.');
xlabel('n’);
title(’kaiser窗函数’);
subplot(2,2,2)
plot(w/pi,20*log(abs(h1)/abs(h1
(1))));
grid;
xlabel(’w/pi');
ylabel(’幅度(dB)’);
title(’kaiser窗函数的频谱');
hn=fir1(N—1,wc,'stop');
[h2,w]=freqz(hn,1,512);
subplot(2,2,3)
stem(n,hn,’。
’);
xlabel('n’);
ylabel('h(n)');
title(’kaiser窗函数的单位脉冲响应');
subplot(2,2,4)
plot(w/pi,20*log(abs(h2)/abs(h2
(1))));
grid;
xlabel('w/pi’);
ylabel(’幅度(dB)’);
title(’kaiser带阻滤波器的幅度特性’);
运行图形如下:
4.2分析
由kaiser窗函数程序可知,N—1=80,所以N=81,通过图形“kaiser带阻滤波器的幅度特性”可知通带内最大衰减约
为0,小于1。
满足指标要求。
当用凯赛窗时,只需在程序中调用函数blackman(N)即可。
其图形如下:
由两图可得,两种设计都符合要求,但要综合过渡带宽度,选择凯塞窗其宽度最小.
5总结与体会
设计带通滤波器时首先要计算出过渡带,然后查表得到不同窗函数所需要的阶数,不同的窗函数所设计的滤波器的形状各有差异,尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比较明显得差别,实际应用中应根据实际情况,折衷处理,兼顾各项指标。
为了这次课程设计,自己自学了数字信号处理领域中窗函数的有关知识。
实际中遇到的离散时间信号总是有限长的,因此不可避免地要遇到数据截断问题。
而在信号处理中,对离散序列的数据截断是通过序列与窗函数相乘来实现的.而且,有关滤波器的设计、功率谱估计等基本概念也要用到窗函数。
本次课程设计对经常用到的下面6窗函数:
矩形窗函数、三角窗函数、汉宁窗函数、哈明窗函函数、布莱克曼窗函数、凯塞窗函窗,先是做了基本概念上的阐释,然后对其MATLAB实现函数做出了说明,最后又结合具体的实例,对这些窗函数的频域特性等进行了介绍。
通过这次学习,我不但掌握了FIR数字滤波器窗函数的基本知识及其实际应用的技巧了,还提高了自己的编程和写报告的能力,收获颇多。
6参考文献
[1]《数字信号处理》(第三版),丁玉美,高西全.西安电子科技大学出版社,2000.
[2]《MATLAB及在电子信息课程中的应用》,陈怀堔,吴大正,高西全.电子工业出版社,2006。
[3]《MATLAB7.0从入门到精通》,求是科技。
人民邮电出版社,2006。
[4]《数字信号处理(第三版)》学习指导,高西全,丁玉美.西安科技大学出版社,2001.