数字信号处理(DSP)课程设计—利用Matlab设计IIR数字滤波器

数字信号处理课程设计报告姓名：一.蒲钇霖二.舒琳学号：一.201021030619二.201021030414学院：微固利用Matlab设计IIR数字滤波器一、前言数字化、智能化和网络化是当代信息技术发展的大趋势，而数字化是智能化和网络化的基础，这也就促进了数字信号处理这门技术的飞速发展，使其不仅自成一门学科，更是以不同形式影响和渗透到其它学科。如今，数字信号处理技术与国民经济息息相关，与国家建设紧密相连，影响着并改变着我们的生产、生活方式，因此它受到了社会的普遍关注。而在数字信号处理当中，数字滤波器又占有相当重要的位置。因为任何需要我们处理的信号都包含有噪声，而滤波是去除噪声的基本手段。数字滤波器简单的说就是一种用来过滤时间离散信号的线性时不变系统，通过对采样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可以分为两类：无限冲激响应(IIR)滤波器和有限冲激响应（FIR）滤波器。数字滤波本质上是一种运算过程，输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器的本质上是一个完成特定运算的数字计算过程。描述离散系统输入和输出关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。二、数字滤波器设计原理在数字滤波过程中，我们主要讨论离散时间序列。设输入序列为x(n)，滤波器对单位采样序列δ(n)的响应为h(n)。因δ(n)在时域离散信号和系统中所起的作用相当于单位冲激函数在时域连续信号和系统中所起的作用。所以数字滤波器的序列y(n)就是以上两个序列的卷积，即：()()()kynhkxnk(1)同样，两个序列卷积的z变换等于各自z变换的乘积，即：()()()YzHzXz(2)用jze代入上式，则可得到时域离散系统的频谱特性：()()()jjjYeHeXe(3)式中()jYe和()jXe分别是数字滤波器的输出序列和输入序列的频谱特性，而()jHe则是单位采样序列h(n)的频特性谱，又称为数字滤波器的频域响应。输入序列()jXe经过滤波后变为()jXe()jHe。因此，只需要按照输入信号频谱()jXe的特点和我们处理信号的目的，选择适当的频域响应()jHe，使得滤波后的输出信号频谱()jYe满足我们设计的要求。三、课程设计目的与内容通过完成本次课程设计，加深我们对于在数字信号处理这门学科的印象，巩固我们在课堂上学习到的理论知识和实验方法，增强我们将理论知识化为实践技巧的能力，并且让我们进一步熟悉数字信号处理中常用软件Matlab的使用方法。在这个过程中同时培养我们发现问题、分析问题以及解决问题的能力。在这次课程设计当中，我们需要分析和研究采用不同设计方法时，滤波器最低阶数与设计指标之间的关系，通过理论分析方式推导建立起相应的简单数学模型。在此基础上，自选设计指标，分别采用双线性变换法和冲激响应不变法进行设计Butterworth低通滤波器，然后通过仿真检验所建模型的正确性和局限性。四、具体设计步骤1.推导滤波器最小阶数与设计指标的关系滤波器阶数就是其系统函数的极点个数。为了避免滤波器的复杂程度与我们的要求不匹配，造成不必要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设计要求的最小阶数。下面就分别以Butterworth滤波器和Chebyshev1型滤波器为例，推导其阶数的数学模型。①对于Butterworth滤波器，若给出滤波器的设计指标为：通带截止频率p，阻带截止频率s，通带最大纹波()pRdB，阻带最小纹波()sRdB。因为滤波器幅频特性为：21(1()NcHj（4）其中c为3dB截止频率，N为滤波器阶数。所以当p以及s时，可得到：2211((1()1()psNpsNccHjHj、(5)然后由()pRdB与()sRdB的关系式可得到：221022101110lg()(101110lg()(10psppRpssRsRHjHjRHjHj(6)由上面的(5)式和(6)式可以联立求得：21010102101()10101()1011()10ppssRpNRpcNRRsNsc(7)通过上面的结果我们就可以求得滤波器阶数N为：1010101lg101lg()psRRpsN(8)所以滤波器的最小阶数就是大于上式所求得的值的最小整数。②对于Chebyshev1型滤波器，推导方式和上面类似，同样的给出滤波器的设计指标为：通带截止频率p，阻带截止频率s，通带最大纹波()pRdB，阻带最小纹波()sRdB。因为滤波器幅频特性为：221()1()cNHjV(9)又因为通带内起伏的等纹波参数/10101pR，系数2010pRA，所以可以得出滤波器阶数为：22lg(1)lg(()1)ssppggN(10)式中22(1)/gA。滤波器的最小阶数就是取大于上式的值的最小整数。2.双线性变换法设计Butterworth低通滤波器双线性变换法是在完成模拟系统设计后，再进行反变换，从模拟系统函数得到数字系统函数。它没有采样过程，不存在频率混叠问题，适应于各类滤波器的变换，但是变换在高频区域会出现强烈的非线性。下面就是基于Butterworth模拟原型滤波器，使用双线性变换法设计数字低通滤波器的过程。设计指标为：通带截止频率p=0.2π,阻带截止频率s=0.3π，通带最大纹波Rp=1dB，阻带最大纹波Rs=15dB，采样频率Fs=20000Hz。同时，为了检测滤波器的性能，我们假定一个输入信号21)()sin(2)0.5cos(2txtftf，其中f1=1000Hz，f2=4000Hz，将该信号与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。Matlab程序如下：Wp=0.2*pi;Ws=0.3*pi;fs=20000;T=1/fs;Rp=1;Rs=15;wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;[N,Wn]=buttord(wp,ws,Rp,Rs,'s');[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,fs);figure(1);[h,w]=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs(h));title('ButterworthLPF幅频特性');xlabel('频率(Hz)');ylabel('幅值(dB)');gridon;subplot(2,1,2),plot(w,angle(h));title('ButterworthLPF相频特性');xlabel('频率(Hz)');ylabel('相位(degree)');gridon;figure(2);f1=1000;f2=4000;N1=100;dt=1/fs;n=0:N1-1;t=n*dt;x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);title('输入信号'),gridon;y=filter(bz,az,x);subplot(2,1,2),plot(t,y,'r-');title('输出信号'),gridon;得到的Butterworth低通滤波器幅频特性和相频特性曲线为：由语句[N,Wn]=buttord(wp,ws,Rp,Rs,'s')我们可以用Matlab得出滤波器的最小阶数N=6。再利用上面我们推导出的式（8），将滤波器各项指标代入可以算出N=4.7433，取大于它的最小整数得到滤波器的最小阶数为5。二者并不一致，但差距只有1，说明我们推导出的计算Butterworth滤波器最小阶数的数学模型存在误差，但误差不大。又因为Matlab计算滤波器的最小阶数也是估算出来的，所以这个这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。输入信号和经过滤波后的输出信号图像为：3.冲激响应不变法设计Butterworth低通滤波器冲激响应不变法是对模拟滤波器的单位冲激响应h(t)等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的系统函数。这次我们设定与上面不同的设计指标：通带截止频率p=2000π，阻带截止频率s=3000π，通带最大纹波Rp=3dB，阻带最小纹波Rs=18dB,采样频率Fs=10000Hz。为了方便比较，我们设定一样的输入信号21)()sin(2)0.5cos(2txtftf。Matlab程序如下：wp=2000*pi;ws=3000*pi;fs=10000;Rp=3;Rs=18;[N,wn]=buttord(wp,ws,Rp,Rs,'s');[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,wn);[bz,az]=impinvar(b,a,fs);figure(1);[h,w]=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs(h));title('ButterworthLPF幅频特性');xlabel('频率(Hz)');ylabel('幅值(dB)');gridon;subplot(2,1,2),plot(w,angle(h));title('ButterworthLPF相频特性');xlabel('频率(Hz)');ylabel('相位(degree)');gridon;figure(2);f1=1000;f2=4000;N1=100;dt=1/fs;n=0:N1-1;t=n*dt;x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);title('输入信号'),gridon;y=filter(bz,az,x);subplot(2,1,2),plot(t,y,'r-');title('输出信号'),gridon;得到的Butterworth低通滤波器幅频特性和相频特性曲线为：通过比较两组幅频特性和相频特性曲线我们可以发现，其实双线性变换法和冲激响应不变法在设计较为简单的Butterworth低通滤波器时差别并不大，并且都拥有不错的设计效果。同样的，我们通过语句[N,wn]=buttord(wp,ws,Rp,Rs,'s')可以得到Matlab估算出的滤波器最小阶数为N=6，和上一个滤波器一样。用式(8)出来的N=5.1925，取大于它的最小整数得到滤波器的最小阶数为6，两者一致，再次验证了该数学模型的可靠性。输入信号和经过滤波后的输出信号图像为：4.双线性变换法设计Chebyshev1型低通滤波器为了和Butterworth滤波器进行对比，我们采用和双线性变换法设计Butterworth低通滤波器时一样的设计指标：通带截止频率p=0.2π,阻带截止频率s=0.3π，通带最大纹波Rp=1dB，阻带最大纹波Rs=15dB，采样频率Fs=20000Hz。输入信号同样为：21)()sin(2)0.5cos(2txtftf。Matlab程序如下：Wp=0.2*pi;Ws=0.3*pi;fs=20000;T=1/fs;Rp=1;Rs=15;wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');[z,p,k]=cheb1ap(N,Rp);[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,fs);figure(1);[h,w]=freqz(bz,az,N,fs);subplot(2,1,1),pl