FIR带阻滤波器的设计

课程设计任务书学生姓名：_ 专业班级： 指导教师： 工作单位：信息工程学院题 目： FIR 带阻滤波器的设计初始条件：具备数字信号处理的理论知识；具备 Matlab 编程能力； 熟悉带阻滤波器的设计原理； 提供编程所需要的计算机一台 要求完成的主要任务：（包括课程设计工作量及其技术要求，以及说明书撰写等 具体要求）1、 设计中心频率为200Hz，带宽为150Hz的FIR数字带阻 滤波器；2、 独立编写程序实现3、完成符合学校要求的设计说明书时间安排：一周，其中 3 天程序设计， 2 天程序调试 指导教师签名： 年 月 日系主任（或责任教师）签名：年月日目录摘要 I1 MATLAB 简介 12 设计原理 22.1 线性相位 FIR 滤波器的特点 22.2 窗函数法的设计原理 22.3 基本窗函数 52.3.1 矩形窗 52.3.2 巴特利特窗 52.3.3 汉宁窗 62.3.4 海明窗 72.3.5 布拉克曼窗 72.3.6 凯泽窗 83 程序设计 94 运行结果及分析 114.1 程序运行结果 114.2 带阻滤波器比较 155 心得体会 16参考文献 16摘要现代图像、语声、数据通信对线性相位的要求是普遍的。

正是此原因，使得 具有线性相位的fir数字滤波器得到大力发展和广泛应用本次课设中要求做的 fir 带阻滤波器也是一种应用极广的滤波器该滤波器 能让一些频率的波段通过，而抑制一些波段的波形在实际的设计中，使用 fir 型的滤波器具有线性相位，因而可以满足一些有特定相位要求的滤波器设计中， 可先用 matlab 软件进行仿真，编程设计中会用到窗函数，达到要求后再通过软 件编程或者硬件来实现关键词：线性相位；fir数字滤波器；窗函数；matlab1 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)是有限长的(OSnSN-1),其z变换为H (z)=艺 h (n) z - n (2.1)显然，H⑵是z-1的(N-1)阶多项式，在有限z平面(OVIzlV^)有(N-1)个零 点，而有(N-1)阶极点全部位于z平面的原点z=0处h(n丿的频率响应H(ej®)可表示为H (ej3)=艺 h( n)e - j^n (2.2)n—0 当h(n)为实序列时，可将H(ej®)表示成一H (eg) = ± H (eg) |ej e(u) — H (o)ej0(u) (2.3 )其中|H(ej3)|是真正的幅度响应，而实函数H@)称为幅度函数，0(®)称为相位函 数有两类准确的线性相位，要求分别满足0(3) — _T3 (2.4)0(3) = B-T3 (2.5)其中T，卩均为常数，表示相位是通过坐标原点®=0或是通过0(0)=卩的斜直线, 二者的群延时都是常数T = -d 0(3)。

一般将满足式(2.4)的相位称为第一类线d3性相位，满足式(2.5)的相位称为第二类线性相位2.2 窗函数法的设计原理数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间 序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形 式变化数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数 字滤波器分为两种，即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR) 滤波器IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配所以IIR 滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法FIR数字 滤波器的单位脉冲响应是有限长序列它的设计问题实质上是确定能满足所要求 的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等因此设计 FIR 滤波器的方法之一可以从时域出发，截取有限长的一段冲击响 应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数只要N足够长， 截取的方法合理，总能满足频域的要求一般这种时域设计、频域检验的方法要 反复几个回合才能成功要设计一个线性相位的 FIR 数字滤波器，首先要求理想 频率响应Hd(ejw)。

Hd(ejw)是w的周期函数，周期为2沢，可以展开成傅氏级数H (ejw) = £ h (n)e-jwn (2.6)dd其中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)在频域上则有je n X (jj. W2n n—n2.8)由此可见，窗函数w(n)不仅仅会影响原信号x(n)在时域上的波形，而且也会 影响到频域内的形状MATLAB信号工具箱主要提供了以下几种窗函数，如表2.1所示 加矩形窗后的频谱和理想频谱可得到以下结论：加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度矩形窗情况下的 过渡带宽是滋/N。

N越大，过渡带越窄、越陡；过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡肩峰幅度取决于窗谱主瓣 和旁瓣面积之比矩形窗情况下是 8.95％，与 N 无关工程上习惯用相对衰耗 来描述滤波器，相对衰耗定义为A( w) = 20lg[| H (ejw )/| H (ej o)|] = 201g[| H (w)/ H (0)|] (2.9) 这样两个肩峰点的相对衰耗分别是0.74dB和-21dB其中（-0.0895）对应的点的 值定义为阻带最小衰耗以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波 器的关键在于改进窗函数表 2.1 MATLAB 窗函数窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗旁瓣峰值衰耗定义为旁瓣峰值衰耗= 201g （第一旁瓣峰值/主瓣峰值） （2.10）为了改善滤波器的性能，需使窗函数谱满足：① 主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带；② 第一副瓣面积相对主瓣面积尽可能小，即能量尽可能集中在主瓣，外泄少，使设计出来的滤波器的肩峰和余振小 但上面两个条件是相互矛盾的，实际应用中，折衷处理，兼顾各项指标 设计滤波器尽量要求窗函数满足以下两项要求：（1） 窗谱主瓣尽可能地窄，以获取较陡的过渡带。

2） 尽量减少窗谱的最大旁瓣的相对幅度也就是能量尽量集中于主瓣，这 样使尖峰和波纹减小，就可增大阻带的衰减但是这两项要求是不能同时满足的当选用主瓣宽度较窄时，虽然得到陡峭 的过渡带，但通带和阻带的波动明显增加；当选用最小的旁瓣幅度时，虽能得到 平坦的幅度响应和较小的阻带波纹，但过渡带加宽，即主瓣会加宽因此，实际 所选用的窗函数往往是它们的折中表 2.2 6 种窗函数的基本参数窗函数类型旁瓣峰值a n/dB过渡带宽度阻带最小衰减a /dBs近似值精确值矩形窗-134兀/N1.8 兀 /N-12巴特利特窗-258兀/N6.1 兀 /N-25汉宁窗-318兀/N6.2 兀 /N-44海明窗-418兀/N6.6 兀 /N-53布拉克曼窗-5712 兀 /N11 兀 /N-74凯泽窗-5710 兀 /N-802.3 基本窗函数数字信号处理领域中所用到的基本窗函数主要有矩形窗函数、三角窗函数、 汉宁窗函数、哈明窗函数、布莱克窗函数和凯塞窗函数下面就对这些窗函数展 开介绍2.3.1 矩形窗矩形窗(Rectangular Window)w函数的时域形式可以表示为2.11)[1, 0 < n < N-1w(n)=W)彳0,其他它的频域特性为2.12)•(sm —12丿Boxcar 函数：生成矩形窗调用方式w = boxcar (N)：输入参数N是窗函数的长度，输出参数w是由窗函数的值组成的 n 阶向量。

从功能上讲，该函数等价于 w = ones(n,1)2.3.2 巴特利特窗巴特利特窗(三角窗)是最简单的频谱函数W(ej®为非负的一种窗函数巴 特利特窗函数的时域形式可以表示为当n为奇数时w(k)= <2(k -1)n -1 ，2 - 2(k-1)n -12.13)

注意：此函数不返回是零点的窗函数的首尾两个元素2) w = hanning(n,'symmetric'):与上面相类似3) w = hanning(n,'periodic') ：此函数返回包括为零点的窗函数的首尾两个元 素2.3.4 海明窗2.17)海明窗函数的时域形式可以表示为( k \w(k) = 0.54一 0.46cos 2n k = 1,2,…,NI N -1 丿它的频域特性为W(①)=0.54W (®) + 0.23RN 一1 丿 R2.18)其中，WRg)为矩形窗函数的幅度频率特性函数R海明窗函数的最大旁瓣值比主瓣值低41dB，但它和汉宁窗函数的主瓣宽度 是一样大的Hamming 函数:生成海明窗调用方式:(1) w = hamming(n):输入参数n是窗函数的长度；输出参数w是由窗函数 的值组成的 n 阶向量2) w = hamming(n,sflag):参数sflag用来控制窗函数首尾的两个元素值；其 取值为 symmetric 或 periodic； 默认值为 symmetric2.3.5 布拉克曼窗布拉克曼窗函数的时域形式可以表示为(k — 1、( k 一1 、2n+ 0.08cos4 n I N -1 丿1 N -1 丿w(k) = 0.42 一 0.5cosk = 1,2,…,N (2」9)它的频域特性为W(w) = 0.42 W 6)+ 0.25R+0.04(W (w 一RI4n、N 一 1丿(+W (w +RI4n讣N -1 丿］2.20)其中，Wr (w)为矩形窗函数的幅度频率特性函数。

R布拉克曼窗函数的最大旁瓣值比主瓣值低57dB，但是主瓣宽度是矩形窗函 数的主瓣宽度的3倍，为12n/NBlackman 函数：生成布拉克曼窗调用方式：(1) w = blackman (n)：输入参数n是窗函数的长度；输出参数w是由窗函数 的值组成的 n 阶向量2) w = blackman (n,sflag) ：参数 sflag 用来控制窗函数首尾的两个元素值； 其取值为 symmetric 或 periodic；默认值为 symmetric2.3.6 凯泽窗上面所讨论的几种窗函数，在获得旁瓣抑制的同时却增加了主瓣的宽度而 凯泽窗定义了一组可调的窗函数，它是由零阶贝塞尔函数构成的，其主瓣能量和 旁瓣能量的比例是近乎最大的而且，这种窗函数可以在主瓣宽度和旁瓣高度之 间自由选择它们的比重，使用户的设计变得非常灵活I PJ1— 1 —0凯泽窗函数的时域形式可表示为0< k < N-1 (2.2DI (P)0其中，10(P)是第1类变形零阶贝塞尔函数，P是窗函数的形状参数，由式2.22)确定.2.22)'0.1102(a- 8.7), a > 50P = \ 0.5482(a — 21)0.4 + 0.07886(a — 21), 21

改变卩的取值， 可以对主瓣宽度和旁瓣衰减进行自由选择P的值越大，窗函数频谱的旁瓣值就 越小，而其主瓣宽度就越宽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程序运行结果0.51020J0£P 砸 InE -E-%Boxcar®函数的单位脉冲响应Boxcar囱函刁-50-100-150-200Boxcar囱函数的频谱0.6 r10Q W 300 ITHz即带阻滤波器的幅频特性100OOOOOZf/H00002AU3~图4.1矩形窗函数经过上述仿真，可以看出矩形窗基本可以满足设计要求，将程序中的相关语 句改变后，即可得到其他窗函数的频谱图，得到的结果如下。