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

1、數(shù)字信號處理課程設計題目： 試實現(xiàn)線性常系數(shù)差分方程的求解 學院： 專業(yè)： 班級： 學號： 組員： 指導教師： 題目：用Matlab實現(xiàn)線性常系數(shù)差分方程求解一 設計要求1 掌握線性常系數(shù)差分方程的求解2 熟練掌握Matlab基本操作和各類函數(shù)調用3 結合Matlab實現(xiàn)線性常系數(shù)差分方程的求解二設計原理1差分與差分方程與連續(xù)時間信號的微分及積分運算相對應，離散時間信號有差分及序列求和運算。設有序列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）二者僅移位不同，沒有原則上的差別，因而它們的性質也相同。此處主要采用后向差分，并簡稱其為差分。由查分的定義，若有序列、和常數(shù)，則 （3.14）這表明差分運算具有線性性質。 二階差分可定義為 （3.15）類似的，可定義三階、四階、n階差分。一般地，n階差分 （3.16）式中 （3.17）為二項式系數(shù)序列f(k)的求和運算為 （3.18）差分方程是包含關于變量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ù)線性差分方程。差分方程是具有遞推關系的代數(shù)方程，若一直初始條件和激勵，利用迭代法渴求的差分方程的數(shù)值解。2. 差分方程的經典解一般而言，如果但輸入單輸出的LTI系統(tǒng)的激勵f(k)，其全響應為y(k)，那么，描述該系統(tǒng)激勵f(k)與響應y(k)之間關系的數(shù)學模型式n階常系數(shù)線性差分方程，它可寫為 （3.110a）式中、都是常

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

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

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

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

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

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

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

11、過程1.用MATLAB求解差分方程 MATLAB信號處理工具箱提供的filter函數(shù)實現(xiàn)線性常系數(shù)差分方程的遞推求解，調用格式如下： yn=filter(B,A.xn) 計算系統(tǒng)對輸入信號向量xn的零狀態(tài)響應輸出信號向量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) 計算系統(tǒng)對輸入信號向量xn的全響應輸出信號yn。所謂全響應，就是由初始狀態(tài)引起的零輸入響應和由輸入信號xn引起的零狀態(tài)響應之和。其中，xi

12、是等效初始條件的輸入序列，所以xi是由初始條件確定的。MATLAB信號處理工具箱提供的filtic就是由初始條件計算xi的函數(shù)，其調用格式如下：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解差分方程y(n)-0.8y(n-1)=x(n)a=0.8；ys=1; %設差分方程系數(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); %由初始條件計算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調用fiter解差分方程，求系統(tǒng)輸出信號y(n)n=0:length(yn)-1;stem(n,yn,.)title(時域波形圖);xlabel(n);ylabel(y(n)程序中取查分方程系數(shù)a=0.8時，得到系統(tǒng)輸出y（n）如圖1.4.1（a）所示，與例1.4.1的解析遞推結果完全相同。如果令初始條件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時，系統(tǒng)輸出時域波形圖，圖（b）為a=0.8，y(-1)=0時，系統(tǒng)輸出時域波形圖。 四設計代碼及結果MATLAB源程序源程序如下%1.m:調用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); %由初始條件計算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調用filter解差分方程，求

15、系統(tǒng)輸出信號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軸分別代表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

16、.(-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 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; %設差分方程系數(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); %由初始條件計算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調用fiter解差分方程，求系統(tǒng)輸出信號y(n)stem(n,yn,.) %畫出時域波形圖title(a);xlabel(n);ylabel(y(n) %x軸、y軸分別代表n，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)結果如下：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(n)波形圖 六比較結果總結由上實驗可知，通過改變初始條件ys的值，得到的輸出波形大小并不一致，即輸出信號y(n)是不相同的；從而我們可以得出，對于同

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

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

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