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

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

2、3.1 1)一阶后向差分定义为f (k) f (k) f (k 1)(3.1 2)式中和称为差分算子。由式(3.1 1)和式(3.1 2)可见，前向差分 与后向差分的关系为f (k) f (k 1)(3.1 3)二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采 用后向差分，并简称其为差分由查分的定义，若有序列fi(k)、f2(k)和常数ai， a2则a,(k) a2f2(k) a,(k) a2f2(k) a,(k 1) &2彳2化 1)afi(k) fi(k 1) &2比化)f2(k 1)(3.1 4)ai fi(k) a? f2(k)这表明差分运算具有线性性

3、质。二阶差分可定义为2f(k) f(k)f(k) f(k 1) f(k) f(k 1)(3f (k) 2f (k 1) f (k 2)类似的，可定义三阶、四阶、n阶差分。一般地，n阶差分f(k)n1f(k)n(1)jj)(3.1 6)式中n n!,j 0,1,2,L ,n (3.1 7) j (n j)! j!为二项式系数序列f(k)的求和运算为kf(i) (3.1 8)i差分方程是包含关于变量k的未知序列y(k)及其各阶差分的方程式，它的一 般形式可写为F k,y(k), y(k)丄，ny(k) 0(3.1 9a)-可编辑修改 -式中差分的最高阶为 n 阶，称为 n 阶差分方程。由式( 3.

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

5、描述该系统激励f(k)与响应y(k)之间关系的数学模型式n阶常系数线性差 分方程，它可写为3.1 10a)bmf(k) bm 1f (k 1) Lb0f (k m)y(k) an 1y(k 1) La0y(k n)式中 aj( j 0,1,L ,n 1)、bi(i 0,1,L ,m) 都是常数。上式可缩写为an j y(kj0j)am i f (ki0i)(式中 an=1)3.110b)与微分方程的经典解类似，上述差分方程的解由齐次解和特解两部分组成齐次解用yh(k)表示，特解用yp(k)表示，即y(k)=yh(k) yp(k)(3.111)a. 齐次解当式(3.1 10)中的f(k)及其各移

6、位项均为零时，齐次方程y(k) an1y(k 1) La0y(k n) 0(3.112)-可编辑修改 -的解称为齐次解。首先分析最简单的一阶差分方程。若一阶差分方程的齐次方程为(3.1 13)y(k) ay(k 1)0它可改写为y(k)y(k 1)y(k)与y(k -1)之比等于一a表明，序列y(k)是一个公比为一a的等比级数，因此y(k)应有如下形式y(k)C( a)k (3.1 14)式中C式常数，有初始条件确定。k对于n阶齐次差分方程，它的齐次解由形式为 C的序列组合而成，将kC 代入到式(3.1 12 )，得an 1Ca1C k n 1 a0C k n 0-可编辑修改 -k n由于cm

7、)，消去c；且入用，以 除上式，得n 1an 1a1ao(3.1 15)上式称为差分方程式(3.1 10)和式(3.1 12)的特征方程，它有n个根 j(j 0,1丄，n)，称为差分方程的特征根。显然，形式为Cj j的序列都满足式 (3.1 12)，因而它们是式(3.1 10)方程的齐次解。依特征根取值的不同，差分方程齐次解的形式见表3 1,其中Cj、Dj、Aj、j等为待定常数表3 1不同特征根所对应的齐次解特征根齐次解yh(k)单实根C kr重实根(Cr 1kr 1 Cr 2kr 2 L C1k C0) k一对共轭复根1,2 a jb e JkCcos( k) Dsin( k)或A kcos

8、( k)其中AeJC JDr重共轭复跟k Av 1kr 1 cos( k r 1) Ar 2kr 2 cos( k r 2) LAo cos( k0)b. 特解特解的函数形式与激励的函数形式有关，表 3 2列出了集中典型的激励f(k)所 对应的特解卩(灯。选定特解后代入原差分方程，求出其待定系数P(A或)等, 就得出方程的特解。表3 2不同激励所对应的特解激励特解yp(k)kmPmkPm lkL pk Po所有特征根均不等于1时PF Pm*"' L Pk Po当有r重等于1时的特征根时k a_. kPa当a不等于特征根时k(Pk Po)a当a是特征单根时Prk Pr1kL P

9、k Poa当a是r重特征根时cos( k)或 sin( k)Pcos( k) Qsin( k)或Acos( k ),其中AeJ P JQ所有特征根均不等于 e Jc. 全解式(3.1 10 )的线性差分方程的全解是齐次解与特解之和。如果方程的特征 根均为单根，则差分方程的全解为y(k)=yh(k) yp(k)n kCj j yp(k)1( 3.1 16)如果特征根1为r重根，而其余n - r个特征根为单根时，差分方程的全解为ry(k)= Ckr j kj 1nCj k yp(k)r 1(3.1 17)式中各系数Cj由初始条件确定如果激励信号是在k=0时接入的，差分方程的解适合于 kX)。对于n

10、阶差 分方程，用给定的n个初始条件y(0),y(1),y(n - 1)就可确定全部待定系数。 如果差分方程的特解都是单根，则方程的全解为式(3.1 16)，将给定的初始 条件y(0),y(1),y(n -1)分别代入到式(3.1 16)，可得y (0)C1 C2 LCn yp(0)y(1)1C12C2LnCn yp(1)L L L L(3.1 18 )y(n 1)1n1C2 1C2 L ；1Cn yp(n 1)由以上方程可求得全部待定系数Cj(j 0,1丄，n)2.1零输入响应系统的激励为零，仅由系统的初始状态引起的响应，称为零输入响应，用 yzi(k)表示。在零输入条件下，式(3.1 10)

11、等号右端为零，化为齐次方程, 即an jYzi(kj 0j)(3.1 25)般设定激励是在k=0时接入系统的，在k v 0时，激励尚未接入，故式(3.125)的几个初始状态满足yzi( yzi( LL1)2) LLy(y(1)2)3.1 26)yzi(n)y(n)式( 3.126)25)和式( 3.1 26)可求得零输入响应 yzi(k)。中的y( 1),y( 2),y( n)为系数的初始状态，由式(3.12.2 零状态响应当系统的初始状态为零， 仅由激励 f(k) 所产生的响应， 称为零状态响应， 用 表 示。在零状态情况下，式( 3.1 10)仍是非齐次方程，其初始状态为零，即零 状态响应

12、满足nman j yzs(kj)am i f (k i)j 0i 0yzs( 1) yzs( 2) Lyzs( n) 0 (3.1 30)的解。若其特征根均为单根，则其零状态响应为nyzs(k)Czsj jkj1yp(k)3.1 31 )式中Czsj为待定常数，yp(k)为特解。yzs( 1),yzs( 2),L ,yzs( n)为零为零,需要指出，零状态响应的初始状态但其初始值 yzs(0), yzs(1),L ,yzs(n 1)不一定等于零。3.线性常系数差分方程3.1 一个 N阶线性常系数差分方程可用下式表示：My(n)bix(ni0i)Nai y(n i)i11.4.1 )或者NMai

13、 y(n i)bix(n i),a0 10i0(1.4.2)式中，X (n)和y(n)分别是系统的输入序列和输出序列，ai和bi均为常系 数，式中 y(n-i) 和 x(n-i) 项只有一次幂，也没有相互交叉相乘项，故称为线性常 系数差分方程。差分方程的阶数是用方程 y(n-i) 项中 i 的最大取值与最小取值之 差确定的。在(1.4.2)式中，y(n-i)项i最大的取值N , i的最小取值为零，因此 称为 N 阶差分方程。4 线性常系数差分方程的求解已知系统的输入序列， 通过求解差分方程可以求出输出序列。 求解差分方程 的基本方法有以下三种：(1 ) 经典解法。这种方法类似于模拟系统中求解微

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

15、.4.1 )式，求 n 时刻的输出，要知道 n 时刻以及 n 时刻以前的输入-可编辑修改 -序列值，还要知道 n 时刻以前的 N 个输出序列值。因此求解差分方程在给定输 入序列的条件下，还需要确定 N 个初始条件。如果求 n0 时刻以后的输出， n0 时刻以前 N 个输出值 y(n0-1 )、 y( n0-2 )、 ?y?(? n、0-N )就构成了初始 条件。(1.4.1 )式表明，已知输入序列和 N 个初始条件，则可以求出 n 时刻的输出； 如果将该公式中的 n 用 n+1 代替，可以求出 n+1 时刻的输出，因此( 1.4.1 ) 式表示的差分方程本身就是一个适合递推法求解的方程。三设计

16、过程1. 用 MATLAB 求解差分方程MATLAB 信号处理工具箱提供的 filter 函数实现线性常系数差分方程的递 推求解，调用格式如下：yn=filter(B,A.xn) 计算系统对输入信号向量 xn 的零状态响应输出信号向 量 yn ，yn 与 xn 长度相等，其中， B 和 A 是( 1.4.2 )式所给差分方程的系数向 量，即B=b0,b1, ?,b?M, A=a0,a1, ?,a?N其中a0=1，如果a0工1，贝U filter用a0对系数向量B和A归一化。yn=filter(B,A.xn,xi) 计算系统对输入信号向量 xn 的全响应输出信号 yn。 所谓全响应，就是由初始状

17、态引起的零输入响应和由输入信号 xn 引起的零状态 响应之和。其中， xi 是等效初始条件的输入序列，所以 xi 是由初始条件确定的。 MATLAB 信号处理工具箱提供的 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。例1.4.1的MATLAB求解程序ep141.m 如下:%1.4.1.m:调用 MATLAB

18、解差分方程 y(n)-0.8y(n-1)=x(n)a=0.8 ; ys=1;y(-1)=1%设差分方程系数a=0.8 ,初始状态:xn=1,zeros(1,30);%x(n)=单位脉冲序列，长度N=31B=1;A=1,-0.8;%差分方程系数xi=filtic(B,A,ys); 列xi%由初始条件计算等效初始条件的输入序yn=filter(B,A,x n,xi);y(n)%调用fiter解差分方程，求系统输出信号n=O:le ngth(y n)-1;stem( n,yntitle('时域波形图)');xlabel('n');ylabel('y(n)

19、9;)程序中取查分方程系数a=0.8时，得到系统输出y (n)如图1.4.1 (a)所 示，与例1.4.1的解析递推结果完全相同。如果令初始条件 y (-1)=0 (仅修改 程序中ys=0 )，则得到系统输出y (n) =h (n)，如图1.4.1 (b)所示。n-可编辑修改-(b)(a)图(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,ze

20、ros(1,30);%x(n)= 单位脉冲序列，长度 N=31B=1;A=1,-0.8;%差分方程系数xi=filtic(B,A,ys);% 由初始条件计算等效初始条件的输入序列 xiyn=filter(B,A,xn,xi);%调用filter解差分方程，求系统输出信号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.八

21、n;%xn=0.5.Anstem(n,xn,'fill'),grid on%画出时域波形图xlabel('n'),ylabel('x(n)'),title('时域波形图')%x轴、y轴分别代表n，x(n)n=-5:5;%n的取值范围a=0.5;% 设 a=0.5xen二 a*a.A n+a4( -n);%xen=a*a.An+a.A(-n)-可编辑修改 -xon二 a*a.A n-a.A( -n);%xo n=a*a.A n-a.A(- n)figure(1); stem(n,xen,'filled'),grid

22、on%画出xe(n)时域波形图title('时域波形图');xlabel('n');ylabel('xe(n)')%x轴、y轴分别代表 n , xe(n)figure(2); stem(n,xon,'filled'),grid 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

23、(-1)=1B=2,0,-1;A=1,-0.7,0.1;%差分方程系数n=-5:5;%n的取值范围xn=0.5.An;%x(n)=0.5.Anxi=filtic(B,A,ys);% 由初始条件计算等效初始条件的输入序列 xiyn=filter(B,A,xn,xi);%调用fiter解差分方程，求系统输出信号y(n)stem(n,yn,'.')%画出时域波形图title('(a)');xlabel('n');ylabel('y(n)') %x轴、y轴分别代表 n , y(n)五程序运行结果如下：-可编辑修改 -可编辑修改-差分序列时

24、域波形图得到二阶线性常系数差分输入参数 a仁-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,得到结果如下:输入初始条件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

25、(n)波形图六.比较结果总结由上实验可知，通过改变初始条件ys的值，得到的输出波形大小并不一致, 即输出信号y(n)是不相同的；从而我们可以得出，对于同一个差分方程和同一个 输入信号，因为初始条件不同，得到的输出信号是不相同的。七.收获与体会本次MATLAB课程设计让我熟悉了该软件的一些功能，但是对于灵活应用 MATLAB，以及掌握各方面的设计思维以及技巧，还需要投入更多的时间。在熟悉MATLAB程序和操作的同时培养了我的独立思考能力，专研精神，解决冋题 能力和动手能力。在此之前了解到MATLAB是一个很重要很有用的工具，但我并没有完全理 解，本课程设计中，通过查阅资料，阅读网上程序并读写程序

26、，对于MATLAB的应用有了更深的了解，同时也认识到 MATLAB功能非常的强大，有着很多方 面的应用，如绘制函数，处理音频，图像数据，创建用户界面等功能，实为一个 功能强大的软件。本次课程设计我完成了基于 MATLAB 的线性常系数差分方程求解的题目， 通过实际操作回顾所学的内容， 强化基础，实践理论知识。 相信在以后的学习中， 还会更加深入的了解 MATLAB ，应用它。随着课程设计报告的基本完成， 本次课程设计终于接近了尾声。 本次课程设 计要求我们利用上学期所学的信号与线性系统分析的知识结合 MATLAB 编程工 具，完成差分方程求解设计的题目，通过实际操作，回顾所学内容，务实基础，

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