课程设计_窗函数法设计FIR_数字滤波器.doc

课程设计（论文）任务书院（系）： 基层教学单位： 学 号学生姓名专业（班级）设计题目利用窗函数法设计FIR数字滤波器设计技术参数窗函数，通带截频，阻带截频，通带最大衰减，阻带最大衰减幅度响应，相位响应等设计要求(1) 给定所要求的频率响应Hd(ej)；(2) 由Hd(ej)的傅立叶反变换求出Hd(n)；(3) 根据过渡带宽及阻带最小衰减的要求,选择窗的形状并估计窗口长度N；(4) 计算滤波器的单位脉冲响应h(n) = Hd (n)(n)；(5)用matlab编程实现，显示其幅频特性曲线，分析并得出结论。工作量共二周，每天工作8小时工作计划第一到三天：查资料并选定题目，完成设计方案第三到七天：根据题目编程并调试第八到十 四天：完成论文参考资料陈后金 数字信号处理(2008年11月第二版)高等教育出版社陈后金 信后与系统(2007年12月1日第1版)高会生 MATLAB原理与工程应用(2006年1月第1版) 电子工业出版社丁玉美 数字信号处理（2000年12月第二版）西安电子科技大学出版社指导教师签字基层教学单位主任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。年 月 日 课程设计评审意见表指导教师评语：成绩： 指导教师： 年 月 日答辩小组评语：成绩： 评阅人： 年 月 日课程设计总成绩：答辩小组成员签字：年 月 日 课 程 设 计 说 明 书利用窗函数法设计FIR数字滤波器Design FIR digital filters with the window function摘要: 在数字信号处理中, 数字滤波器是一种被广泛使用的信号处理部件。分析FIR (有限冲激响应) 数字滤波器的结构特征, 得到了满足系统要求的数字滤波器设计法，结合实际工程所要求的数字滤波器指标, 利用MATLAB 对FIR 数字滤波器进行了设计和仿真, 最后通过讨论,给出窗函数法FIR 滤波器的设计原则。关键词：FIR 数字滤波器;窗函数;Matlab程序;设计Abstract: The digital filter is a widely used signal processing section in digital signal processing,. Analysis of the structural characteristics of the FIR (finite impulse response) digital filters, digital filter design method has been to meet the system requirements, combined with digital filter indicators of practical engineering required, using MATLAB to design and simulation FIR digital filter, and finally given window function FIR filter design principles through discussion.Keywords: FIR digital filter; window function; Matlab program; design一.设计目的1掌握FIR数字滤波器的设计方法和步骤，深入了解其幅度和相位特性，并能在实际中学会选择与应用。2学习Matlab语言工具，并应用Matlab实现数字信号处理的各种方式，理解和熟悉数字信号处理的实际应用。3了解几种窗函数的性质，熟悉其波形图，分析加窗截断对波形的影响。二设计原理1FIR数字滤波器结构特征FIR 数字滤波器的冲激响应只能延续一定的时间,N 阶因果有限冲激响应滤波器的差分方程表达式为: yk=k=0Mhix(k-i)其卷积和是一个有限和,可以直接计算y( k) ,所以其涉及的基本操作就是简单的乘和加,其计算仅需要所给出的初始样本值及所要求的相关样本值。通过精确设计,有限冲激响应滤波器在整个频率范围内均能提供精确的线性相位,由于系统的单位脉冲序列为有限 共 9 页 第 1 页序列,当输入有限时,输出也必然为有限,这样其稳定性总可以独立于滤波器系数之外。因此,在很多情况下,有限冲激响应滤波器成为首选,只要确定能满足要求的转移序列或者脉冲响应的常数,就可以准确地设计出满足要的FIR数字滤波器。通常所采用的设计方法主要有窗函数、频率采样法和等波纹最佳逼近法等,其中窗函数法是从时域进行设计的,其算法简单、物理意义清晰,因此得到了较为广泛的应用。2算法研究为了建立一个具有线性相位和稳定的非递归特性的有限脉冲响应滤波器即 FIR滤波器 , 要考虑两个方 面 : 一是使用有限长的单位取样响应来逼近理想低通，二是单位取样响应对(N -1)/2 对称,从而保证线性相位。根据FIR 数字滤波器输入输出关系的差分方程描述 ,有yx=i=0n-1hix(n-i)上式反映了有限列长的单位取样响应。滤波网络传输函数为:Hz=n=0N-1hd(n)z-n利用傅立叶变换得单位取样响应hd(n)为:hdn=12-Hd(ej)ejnd=sin(wc(n-a)(n-a)wc是截止频率,a =(N-1)/2 ,从而保证线性相位。在窗函数法设计中使用有限列长的 h(n) 逼近hdn，用Hamming窗函数w ( n )将hd(n)截断，进行加权处理后可得hn=hdnw(n)。此时频率响应函数为 ：Hdej=n=0N-1h(n)e-jn Hamming 窗函数wh(n)= 0.54 - 0.46cos2nN-1 RN(n)，Blackman 窗函数wb(n)= 0.42-0. 5cos2nN-1+ 0. 08cos4nN-1。其中RN(n) 是幅度为1 ,长度为N的矩形序列。3窗函数法的原理 可以从时域或频域出发来设计FIR 滤波器,从频域出发的方法称为频率采样设计法,从时域出发的方法称为窗函数法,也称为傅立叶级数法,这种方法应用较广泛。窗函数法的设计思想是按照所要求的理想滤波器频率响应Hdej,设计一个FIR 滤波器,使之频率响应 共 9 页 第 2 页Hej=n=0N-1hne-jn来逼近Hdej。因为设计是在时域中进行的, 需先由Hdej的傅立叶反变换导出序列hdn，即hdn=12-Hdejejnd由于Hdej是矩形频率特性，有hd(n)是一无限长的序列,且是非因果的,而要设计的FIR滤波器的冲激响应序列是有限长的, 所以要用有限长的序列h(n) 来逼近无限长的序列hdn，最有效的方法是截断hdn，或者说用一个有限长度的窗口函数w(n) 序列来截取hdn，即hn=w(n) hdn按照复卷积公式, 在时域中的乘积关系可表示成在频域中的周期性卷积关系, 即可得所设计的FIR滤波器的频率响应Hej=12-HdW( )d其中, W()为截断窗函数的频率特性。由此可见,实际的FIR 数字滤波器的频率响应Hej逼近理想滤波器频率响应Hdej的好坏，完全取决于窗函数的频率特性W()。4窗函数法的设计步骤(1)给定所要求的频率响应Hdej；(2)由Hdej的傅立叶反变换求出hdn；(3)根据过渡带宽及阻带最小衰减的要求,选择窗的形状并估计窗口长度N；(4)计算滤波器的单位脉冲响应hn=w(n) hdn(5)验算技术指标是否满足要求。窗函数设计法设计时要满足以下两个条件:(1)窗谱主辨尽可能窄,以获得较陡的过渡带；(2)尽量减少窗谱的最大旁瓣的相对幅度,使能量尽量集中于主辨,进而增加阻带的衰减。给定的滤波器指标参数一般为通带截止频率p 、阻截止频率s、实际通带波动R p和最小阻带衰减A s。窗函数设计的经验公式为归一化过渡带 =s-p2滤波器阶数 M=A s-7.9514.36当A s50时，=0. 1102(A s-8.7)；当21A s50时，=0. 5842(A s-21)0.4+0. 07886(A s-21) 共 9 页 第 3 页在实际工程中常用的函数有六种, 海明窗(Hamming) 、矩形窗(Boxcar) 、三角窗(Bartlett ) 、汉宁窗( Hanning) 、布拉克曼窗(Blackman) 和凯泽窗( Kaiser) 等。三设计内容1低通滤波器的设计指标。通带截止频率p=0.2,阻带截止频率s=0.3,最小阻带衰减A s=50dB 。2设计分析由以上参数可得过度带宽tr_width = ws - wp , wc = (ws + wp)/2。考虑选择Hamming 窗函数, 则列长N = 8*pi/tr-width , hd(n)=sin (wc*(n-(N-1)/2)/(pi*(n-(N-1)/2), h(n)=hd (n)*wh(n)。利用MATLAB编程时,先计算N、tr-width、hd(n)和h(n) ,并利用MATLAB提供的相关工具计算该低通滤波器的幅度特性值, 最终显示Hamming 窗频谱、实际单位取样响应频谱、理想单位取样响应频谱和该FIR滤波器幅度特性曲线。3编写Matlab程序(1) 实现FIR滤波器的主程序:%dfilter.m文件%clearwp = 0.2*pi; ws = 0.3*pi; %给出通带和阻带的截止频率tr-width = ws - wp %过度带宽度N = ceil(8*pi/tr_width)+1 %计算单位取样响应列长 ,使其偶对称n = 0:1:N-1 ;wc =(ws + wp )/2 %计算低通截止频率hd =ideal_lp(wc ,N); %计算理想低通单位取样响应w_ham =(hamming(N);h =hd.*w_ham ; %应用窗函数加权后的有限长单位取样响应序列db, mag ,pha ,grd ,w = freqz_m(h,1) ; %计算幅值响应 、相位响应和群延迟响应delta_w =2*pi/1000 ;Rp =-( min(db(1:1:wp/delta_w+1)As =-(round(max(db(ws/delta_w+1:1:501)figure(1);subplot(2,2,2);stem(n, w_ham)title(Hamming Window)axis(0 N 0 1.1 );ylabel(w( n) %显示 Ha m mi ng 窗频谱subplot(2,2,3); stem(n,h)title(Actual Impulse Response) 共 9 页 第 4 页axis(0 N - 0.1 0.3); ylabel(h(n) ) %显示实际单位取样响应频谱subplot(2,2,4);plot(w/pi ,db)title(Magnitude Response in db);gridset( gca , XTickMode, manual, XTick,0,0.2,0.3,1)axis(0,1,-100,20);ylabel(Decibels) %显示 F IR 滤波器幅度特性曲线subplot(2,2,1);stem(n,hd)title(Ideal ImpulseRe sponse)axis(0 N -0.1 0.3);ylabel(hd(n) %显示理想单位取样响应频谱(2)计算理想低通取样响应的程序:%ideal_lp.m文件function hd = ideal_lp(wc,N) ;%Ideal Lowpass filter computation% - - - - - - - - - - - - - - - - - - - - - -%hd= ideal_lp(wc,M)%hd =ideal impulse response between 0 to M-1%wc=cutoff frequency in redians%M =length of the ideal filter%alpha=(N - 1)/2 ;n=0:1 :N-1 ;m =(n-alpha+eps);hd =sin(wc*m)./(pi*m);end(3)计算幅值响应、相位响应和群延迟响应% freqz_m.m文件function db,mag,pha,grd,w = freqz_m(b,a) ;%Modified version of freqz subroutine% - - - - - - - - - - - - - - - - - - - - - - - - - -%db,mag,pha,grd,w =freqz_m(b,a)%db = Relative magnitude in dB computed over 0 to pi radians%mag = absolute magnitude computed over 0 to pi radians%p ha = Phase response in radians over 0 to pi radians%grd = Group delay over 0 to pi radians%w = 501 f requency samples between 0 and pi radians%b = nuberrator polynomial of H(z) (for F IR :b = h)%a = denominator polynomial of H(z) (for a = 1 共 9 页 第 5 页%H,w = freqz(b,a,1000 , whole) ;H =(H(1 :501);w =(w(1:501);mag =abs(H);db = 20*log10(mag + eps) ;pha=angle(H) ;grd =grpdelay(b,a,w);end运行后,有关参数:N=81,tr_width=0.3142,wc =0.7854 ,As = 53,相应幅频特性曲线为图1。当改变参数N = 201 时,其幅频特性曲线为图2。使用Blackman 窗函数设计该滤波器, 可将主程序中的w- ham = (hamming (N) ) ; 改为w_ham =(Blackman (N); ,N = 81 ,相应幅频特性曲线为图3 (图14中的子图依次为理想取样响应、窗函数、实际取样响应频谱和低通幅度特性曲线) 。图（1） Hamming N=81图（2）Hamming N=200图（3）Blackman N=81图(1) Hamming N=81 共 9 页 第 6 页图（4）Blackman N=200图(2) Hamming N=200图(3) Blackman N=81 共 9 页 第 7 页四结果分析：在前面的FIR滤波器设计中,一是利用Hamming窗函数并且使用了不同的有限列长值(N) ,如图1和图2 ,二是利用不同的窗函数(Hamming 窗和Blackman 窗) 使用同样的N 值(N = 81) ,如图1 和图3。通过分析图1 和图2,注意到相同的窗函数下,N 值越大,主瓣越窄,过度带越窄(在图中就是p=0.2 处的下降沿的宽度), 阻带的波动频率加快, 但阻带的最小衰减没有改变(均接近-60dB) , 即最大的副瓣相对于主瓣的值没有改变。虽然提高N 值没有加大阻带衰减,但却使过度带变窄了。再来分析图1 和图3。它们一个用的是Hamming窗, 一个用的是Black