数字信号处理 数字滤波器设计.doc

数字信号处理 带通滤波器的设计带通滤波器的设计1 FIR滤波器的原理与技术FIR 滤波器是数字滤波器的一种。数字滤波器是用于修正或改变时域或频域中信号的属性。常见的是数字滤波器是线性时间不变(LineTime-Invariant,LTI)滤波器。1.1 数字滤波器的定义及分类数字滤波器是指输入和输出均为数字信号，且通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的软件或器件。数字滤波器和快速傅里叶变换一样，是数字信号处理学科的重要组成部分，应用非常广泛。数字滤波器，通常是指一种算法，或一种数字处理设备，它的功能是将一组输入的数字序列经过一定的运算后变换为另一组输出的数字序列。因此，数字滤波器既可以是用数字硬件装配成的一台完成给定运算的专用的数字计算机，也可以是所需要的运算编写程序，让通用计算机来执行。对于数字滤波器而言，若系统函数为H（z），其脉冲响应为h（n），输入时间序列为x（n），则它们在时域内的关系是如下： y(n)=h(n)*x(n) 式（1-1）在Z域内，输入和输出存在关系如下： Y(z)=H(z)*X(z) 式（1-2）式中，X(z),Y(z)分别为x(n)和y(n)的Z变换。在频域内，输入和输出则存在如下关系： Y(jw)=H(jw)*X(jw) 式（1-3）式中，H(jw)是数字滤波器的频率特性，X(jw)、Y(jw)分别为x(n)和y(n)的频谱，而w为数字角频率。数字滤波器可以有很多分类方法，但总体上可分为两大类。一类为经典滤波器，即一般的滤波器，特点是输入信号中的有用成分和希望能滤除的成分占用不同的频带，通过合适的选频滤波器可以实现滤波。例如，若输入信号中有干扰，信号中信号和干扰信号的频带互不重叠，即可滤除干扰得到想要的信号。另一类为现代滤波器，如维纳滤波器、卡尔曼滤波器、自适应滤波器等最佳滤波器，其输入信号中有用信号和希望滤除的成分频带重叠，可按照随机信号内部的一些统计分布规律，从干扰中最佳地提取信号。对于一般滤波器，从功能上分类可以分为低通、高通、带通、带阻和全通滤波器；从实现的网络结构或单位冲激响应上分类，可以分为有限长单位冲激响应（FIR）数字滤波器和无限长单位冲激响应（IIR）数字滤波器。对于有限长冲击响应数字滤波器（FIR），其输出y(n)只取决于有限个过去和现在的输入，x(n),x(n-1),x(n-2),x(n-m)，滤波器的输入输出关系可表示为 y(n)= 式（1-4）对于无线长冲击响应数字滤波器（IIR），它的输出不仅取决于过去和现在的输入，而且还取决于过去的输出，其差分方程为 y(n)+ = 式（1-5）该差分方程的单位冲击响应是无限延续的。1.2线性相位FIR数字滤波器的条件和特点本节中主要介绍FIR滤波器具有的线性相位的条件及滤波其器的特点以及其网络结构和特性。1.2.1 线性相位条件对于长度为N的h(n)，传输函数为 式（1-6） 式（1-7） 式中，Hg()称为幅度特性，()称为相位特性。注意，这里Hg()不同于|H(ej)|,Hg()为的实函数，可能取负值，而|H(ej)|总是正值。H(ej)线性相位是指()是的线性函数，即 ()=, 为常数 如果()满足下式： ()=0-, 0是起始相位 严格地说，此时()不具有线性相位，但以上两种情况都满足群时延是一个常数，即式（1-8） 也称这种情况为线性相位。1.2.2 线性相位型 FIR 滤波器的特点按结构划分，FIR 滤波器有四种类型：直接型、级联型、线性相位型和频率采样型。当一个 FIR 滤波器具有线性相位响应时，其脉冲响应呈现某种对称条件。在这种形式中，将利用这些对称关系把相乘的次数减少一半。正是由于此原因，在本文中才会选择线性相位型来进行设计。线性相位响应的优点如下：(1) 设计问题中只有实数运算而没有复数运算，运算简单；(2) 线性相位 FIR 滤波器没有延时失真，仅有某一固定时延，便于设计；(3) 对于长度为 M(或 M-1)的滤波器，其运算次数具有 M/2 量级，提高运算速度，节省资源。1.2.3 线性相位 FIR 滤波器的结构FIR网络结构特点是没有反馈电路，即没有环路，其单位脉冲响应是有限长的。设单位脉冲响应h(n)长度为N。其系统函数H(Z)和差分方程分别为 式（1-8） 式（1-9）为了突出线性相位滤波器的优点，先来看一下滤波器的结构。(1) 直接型结构。其结构如图 1-1 所示图1-1直接型结构(2) 级联型结构。结构如图1-2所示将式（1-4）分解成二阶实系数因子形式 式（1-10）图1-2级联型结构(3) 频率采样结构频率域等间隔采样，相应的时域信号会以采样点数为周期进行周期性延拓，如果在频率域采样点数N大于等于原序列的长度M，则不会引起信号失真，此时原序列的Z变换H(z)与频率采样H(k)满足下面关系式 式（1-11）图1-3 FIR滤波器频率采样结构而线性相位 FIR 滤波器的差分方程为：y( n)= h0 x(n)+h1x(n-1)+.+h1x(n-M+2)+h0x(n-M+1)= h0 x(n)+x(n-M+1)+h1x(n-1)+x(n-M+2)+.当滤波器的阶数为奇数和偶数时，其相应的结构如图所示图1-3 第一类线性相位网络结构图1-4 第二类线性相位网络结构1.2.4 线性相位 FIR 滤波器的性质线性相位的 FIR 滤波器有两种对称的类型，对称和反对称。在对称和反对称的情形下，脉冲响应的特性如下表所述。表1-1 线性相位FIR滤波器的幅度特性与相位特性表1.3 窗函数法设计FIR滤波器FIR的设计方法与IIR滤波器的设计方法有所不同，不能采用有模拟滤波器的设计进行转换的方法，其设计方法有窗函数法和频率采样法进行设计，也可以借助计算机辅助设计软件采用切比雪夫等波纹逼近法进行设计。窗函数设计技术是FIR滤波器设计的主要方法之一，由于其运算方便，物理意义直观，已成为工程实践中应用最广泛的方法。窗函数设计FIR滤波器的基本思想，就是根据给定的滤波器技术指标，选择滤波器长度和窗函数，使其具有最窄的主瓣和最小的旁瓣。其核心是从给定的频率特性通过加窗以确定有限长单位脉冲响应序列h(n)。有限长单位脉冲响应滤波器的系统函数为： 式（1-12） 式（1-13） 式（1-14）相应的单位取样响应hd(n)为 式（1-15）为了构造一个长度为N的线性相位滤波器，只有将h-d(n)截取一段，并保证截取的一段对(N-1)/2对称。设截取的一段用h(n)表示，即h(n)=hd(n)RN(n) 式（1-16）我们实际实现的滤波器的单位取样响应为h(n),长度为N，其系统函数为H(z)。以上就是用窗函数法设计FIR滤波器的思路。另外，我们知道是一个以2为周期的函数，可以展为傅氏级数，即 式（1-17）对（1-8）式进行傅里叶变换，根据复卷积定理，得到：式中，和分别是和的傅里叶变换，即 式（1-18） 式（1-19）对称为矩形窗的幅度函数；将写成下式： 式（1-20）按照(1-6)式，理想低通滤波器的幅度特性为 式（1-21）将和代入(1-13)式，得到：将写成下式： 式（1-22） 设 式（1-23）则实际设计的FIR滤波器的幅频特性为 式（1-24）显然，对实际FIR滤波器的幅频特性有影响的只是窗函数的幅频特性。实际FIR滤波器的幅频特性是理想低通滤波器的幅频特性与窗函数的幅频特性的复卷积。复卷积可用下图说明。图1-5 矩形窗对理想低通幅频特性的影响从图中可以看到(1)=0时的响应H(0),式1-24是图中（a）和（b）两个函数乘积的积分。(2)=时的响应刚好与的一半重叠，卷积值刚好是H(0)的一半，如图1-5(C)所示。(3)时的响应，的全部主瓣都在的通带之内，如图1-5（d）所示。(4)时的响应，的全部主瓣都在的通带之外，如图1-5（e）所示。通过分析可知，对hd(n)加矩形窗处理后，H()和原理想低通Hd()差别有以下两点：(1)在理想特性不连续点=c附近形成过渡带。过渡带的宽度，近似等于RN()主瓣宽度，即4/N。(2)通带内增加了波动，最大的峰值在c-2/N处。阻带内产生了余振，最大的负峰在c+2/N处。以上两点就是对hd(n)用矩形窗截断后，在频域的反映，称为吉布斯效应。在主瓣附近，按照(1-11)式，RN()可近似为 式（1-25）我们常使用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、布莱克曼窗和凯塞窗。(1) 矩形窗(Rectangle Window) 其频率响应和幅度响应分别为：， (2) 三角形窗(Bartlett Window) 式（1-26）其频率响应为： 式（1-27） (3) 汉宁(Hanning)窗，又称升余弦窗 式（1-28）其频率响应和幅度响应分别为：式（1-29）(4) 汉明(Hamming)窗，又称改进的升余弦窗 式（1-30）其幅度响应为： 式（1-31）(5) 布莱克曼(Blankman)窗，又称二阶升余弦窗 式（1-32）其幅度响应为： 式（1-33） (6)凯泽(Kaiser)窗 式（1-34）窗函数旁瓣峰值幅度 dB过度带宽阻带最小衰减dB矩型窗-134/N-21三角窗-258/N-25汉宁窗-318/N-44哈明窗-418/N-53布来克曼窗-5712/N-74凯塞窗（a=7.865）-5710/N-80表1-2六种窗函数的基本参数2 MATLAB软件简介 2.1 MATLAB的发展MATLAB是由美国mathwoks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB 产品族可以用来进行以下各种工作：数值分析，数值和符号计算，工程与科学绘图，控制系统的设计与仿真，通讯系统设计与仿真，财务与金融工程，（1）友好的工作平台和编程环境MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。（2）简单易用的程序语言 MATLAB一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。（3）强大的科学计算机数据处理能力MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而前经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C+ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。（4）出色的图形处理功能MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。（5）应用广泛的模块集合工具箱MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。（6）实用的程序接口和发布平台 新版本的MATLAB可以利用MATLAB编译器和C/C+数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C+代码。允许用户编写可以和MATLAB进行交互的C或C+语言程序。（7）应用软件开发（包括用户界面） 在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等。2.2函数命令 (1)fir1函数设计标准响应FIR滤波器可以使用fir1函数。fir1函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通，带通，高通和带阻滤波器。形式为： b=fir1(n,Wc,ftype,Window)各个参数的含义如下：b:滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为： b(z)=b(1)+b(2)+b(n+1) n:滤波器的阶数。 Wc:截止频率。 ftype:当指定ftype时，可设计高通和带阻滤波器。 Window:窗函数。 (2)freqz函数该函数基于FFT算法计算数字滤波器Z变换频率响应。返回数字滤波器的n点复频域响应 H()= 在简单形式中，b，a为滤波器系数，freqz可得到数字滤波器的n点复频响应,并将这n点保存在w中，相应的频率级联在h中。 (3)设计用的其他函数即命令如下表所示：名称功能名称功能Clear从内存中清除变量和函数Close关闭图形Min取最小值Ceil取整Angle相位角Unwrap相位角展开Figure建立图形窗口Subplot在标定位置上建立坐标系Stem离散序列图Plot线性绘图XlabelX轴标记YlabelY轴标记Title图形标题Axis控制坐标系的刻度和形式 表2-2-1 函数命令表3 带通滤波器的设计、matlab程序和matlab图（1) 设计带通滤波器.参数为采样频率Fs=8000HZ，要求保留8001400Hz频段的频率成分，幅度小于0.5dB;滤除0500Hz和3200Hz以上频段的频率成分，衰减大于30dB.解：用数字滤波器对模拟信号进行带通滤波器设计，先对模拟信号进行A/D变换，再进行数字带通滤波器处理。首先确定数字滤波器的技术指标：Wls=2*fpl/Fs=2*800/8000=0.2Wlp=2*fpu/Fs=2*1400/8000=0.35Whp=2*fsi/Fs=2*500/8000=0.125Whs=2*fsu/Fs=2*3200/8000=0.8最小衰减度p=0.5dB.最大衰减度s=30dB.数字带通的技术指标p=0.5db，s=30db，wls=0.2pi，wlp=0.35pi，whp=0.125pi，whs=0.8pi。逼近通带：wls，wlp。逼近阻带:0,whp,whs,模拟带通的技术指标计算T=1 ls=2/T*(tan（ls/2）) p=0.5 db lp=2/T*(tan（lp/2）) s=30 db模拟低通滤波器的技术指标计算 lp=1/（2/T*(tan（lp/2）) ls=1/（2/T*(tan（lp/2）) 。 p=0.5db，s=30db将lp和ls对3db截止频率c归一化，这里c=lp，p=1， s=c/lp设计归一化模拟低通滤波器G（p）N=-lg（ksp）/lg（sp） ksp=（10*（0.1*p）-1）/（10*（0.1*s）sp=s/p查表得到G（p），去归一化，p=s/c，带入得到G（s） 将模拟低通转换成带通，转换到Z域得H（z）=G（T/2*(1+z*（-1）/（1-z*（-1）调用Matlab信号处理工具箱函数设计数字带通滤波器的程序为： clear; close all; wls=0.2*pi;wlp=0.35*pi; whp=0.125*pi;whs=0.8*pi; delta_w=min(wlp-wls),(whs-whp); wc1=(wls+wlp)/2;wc2=(whp+whs)/2; N1=ceil(1.8*pi/delta_w); hn1=fir1(N1-1,wc1,wc2/pi,boxcar(N1); h1,w1=freqz(hn1,1);%Hamming（汉宁）窗N2=ceil(6.6*pi/delta_w); hn2=fir1(N2-1,wc1,wc2/pi,hamming(N2);h2,w2=freqz(hn2,1);%Blackman （布莱克曼）窗N3=ceil(11*pi/delta_w); hn3=fir1(N3-1,wc1,wc2/pi,blackman(N3);h3,w3=freqz(hn3,1);%Kaiser（凯塞贝塞尔） 窗N4=ceil(10*pi/delta_w); hn4=fir1(N4-1,wc1,wc2/pi,kaiser(N4);h4,w4=freqz(hn4,1);%绘图figure(1) %第一幅图 subplot(3,1,1); n=0:N1-1;stem(n,hn1,.); axis(0,N1-1,-0.4,0.4); xlabel(n);ylabel(h(n);grid on; title(boxcar矩形窗单位冲击响应h(n); subplot(3,1,2);plot(w1/pi,20*log10(abs(h1); axis(0,1,-150,5); xlabel(数字域归一化角频率);ylabel(幅度);grid on; title(boxcar矩形窗幅频响应); subplot(3,1,3);plot(w1/pi,180/pi*unwrap(angle(h1); xlabel(数字域归一化角频率);ylabel(单位：度);grid on;title(boxcar矩形窗相频相应);figure(2) %第二幅图 subplot(3,1,1);n=0:N2-1;stem(n,hn2,.); axis(0,N2-1,-0.2,0.2); xlabel(n);ylabel(h(n);grid on;title(Hamming汉宁窗单位脉冲响应h(n);subplot(3,1,2);plot(w2/pi,20*log10(abs(h2); axis(0,1,-150,5);xlabel(数字域归一化角频率);ylabel(幅度（单位：分贝）);grid on;title(Hamming汉宁窗幅频响应);subplot(3,1,3);plot(w2/pi,180/pi*unwrap(angle(h2); xlabel(数字域归一化角频率);ylabel(单位：度);grid on;title(Hamming汉宁窗相频相应);figure(3) %第三幅图subplot(3,1,1);n=0:N3-1;stem(n,hn3,.);axis(0,N3-1,-0.2,0.2);xlabel(n);ylabel(h(n);grid on;title(Blackman布莱克曼窗单位冲击响应h(n);subplot(3,1,2);plot(w3/pi,20*log10(abs(h3);axis(0,1,-150,5);xlabel(数字域归一化角频率);ylabel(幅度（单位：分贝）);grid on;title(Blackman布莱克曼窗幅频响应);subplot(3,1,3);plot(w3/pi,180/pi*unwrap(angle(h3);xlabel(数字域归一化角频率);ylabel(单位：度);grid on;title(Blackman布莱克曼窗相频相应);figure(4) %第四幅图subplot(3,1,1);n=0:N4-1;stem(n,hn4,.);axis(0,N4-1,-0.2,0.2);xlabel(n);ylabel(h(n);grid on;title(Kaiser凯塞贝塞尔窗单位脉冲响应h(n);subplot(3,1,2);plot(w4/pi,20*log10(abs(h4);axis(0,1,-150,5);xlabel(数字归一化角频率);ylabel(幅度（单位：分贝）);grid on;title(Kaiser凯塞贝塞尔窗幅频响应);subplot(3,1,3);plot(w4/pi,180/pi*unwrap(angle(h4);xlabel(数字域归一化角频率);ylabel(单位：度);grid on;title(Kaiser凯塞贝塞尔窗相频相应);4.结果分析 设计结果中，矩形窗和Hamming 窗的情况结果如图所示，Blackman窗和Kaiser窗的运行结果如图所示。根据理论分析，幅频响应图中，在3db的对应点应该为通带截止频率，在60db处对应的是阻带截止频率。 （1）对于矩形窗：窗宽N=12，h(n)为偶对称，对称中心为n=5.5，由于n为整数，故在n=5和n=6处存在两个极大值；在幅频响应图中，实际设计的低端，高端通带截止频率为0.3262pi和0.6758pi，与技术指标相差6.8%和3.9%。而低端和高端的阻带截止频率为0.167pi和0.833pi，与技术指标相差16.5%和4.13%。其阻带的纹波较大，第一阻带最小衰减27db；由相频特性可以看出，次滤波器是线性相位的。 （2）对于Hamming窗：窗宽N=44，h(n)为偶对称，对称中心为n=21.5，由于n为整数，故在n=21和n=22处存在两个极大值；在幅频响应图中，实际设计的低端，高端通带截止频率为0.29pi和0.705pi，与技术指标相差17.17%和8.47%。而低端和高端的阻带截止频率为0.1953pi和0.8047pi，与技术指标相差2.35%和0.59%。第一阻带最小衰减50db；由相频特性可以看出，次滤波器是线性相位的。 （3）对于Blackman窗：窗宽N=80，h(n)为偶对称，对称中心为n=39.5，由于n为整数，故在n=39和n=40处存在两个极大值；在幅频响应图中，实际设计低端，高端通带截止频率为0.289pi和0.7109pi，与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.207pi和0.793pi，与技术指标相差3.5%和0.875%。第一阻带最小衰减75db；由相频特性可以看出，滤波器为线性相位的。 （4）对于Kaiser窗：窗宽为N=67，h(n)偶对称，对称中心n=33，有用n为整数，故在n=33处存在一个极大值；在幅频响应图中，实际设计的低端和高端通带的截止频率为0.289pi和0.7109pi，与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.205pi和0.7969pi，与技术指标相差2.55%和0.39%。第一阻带最小衰减80db；有相频特性可以看出，次滤波器为线性相位的。5设计心得数字滤波是语音和图象处理、模式识别、谱分析等应用中的一个基本处理算法，有限长单位冲激响应数字滤波( 即 FIR)从理论上以及从实际的有限精度的运算中，都是稳定的，