简 介： 根据学生提出的如何使用C语言来实现系统函数的问题，给出了利用MATLAB，Python对连续时间系统函数进行离散化的方法。继而可以进行C语言实现。
关键词： 系统函数，离散化，MATLAB
?
§01 如何实现传递函数？
??卓晴老师，请教一下，传递函数环节用C需要编写，该怎么编写？在matlab中他能简单的表述，但是在利用C语言，该咋办呢？
??例如写进单片机中。类似于写一个传递函数通式，想实现各种环节。
??在matlab中他有很好的封装，却不好理解matlab怎么做到的，simulink中一个时域信号输入传递模块输出还是时域信号，逐步仿真我认为这个传递模块中封装着一个复杂表达式，想知道用C语言咋能实现，这样以后对信号处理想加个滤波环节，一阶延时，可以利用这种方法直接解决。
??好的，感谢卓大大，最好给出一个G/TS+1这种的一个C代码的例子，然后我就可以仿照思路写出更多高阶的。
?
§02 什么是传递函数？
??要回答前面同学提出的“如何使用C编程语言实现传递函数？”，则需要先了解什么是系统传递函数？为什么引入系统传递函数？
一、线性时不变系统
??下面两个由线性R，C组成的电路网络，它们的共同之处都属于 线性时不变系统 ，不同之处在于左边的电路输入输出之间是一个比例关系， U O ( t ) U I ( t ) = R 12 R 11 + R 12 {{U_O \left( t \right)} \over {U_I \left( t \right)}} = {{R_{12} } \over {R_{11} + R_{12} }} UI?(t)UO?(t)?=R11?+R12?R12??
??这类系统被称为即时系统。
▲ 图2.1.1 左：电阻分压电路；右：阻容低通滤波电路
??右边的低通滤波电路输入输出之间关系则需要通过微分方程才能够描述： R 21 C 11 ? d V O ( t ) d t + V O ( t ) = V I ( t ) R_{21} C_{11} \cdot {{dV_O \left( t \right)} \over {dt}} + V_O \left( t \right) = V_I \left( t \right) R21?C11??dtdVO?(t)?+VO?(t)=VI?(t)
??这种带有储能器件的电路，当前的输出不仅与当前的输入有关系，还和之前时间输入有关系。这类系统被称为动态系统，显然动态系统比即使系统复杂，但功能更强大。即使系统可以看成一类特殊的动态系统。
二、系统函数
??既然两个电路都属于线性时不变系统，可否将它们的输入输出关系都简化成比例关系，这样分析与求解就比较方便了？
??答案是可以的，只不过需要把输入输出信号都进行某种变换，比如 傅里叶变换 ， 拉普拉斯变换 ，这时变换后的系统输入输出之间就是一个比例关系了。这个比值就被称为线性时不变分析对象的系统函数。
??在拉普拉斯变换下，上面的电阻分压电路对应的系统函数为： F 1 ( s ) = U O ( s ) U I ( s ) = R 12 R 11 + R 12 F_1 \left( s \right) = {{U_O \left( s \right)} \over {U_I \left( s \right)}} = {{R_{12} } \over {R_{11} + R_{12} }} F1?(s)=UI?(s)UO?(s)?=R11?+R12?R12??
??阻容低通滤波电路对应的系统函数为：
F
2
(
s
)
=
U
O
(
s
)
U
I
(
s
)
=
1
R
21
C
11
s
+
1
=
1
τ
?
s
+
1
,
?????
τ
=
R
21
C
11
F_2 \left( s \right) = {{U_O \left( s \right)} \over {U_I \left( s \right)}} = {1 \over {R_{21} C_{11} s + 1}} = {1 \over {\tau \cdot s + 1}},\,\,\,\,\,\tau = R_{21} C_{11}
F2?(s)=UI?(s)UO?(s)?=R21?C11?s+11?=τ?s+11?,τ=R21?C11?
??如果抛开这些变换概念背后的数学基础，从工程应用角度也可以将上述系统函数看成对描述线性时不变系统输入输出关系之间的微分方程的另外一种表示方式。
??借助于系统函数的概念，就把线性时不变系统的分析、设计都可以进行简化。
▲ 图2.2.1 在变换域内定义的系统函数等于系统的零状态响应与输入信号的比值
三、离散时间系统函数
??计算机的普及使得离散时间系统应用增加，系统可以通过软件编程来实现。如果实现对于输入信号进行低通滤波，也可以在对信号进行采样之后利用下面离散时间系统来进行处理。
▲ 图2.3.1 一阶低通离散时间系统
??下图中显示的随机噪声的原始信号，经过上面迭代方程滤波之后，就会形成平滑后的数据。其中的滤波器系数 a = 0.99 a = 0.99 a=0.99。
▲ 图2.3.1.1 增加有噪声信号以及一阶低通滤波后的信号
滤波系数a=0.99
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()
??线性时不变动态离散时间系统可以使用差分方程来刻画，经过简单的一项便可成写出的迭代方程形式，即由输入信号以及以前的输出信号计算出当前系统输出值，这是利用软件编程实现的基础。
??同样，在对于系统的输入输出信号进行 Z-变换 ，那么系统的零状态输出信号比上输入信号也是一个常量，定义为离散时间系统的系统函数。上述一阶低通离散时间系统对应的系统函数为： H ( z ) = 1 ? a 1 ? a z ? 1 = ( 1 ? a ) z z ? a H\left( z \right) = {{1 - a} \over {1 - az^{ - 1} }} = {{\left( {1 - a} \right)z} \over {z - a}} H(z)=1?az?11?a?=z?a(1?a)z?
??同样，如果抛开z变换背后的数学原理，仅仅把 z ? 1 z^{ - 1} z?1看成对于信号的延迟操作，那么上面的系统函数也可以看成差分方程另外一种描述方法。可以看到系统函数所对应的有理分式的分子多项式系数对应着差分方程输入信号前的系数；分母多项式的系数对应差分方程输出信号前的系数。
▲ 图2.3.2 离散时间系统的系统方程
??因此，如果已知知道了离散时间系统函数，则按照上面的系统函数与差分方程对应关系可以写出对应的差分方程，进而整理成对应的迭代方程。这样就可以使用编程语言（C,Python）来实现算法了。
?
§03 如何将H(s)转换成H(z)?
??本文一开始同学询问，如果知道了一个连续时间系统的传递函数，怎么使用C语言编程实现系统功能？这个问题就转换成如何将连续时间系统函数 H ( s ) H\left( s \right) H(s)转换成对应的离散时间系统函数 H ( z ) H\left( z \right) H(z)的问题了。这个过程也称为连续时间系统的离散化。
一、连续系统函数离散化
1、离散化方法
??将连续时间系统转换成对应的离散时间系统的过程常常用在模拟设计-离散实现的设计过程中。通常使用的方法包括：
- 冲激响应不变法：这种方法可以保持离散时间系统的单位冲激响应是对应连续时间系统的采样；
- 变量替换法：包括前向差分（Euler方法）、后向差分、双线性法，也称梯形方法；
??除此之外，在 Continuous 中给出了在控制系统中常用到的六种方法及其转换原理。
2、举例
??比如使用双线性方法，它是基于把代表位移算子的
z
z
z所对应的连续系统函数近似展开成有理分式：
z
=
e
s
T
s
≈
1
+
s
T
s
/
2
1
?
s
T
s
/
2
,
???
s
=
2
T
s
?
z
?
1
z
+
1
z = e^{sT_s } \approx {{1 + sT_s /2} \over {1 - sT_s /2}},\,\,\,s = {2 \over {T_s }} \cdot {{z - 1} \over {z + 1}}
z=esTs?≈1?sTs?/21+sTs?/2?,s=Ts?2??z+1z?1?
??以一个一阶低通滤波器为例： H ( s ) = G T s + 1 H\left( s \right) = {G \over {Ts + 1}} H(s)=Ts+1G?
??使用双线性方法转换成离散系统函数：
H
(
z
)
=
G
2
T
T
s
?
z
?
1
z
+
1
+
1
=
G
(
1
+
z
?
1
)
(
2
T
T
s
+
1
)
?
2
T
T
s
z
?
1
H\left( z \right) = {G \over {{{2T} \over {T_s }} \cdot {{z - 1} \over {z + 1}} + 1}} = {{G\left( {1 + z^{ - 1} } \right)} \over {\left( {{{2T} \over {T_s }} + 1} \right) - {{2T} \over {T_s }}z^{ - 1} }}
H(z)=Ts?2T??z+1z?1?+1G?=(Ts?2T?+1)?Ts?2T?z?1G(1+z?1)?
二、使用MATLAB离散化系统函数
??利用MATLAB内部函数 c2d ，可以很方便将连续时间系统函数转换成对应的离散时间系统函数。
1、命令语法
??MATLAB 中给出了一下四种常用到离散系统函数命令形式：
sysd = c2d(sysc,Ts)
sysd = c2d(sysc,Ts,method)
sysd = c2d(sysc,Ts,opts)
[sysd,G] = c2d(___)
??第一种方法利用了零阶保持方法完成系统函数离散化，如果希望指明离散化具体方法，可以使用2,3两种命令形式。
1、转换举例
（1）连续时间系统函数
??将如下带有延时环节的二阶连续时间系统进行离散化：
H ( s ) = e ? 0.25 s 10 s 2 + 3 s + 10 H\left( s \right) = e^{ - 0.25s} {{10} \over {s^2 + 3s + 10}} H(s)=e?0.25ss2+3s+1010?
（2）转换命令
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）系统仿真
??通过以下指令，可以对比连续时间系统与离散时间系统单位冲激响应波形。
step(h,’–’,hd’,’-’)
▲ 图3.2.2.1 系统的单位冲激响应信号对比
?
※ 问题总结 ※
??使用嵌入式系统实现信号处理算法、控制算法，比起利用连续时间系统（或者模拟电路）实现具有很大的优点。在一些模拟设计（仿真）-数字化实现的场景中，将连续时间系统函数转换成对应的离散时间系统函数是编程实现控制算法的第一步。
??离散化系统函数的方法有很多种，在MATLAB中可以利用相关的命令完成这些转换，并进行系统仿真测试。
??手工实现连续系统传递函数的一般过程：
- 将连续时间系统函数转换成离散时间系统函数；
- 将离散时间系统函数整理成关于 z ? 1 z^{ - 1} z?1的有利多项式；
- 根据 z ? 1 z^{ - 1} z?1有理多项式写对应的的差分方程；
- 将差分方程整理成迭代方程形式，便可以通过C语言完成其中的运算了。
??这个过程也可以借助于 MATLAB 代码转成C语言代码 的辅助功能来加速，甚至在一些硬件平台上，生成的C语言代码可以直接部署在嵌入式单片机中来执行。
??在 control.TransferFunction 中给出了利用 Python 中的 control 开发包进行连续时间系统离散化相关指令。
??下面示例给出了给出了将系统函数： H ( s ) = 1 s + 1 H\left( s \right) = {1 \over {s + 1}} H(s)=s+11?转换成对应的离散时间系统函数： H ( z ) = 0.2 z + 0.2 z ? 0.6 , ??? d t = 0.5 H\left( z \right) = {{0.2z + 0.2} \over {z - 0.6}},\,\,\,dt = 0.5 H(z)=z?0.60.2z+0.2?,dt=0.5
import control
sys = control.TransferFunction(1,[1,1])
sysd = sys.sample(0.5, method='bilinear')
■ 相关文献链接:
● 相关图表链接: