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

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

2、画出幅频特性,画出并分析滤波器传输函数的零极点；信号x(t)xi(t)X2(t)sin(2fit)sin(2fzt)经过该滤波器,其中f1300Hz,f2600Hz,滤波器的输出y(t)是什么用Matlab验证你的结论并给出Xi(t),X2(t),x(t),y(t)的图形.任务要求1、掌握用双线性变换法设计原型低通为椭圆型的数字IIR低通滤波器的原理和设计方法.2、求出所设计滤波器的Z变换.3、用MATLAB画出幅频特性图.4、验证所设计的滤波器.奔f文献1、程佩青著,?数字信号处理教程?,清华大学出版社,20012、SanjitK.Mitra著,孙洪,余翔宇译,?数字信号处理实验指导书(MA

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

4、实际工作中被处理的信号都是幅度量化的数字信号,因此,数字滤波器实际上是用有限精度的算法实现一个线性时不变的时域离散系统.数字滤波器的种类很多,分类方法也不同,可以从功能上分类,也可以从实现方法上及设计方法上来分类等等.滤波器在功能上总的可分为四类,即低通(LP)、高通(HP)、带通(BP)、带阻(BS)滤波器等,从实现方法上,由有限长冲激响应的数字滤波器被称为FIR滤波器,具有无限长冲激响应的数字滤波器称作IIR滤波器.椭圆滤波器(日lipticfilter)又称考尔滤波器(Causerfilter),是在通带和阻带等波纹的一种滤波器.椭圆滤波器相比其他类型的滤波器,在阶数相同的条件下有着最小

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

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

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

8、叠现象,映射关系如图1所示图1双线性变换的映射关系为了将S平面的整个虚轴jQ压缩至ijS1平面jQ1轴上的-冗/T到上,可以通过以下的正切变换实现-tanTiT(2.3)式中,T仍是采样间隔.当Q1由-冗/T经过0变化到兀/T时,Q由-00经过0变化到+oo,也即映射了整个j.轴.将式2.3写成将此关系解析延拓到整个ON1T/2j1T/22ee=jiT/2jiT/2Tee(2.4)S平面和S1平面,令jQ=s,jQ1=s1,那么得2eS1T/2sT严SiT/2eSiT/2e0T0T(2.5)再将S1平面通过以下标准变换关系映射到Z平面z=es1T(2.6)从而得到S平面和Z平面的单值映射关系为

9、:21z1(2.(7)1Tsz宜22sT_2sT(2.(8)式2.7与式2.8是S平面与Z平面之间的单值映射关系,这种变换都是两个线性函数之比,因此称为双线性变换式2.6与式2.7的双线性变换符合映射变换应满足的两点要求.首先,把z=e;可得(2.9)21ej.2,s厂jtanjT1ejT2即S平面的虚轴映射到Z平面的单位圆.其次,将s=(r+jQ代入式(2.8),得因此(2.10)(2.11),T由此看出,当(T<0时,|z|<1;当(T>0时,|z|>1.也就是说,S平面的左半平面映射到Z平面的单位圆内,S平面的右半平面映射到Z平面的单位圆外,S平面的虚轴映射到Z平

10、面的单位圆上.因此,稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的.2、设计思路具体实现方法如下:(1)确定数字低通滤波器的技术指标：通带截止频率p、通带衰减Rp>阻带截止频率s、阻带衰减Rs.(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标.如果采用脉冲响不变法,边界频率的转换关系为(2.12)如果采用双线性变换法,边界频率的转换关系为河(2.13)(3)根据模拟低通滤波器的技术指标设计模拟低通滤波器.(4)利用双线性变换法将模拟滤波器Ha(s),从s平面转换到z平面,得到数字低通滤波器系统函数H(z)0(5)数字低通技术指标为:p=0.4兀rad,Rp=1

11、dB;s=0.5兀rad,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通带最小衰减Rp=1dB和阻带最大衰减Rs=40dB;归一化截止角频率p=2pi*Fs/Ft;ws=2pi*Fs/Ft(7)利用模拟椭圆滤波器设计数字滤波器.通带截止频率为：p=0.4*pi；阻带截止频率为：s=0.5*pi;通带最大衰减为：Rp=1;阻带最大衰减为：Rs=40；设定周期为1s；模拟低通滤波器的生成：num1,den1=ellip(N,Rp,Rs,Wn,&

12、#39;s');最后实现输入输出、幅频特性、相频特性的图形.三程序流程图开始确定数字低通滤波器的各项性能指标数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标根据模拟低通滤波器的技术指标设计模拟低通滤波器利用双线性变换法将模拟滤波器转变为数字低通滤波器利用模拟椭圆滤波器设计数字滤波器结束图2滤波器设计流程方框图四程序源代码MATLAB程序如下：clearallcloseallCoeff_LF=1;%低频信号系数F1(300Hz)Coeff_HF=1;%高频信号系数F2(600Hz)Coeff_NOISE=0;%噪声系数,如0(没有噪声),1,sqrt(2),3等等.F1=300;%

13、HzF2=600;%HzFS=2000;%Hz抽样频率%时间定义EndTime=1/FS*1023;Time=0:1/FS:EndTime;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)fi

14、guresubplot(3,1,1)xlabel('时间')xlim(00.5)t田e('时域信号x1(t)')subplot(3,1,2)plot(Time,SinSignal_x1)xlabel('时间')xlim(00.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,F

15、nSin_x1(2:length(FW)holdonplot(FW(2:length(FW)-2*pi)*FS/pi/2,FnSin_x1(2:length(FW)xlim(-FS/2FS/2)xlabel('频率')%绘制x2(t)figuresubplot(3,1,1)plot(Time,SinSignal_x2)xlabel('时间')xlim(00.5)title('时域信号x(2)')subplot(3,1,2)plot(Time,SinSignal_x2)xlabel('时间')xlim(00.05)title(

16、9;时域信号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)holdonplot(FW(2:length(FW)-2*pi)*FS/pi/2,FnSin_x2(2:length(FW)xlim(-FS/2FS/2)xlabel('频率')%绘制信号x(t)figuresubplot(3,1,1)plot(Time,SinSignal)xlabel

17、('时间')xlim(00.5)title('时域信号')subplot(3,1,2)plot(Time,SinSignal)xlabel('时间')xlim(00.05)t田e('时域信号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)holdonplot(FW(2:length(FW)-2*pi)*FS/pi/2,FnSin(2:le

18、ngth(FW)xlim(-FS/2FS/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);%非线性变换OmegaS1=tan(Ws1/2);%非线性变换Eta_P=OmegaP1/OmegaP1;%归一化Eta_S=OmegaS1/OmegaP1

19、;%归一化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(N,Rp,Rs,Wn,'s');num2,den2=lp21P(num1,den1,OmegaP1);num,den=bilinear(num2,den2,0.5);%计算

20、增益响应w=0:pi/255:pi;h=freqz(num,den,w);g=20*log10(abs(h);%绘制增益响应figure;plot(w/pi,g);gridaxis(01-605);xlabel('omega八pi');ylabel('增益/dB');title('椭圆低通滤波器的增益响应');%绘制零点和极点figure;z,p,k=tf2zp(num,den);zplane(z,p);title('零点和极点')figuresubplot(2,1,1)YLP=filter(num,den,SinSignal);

21、plot(Time,YLP)xlim(00.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)holdonplot(FW(2:length(FW)-2*pi)*FS/pi/2,FnSinSum(2:length(FW)xlim(-FS/2FS/2)ylim(01)xlabel('频率')%retur

22、nEndProgram五仿真结果图时域信号x1(t)时间时域信号x1(t)00.0050.010.0150.020.0250.030.0350.040.0450.05时间1CC,ECTiC；0.5-011cJ1-1000-800-600-400-20002004006008001000频率时域信号x(2)时间时域信号x2(t)00.0050.010.0150.020.0250.030.0350.040.0450.05时间10.50-1000-800-600-400-20002004006008001000频率时域信号时间时域信号x(t)91Ln-LiL-1,AJ|JV-!111ii/VjVNI

23、Mf1f1I1,Af乂AAIr/AIfl/Vi/If11j1fa1iiVJJjjj1EIJJrrrrf1rI0-200.0050.010.0150.020.0250.030.0350.040.0450.05时间一LLL51LLI,i1r-J1Jr,0.5-1000-800-600-400-20002004006008001000频率椭圆低通滤波器的增益响应0r111-10-20B1d益-3011增-40-50;1|11j111I11111100.10.20.30.40.50.60.70.80.91/零点和极点0.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.50.510RealPart低通滤波器后的时域信号00.0050.010.