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

《IIR数字滤波器课程设计汇总.pdf》由会员分享，可在线阅读，更多相关《IIR数字滤波器课程设计汇总.pdf（24页珍藏版）》请在三一文库上搜索。

1、1 数字信号处理 课 程 设 计 报 告 基于 MATLAB 的 IIR 数字滤波器设计 专业班级 : 电信工程 1302 班 学号: 311308000626 学生姓名 : 王海龙 指导教师 : 王科平 2016年 7 月 2 目录 摘要3 一、课程设计任务及要求4 1. 本次设计的目的 .4 2. 本次设计的要求 .4 二、课程设计原理 4 1. 脉冲响应不变法原理 .4 2. 双向性变换法原理 .5 三、IIR 数字滤波器设计内容 5 1. 总体方法分析 .5 2. 脉冲相应不变法 .6 3. 双线性变换法 .7 四、IIR 数字滤波器设计过程 .9 1. 设计步骤 .9 2. 程序流程

2、框图 .11 3.MATLAB 程序.11 4. 调试分析过程描述 .19 5 结果分析 19 五、结论 22 六、参考文献 23 3 摘要 在当今社会， 数字信号处理技术飞速发展， 它不但自成一门学科， 更是以不 同的方式影响和渗透到其他学科的研究中，它变得与我们的生活联系越来越紧 密，不断改变着我们的生产生活方式，因此受到人们越来越多的关注。 数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一 种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运 算转变为输出的数字序列。 描述离散系统输出与输入关系的卷积和差分方程只是 给数字信号滤波器提供运算规则，使其按

3、照这个规则完成对输入数据的处理。时 域离散系统的频域特性 :Y(e jw )=X(e jw )H(e jw) , 其中、分别是数字 滤波器的输出序列和输入序列的频域特性（或称为频谱特性）,H(e jw )是数字滤 波器的单位取样响应的频谱，又称为数字滤波器的频域响应。 数字滤波器是具有一定传输选择特性的数字信号处理装置, 其输入、输出均 为数字信号 , 实质上是一个由有限精度算法实现的线性时不变离散系统。IIR 数 字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现, 其差 分方程为：y(n)= N a 0i i x(n-i)+ N i 1 bi y(n-i) 系统函数为： H（

4、z）=( M r br 0 Z -r )/( 1+ N a 0k k Z -k ) 设计IIR 滤波器的任务就是寻求一个物理上可实现的系统函数H(z) ，使其频 率响应 H(z) 满足所希望得到的频域指标。 本次课程设计分别用脉冲响应不变法、双向性变换法设计IIR 低通、高通、 带通、带阻滤波器滤波器。 并在MATLAB 环境下实现了 IIR 数字滤波器的设计和仿 真。其主要内容概括为： 首先对滤波器的原理和设计进行了介绍；接着描述了 IIR 数字滤波器的基本概念，其中包括系统的描述、系统的传递函数、系统的模型； 接着简单介绍 MATLAB ，并对数字滤波器在 MATLAB 环境下如何实现进行

5、了介绍； 重 点描述了 IIR 数字滤波器的设计过程，最后对IIR 滤波器进行仿真。 关键词：数字滤波器频域特性脉冲响应双向性变换法 MATLAB 4 一、课程设计任务及要求 1. 本次设计的目的 1 ）学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法； 2 ）掌握数字信号处理的基本概念、基本理论和基本方法； 3 ）掌握 MATLAB 设计 IIR 滤波器； 4 ）学会用 MATLAB 对信号进行分析和处理。 2. 本次设计的要求 1）分别用脉冲响应不变法、双向性变换法设计IIR 低通、高通、带通、带 阻滤波器滤波器； 2）分别画出其幅频特性、相频特性图； 3）IIR 滤波器的各

6、项指标 : 低通：通带截止频率wc=2s radk /2,阻带截止频率为 8KHZ 通带衰减 pR 小于 3dB，阻带衰减大于15dB,采样频率 20000Hz; 高通：通带截止频率为 2.5KHZ,通带衰减不大于 2dB,阻带上限截止频率为 1.5KHZ，阻带衰减不小于15dB; 带通：中心频率为 p0=0.5,通带截止频率 p1=0.4,p2=0.6;通带 最大衰减 p=3dB;阻带最小衰减 s=15dB;阻带截止频率 s2=0.7; 带阻：抽样频率为10KHZ,在-2dB 衰减处边带频率是1.5KHZ,4KHZ, 在 -13dB 处边带频率为 2KHZ 和 3KHZ 。 二、课程设计原理

7、 1. 脉冲响应不变法原理 脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法，它特 别适合于对滤波器的时域特性有一定要求的场合。 脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应 序列h(n) 模仿模拟滤波器的冲击响应ha(t), 使h(n) 正好等于 ha(t) 的采样值， 即h(n)=ha(nT)T 为采样周期。 如以Ha(s) 及H(z) 分别表示 ha(t) 的拉氏变换及 h(n) 的z变换，即 5 Ha(s)=Lha(t) H(z)=Zh(n) 则根据采样序列 z变换与模拟信号拉氏变换的关系，可知：采用脉冲响应不变法 将模拟滤波器变换为数字滤波器时，它所完

8、成的 S平面到 Z平面的变换，正是以前 讨论的拉氏变换到 Z变换的标准变换关系， 即首先对 Ha(s) 作周期延拓，然后再经 过z=e的映射关系映射到 Z平面上。 脉冲响应不变法映射关系见图2。 2. 双向性变换法： 脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到 Z平面的 标准变换 z=e的多值对应关系导致的 , 为了克服这一缺点，设想变换分为两步。 1) 将整个 S平面压缩到 S1平面的一条横带里。 2)通过标准变换关系将此横带变换到整个Z平面上去。由此建立 S平面与 Z平面一 一对应的单值关系，消除多值性，也就消除了混淆现象。 o 11 Z平面 jImz Rez / T j

9、 1 1 / T S 1平面 S平面 j oo 图1 双线性换法映射关系图 双线性换法的主要优点是S平面与 Z平面一单值对应， S平面的虚轴 ( 整个j ) 对应于 Z平面单位圆的一周， S平面的 =0处对应于 Z平面的 =0处，对应即数字 滤波器的频率响应终止于折迭频率处，所以双线性变换不存在混迭效应。 三、IIR 数字滤波器设计内容 1. 总体方法分析 IIR 数字滤波器是一种离散时间系统，其系统函数为： 6 假设M N， 当M N 时, 系统函数可以看作一个 IIR 的子系统和一个 (M-N)的FIR子系 统的级联。 IIR 数字滤波器的设计实际上是求解滤波器的系数和，它是数学 上的一种

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

11、的关系得 (1-1) 则可看出，脉冲响应不变法将模拟滤波器的S 平面变换成数字滤波器的Z 平面，这个从 s 到 z 的变换 z=e sT 是从 S平面变换到 Z平面的标准变换关系式。 图 2 脉冲响应不变法的映射关系 由（1-1）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为 (1-2) 这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如 采样定理所讨论的， 只有当模拟滤波器的频率响应是限带的，且带限于折叠频率 以内时，即 0)( jHa (1-3) 才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应， 而不产生混叠失真，即 k T jsX T jksX

12、 T zX k as k a ez sT 21 )( 1 )( j 3/ T / T 3/ T / T o o 11 jIm z Rez Z平面 S平面 T k jH T eH k a j 21 )( 2 | s T 7 |w|0 时，| z|1 。也就是说， S 平面的左 半平面映射到 Z 平面的单位圆内， S平面的右半平面映射到Z平面的单位圆外， S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换 后所得的数字滤波器也一定是稳定的。 四、IIR 数字滤波器设计过程 根据以上 IIR 数字滤波器设计方法，下面运用双线性变换法基于MATLAB 设计一个 IIR 带通滤波器，

13、其中带通的中心频率为p0=0.5，;通带截止频率 图 5 部分滤波器设计指标图示 wp1=0.4,p2=0.6;通带最大衰减 p=3dB;阻带最小衰减 s=15dB;阻带截止 频率s2=0.7 1. 设计步骤 （以带通为例） j T j T z 2 2 2 2 2 2 2 2 | T T z 10 (1) 根据任务 , 确定性能指标 : 在设计带通滤波器之前, 首先根据工程实际的需要 确定滤波器的技术指标， 带通滤波器的阻带边界频率关于中心频率p0 几何对称，因此ws1=wp0- (ws2-wp0)=0.3 通带截止频率 wc1=0.4,wc2=0.6 ; 阻带截止频率 wr1=0.3， wr

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

15、2/2)=tan(0.7/2)=1.9626 阻带最小衰减 s=3dB和通带最大衰减 p=15dB; (3) 运用低通到带通频率变换公式=( 2)-( 02)/(B*) 将模拟带通滤 波器指标转换为模拟低通滤波器指标。 B=wc2-wc1=0.6499 normwr1=(wr12)-(w02)/(B*wr1)=2.236 normwr2=(wr22)-(w02)/(B*wr2)=2.236 normwc1=(wc12)-(w02)/(B*wc1)=1 normwc2=(wc22)-(w02)/(B*wc2)=1 得出， normwc=1 ，normwr=2.236 模拟低通滤波器指标： nor

16、mwc=1 ，normwr=2.236，p=3dB ，s=15dB (4) 设计模拟低通原型滤波器。用模拟低通滤波器设计方法得到模拟低通滤波器 的传输函数 Ha(s); 借助巴特沃斯 (Butterworth)滤波器、切比雪夫 (Chebyshev) 滤波器、椭圆 (Cauer) 滤波器、贝塞尔 (Bessel) 滤波器等。 (5) 调用 lp2bp 函数将模拟低通滤波器转化为模拟带通滤波器。 (6) 利用双线性变换法将模拟带通滤波器Ha(s) 转换成数字带通滤波器H(z). 11 图 6 四种数字滤波器 2. 程序流程框图 开始 读入数字滤波器技术指标 将指标转换成归一化模拟低通滤波器的指标

17、 设计归一化的模拟低通滤波器阶数N和 3db 截止频率 模拟域频率变换，将 G(P)变换成模拟带通滤波器H(s) 用双线性变换法将 H(s) 转换成数字带通滤波器H(z) 输入信号后显示相关结果 图 7 程序流程图 3.MATLAB 程序 MATLAB 程序如下 : 带通滤波器： clear 结束 12 wp0=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(wp

18、2/2); wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2); w0=(2/T)*tan(wp0/2); %频率预畸变 B=wc2-wc1; %带通滤波器的通带宽度 normwr1=(wr12)-(w02)/(B*wr1); normwr2=(wr22)-(w02)/(B*wr2); normwc1=(wc12)-(w02)/(B*wc1); normwc2=(wc22)-(w02)/(B*wc2); %带通到低通的频率变换 if abs(normwr1)abs(normwr2) normwr=abs(normwr2) else normwr=abs(norm

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

20、qz(b,a,w); subplot(2,1,2); plot(w,abs(h); grid on 13 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); subpl

21、ot(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); 14 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*

22、tan(Wp1/2),wp2=2*fs*tan(Wp2/2); ws1=2*fs*tan(Ws1/2),ws2=2*fs*tan(Ws2/2); B=ws2-ws1,w0=sqrt(ws1*ws2); ws=1; wp1=(B*wp1)/(-wp12+w02),wp2=(B*wp2)/(-wp22+w02); 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=

23、bilinear(numt,dent,fs); w=linspace(0,pi,1024); H=freqz(numd,dend,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); 15 axis(0,60,min(h),max(h); xlabel(t/s); ylabel(h); title(单位冲击响应 ); % 输出信号 y=filter(numd,dend,x1)

24、; 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/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=ff

25、t(x1,L); 16 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 和双线性变换法 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)

26、; numd,dend=bilinear(numt,dent,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); plot(t1,h); axis(0,max(t1),min(h),max(h); xlabel(t/s); ylabel(h); 17 title(单位冲击响应 ); % 输出信号 y=filter(numd,dend,x1);

27、 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),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

28、*pi*fs; w1=(0:L-1)*ws/pi; subplot(322); stem(0:L-1),abs(X1); 18 xlabel(w/rad); ylabel(X); title(幅度谱 ); % 采用 BW 和双线性变换法 Wp=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=f

29、reqz(numd,dend,w); subplot(324); 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); 19 Y=fft(y,L); ws

30、=2*pi*fs; w=(0:L-1)*ws/pi; subplot(326); stem(0:L-1),abs(Y); 4. 调试分析过程描述 MATLAB 是矩阵实验室（ Matrix Laboratory）之意。除具备卓越的数值计算 能力外，同时它还提供了专业水平的符号计算，文字处理， 可视化建模仿真和实 时控制等功能。它的指令表达式与数学 , 工程中常用的形式十分相似, 故用MATLAB 来解算问题要比用 C,FORTRAN等语言完成相同的事情简捷得多。 用MATLAB 进行模拟原型的数字滤波器的设计，一般步骤如下: （1）按一定规则将给出的数字滤波器的技术指标转换成模拟低通滤波器的技

31、术 指标； （2）根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数 N 和固有频 率Wn ，根据选用的模拟低通滤波器的类型可分别用: buttord,cheblord， cheb2ord，ellipord等函数； （3）运用最小阶数 N产生模拟滤波器原型，模拟低通滤波器的创建函数 有:buttap,cheblap,cheb2ap,ellipap,besselap等； （4）运用固有频率 Wn 把模拟低通滤波器原型转换成模拟低通、高通、带通、带 阻滤波器，可分别用函数lp2lp,lp2hp,lp2bp,lp2bs； （5）运用冲激响应不变法或双线性变换法把模拟滤波器转换成数字滤波器，分

32、别用函数 impinva 和bilinear来实现。 20 图 8 带通滤波器结果图 21 图 9 带阻滤波器结果图 图 10 高通滤波器结 22 图 11 低通滤波器结 程序运行结果： normwr=2.2361 由设计流程计算得normwr=2.236 与运行结果相同。 低通原型的每一个边界频率都映射为带通滤波器两个相应的边界频率。根据 通带截至频率和阻带截至频率与频谱函数曲线比较，满足设计要求。 五总结 本次通过这个数字信号处理课程的课程设计，使得自己的能力和知识得到了 很大的提升。 不但对于上学期学到的数字信号处理内容有了回顾和复习，而且自 己也真正的在我们的实验中体会到了理论知识对我

33、们的实验的指导意义，在自己 的切身的实践动手的过程中， 自己不仅对于上学期学到的内容中设计带通数字滤 波器的整个过程流程的思路有了一定的了解，同时我也通过自己的实际动手操作 发现原来自己的学习效果是如此的差，自己对于我们学到的知识从来都是一知半 解，自己学习不动脑筋，不去思考，尤其是在自己的课程设计遇到差错的时候， 自己没有一点耐心， 十分烦躁， 这一切的原因都是自己学习的不够深入，如果自 23 己学习过程中注重思考， 那么遇到的问题也就不会显得那么突然和无可奈何。同 时在本次的实验中对于我们学到的双线性变换法，巴特沃斯设计模拟滤波器的运 用，也算是比较熟悉了。 同时在对数字带通滤波器的设计过

34、程中，自己也重新学习了MATLAB 的运行 环境，自己也通过复习课本的知识，对于MATLAB 了解的更加深入，同时在上网 学习过程中学习初步掌握了MATLAB 语言在数字信号处理中一些基本库函数的调 用和编写基本程序等的应用； 了解和熟悉了在滤波器设计中， 我们需要具备的设 计思路；掌握到滤波器设计的一般原理， 对滤波器的工作情况以及最终的设计的 结果有了一定的了解； 同时自己也掌握到了数字高通滤波器设计的一般步骤；加 深了对滤波器设计中产生误差的原因以及双线性变换法优缺点的理解和认识。这 些对于自己知识和实践能力的提高有很大的作用。 在这次课程设计中， 自己主要是负责寻找资料， 以及实验报告

35、的排版编写工 作，主要是关于在本次课程设计中需要用到的工作原理和设计思路的流程设计， 自己和自己的队友分工明确， 各取所长。自己也在本次的课程设计中实实在在锻 炼了自己，同时自己也好好补习了了自己所欠缺的MATLAB 语言知识，对于自己 的编程能力有了提高， 同时在复习到自己上学期学习的数字信号处理过程中，使 理论联系了实际，巩固并深化了对课本基本知识的认识和理解，使理论得以升华。 也算是一次不小的收获吧。 六、参考文献 1高西全 , 丁美玉 . 数字信号处理 . 第 1 版. 北京：西安电子科技大学出版社, 2008 2郭仕剑，王宝顺，贺志国，杨可心.MATLAB7.X数字信号处理 . 人民邮电出 版社, 2006 3李正周 .Matlab 数字信号处理与应用工程 . 北京：清华大学出版社，2008 4刘舒帆 . 数字信号处理实验 （MATLAB ） 版. 北京：西安电子科技大学出版社 . 24 5张圣勤 .MATLAB 7.0 使用教程 . 北京：机械工业出版社， 2008