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

卓晴老師，請(qǐng)教一下，傳遞函數(shù)環(huán)節(jié)用C需要編寫，該怎么編寫？在matlab中他能簡(jiǎn)單的表述，但是在利用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)系都簡(jiǎ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)行簡(jiǎ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)過簡(jiǎn)單的一項(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)行離散化：

$$Hleft( s right) = e{ - 0.25s} {{10} over {s2  + 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)的場(chǎng)景中，將連續(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)仿真測(cè)試。

在 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/%E7%BA%BF%E6%80%A7%E6%97%B6%E4%B8%8D%E5%8F%98%E7%B3%BB%E7%BB%9F/2474061

[2]傅里葉變換: https://mathworld.wolfram.com/FourierTransform.html#::text=The%20smoother%20a%20function%20%28i.e.%2C%20the%20larger%20the,Fourier%20transform%20is%20also%20symmetric%20since%20implies%20.

[3]拉普拉斯變換: https://byjus.com/maths/laplace-transform/#::text=The%20Laplace%20transform%20is%20a%20well%20established%20mathematical,to%20calculate%20the%20solution%20to%20the%20given%20problem.

[4]Z-變換: https://faculty.ksu.edu.sa/sites/default/files/dsp_cen352_ch5_z.pdf#:~:text=Definition%20of%20z-transform%3A%20The%20z-transform%20is%20a%20very,filter%20design%20and%20frequency%20analysis%20of%20digital%20signals.

[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