IIR数字滤波器课程设计汇总

1、数字信号处理课程设计报告基于MATLAB勺IIR数字滤波器设计专业班级:电信工程1302班学 号:311308000626学生姓名:王海龙指导教师:王科平2016年7 月目录摘要、课程设计任务及要求 41. 本次设计的目的 42. 本次设计的要求 4、课程设计原理 41. 脉冲响应不变法原理 42. 双向性变换法原理 5三、IIR数字滤波器设计内容51. 总体方法分析 52. 脉冲相应不变法 63. 双线性变换法 7 9911.1119.1922四、IIR 数字滤波器设计过程1. 设计步骤 2. 程序流程框图 3. MATLAB1序4. 调试分析过程描述 5 结果分析 五、结论 六、参考文献

2、23在当今社会，数字信号处理技术飞速发展，它不但自成一门学科，更是以不 同的方式影响和渗透到其他学科的研究中，它变得与我们的生活联系越来越紧 密，不断改变着我们的生产生活方式，因此受到人们越来越多的关注。数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一 种运算过程，实现对信号的运算处理。输入数字信号(数字序列)通过特定的运 算转变为输出的数字序列。描述离散系统输出与输入关系的卷积和差分方程只是 给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时 域离散系统的频域特性：Y(ejw)=X(e jw)H(ejw),其中八、分别是数字滤波器的输出序列和输入序列的频

3、域特性(或称为频谱特性),H(ejw)是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。数字滤波器是具有一定传输选择特性的数字信号处理装置 ，其输入、输出均 为数字信号，实质上是一个由有限精度算法实现的线性时不变离散系统。IIR数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现 ，其差NN分方程为：y(n)= ai X(n-i)+ bi y(n-i)i =0i#MN系统函数为：H(z)=(、brZ)/( 1+akZk)r =0k -0设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数 H(z)，使其频 率响应H(z)满足所希望得到的频域指标。本次课程设计分别用

4、脉冲响应不变法、双向性变换法设计IIR低通、高通、带通、带阻滤波器滤波器。并在MATLA环境下实现了 IIR数字滤波器的设计和仿 真。其主要内容概括为：首先对滤波器的原理和设计进行了介绍； 接着描述了 IIR 数字滤波器的基本概念，其中包括系统的描述、系统的传递函数、系统的模型； 接着简单介绍MATLAB并对数字滤波器在MATLA环境下如何实现进行了介绍；重 点描述了 IIR数字滤波器的设计过程，最后对IIR滤波器进行仿真。关键词：数字滤波器频域特性脉冲响应双向性变换法 MATLAB、课程设计任务及要求1. 本次设计的目的1 )学会MATLA的使用，掌握MATLAB勺程序设计方法；2 )掌握数

5、字信号处理的基本概念、基本理论和基本方法；3 )掌握MATLAB计IIR滤波器；4 )学会用MATLAB寸信号进行分析和处理。2. 本次设计的要求1) 分别用脉冲响应不变法、双向性变换法设计IIR低通、高通、带通、带 阻滤波器滤波器；2) 分别画出其幅频特性、相频特性图；3) IIR滤波器的各项指标：低通：通带截止频率 wc=2n s radk /2,阻带截止频率为8n KHZ!通带衰减pR小于3dB,阻带衰减大于15dB,采样频率20000Hz;高通：通带截止频率为2.5KHZ,通带衰减不大于2dB,阻带上限截止频率为 1.5KHZ，阻带衰减不小于15dB;带通：中心频率为 3 p0=0.5

6、n,通带截止频率 3 p1=0.4n,3 p2=0.6n ;通带 最大衰减a p=3dB；阻带最小衰减a S=15dB;阻带截止频率3 s2=0.7n ;带阻：抽样频率为 10KHZ,在-2dB衰减处边带频率是 1.5KHZ,4KHZ,在 -13dB处边带频率为2KHZ和3KHZ。二、课程设计原理1. 脉冲响应不变法原理脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法，它特 别适合于对滤波器的时域特性有一定要求的场合。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应 序列h(n)模仿模拟滤波器的冲击响应ha(t),使h(n)正好等于ha(t)的采样值， 即h(n)=

7、ha(nT)T为采样周期。如以Ha(s)及H(z)分别表示ha(t)的拉氏变换及h(n)的z变换，即Ha(s)=Lha(t)H(z)=Zh( n)则根据米样序列z变换与模拟信号拉氏变换的关系，可知：米用脉冲响应不变法 将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是以前 讨论的拉氏变换到Z变换的标准变换关系，即首先对Ha(s)作周期延拓，然后再经 过z=e的映射关系映射到Z平面上。脉冲响应不变法映射关系见图2。2. 双向性变换法:脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换z=e的多值对应关系导致的，为了克服这一缺点，设想变换分为两步。1)将

8、整个S平面压缩到S1平面的一条横带里。2)通过标准变换关系将此横带变换到整个 Z平面上去。由此建立S平面与Z平面一s平面一对应的单值关系,图1双线性换法映射关系图双线性换法的主要优点是S平面与Z平面一单值对应，S平面的虚轴(整个j Q ) 对应于Z平面单位圆的一周，S平面的Q =0处对应于Z平面的3 =0处，对应即数字 滤波器的频率响应终 止于折迭频率处，所以双线性变换不存在混迭效应。三、IIR数字滤波器设计内容1.总体方法分析IIR数字滤波器是一种离散时间系统，其系统函数为: 假设 M N,当M N寸,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系 统的级联。IIR数字滤波器

9、的设计实际上是求解滤波器的系数 和.，它是数学 上的一种逼近问题，即在规定意义上(通常采用最小均方误差准则) 去逼近系统 的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就 得到数字滤波器。2.脉冲相应不变法脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应 序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n) 正好等于ha(t)的采样值，满足:h(n)=ha(n!)式中,T是采样周期。如果令H(s)是ha(t)的拉普拉斯变换，H(z)为h(n)的Z变换，利用采样序列 的Z变换与模拟信号的拉普拉斯变换的关系得2 二k(1

10、-1)4 血1 co(X(z)zs-Z Xa(S jkf)=5： Xa S jT kiI k =oa I则可看出，脉冲响应不变法将模拟滤波器的S平面变换成数字滤波器的 Z平面，这个从S到z的变换z=esT是从S平面变换到Z平面的标准变换关系式由(1-1 )式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为(1-2)这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。 正如 采样定理所讨论的，只有当模拟滤波器的频率响应是限带的， 且带限于折叠频率 以内时，即H2 )一k三I j Ha(j0)=O1 C 1 上 T =寸(1-3)才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波

11、器的频率响应，而不产生混叠失真，即(1-4)(e j 驾=T HAa|jT ；但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图7-4所示。这 时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟 滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失 真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。对某一模拟滤波器的单位冲激响应 ha(t)进行采样，采样频率为fs,若使fs 增加，即令采样时间间隔(T=1/fs)减小，则系统频率响应各周期延拓分量之间 相距更远，因

12、而可减小频率响应的混叠效应。3. 双线性变换法脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。 为了克服这一缺点，可以采用非线性频率 压缩方法，将整个频率轴上的频率范围压缩到 -n/Tn/T之间，再用z=esT转 换到Z平面上。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多 值变换性就消除了频谱混叠现象，映射关系如图 4所示。ils平面图4双线性变换的映射关系(1-5)门=-tanT为了将S平面的整个虚轴j Q压缩到S1平面j Q 1轴上的-n /T到n /T段上, 可以通过以下的正切变换实现一 I 2丿式中,T仍是米样间隔。当Q

13、1由-n /T经过0变化到n /T时，Q由-x经过o变化到+x,也即映射了整个j Q轴。将式（1-5）写成.,2 eZ _T JE2 .e/2将此关系解析延拓到整个S平面和S1平面，令j Q =s, j Q 1=s1，贝U得eT/2 e 勿2&T /2 e-S1T /2sT丄.上兰L2 T 1 eT2 1-zs1T 1 z12s 厂再将S1平面通过以下标准变换关系映射到 Z平面 z=es1T从而得到S平面和Z平面的单值映射关系为:(1-6)(1-7)式（1-6）与式（1-7 ）是S平面与Z平面之间的单值映射关系，这种变换都 是两个线性函数之比，因此称为双线性变换式（1-5）与式（1-6）的双线

14、性变换符合映射变换应满足的两点要求。+ tan(1-8)首先,把z=ej ,可得即S平面的虚轴映射到Z平面的单位圆 其次，将s=（T +j Q代入式（1-8 ），得因此由此看出，当（7 0时，|Z|0时，|z|1。也就是说，S平面的左 半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外， S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换 后所得的数字滤波器也一定是稳定的。四、IIR数字滤波器设计过程根据以上IIR数字滤波器设计方法，下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器，其中带通的中心频率为 3 p0=0.5n ,;通带截止频率图5

15、部分滤波器设计指标图示Wp1=0.4n ,3 p2=0.6n ;通带最大衰减a p=3dB;阻带最小衰减a s=15dB;阻带截止频率 3 s2=0.7n1. 设计步骤（以带通为例）(1) 根据任务 ,确定性能指标 :在设计带通滤波器之前 , 首先根据工程实际的需要确定滤波器的技术指标，带通滤波器的阻带边界频率关于中心频率3 p0几何对称，因此 ws1=wpO-(ws2-wp0)=0.3 n通带截止频率 wc1=0.4 n ,wc2=0.6 n ;阻带截止频率 wr1=0.3 n , wr2=0.7 n ; 阻带最小衰减a S=3dB和通带最大衰减a p=15dB；用Q =2/T*tan(w/

16、2)对带通数字滤波器 H(z)的数字边界频率预畸变，得到带 通模拟滤波器 H(s) 的边界频率主要是通带截止频率 3p1, 3p2； 阻带截止频 率3 s1，3 s2 的转换。为了计算简便，对双线性变换法一般 T=2s通带截止频率 wc1=(2/T)*tan(wp1/2)=tan(0.4n /2)=0.7265wc2=(2/T)*tan(wp2/2)=tan(0.6n /2)=1.3764阻带截止频率 wr1=(2/T)*tan(ws1/2)=tan(0.3n /2)=0.5095wr2=(2/T)*tan(ws2/2)=tan(0.7n /2)=1.9626阻带最小衰减a s=3dB和通带最

17、大衰减a p=15dB；(3) 运用低通到带通频率变换公式 入=(Q A2)-( Q 0A2)/(B* Q)将模拟带通滤波器指标转换为模拟低通滤波器指标。B=wc2-wc1=0.6499normwr1=(wr1A2)-(w0A2)/(B*wr1)=2.236normwr2=(wr2A2)-(w0A2)/(B*wr2)=2.236normwc1=(wc1A2)-(w0A2)/(B*wc1)=1normwc2=(wc2A2)-(w0A2)/(B*wc2)=1得出， normwc=1， normwr=2.236模拟低通滤波器指标： normwc=1， normwr=2.236， a p=3dB， a

18、 s=15dB(4) 设计模拟低通原型滤波器。用模拟低通滤波器设计方法得到模拟低通滤波器 的传输函数 Ha(s)； 借助巴特沃斯 (Butterworth) 滤波器、切比雪夫 (Chebyshev) 滤波器、椭圆 (Cauer) 滤波器、贝塞尔 (Bessel) 滤波器等。(5) 调用 lp2bp 函数将模拟低通滤波器转化为模拟带通滤波器。(6) 利用双线性变换法将模拟带通滤波器Ha(s)转换成数字带通滤波器H(z).兀(jQ)恤3)|高通0带阻0Q图6四种数字滤波器2. 程序流程框图开始读入数字滤波器技术指标将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数 N和3db截止

19、频率模拟域频率变换，将g(p)变换成模拟带通滤波器 HS)用双线性变换法将H(s)转换成数字带通滤波器H(z)输入信号后显示相关结果结束图7 程序流程图3. MATLABS 序MATLABg序如下：带通滤波器：clearwp0=0.5*pi;wp1=0.4*pi;wp2=0.6*pi;数字带通滤波器技术指标计算带通滤波器的阻带下截止频率频率预畸变Ap=3;ws2=0.7*pi;As=15;T=2; % ws1=wp0-(ws2-wp0); % wc1=(2/T)*tan(wp1/2);wc2=(2/T)*tan(wp2/2); wr1=(2/T)*tan(ws1/2);wr2=(2/T)*ta

20、n(ws2/2); w0=(2/T)*tan(wp0/2); %将指标转换成归一化模拟低通滤波器的指标%B=wc2-wc1; % 带通滤波器的通带宽度 normwr1=(wr1A2)-(w0A2)/(B*wr1);normwr2=(wr2A2)-(w0A2)/(B*w )； normwc1=(wc1A2)-(w0A2)/(B*wc1); normwc2=(wc2A2)-(w0A2)/(B*wc2)；%if abs(normwr1)abs(normwr2) normwr=abs(normwr2) else normwr=abs(normwr1) end normwc=1； %N=buttord(

21、normwc,normwr,Ap,As,s)； 和 3db 截止频率 bLP,aLP=butter(N,normwc,s)； %G(p)bBP,aBP=lp2bp(bLP,aLP,w0,B)； % 带通滤波器 H(s) b,a=bilinear(bBP,aBP,0.5)； % 通滤波器 H(z) w=linspace (0,2*pi,500)； h=freqz(b,a,w)； subplot(2,1,2)； plot(w,abs(h)； grid on带通到低通的频率变换设计归一化的模拟低通滤波器阶数 N 计算相应的模拟滤波器系统函数模拟域频率变换，将 G(P)变换成模拟用双线性变换法将 H(

22、s) 转换成数字带xlabel(w(rad) ylabel(|H(jw)|) title( 频谱函数 ) subplot(2,2,1); plot(w,20*log10(abs(h); axis(0,2*pi,-120,20); grid on xlabel(w(rad) ylabel(20*lg|H(jw)|(db) title(20*lg|H(jw)|-w)带阻滤波器： %周期方波时域 f1=64; fs=10000; t=0:1/fs:0.1; x1=square(2*pi*f1*t); subplot(321); plot(t,x1); axis(0,0.1,-1.3,1.3); xl

23、abel(t/s); ylabel(x(t); title( 时域谱 ) %周期方波频域 L=1250;X1=fft(x1,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(322); stem(0:L-1),abs(X1);xlabel(w/rad);ylabel(X);title( 幅度谱 );%采用BW和双线性变换法Wp1=0.1*pi,Wp2=0.4*pi;Ws1=0.2*pi,Ws2=0.3*pi;Ap=1;As=50; wp1=2*fs*tan(Wp1/2),wp2=2*fs*tan(Wp2/2); ws1=2*fs*tan(Ws1/2),ws2=2*

24、fs*tan(Ws2/2);B=ws2-ws1,w0=sqrt(ws1*ws2);ws=1;wp1=(B*wp1)/(-wp1A2+w0A2),wp2=(B*wp2)/(-wp2A2+w0A2); wp=max(abs(wp1),abs(wp2);%wp=wp1 wp2,ws=ws1 ws2;N,wc=buttord(wp,ws,Ap,As,s); num,den=butter(N,wc,s);numt,dent=lp2bs(num,den,w0,B); numd,dend=bilinear(numt,dent,fs);w=linspace(0,pi,1024);H=freqz(numd,de

25、nd,w);subplot(324);plot(w/pi,20*log10(abs(H);grid;xlabel(W);ylabel(H/db);title( 幅度响应 );h,t1=impz(numd,dend);subplot(323);plot(t1,h);axis(0,60,min(h),max(h); xlabel(t/s);ylabel(h);title( 单位冲击响应 ); %输出信号 y=filter(numd,dend,x1); subplot(325);plot(t,y); axis(0,0.1,-1.3,1.3); Y=fft(y,L);ws=2*pi*fs; w=(0:

26、L-1)*ws/L; subplot(326); stem(0:L-1),abs(Y);高通滤波器： %周期方波时域 f1=64; fs=10000; t=0:1/fs:0.1; x1=square(2*pi*f1*t); subplot(321); plot(t,x1); axis(0,0.1,-1.3,1.3); xlabel(t/s); ylabel(x(t); title( 时域谱 ) %周期方波频域 L=1250;X1=fft(x1,L);ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(322); stem(0:L-1),abs(X1); xlabel(w/r

27、ad);ylabel(X); title( 幅度谱 );%采用BW和双线性变换法Wp=0.5*pi,Ws=0.1*pi; Ap=1,As=50; wp=2*fs*tan(Wp/2);ws=2*fs*tan(Ws/2); wp1=1/wp,ws1=1/ws;N,wc=buttord(wp1,ws1,Ap,As,s); num,den=butter(N,wc,s); numt,dent=lp2hp(num,den,1); numd,dend=bilinear(numt,dent,fs); w=linspace(0,pi,1024);H=freqz(numd,dend,w); subplot(324

28、); plot(w/pi,abs(H); xlabel(W); ylabel(H);title( 幅度响应 ); h,t1=impz(numd,dend); subplot(323);plot(t1,h); axis(0,max(t1),min(h),max(h); xlabel(t/s);ylabel(h);title( 单位冲击响应 ); %输出信号 y=filter(numd,dend,x1); subplot(325);plot(t,y);axis(0,0.1,-1.3,1.3); Y=fft(y,L);w=(0:L-1)*ws/L; subplot(326); stem(0:L-1)

29、,abs(Y);低通滤波器： %周期方波时域 f1=64;fs=10000;t=0:1/fs:0.1; x1=square(2*pi*f1*t);subplot(321); plot(t,x1);axis(0,0.1,-1.3,1.3); xlabel(t/s);ylabel(x(t);title( 时域谱 ) %周期方波频域L=1250;X1=fft(x1,L); ws=2*pi*fs;w1=(0:L-1)*ws/pi; subplot(322);stem(0:L-1),abs(X1);xlabel(w/rad); ylabel(X); title( 幅度谱 );%采用BW和双线性变换法Wp

30、=0.1*pi;Ws=0.5*pi;Ap=1;As=50; wp=2*fs*tan(Wp/2);ws=2*fs*tan(Ws/2); N,wc=buttord(wp,ws,Ap,As,s); num,den=butter(N,wc,s);numd,dend=bilinear(num,den,fs); w=linspace(0,pi,1024);H=freqz(numd,dend,w);subplot(324); plot(w/pi,abs(H); xlabel(W); ylabel(H); title( 幅度响应 ); h,t1=impz(numd,dend); subplot(323);pl

31、ot(t1,h); axis(0,max(t1),min(h),max(h); xlabel(t/s);ylabel(h);title( 单位冲击响应 ); %输出信号 y=filter(numd,dend,x1); subplot(325); plot(t,y);axis(0,0.1,-1.3,1.3);Y=fft(y,L); ws=2*pi*fs; w=(0:L-1)*ws/pi; subplot(326); stem(0:L-1),abs(Y);4. 调试分析过程描述MATLAB 是矩阵实验室（ Matrix Laboratory ）之意。除具备卓越的数值计算 能力外，同时它还提供了专业

32、水平的符号计算， 文字处理， 可视化建模仿真和实 时控制等功能。它的指令表达式与数学，工程中常用的形式十分相似，故用MATLAB 来解算问题要比用C,FORTRAN语言完成相同的事情简捷得多。用MATLA进行模拟原型的数字滤波器的设计，一般步骤如下 ：（1）按一定规则将给出的数字滤波器的技术指标转换成模拟低通滤波器的技术 指标；（2） 根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数N和固有频 率Wn根据选用的模拟低通滤波器的类型可分别用：buttord,cheblord ， cheb2ord ，ellipord 等函数；（3）运用最小阶数N产生模拟滤波器原型，模拟低通滤波器的创建函数

33、有:buttap,cheblap,cheb2ap,ellipap,besselap等；（4）运用固有频率Wr把模拟低通滤波器原型转换成模拟低通、高通、带通、带阻滤波器，可分别用函数 lp2lp,lp2hp,lp2bp,lp2bs ；（ 5）运用冲激响应不变法或双线性变换法把模拟滤波器转换成数字滤波器，分别用函数impinva和bilinear 来实现。唾IM工EE图8带通滤波器结果图时嘛a1pm* 2 K1, 0M4 -0D10.02 OjOJ OW 0j05006 O.O? 0j0& de 01图9带阻滤波器结果图1.5-iI 一一1IIIII00.1 Q2 0.3040.50.6070.B09120图10高通滤波器结图11低通滤波器结程序运行结果：normwr=2.2361由设计流程计算得normwr=2.236与运行结果相同。低通原型的每一个边界频率都映射为带通滤波器两个相应的边界频率。根据通带截至频率和阻带截至频率与频谱函数曲线比较，