《基于matlab的iir数字滤波器设计.doc》由会员分享，可在线阅读，更多相关《基于matlab的iir数字滤波器设计.doc（19页珍藏版）》请在金锄头文库上搜索。
1、 基于matlab的IIR数字滤波器设计一 IIR数字滤波器介绍1.IIR数字滤波器的基本原理所谓数字滤波器，是指输入，输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的硬件。实质上就是一个由有限精度算法实现的线性时不变离散系统。它的基本工作原理是利用离散系统的特性对系统输入信号进行加工和变换，改变输入序列的频谱或信号波形，让有用的频率分量通过，抑制无用的信号分量输出，因此数字滤波与模拟滤波的概念相同，根据其频率特性同样可以分为低通，高通，带通，带阻，只是信号的形式和实现滤波方式有所不同。如果要处理的信号是模拟信号，就可以通过A/D或者D/A转换，在信号
2、形式上进行匹配转换，同样可以使用数字滤波器对模拟信号进行滤波。数字滤波器滤波的数学表达式：y（n）=x(n)*h(n); 如果滤波器的输入输出信号都是离散信号，那么该滤波器的脉冲响应也一定是离散信号，这样的滤波器就成为了数字滤波器。上面的系统为时域离散系统时，其频域特性为：其中分别是数字滤波器的输出序列和输入序列的频域响应，是数字滤波器的频域响应。可以看见按照输入信号的频谱特点和处理信号的目的适当选择滤波器的频域响应，使得滤波后的输出信号满足设计性能要求，就是滤波器的滤波原理。2IIR数字滤波器传输特性IIR数字滤波器的系统函数可以表示为：H(Z)=，式中H(Z)称为N阶IIR滤波器函数。3.
3、数字滤波器的技术要求.我们通常设计的数字滤波器一般属于选频滤波器，。我们的目的是要设计一个因果可实现的滤波器，另外买也要考虑到成本和复杂性问题，因此实用中通带和阻带都允许一定的误差容限，即通带不一定是完全水平的，阻带也不可能完全衰减到零。而且，通带和阻带之间还要设置一定带宽的过渡带。如下图表示低通滤波器的技术要求： 图中， 分别表示通带截止频率和阻带截止频率，通带频率范围为0w，通带中要求（1-1）|H1，阻带截止频率范围w，再阻带中要求2，从到称为过渡带，在这个频带内，幅度响应从通带平滑的下落到阻带。二 IIR数字滤波器的设计方法1常用的IIR滤波器设计方法 （1）. 以模拟滤波器函数为基础
4、的变换法；即先设计一满足指定条件的模拟滤波器H(s)，再将该模拟滤波器转化为数字滤波器H(z)。 （2）. 直接设计法：在z平面内，根据零、极点对系统特性的影响，调整零极点位置得H(z)。 （3）. 最优化设计法(计算机辅助设计)，在某种最小化误差准则下，建立差分方程系数 a k、b i 对理想特性的逼近方程，使用迭代方法解方程组得到最佳逼近系统。由于此方法计算量大，需要借助于计算机进行设计。下面分别介绍函数设计法和信号处理图形界面来设计IIR数字滤波器，得出最优化的设计方法。 2.基于matlab的函数设计IIR数字滤波器IIR数字滤波器设计的一般方法是先设计低通模拟滤波器，进行频率变换，将
5、其转换为相应的（高通，带通等）模拟滤波器，在转换为高通，带通或带阻数字滤波器，由模拟滤波器设计数字滤波器的方法。这是因为模拟滤波器设计方法已经很成熟，它不仅有完整的设计公式，还有完善的图表供查阅，另外，还有一些典型的滤波器类型可供我们使用。对设计的全过程的各个步骤，MATLAB都提供了了相应的工具箱函数，使IIR数字滤波器设计变得非常简单。2.1. 基于matlab函数的IIR滤波器设计2.1.1.设计方法选择： 程序设计法是基于MATLAB中相应的工具箱函数来实现的，IIR数字滤波器的设计步骤可由图1.2所示的流程图来表示。这个图也清晰的表示了5类20个信号处理工具箱函数的作用，在MATLA
6、B中，模拟滤波器的系统函数H(S)=数字滤波器的系统函数H(Z)= 在实际工程中，需要的设计结果是系数向量B和A，用B和A来综合滤波器的硬件实现结构或软件运算结构，为了直观的看出设计结果，本文的实例均以滤波器幅频响应曲线作为设计结果输出。如果需要滤波器系数，在运行程序后，只要在MATLAB命令窗口键入系数向量名，则相应的系数就显示出来了。2.1.2.程序设计实例分析（a）设计高通和带通Butterworth数字滤波器 我们给出四阶归一化 Butterworth模拟滤波器的系统函数 用双线性变换法从Ha（s）设计四阶带通butterworth数字滤波器,并图示|,设计采样周期T=1s，指标如下
7、现在我们分步进行：建模 由于本例主要涉及三个问题：（1） 由数字滤波器指标求相应的模拟滤波器指标；（2） 模拟滤波器频率变换（因为已给定阶数和模拟滤波器的归一化低通原型）；（3） 由相应的模拟滤波器到数字滤波器（双线性变换法）。由于调用bilinear函数将模拟滤波器转换成数字滤波器非常容易，并且有效抑制频率失真的问题，本例给定了数字滤波器指标，所以首先要设计处与该指标相应的四阶Butterworth模拟滤波器，然后调用bilinear函数将其转换为数字滤波器即可，应当特别注意的是，对于双线性变换法，由数字边界频率求相应的模拟边界频率时，一定要考虑预畸变矫正。只有这样，最终设计结果才能满足所给
8、指标，（）首先按照步骤一的要求 设计高通数字滤波器时，相应的模拟高通滤波器3dB截止频率为 设计带通数字滤波器时，相应的模拟滤波器的3dB截止频率为 （）步骤二的的原理 可调用MATLAB频率变换函数lp2lp，lp2hp，lp2bp，分别实现从模拟低通到模拟低通，高通，带通，带阻的频率变换。【Bt，At】=lp2hp（B，A,wc），将系数向量为B和A的模拟滤波器归一化低通原型（3Db）截止频频为1rad/s），变换成3dB截止频率为wc的高通模拟滤波器，返回高通模拟滤波器系数向量Bt和At。【Bt，At】=lp2bp（B,A，wo，Bw）将系数向量为B和A的模拟滤波器归一化低通原型变换成中
9、心频率为wo，带宽为Bw的带通模拟滤波器，返回带通模拟滤波器的系数向量Bt和At。其中，wo= ，由以上原理我们来编写如下程序：%用双线性变换法设计数字高通和带通滤波器clear;close allT=1;wch=pi/2; %T为采样间隔，wch位数字高通3dB截止频率wlc=0.35*pi;wuc=0.65*pi; %wlc，wuc；数字高通3dB截止频率B=1;A=1,2.6131,3.4142,2.6131,1;h,w=freqs(B,A,512); %求原归一化模拟滤波器的频率响应subplot(3,2,1);plot(w,20*log10(abs(h); %画模拟滤波器幅频特性 g
10、rid;axis(0,10,-90,0)xlabel(w/ );ylabel(模拟低通幅度（dB）)%(1)设计高通omegach=2*tan(wch/2)/T; %预畸变求模拟高通3dB截止频率Bhs,Ahs=lp2hp(B,A,omegach); %模拟域低通转换为高通系数 Bhz,Ahz=bilinear(Bhs,Ahs,1/T); %模拟转换位数字高通系数变量h,w=freqz(Bhz,Ahs,512); %求画出数字滤波器幅频特性Subplot(3,2,3);plot(w/pi,20*log10(abs(h);grid;axis(0,1,-150,0)xlabel(w/ );ylab
11、el(数字滤波器幅度（dB）)%(2)设计带通omegalc=2*tan(wlc/2)/T; %预畸变求滤波器通带低端截止频率omegauc=2*tan(wuc/2)/T; %预畸变求滤波器通带高端截止频率wo=sqrt(omegalc*omegauc);Bw=omegauc-omegalc;Bbs,Abs=lp2bp(B,A,wo,Bw); %模拟域低通转换为带通系数 Bbz,Abz=bilinear(Bbs,Abs,1/T); %模拟转换为数字带通系数变量 h,w=freqz(Bbz,Abz,512); %求并画出数字滤波器幅频特性subplot(3,2,4);plot(w/pi,20*l
12、og10(abs(h);grid;axis(0,1,-150,0)xlabel(w/pi );ylabel(数字滤波器幅度（dB）)程序运行结果如图1.3，1.4，1.5所示分别表示模拟低通幅度，数字高通幅度和数字带通幅度特性。模拟低通幅度 数字高通幅度 图1.3 图1.4 数字带通幅度 图1.5（b）设计滤波器进行图像去噪处理以上是我们运用MATLAB工具箱函数来编程设计IIR数字滤波器并对其频响特性进行分析的实例，下面我们就来引入数字滤波器处理图像的实例具体分析。我们用卷积定理来说明数字滤波器进行图像去噪处理的原理，设图像信号e(x, y)通过线性不变系统h(x, y)的结果是g(x, y
13、)，即r（x,y）=e(x,y)*h(x,y); 由时域卷积，频域乘积定理得R（U，V）=E（U，V）H（U，V）。其中，R(u, v)，E(u, v)，H(u, v)分别是r(x, y)，e（x, y)和h(x, y)的傅里叶变换。实际上，图像的能量大部分集中在幅度频的低频和中频段，而图像的边缘和噪声对应于高频部分。因此，能降低高频成分幅度的滤波器则能过滤噪声，减弱噪声的影响。而Butterworth低通滤波器在物理上是可以实现的滤波器，它的转移函数模的平方为,为N阶，截止频率为。下面我们就来设计Butterwirth滤波器来对加高斯白噪声的图像进行去噪处理I = imread(K:祖国好.jpg);figure, imshow(I),D = imnoise(I,gaussian);figure, imshow(D)D=double(D);F=fft2(D); %傅里叶变换F=fftshift(F); %转换数据矩阵N1,N2=size(F);n=2;d0=40;n1=fix(N1/2);n2=fix(N2/2);for i=1:N1for j=1:N2 d =sqrt(i-n1)2+(j-n2)2); h=1/(1+(d/d0)(2*n); %计算低通转换函数 FD(i,j)=h*F(i,j); %低通滤波e
