基于MATLAB的IIR和FIR滤波器的设计与实现

基于MATLAB的IIR和FIR滤波器的设计与实现陈XX（XXX学院电信XX班）摘要：数字滤波是数字信号处理的重要内容，是由乘法器、加法器和单位延时器组成的一种运算过程，其功能是对输人离散信号进行运算处理，以达到改变信号频谱的目的。数字滤波器根据频域特性可分为低通、高通、带通和带阻四个基本类型。本文用脉冲响应不变法设计的一个满足指标的巴特沃斯IIR滤波器，利用了一种基于Matlab软件的数字滤波器设计方法，完成了低通，高通，带通，帯阻IIR滤波器的设计,文中深入分析了该滤波器系统设计的功能特点、实现原理以及技术关键，阐述了使用MATLAB进行带通滤波器设计及仿真的具体方法。最后把整个设计方案用GUIDE界面制作并演示出来。文章根据IIR滤波器的设计原理，重点介绍巴特沃斯数字滤波器的设计方法和操作步骤，并以实例形式列出设计程序。关键词：信号巴特沃斯MatlabIIR滤波器脉冲响应不变法一、引言在信号处理过程中，所处理的信号往往混有噪音，从接收到的信号中消除或减弱噪音是信号传输和处理中十分重要的问题。根据有用信号和噪音的不同特性，提取有用信号的过程称为滤波，实现滤波功能的系统称为滤波器。在近代电信设备和各类控制系统中，数字滤波器应用极为广泛。数字滤波器精确度高、使用灵活、可靠性高，具有模拟设备所没有的许多优点，已广泛地应用于各个科学技术领域,例如数字电视、语音、通信、雷达、声纳、遥感、图像、生物医学以及许多工程应用领域。随着信息时代数字时代的到来，数字滤波技术已经成为一门极其重要的学科和技术领域。以往的滤波器大多采用模拟电路技术，但是，模拟电路技术存在很多难以解决的问题，例如，模拟电路元件对温度的敏感性，等等。而采用数字技术则避免很多类似的难题，当然数字滤波器在其他方面也有很多突出的优点，在前面部分已经提到，这些都是模拟技术所不能及的，所以采用数字滤波器对信号进行处理是目前的发展方向。二、IIR数字滤波器的设计2.1IIR滤波器的基本结构一个数字滤波器可以用系统函数表示为：01()()()1MkkkNkkkbzYzHzXzaz(2-1)由这样的系统函数可以得到表示系统输入与输出关系的常系数线形差分程为：00()()()NMkkkkynaynkbxnk(2-2)可见数字滤波器的功能就是把输入序列x(n)通过一定的运算变换成输出序列y(n)。不同的运算处理方法决定了滤波器实现结构的不同。无限冲激响应滤波器的单位抽样响应h(n)是无限长的，其差分方程如(2-2)式所示，是递归式的，即结构上存在着输出信号到输入信号的反馈，其系统函数具有(2-1)式的形式，因此在z平面的有限区间(0︱z︱∞)有极点存在。前面已经说明，对于一个给定的线形时不变系统的系统函数，有着各种不同的等效差分方程或网络结构。由于乘法是一种耗时运算，而每个延迟单元都要有一个存储寄存器，因此采用最少常熟乘法器和最少延迟支路的网络结构是通常的选择，以便提高运算速度和减少存储器。然而，当需要考虑有限寄存器长度的影响时，往往也采用并非最少乘法器和延迟单元的结构。IIR滤波器实现的基本结构有：(1)IIR滤波器的直接型结构；优点：延迟线减少一半，变为N个，可节省寄存器或存储单元；缺点：其它缺点同直接I型。通常在实际中很少采用上述两种结构实现高阶系统，而是把高阶变成一系列不同组合的低阶系统(一、二阶)来实现。(2)IIR滤波器的级联型结构；特点：系统实现简单，只需一个二阶节系统通过改变输入系数即可完成；极点位置可单独调整；运算速度快(可并行进行)；各二阶网络的误差互不影响，总的误差小，对字长要求低。缺点：不能直接调整零点，因多个二阶节的零点并不是整个系统函数的零点，当需要准确的传输零点时，级联型最合适。(3)IIR滤波器的并联型结构。优点：简化实现，用一个二阶节，通过变换系数就可实现整个系统；极、零点可单独控制、调整，调整α1i、α2i只单独调整了第i对零点，调整β1i、β2i则单独调整了第i对极点；各二阶节零、极点的搭配可互换位置，优化组合以减小运算误差；可流水线操作。缺点：二阶阶电平难控制，电平大易导致溢出，电平小则使信噪比减小。2.2用脉冲响应不变法设计IIR滤波器2.2.1IIR低通数字滤波器实例数字低通的技术指标为：Wp=0.2πradAp=1dBT=2sWs=0.3πradAs=15dB程序为：T=2;%采样周期fs=1/T;%采样频率为采样周期倒数Wp=0.2.*pi;Ws=0.3.*pi;%设计归一化通带阻带截止频率Ap=1;As=15;%设置通带最大最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'s');%调用butter函数设计巴特沃斯滤波器W=linspace(0,pi,400*pi);%指定一段频率值[D,C]=impinvar(B,A,fs);%调用脉冲不变法¨Hz=freqz(D,C,W);%·返回频率值¦plot(W/pi,abs(Hz)/abs(Hz(1)));%绘出巴特沃斯数字低通滤波器的幅频特性曲线gridon;%显示栅格title('巴特沃斯低通数字滤波器');%标题xlabel('Frequency/Hz');%x轴标签ylabel('Magnitude');%y轴标签gtext(‘100230陈外流’)得出幅频特性如图1所示：图12.2.2IIR高通数字滤波器数字高通的技术指标为：Wp=0.4πradAp=2dBT=2sWs=0.2πradAs=15dB程序为：T=2;%采样周期fs=1/T;%采样频率Wp=0.4.*pi;Ws=0.2.*pi;%设置归一化通带和阻带截止平率Ap=2;As=15;%设置通带最大最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'high','s');%调用butter函数设计巴特沃斯W=linspace(0,pi,400*pi);%指定一段频率值[D,C]=impinvar(B,A,fs);%调用脉冲响应不变法·¨Hz=freqz(D,C,W);%返回频率响应¦plot(W/pi,abs(Hz)/abs(Hz(1)));%绘制巴特沃斯数字高通滤波器的幅频响应曲线gridon;%显示栅格title('巴特沃斯高通数字滤波器');%标题xlabel('Frequency/Hz');%x轴标签ylabel('Magnitude');%y轴标签gtext(‘100230陈外流’)得出幅频特性如图2所示：图22.2.3IIR带通数字滤波器带通滤波器技术指标为：Wp=[0.25π0.35π]radAp=1dBT=2sWs=[0.15π0.40π]radAs=10dB程序为;T=2;%采样周期fs=1/T;%采样频率Wp=[0.25.*pi0.35.*pi];Ws=[0.15.*pi0.4.*pi];%设置归一化通带和阻带截止平率Ap=1;As=10;%设置通带最大最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'bandpass','s');%调用butter函数设计巴特沃斯滤波器W=linspace(0,pi,400*pi);%指定一段频率值[D,C]=impinvar(B,A,fs);%调用脉冲响应不变法¨Hz=freqz(D,C,W);%·返回频率值¦plot(W/pi,abs(Hz));%绘出巴特沃斯数字滤波器的幅频特性曲线gridon;%显示栅格title('巴特沃斯带通滤波器');%标题xlabel('Frequency/Hz');%x轴标签ylabel('Magnitude');%y轴标签gtext(‘100230陈外流’)得出幅频特性如图3所示：图32.2.4用脉冲响应不变法设计IIR带阻数字滤波器实例。带阻滤波器技术指标为：Wp=[0.15π0.40π]radAp=1dBT=2sWs=[0.25π0.35π]radAs=11dB程序为：T=2;%采用周期为2fs=1/T;%采样频率Wp=[0.15.*pi0.40.*pi];Ws=[0.25.*pi0.35.*pi];%设置归一化通带和阻带截止平率Ap=1;As=11;%设置通带最大最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'stop','s');%调用butter函数设计巴特沃斯滤波器÷W=linspace(0,pi,400*pi);%指定一段频率值[D,C]=impinvar(B,A,fs);%调用脉冲响应不变法·¨Hz=freqz(D,C,W);%·返回频率响应plot(W/pi,abs(Hz));%绘出巴特沃斯数字滤波器的幅频特性曲线gridon;%显示栅格title('巴特沃斯带阻滤波器');%显示标题‘巴特沃斯带阻滤波器xlabel('Frequency/Hz');%x轴标签gtext(‘100230陈外流’)得出幅频特性如图4所示：图4结论：经观察滤波器幅频特性图得，巴特沃斯滤波器的特点是通带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。低通滤波器容许低频信号通过,但减弱(或减少)频率高於截止频率的信号的通过。高通滤波器容许高频信号通过,但减弱(或减少)频率低于於截止频率的信号的通过。带通滤波器容许一定频率范围信号通过,但减弱(或减少)频率低于於下限截止频率和高于上限截止频率的信号的通过。带阻滤波器减弱(或减少)一定频率范围信号,但容许频率低于於下限截止频率和高于上限截止频率的信号的通过。对于巴特沃斯滤波器低通和带通效果较好，而高通和帯阻效果较差。2.3双线性变换法设计IIR滤波器2.3.1IIR数字高通滤波器fs=150;fp=250;%模拟技术指标Fs=1000;T=1/Fs;%fs为采样频率，fp为中心频率wp=fp/Fs*2*pi;%滤波器的通带截止频率ws=fs/Fs*2*pi;%滤波器的阻带截止频率Rp=1;As=20;%滤波器的通阻带衰减指标ripple=10^(-Rp/20);%滤波器的通带衰减对应的幅度值Attn=10^(-As/20);%滤波器的阻带衰减对应的幅度值Omgp=(2/T)*tan(wp/2);%原型通带频率的预修正Omgs=(2/T)*tan(ws/2);%原型阻带频率的预修正[n,Omgc]=ellipord(Omgp,Omgs,Rp,As,'s')%计算阶数n和截止频率[z0,p0,k0]=ellipap(n,Rp,As);%设计归一化的椭圆模拟滤波器原型ba=k0*real(poly(z0));%求原型滤波器的系数baa=real(poly(p0));%求原型滤波器的系数a[ba1,aa1]=lp2hp(ba,aa,Omgc);%变换为模拟高通滤波器[bd,ad]=bilinear(ba1,aa1,Fs)%求数字系统的频率特性[H,w]=freqz(bd,ad);%绘制模拟滤波器频响特性dbH=20*log10((abs(H)+eps)/max(abs(H)));subplot(2,1,1),plot(w/2/pi*Fs,abs(H),'k');ylabel('|H|');title('幅度响应');axis([0,Fs/2,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);gridsubplot(2,1,2),pl