基于MATLAB的IIR数字滤波器设计

基于MATLAB的IIR数字滤波器设计肖逸20122301机自05班2016-3-24机自05班20122301肖逸1/13基于MATLAB的IIR数字滤波器设计一概述数字滤波器由数字乘法器、加法器和延时单元组成的一种算法或装置。数字滤波器的功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。数字滤波器对信号滤波的方法是：用数字计算机对数字信号进行处理，处理就是按照预先编制的程序进行计算。数字滤波器的原理如图1所示，它的核心是数字信号处理器。图1数字滤波器的方框图如果采用通用的计算机，随时编写程序就能进行信号处理的工作，但处理的速度较慢。如果采用专用的计算机芯片，它是按运算方法制成的集成电路，连接信号就能进行处理工作，处理的速度飞快，但功能不易更改。如果采用可编程的计算机芯片，那么，装入什么程序机器就能具有什么功能。这种可编程芯片的优点很多，是现代电子产品的首选。如果是对模拟信号进行处理，则需要添加模数转换器和数模转换器。二数字滤波器的特性2.1滤波器基础知识数字滤波器按频率特性也有低通、高通、带通、带阻全通等类型。由于频率响应的周期性，频率变量以数字频率ω来表示（ω=ΩΤ=Ω/fs，Ω位模拟角频率，Τ位抽样时间间隔，fs为抽样频率），所以数字滤波器设计中必须给出抽样频率。图2为各种数字滤波器的理想频率响应（只表示了正频率部分）。在图2中，2π等于数字域抽样频率，即ωs=ΩsT=2πfsT=2πfsfs=2π，T=1fsωs2⁄=π是折叠频率，按照奈奎是特采样定理，频率特性只能限于折叠频率以内，即限于|ω|ωs2⁄=π范围之内。一般情况下，数字滤波器是一个线性移时不变离散时间系统，利用有限精度算法来实现。机自05班20122301肖逸2/132.2数字滤波器的性能要求一般说来，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波器为例，如图3所示，频率响应有通带、过渡带及阻带三个范围（而不是理想的陡截止的通带、阻带两个范围）。在通带内，幅度响应应以误差α1逼近于1，即1−α≤H(ejω)≤1，|ω|≤ωc在阻带中，幅度响应应以误差小于α2而逼近于零，即|H(ejω)|≤α2，ωst≤|ω|≤π其中ωc、ωst分别为通带截止频率和阻带截止平率，它们都是数字域频率。为了逼近理想的低通滤波器特性，还必须有一个非零宽d度（ωst−ωc）的过渡带，在这个过渡带内的频率响应应平滑地从通带降到阻带。2.3表征数字滤波器频率响应的参量2.3.1幅度平方响应当只需要逼近幅度响应而不考虑相位时，例如标准的低通、高通、带通、带阻滤波器的逼近就是这样，这时根据幅度平方响应来进行设计时很方便的。幅度平方响应定义为|H(ejω)|2=H(ejω)H∙(ejω)=H(ejω)H(e−jω)=|H(z)H(z−1)|z=ejω这里由于冲击响应为实函数，故满足H∙(ejω)=H(e−jω)，也就是满足共轭对称条件。2.3.2相位响应由于H(ejω)是复数，可表示成H(ejω)=|H(ejω)|ejβ（ejω)=Re[H(ejω)]+jIm[H(ejω)]所以β(ejω)=arctan⁡{Im[H(ejω)]Re[H(ejω)]}由于H∙(ejω)=H(e−jω)所以又有β(ejω)=12jln[H(ejω)H∙(ejω)]=12jln[H(ejω)H(e−jω)]=12jln[H(z)H(z−1)]z=ejω2.3.3群延迟响应它是滤波器平均延迟的一个度量，定义为相位对角频率的倒数的负值，即τ(ejω)=−dβ（ejω)dω可以化为τ(ejω)=−dβ(z)dzdzdω|z=ejω=−jzdβ(z)dz|z=ejω由于机自05班20122301肖逸3/13ln[H(ejω)]=ln|H(ejω)|+jβ(ejω)所以β(ejω)=Im{ln[H(ejω)]}因而又有τ(ejω)=−Im[ddωlnH(ejω)]同样可以化为τ(ejω)=−Im[d{ln[H(z)]]dzdzdω]z=ejω=−Im[jzd{ln[H(z)]}dz]z=ejω=−Re[zddz{ln[H(z)}]z=ejω=−Re[zdH(z)dz1H(z)]z=ejω当要求滤波器为线性相位响应特性时，则通带内群延迟特性就是常数。三窗函数窗函数是一种以截断信号进行操作的方式，在功率谱估计和滤波器设计中有着重要的应用。加窗可以减少由于无限序列截断带来的吉布斯效应。在选用的时候，需要根据信号类型和处理目的进行选择。常见的窗函数如表1所示。窗函数MATLAB实现函数矩形窗rectwin巴特窗/三角窗Bartlett/triang海明窗hamming汉宁窗hann凯瑟窗kaiser表1常见的窗函数下面结合实例对窗函数做简单的说明。3.1矩形窗建立一个20点长度的矩形窗，得到频率响应曲线如图4，左边为时域窗函数，右边为频域窗函数。机自05班20122301肖逸4/13图420点长度的矩形窗3.2巴特窗/三角窗巴特窗/三角窗是两个矩形窗的卷积。建立不同长度的三角，比较两者的区别。图5、6分别为长度为9、10的三角窗，同样，左边为时域窗函数，右边为频域窗函数。图59点长度的三角窗机自05班20122301肖逸5/13图610点长度的三角窗3.3海明窗创建长度为64点的海明窗，绘制出频率响应曲线，如图7所示。图764点长度的海明窗3.4汉宁窗创建64点长度的汉宁窗，并绘制频率响应曲线，如图8所示。机自05班20122301肖逸6/13图864点长度的汉宁窗3.4凯瑟窗创建64点长度的凯瑟窗，并绘制频率响应曲线，如图9所示图964点长度的凯瑟窗4IIR滤波器在MATLAB中进行信号滤波处理，使用的滤波器主要包括两类：1）无限冲激响应滤波器（IIR）。2）有限冲击响应滤波器（FIR）。4.1IIR滤波器的优点IIR滤波器相较于FIR滤波器，需要使用的阶次更少，却包含了非线性相位。但是由于MATLAB处理信号一般为离散信号，所以在滤波处理中允许一些非因果效应的存在。机自05班20122301肖逸7/134.2经典滤波器的设计过程IIR滤波器的设计方法一般有如下三种：1）简单滤波器的零极点位置累试法2）利用模拟滤波器的设计理论设计IIR数字滤波器3）利用最优化技术来设计滤波器参数其中，利用模拟滤波器的设计理论设计IIR数字滤波器的方法是目前应用最广泛的设计方法，其设计过程为：先设计合适的模拟滤波器，再变换为满足预期目标的数字滤波器。4.3经典法IIR滤波器设计常见的经典法设计的IIR滤波器包括以下四种：1）巴特沃斯滤波器2）切比雪夫I型滤波器3）切比雪夫II型滤波器4）椭圆滤波器4.3.1巴特沃斯滤波器巴特沃斯滤波器可用于设计低通、高通、带通、带阻滤波器和模拟滤波器，其实现过程在MATLAB分为两步：首先根据滤波器要求用buttord函数求出最小滤波器阶数和截止频率，然后求出滤波器系数。其中有几个个较为关键的参量Wp和Ws分别是通带和阻带的拐角归一化频率（截止频率），其取值范围为0~1），当其值为1时，代表采样频率的一半。Rp和Rs分别是通带和阻带的波纹系数。不同类型的滤波器对应的Wp和Ws遵循以下原则：1）高通滤波器：Wp和Ws为一元矢量切WpWs。2）低通滤波器：Wp和Ws为一元矢量切WpWs。3）带通滤波器：Wp和Ws为二元矢量切WpWs，如Wp=[0.20.7]，Ws=[0.10.8]。4）带阻滤波器：Wp和Ws为二元矢量切WpWs，如Wp=[0.10.8]，Ws=[0.20.7]。例如：要求对频率为20HZ和200HZ的单频叠加谐波信号的低通滤波，要求使得输出仅含有20HZ。得到图形10、11、12分别代表原始信号，滤波器情况，和处理后的信号。机自05班20122301肖逸8/13图10原始信号时域和频域图图11滤波器频率响应情况图12输出信号时域和频域图4.3.2切比雪夫I型滤波器切比雪夫I型滤波器的MATLAB设计和实现过程与巴特沃斯滤波器的设计实现过程基本一致。假设要设计一个采样频率为1000Hz的信号，使其在通带内（0~40Hz）波动不超过3dB,在150Hz处至少衰减60dB。则可以得到如图13所示的低通滤波器的频域响应情况。机自05班20122301肖逸9/13图13切比雪夫I型低通滤波器4.3.3切比雪夫II型滤波器其设计过程与前两种类似。以同样的要求设计一个切比雪夫II型低通滤波器，得到如图14所示的低通滤波器的频域响应情况。图14切比雪夫II型低通滤波器4.3.4椭圆滤波器其设计过程与前两种类似。以同样的要求设计一个椭圆低通滤波器，得到如图15所示的低通滤波器的频域响应情况。机自05班20122301肖逸10/13图15椭圆滤波器4.4直接法IIR滤波器设计直接法与前几种设计方法不一样，它直接针对滤波器的频率响应进行设计。例如，设计一个8阶的低通滤波器，将其频率响应与设计目标相比较，如图16所示。图16直接法设计IIR滤波器5总结通过几天的学习，我学会很多关于信号处理的知识，了解了许多滤波的算法，学会了不少的MATLAB中对于滤波器设计的过程，看到了不同滤波方法的不同点和相同点。总而言之，通过这一次的论文写作，我受益良多，同时在课程的学习中也学到了很多，不仅仅只是知识，还有很多的人生经验。机自05班20122301肖逸11/136参考资料[1]韩晶、刘浩编著，《MATLABR2014a完全自学一本通》，电子工业出版社，2015.1[2]程佩青，《数字信号处理教程》，清华大学出版社，2007.2[3]百度文库，百度百科5附录（代码）矩形窗n=20;w=rectwin(n);wvtool(w)三角窗w1=bartlett(9);w2=bartlett(10);wvtool(w1)wvtool(w2)海明窗l=64;wvtool(hamming(l))汉宁窗l=64;wvtool(hann(l))凯瑟窗w=kaiser(200,2.5);wvtool(w)巴氏低通滤波器fs=1200;N=300;n=0:N-1;t=n/fs;fL=20;fH=200;s=cos(2*pi*fL*t)+cos(2*pi*fH*t);subplot(121);plot(t,s)title('输入信号');xlabel('t/s');ylabel('幅度');sfft=fft(s);subplot(122);plot((1:length(sfft)/2)*fs/length(sfft),2*abs(sfft(1:length(sfft)/2))/length(sfft));title('信号频谱');xlabel('频率/Hz');ylabel('幅度');Wp=50/fs;Ws=100/fs;[n,Wn]=buttord(Wp,Ws,1,50);[a,b]=butter(n,Wn);[h,f]=freqz(a,b,'whole',fs);f=(0:length(f)/2)'*fs/length(f);figure(2);plot(f(1:length(f)/2),abs(h(1:length(f)/2)));机自05班20122301肖逸12/13title('巴士低通滤波器');xlabel('频率/Hz');ylabel('幅度');grid;sF=filter(a,b,s);figure(3);subplot(121);plot(t,sF);title('输出信号');xlabel('t/s');ylabel('幅度');SF=fft(sF);subplot(122);plot((1:length(SF)/2)*fs/length(SF),2*abs(SF(1:length(SF)/2))/length(SF));title('低通滤波后的频谱图');xlabel('频率/Hz');ylabel('幅度');切比雪夫I型低通滤波器Wp=40/500;Ws=150/500;Rp=3;Rs=60;[nWp]=cheb1ord(Wp,Ws,Rp,Rs);[