FIR带通滤波器的设计成熟版.doc

目 录1 技术要求12 基本原理12.1 FIR滤波器简介12.2 窗函数法原理23 建立模型描述43.1 MATLAB常用函数43.1.1 矩形窗函数43.1.2 三角窗函数53.1.3 广义余弦窗53.1.4 汉宁窗（升余弦窗）63.1.5 fir1函数73.1.6 freqz函数73.1.7 其他函数与命令83.2 方案设计与论证83.3 程序流程图94 模块功能分析或源程序代码105 调试过程及结论135.1 实验结果135.2 结果分析156 思考题157 心得体会168 参考文献17FIR带通滤波器的设计1 技术要求用窗函数法设计FIR带通滤波器。要求低端阻带截止频率ls=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滤波器的系统差分方程为：因此，FIR滤波器又称为卷积滤波器。FIR滤波器的频率响应表达式为：信号通过FIR滤波器不失真条件是在通带内具有恒定的幅频特性和线性相位特性。理论上可以证明：当FIR滤波器的系数满足下列中心对称条件： 或者 时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个 N 阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。FIR滤波器的设计任务是选择有限长度的h(n)，使传输函数满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本次设计使用窗函数法设计FIR带通滤波器。2.2 窗函数法原理数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是和重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。因此设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲击响应作为H(z)的系数，冲击响应长度N就是系统函数H(z)的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是w的周期函数，周期为，可以展开成傅氏级数： （公式1-1）使用上述的传递函数去逼近，一个理想的频率响应的傅立叶反变换： （公式1-2）其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的h(n)，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用来近似。窗函数的基本思想：先选取一个理想滤波器（它的单位抽样响应是非因果、无限长的），再截取（或加窗）它的单位抽样响应得到线性相位因果FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器。设是一个长序列，是长度为N的窗函数，用截断，得到N点序列，即 （公式1-3）在频域上则有 （公式1-4）由此可见，窗函数不仅仅会影响原信号在时域上的波形，而且也会影响到频域内的形状。MATLAB信号工具箱主要提供了以下几种窗函数，如表1 表1 MATLAB窗函数窗窗 函 数矩形窗Boxcar三角窗Triang 海明窗Hamming汉宁窗Hanning布莱克曼Blackman切比雪夫窗Chebyshev凯塞窗Kaiser加矩形窗后的频谱和理想频谱可得到以下结论：加窗使过渡带变宽，过渡带的带宽取决于窗谱的主瓣宽度。矩形窗情况下的过渡带宽是。N越大，过渡带越窄、越陡；过渡带两旁产生肩峰，肩峰的两侧形成起伏振荡。肩峰幅度取决于窗谱主瓣和旁瓣面积之比。矩形窗情况下是8.95，与N无关。工程上习惯用相对衰耗来描述滤波器，相对衰耗定义为： （公式1-5）这样两个肩峰点的相对衰耗分别是0.74dB和-21dB。其中（-0.0895）对应的点的值定义为阻带最小衰耗。以上的分析可见，滤波器的各种重要指标都是由窗函数决定，因此改进滤波器的关键在于改进窗函数。窗函数谱的两个最重要的指标是：主瓣宽度和旁瓣峰值衰耗。旁瓣峰值衰耗定义为： 旁瓣峰值衰耗20lg(第一旁瓣峰值主瓣峰值) （公式1-6）为了改善滤波器的性能，需使窗函数谱满足：(1)主瓣尽可能窄，以使设计出来的滤波器有较陡的过渡带。(2)尽量减少最大旁瓣的相对幅度，也就是能量集中于主瓣，以减小带内、带外波动的最大幅度，增大阻带衰减。一般来说，以上两点很难同时满足。当选取主瓣宽度很窄时，旁瓣的分量势必增加，从而带内、带外的波动也增加了；当选取最小的旁瓣幅度时，降低了带内、带外的波动，但是过渡带的陡度减小了。所以实际采用的窗函数其特性往往是它们的折中，在保证主瓣宽度达到一定要求的前提下，适当牺牲主瓣宽度来换取旁瓣波动的减小3 建立模型描述3.1 MATLAB常用函数3.1.1 矩形窗函数矩形窗(Rectangular Window)函数的时域形式可以表示为： （公式2-1）它的频域特性为 （公式2-2）Boxcar函数：生成矩形窗调用方式w = boxcar (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。从功能上讲，该函数又等价于w = ones(n,1)。3.1.2 三角窗函数三角窗（Bartlett Window）函数时域形式可表示为： （公式2-3）窗谱为： （公式2-4） 式中，当N远大于1时，此时，窗谱主瓣宽度为8/N。3.1.3 广义余弦窗汉宁窗、海明窗和布莱克曼窗，都可以用一种通用的形式表示，这就是广义余弦窗。这些窗都是广义余弦窗的特例，汉宁窗又被称为余弦平方窗或升余弦窗，海明窗又被称为改进的升余弦窗，而布莱克曼窗又被称为二阶升余弦窗。采用这些窗可以有效地降低旁瓣的高度，但是同时会增加主瓣的宽度。这些窗都是频率为0、2/(N1)和4/(N1)的余弦曲线的合成，其中N为窗的长度。通常采用下面的命令来生成这些窗： （公式2-5） （公式2-6）其中，A、B、C适用于自己定义的常数。根据它们取值的不同，可以形成不同的窗函数，分别是：汉宁窗 A=0.5，B=0.5，C=0；海明窗 A=0.54，B=0.54，C=0；布莱克曼窗 A=0.5，B=0.5，C=0.08；3.1.4 汉宁窗（升余弦窗）汉宁窗（Hanning）函数时域形式可表示为： （公式2-7）利用傅利叶变换的调制特性，由上式可得汉宁窗的平谱函数为： （公式2-8）式中， (公式2-9)当N远大于1时，上式可近似表示为： （公式2-10）这三部分之和使旁瓣互相抵消，能量更集中在主瓣，汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8/N。hanning函数：生成汉宁窗调用方式：(1) w = hanning(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。注意：此函数不返回是零点的窗函数的首尾两个元素。 (2) w = hanning(n,'symmetric')：与上面相类似。(3) w = hanning(n,'periodic')：此函数返回包括为零点的窗函数的首尾两个元素。3.1.5 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=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数；Window窗函数。窗函数的长度应等于FIR滤波器系数个数，即n+1。3.1.6 freqz函数 该函数基于FFT算法计算数字滤波器Z变换频率响应。形式为 h , w = freqz ( b , a , n )返回数字滤波器的n点复频响应在简单形式中，b，a为滤波器系数，freqz可得到数字滤波器的n点复频响应，并将这n点保存在w中，相应的频率记录在h中。3.1.7 其他函数与命令设计所用其他函数及命令如表2所示表2 函数命令及其功能名称功能名称功能Clear从内存中清除变量和函数Close关闭图形Min取最小值Ceil取整Angle相位角Unwrap相位角展开Figure建立图形窗口Subplot在标定位置上建立坐标系Stem离散序列图Plot线性绘图XlabelX轴标记YlabelY轴标记Title图形标题Axis控制坐标系的刻度和形式Grid网格线3.2 方案设计与论证用窗函数法设计一个FIR带通滤波器。指示如下：低端阻带截止频率 wls = 0.2*pi；低端通带截止频率 wlp = 0.35*pi；高端通带截止频率 whp = 0.65*pi；高端阻带截止频率 whs = 0.8*pi；六种窗函数的基本参数如表3：表3 窗函数基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134/N-21三角形窗-258/N-25汉宁窗-318/N-44哈明窗-418/N-53布莱克曼窗-5712/N-74凯塞窗()-5710/N-80以上表格里的参数设置是最佳窗函数设计，根据设计方案的要求，选择一个合适的窗函数进行滤波器的设计，从上表可以看出：最小带阻衰减仅有窗函数决定，不受N的影响，而过渡带的宽度则随窗函数的增加而减小。绘制矩形窗的单位脉冲响应 幅频响应 相频响应绘制Hamming窗的单位脉冲响应 幅频响应 相频响应绘制Kaiser窗的单位脉冲响应 幅频响应 相频响应绘制Blackman窗的单位脉冲响应 幅频响应 相频响应选择Hamming窗，精确过渡带宽6.6/N选Blackman窗，精确过渡带宽11/N选择Kaiser窗，精确过渡带宽10/N选择矩形窗，精确过渡带宽1.8/N使用fir1函数计算通带滤波器特性使用freqz函数计算频率响应Clear语句清除存储空间的变量影响输入参数计算过渡带宽delta_w开始结束3.3 程序流程图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); %求两个过渡带的较小者wc1=(wls+wlp)/2;wc2=(whp+whs)/2; 截止频率取通带阻带边界频率的均值%矩形窗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计算窗宽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 窗N4=ceil(10*pi/delta_w); %根据Kaiser窗技术精确过渡带宽10/N计算窗宽hn4=fir1(N4-1,wc1,wc2/pi,kaiser(N4);h4,w4=freqz(hn4,1);%绘图figure(1) %建立图形窗口subplot(3,1,1); %把窗口分割成3行1列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(3,1,2);plot(w1/pi,20*log10(abs(h1); %绘制矩形窗的幅频特性曲线axis(0,1,-150,5); %设置显示范围xlabel('归一化角频率');ylabel('幅度（单位:分贝）');grid on; %确定x，y坐标title('矩形窗幅频响应'); %添加图形的标题subplot(3,1,3);plot(w1/pi,180/pi*unwrap(angle(h1); %绘制矩形窗的相频特性曲线xlabel('归一化角频率');ylabel('单位：度');grid on;title('矩形窗相频相应');figure(2) %建立图形窗口subplot(3,1,1);n=0:N2-1;stem(n,hn2,'.'); %绘制Hamming窗单位脉冲响应axis(0,N2-1,-0.4,0.4); %确定显示范围xlabel('n');ylabel('h(n)');grid on;title('Hamming窗单位脉冲响应h(n)');subplot(3,1,2);plot(w2/pi,20*log10(abs(h2); %绘制Hamming窗幅频响应axis(0,1,-150,5);xlabel('归一化角频率');ylabel('幅度（单位：分贝）');grid on;title('Hamming窗幅频响应');subplot(3,1,3);plot(w2/pi,180/pi*unwrap(angle(h2); %绘制Hamming窗相频响应xlabel('归一化角频率');ylabel('单位：度');grid on;title('Hamming窗相频相应');figure(3)subplot(3,1,1);n=0:N3-1;stem(n,hn3,'.');axis(0,N3-1,-0.4,0.4);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.4,0.4);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窗相频相应');5 调试过程及结论5.1 实验结果新建m文件然后运行程序，检查无误后，运行结果如下：图1 矩形窗滤波器图2 Hamming窗滤波器图3 Blackman窗滤波器图4 Kaiser窗滤波器5.2 结果分析设计结果中，矩形窗和Hamming 窗的情况结果如图1和2所示，Blackman窗和Kaiser窗的运行结果如图3和4所示。根据理论分析，幅频响应图中，在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；有相频特性可以看出，次滤波器为线性相位的。6 思考题窗函数设计法中，选择窗函数的类型与滤波器阶数对滤波器设计的影响？答：首先明确一个关系式，如果窗宽为N，则滤波器阶数为N-1，所以滤波器的阶数的增大或减小，相当于窗宽的增大或减小。FIR滤波器的设计目标就是尽可能的接近理想滤波器。加窗处理后滤波器的幅度函数为。可见只有当窗谱逼近冲激响应时，才会逼近。但这时，窗宽N趋近于无穷大，相当于没有加窗处理。因此，我们希望是：窗谱主瓣尽可能窄，以便获得较陡过渡带，而主瓣宽度只跟N呈反比，主瓣窄要求窗宽N要大；尽量减小最大旁瓣的相对幅度，以便增大阻带的衰减，而阻带衰减大则要求窗宽N要小。 7 心得体会这已经是我们第三次做课程设计了，可以说是轻车熟路了。这次带通滤波器的设计，从原理上来说，难度并不大。只要把基本参数如通带截止频率，阻带截止频率，通带波动，阻带衰减等掌握其含义，就能很好的完成这次试验。碰到的难点其一就在于Matlab的使用。最初的学习使用MATLAB软件阶段，由于操作不熟练，经常出现函数或者命令输入错误，中英文标点输入没有区分，漏掉棒引号或者分号等情况。虽然都是小错，但是极其容易被忽视而产生错误。通常按照错误要找半天。而慢慢熟练之后只需要调用少量几个函数就能实现设计功能，因此后面的调试过程基本上不存在问题。当我学会使用Matlab后，发现一直让我感到枯燥的编程原来并不困难，起码设计带通滤波器是这样，只要调用几个函数就可以轻松完成。在这之后的难点就是结论分析了。虽然作出的图形效果让我有成就感，但这并不是最终结果。一次实验结束后，通过得到的实验数据可以提出或者验证什么理论，才是实验的目的所在。而这恰恰是我目前所缺乏的能力。一直以来，我们所做的实验设计基本都是验证性的，对与不对，一看结果就知道，这直接导致了我忽视了后续的分析，认为只要结果正确就行。就拿这次来说，虽然用不同的窗实现了同一个功能，但用这么多窗的意义何在？当然是为了找到最好的。那如何找到？任务就落在了后续分析上。而经过分析发现，实际上并没有最好的，只有根据不同情况选择最适合的。尽管现在只是初步学会了简单数字滤波器的设计，离真正掌握还有一定距离，但这段日子确实令我收益匪浅，这将对我今后的学产生积极的影响。8 参考文献1 程佩青.数字信号处理教程（第3版）.清华大学出版社，20072 陈亚勇等.MATLAB信号处理详解.人民邮电出版社，20013 万永革.数字信号处理的MATLAB实现.科学出版社，20074 王力宁.MATLAB与通信仿真.人民邮电出版社，1999