Matlab和双线性变换的滤波器设计

题目:MATLAB课程设计—基于Matlab和双线性变换的滤波器设计初始条件：MATLAB仿真软件数字信号处理与图像处理根底知识双线性变换法的原理和算法巴特沃斯高通滤波器的性能指标要求完成的主要任务:利用MATLAB仿真软件系统结合双线性变换法设计一个数字巴特沃斯高通、低通和带通IIR滤波器。说明：自己确定滤波器的不同参数得到不同的结果，最后自己分析原因。课程设计题目由指导教师提供，每人隶属一组完成任务，每组完成的内容不能雷同〔按学号分组〕时间安排：序号阶段内容所需时间1方案设计1天2软件设计2天3系统调试1天4辩论1天合计5天指导教师签名：年月日系主任〔或责任教师〕签名：年月日目录摘要数字滤波器是数字信号处理的根底，用来对信号进行过滤、检测与参数估计等处理，在通信、图像、语音、雷达等许多领域都有着十分广泛的应用。尤其在图像处理、数据压缩等方面取得了令人瞩目的进展和成就。鉴于此，数字滤波器的设计就显得尤为重要。此报告重点介绍了用双线性不变法设计IIR数字滤波器的根本流程，比拟了各种设计方法的优缺点，总结了模拟滤波器的性能特征。最后以双线性不变法设计了一个高通巴特沃斯FIR数字滤波器，介绍了设计步骤，然后在Matlab环境下进行了仿真与调试，实现了设计目标。关键字：数字滤波器巴特沃斯MATLAB双线性变法AbstractMatlabforresearchandapplicationofdigitalfilterprovidesanintuitive,efficientandconvenientweapon.Hematrixoperationsbasedonthecalculation,visualization,programdesignintegratedintoaninteractiveworkingenvironment.Inparticular,Matlabtoolboxtoenableresearchersinvariousfieldscanbeeasyandintuitiveapplicationofscientificresearchandengineering.ThisreportintroducedwithemphasisofthebasicflowofdesigningtheIIRdigitfilterbythebilinearpoliticalreform,comparedwitheachkindofdesignmethod'sgoodandbadpoints,summarizedanalogfilter'sperformancecharacteristic.FinallydesignonebythebilinearpoliticalreformtopassButterworthhighIIRdigitfilter,introducedthedesignprocedure,thenhascarriedonthesimulationandthedebuggingundertheMatlabenvironment,hasachievedtheprojectobjective.Keyword:digitalfilterButterworthMATLAB1数字滤波器1.1数字滤波器的根本概念数字滤波器是数字信号处理的重要根底，是对信号都是过滤检测与参数估计等处理过程中，它是使用最为广泛的一种线性系统。数字滤波器处理的对象是经由采样期间将模拟信号转换而得到的数字信号。数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散时间线性非时变系统。数字滤波器的输入是一组数字量。它本身既可以是用数字硬件装配而成的一台用于完成给定运算的专用数字计算机，也可以是将所需的运算编写的程序通过计算机来执行。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。随着数字技术的开展，用数字技术实现滤波器的功能愈来愈受到人们的重视，并得到了广泛的应用。数字滤波器的数学运算通常有两种实现方式。一种是频域法，即利用FFT快速运算方法对输入信号进行离散傅里叶变换，分析其频谱，然后再根据所希望的频率特性进行滤波，再利用傅里叶反变换得到时域信号。这种方法具有较好的频域选择特性和灵活性，并且由于信号频率与所希望的频率特性是简单的相乘关系，所以它比计算等价的时域卷积要快得多。另一种方法是时域法，这种方法是通过离散的抽样数据做差分数学运算来到达滤波目的的。1.2数字滤波器的分类数字滤波器按照不同的分类方法，可分为许多种，但总体来讲可以分成两大类。一类称为经典滤波器，即一般滤波器，特点是输入信号中有用的频率成分和希望滤除的频率成分各占有不同的频带，通过一个适宜的选频滤波器到达滤波的目的。但对于一般滤波器，如果信号和干扰的频带互不重叠，那么不能完成对干扰的有效滤除，这是需要采用另一类所谓的现代滤波器，例如维纳滤波器、卡尔曼滤波器、自适应滤波器等最正确滤波器，这些滤波器可按照随机信号内部的一些统计分布规律，从干扰中最正确地提取信号。从功能上，一般数字滤波器可以分为低通、高通、带通、带阻和全通等，此种分类方法是和模拟滤波器一样的。1.3数字滤波器的MATLAB实现MATLAB工具箱未滤波器的设计应用提供了丰富而简便的方法，如函数方法和图形工具方法等，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设计和实现开辟了广阔的天地。数字滤波器的一般设计过程为：按照实际需要，确定滤波器的性能要求；用一个因果的、稳定的离散线性时不变系统，去逼近这一性能指标；用有限精度的运算实现所设计的系统；通过模拟，验证所设计的系统是否符合给定性能要求。Dd2双线性变换法2.1双线性变换法知识简介脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Ｚ平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。由图2-1看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频局部超过折叠频率而混淆到低频局部去的现象，从而消除了频率混叠现象。图2-12.2双线性变换法设计数字滤波器原理将S平面j轴压缩变换到s1平面j轴上的-/T到/T一段，可以采用以下变换关系：这样变换到，变换到，可将上式写成令，，解析延拓到整个s平面和s1平面，可得再将平面通过以下变换关系映射到z平面，即从而得到s平面和z平面的单值映射关系为一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应关系，可引入待定常数c，将代入到上式，可得在MATLAB中，双线性Z变换可以通过bilinear函数实现，其调用格式为：[Bz，Az]＝bilinear(B，A，Fs)；其中B，A为模拟滤波器传递函数G(s)的分子分母多项式的系数向量，而Bz，Az为数字滤波器的传递函数H(z)的分子分母多项式的系数向量。3设计任务及方案选择3.1设计任务与要求利用MATLAB仿真软件系统结合双线性变换法设计一个数字巴特沃斯高通、低通和带通FIR滤波器。说明：自己确定滤波器的不同参数得到不同的结果，最后自己分析原因。假设一个信号：其中f1=100Hz，f2=550Hz。试将该信号分别通过低通、高通、带通滤波器，对输出信号进行分析比拟。Step1:Step2:Step3:Step4:3.2方案设计论证与选择本次滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，那么通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：方案一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻模拟域双线性变换法归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻频率变换图3-1先频率变换再离散方案二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻数字原型低通数字域归一化模拟低通原型数字高,带通或带阻数字原型低通双线性变换法频率变换图3-2先离散再频率变换以上两种方案都可以，我最后选择了第一种方案进行设计，即先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。因为高通，带通滤波器的设计方法是先将要设计的滤波器的技术特性指标通过频率转换成模拟低通滤波器的技术指标，再根据这些性能指标设计出低通滤波器的传递函数，为了设计程序的简洁，应选择了方案一。4程序设计与调试4.1设计要求假设实验设计要求得巴特沃斯高通滤波器的性能指标为：通带边界频率fp：fp=400Hz；通带波纹Rp：Rp=1dB；阻带边界频率fs：fs=300Hz；阻带衰减Rs：Rs=30dB；采样频率Fs：Fs=2000Hz；要求巴特沃斯低通滤波器通带边界频率与阻带边界频率与高通的正好相反。模拟带通滤波器中心频率为f0=200Hz，带宽为Bw=400Hz，4.2程序设计与说明4.2.1模拟低通滤波首先应用双线性变换作为数字低通滤波器到模拟高通滤波器的频率预畸变处理：利用上式确定的技术指标和的参数Rp、Rs，就可以进行模拟原型低通滤波器设计了。该局部的程序设计与注释如下：Fs=2000;Ts=1/Fs;%采样频率和采样周期wp=400*2*pi;ws=300*2*pi;%低通、高通数字滤波器截止频率Rp=1;Rs=30;%通带波纹和阻带衰减wpl=Fs*2*tan(ws*Ts/2);%按频率转换公式进行转换wsl=Fs*2*tan(wp*Ts/2);[N,Wn]=buttord(wpl,wsl,Rp,Rs,'s');%计算模拟滤波器的最小阶数[z,p,k]=buttap(N);%设计模拟原型低通滤波器[b,a]=zp2tf(z,p,k);%变零极点增益形式为系统传递函数形式[H,w]=freqs(b,a);%求模拟低通滤波器的频率响应figure;plot(w,abs(H));grid;%绘制低通频率特性title('模拟低通原型');xlabel('w');ylabel('幅度');4.2.2模拟到数字的转换〔1〕从低通向低通的转换[b,a]=lp2lp(Bap,Aap,Wn);功能是把模拟滤波器原型转换成截至频率为Wn的低通滤波器。再由得到的模拟低通滤波器通过函数[Bz，Az]＝bilinear(B，A，Fs)得到数字低通滤波器。该局部的程序设计如下：[e1,f1]=lp2lp(b,a,300*2*pi);%原型低通转换为低通模拟滤波器[h1,w1]=freqs(e1,f1);%求模拟低通滤波器的频率响应figure;plot(w1/(2*pi),abs(h1));%绘制高通频率特性grid;title('模拟低通滤波器');xlabel('频率Hz');ylabel('幅度');[r1,q1]=bilinear(e1,f1,1000);%模拟低通滤波器转换为数字滤波器[s1,w2]=freqz(r1,q1);%求数字低通滤波器的频率响应figure;plot(w2/pi,20*log10(abs(s1)));grid;%绘制数字低通滤波器的频率响应title('数字低通滤波器');xlabel('w/pi');ylabel('幅度(dB)');〔2〕从低通向高通的转换[b,a]=lp2hp(Bap,Aap,Wn);功能是把模拟滤波器原型转换成截至频率为Wn的模拟高通滤波器，再进一步转换成数字高通滤波器。程序设计如下：[e2,f2]=lp2hp(b,a,300*2*pi);%原型低通模拟滤波器转换为高通模拟滤波器[h2,w3]=freqs(e2,f2);%求模拟高通滤波器的频率响应figure;plot(w3/(2*pi),abs(h2));%绘制高通频率特性grid;title('模拟高通滤波器');xlabel('频率Hz');ylabel('幅度');[r2,q2]=bilinear(e2,f2,1000);%模拟高通滤波器转换为数字滤波器[s2,w4]=freqz(r2,q2);%求数字高通滤波器的频率响应figure;plot(w4/pi,20*log10(abs(s2)));grid;%绘制数字高通滤波器的频率响应title('数字高通滤波器');xlabel('w/pi');ylabel('幅度(dB)');〔3〕从低通向带通的转换[b,a]=lp2bp(Bap,Aap,W0，Bw);功能是把模拟滤波器原型转换成中心频率为W0、带宽为Bw的模拟带通滤波器，在进一步转换成数字带通滤波器。程序设计如下：W0=200*2*pi;Bw=200;[e3,f3]=lp2bp(b,a,W0,Bw);%原型低通转换为带通模拟滤波器[h3,w5]=freqs(e3,f3);%求模拟带通滤波器的频率响应figure;plot(w5/(2*pi),abs(h3));%绘制带通频率特性grid;title('模拟带通滤波器');xlabel('频率Hz');ylabel('幅度');[r3,q3]=bilinear(e3,f3,1000);%模拟滤波器转换为数字滤波器[s3,w6]=freqz(r3,q3);%求数字带通滤波器的频率响应figure;plot(w6/pi,20*log10(abs(s3)));grid;%绘制数字带通滤波器的频率响应title('数字带通滤波器');xlabel('w/pi');ylabel('幅度(dB)');4.1.3对输入信号进行滤波读取输入信号程序：n=0:180;dt=1/Fs;t=n*dt;%时间序列f1=100;f2=550;%输入信号频率x=sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*1*t);%输入信号figure;plot(t,x);grid;%绘制输入信号title('输入信号');xlabel('时间t');ylabel('幅度');将输入信号分别通过得到的低通、高通、带通滤波器，调用函数y=filter(bz,az,x)对输入信号进行滤波。y=filter(r1,q1,x);%对输入信号进行低通滤波figure;plot(t,y);grid;%绘制输出信号title('低通滤波输出信号');xlabel('时间t');ylabel('幅度')再进行高通、带通滤波只需将调用函数filter中参量改为r2、q2和r3、q3即可。4.3仿真结果与分析将设计好的程序在matlab软件下进行仿真，得到仿真波形。仿真结果记录图4-1模拟低通原型幅频响应〔a〕低通模拟滤波器的频率特性曲线〔b〕巴特沃斯低通数字滤波器的频率特性曲线图4-2〔a〕高通模拟滤波器的频率特性曲线〔b〕巴特沃斯高通数字滤波器的频率特性曲线图4-3〔a〕带通模拟滤波器的频率特性曲线〔b〕巴特沃斯带通数字滤波器的频率特性曲线图4-4图4-5输入信号波形图4-6经过巴特沃斯低通滤波器滤波后的输出波形图4-7经过巴特沃斯高通滤波器滤波后的输出波形图4-8经过巴特沃斯高通滤波器滤波后的输出波形4.3.2仿真结果分析由仿真结果可知，设计出的数字低通，高通，带通的幅频响应均符合各自的特性，从模拟滤波器特性曲线可以看出，所设计的高通滤波器在大于400Hz为通带，其衰减均小于1dB；小于300Hz为阻带，其衰减大于30dB，而设计的低通滤波器情况正好相反，符合本次设计的要求。当滤波器输入100Hz和550Hz两种信号后，低通滤波器可以滤除频率为550Hz的高频信号，而低频信号100Hz的信号可以通过滤波器，验证了该低通滤波器的功能。而高通滤波器可以滤除频率为100Hz的信号，而高频信号550Hz的信号可以通过滤波器。4.3.3程序调试中出现的问题带通滤波器的仿真结果虽然符合其带通特性，但要求得到的性能指标存在较大误差，由于时间紧迫，将在课后进一步寻找原因及解决方案。一个线性是不变系统的典型表达式有状态空间型、传递函数型、零极增益型、极点留数型共四种，在开始的设计过程中我发现在把零极点增益型向传递函数型的转换过程中多用了一步，即先把零极点增益型转换为状态空间型，进行相应的频率变换和双线性变换后，再由状态空间型转换为传递函数型，进而绘制出其频率响应图，但是经仿真分析后发现，性能并不理想，所以并没有采用。5心得体会本次课程设计是基于数字滤波器的MATLAB软件应用，在本学期中，我在数字信号处理课程中学习了关于本次课程设计的一些内容，对于MATLAB软件在数字滤波器的设计中的一些函数也有所认识，通过着一周的时间更是让自己在这个方面的学习更进一步，使自己在这方面的学习深度有所加深。本次用双线性变换法对巴特沃斯低通，高通及带通滤波器的设计，要求我能够用所学理论知识对巴特沃斯高通滤波器进行设计，并能够对设计的滤波器的频率特性进行分析。在设计过程中，首先，我通过查找资料等途径，对关于所要设计的题目进行更加全面深入的学习，认真分析其中的原理，弄明白设计的思路；其次，根据自己选择的性能指标要求，结合自己通过学习掌握的知识和资料进行程序设计，根据自己提出的方案利用软件实现；最后，运行程序，分析所得结果是否符合自己的设计要求，对结果进行分析，实现由理论到实践，再由实践结果验证理论的过程。总的来说，这次课程设计让我对MATLAB有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的根本方法，提高了理论用于实践的能力，掌握了更多专业相关的使用知识与技能。同时，也暴露了我很多的缺乏，在以后的学习中，将进一步发扬有点，克服缺点。参考文献[1]袁润章，自蔓延高温合成技术研究进展.武汉:武汉工业大学出版社,1994[2]SchachtE.Industrialpolysaccharides.Amsterdam:ElsevierScience,1987[3]郭景坤.陶瓷材料的强化与增韧新途径的探索.无机材料学报,1998,13910;23-26[4]SpriggsGE.Ahistoryoffinegrainedhardmetal.IntJofRefractoryMetalandHardMaterial,1995,13:241-255[5]王为民.二硼化钛陶瓷的自蔓延高温合成与加工:[博士论文].武汉工业大学,1988[6]吴中庆,张基田.SiC颗粒与液体铝的界石层.广州:92秋季中国材料科学研讨会,C-MRS,1992.603-609.附录本次课设总程序如下：Fs=2000;Ts=1/Fs;%采样频率和采样周期wp=400*2*pi;ws=300*2*pi;%低通、高通数字滤波器截止频率Rp=1;Rs=30;%通带波纹和阻带衰减wpl=Fs*2*tan(ws*Ts/2);%按频率转换公式进行转换wsl=Fs*2*tan(wp*Ts/2);[N,Wn]=buttord(wpl,wsl,Rp,Rs,'s');%计算模拟滤波器的最小阶数[z,p,k]=buttap(N);%设计模拟原型低通滤波器[b,a]=zp2tf(z,p,k);%变零极点增益形式为系统传递函数形式[H,w]=freqs(b,a);%求模拟低通滤波器的频率响应figure;plot(w,abs(H));grid;%绘制低通频率特性title('模拟低通原型');xlabel('w');ylabel('幅度');[e1,f1]=lp2lp(b,a,300*2*pi);%原型低通转换为低通模拟滤波器[h1,w1]=freqs(e1,f1);%求模拟低通滤波器的频率响应figure;plot(w1/(2*pi),abs(h1));%绘制高通频率特性grid;title('模拟低通滤波器');xlabel('频率Hz');ylabel('幅度');[r1,q1]=bilinear(e1,f1,1000);%模拟低通滤波器转换为数字滤波器[s1,w2]=freqz(r1,q1);%求数字低通滤波器的频率响应figure;plot(w2/pi,20*log10(abs(s1)));grid;%绘制数字低通滤波器的频率响应title('数字低通滤波器');xlabel('w/pi');ylabel('幅度(dB)');[e2,f2]=lp2hp(b,a,300*2*pi);%原型低通模拟滤波器转换为高通模拟滤波器[h2,w3]=freqs(e2,f2);%求模拟高通滤波器的频率响应figure;plot(w3/(2*pi),abs(h2));%绘制高通频率特性grid;title('模拟高通滤波器');xlabel('频率Hz');ylabel('幅度');[r2,q2]=bilinear(e2,f2,1000);%模拟高通滤波器转换为数字滤波器[s2,w4]=freqz(r2,q2);%求数字高通滤波器的频率响应figure;plot(w4/pi,20*log10(abs(s2)));grid;%绘制数字高通滤波器的频率响应title('数字高通滤波器');xlabel('w/pi');ylabel('幅度(dB)');W0=200*2*pi;Bw=200;[e3,f3]=lp2bp(b,a,W0,Bw);%原型低通转换为带通模拟滤波器[h3,w5]=freqs(e3,f3);%求模拟带通滤波器的频率响应figure;plot(w5/(2*pi),abs(h3));%绘制带通频率特性grid;title('模拟带通滤波器');xlabel('频率Hz');ylabel('幅度');[r3,q3]=bilinear(e3,f3,1000);%模拟滤波器转换为数字滤波器[s3,w6]=freqz(r3,q3);%求数字带通滤波器的频率响应figure;plot(w6/pi,20*log10(abs(s3)));grid;%绘制数字带通滤波器的频率响应n=0:180;dt=1/Fs;t=n*dt;%时间序列f1=100;f2=550;%输入信号频率x=sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*1*t);%输入信号figure;plot(t,x);grid;%绘制输入信号title('输入信号');xlabel('时间t');ylabel('幅度');y=filter(r1,q1,x);%对输入信号进行低通滤波figure;plot(t,y);grid;%绘制输出信号title('低通滤波输出信号');xlabel('时间t');ylabel('幅度')y=filter(r2,q2,x);%对输入信号进行高通滤波figure;plot(t,y);grid;