离散时间系统入门 | 硬件学习笔记

19.2 离散时间系统入门

上一节，我们揭开了数字环路的面纱——量化像台阶一样粗糙，延迟像幽灵一样存在。这些非理想因素把原本连续、平滑的模拟控制切得支离破碎。为了在这个破碎的世界里重建稳定，我们手里需要一把新的锤子。这把锤子，就是离散时间系统理论。

本节的任务，就是把这把锤子递到你的手里。我们将看到如何在离散的数字域里，用加减乘除和单元延迟，去重建那些我们在模拟域里熟稔的积分、比例和微分运算。

19.2.1 连续积分与离散积分

让我们从最熟悉的场景开始。一个标准的模拟控制环路你早就见过：传感器采回来输出电压，减去参考值得到误差，再丢进那个连续域的补偿器 $G_{c} (s)$ 。我们在第 9 章花了大量篇幅学习怎么设计它。

考虑一个最简单的纯积分补偿器（公式 19.15）：

G_{c} (s) = \frac{v_{c} (s)}{v_{e} (s)} = \frac{ω_{o}}{s}

在 $s$ 域里，它是一个位于原点的极点。这意味着在时域上，输出 $v_{c} (t)$ 是输入 $v_{e} (t)$ 的积分（公式 19.16）：

v_{c} (t) = v_{c} (0) + ω_{o} \int_{0}^{t} v_{e} (τ) d τ

想象一下这两个信号的波形：在这个例子里，输入信号 $v_{e} (t)$ 是一个正弦波，频率为 $f_{s} / 10$ （ $f_{s}$ 是采样频率），输出 $v_{c} (t)$ 则是它的积分——一条带斜率的、被「抹平」了的曲线。

现在的关键问题是：在我们那个数字控制器里，这个积分到底该怎么算？ 我们手里只有一堆离散的采样点 $v_{e} [n] = v_{e} (n T_{s})$ 。我们要算出离散的输出点 $v_{c} [n]$ ，让它尽可能接近连续积分器在同一时刻的值 $v_{c} (n T_{s})$ 。

【模式 B：侯捷紧密模式】 先把连续积分公式 (19.16) 变个形。我们把积分区间切开，切成一个个长度为 $T_{s}$ 的小段。

v_{c} (t) = v_{c} (t - T_{s}) + ω_{o} \int_{t - T_{s}}^{t} v_{e} (τ) d τ (19.17)

要想在数字域精确复现公式 (19.17)，补偿器得算出这个积分：也就是算出 $v_{e} (t)$ 在 $t = (n - 1) T_{s}$ 到 $t = n T_{s}$ 这段时间曲线下的面积（公式 19.18）：

v_{c} [n] = v_{c} [n - 1] + ω_{o} \int_{(n - 1) T_{s}}^{n T_{s}} v_{e} (τ) d τ

但这里有个硬伤：我们根本不知道曲线长什么样，只有头尾两个采样点 $v_{e} [n]$ 和 $v_{e} [n - 1]$ 。所以，精确积分是不可能的。我们只能近似。

最常用的近似方法，是梯形近似（Trapezoidal Approximation）。思路很简单：把相邻两个采样点连起来，用这俩点撑起的梯形面积，去代替曲线底下那一段真实的面积。

v_{c} [n] = v_{c} [n - 1] + ω_{o} T_{s} \frac{v_{e} [n] + v_{e} [n - 1]}{2} (19.19)

这个公式非常有价值。它把微积分变成了简单的算术：一个加法（取平均），一个乘法（乘系数），再加到上一时刻的结果上。这正是 CPU 或 FPGA 擅长的事情。

把离散积分器算出来的 $v_{c} [n]$ （一个个采样点）和模拟积分器的真实值 $v_{c} (t)$ （连续曲线）放在一起比，你会发现它们贴得很近，但不完全重合。随着采样频率 $f_{s}$ 升高、 $T_{s}$ 变小，这个误差会迅速减小。

当然，梯形近似不是唯一的路。如果你把面积近似成矩形，还有两种常见的“偷懒”方法：

后向欧拉法（公式 19.20）：只看上一个时刻的高度。 $v_{c} [n] = v_{c} [n - 1] + ω_{o} T_{s} v_{e} [n - 1]$
前向欧拉法（公式 19.21）：只看当前时刻的高度。 $v_{c} [n] = v_{c} [n - 1] + ω_{o} T_{s} v_{e} [n]$

这三种方法都有人用，但通常来说，梯形近似是最准的。为什么？因为它抓住了“平均值”的本质。

19.2.2 Z 变换与频率响应

公式 (19.19), (19.20), (19.21) 给出了时域里的递推公式。但在前几章里，我们已经习惯了用拉普拉斯变换、传递函数和波特图来思考问题。在离散域里，对应的家伙是谁？是 Z 变换。

它是离散时间系统的“拉普拉斯变换”。

对于一个离散信号 $x [n]$ ，Z 变换定义为（公式 19.22）：

Z {x [n]} = X (z) = \sum_{n = - \infty}^{+ \infty} x [n] z^{- n}

它和拉普拉斯变换一样，也是线性的（公式 19.23）：

Z {a x [n] + b y [n]} = a X (z) + b Y (z)

最核心的性质来了（公式 19.24）。如果你把信号 $x [n]$ 延迟一个采样周期，变成 $x [n - 1]$ ，它的 Z 变换会怎样？

Z {x [n - 1]} = \sum_{n = - \infty}^{+ \infty} x [n - 1] z^{- n}

做变量替换 $k = n - 1$ ：

= z^{- 1} \sum_{k = - \infty}^{+ \infty} x [k] z^{- k} = z^{- 1} X (z)

结论很简单，但极其重要：时域延迟一个周期 $T_{s}$ ，等于在 Z 域乘以 $z^{- 1}$ 。换句话说， $z^{- 1}$ 就是离散域里的“算子”，代表单位延迟。

好了，有了这个武器，我们回过头看那个梯形近似的积分器公式 (19.19)。两边做 Z 变换（利用线性性质和延迟性质）：

V_{c} (z) = z^{- 1} V_{c} (z) + ω_{o} T_{s} \frac{V_{e} (z) + z^{- 1} V_{e} (z)}{2} (19.25)

把它整理成传递函数的形式 $G_{c d} (z) = V_{c} (z) / V_{e} (z)$ ，我们就得到了离散积分器的 Z 域表达式（公式 19.26）：

G_{c d} (z) = \frac{V_{c} (z)}{V_{e} (z)} = \frac{ω_{o} T_{s}}{2} \frac{1 + z^{- 1}}{1 - z^{- 1}} = ω_{o} T_{s} \frac{z + 1}{2 (z - 1)}

我们把这个结果列在表 19.2 里，你可以看到前向欧拉和后向欧拉的结果也在那儿。

表 19.2 离散积分器传递函数对比

近似方法	$G_{c d} (z)$
梯形近似	$\frac{ω_{o} T_{s}}{2} \frac{z + 1}{z - 1}$
后向欧拉	$\frac{ω_{o} T_{s}}{z - 1}$
前向欧拉	$\frac{ω_{o} T_{s} z}{z - 1}$

【模式 A：SICP 紧绷模式】 现在我们站在了一个分岔路口。我们有了离散域的传递函数，但怎么评估它的稳定性？在模拟域，我们把 $s$ 换成 $j ω$ 来看频率响应。在 Z 域，我们该怎么做？

回想一下刚才的结论：

Z 域的 $z^{- 1}$ 是单位延迟。
拉普拉斯域里，延迟 $T_{s}$ 对应的因子是 $e^{- s T_{s}}$ （公式 19.27）。

这意味着，我们可以建立这样一个映射：

z^{- 1} \leftrightarrow e^{- s T_{s}} ⟹ z \leftrightarrow e^{s T_{s}}

所以，求 $G_{c d} (z)$ 的频率响应，方法就是（公式 19.28）：

G_{c d} (j ω) = G_{c d} (z) |_{z \to e^{j ω T_{s}}}

我们来试一下。把那个梯形积分器的公式 (19.26) 代进去（公式 19.29）：

G_{c d} (j ω) = ω_{o} T_{s} \frac{e^{j ω T_{s} / 2} + e^{- j ω T_{s} / 2}}{2 (e^{j ω T_{s} / 2} - e^{- j ω T_{s} / 2})}

利用欧拉公式 $e^{j x} = \cos x + j \sin x$ ，上下同除，最后化简得到（公式 19.30）：

G_{c d} (j ω) = - j ω_{o} T_{s} \frac{\cos (ω T_{s} / 2)}{\sin (ω T_{s} / 2)}

这东西看起来很怪。我们来对比一下原版模拟积分器的响应（公式 19.31）：

G_{c} (j ω) = - j \frac{ω_{o}}{ω}

惊人的事情发生了。如果你看相位，公式 (19.30) 在所有频率下的相位都是 $- 90^{\circ}$ 。它和模拟积分器完全一致。注意：这只是梯形近似的特权。前向和后向欧拉法做不到这一点。

我们再来看看幅值。假设频率很低，满足 $ω T_{s} / 2 ≪ 1$ （也就是说信号频率远小于采样频率 $f_{s}$ ）。这时候 $\tan (x) \approx x$ ，公式 (19.30) 就变成了（公式 19.32）：

| G_{c d} (j ω) | \approx ω_{o} T_{s} \frac{1}{ω T_{s} / 2} = \frac{ω_{o}}{ω} = | G_{c} (j ω) |

这说明：只要采样频率足够高，在低频段，离散积分器几乎完美复现了模拟积分器。

但是，随着频率升高，接近奈奎斯特频率（ $f_{s} / 2$ ）时，差距会越来越大。对比模拟积分器和离散积分器的幅频曲线（模拟那条是平滑的，离散那条是周期波动的），你会看到离散曲线在高频段越走越歪。更致命的是，离散积分器在某些频率点的增益会掉到零（公式 19.33）：

| G_{c d} (j ω) | = 0, 当 f = \frac{f_{s}}{2}, \frac{3 f_{s}}{2}, \dots

这在模拟积分器里是不可想象的—— $s = 0$ 处是极点，增益应该是无穷大才对。但在离散域，周期性导致了零点的出现。

【阶段收尾 - 喘息】 这里我们看到了离散域的第一个代价：高频特性的畸变。但在控制环路里，我们的穿越频率通常远小于 $f_{s} / 2$ 。所以，只要别去触碰那个高频禁区，我们依然可以用 Z 变换来设计。

💡 一句话记住梯形近似为什么「最准」：前向、后向欧拉各自只偷看一个采样点（左端点或右端点），本质是拿矩形去糊弄曲线面积；梯形近似左右两点都看，等价于把那段曲线当线性插值——只要采样够密，曲线段本来就近似直线，所以梯形几乎是「白捡」的精度。这也是为什么双线性映射在相位上能撞上 $- 90^{\circ}$ 、而两种欧拉法做不到——精度不是天上掉下来的，是「多用一个采样点信息」换的。

接下来的问题是，对于一个更复杂的补偿器（比如 PID），我们怎么从 $s$ 域优雅地跳到 $z$ 域？

19.3.3 连续域到离散域的映射

我们在上一节里，用梯形积分的例子推导出了一个从 $s$ 到 $z$ 的映射关系。从公式 (19.34) 可以看出这个端倪：

G_{c} (s) = \frac{ω_{o}}{s} \overset{离散化}{\to} G_{c d} (z) = ω_{o} T_{s} \frac{z + 1}{2 (z - 1)}

如果我们把这个 $s$ 和 $z$ 的关系解出来，就会发现一个通用的置换公式（公式 19.35）：

s \to \frac{2}{T_{s}} \frac{z - 1}{z + 1}

这就是大名鼎鼎的双线性映射，或者叫 Tustin 映射。只要把你设计好的模拟补偿器 $G_{c} (s)$ 里的 $s$ 全部换成这个关于 $z$ 的分式，你就得到了离散形式的 $G_{c d} (z)$ （公式 19.36）。

G_{c d} (z) = G_{c} (s) |_{s \to \frac{2}{T_{s}} \frac{z - 1}{z + 1}}

这个映射到底干了什么？它的秘密藏在 $s$ 平面和 $z$ 平面的几何对应里。把公式 (19.35) 变形一下，解出 $z$ （公式 19.37）：

z = \frac{1 + s T_{s} / 2}{1 - s T_{s} / 2}

$s$ 平面的原点（ $s = 0$ ，积分器极点），跑到了 $z$ 平面的 $z = 1$ 。
$s$ 平面的虚轴（ $s = j ω$ ），跑到了 $z$ 平面的单位圆上（ $| z | = 1$ ）。
$s$ 平面的左半平面（稳定区），跑到了 $z$ 平面的单位圆内部。

这太完美了。它保证了稳定的模拟系统映射后还是稳定的离散系统。

【模式 B：侯捷紧密模式】 我们来实战一下，把我们在第 9 章熟悉的 PI 补偿器搬过来。模拟 PI 补偿器长这样（公式 19.38）：

G_{c} (s) = G_{c \infty} (1 + \frac{ω_{L}}{s})

用双线性映射公式代入 $s$ 。经过一番代数化简（公式 19.39），我们可以把它写成极零点的形式（公式 19.40）：

G_{c d} (z) = G_{c \infty} \frac{z - (1 - ω_{L} T_{s} / 2)}{z - 1} \dots

通常 PI 补偿器的零点频率 $f_{L}$ 很低，满足 $ω_{L} T_{s} / 2 ≪ 1$ 。公式可以简化为（公式 19.41）：

G_{c d} (z) \approx G_{c \infty} \frac{z - (1 - ω_{L} T_{s})}{z - 1}

这里看到了一个非常有意思的现象： 离散零点跑到了 $z \approx 1$ 的位置。当 $f_{L}$ 很低的时候，零点甚至非常接近极点（ $z = 1$ ）。在 $z$ 平面上，两个点凑得太近会带来数值计算的麻烦，这在后面的实现部分（19.4 节）会是个坑。

在低频段（远低于 $f_{s}$ ），模拟的 $G_{c} (s)$ 和离散的 $G_{c d} (z)$ 的幅值和相位几乎重合。这就是双线性映射的魅力。

再来试试 PD 补偿器（公式 19.43）：

G_{c} (s) = G_{c 0} \frac{1 + s / ω_{z}}{1 + s / ω_{p}}

同样代入双线性映射公式。我们会发现，零点和极点都被映射到了 $z$ 平面的实轴上。但是，如果原来的零极点频率很高（接近 $f_{s}$ ），麻烦就来了。

⚠️ 踩坑预警 如果 $f_{z}$ 和 $f_{p}$ 并不远小于 $f_{s}$ ，双线性映射带来的频率畸变会非常明显。离散那条曲线在高频段会明显偏离模拟那条——相位超前刚爬到峰值，紧接着就一头栽下去。这时候，我们得用一招叫频率预畸变（Frequency Prewarping）的高级技巧。

它的核心思想是：既然映射会把频率搞弯，那我就事先把那个关键的频率点“弯”回去。我们在映射公式里加个系数 $k_{p r e w a r p}$ （公式 19.46, 19.47）：

s \to k_{p r e w a r p} \frac{2}{T_{s}} \frac{z - 1}{z + 1}

这个系数是精心设计的，目的是保证在某个特定频率 $ω_{p r e w a r p}$ （通常是穿越频率）上，离散和模拟的响应完全一致（公式 19.48, 19.49）。

用了预畸变之后，你会看到在 $f_{p r e w a r p} = 200 kHz$ 处，离散和模拟两条曲线被「焊」在一起，分毫不差。代价是低频段的误差稍微大了一点点，但这通常是值得的。

最后，作为一个总结，我们来看看 PID 补偿器（公式 19.50）。在模拟域，PID 通常有一个高频极点 $f_{p 2}$ 用来衰减开关纹波。但在数字域，A/D 采样本身就是个低通过程，开关纹波已经被滤掉了。所以，我们在设计 $G_{c d} (z)$ 的时候，通常不需要把那个高频极点 $f_{p 2}$ 映射进来。那个任务交给模拟侧的抗混叠滤波器 $H (s)$ 去做就行。

直接用带预畸变的双线性映射，得到的离散 PID 传递函数 $G_{c d}^{*} (z)$ 形如（公式 19.51）：

G_{c d}^{*} (z) = G_{d} \frac{(z - z_{L}) (z - z_{z})}{(z - 1) (z - z_{p})}

这里的系数 $z_{L}, z_{z}, z_{p}$ 都是给定的公式算出来的（公式 19.52 - 19.54），看着复杂，其实就是映射后的极零点位置。

【模式 C：费曼流动模式】 说实话，手动去算这些 $z$ 域的极零点简直是灾难。好在现在的工具都很聪明。MATLAB 里的 c2d 函数（表 19.3）一键就能搞定这些繁琐的代数运算。

表 19.3 离散化映射工具

映射方法	替换公式	MATLAB 函数
双线性	$s \to \frac{2}{T_{s}} \frac{z - 1}{z + 1}$	`c2d(Gc, Ts, 'tustin')`
带预畸变的双线性	$s \to k_{p r e w a r p} \frac{2}{T_{s}} \frac{z - 1}{z + 1}$	`c2d(Gc, Ts, 'prewarp', wprewarp)`

至此，我们完成了从理论到工具的准备。我们已经知道怎么把模拟补偿器“翻译”成数字算法。但“翻译”仅仅是第一步。当你把这行代码真的塞进单片机里跑的时候，你会发现有些坑是数学公式里看不出来的—— 比如，当积分值大到连 uint32_t 都装不下时会发生什么？比如，当 DPWM 的分辨率不够，你想调 1% 占空比却只能调 2% 时会怎样？这就是下一节（19.3 和 19.4）我们要面对的残酷现实。

参考说明：参考自 geqianQWQ 同学阅读《Fundamentals of Power Electronics》的笔记，仅作理解线索；本文为结合自己理解重新整理的学习笔记，不涉及对原书的复制或翻译。

19.2 离散时间系统入门 ​

19.2.1 连续积分与离散积分 ​

19.2.2 Z 变换与频率响应 ​

19.3.3 连续域到离散域的映射 ​

19.2 离散时间系统入门

19.2.1 连续积分与离散积分

19.2.2 Z 变换与频率响应

19.3.3 连续域到离散域的映射