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

1、數(shù)字信號處理課程設(shè)計題目試實現(xiàn)線性常系數(shù)差分方程的求解學(xué)院：專業(yè)： 班級：學(xué)號：組員：指導(dǎo)教師：題目：用Matlab實現(xiàn)線性常系數(shù)差分方程求解一. 設(shè)計要求1.掌握線性常系數(shù)差分方程的求解2 . 熟練掌握Matlab基本操作和各類函數(shù)調(diào)用3 . 結(jié)合Matlab實現(xiàn)線性常系數(shù)差分方程的求解二. 設(shè)計原理1 .差分與差分方程與連續(xù)時間信號的微分及積分運算相對應(yīng)，離散時間信號有差分及序列求和 運算。設(shè)有序列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)可見，前向差分 與后向差分的關(guān)系為f (k) f (k 1)(3.1 3)二者僅移位不同，沒有原則上的差別，因而它們的性質(zhì)也相同。此處主要采 用后向差分，并簡稱其為差分由查分的定義，若有序列fi(k)、f2(k)和常數(shù)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、質(zhì)。二階差分可定義為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!為二項式系數(shù)序列f(k)的求和運算為kf(i) (3.1 8)i差分方程是包含關(guān)于變量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)及其各移位序列均為常數(shù)，就稱其為常系數(shù)差分 方程；如果某些系數(shù)是變量 k 的函數(shù)，就稱其為變系數(shù)差分方程。描述 LTI 離散 系統(tǒng)的是常系數(shù)線性差分方程。差分方程是具有遞推關(guān)系的代數(shù)方程， 若一直初始條件和激勵， 利用迭代法 渴求的差分方程的數(shù)值解。2. 差分方程的經(jīng)典解一般而言，如果但輸入一單輸出的 LTI系統(tǒng)的激勵f(k)，其全響應(yīng)為y(k), 那么，

5、描述該系統(tǒng)激勵f(k)與響應(yīng)y(k)之間關(guān)系的數(shù)學(xué)模型式n階常系數(shù)線性差 分方程，它可寫為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) 都是常數(shù)。上式可縮寫為an j y(kj0j)am i f (ki0i)(式中 an=1)3.110b)與微分方程的經(jīng)典解類似，上述差分方程的解由齊次解和特解兩部分組成齊次解用yh(k)表示，特解用yp(k)表示，即y(k)=yh(k) yp(k)(3.111)a. 齊次解當(dāng)式(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的等比級數(shù)，因此y(k)應(yīng)有如下形式y(tǒng)(k)C( a)k (3.1 14)式中C式常數(shù)，有初始條件確定。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等為待定常數(shù)表3 1不同特征根所對應(yīng)的齊次解特征根齊次解yh(k)單實根C kr重實根(Cr 1kr 1 Cr 2kr 2 L C1k C0) k一對共軛復(fù)根1,2 a jb e JkCcos( k) Dsin( k)或A kcos

8、( k)其中AeJC JDr重共軛復(fù)跟k Av 1kr 1 cos( k r 1) Ar 2kr 2 cos( k r 2) LAo cos( k0)b. 特解特解的函數(shù)形式與激勵的函數(shù)形式有關(guān)，表 3 2列出了集中典型的激勵f(k)所 對應(yīng)的特解卩(燈。選定特解后代入原差分方程，求出其待定系數(shù)P(A或)等, 就得出方程的特解。表3 2不同激勵所對應(yīng)的特解激勵特解yp(k)kmPmkPm lkL pk Po所有特征根均不等于1時PF Pm*"' L Pk Po當(dāng)有r重等于1時的特征根時k a_. kPa當(dāng)a不等于特征根時k(Pk Po)a當(dāng)a是特征單根時Prk Pr1kL P

9、k Poa當(dāng)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)式中各系數(shù)Cj由初始條件確定如果激勵信號是在k=0時接入的，差分方程的解適合于 kX)。對于n

10、階差 分方程，用給定的n個初始條件y(0),y(1),y(n - 1)就可確定全部待定系數(shù)。 如果差分方程的特解都是單根，則方程的全解為式(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)由以上方程可求得全部待定系數(shù)Cj(j 0,1丄，n)2.1零輸入響應(yīng)系統(tǒng)的激勵為零，僅由系統(tǒng)的初始狀態(tài)引起的響應(yīng)，稱為零輸入響應(yīng)，用 yzi(k)表示。在零輸入條件下，式(3.1 10)

11、等號右端為零，化為齊次方程, 即an jYzi(kj 0j)(3.1 25)般設(shè)定激勵是在k=0時接入系統(tǒng)的，在k v 0時，激勵尚未接入，故式(3.125)的幾個初始狀態(tài)滿足yzi( yzi( LL1)2) LLy(y(1)2)3.1 26)yzi(n)y(n)式( 3.126)25)和式( 3.1 26)可求得零輸入響應(yīng) yzi(k)。中的y( 1),y( 2),y( n)為系數(shù)的初始狀態(tài)，由式(3.12.2 零狀態(tài)響應(yīng)當(dāng)系統(tǒng)的初始狀態(tài)為零， 僅由激勵 f(k) 所產(chǎn)生的響應(yīng)， 稱為零狀態(tài)響應(yīng)， 用 表 示。在零狀態(tài)情況下，式( 3.1 10)仍是非齊次方程，其初始狀態(tài)為零，即零 狀態(tài)響應(yīng)

12、滿足nman j yzs(kj)am i f (k i)j 0i 0yzs( 1) yzs( 2) Lyzs( n) 0 (3.1 30)的解。若其特征根均為單根，則其零狀態(tài)響應(yīng)為nyzs(k)Czsj jkj1yp(k)3.1 31 )式中Czsj為待定常數(shù)，yp(k)為特解。yzs( 1),yzs( 2),L ,yzs( n)為零為零,需要指出，零狀態(tài)響應(yīng)的初始狀態(tài)但其初始值 yzs(0), yzs(1),L ,yzs(n 1)不一定等于零。3.線性常系數(shù)差分方程3.1 一個 N階線性常系數(shù)差分方程可用下式表示：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)分別是系統(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)的輸入序列， 通過求解差分方程可以求出輸出序列。 求解差分方程 的基本方法有以下三種：(1 ) 經(jīng)典解法。這種方法類似于模擬系統(tǒng)中求解微

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

15、.4.1 )式，求 n 時刻的輸出，要知道 n 時刻以及 n 時刻以前的輸入-可編輯修改 -序列值，還要知道 n 時刻以前的 N 個輸出序列值。因此求解差分方程在給定輸 入序列的條件下，還需要確定 N 個初始條件。如果求 n0 時刻以后的輸出， n0 時刻以前 N 個輸出值 y(n0-1 )、 y( n0-2 )、 ?y?(? n、0-N )就構(gòu)成了初始 條件。(1.4.1 )式表明，已知輸入序列和 N 個初始條件，則可以求出 n 時刻的輸出； 如果將該公式中的 n 用 n+1 代替，可以求出 n+1 時刻的輸出，因此( 1.4.1 ) 式表示的差分方程本身就是一個適合遞推法求解的方程。三設(shè)計

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

17、態(tài)引起的零輸入響應(yīng)和由輸入信號 xn 引起的零狀態(tài) 響應(yīng)之和。其中， xi 是等效初始條件的輸入序列，所以 xi 是由初始條件確定的。 MATLAB 信號處理工具箱提供的 filtic 就是由初始條件計算 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)用時可缺省 xs。例1.4.1的MATLAB求解程序ep141.m 如下:%1.4.1.m:調(diào)用 MATLAB

18、解差分方程 y(n)-0.8y(n-1)=x(n)a=0.8 ; ys=1;y(-1)=1%設(shè)差分方程系數(shù)a=0.8 ,初始狀態(tài):xn=1,zeros(1,30);%x(n)=單位脈沖序列，長度N=31B=1;A=1,-0.8;%差分方程系數(shù)xi=filtic(B,A,ys); 列xi%由初始條件計算等效初始條件的輸入序yn=filter(B,A,x n,xi);y(n)%調(diào)用fiter解差分方程，求系統(tǒng)輸出信號n=O:le ngth(y n)-1;stem( n,yntitle('時域波形圖)');xlabel('n');ylabel('y(n)

19、9;)程序中取查分方程系數(shù)a=0.8時，得到系統(tǒng)輸出y (n)如圖1.4.1 (a)所 示，與例1.4.1的解析遞推結(jié)果完全相同。如果令初始條件 y (-1)=0 (僅修改 程序中ys=0 )，則得到系統(tǒng)輸出y (n) =h (n)，如圖1.4.1 (b)所示。n-可編輯修改-(b)(a)圖(a)為a=0.8，y (-1 ) =1時，系統(tǒng)輸出時域波形圖，圖(b)為a=0.8， y(-1)=0時，系統(tǒng)輸出時域波形圖。四.設(shè)計代碼及結(jié)果MATLAB 源程序源程序如下%1.m:調(diào)用 MATLAB 解差分方程 y(n)-0.8y(n-1)=x(n)ys=1;% 初始狀態(tài)： y(-1)=1xn=1,ze

20、ros(1,30);%x(n)= 單位脈沖序列，長度 N=31B=1;A=1,-0.8;%差分方程系數(shù)xi=filtic(B,A,ys);% 由初始條件計算等效初始條件的輸入序列 xiyn=filter(B,A,xn,xi);%調(diào)用filter解差分方程，求系統(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.八

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;% 設(shè) 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;%設(shè)差分方程系數(shù)，初始狀態(tài)：y

23、(-1)=1B=2,0,-1;A=1,-0.7,0.1;%差分方程系數(shù)n=-5:5;%n的取值范圍xn=0.5.An;%x(n)=0.5.Anxi=filtic(B,A,ys);% 由初始條件計算等效初始條件的輸入序列 xiyn=filter(B,A,xn,xi);%調(diào)用fiter解差分方程，求系統(tǒng)輸出信號y(n)stem(n,yn,'.')%畫出時域波形圖title('(a)');xlabel('n');ylabel('y(n)') %x軸、y軸分別代表 n , y(n)五程序運行結(jié)果如下：-可編輯修改 -可編輯修改-差分序列時

24、域波形圖得到二階線性常系數(shù)差分輸入?yún)?shù) 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)結(jié)果如下：x (n)的時域波形圖共軛對稱分量xe (n)的時域波形圖共軛對稱分量xo (n)的時域波形圖輸入初始條件ys=y(-1)=1,得到結(jié)果如下:輸入初始條件ys=y(-1)=1 時，輸出y(n)波形圖改變初始條件ys=y(-1)=0,得到結(jié)果如下:輸入初始條件ys=y(-1)=0 時，輸出y(n)波形圖改變初始條件ys=y(-1)=100,得到結(jié)果如下：輸入初始條件ys=y(-1)=100 時，輸出y

25、(n)波形圖六.比較結(jié)果總結(jié)由上實驗可知，通過改變初始條件ys的值，得到的輸出波形大小并不一致, 即輸出信號y(n)是不相同的；從而我們可以得出，對于同一個差分方程和同一個 輸入信號，因為初始條件不同，得到的輸出信號是不相同的。七.收獲與體會本次MATLAB課程設(shè)計讓我熟悉了該軟件的一些功能，但是對于靈活應(yīng)用 MATLAB，以及掌握各方面的設(shè)計思維以及技巧，還需要投入更多的時間。在熟悉MATLAB程序和操作的同時培養(yǎng)了我的獨立思考能力，專研精神，解決冋題 能力和動手能力。在此之前了解到MATLAB是一個很重要很有用的工具，但我并沒有完全理 解，本課程設(shè)計中，通過查閱資料，閱讀網(wǎng)上程序并讀寫程序

26、，對于MATLAB的應(yīng)用有了更深的了解，同時也認識到 MATLAB功能非常的強大，有著很多方 面的應(yīng)用，如繪制函數(shù)，處理音頻，圖像數(shù)據(jù)，創(chuàng)建用戶界面等功能，實為一個 功能強大的軟件。本次課程設(shè)計我完成了基于 MATLAB 的線性常系數(shù)差分方程求解的題目， 通過實際操作回顧所學(xué)的內(nèi)容， 強化基礎(chǔ)，實踐理論知識。 相信在以后的學(xué)習(xí)中， 還會更加深入的了解 MATLAB ，應(yīng)用它。隨著課程設(shè)計報告的基本完成， 本次課程設(shè)計終于接近了尾聲。 本次課程設(shè) 計要求我們利用上學(xué)期所學(xué)的信號與線性系統(tǒng)分析的知識結(jié)合 MATLAB 編程工 具，完成差分方程求解設(shè)計的題目，通過實際操作，回顧所學(xué)內(nèi)容，務(wù)實基礎(chǔ)，

27、強化理論知識，并體驗理論與實際相結(jié)合的過程。設(shè)計過程中遇到的第一個問題便是對于 MATLAB 語言的不熟悉，其實現(xiàn)在 想想這個問題不應(yīng)該成為問題。 畢竟本專業(yè)曾開設(shè)過 MATLAB 程序設(shè)計這門課， 而且老師還特別提醒過課程設(shè)計會用到 MATLAB 語言。可惜的是，老師的話并 沒有引起我的足夠重視， 才導(dǎo)致要真正設(shè)計的時候一頭霧水， 無從下手。 通過這 件事，我明白了對于一件事情， 想要做得很好， 提前做功課和準備是必不可少的， 機會只垂青有準備的人。當(dāng)然，通過本次課程設(shè)計，我還是基本熟悉了一些 MATLAB 模塊以及與本 課程有關(guān)的一些函數(shù)的用法。 但是由于上學(xué)期信號與線性系統(tǒng)分析課程學(xué)的不是 很好，所以也就是在懂得同學(xué)的指導(dǎo)下按部就班的寫了一些代碼， 也沒有什么創(chuàng) 新，但還是很好的回顧了差分方程求解的內(nèi)容。 老師在任務(wù)書中提到的目的要求， 我自認為多多少少完成了一些。 但不可否認的是還有很多沒有達到的。 總之，很 多東西還是需要自己在不斷摸索中找到答案。