用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器

1、课 程 设 计课程设计名称：用双线性变换法设计原型低通为 椭圆型的数字IIR带通滤波器 专 业 班 级 ： xxxxxx 学 生 姓 名 ： xxx 学 号 ： xxxxxxxxxxxx 指 导 教 师 ： xxx 课程设计时间： 2014.6.16至2014.6.20 电子信息工程 专业课程设计任务书学生姓名XXX专业班级XXXX 学号XXXXXX题 目用双线性变换法设计原型低通为椭圆型的数字IIR低通滤波器课题性质其他课题来源自拟指导教师XX同组姓名主要内容用双线性变换法设计原型低通为椭圆型的数字IIR低通滤波器，要求通带边界频率为400Hz，阻带边界频率分别为500Hz，通带最大衰减1d

2、B，阻带最小衰减40dB，抽样频率为2000Hz，用MATLAB画出幅频特性，画出并分析滤波器传输函数的零极点；信号经过该滤波器，其中300Hz，600Hz，滤波器的输出是什么？用Matlab验证你的结论并给出的图形。任务要求1、掌握用双线性变换法设计原型低通为椭圆型的数字IIR低通滤波器的原理和设计方法。2、求出所设计滤波器的Z变换。3、用MATLAB画出幅频特性图。4、验证所设计的滤波器。参考文献1、程佩青著，数字信号处理教程，清华大学出版社，20012、Sanjit K. Mitra著，孙洪，余翔宇译，数字信号处理实验指导书（MATLAB版），电子工业出版社，2005年1月3、郭仕剑等，

3、MATLAB 7.x数字信号处理，人民邮电出版社，2006年4、胡广书，数字信号处理 理论算法与实现，清华大学出版社，2003年审查意见指导教师签字：教研室主任签字： 2014 年06 月 12 日 说明：本表由指导教师填写，由教研室主任审核后下达给选题学生，装订在设计（论文）首页一需求分析和技术要求 1、需求分析滤波器从广义上来说对特定的频点或频点以外的频率进行有效滤波的电路，这种电路保留输入信号中的有用信息，滤除不需要的信息，从而达到信号的检测、提取、识别等不同的目的。如果处理的信号是时域离散信号，那么相应的处理系统就称为数字滤波器，由于在实际工作中被处理的信号都是幅度量化的数字信号，因此

4、，数字滤波器实际上是用有限精度的算法实现一个线性时不变的时域离散系统。数字滤波器的种类很多，分类方法也不同，可以从功能上分类，也可以从实现方法上及设计方法上来分类等等。滤波器在功能上总的可分为四类，即低通(LP)、高通(HP)、带通(BP)、带阻(BS)滤波器等，从实现方法上，由有限长冲激响应的数字滤波器被称为FIR滤波器，具有无限长冲激响应的数字滤波器称作IIR滤波器。椭圆滤波器（Elliptic filter）又称考尔滤波器（Causer filter），是在通带和阻带等波纹的一种滤波器。椭圆滤波器相比其他类型的滤波器，在阶数相同的条件下有着最小的通带和阻带波动。它在通带和阻带的波动相同，

5、这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及阻带平坦、通带等波纹的切比雪夫滤波器。利用双线性变换法将模拟传输信号Ha(s)变换为数字传输函数G(z)，从而是z域的数字传输函数保留s域的模拟传输函数的基本性质。设计成的IIR数字低通滤波器能够去掉信号中不必要的高频成分，降低采样频率，避免频率混淆，去掉高频干扰。2、技术要求 1、掌握用双线性变换法设计原型低通为椭圆型的数字IIR低通滤波器的原理和设计方法。2、求出所设计滤波器的Z变换。3、用MATLAB画出幅频特性图。4、验证所设计的滤波器。二设计原理及设计思路1、设计原理数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行

6、数学处理来达到频域滤波的目的。数字滤波器也是具有一定传输选择特性的数字信号处理装置，其输入、输出均为数字信号，实质上是一个由有限精度算法实现的线性时不变离散系统。IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为： （2.1） （2.2）设计IIR滤

7、波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。本次课程设计采用双线性变换法。双线性变换法设计IIR数字滤波器双线性变换法主要是采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1所示。 图

8、1双线性变换的映射关系为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上，可以通过以下的正切变换实现 (2.3)式中,T仍是采样间隔。当1由-/T经过0变化到/T时，由-经过0变化到+，也即映射了整个j轴。将式（2.3）写成(2.4)将此关系解析延拓到整个S平面和S1平面，令j=s，j1=s1，则得 （2.5）再将S1平面通过以下标准变换关系映射到Z平面z=es1T (2.6) 从而得到S平面和Z平面的单值映射关系为： （2.7） （2.8） 式（2.7）与式（2.8）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（2.6）与式（2.7）的

9、双线性变换符合映射变换应满足的两点要求。首先,把z=ej，可得 （2.9）即S平面的虚轴映射到Z平面的单位圆。其次，将s=+j代入式（2.8），得 (2.10)因此(2.11)由此看出，当<0时，|z|<1；当>0时，|z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。2、设计思路具体实现方法如下：（1）确定数字低通滤波器的技术指标：通带截止频率p、通带衰减Rp、阻带截止频率s、阻带衰减Rs。（2）将数字低通滤波器的技

10、术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为: (2.12)如果采用双线性变换法，边界频率的转换关系为(2.13)（3）按照模拟低通滤波器的技术指标设计模拟低通滤波器。（4）利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。（5）数字低通技术指标为：p=0.4rad， Rp=1dB； s=0.5rad， Rs=40dB（6）模拟低通的技术指标为：通带截止频率p=(2/T)*tan(p/2)=tan(0.392)=0.0068阻带截止频率s=(2/T)*tan(s/2)=tan(0.628)=0.011通带最小衰

11、减Rp=1dB和阻带最大衰减Rs=40dB；归一化截止角频率p=2pi*Fs/Ft； s=2pi*Fs/Ft（7） 利用模拟椭圆滤波器设计数字滤波器。通带截止频率为：p=0.4*pi；阻带截止频率为：s=0.5*pi；通带最大衰减为：Rp=1；阻带最大衰减为：Rs=40；设定周期为1s；模拟低通滤波器的生成：num1,den1 =ellip(N,Rp,Rs,Wn,'s')；最后实现输入输出、幅频特性、相频特性的图形。三程序流程图确定数字低通滤波器的各项性能指标开始数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标按照模拟低通滤波器的技术指标设计模拟低通滤波器利用双线性变换法

12、将模拟滤波器转变为数字低通滤波器利用模拟椭圆滤波器设计数字滤波器结束图2 滤波器设计流程方框图四程序源代码MATLAB 程序如下：clear all close all Coeff_LF = 1; % 低频信号系数 F1(300 Hz)Coeff_HF = 1; % 高频信号系数 F2(600 Hz) Coeff_NOISE = 0; % 噪声系数, 如0(没有噪声), 1, sqrt(2), 3等等. F1 = 300; % HzF2 = 600; % HzFS = 2000; % Hz 抽样频率 % 时间定义EndTime = 1/FS * 1023;Time = 0:1/FS:EndTi

13、me;N = length(Time); % 产生正弦信号 x1(t),x2(t)和 x(t)SinSignal_x1= Coeff_LF*sin(2*pi*F1*Time)% x1(t)SinSignal_x2= Coeff_HF*sin(2*pi*F2*Time)% x2(t)SinSignalSum = SinSignal_x1 + SinSignal_x2+ Coeff_NOISE*randn(1, N);SinSignal = SinSignalSum;% x(t) % 绘制信号 x1(t)figuresubplot(3,1,1)plot(Time, SinSignal_x1)xla

14、bel('时间')xlim(0 0.5)title('时域信号 x1(t)') subplot(3,1,2)plot(Time, SinSignal_x1)xlabel('时间')xlim(0 0.05)title('时域信号 x1(t)') FFTSinSignal_x1 = 2/N*fft(SinSignal_x1)subplot(3,1,3)FnSin_x1 = abs(FFTSinSignal_x1); FW = 2*pi*(Time)/EndTime; plot(FW(2:length(FW)*FS/pi/2, FnSi

15、n_x1(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin_x1(2:length(FW)xlim(-FS/2 FS/2) xlabel('频率') % 绘制 x2(t)figuresubplot(3,1,1)plot(Time, SinSignal_x2)xlabel('时间')xlim(0 0.5)title('时域信号 x(2)') subplot(3,1,2)plot(Time, SinSignal_x2)xlabel('时间')xlim(0 0.05

16、)title('时域信号 x2(t)') FFTSinSignal_x2 = 2/N*fft(SinSignal_x2)subplot(3,1,3)FnSin_x2 = abs(FFTSinSignal_x2); plot(FW(2:length(FW)*FS/pi/2, FnSin_x2(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin_x2(2:length(FW)xlim(-FS/2 FS/2) xlabel('频率') % 绘制信号 x(t)figuresubplot(3,1,1)p

17、lot(Time, SinSignal)xlabel('时间')xlim(0 0.5)title('时域信号') subplot(3,1,2)plot(Time, SinSignal)xlabel('时间')xlim(0 0.05)title('时域信号 x(t)') FFTSinSignal = 2/N*fft(SinSignal)subplot(3,1,3)FnSin = abs(FFTSinSignal); plot(FW(2:length(FW)*FS/pi/2, FnSin(2:length(FW)hold onplot

18、(FW(2:length(FW)-2*pi)*FS/pi/2, FnSin(2:length(FW)xlim(-FS/2 FS/2) xlabel('频率') %-设计低通滤波器 - Rp = 1; % 通带衰减 /dBRs = 40; % 阻带衰减 /dBOmegaP1_1=400; % 通带边界频率OmegaS1_1=500; % 阻带边界频率FS=2000; % 抽样频率 Wp1=2*pi*OmegaP1_1/FS; % 将模拟频率变为数字角频率Ws1=2*pi*OmegaS1_1/FS; % 将模拟频率变为数字角频率OmegaP1=tan(Wp1/2); % 非线性变换

19、OmegaS1=tan(Ws1/2); % 非线性变换Eta_P=OmegaP1/OmegaP1; % 归一化Eta_S=OmegaS1/OmegaP1; % 归一化 Lemta_P_EquivalentLowPass=1/Eta_P; % 变为等效低通参数Lemta_S_EquivalentLowPass=1/Eta_S; % 变为等效低通参数 % 估计滤波器阶数N, Wn = Ellipord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass,Rp, Rs,'s'); % 时间滤波器num1,den1 = ellip(

20、N,Rp,Rs,Wn,'s');num2,den2=lp2lp(num1,den1,OmegaP1);num,den=bilinear(num2,den2,0.5); % 计算增益响应w = 0:pi/255:pi;h = freqz(num,den,w);g = 20*log10(abs(h); % 绘制增益响应figure;plot(w/pi,g);grid axis(0 1 -60 5);xlabel('omega /pi'); ylabel('增益 /dB');title('椭圆低通滤波器的增益响应'); % 绘制零点和极

21、点figure;z,p,k=tf2zp(num,den);zplane(z,p);title('零点和极点') figuresubplot(2,1,1)YLP = filter(num,den, SinSignal);plot(Time, YLP)xlim(0 0.05)xlabel('时间')title('低通滤波后的时域信号') FFTSinSignalSum = 2 / N * fft(YLP)subplot(2,1,2)FnSinSum = abs(FFTSinSignalSum); plot(FW(2:length(FW)*FS/pi/2, FnSinSum(2:length(FW)hold onplot(FW(2:length(FW)-2*pi)*FS/pi/2, FnSinSum(2:length(FW)xlim(-FS/2 FS/2) ylim(0 1)xlabel('频率') %return%-End Program-五仿真结果图 六 参考资料1程佩青著.数字信号处理教程