利用MATLAB结合频率取样法设计数字高通FIR滤波器Word下载.docx

1、它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。本文介绍了如何利用MATLAB仿真软件系统及数字信号处理所学知识利用频率采样法设计一个数字高通滤波器。以此来巩固课堂理论学习，并能用所学理论知识正确分析信号处理的基本问题和解释信号处理的基本现象。关键字： MATLAB； 数字信号处理； 数字滤波器； 频率采样法AbstractMAT

2、LAB is released by the United States mathworks mainly for scientific computing, visualization and interactive program designed high-tech computing environment. It numerical analysis, matrix computation, scientific data visualization as well as non-linear dynamic systems modeling and simulation, and

3、many other powerful integrated in an easy-to-use Windows environment, scientific research, engineering design and the need for effective numerical the edit mode many scientific fields provides a comprehensive solution, and in large part to get rid of the traditional non-interactive programming langu

4、age （such as C, Fortran）, on behalf of the advanced level of todays international scientific computing software.This article describes how to use MATLAB simulation software systems and digital signal processing learned knowledge using frequency sampling method to design a digital high-pass filter. I

5、n order to consolidate the theoretical classroom learning, and basic questions and explain basic signal processing phenomenon can be learned theoretical knowledge to correctly analyze the signal processing.Keyword: digital signal processing； digital filter； frequency sampling method1 FIR数字滤波器1.1 FIR

6、滤波器的特点FIR滤波器的脉冲响应h（n）是有限长的（0nN-1），其z变换为： （式1）它是z-1的（N-1）阶多项式，在有限z平面（0n）上有（N-1）个零点，而极点位于z平面原点z=0处，且有（N-1）阶。FIR滤波器的基本结构可以理解为一个分节的延时线，把每一节的输出加权累加，可得到滤波器的输出，FIR滤波器的冲激响应h（n）是有限长的，数学上M阶FIR滤波器可以表示为： y（n）= （式2）其系统函数为： H（z）= （式3） 普通的直接型FIR 滤波器结构如图1 所示。图1 FIR滤波器的直接型结构FIR滤波器最突出的优点有2个：一是只要对h（n）附加一定的条件，很容易获得严格的线

7、性相位特性；二是由于H（z）的极点位于原点z=0处，始终满足稳定条件，所以FIR滤波器永远稳定。三是FIR滤波器由于单位脉冲响应是有限长的，因而可以用快速傅里叶变换（FFT）算法来实现过滤信号，从而可大大提高运算效率。但是，要取得很好的衰减特性，FIR滤波器H（z）的阶次比IIR滤波的要高。1.2 FIR数字滤波器设计方法IIR滤波器设计中的各种变换法对FIR滤波器设计是不适用的，这是因为那里是利用有理分式的系统函数，而FIR滤波器的系统函数只是z-1的多项式。FIR的设计任务是选择有限长度的脉冲响应h（n），得到系统函数H（z），使幅频特性满足技术指标要求，同时使相频特性达到线性相位。 常用

8、设计方法：（1）窗函数法（2）频率采样法（3）切比雪夫等波纹逼近法。人们最感兴趣的是FIR滤波器具有线性相位的相频特性。对非线性相位的FIR滤波器，一般可以用IIR滤波器来代替，因为同样幅度特性，IIR滤波器所需阶数比FIR滤波器的阶数要少得多。1.3 线性相位FIR数字滤波器的条件和特点1.3.1 线性相位条件对于长度为N的h（n），传输函数为 （式4）H（ej）=Hg（）ej（） （式5）式中，Hg（）称为幅度特性，（）称为相位特性。注意，这里Hg（）不同于|H（ej）|，Hg（）为的实函数，可能取负值，而|H（ej）|总是正值。H（ej）线性相位是指（）是的线性函数，即 （）= - ，为

9、常数 （式6）如果（）满足 （）= 0- ，0是起始相位 严格地说，此时（）不具有线性相位，但以上两种情况都满足群时延是一个常数，即 （式7）也称这种情况为线性相位。1.3.2 线性相位FIR滤波器的幅度特性与相位特性线性相位FIR滤波器的幅度特性与相位特性如下图：图2线性相位FIR滤波器的幅度特性与相位特性一览表在设计时，要注意选择合适的h（n）对称形式（奇或偶）和h（n）长度N（奇数或偶数）。如要设计高通滤波器，只能选情况1和情况4；要设计低通滤波器，只能选情况1和情况2。2 利用频率采样法设计FIR滤波器2.1 用频率采样法设计滤波器的基本原理 待设计的滤波器的传输函数用Hd（ej）表示

10、，可按下列思路进行设计：1 它在=0到2之间等间隔采样N点，得到Hd（k） （式8）2 N点Hd（k）进行IDFT，得到h（n） （式9）式中，h（n）作为所设计的滤波器的单位取样响应。3 h（n）求系统函数H（z） （式10）将插值公式重写如下 （式11）此式就是直接利用频率采样值Hd（k）形成滤波器的系统函数。用频率采样法设计线性相位滤波器的条件 ：FIR滤波器具有线性相位的条件是h（n）是实序列，且满足h（n）= h（N1n），其传输函数应满足的条件是 （式12） （式13） （式14） （式15）且Hg（）=0 。在=02之间等间隔采样N点，将=k代入式（47）中，并写成k的函数： （

11、式16） （式17） ，N为奇数（式18） ，N为偶数且 （式19） （式20）说明N等于奇数时Hg（k）对（N1）/2偶对称，N等于偶数时, Hg（k）对N/2奇对称，且Hg（N/2）=0。对于高通滤波器，这里N只能取奇数。截止频率为c，采样点数N，Hg（k）和（k）用下面公式计算 （式21）以上是用频率采样法设计滤波器的基本原理。2.2 线性相位的约束条件以h（n）为偶对称，N为奇数的情况进行分析。1）FIR的频响具有线性相位的一般表达式当h（n）为偶对称，N为奇数时，则 （式22） 而且幅度函数H（w）应为偶对称，即 （式23） 2）采样值H（k）具有线性相位的约束 （式24） 其中，

12、表示采样值的模（纯标量），表示其相角。因此，在采样点上具有线性相位的条件应为： （式25）而且，必须满足偶对称，即： （式26） 实际滤波器的传输函数，与理想的传输函数Hd（ej）间存在误差，如图2图3频率采样的响应需要讨论逼近误差问题及其改进措施。2.3 逼近误差及其改进措施2.3.1 产生误差的原因从图3可看出，实际的H（ej）与理想的Hd（ej）相比，误差主要体现在一是通带和阻带出现波动，二是过渡带加宽，与窗函数设计法情况类似，产生误差的原因可从时域和频域两方面进行分析。从时域分析：如果Hd（ej）有间断点，那么相应单位取样响应hd（n）应是无限长的。这样，由于时域混叠，引起所设计的h（

13、n）和hd（n）有偏差。为此，希望在频域的采样点数N加大。N愈大，设计出的滤波器愈逼近待设计的滤波器Hd（ej）。从频域分析：在采样点=2k，k=0,1,2,N-1,（-2k/N）=1，因此，采样点处H （ejk） （k=2k/N）与H（k）相等，逼近误差为0。在采样点之间，H（ej）由有限项的H（k）（-2k/N）之和形成。其误差和Hd（ej）特性的平滑程度有关，特性愈平滑的区域，误差愈小；特性曲线间断点处，误差最大。表现形式为间断点用倾斜线取代，且间断点附近形成振荡特性，使阻衰减减小，往往不能满足技术要求。2.3.2 减小误差的方法 最直观的想法是增加采样点数，即加大N值，由于过渡带就等于

14、采样间隔（参看图3），即 （式27）所以加大N，可使过渡带变窄，但增加要适当，否则会增加滤波器体积与成本。但是，增加N并不会改善滤波器的阻带衰减特性，因为Hd（ej）是理想矩形, 无论怎样增多频率采样的点数，在通、阻带交界处，幅值总是从1突变到0，会引起较大的起伏振荡。为使逼近误差更小，和窗口法的平滑截断一样，通过在理想频率响应的不连续点的边缘上加一些过渡的抽样点，减小频带边缘的突变，也就减小了起伏振荡，增大了阻带最小衰减。一般过渡带取一、二、三点抽样值即可得到满意结果。如在低通设计中，不加过渡点时，阻带最小衰减为-20dB，加三个过渡点（最优设计）则可达-80dB到-95dB左右。加过渡点的

15、示意如图4所示。图4理想低通滤波器增加过渡点增加过渡点，可使阻带衰减明显提高，但付出的代价是过渡带加宽，可通过下式加大N来调整。 m=0，1，2，3 （式28）2.4 频率采样法的特点 频率采样法设计滤波器最大的优点是直接从频率域进行设计，比较直观，也适合于设计具有任意幅度特性的滤波器。但边界频率不易控制。如果增加采样点数N，对确定边界频率有好处，但会增加滤波器的成本。因此，它适合于窄带滤波器的设计。3 频率取样法的数字高通滤波器的实现3.1 MATLAB的介绍MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数

16、据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。3.2 设计条件ws=0.6pi, wp=0.8pi, 通带波动1dB, 阻带衰减50dB，M=33。3.3 设计程序% 频率采样技术: 高通% ws=0.6pi, wp=0.8pi, Rp=1dB, As=50dB% M=33, T1 = 0.1095; T2 = 0.598;M = 33; alpha = （M-1）/2; l = 0:M-1; wl = （2*pi/M）*l;T1 = 0.1095;Hrs = zeros（1,11）,T1,T2,ones（1,8）,T2,T1,z

17、eros（1,10）;Hdr = 0,0,1,1; wdl = 0,0.6,0.8,1;k1 = 0:floor（M-1）/2）; k2 = floor（M-1）/2）+1:angH = -alpha*（2*pi）/M*k1, alpha*（2*pi）/M*（M-k2）;H = Hrs.*exp（j*angH）;h = real（ifft（H,M）;db,mag,pha,grd,w = freqz_m（h,1）;Hr,ww,a,L = Hr_Type1（h）;subplot（1,1,1）subplot（2,2,1）;plot（wl（1:17）/pi,Hrs（1:17）,o,wdl,Hdr）;a

18、xis（0,1,-0.1,1.1）; title（高通: M=33,T1=0.1095,T2=0.598）xlabel（）; ylabel（Hr（k）set（gca,XTickMode,manualXTick,0;.6;.8;1）XTickLabelModeXTickLabels, 0;.6.8 1）YTickModeYTick,0,0.109,0.59,1）; gridsubplot（2,2,2）; stem（l,h）; axis（-1,M,-0.4,0.4）title（脉冲响应h（n）text（M+1,-0.4,nsubplot（2,2,3）; plot（ww/pi,Hr,wl（1:振幅响

19、应频率（单位：pi）Hr（w）subplot（2,2,4）;plot（w/pi,db）; axis（0,1,-100,10）;幅度响应 xlabel（ylabel（分贝数Manual,-50;0）;YTickLabelModeYTickLabels50%3.4 调试结果 图5 频率采样技术：高通，最优法结果分析：第一幅图为要高通滤波器原型，可以看到它在过渡带添加了两个采样点，以增加阻带衰减；第二幅图为系统函数单位脉冲响应图形，可以看出，它以中点成偶对称，由于采样点数为奇数，故在对称轴处有取值；第三幅图（左下）为根据频率取样法设计出的滤波器振幅响应，可以看出它在采样点处的取值与原高通滤波器精确一

20、致，在其他点处与原高通滤波器取值逼近有上下波动；第四幅图为用分贝数表示的幅度响应，可以看到采用线性最优法设计的高通滤波器的阻带衰减大于50db。设计取得了良好的效果。4 心得体会Matlab的课程设计做到现在已经基本接近尾声了，既然学习一门课程，简单的总结是必须要有的。以前在信号与系统和数字信号处理的实验中已经接触过matlab，所以上手并不是很难，不过在设计的时候还是遇到了不少问题，首先是对频率取样法掌握的不到位，重新学习了频率取样法后，发现如何利用程序实现频率取样法成了一个问题。通过自己在网上查找资料，看从图书馆借来的书以及对照着老师的PPT，不断的调试，终于做出了成果。课程设计虽然做完了

21、，但现在学的这点知识还远远不够，特别是这个软件的函数非常多，要能够熟练运用我们还有很多要学习。不过我觉得Matlab的函数设计都比较合理，她总是从函数本身的意义出发命名，这使我们记不会很难。 总之这次课程设计完成的还算顺利，虽然也遇到过一些问题，但通过和同学讨论一起学习都能解决。当然，我们也都明白matlab的确是一个很实用的工具，在今后的学习中我们会不断的边学边运用它，而且我们还可以将它用在我们专业的学习中。5 参考文献1 刘 泉，数字信号处理原理与实现，电子工业出版社，20092 郭仕剑，MATLAB7.x数字信号处理，人民邮电出版社，20073 陈怀琛，MATLAB及在电子信息课程中的应

22、用，电子工业出版社，20064 高会生，MATLAB实用教程（第2版），电子工业出版社，2010 5 陈怀琛，数字信号处理教程MATLAB释义与实现，电子工业出版社，2004附录 辅助函数1function Hr,w,a,L = Hr_Type1（h）;% 计算第一种低通滤波器设计的振幅响应Hr（w） % -% Hr,w,a,L = Hr_Type1（h）% Hr = 振幅响应% w = 在0 pi 区间内计算Hr 的500个频率点% a = 第一种低通滤波器的系数% L = Hr的阶次% h = 第一种低通滤波器的频率响应 M = length（h）; L = （M-1）/2; a = h（

23、L+1） 2*h（L:-1:1）; % 1乘（L+1） 行向量row vector n = 0:1:L; % （L+1）乘1 列向量 w = 0:500*pi/500;Hr = cos（w*n）*a辅助函数2function db,mag,pha,grd,w = freqz_m（b,a）;% freqz 子程序的改进版本% -% db,mag,pha,grd,w = freqz_m（b,a）;% db = 0 到pi弧度区间内的相对振幅（db）% mag = 0 到pi弧度区间内的绝对振幅% pha = 0 到pi弧度区间内的相位响应% grd = 0 到pi弧度区间内的群迟延% w = 0

24、到pi弧度区间内的501个频率样本向量% b = Ha（z）的分子多项式系数（对FIR b=h）% a = Ha（z）的分母多项式系数（对 FIR: a=1）H,w = freqz（b,a,1000,whole H = （H（1:501） w = （w（1: mag = abs（H）; db = 20*log10（mag+eps）/max（mag）; pha = angle（H）;% pha = unwrap（angle（H）; grd = grpdelay（b,a,w）;% grd = diff（pha）;% grd = grd（1） grd;% grd = 0 grd（1:500）; grd; grd（2:501） 0;% grd = median（grd）*500/pi;