用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器.doc

课 程 设 计课程设计名称：数字信号处理课程设计专 业 班 级 ： 电信1108班 学 生 姓 名 ： 学 号 ： 指 导 教 师 ： 乔丽红 课程设计时间： 20140616-20140618 电子信息工程专业课程设计任务书学生姓名专业班级学号题 目用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器课题性质其他课题来源自拟课题指导教师乔丽红同组姓名主要内容用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器，要求通带边界频率为400Hz，500Hz，阻带边界频率分别为350Hz，550Hz，通带最大衰减1dB，阻带最小衰减40dB，抽样频率为2000Hz，用MATLAB画出幅频特性，画出并分析滤波器传输函数的零极点；信号经过该滤波器，其中450Hz，600Hz，滤波器的输出是什么？用Matlab验证你的结论并给出的图形。任务要求1、掌握用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器的原理和设计方法。2、求出所设计滤波器的Z变换。3、用MATLAB画出幅频特性图。4、验证所设计的滤波器。参考文献1、程佩青著，数字信号处理教程，清华大学出版社，20012、Sanjit K. Mitra著，孙洪，余翔宇译，数字信号处理实验指导书（MATLAB版），电子工业出版社，2005年1月3、郭仕剑等，MATLAB 7.x数字信号处理，人民邮电出版社，2006年4、胡广书，数字信号处理 理论算法与实现，清华大学出版社，2003年审查意见指导教师签字：教研室主任签字： 年 月 日 一需求分析和设计内容 1、数字滤波器简介当前我们正处于数字化时代，数字信号处理技术受到了人们的广泛关注，其理论及算法随着计算机技术和微电子技术的发展得到了飞速的发展，被广泛应用于语音图象处理、数字通信、谱分析、模式识别、自动控制等领域。数字滤波器是数字信号处理中最重要的组成部分之一，几乎出现在所有的数字信号处理系统中。数字滤波是语音和图象处理、模式识别、谱分析等应用中的一个基本处理算法，在数字信号处理中占有极其重要的地位。数字滤波器是数字信号处理中最重要的组成部分之一。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。可以设计系统的频率响应，让它满足一定的要求，从而对通过该系统的信号的某些特定的频率成分进行过滤，这就是滤波器的基本原理。而数字带通滤波器就是让满足要求的某一特定频段通过，其他频率分量全部滤除。IIR数字滤波器采用递归型结构，即结构上带有反馈环路，而且相位特性不好控制，对相位要求较高时，需加相位校准网络。另外，IIR数字滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。而且IIR幅频特性精度很高，不是线性相位的，可以应用于对相位信息不敏感的音频信号上。2、设计内容简介分析要求设计一个通带边界频率为400Hz，500Hz，阻带边界频率分别为350Hz，550Hz，通带最大衰减1dB，阻带最小衰减40dB，抽样频率为2000Hz的数字IIR带通滤波器。并且实现方法要求使用双线性法而且低通原型为巴特沃兹型。用MATLAB画出所设计的数字IIR带通滤波器的幅频特性，画出并分析其传输函数的零极点。求出所设计的滤波器的Z变换，并用给定的信号验证该滤波器的性能，画出给定信号的图像以及输出信号的图像。二设计原理及设计思路1、设计原理2、 设计思路 利用下式（1）预计所求数字滤波器Gbp(z)的数字频率指标，再利用式（2）得到一个等价的模拟带通滤波器的频率指标。 ，（1） ，（2） 采取归一化方法，利用下式（3）和（4）进行频率变换，将模拟带通滤波器的频率指标转换成原型低通滤波器的频率指标。 ，（3） ，（4）利用巴特沃兹滤波器的设计方法设计模拟低通滤波器，在此基础上得到模拟带通滤波器，最后再利用双线性变换法得到所要的数字IIR带通滤波器。所需用到的MATLAB函数如下：N, Wn =buttord( Wp ,Ws ,Rp, Rs,s); num1,den1 = butter(N,Wn,s);num2,den2=lp2bp(num1,den1,OmegaP0,Bw); num,den=bilinear(num2,den2,0.5); ff=filter(num,den,x); h=freqz(num,den,w); zplane(num,den); 最后根据滤波器的幅频特性、零极点分布图以及对给定信号的滤波处理效果来分析所设计的数字IIR带通滤波器的性能。三程序流程图指定所要设计滤波器的参数指标计算所求数字滤波器的数字频率指标归一化，频率变换得到转换成原型低通滤波器的频率指标求得等价带通滤波器的频率指标利用M文件buttord，求得巴特沃兹模拟传输函数的阶数与3dB截止角频率利用M文件butter设计一个指定3dB截止角频率Wn的 N 阶低通传输函数利用M文件bilinear实现双线性变换，得到数字IIR带通滤波器Gbp(z)利用M文件lp2bp实现低通滤波器到带通滤波器的转换，确定模拟带通滤波器Hbp(s)用MATLAB设定输入信号x(t),并画出输入信号的图像用所设计的滤波器对给定信号进行滤波处理，并画出滤波后的信号y(t)的图像，验证其性能将输出信号y（t）的图像与输入信号频率分量x1(t)图像相对照,并画出输入输出信号的频谱图画出滤波器的幅频特性以及零极点分布图。并结合信号处理效果对滤波器性能进行分析四程序源代码% Program % 用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器clc;clear all;Rp = 1; % 通带最大衰减 dBRs = 40; % 阻带最小衰减 dBOmegaP1_1=400; % 通带边界频率OmegaP1_2=500; % 通带边界频率OmegaS1_1=350; % 阻带边界频率OmegaS1_2=550; % 阻带边界频率Fp=2000; % 抽样频率Wp1=2*pi*OmegaP1_1/Fp % 转换成数字信号时的带通截止角频率Wp2=2*pi*OmegaP1_2/Fp % 转换成数字信号时的带通截止角频率Ws1=2*pi*OmegaS1_1/Fp % 转换成数字信号时的带阻截止角频率Ws2=2*pi*OmegaS1_2/Fp % 转换成数字信号时的带阻截止角频率OmegaP1=tan(Wp1/2); % 模拟带通滤波器的截止角频率OmegaP2=tan(Wp2/2); % 模拟带通滤波器的截止角频率OmegaS1=tan(Ws1/2); % 模拟带通滤波器的截止角频率OmegaS2=tan(Ws2/2); % 模拟带通滤波器的截止角频率OmegaP0=sqrt(OmegaP1*OmegaP2);% 通带中心频率Bw=OmegaP2-OmegaP1; % 带通滤波器通带带宽Eta_P0=OmegaP0/Bw; % 归一化Eta_P1=OmegaP1/Bw; % 归一化Eta_P2=OmegaP2/Bw; % 归一化Eta_S1=OmegaS1/Bw; % 归一化Eta_S2=OmegaS2/Bw; % 归一化%模拟低通滤波器的指标Lemta_P_EquivalentLowPass=(Eta_P22-Eta_P02)/Eta_P2; % 归一化通带截止角频率Lemta_S1_EquivalentLowPass=-(Eta_S12-Eta_P02)/Eta_S1; % 归一化阻带截止角频率Lemta_S2_EquivalentLowPass=(Eta_S22-Eta_P02)/Eta_S2; % 归一化阻带截止角频率Lemta_S_EquivalentLowPass=min(Lemta_S1_EquivalentLowPass,Lemta_S2_EquivalentLowPass);% 取归一化阻带截止频率中的较小值作为阻带截止角频率N, Wn = buttord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass, Rp, Rs,s); % 确定低通滤波器的阶数与3dB截止频率num1,den1 = butter(N,Wn,s); % 确定低通滤波器原型的传输函数Hlp(s)num2,den2=lp2bp(num1,den1,OmegaP0,Bw); % 确定带通滤波器传输函数Hbp(s)num,den=bilinear(num2,den2,0.5); % 数字IIR带通滤波器的传输函数Gbp(z)% Compute the gain responsew = 0:pi/255:pi;h = freqz(num,den,w); % 数字带通滤波器的频率响应g = 20*log10(abs(h); % 转化为 dB% Plot the gain response%输入信号 x(t)%t=0:pi/(255*255):8/900;t=(0:9999)/2000;x1=sin(2*pi*450*t);x2=sin(2*pi*600*t);x=x1+x2;%subplot(1,3,1),figure;plot(x1);axis(100 200 -1 1);grid;title(x1(t);%subplot(1,3,2),figure;plot(x2);axis(100 200 -1 1);grid;title(x2(t);%subplot(1,3,3),figure;plot(x);axis(100 200 -2 2);grid;title(x(t);ff=filter(num,den,x); % 滤波器滤波figure;plot(ff),axis(100 200 -1 1);grid;title(滤波处理后输出信号y(t)figure;plot(x1);axis(100 200 -1 1);hold on;plot(ff,r),axis(100 200 -1 1);title(蓝色线为x1(t)图像，红色线为滤波器输出信号y(t)图像);Df1=ff-x1;figure;plot(Df1);axis(100 200 -0.2 0.2);grid;title(将输出信号y(t)与输入频率分量x1（t）相减后的效果);% 输入输出信号频率分布图 fs=2000;N2=1024;n=0:N2-1;t=n/fs;Y1=fft(x,N2);Y2=fft(ff,N2);mag1=abs(Y1);mag2=abs(Y2);f=n*fs/N2;figure;subplot(2,1,1);plot(f(1:N2/2),mag1(1:N2/2);title(输入信号的频谱图);xlabel(频率/HZ);ylabel(振幅);grid on;subplot(2,1,2);plot(f(1:N2/2),mag2(1:N2/2); title(输出信号的频谱图);xlabel(频率/HZ);ylabel(振幅);grid on; figure;plot(w/pi,g),grid; % 滤波器的幅频特性axis(0 1 -60 5);xlabel(omega /pi); ylabel(Gain in dB);title(Gain Response of a Butterworth Bandpass Filter);% 求IIR数字带通滤波器传输函数的零极点figure;zplane(num,den),grid;title(滤波器传输函数零极点);五仿真结果图1、仿真过程截图 图（1） 原输入信号频率分量为450Hz的x1(t)=sin(2*pi*f1*t)的图像 图（2） 原输入信号频率分量为600Hz的x2(t)=sin(2*pi*f2*t)的图像 图（3） 原输入信号x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t)的图像 图（4） 经过滤波器处理后的输出信号y(t)的图像 图（5） 输出信号y(t)与450Hz频率分量x1(t)的对比图像 图（6） 输出信号y（t）与450Hz频率分量x1(t)相减输出 图（7） 输入输出信号频率分布图 图（8） 滤波器幅频特性 图（9） 滤波器传输函数的零极点分布图 2、仿真效果以及滤波器性能分析 由图像（5）输出信号y(t)与原信号频率分量x1(t)的对比图像以及图像（6）可以明显地看出滤波后的信号图像只剩下了频率为450Hz的原输入信号分量x1(t)；再结合图像（7）可以看出被滤除的信号为频率分量为600Hz的x2(t)，即频率为600Hz的原输入信号分量x2(t)已经被基本滤除。由此可看出本数字IIR带通滤波器在满足设计要求下能够正常工作。而且由程序算出滤波器的数字频率指标如下： Wpl = 1.2566 = 0.4,Wph =1.5708=0.5, Wsl=1.0996=0.35, Wsh=1.7279=0.55。与图（8）滤波器的幅频特性图像相对照并对比滤波器的阻带衰减与通带衰减可以知道本次所设计的滤波器应该说基本达到了设计指标。 由滤波器的零极点分布图（9）可以看出：该滤波器是线性滤波器，他们的零点都是关于单位圆镜像对称，并且关于实轴对称。混合相位传输函数，其零点既在单位圆外也在单位圆内。关于极点，其描述的是滤波器的稳定性。本滤波器传输函数为因果序列，而且其极点全部位于单位圆内，因此本滤波器BIBO稳定。六 参考资料1、杨大柱，MATLAB在IIR滤波器设计中的应用，现代机械期刊，2006年第五期2、Sanjit K. Mitra著，孙洪，余翔宇译，数字信号处理实验指导书（MATLAB版），电子工业出版社，2005年1月3、Sanjit K. Mitra著，孙洪译，数字信号处理-基于计算机的方法（第三版），电子工业出版社，2006年6月七设计心得 原先在学习数字信号处理时，看着课本上各种各样的分析公式，以及各种各样的滤波器设计方