课 程 设 计课程设计名称: 数字信号处理 数字信号处理 专业课程设计任务书题 目 用双线性变换法设计原型低通为巴特沃兹型的数字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.“课题性质”一栏:A.工程设计;B.工程技术研究;C.软件工程(如CAI课题等);D.文献型综述;E.其它。2.“课题来源”一栏:A.自然科学基金与部、省、市级以上科研课题;B.企、事业单位委托课题;C.校、院(系、部)级基金课题;D.自拟课题。1 需求分析本课程设计要求用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器,给定的设计参数为通带边界频率为400Hz,500Hz,阻带边界频率分别为350Hz,550Hz,通带最大衰减1dB,阻带最小衰减40dB,抽样频率为2000Hz。要求采用双线性变换法,使得s平面与z平面是单值的一一对应关系,不存在频谱混淆的问题。由于数字频域和模拟频域的频率不是线性关系,这种非线性关系使得通带截止频率、过渡带的边缘频率的相对位置都发生了非线性畸变。设计出巴特沃兹型的带通滤波器之后,让信号 经过该滤波器,其中 450Hz, 600Hz,分析滤波器的输出 是什么,用Matlab验证结论并给出 的图形。2 概要设计1>用双线性变换法设计IIR数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点,可以采用非线性频率压缩方法,将整个频率轴上的频率范围压缩到-π/T~π/T之间,再用z=esT转换到Z平面上。也就是说,第一步先将整个S平面压缩映射到S1平面的-π/T~π/T一条横带里;第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系,消除了多值变换性,也就消除了频谱混叠现象,映射关系如图1-3所示。图1-3双线性变换的映射关系为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上,可以通过以下的正切变换实现(1-5)式中,T仍是采样间隔。当Ω1由-π/T经过0变化到π/T时,Ω由-∞经过0变化到+∞,也即映射了整个jΩ轴。将式(1-5)写成将此关系解析延拓到整个S平面和S1平面,令jΩ=s,jΩ1=s1,则得再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为:(1-6)(1-7)式(1-6)与式(1-7)是S平面与Z平面之间的单值映射关系,这种变换都是两个线性函数之比,因此称为双线性变换式(1-5)与式(1-6)的双线性变换符合映射变换应满足的两点要求。首先,把z=ejω,可得(1-8)即S平面的虚轴映射到Z平面的单位圆。其次,将s=σ+jΩ代入式(1-8),得因此由此看出,当σ<0时,|z|<1;当σ>0时,|z|>1。也就是说,S平面的左半平面映射到Z平面的单位圆内,S平面的右半平面映射到Z平面的单位圆外,S平面的虚轴映射到Z平面的单位圆上。因此,稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。2>双线性变换法优缺点双线性变换法与脉冲响应不变法相比,其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单值地对应于Z平面单位圆一周,即频率轴是单值变换关系。这个关系如式(1-8)所示,重写如下:上式表明,S平面上Ω与Z平面的ω成非线性的正切关系,如图7-7所示。由图7-7看出,在零频率附近,模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系;但当Ω进一步增加时,ω增长得越来越慢,最后当Ω→∞时,ω终止在折叠频率ω=π处,因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象,从而消除了频率混叠现象。图1-4双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的,如式(1-8)及图1-4所示。由于这种频率之间的非线性变换关系,就产生了新的问题。首先,一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器,不再保持原有的线性相位了;其次,这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的,即某一频率段的幅频响应近似等于某一常数(这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性),不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变,如图1-5所示。图1-5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器,双线性变换后,仍得到幅频特性为分段常数的滤波器,但是各个分段边缘的临界频率点产生了畸变,这种频率的畸变,可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变,然后经变换后正好映射到所需要的数字频率上。根据以上IIR数字滤波器设计方法,下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器,通带截止频率fp1=500HZ,fp2=400HZ;阻带截止频率fs1=350HZ,fs2=550HZ;通带最大衰减αp=1dB;阻带最小衰减αs=40dB,抽样频率为2000HZ。I.设计步骤:(1)根据任务,确定性能指标:通带截止频率fp1=500HZ,fp2=400HZ;阻带截止频率fs1=350HZ,fs2=550HZ;通带最大衰减αp=1dB;阻带最小衰减αs=40dB;抽样频率为2000HZ。 (2)用Ω=(2/T)tan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变,得到带通模拟滤波器H(s)的边界频率主要是通带截止频率fp1,fp2;阻带截止频率fs1,fs2的转换。为了计算简便,对双线性变换法一般T=2s通带截止频率fc1=(2/T)*tan(2*pi*fp1/2)fc2=(2/T)*tan(2*pi*fp2/2)阻带截止频率fr1=(2/T)*tan(2*pi*fs1/2)fr2=(2/T)*tan(2*pi*fs2/2)阻带最小衰减αs=1dB和通带最大衰减αp=40dB;(3)运用低通到带通频率变换公式λ=(((f^2)-(f0^2))/ f)将模拟带通滤波器指标转换为模拟低通滤波器指标。B=fc2-fc1normwr1=(((fc1*fc2)-(fc1^2))/fc1) normwr2=(((fc1*fc2)-(fc2^2))/fc2)normwc1=(((fc1*fc2)-(fr1^2))/fr1)normwc2=(((fc1*fc2)-(fr1^2))/fr1)模拟低通滤波器指标:normfc,normfr,αp,αs(4)设计模拟低通原型滤波器。借助巴特沃斯(Butterworth)滤波器,用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数H(s);(5)调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。(6)利用双线性变换法将模拟带通滤波器H(s)转换成数字带通滤波器H(z).II.程序流程框图:开始↓读入数字滤波器技术指标↓将指标转换成归一化模拟低通滤波器的指标↓设计归一化的模拟低通滤波器阶数N和1db截止频率↓模拟域频率变换,将G(P)变换成模拟带通滤波器H(s)↓用双线性变换法将H(s)转换成数字带通滤波器H(z)↓输入信号后显示相关结果↓结束3 运行环境PC机 MATLAB4 开发工具和编程语言MATLAB语言5 详细设计clc;clear all;%设计滤波器%所需设计的带通滤波器的参数设置Fres=2000;Ts=1/Fres;Omgap1=500Omgap2=400Omgas1=350Omgas2=550%进行双线性变换,使得s平面与z平面是单值的一一对应关系Omgap=(Omgap1-Omgap2)*TsOmgap3=tan(Omgap1*Ts*2*pi/2)Omgap4=tan(Omgap2*Ts*2*pi/2)Omgas3=tan(Omgas1*Ts*2*pi/2)Omgas4=tan(Omgas2*Ts*2*pi/2)%对参数归一化ap1=Omgap3/Omgapap2=Omgap4/Omgapas1=Omgas3/Omgapas2=Omgas4/Omgap%将模拟带通滤波器指标转化成模拟低通滤波器指标I1=(ap1*ap2-as1*as1)/as1I2=(ap1*ap2-as2*as2)/as2I3=(ap1*ap2-ap1*ap1)/ap1I4=(ap1*ap2-ap2*ap2)/ap2I5=min(I1,-I2)alfpp=1;alfps=40;%设计巴特沃兹型低通滤波器HL(s)[N,Wn]=buttord(I4,I5,alfpp,alfps,'s')[Num,Den]=butter(N,Wn,'s')%将设计的模拟低通滤波器传递函数HL(s)转化成模拟带通滤波器H(s)[Num2,Den2]=lp2bp(Num,Den,sqrt(Omgap3*Omgap4),Omgap);%将设计的模拟带通滤波器H(s)转化成数字带通滤波器H(Z)[Num3,Den3]=bilinear(Num2,Den2,0.5);%画出H(Z)幅频特性曲线w=0:pi/255:pi;h=freqz(Num3,Den3,w);H=20*log10(abs(h));plot(w,H);%设计信号函数f1=450;f2=600;t=0:1/2000:40/f2;f=sin(2*pi*f1*t)+sin(2*pi*f2*t);plot(f);%信号通过带通滤波器g=invfreqz(h,w,40,50)g1=conv(g,f)plot(g1)6 调试分析在设计滤波器的时候,计算的是幅频特性。而相频特性没有进行频谱分析。而信号和信号通过滤波器的图形分析都为时与分析,没有进行频域分析,是不知道该用哪个函数对其进行傅里叶变换。试过freqz,fft等函数,但是都没有出来结果。说明调用格式不正确。本次试验用的是巴特沃兹型低通滤波器设计的方法,还可以利用切比雪夫型低通滤波器设计,或者椭圆型等。在本次试验中,我们只用了两种频率简单的叠加作为输入信号,但实际信号是多种频率的组合,可以再增加一些频率。现实中还应该有噪声的影响,本实验中没有考虑。可以再加上噪声信号。7 测试结果Fres = 2000Ts = 5.0000e-004Omgap1 = 500Omgap2 = 400Omgas1 = 350Omgas2 = 550Omgap = 0.0500Omgap3 = 1.0000Omgap4 = 0.7265Omgas3 =0.6128Omgas4 = 1.1708ap1 = 20.0000ap2 = 14.5309as1 = 12.2560as2 = 23.4170I1 = 11.4562I2 = -11.0065I3 = -5.4691I4 = 5.4691I5 = 11.0065alfpp = 1alfps = 40N = 8Wn = 6.1894Num = 1.0e+006 * Columns 1 through 8 0 0 0 0 0 0 0 0 Column 9 2.1538Den = 1.0e+006 * Columns 1 through 8 0.0000 0.0000 0.0005 0.0052 0.0377 0.1984 0.7386 1.7837 Column 9 2.1538Num2 = 1.0e-004 * Columns 1 through 8 0.8413 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Column 9 -0.0000Den2 = Columns 1 through 8 1.0000 1.5863 7.0705 8.7151 20.5005 19.9985 32.1353 24.8474 Columns 9 through 16 29.9185 18.0527 16.9631 7.6698 5.7123 1.7643 1.0400 0.1695 Column 17 0.0776Num3 = 1.0e-004 * Columns 1 through 8 0.0043 0.0000 -0.0341 0.0000 0.1194 0.0000 -0.2389 0.0000 Columns 9 through 16 0.2986 0.0000 -0.2389 0.0000 0.1194 0.0000 -0.0341 0.0000 Column 17 0.0043Den3 = Columns 1 through 8 1.0000 -2.2463 8.3946 -13.4519 27.6744 -33.6694 48.3386 -45.7295 Columns 9 through 16 49.5687 -36.4140 30.6576 -16.9904 11.1207 -4.2944 2.1332 -0.4524 Column 17 0.1603h = Columns 1 through 5 -0.0000 -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i …… Columns 251 through 256 -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000iH = Columns 1 through 8 -279.2737 -279.2732 -279.2723 -279.2813 -279.3855 -280.0020 -283.0382 -292.4507…… Columns 249 through 256 -281.0032 -280.3352 -280.1410 -280.0979 -280.0943 -280.0973 -280.0996 -280.1004f1 = 450f2 = 600t = Columns 1 through 8 0 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 …… Columns 129 through 134 0.0640 0.0645 0.0650 0.0655 0.0660 0.0665f = Columns 1 through 8 0 1.9387 -0.2788 -1.4788 0.3633 0.7071 -0.1420 0.1338…… Columns 129 through 134 -0.3633 -0.7946 1.0000 1.1075 -1.5388 -1.0418g = Columns 1 through 8 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0001 -0.0002…… Columns 33 through 41 -0.0002 0.0000 0.0001 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000g1 = Columns 1 through 8 0 0.0000 0.0000 -0.0000 -0.0000 0.0001 0.0001 -0.0003…… Columns 169 through 174 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000综合上述分析,数据基本正确。上图为带通滤波器幅频响应 、上图为传输函数零极点分析 上图为输入信号时域上图为输出波形的时域分析参考文献 [1] 吴大正 《信号与线性系统分析》第四版 高等教育出版社[2] 郑君里 《信号与系统》第二版 高等教育出版社[3] Sanjit K. Mitra 《数字信号处理—基于计算机的方法》第三版 清华大学出版社[4] 余成波 《数字信号处理及MATLAB实现》清华大学出版社[5] 周利清 《数字信号处理基础》 北京邮电大学出版社[6] 美国莱昂 《数字信号处理》英文第二版 机械工业出版社心得体会................0.0