離散LTI系統(tǒng)數(shù)值解法實用程序設計

畢業(yè)設計 （ 論文 ） 題 目 離散 LTI 系統(tǒng)數(shù)值解法實用程序設計 學生姓名 張丹 學號 1010064047 專業(yè)班級 電子信息科學與技術 102 指導教師 龍姝明 完成地點 陜西理工學院 物電學院 505 實驗室 2014 年 06 月 12 日 陜西理工學院畢業(yè)設計論文 離散離散 LTILTI 系統(tǒng)數(shù)值解法實用程序設計系統(tǒng)數(shù)值解法實用程序設計 張丹 （陜理工 物理與電信工程學院 電子信息科學與技術專業(yè) 2010 級 2 班,陜西 漢中,723000） 指導老師：龍姝明 摘要 探索連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為離散 LTI 系統(tǒng)的有效方法及初值條件轉(zhuǎn)化的有效方法，對 LTI 離散系統(tǒng)用迭 代法求解指定時區(qū)上的輸入響應和零狀態(tài)響應，編寫 LTI 連續(xù)系統(tǒng)轉(zhuǎn)化為 LTI 離散系統(tǒng)并用迭代方法求解系統(tǒng)響應 的 Mathematica 程序，以 LTI 連續(xù)系統(tǒng)分析實例展示程序的高效性和可靠性。 關鍵詞 離散系統(tǒng)；迭代法；Mathematica 程序；數(shù)值解 Numerical Solution of discrete LTI system utility designed Zhang Dan (Grade10,Class2,Major Electronic Information Science and Technology School of Physics and Telecommunication Engineering, Shannxi University of Technology,Hanzhong,723000) Tutor:Long Shuming Abstract Explore effective method of continuous LTI system into discrete LTI system effective methods and initial conditions for transformation, for LTI discrete systems with iterative method specified input response and zero status area on the response, the preparation of LTI continuous system into LTI discrete system with iterative method for solving the system response Mathematica program to LTI continuous systems analysis examples demonstrate the efficiency and reliability of the program. Keywords Discrete systems; iterative method;Mathematica program; numerical solution 陜西理工學院畢業(yè)設計論文 目錄 引言 1 1 連續(xù) LTI 系統(tǒng)的解析解和數(shù)值解1 1.1 求解連續(xù)低階 LTI 系統(tǒng)解析解的特點1 1.2 求解連續(xù)高階 LTI 系統(tǒng)解析解存在的問題2 1.3 求解連續(xù) LTI 系統(tǒng)數(shù)值解的特點3 2 連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為離散 LTI 系統(tǒng)的方法4 3 連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為 LTI 離散系統(tǒng)初值條件的方法5 4 迭代法求解指定時區(qū)上的 LTI 連續(xù)系統(tǒng)的編程技巧5 4.1 求解指定時區(qū)上零輸入響應的編程技巧 .5 4.2 求解指定時區(qū)上零狀態(tài)響應的編程技巧 .6 4.3 求解指定時區(qū)上全響應的編程技巧 .6 5 迭代法求解 LTI 系統(tǒng)的 Mathematica 程序設計思路7 6 程序應用實例 8 6.1 用 Mathematica 軟件編程求解連續(xù) LTI 系統(tǒng)數(shù)值解的思 路8 6.2 用 Mathematica 軟件編程求解連續(xù) LTI 系統(tǒng)的實 例8 7 結語 9 參考文獻 9 附錄 A .9 附錄 B .10 附錄 C .10 陜西理工學院畢業(yè)設計論文 第 0 頁 共 15 頁 引言引言 隨著數(shù)字技術以及計算機技術的飛速發(fā)展，鑒于離散系統(tǒng)在精度、抗干擾能力和可集成化等方 面，比連續(xù)系統(tǒng)具有更大的優(yōu)越性，因此原來對連續(xù)信號和系統(tǒng)的研究問題，越來越多地轉(zhuǎn)化為對 離散信號和系統(tǒng)的處理問題。通信和計算機設備等數(shù)字化的高科技產(chǎn)品滲透于人們的生活、學習、 工作等諸多方面，這樣，對于離散系統(tǒng)的分析、研究和改進成為了必不可少的課題1。 離散系統(tǒng)的響應問題是求解及分析離散系統(tǒng)的基礎理論問題，是我們深入分析線性時不變離散 系統(tǒng)的基礎。離散 LTI 系統(tǒng)的求解方法有多種，有時域分析法和變換域分析法，其中時域分析法有 利用差分方程直接求解和利用用迭代法求解。本課題所研究的是用時域分析法中的迭代法求解離散 LTI 系統(tǒng)的數(shù)值解2。 利用迭代的方法分析不借助任何變換而直接求解，直觀且準確。根據(jù)差分方程，用迭代的方法 可以求解出零輸入響應 yzi(k)和零狀態(tài)響應 yzs(k)，也可整體求出系統(tǒng)的全響應。迭代法是利用前 n-1 個時點的的響應值和第 n 個時點的輸入信號響應的值求出第 n 個時點的系統(tǒng)響應值，這種方法是逐 次求解直到求出各點的響應值，方法簡單，概念清楚，對于低階的系統(tǒng)手工操作就可以求解出，但 對于高階系統(tǒng)計算量比較大時，利用計算機運算速度快、適合做重復性操作的特點，用 Mathematica 軟件編程實現(xiàn)這一過程，則更方便快捷。迭代法不僅可求解線性系統(tǒng)在指定區(qū)間上的 值，還可求非線性系統(tǒng)在指定區(qū)間上的值。 作為理論上的研究，此課題雖然簡單，但在實用技術上有廣泛的應用，為進一步深入研究奠定 基礎。例如在通信、計算機、自動化等很多領域都離不開對各類離散系統(tǒng)的分析處理，其中必定涉 及高階系統(tǒng)和非線性系統(tǒng)，對于它們的求解，迭代解法是最有效的分析方法。比如我們在聲音的處 理過程中，就是經(jīng)過濾波器，將聲音信號轉(zhuǎn)換為離散的差分方程來處理的，再如我們在處理圖像時， 也是將其轉(zhuǎn)化為離散的差分方程來求解。在未來的“數(shù)字化”工業(yè)發(fā)展進程中，此課題的研究方法 將有更加廣泛和深入的應用3。 本課題研究的方法是時域分析法中的迭代法，研究的工具軟件是 Mathematica。 Mathematica 是一款科學計算軟件，它有強大的數(shù)值運算和符號運算，并且能與其他應用程序 連接，它的很多功能在科學研究領域中處于世界領先地位，截至 2009 年，它已是使用最廣泛的理 工軟件之一。Mathematica 的發(fā)布標志著現(xiàn)代科技計算的開始，Mathematica 是世界上通用計算系統(tǒng) 中最強大的系統(tǒng)，自從 1988 發(fā)布以來，它已經(jīng)對如何在科技領域運用計算機技術產(chǎn)生了深刻的、 巨大的和廣泛的影響。 Mathematica 和 MATLAB、Maple 統(tǒng)稱為稱為三大數(shù)學軟件。本課題中用 Mathematica 軟件用迭代 的方法進行編程求解系統(tǒng)的數(shù)值解。 1 連續(xù)連續(xù) LTI 系統(tǒng)系統(tǒng)的解析解和數(shù)值解的解析解和數(shù)值解 1.1 求解連續(xù)低階求解連續(xù)低階 LTI 系統(tǒng)解析解的特點系統(tǒng)解析解的特點 當用數(shù)字計算機求解 LTI 連續(xù)系統(tǒng)的解析解時，或直接在系統(tǒng)中采用數(shù)字計算機進行求解時， 對于連續(xù)低階系統(tǒng)，可以通過 Mathematica 軟件編程來實現(xiàn)，例如求解一個連續(xù)低階系統(tǒng)，其微分 方程為 yt+5yt+6yt=-10Cos20t （1.1-1） y0=1，y0=-0.5 用 Mathematica 軟件進行編程求解，程序如下： Cleary; eq=yt+5yt+6yt=10Cos20t, y0=1,y0=-0.5; sol=DSolveeq,yt,t1/Expand; yt_=yt/.sol; Plotyt,t,0,5 陜西理工學院畢業(yè)設計論文 第 1 頁 共 11 頁 所得的結果如圖 1.1 所示： 12345 0.2 0.4 0.6 0.8 1.0 圖 1.1 輸出信號響應圖像 由上面這個例子可知，對于連續(xù) LTI 低階系統(tǒng)，用 Mathematica 軟件編程可以求解出它的解析 解，解的過程很容易而且沒有出現(xiàn)錯誤。 1.2 求解連續(xù)高階求解連續(xù)高階 LTI 系統(tǒng)解析解存在的問題系統(tǒng)解析解存在的問題 求連續(xù)高階 LTI 系統(tǒng)的解析解，就不會像求解連續(xù)低階 LTI 系統(tǒng)那么簡單。例如求解高階微分 方程 yt+150000yt-12yt-5.67yt+123yt=100Sin15.7t （1.2- 1） y0=0，y0=-4，y0=-1.86，y0=12 用 Mathematica 軟件進行編程求解，程序如下： eq=yt+150000yt-12yt-5.67yt+123yt =100Sin15.7t, y0=0,y0=-4,y0=-1.86,y0=12; sol=DSolveeq,yt,t1/Expand; yt_=yt/.sol; Plotyt,t,0,50; 用這種方法求解，在解的過程中出現(xiàn)了問題，如下面所示： RowReduce:luc：病態(tài)矩陣1.+0.i,0.+0.i,1.+0.i,1. +0.i,1.722710-7+1.6311910-21i, (-3.5556410-15-3.9283710-35i)-150000.t- (56.4103-1.5275210-17i)-0.093707 t+ (56.4103-1.5276810-17i)0.0468935tCos0.0809425t- (0.00161326+0.i)Cos15.7t+ (0.00161344_+1.6311910-21i)Cos0.0809425t2 Cos15.7t+(1.1991810-23-2.7377410-41i) Cos0.161885t Cos15.7t- (147.405_-2.6533610-17i) 0.0468935tSin0.0809425 t+ (4.3368110-19+2.8332410-21 i) Cos0.0809425tCos15.7t Sin0.0809425t+ (0.00161344_+0.i)Cos15.7t Sin0.0809425 t2-(9.7419310-20+0.i) Cos15.7tSin0.161885t+ (9.6289210-6+3.6734210-40i)Sin15.7t- (9.628910-6+9.5523410-24i) Cos0.0809425t2 Sin15.7t+(1.0044310-21-2.2931310-39i) 陜西理工學院畢業(yè)設計論文 第 2 頁 共 11 頁 Cos0.161885 tSin15.7 t-(3.3881310-21-1.6555810-23i) Cos0.0809425 tSin0.0809425tSin15.7t- (9.628910-6+1.1020310-39i) Sin0.0809425t2Sin15.7t+ (5.815810-22+4.0783210-56i) Sin0.161885tSin15.7t 解的結果如圖 1.2 所示： 12345 0.2 0.4 0.6 0.8 1.0 圖 1.2 輸出信號響應圖像 通過這個例子，我們發(fā)現(xiàn)要求解此類連續(xù)高階系統(tǒng)的解析解是比較困難的，雖然也可以得到結 果，但是卻耗費較長時間，而且求出的結果表達較為復雜，不僅含有實函數(shù)，而且含有復函數(shù)。當 給出一個系數(shù)很大很復雜的連續(xù)高階 LTI 系統(tǒng)，我們就解不出來了，因為求解一個復雜的連續(xù)高階 LTI 系統(tǒng)，后臺實際上是解一個一元高次代數(shù)方程，在用 Mathematica 解這個系數(shù)很復雜的一元高 次代數(shù)方程時，會有很大的誤差，故在用 Mathematica 進行編程求解此類系統(tǒng)時，就解不出來了。 例如求解微分方程 yt+18000000yt-12yt-5.67yt+123yt=100Sin15.7t （1.2- 2） y0=0,y0=-4，y0=-1.86，y0=12 程序如下： Cleary; eq=yt+18000000yt-12yt-5.67yt+123yt =100Sin15.7t,y00,y0=-4,y0=-1.86,y0=12; sol=DSolveeq,yt,t1/Expand; yt_=yt/.sol Plotyt,t,0,50 對于上面這個微分方程，運行了好長時間，仍然沒有求出結果。故對于此類連續(xù)高階系統(tǒng)，就 需要將它化為離散系統(tǒng)，利用 Mathematica 來求解它在指定區(qū)間上的數(shù)值解。這樣既簡便又可靠， 而且不會出現(xiàn)的錯誤。 1.3 求解連續(xù)求解連續(xù) LTI 系統(tǒng)數(shù)值解的特點系統(tǒng)數(shù)值解的特點 對于連續(xù) LTI 系統(tǒng)的微分方程 yt+150000yt-12yt-5.67yt+123yt=100Sin15.7t （1.3- 1） y0=0，y0=-4，y0=-1.86，y0=12 由于求它的解析解既耗時結果又復雜，故我們求它的數(shù)值解，用 Mathematica 進行編程，求解它在 一個 0-50s 這個時間段上的數(shù)值解，程序如下: Cleary; eq=yt+150000yt-12yt-5.67yt+123yt= 100Sin15.7t,y0=0,y0=-4,y0=-1.86,y0=12; 陜西理工學院畢業(yè)設計論文 第 3 頁 共 11 頁 sol=NDSolveeq,yt,t,0,50,MaxSteps10000001/Expand; yt_=yt/.sol Plotyt,t,0,50 所得到的結果如圖 1.3 所示： InterpolatingFunction0.,50.,t 1020304050 600 400 200 200 400 600 800 圖 1.3 輸出信號響應圖像 通過上面的例子，我們知道用 Mathematica 進行編程求解此類方程在指定區(qū)間上的數(shù)值解，非 常簡便，而且很實用。 2 連續(xù)連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為離散系統(tǒng)轉(zhuǎn)化為離散 LTI 系統(tǒng)的方法系統(tǒng)的方法 由于用離散系統(tǒng)求解數(shù)值解更為簡便、快捷和準確，故通常我們先將連續(xù)系統(tǒng)轉(zhuǎn)化為離散系統(tǒng)， 離散化就是導出能在采樣時刻上與連續(xù)系統(tǒng)狀態(tài)等價的離散化狀態(tài)方程。連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為離散 LTI 系統(tǒng)，包含自變量的離散化、導函數(shù)的離散化和方程的離散化。 例如對于連續(xù) LTI 方程 yt+ayt+byt=cft （2- 1）y0_=y0，y0_=yp0 要將它離散化為離散 LTI 系統(tǒng)，首先是確定求解的時間區(qū)間，將自變量也就是時間離散化，用 Mathematica 進行編程，程序如下： ta=0;tb=10; n=100;Ts=(tb-ta)/n; ts=Range0,n*Ts; 其中，ta 和 tb 是初始時刻和終止時刻，為采樣數(shù)，Ts 為采樣間隔，ts 為采樣時間。然后是輸 入函數(shù)的離散化，將 f(t)f(ts(k)，f(ts(k)f(k)，f(t)f(k)。再是導函數(shù)的離散化， 程序如下： y(k)(y(k)-y(k-1)/Ts; y(k)(y(k)-2y(k-1)+y(k-2)/Ts2; y(k)(y(k)-3y(k-1)+3y(k-2)-y(k-3)/Ts3; y(k)(y(k)-4y(k-1)+6y(k-2)-4y(k-3)+y(k-4)/Ts4; .; 其中 y(k )是一階導函數(shù)的離散化，y(k)是二階導函數(shù)的離散化，y(k)是三階導函 數(shù)的離散化，y(k) 是四階導函數(shù)的離散化，依次類推，將 y(m)(t)轉(zhuǎn)化為 y(m)(k)。最后 是方程的離散化，對于對于方程 yt+ayt+ byt=cft，將導函數(shù)離散化的結果代入 此方程，整理后，將各項的系數(shù)分別用字母表示便可得到離散化的方程，用 Mathematica 進行編程， 程序如下： y(k)-2y(k-1)+y(k-2)+aTs(y(k)-y(k-1)+bTs2yk=cTs2fk; (1+a.Ts+b.Ts2)yk+(-2-a.Ts)yk-1+yk-2=cTs2fk; q=1/(1+a.Ts+b.Ts2); yk+(-2-a Ts)qyk-1+qyk-2=cTs2qfk; 陜西理工學院畢業(yè)設計論文 第 4 頁 共 11 頁 p=(-2-aTs)q; =cTs2q; yk+pyk-1+qyk-2=fk; 3 連續(xù)連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為系統(tǒng)轉(zhuǎn)化為 LTI 離散系統(tǒng)初值條件的方法離散系統(tǒng)初值條件的方法 探索連續(xù) LTI 系統(tǒng)初值條件轉(zhuǎn)化為 LTI 離散系統(tǒng)初值條件的有效方法，首先我們是對上面導函 數(shù)離散化的方程進行變形，然后賦值就可得到，具體做法如下： 對y(k)(y(k)-y(k-1)/Ts，進行變形，得到 Tsy(k)y(k)-y(k-1)，然后給 k 賦值為-，將 y(-1)=y0_帶入，再移項就可得到 y(-2)=y(-1)-Tsy0_;對 y(k) (y(k)-2y(k-1)+y(k-2)/Ts2進行變形，得到 Ts2y(k)y(k)-2y(k-1)+y(k-2)， 然后給 k 賦值為-，將 y(-1)=y0_帶入，再移項就可得到 y(-)=y(-)y(-)+ Ts2y0_，用 Mathematica 進行編程，程序如下： y(-1)=y0_; y(-2)=y(-1)-Tsy0_; y(-3)=2y(-2)y(-1)+Ts2y0_; ; 4 迭代法求解指定時區(qū)上的迭代法求解指定時區(qū)上的 LTI 連續(xù)系統(tǒng)的編程技巧連續(xù)系統(tǒng)的編程技巧 用迭代解法可以分別求出零輸入響應和零狀態(tài)響應，也只有用迭代解法可以一次性求出系統(tǒng)的 全響應，其它的方法比如狀態(tài)空間和傳遞函數(shù)都只能分別求系統(tǒng)的完全響應。 迭代解法不僅可以求線性系統(tǒng)的完全響應，還可以求非線性系統(tǒng)模型的完全響應，至于求出的 結果是否有用，取決于非線性系統(tǒng)的模型能不能很好的反應實際系統(tǒng)。而非線性系統(tǒng)的模型能不能 很好的反應實際系統(tǒng)，則取決于構建的描述系統(tǒng)的差分方程是否正確，如果正確，那解出的結果就 能很好的反應該實際系統(tǒng)。 4.1 求解指定時區(qū)上零輸入響應的編程技巧求解指定時區(qū)上零輸入響應的編程技巧 零輸入響應是在沒有外加激勵時，僅有 0 時刻的非零初始狀態(tài)引起的響應。它是取決于初始狀 態(tài)和電路特性，這種響應是隨時間按指數(shù)規(guī)律衰減變化的4。 在用 Mathematica 進行編程的過程中，連續(xù) LTI 系統(tǒng)轉(zhuǎn)化為離散 LTI 系統(tǒng)的方法和連續(xù) LTI 系 統(tǒng)初值條件轉(zhuǎn)化為 LTI 離散系統(tǒng)初值條件的方法上面我們提到，用迭代法求解指定時區(qū)上的 LTI 離 散系統(tǒng)零輸入響應,那我們只需要將輸入信號設為 0，進而求出系統(tǒng)的響應即可。 例如對于方程 yt+ayt+byt=cft （4.1- 1） y(-)=m，y(-2)=n 將上文我們研究的導函數(shù)離散化的結果代入（4.1-1） ，再整理便可得到離散化的方程： y(k)-2y(k-1)+y(k-2)+aTs(y(k)-y(k-1)+bTs2yk=cTs2fk （4.1- 2）要用迭代法求解指定時區(qū)上的 LTI 離散系統(tǒng)零輸入響應,那我們可以將上式中的 fk取值為 0，同 類的合并后，把 y(k)保留在等式的左邊，其余項移到等式的左邊，即是 (1+a.Ts+b.Ts2)yk+(-2-a.Ts)yk-1+yk-2=0 （4.1- 3） 然后我們可以依次給 k 賦值為 0、1、2就可以依次得到 y(0)、y(1)、y(2)的值，即是該系統(tǒng) 函數(shù)在指定區(qū)間上的零輸入響應。 例如對于方程 yt+260300yt-12yt-6yt+123yt120Sin17t （4.1- 4） 陜西理工學院畢業(yè)設計論文 第 5 頁 共 11 頁 y00，y0-4，y0-1.9，y012 將它離散化后，用迭代法求解它在 0-50s 的零輸入響應。即就是輸入信號為的情況下，通過已 知的前四個時刻的值依次求出下一個時刻的值，直到求出指定區(qū)間上的所有的時刻的零輸入響應。 將此微分方程中輸入信號取為 0。則微分方程變?yōu)?yt+260300yt-12yt-6yt+123yt （4.1-5）再用 Mathematica 軟件進行編程求解它的零輸入 響應，程序見附錄 A。 結果如圖 4.1 所示： 1020304050 800 600 400 200 0 圖 4.1 零輸入響應的圖像 4.2 求解指定時區(qū)上零狀態(tài)響應的編程技巧求解指定時區(qū)上零狀態(tài)響應的編程技巧 零狀態(tài)響應就是電路的儲能元器件（如電容、電感類元器件）無初始儲能，僅由外部激勵作用 而產(chǎn)生的響應5。 在用 Mathematica 進行編程的過程中，用迭代法求解指定時區(qū)上的 LTI 離散系統(tǒng)零狀態(tài)響應, 那我們只需要將初始狀態(tài)設為 0，進而求出系統(tǒng)的響應即可。 例如對于方程 yt+ayt+byt=cft （4.2- 1）y(-1)=m，y(-2)=n 將上文研究的導函數(shù)離散化的結果代入（4.6） ，再整理便可得到離散化的方程： y(k)-2y(k-1)+y(k-2)+aTs(y(k)-y(k-1)+bTs2yk=cTs2fk （4.2- 2） 合并，把 y(k)保留在等式的左邊，其余項移到等式的左邊，即是 (1+a.Ts+b.Ts2)yk+(-2-a.Ts)yk-1+yk-2=cTs2fk （4.2- 3） 用迭代法求解指定時區(qū)上的 LTI 離散系統(tǒng)零狀態(tài)響應,那我們可以將上式中的 y(-)和 y(-2)取 值為 0，然后我們可以依次給 k 賦值為 0、1、2就可以依次得到 y(0)、y(1)、y(2) 的值， 即是該系統(tǒng)函數(shù)在指定區(qū)間上的零狀態(tài)響應。 同樣對于方程 yt+260300yt-12yt-6yt+123yt120Sin17t （4.2- 4） y00，y0-4，y0-1.，y012 用迭代法求解它在 0-100s 的零狀態(tài)響應。即就是初始狀態(tài)都為的情況下通過已知的前四個 時刻的值依次求出下一個時刻的值，直到求出指定區(qū)間上的所有的時刻的零狀態(tài)響應。則 y0 0,y0，y0，y0，用 Mathematica 軟件進行編程求解，程序見附 B。 結果如圖 4.2 所示： 陜西理工學院畢業(yè)設計論文 第 6 頁 共 11 頁 20406080100 0.20 0.15 0.10 0.05 0.05 圖 4.2 零狀態(tài)響應的圖像 4.3 求解指定時區(qū)上零狀態(tài)響應的編程技巧求解指定時區(qū)上零狀態(tài)響應的編程技巧 全響應就是線性系統(tǒng)或電路在激勵作用下產(chǎn)生的零狀態(tài)響應與零輸入響應之和。它是系統(tǒng)或電 路在輸入和初始條件共同作用下的響應。是零輸入響應和零狀態(tài)響應疊加的結果,也體現(xiàn)了線性電 路的疊加性6。 用迭代法求解指定時區(qū)上零狀態(tài)響應的，即是綜合了零輸入零輸入響應和零狀態(tài)響應，輸入狀 態(tài)和初始狀態(tài)均不為 0，其他的編程技巧均不變。 5 迭代法求解迭代法求解 LTI 系統(tǒng)的系統(tǒng)的 Mathematica 程序程序設計思路設計思路 例如用迭代解法求解方程 yt+150000yt-12yt-5.67yt+123yt100Sin15.7t （5.1- 1） y00,y0-4，y0-1.86，y012 的數(shù)值解，并用 Mathematica 進行編程求解它在 0-50s 的解。 在離散化的過程中，將導函數(shù)的離散化結果帶入到微分方程（5.1）中，合并同類項，將 yk 的系數(shù)化為 1，其它項的系數(shù)分別用字母表示，這樣就得到離散化的差分方程為： yk+yk-1+yk-2+yk-3+qyk-4=qTs4fnk （5.1- 2） 求解的具體過程，用 Mathematica 編程，程序如下： Cleary,u; ta=0;tb=50;n=200; Ts=(tb-ta)/n; 自變量的離散化編程： ts=Range0,n*Ts; 輸入函數(shù)的離散化編程： fn=100*Sin15.7ts; 方程的離散化編程： q=1/(1+150000*Ts-12*Ts2-5.67*Ts3+123Ts4); =q*(-4-3*150000*Ts+2*12*Ts2+5.67*Ts3); =q*(6+3*150000*Ts-12*Ts2); =q*(-4-150000*Ts); 初值的離散化編程： y-1=0; y-2=4Ts; y-3=-1.86Ts2+8Ts; y-4=-3*1.86Ts2+12Ts-12Ts3; 將數(shù)值解存入數(shù)組編程： u=ConstantArray0,n+1+4; 陜西理工學院畢業(yè)設計論文 第 7 頁 共 11 頁 u1;4=y-4,y-3,y-2,y-1; Forj=5,jn+5,j+,uj=-A*uj-1-B*uj-2-C*uj-3- q*uj-4+q*Ts4 fnj-4; 用圖像表示解的結果： data=Transposets,u5;-1; ListPlotdata,JoinedTrue 畫出它的圖形，如圖.所示： 1020304050 600 400 200 200 400 600 800 圖 5.1 輸出信號響應波形圖 6 程序應用實例程序應用實例 6.1 用用 Mathematica 軟件編程求解連續(xù)軟件編程求解連續(xù) LTI 系統(tǒng)數(shù)值解的思路系統(tǒng)數(shù)值解的思路 （1）先給定一個連續(xù)系統(tǒng)，比如一個三階的 RLC 電路，對電路進行分析，得到它的系統(tǒng)函數(shù)； （2）根據(jù)系統(tǒng)函數(shù)寫出它的微分方程，再將微分方程離散化為差分方程，連續(xù)系統(tǒng)的初值也 離散化為差分方程的初值； （3）利用 Mathematica 軟件和迭代法求出差分方程的數(shù)值解。 6.2 用用 Mathematica 軟件編程求解連續(xù)軟件編程求解連續(xù) LTI 系統(tǒng)的實例系統(tǒng)的實例 以 RLC 三階電路電路為例7，具有電阻電感電容的二端網(wǎng)絡如圖 5.1 所示，其中：R1=2 ，R2=80，L1=0.4H, L2=0.2H ,C=5*10-3F。電壓為輸入，電壓=150Cos18t為輸 1( ) u t 2( ) u t 出，求該三階電路系統(tǒng)的離散化數(shù)值解。 圖 6.1 RLC 時域電路 利用拉普拉斯變換進行分析，建立復頻域代數(shù)方程8： L2 和 R2 的串聯(lián)阻抗為 （6.2-122zL sR 1） C 和 L2 與 R2 的并聯(lián)阻抗為 （6.2- 1 2 11z z z sc 2）總阻抗為 （6.2-3211zzRL s 3） 系統(tǒng)函數(shù)為 陜西理工學院畢業(yè)設計論文 第 8 頁 共 11 頁 23 2 12121 2 (s 2 11 21 ) 2 R RRLLcR Rsc L H RL RscL L s (6.2-4) 代入 R,L,C 的值，得系統(tǒng)函數(shù)9： （6.2- 23 200000 2050003500 ( 4 s) 05 H sss 5） 根據(jù)系統(tǒng)函數(shù)可以列出微分方程： yt+405yt+3500yt+205000yt=200000ft y00，y00，y00 （6.2- 6） 用上文的方法將（6.2-6）離散化為以下形式： yk+yk-1+yk-2+yk-3=200000qTs3fnk （6.2- 7） y-1=0，y-2=0，y-3=0 用 Mathematica 進行編程求解，程序見附錄 C10。 結果如圖.所示： 12345 400 200 200 400 圖 6.3 輸出系統(tǒng)響應波形圖 7 結語結語 畢業(yè)設計是對大學四年學習成果的一次大檢閱，平時課堂上學到的知識很難以融會貫通，通過 本次畢業(yè)設計，讓我在平時學習的知識得到了進一步鞏固和加強，通過畢業(yè)設計還可以將平時所學 的一些知識應用到實際的設計中。 設計剛開始時，由于對編程不是很熟悉，出現(xiàn)了許多錯誤，造成了多次的返工。但是，正是這 一次次的嘗試磨練了我的耐性，并提高了我對軟件的操作水平。在這次設計中，我不僅收獲了專業(yè) 知識，還在與同學的溝通交流中增長了很多的見識，特別是要非常感謝龍老師細心和認真的指導， 正因為有了老師和同學的幫助，我的畢業(yè)設計才能順利完成。這次畢業(yè)設計為我未來踏上社會、步 入工作崗位打下了良好基礎。 這個畢業(yè)設計，讓我深深地體會到這是一個連接學習和工作的橋梁。畢業(yè)設計的完成標志著大 學生活的結束，今后迎接我的是更多的挑戰(zhàn)，通過畢業(yè)設計的磨練，我相信我能夠更好的面對這些， 把握機遇，創(chuàng)造未來。在大學里我得到了最好的鍛煉，我要將學到的知識轉(zhuǎn)換成力量，為了自己的 夢想而努力奮斗。 陜西理工學院畢業(yè)設計論文 第 9 頁 共 11 頁 參考文獻參考文獻 1錢琳琳,牛瑞燕,李秀麗.離散 LTI 系統(tǒng)單位脈沖響應求解方法研究N.電氣電子學報,2010(01). 2 Alan VOppenheim.Discrete-Time Signal Processing（Third Edition）M北京:電子工業(yè)出版社,2011:9-70. 3張正文,鐘東.基于 Matlab 的離散時間系統(tǒng)分析N.咸寧學院學報,2007(06). 4楊忠根,任蕾.陳紅亮因果周期信號通過 LTI 系統(tǒng)的零狀態(tài)響應N.電氣電子教學學報,2011(03). 5柴黎,姚秀芳.基于復頻域的連續(xù)時間 LTI 系統(tǒng)時域分析J,北京:電子技術,2013(04):1-3. 6侯靜怡,劉志杰.離散時間 LTI 系統(tǒng)的響應求解與其 Matlab 的實現(xiàn)J,北京:中國科技信息 2012(2