数字信号处理课程设计(1)MATLAB实现

PAGE数字信号处理课程设计报告选题名称:IIR滤波器分析与设计 系（院）:计算 机 工 程 学院 专业: 通信工程 班级: 姓名:学号: 指导教师: 学年学期: 2010~2011学年第1学期2011年01

摘要：随着社会的进步，数字信号处理技术也在飞跃的发展，作为通信工程的学生要对数字信号有更深一层的理解，本次要做的课程设计就是用MATLAB做滤波器的设计，这次课程设计是检验本学期学习的数字信号处理这门课，加深了对滤波器设计原理以及实现方法的理解。常用的设计低通滤波器的方法有脉冲响应不变法以及双线性法，这两种方法都有各自的优缺点，可以通过编写MATLAB程序观察他们的幅度特性，相位特性，以及更高一层次的就是用自己设计的滤波器来实现数字滤波。现在也有相当成熟的滤波器设计模型，如巴特沃斯，以及切比雪夫滤波器。通过设计界面，编写函数，最终还可以达到可以将自己录制的声音加载到设计的滤波器中，听滤波前后的声音的变化，这种更加直观的感受滤波器的效果的方法，而此次设计的滤波器在现实当中以用也是相当广范的，这也是学习数字信号处理这门课程，以及掌握滤波器设计方法的现实所在。关键词：数字信号；脉冲响应不变法；双线性法；数字滤波

目 录TOC\o"1-2"\u1课题综述 11.1课程设计的意义 11.2课程设计的目的 11.3课程设计面对的问题以及关键技术 12系统分析 22.1整体分析 22.2脉冲响应不变法 22.3双线性变换法的基本原理 32.4巴特沃斯滤波器 52.5切比雪夫滤波器 52.6总体设计图 53代码编写 63.1显示幅度谱函数IIR_DispMag() 63.2创建滤波器IIR_Createfilter() 73.3显示相位谱函数IIR_DispAngle() 83.4显示零极点图函数IIR_DispZplane() 83.5选择波形文件IIR_GetWAVFile() 83.6播放原始文件IIR_Playorg 83.7播放时滤波后文件IIR_Playmod 94代码的运行与调试 95软件运行与测试 9总结 12致谢 13参考文献 14设计任务书课题名称IIR滤波器分析与设计设计目的理解并掌握无限脉冲响应数字滤波器（IIR）的机理，分析IIR滤波器的结构特性，观察IIR滤波器的频域特性；学习IIR数字滤波器的三种典型的设计方法,并观察、比较设计结果；熟悉IIR数字滤波器的计算机编程；加深对课堂知识的理解与灵活应用，学习开发资料的收集与整理，学会撰写课程设计报告。实验环境微型电子计算机（PC）；MATLAB6.5或VisualC++6.0等开发环境。任务要求利用课余时间查阅课题相关资料，深入理解课题含义及设计要求；在第18周完成预设计，并请指导教师审查，通过后方可进行下一步工作；按指导书要求设计软件；要求形成稳定的程序软件，可以运行，方可申请答辩；结束后，及时提交课程设计报告（含纸质稿、电子稿）。工作进度计划序号起止日期工作内容12011.01.03~2011.01.03在预设计的基础上，进一步查阅资料，完善设计方案，形成书面材料。22011.01.04~2011.01.05设计总体方案，绘制流程框图，编写代码，上机调试。32011.01.06~2011.01.07测试程序，优化代码，增强功能，撰写设计报告。42011.01.08~2011.01.08提交软件代码、设计报告，参加答辩，根据教师反馈意见，修改、完善设计报告。指导教师（签章）：年月日数字信号处理课程设计报告PAGE141课题综述数字滤波器是数字信号处理理论的一部分。数字信号处理主要是研究用数字或符号的序列来表示信号波形，并用数字的方式去处理这些序列，把它们改变成在某种意义上更为有希望的形式，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。具体来说，凡是用数字方式对信号进行滤波、变换、调制、解调、均衡、增强、压缩、固定、识别、产生等加工处理，都可纳入数字信号处理领域。1.1课程设计的意义《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2课程设计的目的1.理解并掌握无限脉冲响应数字滤波器（IIR）的机理，分析IIR滤波器的结构特性，观察IIR滤波器的频域特性；2.学习IIR数字滤波器的三种典型的设计方法,并观察、比较设计结果；3.熟悉IIR数字滤波器的计算机编程；4.加深对课堂知识的理解与灵活应用，学习开发资料的收集与整理，学会撰写课程设计报告。1.3课程设计面对的问题以及关键技术本次课程设计是用MATLAB软件来实现的，所以第一个需要解决的问题就是如何把自己的思想转化为MATLAB语言，但是在这次的课程设计中要感谢我们的老师，因为他已经将大部分的代码提供给我们了。本次设计的关键技术就是脉冲响应不变法以及双线性法的函数调用以及求极点，零点的Z变换函数调用。最后就是如何实现语音信号的数字滤波。2系统分析2.1整体分析利用模拟滤波器设计数字滤波器，就是从已知的模拟滤波器系统函数Ha(s)设计数字滤波器系统函数H(z)，这归根到底是一个由S平面到Z平面的变换，这种映射变换应遵循两个基本原则： 1）H（z）的频响要能与Ha(s)的频响保持一致，即S平面的虚轴应映射到Z平面的单位圆上。 2）Ha(s)的因果稳定性映射成H（z）后保持不变，即S平面的左半平面Re{S}＜0应映射到Z平面的单位圆以内|Z|<1。所以滤波器的设计就是在给定了滤波器的技术指标后，确定滤波器的阶数N和系数和系数{ai，bi}。在满足技术指标的条件下，滤波器的阶数应尽可能低，因为滤波器的阶数越低，实现滤波器的成本就越低。2.2脉冲响应不变法脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)正好等于模拟滤波器的冲激响应ha(t)的采样值，即h(n)=ha(nT),T为采样周期。如以Ha(s)及H（z）分别表示ha(t)的拉氏变换及h(n)的Z变换，即Ha(s)=L[ha(t)],H(z)=Z[h(n)]计算H(Z):脉冲响应不变法特别适用于用部分分式表达系统函数，模拟滤波器的系统函数若只有单阶极点，且分母的阶数高于分子阶数N＞M，则其拉氏反变换为：（2.2.1）单位阶跃对ha(t)采样得到数字滤波器的单位脉冲响应序列（2.2.2）再对h(n)取Z变换，得到数字滤波器的传递函数：（2.2.3）第二个求和为等比级数之和，要收敛的话，必有所以有（2.2.4）根据理想采样序列拉氏变换与模拟信号拉氏变换的关系（2.2.5）以上表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是拉氏变换到Z变换的标准变换关系，即首先对Ha(s)作周期延拓，然后再经过z=est的映射关系映射到Z平面上。稳定性：如果模拟滤波器是稳定的，则所有极点si都在s左半平面，即Re[si]＜0,那么变换后H(Z)的极点est也都在单位圆以内，即：因此数字滤波器保持稳定。2.3双线性变换法的基本原理脉冲响应不变法使得数字滤波器在时域上能够较好的模仿模拟滤波器，但是由于从S平面到Z平面的映射具有多值性，使得设计出来的数字滤波器不可避免的出现频谱混迭现象。为了克服脉冲响应不变法可能产生的频谱混跌效应的缺点，我们使用一种新的变换双线性变换。双线性变换法可认为是基于对微分方程的积分，利用对积分的数值逼近的道德。仿真滤波器的传递函数为(2.2.5)将展开为部份分式的形式，并假设无重复几点，则(2.2.6)那么，对于上述函数所表达的数字信号处理系统来讲，其仿真输入和模拟输出有如下关系(2.2.7)利用差分方程来代替导数，即(2.2.8)同时令这样，便可将上面的微分方程写为对应的差分方程形式(2.2.9)两边分别取变换，可得(2.2.10)这样，通过上述过程，就可得到双线性变换中的基本关系，如下所示所谓的双线性变换，仅是指变换公式中与的关系无论是分子部份还是分母部份都是线性的。MATLAB中设计IIR数字滤波器的具体步骤如下：1.把给出的数字滤波器的性能指标转换为模拟低通滤波器的性能指标；2.根据转换后的性能指标，通过滤波器结束选择函数，来确定滤波器的最小阶数n和固有频率wn；3.由最小阶数n得到低通滤波器原型；4.由固有频率wn把模拟低通滤波器转换为模拟低通、高通、带通或带阻滤波器；运用双线性变换法把模拟滤波器转换成数字滤波器。2.4巴特沃斯滤波器巴特沃斯滤波器是最常用的滤波器，而巴特沃斯滤波器又可以采用脉冲响应不变法和双线性法来设计，本文主要讲述的是脉冲相应不变法的的设计。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满足：h(n)=ha(nT)，式中,T是采样周期。这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内，才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真。但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。2.5切比雪夫滤波器切比雪夫滤波器的幅度特性具有等波纹特性，它有两种形式，这里仅描述切比雪夫I滤波器，其振幅特性在通带内是等波纹的、在阻带内是单调递减的。切比雪夫滤波器特点是在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。2.6总体设计图设计的步骤流程：求出模拟低通滤波器指标求出模拟低通滤波器指标设计模拟低通滤波器将模拟滤波器转换为数字滤波器利用MATLAB编程实现滤波器得到结果开始图2.1设计流程图MATLAB数字滤波器设计整体图：图2.2整体界面3代码编写3.1显示幅度谱函数IIR_DispMag()functionIIR_DispMag()%显示幅度谱函数IIR_DispMag[numdenCYPLYYYFlag]=IIR_Createfilter;%参数调用以设计好的滤波器参数if(~YYYFlag)return;endfigure(1);[hf]=freqz(num,den,512,CYPL);%转换成数字滤波器plot(f,abs(h));%回执幅频曲线gridon;zoomon;3.2创建滤波器IIR_Createfilter()3.2.1Wp=2*pi*f2/f1;Ws=2*pi*f3/f1;OmigaP=Wp*f1;脉冲相应不变法的数字指标与模拟指标频率转换OmigaS=Ws*f1;[nW]=buttord(OmigaP,OmigaS,Rp,Rs,'s');%调用巴特沃斯滤波器函数[BA]=butter(n,W,'s');%默认为低通滤波器[num,den]=impinvar(B,A,f1);3.2.2Wp=2*pi*f2/f1;Ws=2*pi*f3/f1;OmigaP=Wp*f1;OmigaS=Ws*f1;[nW]=buttord(OmigaP,OmigaS,Rp,Rs,'s');[BA]=butter(n,W,'high','s');%高通滤波器设计[num,den]=impinvar(B,A,f1);3.2.Wp=2*pi*f2/f1;Ws=2*pi*f3/f1;wp1=Wp*f1;%频率转换ws1=Ws*f1;[nW]=cheb1ord(wp1,ws1,Rp,Rs,'s');%切比雪夫I型函数调用[BA]=cheby1(n,1,W,'low','s');[num,den]=bilinear(B,A,f1);3.2.4切比雪夫I型Wp=2*pi*f2/f1;Ws=2*pi*f3/f1;wp1=Wp*f1;%频率转换ws1=Ws*f1;[nW]=cheb1ord(wp1,ws1,Rp,Rs,'s');%切比雪夫I型函数调用[BA]=cheby1(n,1,W,'high','s');[num,den]=bilinear(B,A,f1);3.3显示相位谱函数IIR_DispAngle()functionIIR_DispAngle()[numdenCYPLYYYFlag]=IIR_Createfilter;if(~YYYFlag)return;endfigure(2);[hf]=freqz(num,den,512,CYPL);plot(f,(angle(h)));gridon;zoomon;3.4显示零极点图函数IIR_DispZplane()functionIIR_DispZplane()[numdenCYPLYYYFlag]=IIR_Createfilter;if(~YYYFlag)return;endfigure(3);zplane(num,den);%Z域变换axis([-5,5,-2,2]);3.5选择波形文件IIR_GetWAVFile()functionsig=IIR_GetWAVFile()[name,path]=uigetfile('*.wav','请选择一个.wav波形文件');%选取文件的函数调用file_path=[path,name];if(~(isstr(name))|~(isstr(path)))sig=[];return;end[datafreqbits]=wavread(file_path);time=length(data)/freq;%info=['源信号-'name';'str2num(time)'秒;fs='str2num(freq)'Hz'];s=[];s.data=data;s.freq=freq;s.bits=bits;%=info;s.time=time;data_hndl=findobj(gcbf,'Tag','Button_WAV');set(data_hndl,'Userdata',s);3.6播放原始文件IIR_PlayorgfunctionIIR_Playorg()data_hndl=findobj(gcbf,'Tag','Button_WAV');sig_info=get(data_hndl,'Userdata');data=sig_info.data;fs=sig_info.freq;sound(data,fs);3.7播放时滤波后文件IIR_PlaymodfunctionIIR_Playmoddata_hndl=findobj(gcbf,'Tag','Button_WAV');sig_info=get(data_hndl,'Userdata');data=sig_info.data;fs=sig_info.freq;[BAsuc1suc2]=IIR_Createfilter;if(~suc2)return;endout=filter(B,A,data);%音频文件通过选择的滤波器sound(out,fs);4代码的运行与调试本次试验的大部分代码是书上给出的，但是有点小问题，就是在自编函数哪有错误，原句是：if(strcmp(fil_name,’错误代码’))，修改过后的if(strcmp(fil_name,’巴特沃斯’))，这样程序就可以运行了。再有就是要录一段自己的音频文件，并且格式是WMA的，最好是有高音和低音区别的，因为这样我们就可以更清楚的辨别滤波器的滤音效果。5软件运行与测试在巴特沃斯低通滤波器环境情况下，设置通带截止频率与阻带截止频率分别为100，500；得到的幅度、相位，零极点图如图示：幅度响应(a)相位响应(b)零极点(c)图5.1低通的幅度、相位，零极点图在巴特沃斯高通滤波器环境情况下，设置通带截止频率与阻带截止频率分别为500，100；得到的幅度、相位，零极点图如图示：幅度响应(a)相位响应(b)零极点(c)图5.2高通的幅度、相位，零极点图在切比沃雪夫I高通滤波器环境下，设置通带截止频率与阻带截止频率分别为200，100；得到幅度、相位，零极点图如图5.3所示。幅度（a）相位（b）零极点（c）图5.3幅度、相位，零极点图在切比沃雪夫I低通滤波器环境下，设置通带截止频率与阻带截止频率分别为200，100；得到幅度、相位，零极点图如图5.4所示。幅度（a）相位（b）零极点（c）图5.4幅度、相位，零极点图点击选择音频文件按钮，会出现选择文件的信息框，如图5.4。图5.4读取音频文件选择音频文件李栋.wav，先播放原音频文件，然后播放滤波后的音频文件，二者一对比明显可以听出滤波后音频文件的效果了，高音部分是被滤除的，只保留了满足低通条件的那部分音频，听起来声音低沉，还有点延时。

总 结这次课程设计是相当纠结的，因为我们还有考试，时间对我来说是宝贵的不行了，每天都是两点多睡觉，看书的同时还要把课程设计的界面，代码，以及原理进行详细的阅读分析。不过虽然是疲惫的，但是还是很充实的。而就是前一个星期，我们刚考完数字信号处理，所以对数字滤波器的设计还是很了解的，不管是脉冲响应不变法还是双线性法，我对原理部分还是比较懂的，设计的步骤以及关键技术如何用MATLAB语言实现我也都可以独立完成，但是本次是课程设计是一个综合的滤波器，不仅有巴特沃斯滤波器，同时还有切比雪夫滤波器，但是原理大致还是相通的。而书上都已经给出了用脉冲响应不变法以及双线性法的设计步骤，我们要做的就是把代码同我们设计的滤波器界面进行一一对应，这个工作是需要耐心的，因为只要有控件对应不上，那么这个功能就实现不了。在这次的数字滤波器是设计中，我们学到了课上老师没怎么讲的关于语音信号的滤波功能，其实这也只是数字滤波器的其中一个小小的功能，其实滤波器还可以对图片，视频信号等进行滤波，经过处理过后的信号才是我们想要的。就我们本次设计的IIR数字滤波器，我们可以知道，IIR滤波器的优点就是阶数少，也就是说延时少，还可以借助成熟的滤波器的成功，比如巴特沃斯，椭圆滤波器等，但同时不可避免的，IIR滤波器也有着它自己的缺点，主要就是非线性相位，易出现震荡，不稳定，主要是因为用的递归计算。这次课程设计由于时间的原因，我做的还是不怎么好的，其中对于代码的研究也没那么好，很多的用法，调用等，我都是照着课程设计指导书打上去的，也没多问几个为什么，这是我以后在学习上需要注意的地方。数字信号这门课虽然以及结束了，但是以后还将会结束很多关于数字信号处理的地方，都需要我们掌握扎实的基本功，这样不管遇到什么问题，只要我们把原理搞清楚，以及问题的思想，我们就可以运用自己所学的知识解决这些问题。

致 谢在这次的课程设计中，我首先要感谢淮阴工学院，计算机工程学院给我这次十分宝贵的实践机会，以及为我们准备实验而勤勤恳恳工作的机房工作人员，还有我的那帮乐于助人的同学们，没有他们的帮助我就没办法按时完成课程设计。课程设计的时间是短暂的，只有短短的