[導(dǎo)讀]簡介：根據(jù)學(xué)生提出的如何使用C語言來實(shí)現(xiàn)系統(tǒng)函數(shù)的問題，給出了利用MATLAB，Python對(duì)連續(xù)時(shí)間系統(tǒng)函數(shù)進(jìn)行離散化的方法。繼而可以進(jìn)行C語言實(shí)現(xiàn)。關(guān)鍵詞：系統(tǒng)函數(shù)，離散化，MATLAB01如何實(shí)現(xiàn)傳遞函數(shù)？??卓晴老師，請(qǐng)教一下，傳遞函數(shù)環(huán)節(jié)用C需要編寫，該怎么編寫？在ma...

簡介：根據(jù)學(xué)生提出的如何使用C語言來實(shí)現(xiàn)系統(tǒng)函數(shù)的問題，給出了利用MATLAB，Python對(duì)連續(xù)時(shí)間系統(tǒng)函數(shù)進(jìn)行離散化的方法。繼而可以進(jìn)行C語言實(shí)現(xiàn)。
關(guān)鍵詞： 系統(tǒng)函數(shù)，離散化，MATLAB

01 如何實(shí)現(xiàn)傳遞函數(shù)？

??卓晴老師，請(qǐng)教一下，傳遞函數(shù)環(huán)節(jié)用C需要編寫，該怎么編寫？在matlab中他能簡單的表述，但是在利用C語言，該咋辦呢？

??例如寫進(jìn)單片機(jī)中。類似于寫一個(gè)傳遞函數(shù)通式，想實(shí)現(xiàn)各種環(huán)節(jié)。

??在matlab中他有很好的封裝，卻不好理解matlab怎么做到的，simulink中一個(gè)時(shí)域信號(hào)輸入傳遞模塊輸出還是時(shí)域信號(hào)，逐步仿真我認(rèn)為這個(gè)傳遞模塊中封裝著一個(gè)復(fù)雜表達(dá)式，想知道用C語言咋能實(shí)現(xiàn)，這樣以后對(duì)信號(hào)處理想加個(gè)濾波環(huán)節(jié)，一階延時(shí)，可以利用這種方法直接解決。

??好的，感謝卓大大，最好給出一個(gè)G/TS 1這種的一個(gè)C代碼的例子，然后我就可以仿照思路寫出更多高階的。

02 什么是傳遞函數(shù)？

??要回答前面同學(xué)提出的“如何使用C編程語言實(shí)現(xiàn)傳遞函數(shù)？”，則需要先了解什么是系統(tǒng)傳遞函數(shù)？為什么引入系統(tǒng)傳遞函數(shù)？

一、線性時(shí)不變系統(tǒng)

??下面兩個(gè)由線性R，C組成的電路網(wǎng)絡(luò)，它們的共同之處都屬于 線性時(shí)不變系統(tǒng)^[1] ，不同之處在于左邊的電路輸入輸出之間是一個(gè)比例關(guān)系，

??這類系統(tǒng)被稱為即時(shí)系統(tǒng)。

▲ 圖2.1.1 左：電阻分壓電路；右：阻容低通濾波電路

??右邊的低通濾波電路輸入輸出之間關(guān)系則需要通過微分方程才能夠描述：

??這種帶有儲(chǔ)能器件的電路，當(dāng)前的輸出不僅與當(dāng)前的輸入有關(guān)系，還和之前時(shí)間輸入有關(guān)系。這類系統(tǒng)被稱為動(dòng)態(tài)系統(tǒng)，顯然動(dòng)態(tài)系統(tǒng)比即使系統(tǒng)復(fù)雜，但功能更強(qiáng)大。即使系統(tǒng)可以看成一類特殊的動(dòng)態(tài)系統(tǒng)。

二、系統(tǒng)函數(shù)

??既然兩個(gè)電路都屬于線性時(shí)不變系統(tǒng)，可否將它們的輸入輸出關(guān)系都簡化成比例關(guān)系，這樣分析與求解就比較方便了？

??答案是可以的，只不過需要把輸入輸出信號(hào)都進(jìn)行某種變換，比如 傅里葉變換^[2] ， 拉普拉斯變換^[3] ，這時(shí)變換后的系統(tǒng)輸入輸出之間就是一個(gè)比例關(guān)系了。這個(gè)比值就被稱為線性時(shí)不變分析對(duì)象的系統(tǒng)函數(shù)。

??在拉普拉斯變換下，上面的電阻分壓電路對(duì)應(yīng)的系統(tǒng)函數(shù)為：

??阻容低通濾波電路對(duì)應(yīng)的系統(tǒng)函數(shù)為：

??如果拋開這些變換概念背后的數(shù)學(xué)基礎(chǔ)，從工程應(yīng)用角度也可以將上述系統(tǒng)函數(shù)看成對(duì)描述線性時(shí)不變系統(tǒng)輸入輸出關(guān)系之間的微分方程的另外一種表示方式。

??借助于系統(tǒng)函數(shù)的概念，就把線性時(shí)不變系統(tǒng)的分析、設(shè)計(jì)都可以進(jìn)行簡化。

▲ 圖2.2.1 在變換域內(nèi)定義的系統(tǒng)函數(shù)等于系統(tǒng)的零狀態(tài)響應(yīng)與輸入信號(hào)的比值

三、離散時(shí)間系統(tǒng)函數(shù)

??計(jì)算機(jī)的普及使得離散時(shí)間系統(tǒng)應(yīng)用增加，系統(tǒng)可以通過軟件編程來實(shí)現(xiàn)。如果實(shí)現(xiàn)對(duì)于輸入信號(hào)進(jìn)行低通濾波，也可以在對(duì)信號(hào)進(jìn)行采樣之后利用下面離散時(shí)間系統(tǒng)來進(jìn)行處理。

▲ 圖2.3.1.1 一階低通離散時(shí)間系統(tǒng)

??下圖中顯示的隨機(jī)噪聲的原始信號(hào)，經(jīng)過上面迭代方程濾波之后，就會(huì)形成平滑后的數(shù)據(jù)。其中的濾波器系數(shù)。

▲ 圖2.3.3 增加有噪聲信號(hào)以及一階低通濾波后的信號(hào)

from?headm?import?*

t?=?linspace(0,?10,?1000)
sint?=?sin(t)

data?=?sint? ?random.rand(len(t))

def?filter_1(data,?a):
????yn?=?data[0]
????newd?=?[yn]

????for?d?in?data[1:]:
????????yn?=?d*(1-a)? ?a?*?yn
????????newd.append(yn)

????return?newd

data_filter?=?filter_1(data,?0.99)

plt.plot(t,?data,?label='origin')
plt.plot(t,?data_filter,?label='filter?data')
plt.xlabel("t")
plt.ylabel("data")
plt.grid(True)
plt.legend(loc='upper?right')
plt.tight_layout()
plt.show()

??線性時(shí)不變動(dòng)態(tài)離散時(shí)間系統(tǒng)可以使用差分方程來刻畫，經(jīng)過簡單的一項(xiàng)便可成寫出的迭代方程形式，即由輸入信號(hào)以及以前的輸出信號(hào)計(jì)算出當(dāng)前系統(tǒng)輸出值，這是利用軟件編程實(shí)現(xiàn)的基礎(chǔ)。

??同樣，在對(duì)于系統(tǒng)的輸入輸出信號(hào)進(jìn)行 Z-變換^[4] ，那么系統(tǒng)的零狀態(tài)輸出信號(hào)比上輸入信號(hào)也是一個(gè)常量，定義為離散時(shí)間系統(tǒng)的系統(tǒng)函數(shù)。上述一階低通離散時(shí)間系統(tǒng)對(duì)應(yīng)的系統(tǒng)函數(shù)為：

??同樣，如果拋開z變換背后的數(shù)學(xué)原理，僅僅把看成對(duì)于信號(hào)的延遲操作，那么上面的系統(tǒng)函數(shù)也可以看成差分方程另外一種描述方法?？梢钥吹较到y(tǒng)函數(shù)所對(duì)應(yīng)的有理分式的分子多項(xiàng)式系數(shù)對(duì)應(yīng)著差分方程輸入信號(hào)前的系數(shù)；分母多項(xiàng)式的系數(shù)對(duì)應(yīng)差分方程輸出信號(hào)前的系數(shù)。

▲ 圖2.3.2 離散時(shí)間系統(tǒng)的系統(tǒng)方程

??因此，如果已知知道了離散時(shí)間系統(tǒng)函數(shù)，則按照上面的系統(tǒng)函數(shù)與差分方程對(duì)應(yīng)關(guān)系可以寫出對(duì)應(yīng)的差分方程，進(jìn)而整理成對(duì)應(yīng)的迭代方程。這樣就可以使用編程語言（C,Python）來實(shí)現(xiàn)算法了。

03 如何將H(s)轉(zhuǎn)換成H(z)?

??本文一開始同學(xué)詢問，如果知道了一個(gè)連續(xù)時(shí)間系統(tǒng)的傳遞函數(shù)，怎么使用C語言編程實(shí)現(xiàn)系統(tǒng)功能？這個(gè)問題就轉(zhuǎn)換成如何將連續(xù)時(shí)間系統(tǒng)函數(shù)轉(zhuǎn)換成對(duì)應(yīng)的離散時(shí)間系統(tǒng)函數(shù)的問題了。這個(gè)過程也稱為連續(xù)時(shí)間系統(tǒng)的離散化。

一、連續(xù)系統(tǒng)函數(shù)離散化

1、離散化方法

??將連續(xù)時(shí)間系統(tǒng)轉(zhuǎn)換成對(duì)應(yīng)的離散時(shí)間系統(tǒng)的過程常常用在模擬設(shè)計(jì)-離散實(shí)現(xiàn)的設(shè)計(jì)過程中。通常使用的方法包括：

沖激響應(yīng)不變法：這種方法可以保持離散時(shí)間系統(tǒng)的單位沖激響應(yīng)是對(duì)應(yīng)連續(xù)時(shí)間系統(tǒng)的采樣；
變量替換法：包括前向差分（Euler方法）、后向差分、雙線性法，也稱梯形方法；

??除此之外，在 Continuous^[5] 中給出了在控制系統(tǒng)中常用到的六種方法及其轉(zhuǎn)換原理。

2、舉例

??比如使用雙線性方法，它是基于把代表位移算子的所對(duì)應(yīng)的連續(xù)系統(tǒng)函數(shù)近似展開成有理分式：

??以一個(gè)一階低通濾波器為例：

??使用雙線性方法轉(zhuǎn)換成離散系統(tǒng)函數(shù)：

^{二、使用MATLAB離散化系統(tǒng)函數(shù)
??利用MATLAB內(nèi)部函數(shù) ?c2d^[6] ，可以很方便將連續(xù)時(shí)間系統(tǒng)函數(shù)轉(zhuǎn)換成對(duì)應(yīng)的離散時(shí)間系統(tǒng)函數(shù)。
1、命令語法
??MATLAB 中給出了一下四種常用到離散系統(tǒng)函數(shù)命令形式：
sysd?=?c2d(sysc,Ts)
sysd?=?c2d(sysc,Ts,method)
sysd?=?c2d(sysc,Ts,opts)
[sysd,G]?=?c2d(___)
??第一種方法利用了零階保持方法完成系統(tǒng)函數(shù)離散化，如果希望指明離散化具體方法，可以使用2,3兩種命令形式。
1、轉(zhuǎn)換舉例
（1）連續(xù)時(shí)間系統(tǒng)函數(shù)
??將如下帶有延時(shí)環(huán)節(jié)的二階連續(xù)時(shí)間系統(tǒng)進(jìn)行離散化：
$$H\left( s \right) = e^{{ - 0.25s} {{10} \over {s}2 ? 3s 10}} $$（2）轉(zhuǎn)換命令
h?=?tf(10,[1?3?10],'IODelay',0.25);?
hd?=?c2d(h,0.1)
hd?=
?
???????????0.01187?z^2? ?0.06408?z? ?0.009721
??z^(-3)?*?----------------------------------
?????????????????z^2?-?1.655?z? ?0.7408
?
Sample?time:?0.1?seconds
Discrete-time?transfer?function.
（3）系統(tǒng)仿真
??通過以下指令，可以對(duì)比連續(xù)時(shí)間系統(tǒng)與離散時(shí)間系統(tǒng)單位沖激響應(yīng)波形。
step(h,'--',hd','-')
▲ 圖3.2.2.1 系統(tǒng)的單位沖激響應(yīng)信號(hào)對(duì)比
※ 問題總結(jié) ※
??使用嵌入式系統(tǒng)實(shí)現(xiàn)信號(hào)處理算法、控制算法，比起利用連續(xù)時(shí)間系統(tǒng)（或者模擬電路）實(shí)現(xiàn)具有很大的優(yōu)點(diǎn)。在一些模擬設(shè)計(jì)（仿真）-數(shù)字化實(shí)現(xiàn)的場景中，將連續(xù)時(shí)間系統(tǒng)函數(shù)轉(zhuǎn)換成對(duì)應(yīng)的離散時(shí)間系統(tǒng)函數(shù)是編程實(shí)現(xiàn)控制算法的第一步。
??離散化系統(tǒng)函數(shù)的方法有很多種，在MATLAB中可以利用相關(guān)的命令完成這些轉(zhuǎn)換，并進(jìn)行系統(tǒng)仿真測試。
??在 control.TransferFunction^[7] 中給出了利用 Python 中的 control 開發(fā)包進(jìn)行連續(xù)時(shí)間系統(tǒng)離散化相關(guān)指令。
??下面示例給出了給出了將系統(tǒng)函數(shù)：轉(zhuǎn)換成對(duì)應(yīng)的離散時(shí)間系統(tǒng)函數(shù)：

import?control
sys?=?control.TransferFunction(1,[1,1])
sysd?=?sys.sample(0.5,?method='bilinear')
參考資料
[1]線性時(shí)不變系統(tǒng): https://baike.baidu.com/item/線性時(shí)不變系統(tǒng)/2474061[2]傅里葉變換: https://mathworld.wolfram.com/FourierTransform.html#:_{:text=The smoother a function (i.e., the larger the,Fourier transform is also symmetric since implies .}[3]拉普拉斯變換: https://byjus.com/maths/laplace-transform/#::text=The Laplace transform is a well established mathematical,to calculate the solution to the given problem.[4]Z-變換: https://faculty.ksu.edu.sa/sites/default/files/dsp_cen352_ch5_z.pdf#:~:text=Definition of z-transform: The z-transform is a very,filter design and frequency analysis of digital signals.[5]Continuous: https://ww2.mathworks.cn/help/control/ug/continuous-discrete-conversion-methods.html[6]c2d: https://ww2.mathworks.cn/help/control/ref/lti.c2d.html#mw_53fc4689-2099-41d0-93b3-de1e51a174c1[7]control.TransferFunction: https://python-control.readthedocs.io/en/0.8.3/generated/control.TransferFunction.html}

新一屆培訓(xùn)開始啦

新一屆的培訓(xùn)開始啦，希望卓大大規(guī)則也可以今早出臺(tái)奧。

為祖國獻(xiàn)禮，祝我們的祖國生日快樂，繁榮富強(qiáng)！