利用MATLAB仿真软件系统结合窗函数法设计一个数字带通FIR滤波器.doc

武汉理工大学MATLAB课程设计报告课程设计任务书学生姓名： 胡俊 专业班级： 信息sy0901 指导教师： 魏洪涛 工作单位： 信息工程学院 题 目:利用Matlab仿真软件系统结合窗函数法设计一个数字带通FIR滤波器初始条件： 1.数字信号处理基本理论知识 2.信号与系统基本理论知识3.MATLAB编程基础知识4.装有MATLAB的PC机要求完成的主要任务: 分别使用矩形窗、三角形窗、汉明窗、布莱克曼窗、凯泽窗进行设计，并输出滤波器的频率特性。时间安排：第20周理论设计、实验室安装调试，地点： 鉴主15楼通信实验室一序号阶段内容所需时间1方案设计1天2软件设计2天3系统仿真1天4答辩1天总计5天指导教师签名： 年 月 日系主任（或责任教师）签名： 年 月 日武汉理工大学MATLAB课程设计报告目录中文摘要1Abstract21 MATLAB的概况32 MATLAB结合窗函数设计法原理43 各类窗函数简介73.1 基本窗函数73.1.1 矩形窗函数73.1.2三角窗73.2 广义余弦窗84 方案设计与论证94.1 程序设计法105 窗函数仿真结果分析146总结与体会177参考文献18武汉理工大学MATLAB课程设计报告中文摘要现代图像、语声、数据通信对线性相位的要求是普遍的。正是此原因，使得具有线性相位的FIR数字滤波器得到大力发展和广泛应用。在实际进行数字信号处理时，往往需要把信号的观察时间限制在一定的时间间隔内，只需要选择一段时间信号对其进行分析。这样，取用有限个数据，即将信号数据截断的过程，就等于将信号进行加窗函数操作。而这样操作以后，常常会发生频谱分量从其正常频谱扩展开来的现象，即所谓的“频谱泄漏”。当进行离散傅立叶变换时，时域中的截断是必需的，因此泄漏效应也是离散傅立叶变换所固有的，必须进行抑制。而要对频谱泄漏进行抑制，可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。而在后面的FIR滤波器的设计中，为获得有限长单位取样响应，需要用窗函数截断无限长单位取样响应序列。另外，在功率谱估计中也要遇到窗函数加权问题。由此可见，窗函数加权技术在数字信号处理中的重要地位。传统的数字滤波器的设计过程复杂，计算工作量大，滤波特性调整困难，影响了它的应用。本文介绍了一种利用MATLAB信号处理工具箱（Signal Processing Toolbox）快速有效的设计由软件组成的常规数字滤波器的设计方法。本文给出了使用MATLAB语言进行程序设计和利用信号处理工具箱的FDA Tool工具进行界面设计的方法。利用MATLAB设计滤波器，可以随时对比设计要求和滤波器特性调整参数，直观简便，极大的减轻了工作量，有利于滤波器设计的最优化。AbstractA modern image, Voice and data communications to the linear phase requirement is universal. It is for this reason, making a linear phase FIR digital filters are to develop and widely used.In practice, digital signal processing, they often need to observe the time limit for the signal in a certain time interval, only need to select a period of time to analyze the signals. In this way, access to a finite number of data about the process of signal data truncation is equivalent to the signal, the additional window function operation. And this after the operation, often occur from the normal component of the spectrum spread spectrum open to the phenomenon of so-called spectral leakage. When it comes to discrete Fourier transform, the time domain truncation is necessary, therefore, leakage effects are also inherent in discrete Fourier transform, and must be inhibited. But would like to suppress the spectral leakage can be inhibited by the weighted window function equivalent DFT filter amplitude characteristics of side-lobe, or the weighted window function so that the limited length of the extension of the input signal cycle, after the border to minimize the degree of discontinuity The method of implementation. In the back of the FIR filter design, in order to obtain finite sampling units to respond to the need to truncate the infinite length window function with unit sampling response sequence. In addition, power spectrum estimation problem also encountered in the weighted window function. Thus, the weighted windows function in digital signal processing in an important position. Conventional digital filter design process is complex, computing workload, difficulties in adjustment of filter characteristics, affecting its application. This paper introduces a MATLAB Signal Processing Toolbox (Signal Processing Toolbox) quickly and efficiently design formed by the software of conventional digital filter design methods. This article presents use MATLAB language for programming and utilization of signal processing Toolbox FDA Tool for interface design. A filter was designed using MATLAB. You can always compare the design requirements and adjust the parameters of filter characteristics, visual simple and greatly reduced the workload will help filter design optimized.2武汉理工大学MATLAB课程设计报告1 MATLAB的概况MATLAB是矩阵实验室（MatrixLaboratory）之意。除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。MATLAB的基本数据单位是矩阵，它的指令表达式与数学,工程中常用的形式十分相似,故用MATLAB来解算问题要比用C,FORTRAN等语言完相同的事情简捷得多.应用MATLAB这一软件来设计滤波器,与传统的设计思路相比是非常的简便的.具体来说主要有以下几点:1) MATLAB的GUID可以提供给使用者一个进行人机交换的环境,在此环境中,使用者没有必要了解程序内部具体的编写情况,在这个界面上他只要输入他所期望的滤波器的性能指标就可以了。2) 利用MATLAB可以进行系统仿真,也就是说,使用者可以通过计算机进行对数字滤波器的频率响应图据的处理和波形的检测,不像以前必须借助一定的实验器材才可以得到。3) MATLAB内部有丰富的函数可供调用,使用者只需要根据自己的需要查到所需要的函数名,那就可以直接调用使用了.不需要自己把函数详细的编写出来.使编写内部程序变的简单。 MATLAB包括拥有数百个内部函数的主包和三十几种工具包(Toolbox).工具包又可以分为功能性工具包和学科工具包.功能工具包用来扩充MATLAB的符号计算,可视化建模仿真,文字处理及实时控制等功能.学科工具包是专业性比较强的工具包,控制工具包,信号处理工具包,通信工具包等都属于此类. 开放性使MATLAB广受用户欢迎.除内部函数外,所有MATLAB主包文件和各种工具包都是可读可修改的文件,用户通过对源程序的修改或加入自己编写程序构造新的专用工具包.2 MATLAB结合窗函数设计法原理数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是和重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。因此设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲击响应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数： （公式2-1）使用上述的传递函数去逼近，一个理想的频率响应的傅立叶反变换： （公式2-2）其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用来近似。窗函数的基本思想：先选取一个理想滤波器（它的单位抽样响应是非因果、无限长的），再截取（或加窗）它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。设是一个长序列，是长度为N的窗函数，用截断，得到N点序列，即 （公式2-3）在频域上则有 （公式2-4）由此可见，窗函数不仅仅会影响原信号在时域上的波形，而且也会影响到频域内的形状。MATLAB信号工具箱主要提供了以下几种窗函数，如表1-1所示：窗窗 函 数矩形窗Boxcar三角窗Triang 汉明窗Hamming布莱克曼Blackman凯塞窗Kaiser表1-1 MATLAB窗函数加矩形窗后的频谱和理想频谱可得到以下结论：加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。N越大，过渡带越窄、越陡；过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95，与N无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为： （公式2-5）这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中（-0.0895）对应的点的值定义为阻带最小衰耗。以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波器的关键在于改进窗函数。窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为： 旁瓣峰值衰耗20lg(第一旁瓣峰值主瓣峰值) （公式2-6）为了改善滤波器的性能，需使窗函数谱满足：(1)主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带。(2)尽量减少最大旁瓣的相对幅度，也就是能量集中于主瓣，以减小带内、带外波动的最大幅度，增大阻带衰减。一般来说，以上两点很难同时满足。当选取主瓣宽度很窄时，旁瓣的分量势必增加，从而带内、带外的波动也增加了；当选取最小的旁瓣幅度时，降低了带内、带外的波动，但是过渡带的陡度减小了。所以实际采用的窗函数其特性往往是它们的折中，在保证主瓣宽度达到一定要求的前提下，适当牺牲主瓣宽度来换取旁瓣波动的减小。3 各类窗函数简介数字信号处理领域中所用到的窗函数主要分为两大类，分别是：基本窗函数与广义窗函数。3.1 基本窗函数基本窗函数包括以下三类，主要有：矩形窗函数、三角窗函数和巴特利特窗函数。下面就对这些窗函数展开介绍。3.1.1 矩形窗函数矩形窗(Rectangular Window)函数的时域形式可以表示为： （公式3-1）它的频域特性为 （公式3-2）Boxcar函数：生成矩形窗调用方式w = boxcar (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。从功能上讲，该函数又等价于w = ones(n,1)。3.1.2三角窗三角窗（Triangular Window）函数时域形式可表示为： （公式3-3）窗谱为： （公式3-4） 式中，当N远大于1时，此时，窗谱主瓣宽度为8/N。3.2 广义余弦窗汉明窗和布莱克曼窗，都可以用一种通用的形式表示，这就是广义余弦窗。这些窗都是广义余弦窗的特例，汉明窗又被称为改进的升余弦窗，而布莱克曼窗又被称为二阶升余弦窗。采用这些窗可以有效地降低旁瓣的高度，但是同时会增加主瓣的宽度。这些窗都是频率为0、2/(N1)和4/(N1)的余弦曲线的合成，其中N为窗的长度。通常采用下面的命令来生成这些窗： （公式3-5） （公式3-6）其中，A、B、C适用于自己定义的常数。根据它们取值的不同，可以形成不同的窗函数，分别是：汉明窗 A=0.54，B=0.54，C=0；布莱克曼窗 A=0.5，B=0.5，C=0.08；4 方案设计与论证用窗函数法设计一个FIR带通滤波器。指示如下：低端阻带截止频率 wls = 0.2*pi；低端通带截止频率 wlp = 0.35*pi；高端通带截止频率 whp = 0.65*pi；高端阻带截止频率 whs = 0.8*pi；六种窗函数的基本参数如表1-2：窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134/N-21三角形窗-258/N-25汉明窗-418/N-53布莱克曼窗-5712/N-74凯泽窗-5710/N-80表1-2 窗函数基本参数以上表格里的参数设置是最佳窗函数设计，根据设计方案的要求，选择一个合适的窗函数进行滤波器的设计，从上表可以看出：最小带阻衰减仅有窗函数决定，不受N的影响，而过渡带的宽度则随窗函数的增加而减小。选定窗函数类型和长度N之后，求出单位脉冲响应，并按照公式2-1求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足 （公式4-1）根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。4.1 程序设计法主程序框图如图4.1所示。其中幅度特性要求用dB表示。开始读入窗口长度N计算hd(n)调用窗函数子程序求w(n)调用子程序（函数）计算H(k)=DFTh(n)调用绘图子程序（函数）绘制H(k)幅度相位曲线结束图4-1 主程序框图计算h(n)= hd(n) w(n)下面为matlab主程序：%用窗函数法设计FIR数字滤波器b=1;close all;i=0;while(b); temp=menu(选择窗函数长度N,N=10,N=15,N=20,N=25,N=30,N=33,N=35,N=40,N=45,N=50,N=55,N=60,N=64); menu1=10,15,20,25,30,33,35,40,45,50,55,60,64; N=menu1(temp); temp=menu(选择逼近理想低通滤波器截止频率Wc,Wc=pi/4,Wc=pi/2,Wc=3*pi/4,Wc=pi,Wc=0.5,Wc=1.0,Wc=1.5,Wc=2.0,Wc=2.5,Wc=3.0); menu2=pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3; w=menu2(temp); n=0:(N-1); alpha=(N-1)/2; m=n-alpha+eps; hd=sin(w*m)./(pi*m) %得到理想低通滤波器 k=menu(请选择窗口类型:,boxcar,Triang,Hamming,blackman,Kaiser); if k=1 B=boxcar(N); string=Boxcar,N=,num2str(N); else if k=2 B=Triang(N); string=Triang,N=,num2str(N); else if k=3 B=Hamming(N); string=Hamming,N=,num2str(N); else if k=4 B=blackman(N); string=Blackman,N=,num2str(N); else if k=5 B=Kaiser(N); string=Kaiser,N=,num2str(N); end end end end end h=hd.*(B); %得到FIR数字滤波器 H,m=freqz(h,1,1024,whole); %求其频率响应 mag=abs(H); %得到幅值 db=20*log10(mag+eps)/max(mag); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1; stem(n,h,.); axis(0,N-1,-0.1,0.3); hold on; n=0:N-1; x=zeros(N); plot(n,x,-); xlabel(n); ylabel(h(n); title(实际低通滤波器的h(n); text(0.3*N),0.27,string); hold off; subplot(2,2,2); plot(m/pi,db); axis(0,1,-100,0); xlabel(w/pi); ylabel(dB); title(衰减特性(dB); grid; subplot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,-); title(相频特性); xlabel(频率(rad); ylabel(相位(rad); axis(0,3.15,-4,4); subplot(2,2,4); plot(m,mag);