用matlab實(shí)現(xiàn)線性常系數(shù)差分方程的求解

1、數(shù)字信號處理課程設(shè)計(jì)題目： 試實(shí)現(xiàn)線性常系數(shù)差分方程的求解 學(xué)院： 專業(yè)： 班級： 學(xué)號： 組員： 指導(dǎo)教師： 題目：用Matlab實(shí)現(xiàn)線性常系數(shù)差分方程求解一 設(shè)計(jì)要求1 掌握線性常系數(shù)差分方程的求解2 熟練掌握Matlab基本操作和各類函數(shù)調(diào)用3 結(jié)合Matlab實(shí)現(xiàn)線性常系數(shù)差分方程的求解二設(shè)計(jì)原理1差分與差分方程與連續(xù)時(shí)間信號的微分及積分運(yùn)算相對應(yīng)，離散時(shí)間信號有差分及序列求和運(yùn)算。設(shè)有序列f(k),則稱,f(k+2),f(k+1),f(k1),f(k2),為f(k)的移位序列。序列的差分可以分為前向差分和后向差分。一階前向差分定義為 （3.11）一階后向差分定義為 （3.12）式中和

2、稱為差分算子。由式（3.11）和式（3.12）可見，前向差分與后向差分的關(guān)系為 （3.13）二者僅移位不同，沒有原則上的差別，因而它們的性質(zhì)也相同。此處主要采用后向差分，并簡稱其為差分。由查分的定義，若有序列、和常數(shù)，則 （3.14）這表明差分運(yùn)算具有線性性質(zhì)。 二階差分可定義為 （3.15）類似的，可定義三階、四階、n階差分。一般地，n階差分 （3.16）式中 （3.17）為二項(xiàng)式系數(shù)序列f(k)的求和運(yùn)算為 （3.18）差分方程是包含關(guān)于變量k的未知序列y(k)及其各階差分的方程式，它的一般形式可寫為 （3.19a）式中差分的最高階為n階，稱為n階差分方程。由式（3.16）可知，各階差分均

3、可寫為y(k)及其各移位序列的線性組合，故上式常寫為 （3.19b）通常所說的差分方程是指式（3.19b）形式的方程。若式（3.19b）中，y(k)及其各移位序列均為常數(shù)，就稱其為常系數(shù)差分方程；如果某些系數(shù)是變量k的函數(shù)，就稱其為變系數(shù)差分方程。描述LTI離散系統(tǒng)的是常系數(shù)線性差分方程。差分方程是具有遞推關(guān)系的代數(shù)方程，若一直初始條件和激勵(lì)，利用迭代法渴求的差分方程的數(shù)值解。2. 差分方程的經(jīng)典解一般而言，如果但輸入單輸出的LTI系統(tǒng)的激勵(lì)f(k)，其全響應(yīng)為y(k)，那么，描述該系統(tǒng)激勵(lì)f(k)與響應(yīng)y(k)之間關(guān)系的數(shù)學(xué)模型式n階常系數(shù)線性差分方程，它可寫為 （3.110a）式中、都是常

4、數(shù)。上式可縮寫為 （3.110b）與微分方程的經(jīng)典解類似，上述差分方程的解由齊次解和特解兩部分組成。齊次解用表示，特解用表示，即 （3.111）a.齊次解當(dāng)式（3.110)中的f(k)及其各移位項(xiàng)均為零時(shí)，齊次方程 （3.112）的解稱為齊次解。首先分析最簡單的一階差分方程。若一階差分方程的齊次方程為 （3.113）它可改寫為y(k)與y(k1)之比等于a表明，序列y(k)是一個(gè)公比為a的等比級數(shù)，因此y(k)應(yīng)有如下形式 （3.114）式中C式常數(shù)，有初始條件確定。對于n階齊次差分方程，它的齊次解由形式為的序列組合而成，將代入到式（3.112），得由于C0，消去C；且0，以除上式，得（3.1

5、15）上式稱為差分方程式（3.110）和式（3.112）的特征方程，它有n個(gè)根，稱為差分方程的特征根。顯然，形式為的序列都滿足式（3.112），因而它們是式（3.110）方程的齊次解。依特征根取值的不同，差分方程齊次解的形式見表31,其中、等為待定常數(shù)表31 不同特征根所對應(yīng)的齊次解特征根齊次解單實(shí)根重實(shí)根一對共軛復(fù)根重共軛復(fù)跟b.特解特解的函數(shù)形式與激勵(lì)的函數(shù)形式有關(guān)，表32列出了集中典型的激勵(lì)f(k)所對應(yīng)的特解。選定特解后代入原差分方程，求出其待定系數(shù)等，就得出方程的特解。表32 不同激勵(lì)所對應(yīng)的特解激勵(lì)特解 所有特征根均不等于1時(shí) 當(dāng)有重等于1時(shí)的特征根時(shí) 當(dāng)不等于特征根時(shí) 當(dāng)是特征單

6、根時(shí) 當(dāng)是重特征根時(shí)或所有特征根均不等于c.全解 式（3.110）的線性差分方程的全解是齊次解與特解之和。如果方程的特征根均為單根，則差分方程的全解為（3.116）如果特征根為重根，而其余nr個(gè)特征根為單根時(shí)，差分方程的全解為（3.117）式中各系數(shù)由初始條件確定。如果激勵(lì)信號是在k=0時(shí)接入的，差分方程的解適合于k0。對于n階差分方程，用給定的n個(gè)初始條件y(0),y(1),y(n1)就可確定全部待定系數(shù)。如果差分方程的特解都是單根，則方程的全解為式（3.116），將給定的初始條件y(0),y(1),y(n1)分別代入到式（3.116），可得（3.118）由以上方程可求得全部待定系數(shù)。2.1

7、零輸入響應(yīng)系統(tǒng)的激勵(lì)為零，僅由系統(tǒng)的初始狀態(tài)引起的響應(yīng)，稱為零輸入響應(yīng)，用表示。在零輸入條件下，式（3.110）等號右端為零，化為齊次方程，即 （3.125）一般設(shè)定激勵(lì)是在k=0時(shí)接入系統(tǒng)的，在k0時(shí)，激勵(lì)尚未接入，故式（3.125）的幾個(gè)初始狀態(tài)滿足 （3.126）式（3.126）中的y(1),y(2),y(n)為系數(shù)的初始狀態(tài)，由式（3.125）和式（3.126）可求得零輸入響應(yīng)。2.2零狀態(tài)響應(yīng)當(dāng)系統(tǒng)的初始狀態(tài)為零，僅由激勵(lì)f(k)所產(chǎn)生的響應(yīng)，稱為零狀態(tài)響應(yīng)，用 表示。在零狀態(tài)情況下，式（3.110）仍是非齊次方程，其初始狀態(tài)為零，即零狀態(tài)響應(yīng)滿足（3.130）的解。若其特征根均為單

8、根，則其零狀態(tài)響應(yīng)為 （3.131）式中為待定常數(shù)，為特解。需要指出，零狀態(tài)響應(yīng)的初始狀態(tài)為零，但其初始值不一定等于零。3.線性常系數(shù)差分方程3.1一個(gè)N 階線性常系數(shù)差分方程可用下式表示： （1.4.1）或者 （1.4.2）式中，x（n）和y(n)分別是系統(tǒng)的輸入序列和輸出序列，ai和bi均為常系數(shù)，式中y(n-i)和x(n-i)項(xiàng)只有一次冪，也沒有相互交叉相乘項(xiàng)，故稱為線性常系數(shù)差分方程。差分方程的階數(shù)是用方程y(n-i)項(xiàng)中i的最大取值與最小取值之差確定的。在（1.4.2）式中，y(n-i)項(xiàng)i最大的取值N，i的最小取值為零，因此稱為N階差分方程。4 線性常系數(shù)差分方程的求解已知系統(tǒng)的輸

9、入序列，通過求解差分方程可以求出輸出序列。求解差分方程的基本方法有以下三種：（1） 經(jīng)典解法。這種方法類似于模擬系統(tǒng)中求解微分方程的方法，它包括齊次解與特解，由邊界條件求待定系數(shù)，上節(jié)已作簡單介紹，這里不作介紹。（2） 遞推解法。這種方法簡單，且適合用計(jì)算機(jī)求解，但只能得到數(shù)值解，對于階次較高的線性常系數(shù)差分方程不容易得到封閉式（公式）解答。（3） 變換域方法。這種方法是將差分方程變換到z域進(jìn)行求解，方法簡便有效。當(dāng)然還可以不直接求解差分方程，而是先由差分方程求出系統(tǒng)的單位脈沖響應(yīng)，再與已知的輸入序列進(jìn)行卷積運(yùn)算，得到系統(tǒng)輸出。但是系統(tǒng)的單位脈沖響應(yīng)如果不是預(yù)先知道，仍然需要求解差分方程，求其

10、零狀態(tài)響應(yīng)解。（4） 卷積法：由差分方程求出系統(tǒng)的h(n)，再與已知的x(n) 進(jìn)行卷積，得到y(tǒng)(n)。觀察（1.4.1）式，求n時(shí)刻的輸出，要知道n時(shí)刻以及n時(shí)刻以前的輸入序列值，還要知道n時(shí)刻以前的N個(gè)輸出序列值。因此求解差分方程在給定輸入序列的條件下，還需要確定N個(gè)初始條件。如果求n0時(shí)刻以后的輸出，n0時(shí)刻以前N個(gè)輸出值y（n0-1）、y（n0-2）、y（n0-N）就構(gòu)成了初始條件。（1.4.1）式表明，已知輸入序列和N個(gè)初始條件，則可以求出n時(shí)刻的輸出；如果將該公式中的n用n+1代替，可以求出n+1時(shí)刻的輸出，因此（1.4.1）式表示的差分方程本身就是一個(gè)適合遞推法求解的方程。三設(shè)計(jì)

11、過程1.用MATLAB求解差分方程 MATLAB信號處理工具箱提供的filter函數(shù)實(shí)現(xiàn)線性常系數(shù)差分方程的遞推求解，調(diào)用格式如下： yn=filter(B,A.xn) 計(jì)算系統(tǒng)對輸入信號向量xn的零狀態(tài)響應(yīng)輸出信號向量yn，yn與xn長度相等，其中，B和A是（1.4.2）式所給差分方程的系數(shù)向量，即 B=b0,b1,bM, A=a0,a1,aN其中a0=1，如果a01，則filter用a0對系數(shù)向量B和A歸一化。 yn=filter(B,A.xn,xi) 計(jì)算系統(tǒng)對輸入信號向量xn的全響應(yīng)輸出信號yn。所謂全響應(yīng)，就是由初始狀態(tài)引起的零輸入響應(yīng)和由輸入信號xn引起的零狀態(tài)響應(yīng)之和。其中，xi

12、是等效初始條件的輸入序列，所以xi是由初始條件確定的。MATLAB信號處理工具箱提供的filtic就是由初始條件計(jì)算xi的函數(shù)，其調(diào)用格式如下：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.調(diào)用時(shí)可缺省xs。 例1.4.1的MATLAB求解程序ep141.m如下：%1.4.1.m:調(diào)用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)a=0.8；ys=1; %設(shè)差分方程系數(shù)a=0.8，初始狀態(tài)：y(-1)=1xn=1,z

13、eros(1,30); %x(n)=單位脈沖序列，長度N=31B=1;A=1,-0.8; %差分方程系數(shù)xi=filtic(B,A,ys); %由初始條件計(jì)算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調(diào)用fiter解差分方程，求系統(tǒng)輸出信號y(n)n=0:length(yn)-1;stem(n,yn,.)title(時(shí)域波形圖);xlabel(n);ylabel(y(n)程序中取查分方程系數(shù)a=0.8時(shí)，得到系統(tǒng)輸出y（n）如圖1.4.1（a）所示，與例1.4.1的解析遞推結(jié)果完全相同。如果令初始條件y（-1）=0（僅修改程序中ys=0），則得到系統(tǒng)輸出y（n）=

14、h（n），如圖1.4.1（b）所示。(a) (b)圖（a）為a=0.8，y（-1）=1時(shí)，系統(tǒng)輸出時(shí)域波形圖，圖（b）為a=0.8，y(-1)=0時(shí)，系統(tǒng)輸出時(shí)域波形圖。 四設(shè)計(jì)代碼及結(jié)果MATLAB源程序源程序如下%1.m:調(diào)用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)ys=1; %初始狀態(tài)：y(-1)=1xn=1,zeros(1,30); %x(n)=單位脈沖序列，長度N=31B=1;A=1,-0.8; %差分方程系數(shù)xi=filtic(B,A,ys); %由初始條件計(jì)算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調(diào)用filter解差分方程，求

15、系統(tǒng)輸出信號y(n)n=0:length(yn)-1; %n的取值范圍stem(n,yn,.) %畫出時(shí)域波形圖title(時(shí)域波形圖);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 %畫出時(shí)域波形圖xlabel(n),ylabel(x(n), title(時(shí)域波形圖) %x軸、y軸分別代表n，x(n)n=-5:5; %n的取值范圍a=0.5; %設(shè)a=0.5xen=a*a.n+a.(-n); %xen=a*a.n+a.(-n)xon=a*a.n-a

16、.(-n); %xon=a*a.n-a.(-n)figure(1); stem(n,xen,filled),grid on %畫出xe(n)時(shí)域波形圖 title(時(shí)域波形圖);xlabel(n);ylabel(xe(n) %x軸、y軸分別代表n，xe(n)figure(2); stem(n,xon,filled),grid on %畫出xo(n)時(shí)域波形圖 title(時(shí)域波形圖);xlabel(n);ylabel(xo(n) %x軸、y軸分別代表n，xo(n)b0=2;b2=0;b2=-1;a1=-0.7;a2=0.1;ys=0; %設(shè)差分方程系數(shù)，初始狀態(tài)：y(-1)=1B=2,0,-1

17、;A=1,-0.7,0.1; %差分方程系數(shù)n=-5:5; %n的取值范圍xn=0.5.n; %x(n)=0.5.nxi=filtic(B,A,ys); %由初始條件計(jì)算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調(diào)用fiter解差分方程，求系統(tǒng)輸出信號y(n)stem(n,yn,.) %畫出時(shí)域波形圖title(a);xlabel(n);ylabel(y(n) %x軸、y軸分別代表n，y(n)五程序運(yùn)行結(jié)果如下： 差分序列時(shí)域波形圖輸入?yún)?shù)a1=-0.7,a2=0.1,b0=2,b1=0,b2=-1 得到二階線性常系數(shù)差分方程為 y(n)-0.7y(n-1)+0.

18、1y(n-2)=2x(n)-x(n-2)結(jié)果如下：x（n）的時(shí)域波形圖 共軛對稱分量xe（n）的時(shí)域波形圖 共軛對稱分量xo（n）的時(shí)域波形圖 輸入初始條件ys=y(-1)=1,得到結(jié)果如下：輸入初始條件ys=y(-1)=1時(shí)，輸出y(n)波形圖改變初始條件ys=y(-1)=0,得到結(jié)果如下：輸入初始條件ys=y(-1)=0時(shí)，輸出y(n)波形圖 改變初始條件ys=y(-1)=100,得到結(jié)果如下: 輸入初始條件ys=y(-1)=100時(shí)，輸出y(n)波形圖 六比較結(jié)果總結(jié)由上實(shí)驗(yàn)可知，通過改變初始條件ys的值，得到的輸出波形大小并不一致，即輸出信號y(n)是不相同的；從而我們可以得出，對于同

19、一個(gè)差分方程和同一個(gè)輸入信號，因?yàn)槌跏紬l件不同，得到的輸出信號是不相同的。七收獲與體會 本次MATLAB課程設(shè)計(jì)讓我熟悉了該軟件的一些功能，但是對于靈活應(yīng)用MATLAB，以及掌握各方面的設(shè)計(jì)思維以及技巧，還需要投入更多的時(shí)間。在熟悉MATLAB程序和操作的同時(shí)培養(yǎng)了我的獨(dú)立思考能力，專研精神，解決問題能力和動手能力。在此之前了解到MATLAB是一個(gè)很重要很有用的工具，但我并沒有完全理解，本課程設(shè)計(jì)中，通過查閱資料，閱讀網(wǎng)上程序并讀寫程序，對于MATLAB的應(yīng)用有了更深的了解，同時(shí)也認(rèn)識到MATLAB功能非常的強(qiáng)大，有著很多方面的應(yīng)用，如繪制函數(shù)，處理音頻，圖像數(shù)據(jù)，創(chuàng)建用戶界面等功能，實(shí)為一個(gè)

20、功能強(qiáng)大的軟件。本次課程設(shè)計(jì)我完成了基于MATLAB的線性常系數(shù)差分方程求解的題目，通過實(shí)際操作回顧所學(xué)的內(nèi)容，強(qiáng)化基礎(chǔ)，實(shí)踐理論知識。相信在以后的學(xué)習(xí)中，還會更加深入的了解MATLAB，應(yīng)用它。隨著課程設(shè)計(jì)報(bào)告的基本完成，本次課程設(shè)計(jì)終于接近了尾聲。本次課程設(shè)計(jì)要求我們利用上學(xué)期所學(xué)的信號與線性系統(tǒng)分析的知識結(jié)合MATLAB編程工具，完成差分方程求解設(shè)計(jì)的題目，通過實(shí)際操作，回顧所學(xué)內(nèi)容，務(wù)實(shí)基礎(chǔ)，強(qiáng)化理論知識，并體驗(yàn)理論與實(shí)際相結(jié)合的過程。設(shè)計(jì)過程中遇到的第一個(gè)問題便是對于MATLAB語言的不熟悉，其實(shí)現(xiàn)在想想這個(gè)問題不應(yīng)該成為問題。畢竟本專業(yè)曾開設(shè)過MATLAB程序設(shè)計(jì)這門課，而且老師還

21、特別提醒過課程設(shè)計(jì)會用到MATLAB語言?？上У氖?，老師的話并沒有引起我的足夠重視，才導(dǎo)致要真正設(shè)計(jì)的時(shí)候一頭霧水，無從下手。通過這件事，我明白了對于一件事情，想要做得很好，提前做功課和準(zhǔn)備是必不可少的，機(jī)會只垂青有準(zhǔn)備的人。當(dāng)然，通過本次課程設(shè)計(jì)，我還是基本熟悉了一些MATLAB模塊以及與本課程有關(guān)的一些函數(shù)的用法。但是由于上學(xué)期信號與線性系統(tǒng)分析課程學(xué)的不是很好，所以也就是在懂得同學(xué)的指導(dǎo)下按部就班的寫了一些代碼，也沒有什么創(chuàng)新，但還是很好的回顧了差分方程求解的內(nèi)容。老師在任務(wù)書中提到的目的要求，我自認(rèn)為多多少少完成了一些。但不可否認(rèn)的是還有很多沒有達(dá)到的。總之，很多東西還是需要自己在不斷摸索中找到答案。同時(shí)，通過本次數(shù)字信號處理課程設(shè)計(jì)，使我全面系統(tǒng)的了解了差分方程求解的原理。通過結(jié)合課本的知識去完成課程設(shè)計(jì)也讓我明白了理論與實(shí)踐的相結(jié)合的重要性。只有