数字信号处理课程论文-基于MATLAB的滤波器设计.doc

数字信号处理课程论文基于MATLAB的滤波器设计姓名： 学号：10年级：2010级 学院：物理电气信息学院 专业：通信工程指导老师： 完成日期：2112年12月9日基于MATLAB的滤波器设计摘 要：利用MATLAB 设计滤波器,可以按照设计要求非常方便地调整设计参数,极大地减轻了设计的工作量,有利于滤波器设计的最优化。Matlab因其强大的数据处理功能被广泛应用于工程计算,其丰富的工具箱为工程计算提供了便利,利用Matlab信号处理工具箱可以快速有效地设计各种数字滤波器,设计简单方便。本文介绍了在MATLAB R2011a 环境下滤波器设计的方法和步骤。关键词：滤波器，matlab，FIR，IIR Abstract：By using MATLAB , we can design filters and modify the filtersparameters conveniently according to our demands. This relieves greatly design work loads and makes for optimization of filter designing. Matlab can be widely used in engineering calculations because of its powerful functions of data processing. Its rich toolbox makes the calculations easy. With Matlab signal processing toolbox, various digital filters can be designed effectively in simple way. This article introduce the methods and processes in the circumstance of MATLAB R2011a.Keywords：filter，matlab，fdatool1.滤波器的原理凡是可以使信号中特定的频率成分通过，而极大地衰减或抑制其他频率成分的装置或系统都称之为滤波器，相当于频率“筛子”。 滤波器的功能就是允许某一部分频率的信号顺利的通过，而另外一部分频率的信号则受到较大的抑制，它实质上是一个选频电路。滤波器中，把信号能够通过的频率范围，称为通频带或通带；反之，信号受到很大衰减或完全被抑制的频率范围称为阻带；通带和阻带之间的分界频率称为截止频率；理想滤波器在通带内的电压增益为常数，在阻带内的电压增益为零；实际滤波器的通带和阻带之间存在一定频率范围的过渡带。a理想滤波器的频率特性理想滤波器:使通带内信号的幅值和相位都不失真,阻喧内的频率成分都衰减为零的滤波器，其通带和阻带之间有明显的分界线。如理想低通滤波器的频率响应函数为 或 理想滤波器实际上并不存在。b实际滤波器 实际滤波器的特性需要以下参数描述：1）恒部平均值A0：描述通带内的幅频特性；波纹幅度：d。2）上、下截止频率：以幅频特性值为A0/2时的相应频率值WC1，WC2作为带通滤波器的上、下截止频率。带宽。因为 所以 也称“-3dB”带宽3）选择性：实际滤波器过渡带幅频曲线的倾斜程度表达了滤波器对通带外频率成分的衰减能力，用信频程选择性和滤波器因素描述。信频程选择性：与上、下截止频率处相比，频率变化一倍频程时幅频特性的衰减量，即倍频程选择性：= 或 =信频程选择性总是小于等于零，显然，计算信量的衰减量越大，选择性越好。滤波器因素：60dB处的带宽与3dB处的带宽之比值，即越小,选择性越好.分辨力：即分离信号中相邻频率成分的能力，用品质因素Q描述。，Q越大，分辨率越高。2. MTALAB的在DSP应用上的功能MATLAB中的信号处理工具箱内容丰富，使用简便。在数字信号处理中常用的算法，如FFT，卷积，相关，滤波器设计，参数模型等，几乎都只用一条语句即可调用。数字信号处理常用的函数有波形的产生、滤波器的分析和设计、傅里叶变换等，如：波形产生：sawtooth（锯齿波或三角波）Diric（Dirichlet或周期sinc函数）rand（白噪声信号波形）square（方波）sinc（sinc或 函数）chirp（chirp信号波形）滤波器的分析：abs（求幅值）angle（求相角）conv（求卷积）freqz（数字滤波器频率响应）impz（数字滤波器的冲击响应）zplane（数字系统零极点图）IIR滤波器设计：butter（巴特沃思数字滤波器）cheby1（切比雪夫I型）cheby2（切比雪夫II型）maxflat（最平滤波器）ellip（椭圆滤波器）yulewalk（递归数字滤波器）bilinear（双线性变换）impinvar（冲激响应不变法）FIR滤波器设计：triang（三角窗）blackman（布莱克曼窗）boxcar（矩形窗）hamming（海明窗）hanning（汉宁窗）kaiser（凯塞窗）fir1（基于窗函数法）fir2（基于频率抽样法）firrcos（上升余弦FIR滤波器设计法）intfilt（内插FIR滤波器设计法）kaiserord（用Kaiser窗设计FIR滤波器的参数估计3.滤波器的设计3.1巴特沃斯IIR滤波器的设计在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。butter函数的用法为：b,a=butter(n,Wn,/ftype/)其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：n,Wn= buttord(Wp,Ws,Rp,Rs)其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则：a高通滤波器：Wp和Ws为一元矢量且WpWs；b低通滤波器：Wp和Ws为一元矢量且WpWs；c带通滤波器：Wp和Ws为二元矢量且WpWs，如Wp=0.1,0.8,Ws=0.2,0.7。例如，采样频率 10Hz，通带截止频率 fp=3Hz，阻带截止频率 fs=4Hz，通带衰减小于 1dB，阻带衰减大于 20dB，使用双线性变换法由模拟滤波器原型设计数字滤波器。MATLAB程序见附录，如图1所示。图1 ButterWorth 低通滤波器3.2 契比雪夫I型IIR滤波器的设计在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。cheby1函数的用法为：b,a=cheby1(n,Rp,Wn,/ftype/)在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。cheblord函数的用法为：n,Wn=cheblord(Wp,Ws,Rp,Rs)其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。譬如，采样频率为10Hz，设计一个数字低通滤波器，要求其通带临界频率 fp= 3Hz ，通带p= 1dB )，阻带临界频率fs= 4Hz ，阻带内衰减大于 15dB（s= 15）内衰减小于 1dB(s= 1dB)。其程序见附录，设计出的滤波器如图2图2 Chebyshev 低通（I 型）滤波器3.3数字滤波器的设计数字滤波器(digital filter)是由数字乘法器、加法器和延时单元组成的一种装置。其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。FIR数字滤波器的单位脉冲响应是有限长序列。它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题，设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。在对滤波器实际设计时，整个过程的运算量是很大的。例如利用窗函数法设计M阶FIR低通滤波器时，首先要根据（1）式计算出理想低通滤波器的单位冲激响应序列，然后根据（2）式计算出M个滤波器系数。当滤波器阶数比较高时，计算量比较大，设计过程中改变参数或滤波器类型时都要重新计算。（1）（2）设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频相频响应特性，运算量也是很大的。我们平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很多时候都是要根据设计要求和滤波效果不断的调整，以达到设计的最优化。在这种情况下，滤波器的设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内完成设计。利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效的设计数字滤波器，大大的简化了计算量，直观简便。FIR数字滤波器的窗函数设计法设理想带阻滤波器频率响应为 例如，利用Kaiser窗函数，设计长度为55的阻带滤波器，使阻带衰减为60dB。参数可由式 =0.1102（-8.7）其中=60，确定MATLAB程序见附录，如图3所示。与其他高级语言的程序设计相比,MA TLAB环境下可以更方便、快捷地设计出具有严格线性相位的FIR 滤波器,节省大量的编程时间,提高编程效率,且参数的修改也十分方便. 还可以进一步进行优化设计,相信随着版本的不断更新,MA TLAB在数字滤波器技术中必将发挥更大的作用.IIR数字滤波器的设计方法IIR数字滤波器最大的优点是给定一组指标时，他的阶数要比相同组的FIR滤波器的低得多。基于Matlab的IIR滤波器的设计方法主要有完全设计法、模拟原型设计法和直接设计法等。a. 模拟原型法采用经典低通滤波器作为连续域上的设计模型，通过频域变换得到IIR数字滤波器，最后还要进行离散化处理。Matlab提供的低通模拟滤波器原型函数包括：besselap，buttap，cheb1lp，cheb2ap，ellipap；频域变换函数包括：lp2bp，lp2bs，lp2hp，lp2lp；离散化处理函数有bilinear和impinvar。图3带阻滤波器幅频特性b. 完全设计法Matlab信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，直接调用就可以设计滤波器，这为设计通用滤波器提供了方便。设计Butterworth滤波器用函数butter()，可以设计低通、高通、带通和带阻的数字和模拟滤波器，其特性是通带内的幅度响应最大限度的平滑，但损失了截止频率处的下降斜度。设计Chebyshev I型滤波器用函数chebyl()。可以设计低通、高通、带通和带阻的数字和模拟ChebyshevI型滤披器，其通带内为等波纹，阻带内为单调。Chebyshev I型滤波器的下降斜度比II型大，但其代价是通带内波纹较大。设计Chebyshev II型滤波器用函数cheby2()。可以设计低通、高通、带通和带阻的数字和模拟Chebyshev II型滤波器，其通带内为单调，阻带内等波纹。Chebyshev II型滤波器的下降斜度比I型小，但其阻带内波纹较大。设计椭圆滤波器用函数ellip()，与cheby1，cheby2类似，可以设计低通、高通、带通和带阻的数字和模拟滤波器。与Butterworth和chebyshev滤波器相比，ellip函数可以得到下降斜度更大的滤波器，得到通带和阻带均为等波纹。一般情况下，椭圆滤波器能以最低的阶实现指定的性能指标。c. 直接设计法直接设计方法的思想是基于给定的滤波器参数直接在离散域上寻找合适的数字滤波器，他不限于常规的滤波器类型，如低通、高通、带通和带阻等。这种方法甚至可以设计多带的频率响应，Matlab提供yulewalk函数用于辅助设计。设计一个数字低通滤波器满足以下参数，采样频率为 1Hz，通带临界频率 fp =0.2Hz，通带内衰减小于 1dB（p=1）；阻带临界频率 fs=0.3Hz，阻带内衰减大于 25dB（s=25）。Matlab 使用归一化的频率参数（临界频率除以采样频率的 1/2）。这样临界频率参数的取值范围在 0 和 1 之间，1 代表 Fs/2（用角频率表示的时候对应）。MATLAB源程序见附录，如图4所示。图4 直接设计的数字高通滤波器d. 通用Butterworth设计方法使用这种方法设计的Butterworth数字滤波器可以有不同数目的零点和极点，Matlab提供的maxflat函数实现了这一功能。这个函数与butter函数很相似，但他可以指定两个阶参数，其中归一化和非归一化各一个。如果这两个参数的值相同，那么他与butter函数的结果就是相同的。例如butter()函数设计的滤波器，其MATLAB程序见附录，增益图如图5所示。e. 参数建模法寻找接近于所需要设计的滤波器的通用模型，时域上的建模函数为lpc，prony，Stmcb；频域上的建模函数有invfreqs和invfreqz。典型IIR数字滤波器的设双线性变换法为了克服冲激响应不变法的频率混叠现象，需要使s平面与z平面建立一一对应的单值映射关系，即求出s=f（z），然后将它带入H（s），就可以求得H（z），即 H（z）=H（s）|s=f（z） 为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，其次再通过上面讨论过的标准变换关系z=es1T将此横带变换到整个z平面上去，这样就使s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。例如，双线性变换法设计一个Chebyshev高通滤波器，使其幅频特性逼近一个具有以下技术指标的模拟Chebyshev高通滤波器：Ws=2*pi*1kHz，Wp=2*pi*1.4kHz，在Ws处的最小衰减为15dB，在Wp处的最大衰减不超过0.3dB，抽样频率为20kHz。其MATLAB程序见附录，频率响应如图6图5 Butterworth带阻滤波器增益图图6 Chebyshev高通滤波器的频率响应可以看到经过离散采样、数字滤波后分离出了83.3hz的频率分量。之所以选取上面的叠加信号作为原始信号，是由于在实际工作中是要对已经经过差分滤波的信号进一步做带通滤波，信号的各分量基本同一致，可以反映实际的情况。本例设计的滤波器已在实际工作中应用，取得了不错的效果。结论采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是Matlab工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。MATLAB 信号处理工具箱为滤波器设计及分析提供了非常优秀的辅助设计工具, 在设计数字滤波器时, 善于应用MATLAB进行辅助设计, 能够大大提高设计效率。参 考 文 献1 陈德树. 计算机继电保护原理与技术.北京：水利电力出版社，1992.2 蒋志凯. 数字滤波与卡尔曼滤波.北京：中国科学技术出版社，19933 楼顺天、李博菡. 基于MATLAB的系统分析与设计信号处理.西安：西安电子4 董长虹等. MATLAB信号处理与应用.北京：国防工业出版社，2005 5 美 M.H.海因斯 著，张建华等 译.数字信号处理.北京：科学出版社，20026 张葛祥，李 娜. MATLAB仿真技术与应用.北京：清华大学出版社，2002附录1、Chebyshev高通滤波器的频率响应ws=2*pi*1000;ws1=ws*2*pi;wp=2*pi*1400;wp1=wp*2*pi;rp=0.3;rs=15;fs=20000;N,Wn=cheb2ord(wp1,ws1,rp,rs,s);z,p,k=cheb2ap(N,rs);A,B,C,D=zp2ss(z,p,k);At,Bt,Ct,Dt=lp2hp(A,B,C,D,Wn);At1,Bt1,Ct1,Dt1=bilinear(At,Bt,Ct,Dt,fs);num,den=ss2tf(At1,Bt1,Ct1,Dt1);freqz(num,den);H,W=freqz(num,den);plot(W*fs/(2*pi),abs(H);grid;xlabel(频率/Hz);ylabel(幅值);2、带阻滤波器幅频特性实现n=55-1;w=0.4 0.6;beta=0.1102*(60-8.7);kaiw=Kaiser(n+1,beta);b=fir1(n,w,stop,kaiw);h,w=freqz(b,1,512,2);plot(w,20*log10(abs(h);grid;xlabel(频率（归一化）);ylabel(幅度dB);3、Chebyshev 低通（I 型）滤波器T=0.1; FS=1/T; fp=3;fs=4;Rp = 1; As = 15;wp=fp/FS*2*pi; ws=fs/FS*2*pi; OmegaP = (2/T)*tan(wp/2);OmegaS = (2/T)*tan(ws/2);ep = sqrt(10(Rp/10)-1);A = 10(As/20);OmegaC = OmegaP;OmegaR = OmegaS/OmegaP;g = sqrt(A*A-1)/ep;N = ceil(log10(g+sqrt(g*g-1)/log10(OmegaR+sqrt(OmegaR*OmegaR-1);z,p,k = cheb1ap(N,Rp); a = real(poly(p);aNn = a(N+1);p = p*OmegaC;a = real(poly(p);aNu = a(N+1); k = k*aNu/aNn;b0 = k;B = real(poly(z);b = k*B;bz,az = bilinear(b,a,FS);%freqz(bz,az,200,FS,whole);H=freqz(bz,az,200,whole);plot(abs(H);xlabel(频率/Hz);ylabel(幅值);4、直接设计数字滤波器FS=1;n,Wn=buttord(0.3/(FS/2),0.2/( FS /2),1,25);b,a=butter(n,Wn,high);freqz(b,a,512, FS);5、ButterWorth 低通滤波器T=0.1; FS=1/T; fp=3;fs=4;wp=fp/FS*2*pi; ws=fs/FS*2*pi; Rp = 1; As = 20; OmegaP = (2/T)*tan(wp/2);OmegaS = (2/T)*tan(ws/2); N = ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(OmegaP/OmegaS);OmegaC = OmegaP/(10(Rp/10)-1)(1/(2*N);z,p,k = buttap(N); p = p*OmegaC;k = k*OmegaCN; B = real(poly(z);b0 = k;cs = k*B;ds = real(poly(p);b,a = bilinear