课程设计任务书
学生姓名: 专业班级:
指导教师: 工作单位: 信息工程学院
题 目:FIR带阻滤波器的设计 初始条件:
具备数字信号处理的理论知识;
具备Matlab编程能力; 熟悉带阻滤波器的设计原理; 提供编程所需要的计算机一台
要求完成的主要任务:(包括课程设计工作量及其技术要求,以及说明书撰写等
具体要求)
1、设计中心频率为200Hz,带宽为150Hz的FIR数字带阻滤波器;
2、独立编写程序实现
3、完成符合学校要求的设计说明书 时间安排:
一周,其中3天程序设计,2天程序调试
指导教师签名: 年 月 日
系主任(或责任教师)签名: 年 月 日
目 录
摘要 ............................................................................................................................... I 1 MATLAB简介 ......................................................................................................... 1 2 设计原理 .................................................................................................................. 2
2.1 线性相位FIR滤波器的特点 ....................................................................... 2 2.2 窗函数法的设计原理 .................................................................................... 2 2.3 基本窗函数 .................................................................................................... 5
2.3.1 矩形窗 .................................................................................................. 5 2.3.2 巴特利特窗 .......................................................................................... 5 2.3.3 汉宁窗 .................................................................................................. 6 2.3.4 海明窗 .................................................................................................. 7 2.3.5 布拉克曼窗 .......................................................................................... 7 2.3.6 凯泽窗 .................................................................................................. 8
3 程序设计 .................................................................................................................. 9 4 运行结果及分析 .................................................................................................... 11
4.1 程序运行结果 .............................................................................................. 11 4.2 带阻滤波器比较 .......................................................................................... 15 5 心得体会 ................................................................................................................ 16 参考文献 .................................................................................................................... 16
摘 要
现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因,使得具有线性相位的fir数字滤波器得到大力发展和广泛应用。
本次课设中要求做的fir带阻滤波器也是一种应用极广的滤波器。该滤波器能让一些频率的波段通过,而抑制一些波段的波形。在实际的设计中,使用fir型的滤波器具有线性相位,因而可以满足一些有特定相位要求的滤波器。设计中,可先用matlab软件进行仿真,编程设计中会用到窗函数,达到要求后再通过软件编程或者硬件来实现。
关键词:线性相位;fir数字滤波器;窗函数;matlab。
1 MATLAB简介
MATLAB是“矩阵实验室”(Matrix Laboratory)的缩写,是一种科学计算软件,主要适用于矩阵运算及控制和信息处理领域的分析设计。它使用方便,输入简捷,运算高效,内容丰富,因此很多专家在自己擅长的领域用它编写了许多专门的MATLAB工具包。由于MATLAB功能的不断扩展,所以是科学研究中最常用必不可少的工具。
MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件,其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级,MATLAB的用户界面也越来越精致,更加接近Windows的标准界面,人机交互性更强,操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统,极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统,程序不必经过编译就可以直接运行,而且能够及时地报告出现的错误及进行出错原因分析。
MATLAB是一种高级的矩阵/阵列语言,它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步,也可以先编写好一个较大的复杂的应用程序(M文件)后再一起运行。新版本的MATLAB语言是基于最为流行的C++语言基础上的,因此语法特征与C++语言极为相似,而且更加简单,更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强,这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。 MATLAB 产品族可以用来进行数值分析、数值和符号计算、工程与科学绘图、控制系统的设计与仿真、数字图像处理技术、数字信号处理技术、通讯系统设计与仿真、财务与金融工程等各种工作。
MATLAB 的应用范围非常广,包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱(单独提供的专用MATLAB 函数集)扩展了MATLAB 环境,以解决这些应用领域内特定类型的问题。
2 设计原理
2.1 线性相位FIR滤波器的特点
FIR系统最大的特点之一就是能够做成严格的线性相位。
FIR滤波器的单位冲激响应h(n)是有限长的(0≤n≤N-1),其z变换为
N-1
H(z)=
-1
∑h(n)z (2.1)
-n
n=0
显然,H(z)是z的(N-1)阶多项式,在有限z平面(0<|z|<∞)有(N-1)个零点,而有(N-1)阶极点全部位于z平面的原点z=0处。 h(n)的频率响应H(ejω)可表示为
N-1
H(e
jω
jω
jω
)=
∑h(n)e
n=0jω
-jωn
(2.2)
当h(n)为实序列时,可将H(e)表示成
H(e
)=±H(e
)e
jθ(ω)
=H(ω)e
jθ(ω)
(2.3)
其中H(ejω)是真正的幅度响应,而实函数H(ω)称为幅度函数,θ(ω)称为相位函数。有两类准确的线性相位,要求分别满足
θ(ω)=-τω
(2.4) (2.5)
θ(ω)=β-τω
其中τ,β均为常数,表示相位是通过坐标原点ω=0或是通过θ(0)=β的斜直线, 二者的群延时都是常数τ=-dθ(ω)dω
。一般将满足式(2.4)的相位称为第一类线
性相位,满足式(2.5)的相位称为第二类线性相位。
2.2 窗函数法的设计原理
数字滤波器可以理解为是一个计算程序或算法,将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列,并在转化过程中,使信号按预定的形式变化。数字滤波器有多种分类,根据数字滤波器冲激响应的时域特征,可将数字滤波器分为两种,即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。IIR数字滤波器具有无限宽的冲激响应,与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题,设计方法主要有窗函数法、频率采样法和等
波纹最佳逼近法等。
因此设计FIR滤波器的方法之一可以从时域出发,截取有限长的一段冲击响应作为H(z)的系数,冲击响应长度N就是系统函数H(z)的阶数。只要N足够长,截取的方法合理,总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器,首先要求理想频率响应
Hd(e
jw
)
。
Hd(e
jw
)
是w的周期函数,周期为2π,可以展开成傅氏级数
∞
Hd(e
jw
)=
∑h
n=-∞
d
(n)e
-jwn
(2.6)
其中
hd(n)
是与理想频响对应的理想单位抽样响应序列。但不能用来作为设
计FIR-DF用的h(n),因为hd(n)一般都是无限长、非因果的,物理上无法实现。为了设计出频响类似于理想频响的滤波器,可以考虑用h(n)来近似
hd(n)
。
窗函数的基本思想:先选取一个理想滤波器(它的单位抽样响应是非因果、无限长的),再截取(或加窗)它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。
设x(n)是一个长序列,w(n)是长度为N的窗函数,用w(n)截断x(n),得到N点序列an(n),即
an(n) = x(n)w(n) (2.7) 在频域上则有
X
N
(e)=
jω
⎰2π
1
π-π
Xe
()⋅W(e(
jθ
jω-θ
)
)dθ (2.8)
由此可见,窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形,而且也会影响到频域内的形状。
MATLAB信号工具箱主要提供了以下几种窗函数,如表2.1所示。 加矩形窗后的频谱和理想频谱可得到以下结论:
加窗使过渡带变宽,过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是4π/N。N越大,过渡带越窄、越陡;
过渡带两旁产生肩峰,肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95%,与N无关。工程上习惯用相对衰耗来描述滤波器,相对衰耗定义为
A(w)=20lg[H(ejw)/H(ej0)]=20lgH(w)/H(0)] (2.9)
这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中(-0.0895)对应的点的值定义为阻带最小衰耗。
以上的分析可见,滤波器的各种重要指标都是由窗函数决定,因此改进滤波器的关键在于改进窗函数。
表2.1 MATLAB窗函数
窗函数谱的两个最重要的指标是:主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为
旁瓣峰值衰耗=20lg(第一旁瓣峰值/主瓣峰值) (2.10) 为了改善滤波器的性能,需使窗函数谱满足:
①主瓣尽可能窄,以使设计出来的滤波器有较陡的过渡带;
②第一副瓣面积相对主瓣面积尽可能小,即能量尽可能集中在主瓣,外泄少,使设计出来的滤波器的肩峰和余振小。
但上面两个条件是相互矛盾的,实际应用中,折衷处理,兼顾各项指标。 设计滤波器尽量要求窗函数满足以下两项要求: (1)窗谱主瓣尽可能地窄,以获取较陡的过渡带。
(2)尽量减少窗谱的最大旁瓣的相对幅度。也就是能量尽量集中于主瓣,这样使尖峰和波纹减小,就可增大阻带的衰减。
但是这两项要求是不能同时满足的。当选用主瓣宽度较窄时,虽然得到陡峭的过渡带,但通带和阻带的波动明显增加;当选用最小的旁瓣幅度时,虽能得到平坦的幅度响应和较小的阻带波纹,但过渡带加宽,即主瓣会加宽。因此,实际所选用的窗函数往往是它们的折中。
表2.2 6种窗函数的基本参数
2.3 基本窗函数
数字信号处理领域中所用到的基本窗函数主要有矩形窗函数、三角窗函数、汉宁窗函数、哈明窗函数、布莱克窗函数和凯塞窗函数。下面就对这些窗函数展开介绍。
2.3.1 矩形窗
矩形窗(Rectangular Window)w函数的时域形式可以表示为
它的频域特性为
WRe
⎧1,
w(n)=RN(n)=⎨
⎩0,
0≤n≤N-1
其他
(2.11)
()=e
jω
⎛N-1⎫-j ⎪ω⎝2⎭
⎛ωN⎫sin ⎪⎝2⎭⎛ω⎫sin ⎪⎝2⎭
(2.12)
Boxcar函数:生成矩形窗
调用方式w = boxcar (N):输入参数N是窗函数的长度,输出参数w是由窗函数的值组成的n阶向量。从功能上讲,该函数等价于w = ones(n,1)。
2.3.2 巴特利特窗
巴特利特窗(三角窗)是最简单的频谱函数W(e特利特窗函数的时域形式可以表示为
jω
)
为非负的一种窗函数。巴
当n为奇数时
n+1⎧2(k-1)
,1≤k≤⎪2w(k)=⎨n-1
2(k-1)n+1⎪2-,≤k≤n
n-12⎩
(2.13)
当n为偶数时
n⎧2(k-1)
,1≤k≤⎪2w(k)=⎨n-1
2(n-k)n⎪,≤k≤n
2⎩n-1
(2.14)
Bartlett函数:生成巴特利特窗 调用方式w = bartlett(n): (1) 输入参数n是窗函数的长度;
(2) 输出参数w是由窗函数的值组成的n阶向量。 (3) 巴特利特窗也是两个矩形窗的卷积。
巴特利特窗函数的首尾两个数值通常是不为零的。当n是偶数时,巴特利特的傅立叶变换总是非负数。
2.3.3 汉宁窗
汉宁窗(升余弦窗)函数的时域形式可以表示为
⎛k⎫⎫⎛ w(k)=0.5 2π⎪⎪ 1-cos⎪
n+1⎭⎭⎝⎝
k=1,2, ,N
(2.15)
它的频域特性为
⎧⎡2π⎫2π⎫⎤⎫-jω ⎛⎛⎝
W(ω)=⎨0.5WR(ω)+0.25⎢WR ω-⎪+WR ω+⎪⎥⎬e
N-1⎭N-1⎭⎦⎭⎝⎝⎣⎩
⎛N-1⎫
⎪2⎭
(2.16)
其中,WR(ω)为矩形窗函数的幅度频率特性函数。
汉宁窗函数的最大旁瓣值比主瓣值低31dB,但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍,为8π/N。
Hanning函数:生成汉宁窗 调用方式:
(1) w = hanning(n):输入参数n是窗函数的长度,输出参数w是由窗函数的值组成的n阶向量。
注意:此函数不返回是零点的窗函数的首尾两个元素。
(2) w = hanning(n,'symmetric'):与上面相类似。
(3) w = hanning(n,'periodic'):此函数返回包括为零点的窗函数的首尾两个元素。
2.3.4 海明窗
海明窗函数的时域形式可以表示为
k⎫⎛
w(k)=0.54-0.46cos 2π⎪k=1,2, ,N
N-1⎝⎭
(2.17)
它的频域特性为 其中,W
R
⎡2π⎫2π⎫⎤⎛⎛
W(ω)=0.54WR(ω)+0.23⎢WR ω-⎪+WR ω+⎪⎥
N-1N-1⎝⎭⎝⎭⎦⎣
(2.18)
(ω)
为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB,但它和汉宁窗函数的主瓣宽度是一样大的。
Hamming函数:生成海明窗 调用方式:
(1) w = hamming(n):输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2) w = hamming(n,sflag):参数sflag用来控制窗函数首尾的两个元素值;其取值为symmetric或periodic;默认值为symmetric。
2.3.5 布拉克曼窗
布拉克曼窗函数的时域形式可以表示为
k-1⎫k-1⎫⎛⎛
w(k)=0.42-0.5cos 2π⎪+0.08cos 4π⎪
N-1⎭N-1⎭⎝⎝
k=1,2, ,N
(2.19)
它的频域特性为
⎛
W(ω)=0.42WR(ω)+0.25WR ω-⎢
⎣
⎝⎡
2π⎫⎤⎫⎛
⎪+WR ω+⎪⎥N-1⎭N-1⎭⎦⎝2π
⎡4π⎫4π⎫⎤⎛⎛
+0.04⎢WR ω-⎪+WR ω+⎪⎥
N-1⎭N-1⎭⎦⎝⎝⎣
(2.20)
其中,W
R
(ω)
为矩形窗函数的幅度频率特性函数。
布拉克曼窗函数的最大旁瓣值比主瓣值低57dB,但是主瓣宽度是矩形窗函数的主瓣宽度的3倍,为12π/N。
Blackman函数:生成布拉克曼窗 调用方式:
(1) w = blackman (n):输入参数n是窗函数的长度;输出参数w是由窗函数的值组成的n阶向量。
(2) w = blackman (n,sflag):参数sflag用来控制窗函数首尾的两个元素值;其取值为symmetric或periodic;默认值为symmetric。
2.3.6 凯泽窗
上面所讨论的几种窗函数,在获得旁瓣抑制的同时却增加了主瓣的宽度。而凯泽窗定义了一组可调的窗函数,它是由零阶贝塞尔函数构成的,其主瓣能量和旁瓣能量的比例是近乎最大的。而且,这种窗函数可以在主瓣宽度和旁瓣高度之间自由选择它们的比重,使用户的设计变得非常灵活。
凯泽窗函数的时域形式可表示为
2⎤⎡
2k⎫⎛
I0⎢β- 1-⎪⎥
N-1⎢⎝⎭⎥
⎣⎦
w(k)=
I0(β)
0≤k≤N-1 (2.21)
其中,I0(β)是第1类变形零阶贝塞尔函数,β是窗函数的形状参数,由式(2.22)确定.
α>50⎧0.1102(α-8.7),
⎪0.4
β=⎨0.5482(α-21)+0.07886(α-21),21≤α≤50
⎪0,α
(2.22)
其中,α为凯泽窗函数的主瓣值和旁瓣值之间的差值(dB)。改变β的取值,可以对主瓣宽度和旁瓣衰减进行自由选择。β的值越大,窗函数频谱的旁瓣值就越小,而其主瓣宽度就越宽。
Kaiser函数:生成凯泽窗 调用方式:
w = kaiser(n,beta):输入参数n是窗函数的长度;输入参数beta用于控制旁瓣的高度;输出参数w是由窗函数的值组成的n阶向量。n一定时,beta越大,其频谱的旁瓣就越小,但主瓣宽度相应的增加;当beta一定时,n发生变化,其旁瓣高度不会发生变化。
3 程序设计
任务书中给出的要求为中心频率200Hz,带宽150Hz。故设上通带截止频率为110Hz,下通带截止频率290Hz,阻带上限频率140Hz,阻带下限频率260Hz。
此处仅以Boxcar窗为示例,其他窗函数的程序代码基本相同,只是在window=Boxcar(N)、N=ceil(1.8*pi/delta_w)两处作出各个窗函数相应的修改即可。 MATLAB程序:
flp=110; fhp=290; fls=140; fhs=260; fs=600; wlp=2*pi*flp/fs; whp=2*pi*fhp/fs; wls=2*pi*fls/fs; whs=2*pi*fhs/fs;
wc=[(wlp+wls)/(2*pi),(whp+whs)/(2*pi)]; delta1=wls-wlp; delta2=whp-whs;
delta_w=min(delta1,delta2);
N=ceil(1.8*pi/delta_w); //不同的窗要选择系数不同// N=N+rem(N,2); n=0:N-1;
window=Boxcar(N+1); //选择窗函数// [h1,w]=freqz(window,1); subplot(2,2,1) stem(window,'.'); xlabel('n');
title('Boxcar窗函数');
subplot(2,2,2)
plot(w*fs/(2*pi),20*log(abs(h1)/abs(h1(1)))); grid; xlabel('f/Hz'); ylabel('幅度(dB)');
title('Boxcar窗函数的频谱'); hn=fir1(N,wc,'stop',window); [h2,w]=freqz(hn,1,512); subplot(2,2,3) stem(hn,'.'); xlabel('n'); ylabel('h(n)');
title('Boxcar窗函数的单位脉冲响应'); subplot(2,2,4)
plot(w*fs/(2*pi),20*log(abs(h2)/abs(h2(1)))); grid; xlabel('f/Hz'); ylabel('幅度(dB)');
title('Boxcar带阻滤波器的幅频特性');
4 运行结果及分析
4.1 程序运行结果
图4.1 矩形窗函数
经过上述仿真,可以看出矩形窗基本可以满足设计要求,将程序中的相关语句改变后,即可得到其他窗函数的频谱图,得到的结果如下。
图4.2 巴特利特窗
图4.3 汉宁窗
图4.4 海明窗
图4.5 布拉克曼窗
图4.6 凯泽窗
4.2 fir带阻滤波器比较
以上比较可知,窗函数不同,频谱也就不一样,综合比较,Blackman窗最合适。
5 心得体会
从考试结束到一月六号,前后一个星期的时间里,我把时间都投入到了两
个课设上面,时间有点仓促,课设也有一定的难度,但最后在老师和同学们的共同帮助下,终于完成了。回顾此次课设,自己有不小的收获。
首先,这次课设我们都要用到matlab这种仿真软件,虽然很早就接触到了这种软件,但只是会一点皮毛。这个课设要用到很多新的函数,运用起来有一定的困难,不过最后通过查询一些资料,能较好地掌握这些知识。主要的困难在后面的理论部分。
为了完成此次课设,我再次翻阅了所学的理论知识,对题目有了一定的理解后,开始相关的设计。设计带阻滤波器时首先要计算出过渡带,然后查表得到不同窗函数所需要的阶数,不同的窗函数所设计的滤波器的形状各有差异,尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比较明显得差别,实际应用中应根据实际情况,折衷处理,兼顾各项指标。
为了这次课程设计,自己自学了数字信号处理领域中窗函数的有关知识。实际中遇到的离散时间信号总是有限长的,因此不可避免地要遇到数据截断问题。而在信号处理中,对离散序列的数据截断是通过序列与窗函数相乘来实现的。而且,有关滤波器的设计、功率谱估计等基本概念也要用到窗函数。本次课程设计对经常用到的下面6种窗函数:矩形窗函数、巴特利特函数、汉宁窗函数、海明窗函数、布拉克曼窗函数、凯泽窗函窗,先是做了基本概念上的阐释,然后对其MATLAB实现函数做出了说明,最后又结合具体的实例,对这些窗函数的频域特性等进行了介绍。
通过此次课设,我对fir滤波器的有关知识有了更深入的认识,同时提高了读程序和写程序的能力,理论和实践都有了提高。课设做完后,也发现了自己的一些不足,平时很少自己动手写程序,以至于用的时候有很多困难,在以后的时间里,我会多看看他人写的程序,然后自己动手写一部分,提高自己的动手实践能力。
参考文献
[1]刘泉,厥大顺,郭志强编. 数字信号处理原理与实现. 北京:电子工业出版
社,2009.6.
[2]董长虹. MATLAB信号处理与应用. 北京:国防工业出版社,2005.4. [3]王济. MATLAB在振动信号处理中的应用. 北京:中国水利水电出版社、知识
产权出版社,2006.9.
[4]张志涌. 精通MATLAB6. 5版[M]. 北京:北京航空航天大学出版社,2004.1. [5]候正信译. 数字信号处理基础. 北京:电子工业出版社,2003.8.
本科生课程设计成绩评定表
指导教师签字:
年 月 日