带通滤波器的设计.doc

陕西理工学院课程设计带通滤波器设计作者：汤美玲陕西理工学院（物电学院）电子信息科学与技术专业2008级 陕西 汉中 723000指导教师：蒋媛摘要：带通滤波器（bandpass filter）是从滤波器的特性上划分的，带通滤波器是指能通过某一频率范围内的频率分量、但将其他范围的频率分量衰减到极低水平的滤波器，与带阻滤波器的概念相对。从实现的网络结构或者从单位脉冲响应长度分类，可以分为无限长单位脉冲响应（IIR）滤波器和有限长单位脉冲响应（FIR）滤波器。IIR数字滤波器的设计方法是利用模拟滤波器成熟的理论及设计图表进行设计的，因而保留了一些经典模拟滤波器优良的幅度特性。但设计中只考虑了幅度特性，没考虑相位特性，所设计的滤波器一般是某种确定的非线性相位特性。为了得到线性相位特性，对IIR滤波器必须另外增加相位相位校正网络，是滤波器设计变得复杂，成本也高，又难以得到严格的线性相位特性。FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。两者各有优点，择其而取之。后面的FIR滤波器的设计中，为获得有限长单位取样响应，需要用窗函数截断无限长单位取样响应序列。另外，在功率谱估计中也要遇到窗函数加权问题。由此可见，窗函数加权技术在数字信号处理中的重要地位。关键词：带通滤波器,模拟，数字,IIR,FIR,MATLAB软件Abstract：Bandpass filter (bandpass filter) from the characteristics of the classification of the filter, belt filter is to point to by a frequency can within the scope of the frequency component, but will other range of frequency components to a very low level of attenuation filter, belt and the concept of elimination filter relative. From the network structure or realize from the unit impulse response length classification, can be divided into an infinite long unit impulse response (IIR) filter and limited long unit impulse response (FIR filter. IIR the design of the digital filter method is to use the filter mature theory and simulation design charts for design, so keep some classic simulation filter excellent range characteristics. But design only considered the range characteristics, didnt consider phase characteristic, the design is a certain general filter nonlinear phase characteristic. In order to get the linear phase characteristic, for an additional filter must IIR phase phase correction network, is filter design complicated, the cost is high, and hard to get the strict linear phase characteristic. FIR filter in the guarantee range characteristics to meet technical requirements at the same time, very easy to do have the strict linear phase characteristic. Both have their advantages, pick the and of the take. The back of the FIR filters design, to acquire limited long unit sampling response, need to use the window function truncation infinite long unit sampling response sequence. In addition, in the power spectrum estimation to meet a window function and weighted problem. This shows, window function weighted technology in the digital signal processing to the important position. Key words：Bandpass filter, simulation , digital , IIR , FIR , MATLAB software 一.任务 1基于IIR模拟带通滤波器的设计.2基于IIR数字带通滤波器的设计. 3基于窗函数的FIR带通滤波器的设计二.要求 1.基本要求 1.1可显示任何汉字字符. 1.2 可实现花样显示. 2.发挥部分 2.1 不需要使用专门的字模软件提取固定汉字字模. 2.2 可人性化设置.三.说明3.1 时间要求：11月12日到11月24日.3.2 完成matlab设计程序、仿真,总结报告.四. 带通滤波器的设计原理、指标及方法步骤1. 带通滤波器的设计原理及窗函数法的MATLAB设计函数简介1.1. 带通滤波器的设计原理一个理想的滤波器应该有一个完全平坦的通带，例如在通带内没有增益或者衰减，并且在通带之外所有频率都被完全衰减掉，另外，通带外的转换在极小的频率范围完成。实际上，并不存在理想的带通滤波器。滤波器并不能够将期望频率范围外的所有频率完全衰减掉，尤其是在所要的通带外还有一个被衰减但是没有被隔离的范围。这通常称为滤波器的滚降现象，并且使用每十倍频的衰减幅度dB来表示。通常，滤波器的设计尽量保证滚降范围越窄越好，这样滤波器的性能就与设计更加接近。然而，随着滚降范围越来越小，通带就变得不再平坦开始出现“波纹”。这种现象在通带的边缘处尤其明显，这种效应称为吉布斯现象。除了电子学和信号处理领域之外，带通滤波器应用的一个例子是在大气科学领域，很常见的例子是使用带通滤波器过滤最近3到10天时间范围内的天气数据，这样在数据域中就只保留了作为扰动的气旋。在频带较低的剪切频率f1和较高的剪切频率f2之间是共振频率，这里滤波器的增益最大，滤波器的带宽就是f2和f1之间的差值。下图1.1为理想模拟带通滤波器幅频特性： 图4.1 带通滤波器的幅频特性曲线1.2. 窗函数法的MATLAB设计函数简介（1）窗函数法设计原理 设希望逼近的滤波器频率响应函数为，其单位脉冲响应是。如果能够由已知的求出，经过变换可得到滤波器的系统函数。但通常以理想滤波器作为，其幅度特性逐段恒定，在边界频率处有不连续点，因而是无限时宽的，且是非因果序列。为了构造一个长度为N的第一类线性相位FIR滤波器，只有将截取一段，并保证截取的一段关于偶对称。设截取的一段用表示，即式中，是一个矩形序列，长度为N，当取值为时，截取的一段关于偶对称，保证所涉及的滤波器具有线性相位。我们实际所设计的滤波器的单位脉冲响应为，长度为N，其系统函数为，即这样用一个有限长的序列去代替，肯定会引起误差，表现在频域就是通常所说的吉布斯效应。该效应引起过渡带加宽及通带和阻带内的波动，尤其使阻带的衰减减小，从而不满足技术上的要求，这种吉布斯效应是由于将直接截断引起的，也称为截断效应。通过构造窗函数，用来减少截断效应，从而设计一个能满足技术要求的FIR线性相位滤波器。以上即为窗函数法设计FIR滤波器的思想。总结为以下步骤：通过傅里叶逆变换获得理想滤波器的单位脉冲响应。从时域出发，截取有限长的一段冲击响应作为的系数，冲击响应长度N就是系统函数的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。一般这种时域设计、频域检验的方法要反复几个回合才能成功。要设计一个线性相位的FIR数字滤波器，首先要求理想频率响应。是的周期函数，周期为，可以展开成傅氏级数： 其中是与理想频响对应的理想单位抽样响应序列。但不能用来作为设计FIR DF用的，因为一般都是无限长、非因果的，物理上无法实现。为了设计出频响类似于理想频响的滤波器，可以考虑用来近似。由性能指标确定窗函数和窗口长度N。设是一个长序列，是长度为N的窗函数，用截断，得到N点序列，在频域上则有求得实际滤波器的单位脉冲响应，即为所设计FIR滤波器系数向量。 由此可见，窗函数不仅仅会影响原信号在时域上的波形，而且也会影响到频域内的形状。（2）MATLAB窗函数MATLAB软件工具箱主要提供了以下五种窗函数，如表4-1所示为五种窗函数的定义式，表4-2为五种窗函数的性能。表4-1 MATLAB窗函数窗 函 数 定义式矩形窗（Boxcar）三角窗(Triang)海明窗(Hamming) 0nM-1汉宁窗(Hanning) 0nM-1巴特利特窗（Bartlett）表4-2 MATLAB窗函数窗类型旁瓣峰值主瓣峰值最小阻带衰减矩形窗13dB4/M21dB三角窗25dB8/M25dB汉宁窗31dB8/M44dB海明窗41dB8/M53dB凯塞窗57dB12/M74dB（3）窗函数法的MATLAB设计函数简介实际设计一般用MATLAB工具箱函数。可调用工具箱函数fir1实现窗函数法设计步骤。firl是用窗函数法设计线性相位FIR数字滤波器的工具箱函数，实现线性相位FIR数字滤波器的标准窗函数法设计。fir2为任意形状幅度特性的窗函数法设计函数，用fir2设计时，可以指定任意形状的，它实质是一种频率采样法与窗函数法的综合设计函数。主要用于设计幅度特性形状特殊的滤波器。用help命令查阅其调用格式及调用参数的含义。2. 带通滤波器的设计指标逼近通带：,通带最大衰减：dB；逼近阻带：0,阻带最小衰减： dB。remezord调用参数： 。3. 带通滤波器的设计方法及步骤3.1. 模拟带通滤波器的设计步骤（1）通过映射关系式，将希望设计的带通滤波器指标转换为相应的低通圆形滤波器的指标。（2）设计相应的归一化低通系统函数。（3）用式将转换成所希望设计得到的带通滤波器系统函数。3.2. 数字带通滤波器的设计步骤（1）确定数字带通滤波器的技术指标：通带边界频率、通带最大衰减、阻带截止频率、阻带最小衰减。（2）将数字带通滤波器的技术指标转换成相应的模拟带通滤波器的技术指标。这里主要是边界频率和的转换，和指标不变。如果采用脉冲响应不变法，边界频率的转换关系为如果采用双线性变换法，边界频率的转换关系为（3）按照模拟带通滤波器的技术指标设计过渡模拟带通滤波器。（4）用所选的转换方法，将模拟滤波器转换成数字带通滤波器系统函数。3.3. 窗函数法设计FIR带通滤波器的设计步骤（1）根据对阻带衰减及过渡带的指标要求，选择窗函数的类型，并估计窗口长度N。（2）构造希望逼近的频率响应函数，即所谓的“标准窗函数法”，就是选择为线性相位理想带通滤波器。（3）计算。如果给出待求滤波器的频率响应函数为，那么单位脉冲响应可用下式求出：（4）加窗得到设计结果：。五.带通滤波器的设计1. IIR和FIR滤波器的比较从性能上说，IIR滤波器系统函数的极点可位于单位圆内的任何地方，因此零点和极点相结合，可用较低的阶数获得较高的选择性，所用的存储单元少，计算量小，所以经济高效。但这个高效是以相位的非线性为代价的。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器系统函数的极点固定在原点，因而只能用较高的阶数达到高的选择性；对于同样的滤波器幅频特性指标,FIR滤波器所要求的阶数一般比IIR滤波器高510倍，成本较高，信号延时也较大。下图5.1为IIR数字滤波器的设计方法分类框图：IIR数字滤波器的设计方法 直接设计IIR数字滤波器从模拟滤波器设计IIR滤波器时域逼近法频域逼近法双线性不变法脉冲响应不变法零极点累试法图5.1 IIR数字滤波器的设计方法分类从结构上看，IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。另外，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引起寄生振荡。相反，FIR滤波器主要采用非递归结构，不论在理论上还是实际的有限精度运算中都不存在稳定性问题，运算误差引起的输出信号噪声功率也较小。FIR滤波器可采用FFT算法实现，在相同阶数条件下，运算速度大大提高。从设计工具看，IIR滤波器可以借助成熟模拟滤波器设计成果，因此一般都有封闭形式的设计公式可供准确计算，计算量小，对工具要求低。FIR滤波器计算通带和阻带衰减等仍无显示表达式，其边界频率也不易精确控制。FIR滤波器设计只有计算程序可循，因此对计算工具要求高。IIR滤波器虽然设计简单，但主要适用于设计具有片段常数特性的选频型滤波器，如低通、高通、带通、带阻等，往往脱离不了几种典型模拟滤波器的频响特性的约束。而FIR滤波器 则要灵活的多，易于适应某些特殊的应用，两者各有所长，所以在实际应用时应该全面考虑加以选择。2. 基于IIR模拟带通滤波器的设计（1）利用MATLAB设计巴特沃斯带通滤波器：利用MATLAB函数buttap(语法Z,P,K=buttap(N)，Z,P,K分别表示滤波器的零点，极点，增益)设计，验证极点位置是否正确。技术指标为：通带截止频率，通带的最大衰减，阻带截止频率阻带的最小衰减。MATLAB程序如下：a=10000;b=40000;Ap=3;As=35; %技术指标T1=sqrt(10(0.1*P)-1);T2=sqrt(10(0.1*S)-1);N=ceil(log10(T2/T1)/log10(b/a); %求阶数Z,P,K=buttap(N); %求零、极点，增益syms rad; %定义基本符号对象hs1=K/(i*rad/a-P(1)/(i*rad/a-P(2)/(i*rad/a-P(3);hs2=10*log10(abs(hs1)2);ezplot(hs2,-60000,60000); %绘制符号表达式的二维曲线grid on; %画坐标网格线仿真图形如图5.2所示： 图5.2 巴特沃斯带通滤波器频率响应曲线（2）利用MATLAB设计切比雪夫带通滤波器：现在以MATLAB的函数cheblap(语法：z,p,k=cheblap(n，rp)，其中n为阶数，rp为通带的幅度差，z,p,k分别代表零点、极点、增益)为工具设计切比雪夫滤波器。设计技术指标为：通带截止频率，通带的最大衰减，阻带截止频率，阻带的最小衰减。程序如下：a=10000;b=40000;Ap=3;As=35; %技术指标T1=sqrt(10(0.1*AP)-1);T2=sqrt(10(0.1*AS)-1);N=ceil(acosh(T2/T1)/acosh(b/a); %求阶数Z,P,K=cheb1ap(N,AP); %求零、极点，增益syms rad; %定义基本符号对象hs1=K/(i*rad-P(1)/(i*rad-P(2)/(i*rad-P(3);hs2=10*log10(abs(hs1)2);ezplot(hs2,-15,15); %绘制符号表达式的二维曲线grid on; %画坐标网格线仿真图形如图5.3所示： 图5.3 切比雪夫型带通滤波器响应曲线（3）利用MATLAB设计椭圆带通滤波器：用ellipap函数设计 (语法：z,p,k=ellipap(n，rp，rs)，其中n为阶数，rp为通带的幅度差，rs为阻带内的最小衰减，z,p,k分别代表零点，极点，增益)。技术指标为：通带截止频率，通带的最大衰减，阻带截止频率阻带的最小衰减。程序如下：pd=1;sd=4;pn=1;sn=35; %技术指标rad0=sqrt(pd*sd); %求归一化基准频率k=pd/sd; %求选择因数nd=pd/rad0; %求通带截止频率归一化值nad=sd/rad0; %求阻带截止频率归一化值q0=(1-(1-k*k)0.25)/(1+(1-k*k)0.25)/2;q=q0+2*q05+15*q09+15*q013;b=10(0.1*sn);e=10(0.1*pn);n=ceil(log10(16*(b-1)/(e-1)/log10(1/q); %求阶数z,p,k=ellipap(n,pn,sn); %求零、极点，增益syms rad; %定义基本符号对象hs1=k*(i*rad/rad0-z(1)*(i*rad/rad0-z(2)/(i*rad/rad0-p(1)/(i*rad/rad0-p(2)/(i*rad/rad0-p(3);hs2=10*log10(abs(hs1)2);ezplot(hs2,-10,10); %绘制符号表达式的二维曲线grid on; %画坐标网格线仿真图形如图5.4所示：图5.4 椭圆带通滤波器频率响应曲线从上面的IIR带通滤波器设计特性曲线图可以看出巴特沃斯滤波器的频率特性在通带和阻带内部都是随着频率的变化而单调下降的，在通带内部能够完成设计出超过预定指标的要求，但会造成滤波器的阶数N比较高；切比雪夫滤波器在通带范围内具有良好的等波纹特性，能将指标的精度均匀分布在整个通带内，或者均匀分布在阻带内，但在通带范围外的幅频特性是单调递减的；椭圆滤波器的系统函数在截止频率的附近含有有限个零点和极点，能够更好的逼近理想低通滤波器的特性。在三种常用滤波器的设计方法中，选择哪种形式的滤波器可以根据具体应用情况而定。通常椭圆滤波器的阶数低，切比雪夫滤波器次之，巴特沃斯滤波器最高，而参数量化的灵敏度与此相反，巴特沃斯滤波器最佳，切比雪夫滤波器次之，椭圆滤波器对量化最为敏感。2. 基于IIR数字带通滤波器的设计用MATLAB设计一个IIR数字带通滤波器：要求：通带中心频率为，通带截止频率为，通带最大衰减为；阻带最小衰减为，阻带截止频率。设计步骤：根据要求，采用切比雪夫型数字滤波器涉及该滤波器；设采样频率Fs=2000Hz,则可计算出fp0=500Hz, fp1=400Hz, fp2=600Hz, fs1=300 Hz, fs2=700Hz；给出一测试信号，测试该带通滤波器的功能。本例调用MATLAB函数直接设计数字带通滤波器。设计程序ep2.m如下：%设计数字带通滤波器程序：ep2.mfunction y=bandp(x,fp1,fp2,fs1,fs2,rp,rs,Fs) %带通滤波%使用注意事项：通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半，%即，fp1,fp2,fs1,fs2的值小于Fs/2%x：需要带通滤波的序列%fp1：通带左边界%fp2：通带右边界%fs1：衰减截止左边界%fs2：衰减截止右边界%rp：边带区衰减DB数设置%rs：截止区衰减数DB设置%Fs：序列x的采样频率fp1-400;fp2-600;fs1-300;fs2-700;rp=3;rs=15;Fs=2000;wp1=0.4*pi;wp2=0.6*pi;ws1=0.3*pi;ws2=0.7*pi;wp=wp1 wp3;ws=ws1 ws2;%设计切比雪夫滤波器;n,wn=cheblord(ws/pi,rp,rs);bz1,az1=cheby1(n,rp,wp/pi);%查看滤波器的曲线h,w=freqz(bz1,az1,256,Fs);h=20*log10(abs(h);figure;plot(w,h);title(所设计滤波器的通带曲线);grid on;y=filter(hz1,az1,x);end %带通滤波器测试程序Fs=2000;t=(1:Fs)/Fs;ff1=100;ff2=400;ff3=700;x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t)+sin(2*pi*ff3*t);figure;subplot(211);plot(t,x);subplot(212);hua_fft(y,Fs,1); function hua_fft(y,style,varargin)%当style=1，画幅值谱;当style=2，画功率谱；当style=其他的，那么画幅值谱和功率谱%当style=1时，还可以多输入2个可选参数%可选输入参数是用来控制需要查看的频率段的%第一个是需要查看的频率段起点%第一个是需要查看的频率段的终点%其他style不具备可选输入参数，如果输入发生位置错误Nfft=2nextpow2(length(y); %找出大于y的个数的最大的2的指数值（自动运算最佳FFT步长nfft）%nfft=1024;%人为设置FFT的步长nffty=y-mean(y);%去除直流分量y_ft=fft(y,nfft);%对y信号进行DFT，得到频率的幅值分量y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是它本身。y_f=fs*(0:nfft/2-1)/nfft;if style=1if nargin=3plot(y_f,2*abs(y_ft(1:nfft/2)/length(y);%matlab的帮助里画FFT的方法 %ylabel(幅值);xlabel(频率);title(信号幅值谱);%plot(y_f,abs(y_ft(1:nfft/2);%论坛上画FFT的方法else f1=varargin1; fn=varargin2; ni=round(f1*nfft/fs+1); na=round(fn*nfft/fs+1); plot(y_f(ni:na),abs(y_ft(ni:na)*2/nfft);endelseif style=2 plot(y_f,y_p(1:nfft/2); %ylabel(功率谱密度);xlabel(频率);title(信号功率谱);else Subplot(211);plot(y_ft(1:nfft/2)/length(y); ylabel(幅值);xlabel(频率);title(信号幅值谱); Subplot(212);plot(y_f，y_p(1:nfft/2); ylabel(功率谱密度);xlabel(频率);title(信号功率谱);endend仿真结果如下图5.3为该带通滤波器频谱图形，图5.4为经过带通滤波后的信号时及频谱图形，图5.5原始输入信号时域及频谱波形。 图5.3 带通滤波器频谱图图5.4 经带通滤波后的信号时及频谱图图5.5 原始输入信号时域及频谱波形经图示可看出该滤波器具有带通特性，它将通过它的信号在带通范围内的频率成分保留，阻带范围内的频率成分截止，从而达到滤波作用。3. 基于窗函数的FIR带通滤波器的设计（1）利用矩形窗进行设计程序如下：fs=20000; %设定采样频率fp1=4500;fp2=6500; %第一截止频率fs1=3000;fs2=7500; %第二截止频率As=40; %最小阻带衰减Ws1=(fp1+fs1)/fs;Ws2=(fp2+fs2)/fs; %截止频率归一化处理w=(fp1-fs1)/fs; %求归一化过渡带M=ceil(As-7.95)/(14.36*w) %计算所需滤波器的阶数juxing=boxcar(M+1); %生成长度为M+1的矩形窗boxb=fir1(M,Ws1,Ws2,juxing); %生成矩形窗设计的fir滤波器freqz(boxb,1,fs,fs); %绘制幅频和相频响应曲线运行结果： 图5.6 矩形窗fir滤波器幅频和相频响应曲线从幅频响应上看，通带基本无波纹，阻带中波纹较大，因而阻带较不理想，相频响应曲线在通带内为直线，效果较好，信号失真小。 （2）利用三角窗进行设计利用三角窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用triang（）函数生成三角窗，程序运行结果如下： 图5.7 三角窗fir滤波器幅频和相频响应曲线（3）利用凯塞窗进行设计程序如下：fs=20000; %设定采样频率fp1=4500;fp2=6500; %第一截止频率fs1=3000;fs2=7500; %第二截止频率n,wn,bta,ftype=kaiserord(fs1,fp1,fp2,fs2;0,1,0,0.01 0.1087 0.01,fs) %求滤波器参数b=fir1(n,wn,ftype,kaiser(n+1,bta); %生成fir滤波器freqz(b,1,fs,fs); %绘制幅频和相频响应曲线程序运行结果：图5.8 凯塞窗fir滤波器幅频和相频响应曲线（4）利用汉宁窗进行设计利用汉宁窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用hanning（）函数生成三角窗，程序运行结果如下： 图5.9 汉宁窗fir滤波器幅频和相频响应曲线（5）利用海明窗进行设计利用海明窗进行设计时，原理与矩形窗基本相同，只不过生成窗函数时采用hamming（）函数生成三角窗，程序运行结果如下：图5.10 海明窗fir滤波器幅频和相频响应曲线（6）滤波器的性能测试：由MATLAB模拟生成含有不同频率的数字信号，然后利用设计的滤波器对数字信号进行滤波，为方便观察，模拟生成的信号只含有包含在阻带的两个频率（2000Hz，8000Hz）一个包含于通带的频率（5000Hz）。程序如下：fs=20000;t=0:1/fs:2; x=sin(2*pi*2000*t)+sin(2*pi*5000*t).+sin(2*pi*8000*t); %生成混合信号xo=filter(b,2,xn); %用滤波器对信号进行滤波figure;nn=5000:5100; %取一段信号subplot(211);tt=nn/fs;plot(tt,x(nn); %绘制原始信号axis(0.25,0.255,-4,4);ylabel(原始信号);xlabel(时间);subplot(212);plot(tt,xo(nn); %绘制滤波后的信号axis(0.25,0.