FIR带通滤泼器的设计.doc

1、数字信号处理课程设计报告题 目： FIR带通滤泼器的设计 学 院： 沈阳工业大学工程学院 专 业： 测控技术与仪器 班 级： 0901 姓 名： 指导教师： 2012年 6月24日目 录1 课程设计要求32 基本原理3 2.1 FIR滤泼器简介3 2.2 窗函数法原理43 建立模型概述5 3.1 MATLB 常用函数5 3.1.1 窗函数5 3.1.2 fir1函数6 3.1.3 freqz函数7 3.1.4 其他函数与命令7 3.2 程序流程图84 模块功能分析95调试过程及结论136 心得体会177 参考文献17 1 课程设计要求 用窗函数法设计FIR带通滤波器。要求低端阻带截止频率ls=

2、0.2，低端通带截止频率lp=0.35，高端通带截止频率hp=0.65，高端阻带截止频率hs=0.8。绘出h(n)及其幅频响应特性曲线。2 基本原理2.1 FIR滤泼器简介 数字滤波器包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。在现代信号处理技术中，例如数据传输、雷达接收以及一些要求较高的电子系统，都越来越多地要求信道具有线性的相位特性。在这方面，FIR滤波器具有独到的优点，它可以在幅度特性随意设计的同时，保证精确、严格的线性相位特性。FIR滤波器的单位脉冲响应h（n）是有限长的（0nN-1），其z变换为的（N-1）阶多项式：可得FIR滤波器的系统差分方程为：因

3、此，FIR滤波器又称为卷积滤波器。FIR滤波器的频率响应表达式为：信号通过FIR滤波器不失真条件是在通带内具有恒定的幅频特性和线性相位特性。理论上可以证明：当FIR滤波器的系数满足下列中心对称条件： 或者 时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个 N 阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。FIR滤波器的设计任务是选择有限长度的h(n)，使传输函数满足技术要求。FIR滤波器的设计方法有多

4、种，如窗函数法、频率采样法及其它各种优化设计方法，本次设计使用窗函数法设计FIR带通滤波器。2.2 窗函数法原理 如果所希望的滤波器的理想频率响应函数为，则其对应的单位脉冲响应为：窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼近hd(n)。由于hd(n)往往是无限长序列，且是非因果的，变成物理上可实现呢？一个自然的想法是只取其中的某些项，即只截取hd(n)中的一部分，比如n=0，N-1，N为正整数。这种处理相当于将hd(n)与函数w(n)相乘，w(n)具有下列形式：w(n)相当于一个矩形，我们称之为矩形窗。即我们可采用矩形窗函数w(n)将无限脉冲响应hd(n)截取一段h(n)来近似

5、为hd(n)，这种截取在数学上表示为：这里应该强调的是，加窗函数不是可有可无的，而是将设计变为物理可实现所必须的。h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数Hd(ej)为： 式中，N为所选窗函数w(n)的长度。经过加矩形窗后所得的滤波器实际频率响应能否很好地逼近理想频率响应呢？图1给出了理想滤波器加矩形窗后的情况。理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图。时间域内的乘积要求实际频率响应为这两个频率响应函数在频域内的卷积（卷积定理），即得到图形为图中右边的结果。图1 FIR滤波器理想与实际频率响应用窗函数法设计的滤波器性能取决于窗函数w(

6、n)的类型及窗口长度N的取值。因此，这种方法的重点在于选择某种恰当的窗函数和一种合适的理想滤波器。3 建立模型概述3.1 MATLB常用函数3.1.1 窗函数矩形窗：w=boxcar(N) 海明（Hamming）窗：w=hamming(N)布拉克曼（Blackman）窗：w=balckman(N)凯泽（Kaiser）窗：w=kaiser(N)其中是第一类变形修正零阶贝塞尔函数这些窗函数的基本参数如表1所示：表1 窗函数基本参数名称旁瓣峰值/dB近似过渡带宽精确过渡带宽最小阻带衰减/dB矩形窗-134/N1.8/N-21海明窗-418/N6.6/N-53布拉克曼窗-5712/N11/N-74凯泽

7、窗(=7.865)-5710/N-803.1.2 fir1函数设计标准响应FIR滤波器可使用firl函数。fir1函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通，带通，高通和带阻滤波器。形式为：b=fir1 (n,Wc,ftype,Window)各个参数的含义如下：b滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：n滤波器阶数；Wc截止频率，0Wc1，Wc=1对应于采样频率的一半。当设计带通滤波器时，Wc=Wc1 Wc2,Wc1Wc2；ftype当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=s

8、top时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数；Window窗函数。窗函数的长度应等于FIR滤波器系数个数，即n+1。3.1.3 freqz函数 该函数基于FFT算法计算数字滤波器Z变换频率响应。形式为 h , w = freqz ( b , a , n )返回数字滤波器的n点复频响应在简单形式中，b，a为滤波器系数，freqz可得到数字滤波器的n点复频响应，并将这n点保存在w中，相应的频率记录在h中。3.1.4 其他函数与命令设计所用其他函数及命令如表2所示表2 函数命令及其功能名称功能名称功能Clear从内存中清除变量和函数Close关闭图形Min取最小值Cei

9、l取整Angle相位角Unwrap相位角展开Figure建立图形窗口Subplot在标定位置上建立坐标系Stem离散序列图Plot线性绘图XlabelX轴标记YlabelY轴标记Title图形标题Axis控制坐标系的刻度和形式Grid网格线3.2 程序流程图Clear语句清除存储空间的变量,以免对下面的程序运行产生影响，Close all语句关闭所有已经打开的图像窗口使用fir1函数分别计算具有标准通带滤波器特性的四种不同窗的单位脉冲响应hn1,hn2,hn3,hn4；使用freqz函数计算频率响应h1,h2,h3,h4参数输入：wls,wlp,whp,whs计算过渡带宽delta_w,由于有

10、两个过渡带，使用min函数取小者计算截止频率wc1,wc2,其值取通带阻带边界频率的平均值选择矩形窗，精确过渡带宽1.8/N,使用ceil函数求满足条件的最小窗宽N1选择Blackman窗，精确过渡带宽11/N,使用ceil函数求满足条件的最小窗宽N3选择Hamming窗，精确过渡带宽6.6/N,使用ceil函数求满足条件的最小窗宽N2选择Kaiser窗(=7.865)，精确过渡带宽10/N,使用ceil函数求满足条件的最小窗宽N4取定N=100绘制矩形窗的单位脉冲响应stem (n, hn1); 幅频响应（频率归一化，幅值化分贝）plot(w1/pi,20*log10(abs(h1)； 相频

11、响应（频率归一化，angle函数计算相位，用unwrap函数按相位角展开）plot(w1/pi,180/pi*unwrap(angle(h1)绘制Blackman窗的单位脉冲响应stem(n, hn3); 幅频响应plot(w3/pi,20*log10(abs(h3) 相频响应plot(w3/pi,180/pi*unwrap(angle(h3)绘Kaiser窗的单位脉冲响应stem(n,hn4); 幅频响应plot(w4/pi,20*log10(abs(h4) 相频响应plot(w4/pi,180/pi*unwrap(angle(h4)绘制Hamming窗的单位脉冲响应stem(n, hn2)

12、; 幅频响应plot(w2/pi,20*log10(abs(h2) 相频响应plot(w2/pi,180/pi*unwrap(angle(h2)使用plot函数在同一幅图上绘制四种窗的幅频响应以便加以比较结 束使用fir1函数以及freqz函数计算四种窗频率响应h1,h2,h3,h4开 始4 模块功能分析clear; %清除掉工作空间中原有的变量，避免残留数据对本程序造成影响close all; %关闭所有已经打开的图像窗口 wls=0.2*pi;wlp=0.35*pi; %参数设置whp=0.65*pi;whs=0.8*pi;delta_w= min(wlp-wls),(whs-whp);

13、% 求两个过渡带中的小者wc1 = (wls+wlp)/2; wc2 = (whp+whs)/2; % 截止频率取通带阻带边界频率的平均值%由于N未给定，故根据各种不同窗满足技术参数所需要的最小窗宽作图%矩形窗N1=ceil(1.8*pi/delta_w); %根据矩形窗精确过渡带宽1.8/N计算窗宽hn1=fir1(N1-1,wc1,wc2/pi,boxcar(N1); % 检验设计出的滤波器单位脉冲响应h1,w1=freqz(hn1,1); % 检验设计出的滤波器的频率响应 %Hamming窗N2=ceil(6.6*pi/delta_w); %根据Hamming窗精确过渡带宽6.6/N计算

14、窗宽hn2=fir1(N2-1,wc1,wc2/pi,hamming(N2); h2,w2=freqz(hn2,1); %Blackman窗N3=ceil(11*pi/delta_w); %根据Blackman窗精确过渡带宽11/N计算窗宽hn3=fir1(N3-1,wc1,wc2/pi,blackman(N3); h3,w3=freqz(hn3,1); %Kaiser窗（=7.865)N4=ceil(10*pi/delta_w); %根据Kaiser窗技术精确过渡带宽10/N窗宽hn4=fir1(N4-1,wc1,wc2/pi,kaiser(N4,7.865); h4,w4=freqz(hn

15、4,1);%绘图figure(1) %建立图形窗口subplot(2,3,1); %分割窗口成2行3列n=0:N1-1; stem(n,hn1,.); %绘制矩形窗的单位脉冲响应axis(0 N1-1 -0.4 0.4); %设置显示范围xlabel(n);ylabel(h(n);grid on; %确定x，y轴坐标名称，加网格title(矩形窗单位脉冲响应h(n); %添加图形标题subplot(2,3,2); plot(w1/pi,20*log10(abs(h1); %绘制矩形窗的幅频特性axis(0 1 -150 5);xlabel(归一化角频率); ylabel(幅度（单位：分贝）);

16、grid on;title(矩形窗幅频响应);subplot(2,3,3); plot(w1/pi,180/pi*unwrap(angle(h1); %绘制矩形窗的相频特性xlabel(归一化角频率); ylabel(相位（单位：度）);grid on;title(矩形窗相频响应);subplot(2,3,4);n=0:N2-1; stem(n,hn2,.); axis(0 N2-1 -0.4 0.4);xlabel(n);ylabel(h(n);grid on;title(Hamming窗单位脉冲响应h(n);subplot(2,3,5); plot(w2/pi,20*log10(abs(h

17、2);axis(0 1 -150 5);xlabel(归一化角频率); ylabel(幅度（单位：分贝）);grid on;title(Hanmming窗幅频响应);subplot(2,3,6); plot(w2/pi,180/pi*unwrap(angle(h2);xlabel(归一化角频率); ylabel(相位（单位：度）);grid on;title(Hamming窗相频响应);figure(2)subplot(2,3,1);n=0:N3-1; stem(n,hn3,.); axis(0 N3-1 -0.4 0.4);xlabel(n);ylabel(h(n);grid on;titl

18、e(Blackman窗单位脉冲响应h(n);subplot(2,3,2); plot(w3/pi,20*log10(abs(h3);axis(0 1 -150 5);xlabel(归一化角频率); ylabel(幅度（单位：分贝）);grid on;title(Blackman窗幅频响应);subplot(2,3,3); plot(w3/pi,180/pi*unwrap(angle(h3);xlabel(归一化角频率); ylabel(相位（单位：度）);grid on;title(Blackman窗相频响应);subplot(2,3,4);n=0:N4-1; stem(n,hn4,.); a

19、xis(0 N4-1 -0.4 0.5);xlabel(n);ylabel(h(n);grid on;title(Kaiser窗单位脉冲响应h(n);subplot(2,3,5); plot(w4/pi,20*log10(abs(h4);axis(0 1 -150 5);xlabel(归一化角频率); ylabel(幅度（单位：分贝）);grid on;title(Kaiser窗幅频响应);subplot(2,3,6); plot(w4/pi,180/pi*unwrap(angle(h4);xlabel(归一化角频率); ylabel(相位（单位：度）);grid on;title(Kaise

20、r窗相频响应);%假定当N=100时，各种窗对比%矩形窗N=150;w1=boxcar(N);hn1=fir1(N-1,wc1,wc2/pi,boxcar(N); % 检验设计出的滤波器单位脉冲响应h1,w=freqz(hn1,1); % 检验设计出的滤波器的频率响应 %Hamming窗w2=hamming(N);hn2=fir1(N-1,wc1,wc2/pi,hamming(N); h2,w=freqz(hn2,1); %Blackman窗w3=blackman(N);hn3=fir1(N-1,wc1,wc2/pi,blackman(N); h3,w=freqz(hn3,1); %Kaise

21、r窗(=7.865)w4=kaiser(N,7.865);hn4=fir1(N-1,wc1,wc2/pi,kaiser(N,7.865); h4,w=freqz(hn4,1);figure(3)n=0:N-1;plot(n,w1,:b,n,w2,-.k,n,w3,-r,n,w4,-k); %四种窗函数在同一幅图上显示axis(-10 160 0 1.1);xlabel(n); ylabel(w(n);title(各种窗函数(N=100);legend(矩形窗带通滤波器,Hamming窗带通滤波器,Blackman窗带通滤波器,Kaiser窗带通滤波器);figure(4)plot(w/pi,2

22、0*log10(abs(h1),:b,w/pi,20*log10(abs(h2),-.g,w/pi,20*log10(abs(h3),-r,w/pi,20*log10(abs(h4),-k);axis(0 1 -150 5);xlabel(归一化角频率); ylabel(幅度（单位：分贝）);grid on;title(各种窗函数幅频特性比较(N=100);legend(矩形窗带通滤波器,Hamming窗带通滤波器,Blackman窗带通滤波器,Kaiser窗带通滤波器);5 调试过程及结论最初的学习使用Matlab阶段，由于操作不熟练，经常出现函数或者命令输入错误，中英文标点输入没有区分，漏

23、掉半个引号或分号等情况。虽然都是小错，但极其容易忽视。通常需要对照错误提示查找半天。而熟练后，因为只需调用少量几个函数即可实现设计功能，因此后面的调试过程基本不存在问题。设计结果中，矩形窗和Hamming窗情况如图2所示，Blackman窗和Kaiser窗情况如图3所示，当N=100时，四种窗形状如图4所示，其频谱特性比较如图5所示。根据理论，幅频响应图中，在3dB的对应点应该为通带截止频率，在60dB处对应阻带截止频率。图2 矩形窗和Hamming窗 由图2可以看出：对于矩形窗：窗宽N=12，h(n)为偶对称，对称中心n=5.5，由于n为整数，故在n=5和n=6位置存在两个极大值；在幅频响应

24、图中，实际设计的低端，高端通带截止频率为0.3262pi和0.6758pi，与技术指标相差6.8%和3.96%。而低端和高端的阻带截止频率为0.167pi和0.833pi，与技术指标相差16.5%和4.13%。其阻带波纹较大，第一阻带最小衰减27dB；由相频特性可以看出，此滤波器是线性相位的。对于Hamming窗：窗宽N=44，h(n)为偶对称，对称中心n=21.5，由于n为整数，故在n=21和n=22位置存在两个极大值；在幅频响应图中，实际设计的低端，高端通带截止频率为0.29pi和0.7051pi，与技术指标相差17.14%和8.47%。而低端和高端的阻带截止频率为0.1953pi和0.8

25、047pi，与技术指标相差2.35%和0.59%。第一阻带最小衰减50dB；由相频特性可以看出，此滤波器是线性相位的。图3 Blackman窗和Kaiser窗 由图3可以看出：对于Blackman窗：窗宽N=80，h(n)为偶对称，对称中心n=39.5，由于n为整数，故在n=39和n=40位置存在两个极大值；在幅频响应图中，实际设计的低端，高端通带截止频率为0.2891pi和0.7109pi，与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.207pi和0.793pi，与技术指标相差3.5%和0.875%。第一阻带最小衰减75dB；由相频特性可以看出，此滤波器是线性相位的。

26、对于Kaiser窗：窗宽N=67，h(n)为偶对称，对称中心n=33，由于n为整数，故在n=33存在一个极大值；在幅频响应图中，实际设计的低端，高端通带截止频率为0.2891pi和0.7109pi，与技术指标相差17.4%和9.37%。而低端和高端的阻带截止频率为0.2051pi和0.7969pi，与技术指标相差2.55%和0.39%。第一阻带最小衰减80dB；由相频特性可以看出，此滤波器是线性相位的。图4 四种窗函数图5 四种窗幅频特性比较 由图4，图5可以看出，当N一定时，Blackman窗的窗形状最窄，阻带衰减效果最好，但代价是主瓣宽度最大；而矩形窗的窗形状最宽，阻带衰减效果最差，但其主

27、瓣宽度最小。综合上述分析，可得：并没有哪种窗是最好的，往往某种窗在主瓣及过渡带宽方面良好，而在阻带衰减上表现不佳，反之亦如此。因此只有根据具体条件和实际需求选取最合适的。由于总的来说，滤波器主要还是强调滤波效果，即阻带衰减，因此使用Blackman或Kaiser窗效果较好。6 心得体会 这已经是我们的第三次课程设计了，虽然在闷热的夏天心里感觉很是烦躁，但相比之下，感觉没那么困惑，有想法思路计划。这次带通滤波器的设计，从原理上来说，难度并不大。只要把基本参数如通带截止频率，阻带截止频率，通带波动，阻带衰减等掌握其含义，应该不难明白滤波器的原理。 碰到的难点其一就在于Matlab的使用。本来对Ma

28、tlab的基本操作就不熟练，加上这里要用到一些信号处理的专用函数，基本上是一窍不通，导致了明白原理难下手的情况。当然，这方面的书籍比较多，在查阅资料后，就本次设计而言已经不是问题了。而且可以说这是第一次比较系统的学习了Matlab的使用，就Matlab在信息领域的广泛应用而言，对我以后不论工作或是深造，都打下了良好的基础。总的来说，每经过一次课程设计，都会对设计所用到的只是有一个真正深入的了解与应用，这样训练得到的，才是经过消化吸收的精华所在，才能在以后工作实践中灵活运用的。当然也给了和同学们相互沟通讨论的学习机会，让我从中受益匪浅。7 参考文献1 程佩青.数字信号处理教程（第3版）.清华大学出版社，20072 陈亚勇等.MATLAB信号处理详解.人民邮电出版社，20013 万永革.数字信号处理的MATLAB实现.科学出版社，20074 王力宁.MATLAB与通信仿真.人民邮电出版社，1999- 17 -