用matlab实现线性常系数差分方程的求解

1、数字信号处理课程设计题目： 试实现线性常系数差分方程的求解 学院： 专业： 班级： 学号： 组员： 指导教师： 题目：用Matlab实现线性常系数差分方程求解一 设计要求1 掌握线性常系数差分方程的求解2 熟练掌握Matlab基本操作和各类函数调用3 结合Matlab实现线性常系数差分方程的求解二设计原理1差分与差分方程与连续时间信号的微分及积分运算相对应，离散时间信号有差分及序列求和运算。设有序列f(k),则称,f(k+2),f(k+1),f(k1),f(k2),为f(k)的移位序列。序列的差分可以分为前向差分和后向差分。一阶前向差分定义为 （3.11）一阶后向差分定义为 （3.12）式中和

2、称为差分算子。由式（3.11）和式（3.12）可见，前向差分与后向差分的关系为 （3.13）二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采用后向差分，并简称其为差分。由查分的定义，若有序列、和常数，则 （3.14）这表明差分运算具有线性性质。 二阶差分可定义为 （3.15）类似的，可定义三阶、四阶、n阶差分。一般地，n阶差分 （3.16）式中 （3.17）为二项式系数序列f(k)的求和运算为 （3.18）差分方程是包含关于变量k的未知序列y(k)及其各阶差分的方程式，它的一般形式可写为 （3.19a）式中差分的最高阶为n阶，称为n阶差分方程。由式（3.16）可知，各阶差分均

3、可写为y(k)及其各移位序列的线性组合，故上式常写为 （3.19b）通常所说的差分方程是指式（3.19b）形式的方程。若式（3.19b）中，y(k)及其各移位序列均为常数，就称其为常系数差分方程；如果某些系数是变量k的函数，就称其为变系数差分方程。描述LTI离散系统的是常系数线性差分方程。差分方程是具有递推关系的代数方程，若一直初始条件和激励，利用迭代法渴求的差分方程的数值解。2. 差分方程的经典解一般而言，如果但输入单输出的LTI系统的激励f(k)，其全响应为y(k)，那么，描述该系统激励f(k)与响应y(k)之间关系的数学模型式n阶常系数线性差分方程，它可写为 （3.110a）式中、都是常

4、数。上式可缩写为 （3.110b）与微分方程的经典解类似，上述差分方程的解由齐次解和特解两部分组成。齐次解用表示，特解用表示，即 （3.111）a.齐次解当式（3.110)中的f(k)及其各移位项均为零时，齐次方程 （3.112）的解称为齐次解。首先分析最简单的一阶差分方程。若一阶差分方程的齐次方程为 （3.113）它可改写为y(k)与y(k1)之比等于a表明，序列y(k)是一个公比为a的等比级数，因此y(k)应有如下形式 （3.114）式中C式常数，有初始条件确定。对于n阶齐次差分方程，它的齐次解由形式为的序列组合而成，将代入到式（3.112），得由于C0，消去C；且0，以除上式，得（3.1

5、15）上式称为差分方程式（3.110）和式（3.112）的特征方程，它有n个根，称为差分方程的特征根。显然，形式为的序列都满足式（3.112），因而它们是式（3.110）方程的齐次解。依特征根取值的不同，差分方程齐次解的形式见表31,其中、等为待定常数表31 不同特征根所对应的齐次解特征根齐次解单实根重实根一对共轭复根重共轭复跟b.特解特解的函数形式与激励的函数形式有关，表32列出了集中典型的激励f(k)所对应的特解。选定特解后代入原差分方程，求出其待定系数等，就得出方程的特解。表32 不同激励所对应的特解激励特解 所有特征根均不等于1时 当有重等于1时的特征根时 当不等于特征根时 当是特征单

6、根时 当是重特征根时或所有特征根均不等于c.全解 式（3.110）的线性差分方程的全解是齐次解与特解之和。如果方程的特征根均为单根，则差分方程的全解为（3.116）如果特征根为重根，而其余nr个特征根为单根时，差分方程的全解为（3.117）式中各系数由初始条件确定。如果激励信号是在k=0时接入的，差分方程的解适合于k0。对于n阶差分方程，用给定的n个初始条件y(0),y(1),y(n1)就可确定全部待定系数。如果差分方程的特解都是单根，则方程的全解为式（3.116），将给定的初始条件y(0),y(1),y(n1)分别代入到式（3.116），可得（3.118）由以上方程可求得全部待定系数。2.1

7、零输入响应系统的激励为零，仅由系统的初始状态引起的响应，称为零输入响应，用表示。在零输入条件下，式（3.110）等号右端为零，化为齐次方程，即 （3.125）一般设定激励是在k=0时接入系统的，在k0时，激励尚未接入，故式（3.125）的几个初始状态满足 （3.126）式（3.126）中的y(1),y(2),y(n)为系数的初始状态，由式（3.125）和式（3.126）可求得零输入响应。2.2零状态响应当系统的初始状态为零，仅由激励f(k)所产生的响应，称为零状态响应，用 表示。在零状态情况下，式（3.110）仍是非齐次方程，其初始状态为零，即零状态响应满足（3.130）的解。若其特征根均为单

8、根，则其零状态响应为 （3.131）式中为待定常数，为特解。需要指出，零状态响应的初始状态为零，但其初始值不一定等于零。3.线性常系数差分方程3.1一个N 阶线性常系数差分方程可用下式表示： （）或者 （）式中，x（n）和y(n)分别是系统的输入序列和输出序列，ai和bi均为常系数，式中y(n-i)和x(n-i)项只有一次幂，也没有相互交叉相乘项，故称为线性常系数差分方程。差分方程的阶数是用方程y(n-i)项中i的最大取值与最小取值之差确定的。在（）式中，y(n-i)项i最大的取值N，i的最小取值为零，因此称为N阶差分方程。4 线性常系数差分方程的求解已知系统的输入序列，通过求解差分方程可以求

9、出输出序列。求解差分方程的基本方法有以下三种：（1） 经典解法。这种方法类似于模拟系统中求解微分方程的方法，它包括齐次解与特解，由边界条件求待定系数，上节已作简单介绍，这里不作介绍。（2） 递推解法。这种方法简单，且适合用计算机求解，但只能得到数值解，对于阶次较高的线性常系数差分方程不容易得到封闭式（公式）解答。（3） 变换域方法。这种方法是将差分方程变换到z域进行求解，方法简便有效。当然还可以不直接求解差分方程，而是先由差分方程求出系统的单位脉冲响应，再与已知的输入序列进行卷积运算，得到系统输出。但是系统的单位脉冲响应如果不是预先知道，仍然需要求解差分方程，求其零状态响应解。（4） 卷积法：

10、由差分方程求出系统的h(n)，再与已知的x(n) 进行卷积，得到y(n)。观察（）式，求n时刻的输出，要知道n时刻以及n时刻以前的输入序列值，还要知道n时刻以前的N个输出序列值。因此求解差分方程在给定输入序列的条件下，还需要确定N个初始条件。如果求n0时刻以后的输出，n0时刻以前N个输出值y（n0-1）、y（n0-2）、y（n0-N）就构成了初始条件。（）式表明，已知输入序列和N个初始条件，则可以求出n时刻的输出；如果将该公式中的n用n+1代替，可以求出n+1时刻的输出，因此（1.4.1）式表示的差分方程本身就是一个适合递推法求解的方程。三设计过程1.用MATLAB求解差分方程 MATLAB信

11、号处理工具箱提供的filter函数实现线性常系数差分方程的递推求解，调用格式如下： yn=filter(B,A.xn) 计算系统对输入信号向量xn的零状态响应输出信号向量yn，yn与xn长度相等，其中，B和A是（）式所给差分方程的系数向量，即 B=b0,b1,bM, A=a0,a1,aN其中a0=1，如果a01，则filter用a0对系数向量B和A归一化。 yn=filter(B,A.xn,xi) 计算系统对输入信号向量xn的全响应输出信号yn。所谓全响应，就是由初始状态引起的零输入响应和由输入信号xn引起的零状态响应之和。其中，xi是等效初始条件的输入序列，所以xi是由初始条件确定的。MAT

12、LAB信号处理工具箱提供的filtic就是由初始条件计算xi的函数，其调用格式如下：xi=filtic（B,A,ys,xs）其中，ys和xs是初始条件向量：ys=y(-1),y(-2),y(-3),y(-N),xs=x(-1),x(-2),x(-3),x(-M)。如果xn是因果序列，则xs=0.调用时可缺省xs。 例的MATLAB求解程序ep141.m如下：%1.4.1.m:调用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)a=0.8；ys=1; %设差分方程系数a=0.8，初始状态：y(-1)=1xn=1,zeros(1,30); %x(n)=单位脉冲序列，长度N=31B=1;

13、A=1,-0.8; %差分方程系数xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xiyn=filter(B,A,xn,xi); %调用fiter解差分方程，求系统输出信号y(n)n=0:length(yn)-1;stem(n,yn,'.')title('时域波形图)');xlabel('n');ylabel('y(n)')程序中取查分方程系数a=0.8时，得到系统输出y（n）如图（a）所示，与例1.4.1的解析递推结果完全相同。如果令初始条件y（-1）=0（仅修改程序中ys=0），则得到系统输出y（n

14、）=h（n），如图1.4.1（b）所示。(a) (b)图（a）为a=0.8，y（-1）=1时，系统输出时域波形图，图（b）为a=0.8，y(-1)=0时，系统输出时域波形图。 四设计代码及结果MATLAB源程序源程序如下%1.m:调用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)ys=1; %初始状态：y(-1)=1xn=1,zeros(1,30); %x(n)=单位脉冲序列，长度N=31B=1;A=1,-0.8; %差分方程系数xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xiyn=filter(B,A,xn,xi); %调用filter解差分方程

15、，求系统输出信号y(n)n=0:length(yn)-1; %n的取值范围stem(n,yn,'.') %画出时域波形图title('时域波形图)');xlabel('n');ylabel('y(n)') %x轴、y轴分别代表n，x(n)n=-5:5; %n的取值范围xn=0.5.n; %xn=0.5.nstem(n,xn,'fill'),grid on %画出时域波形图xlabel('n'),ylabel('x(n)'), title('时域波形图') %x轴、y轴

16、分别代表n，x(n)n=-5:5; %n的取值范围a=0.5; %设a=0.5xen=a*a.n+a.(-n); %xen=a*a.n+a.(-n)xon=a*a.n-a.(-n); %xon=a*a.n-a.(-n)figure(1); stem(n,xen,'filled'),grid on %画出xe(n)时域波形图 title('时域波形图');xlabel('n');ylabel('xe(n)') %x轴、y轴分别代表n，xe(n)figure(2); stem(n,xon,'filled'),grid

17、on %画出xo(n)时域波形图 title('时域波形图');xlabel('n');ylabel('xo(n)') %x轴、y轴分别代表n，xo(n)b0=2;b2=0;b2=-1;a1=-0.7;a2=0.1;ys=0; %设差分方程系数，初始状态：y(-1)=1B=2,0,-1;A=1,-0.7,0.1; %差分方程系数n=-5:5; %n的取值范围xn=0.5.n; %x(n)=0.5.nxi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xiyn=filter(B,A,xn,xi); %调用fiter解差分方程

18、，求系统输出信号y(n)stem(n,yn,'.') %画出时域波形图title('(a)');xlabel('n');ylabel('y(n)') %x轴、y轴分别代表n，y(n)五程序运行结果如下： 差分序列时域波形图输入参数a1=-0.7,a2=0.1,b0=2,b1=0,b2=-1 得到二阶线性常系数差分方程为 y(n)-0.7y(n-1)+0.1y(n-2)=2x(n)-x(n-2)结果如下：x（n）的时域波形图 共轭对称分量xe（n）的时域波形图 共轭对称分量xo（n）的时域波形图 输入初始条件ys=y(-1)=1,得

19、到结果如下：输入初始条件ys=y(-1)=1时，输出y(n)波形图改变初始条件ys=y(-1)=0,得到结果如下：输入初始条件ys=y(-1)=0时，输出y(n)波形图 改变初始条件ys=y(-1)=100,得到结果如下: 输入初始条件ys=y(-1)=100时，输出y(n)波形图 六比较结果总结由上实验可知，通过改变初始条件ys的值，得到的输出波形大小并不一致，即输出信号y(n)是不相同的；从而我们可以得出，对于同一个差分方程和同一个输入信号，因为初始条件不同，得到的输出信号是不相同的。七收获与体会 本次MATLAB课程设计让我熟悉了该软件的一些功能，但是对于灵活应用MATLAB，以及掌握各

20、方面的设计思维以及技巧，还需要投入更多的时间。在熟悉MATLAB程序和操作的同时培养了我的独立思考能力，专研精神，解决问题能力和动手能力。在此之前了解到MATLAB是一个很重要很有用的工具，但我并没有完全理解，本课程设计中，通过查阅资料，阅读网上程序并读写程序，对于MATLAB的应用有了更深的了解，同时也认识到MATLAB功能非常的强大，有着很多方面的应用，如绘制函数，处理音频，图像数据，创建用户界面等功能，实为一个功能强大的软件。本次课程设计我完成了基于MATLAB的线性常系数差分方程求解的题目，通过实际操作回顾所学的内容，强化基础，实践理论知识。相信在以后的学习中，还会更加深入的了解MAT

21、LAB，应用它。随着课程设计报告的基本完成，本次课程设计终于接近了尾声。本次课程设计要求我们利用上学期所学的信号与线性系统分析的知识结合MATLAB编程工具，完成差分方程求解设计的题目，通过实际操作，回顾所学内容，务实基础，强化理论知识，并体验理论与实际相结合的过程。设计过程中遇到的第一个问题便是对于MATLAB语言的不熟悉，其实现在想想这个问题不应该成为问题。毕竟本专业曾开设过MATLAB程序设计这门课，而且老师还特别提醒过课程设计会用到MATLAB语言。可惜的是，老师的话并没有引起我的足够重视，才导致要真正设计的时候一头雾水，无从下手。通过这件事，我明白了对于一件事情，想要做得很好，提前做功课和准备是必不可少的，机会只垂青有准备的人。当然，通过本次课程设计，我还是基本熟悉了一些MATLAB模块以及与本课程有关的一些函数的用法。但是由于上学期信号与线性系统分析课程学的不是很好，所以也就是在懂得同学的指导下按部就班的写了一些代码，也没有什么创新，但还是很好的回顾了差分方程求解的内容。老师在任务书中提到的目的要求，我自认为多多少少完成了一些。但不可否认的是还有很多没有达到的。总之，很多东西还是需要自己在不断摸索中找到答案。同时，通过本次数字信号处理课程设计，使我全面系统的了解了差分方程求解的原理。通过结合课本的知识去完成课程设