手里攥着一个s域传递函数或者一版在频域里设计好的滤波器参数结果往单片机上一放发现根本跑不起来——做控制、做信号处理的人多多少少都撞过这堵墙。问题不在算法本身而在连续域那套微积分语言芯片压根不认它只认加减乘和移位只认一个接一个在时间轴上排好队的采样点。把连续世界的方程翻译成芯片能执行的迭代式子中间搭的桥叫Z变换桥那头落地的成品叫差分方程。这篇文章要讲的正是这条链路Z变换方程怎么一步步变成差分方程系数到底从哪来位置式PID为什么常被写成离散差分形式采样周期该按什么依据定代码里怎么写才不会发散。不管你是刚学完信号与系统还对着z的负幂发愣的学生还是手头有个实际控制器等着上板子的工程师这套推导、这些参数和这些坑都值得从头到尾走一遍。1. 为什么非得把Z变换方程翻译成差分方程1.1 一个真实场景算出来的公式跑不进单片机拿电机电流环里最常见的一阶模型举例连续传递函数写成 G(s) 1/(0.05s 1)。在纸面上分析它特别舒服波特图、相位裕度、阶跃响应全都能算得漂漂亮亮。可单片机的处境完全不同它手里只有一个定时器中断每1毫秒触发一次进来读一次ADC算一次输出然后退出下一次中断再从零开始。它没法表达ds没法表达无穷小的时间也没法表达连续时间t。芯片能操作的只有离散序号n以及第n拍、第n-1拍存下来的数值。所以必须做一次坐标转换把连续时间t换成离散序号n把微分算子换成差分把拉普拉斯算子s换成z。这个转换做完纸面上的G(s)就变成了一个只含乘法和加法的递推关系芯片每一拍照着算就行。这里有个容易被忽略的点这不是数学上的近似游戏而是一个采样-保持的物理过程在建模。DAC输出后一般保持一个周期不变零阶保持ADC也是每隔T秒抓一次瞬时值Z变换恰好就是描述这种离散取样保持行为的数学工具。理解了这一层后面z的负幂为什么代表延迟、为什么要挑采样周期就都顺理成章了。1.2 差分方程的本质是延迟加乘加差分方程长什么样典型形式是 y[n] 等于一堆历史值的线性组合比如 y[n] a1·y[n-1] a2·y[n-2] b0·x[n] b1·x[n-1] b2·x[n-2]。仔细看这个式子右边出现的全是已知量历史输出a1、a2对应的y值是上一拍、上上拍已经算完存下来的输入x的当前值和历史值也都在数组里唯一未知的只有左边这一拍的y[n]。这意味着每拍的计算量就是固定几次乘加不用解方程、不用矩阵求逆、不用迭代收敛。这正是它能上硬件的根本原因。MCU里的单周期乘加指令、DSP里的MAC单元、FPGA里的DSP48硬核天生就是干乘一下累加一下这件事的。所以无论是数字滤波器、PID控制器还是观测器最终几乎都会被整理成差分方程这种统一形式。反过来说如果你拿到的是一个Z变换表达式却没有把它整理成当前输出 历史输出和历史输入的加权和那它离能写进代码还差最后也是最关键的一步而这一步恰恰是很多人卡住的地方。1.3 三条离散化路径怎么选从s域到z域不是只有一条路常见的有前向差分、后向差分和双线性变换也叫梯形法、Tustin法。选哪条直接决定了最终系统的稳定性、频率特性和实现难度不能随手抓一个就用。下面这张表把三条路的关键差异摆在一起方便对照着挑。方法s的替换式稳定性频率映射实现难度典型用途前向差分s (1 - z) / T可能把稳定系统变成不稳定畸变较大最简单教学演示实际少用后向差分s (1 - z⁻¹) / T稳定域映射后仍稳定低频准高频压缩简单快速落地、对精度要求一般双线性s (2/T)·(1 - z⁻¹)/(1 z⁻¹)稳定域一一对应需预畸变映射好中等滤波器、精度要求高的控制后向差分的最大好处是无条件稳定只要原来的连续系统是稳定的离散化之后一定还稳定这对工程落地非常友好代价是频率轴被非线性压缩高频段误差明显。双线性变换把整个左半s平面映射到z平面单位圆内稳定性和频率对应都很好缺点是频率轴有畸变需要做预畸变处理而且会引入一个前向的零点阶跃响应可能带一点小过冲。前向差分因为可能把一个稳定极点推到单位圆外实际工程里基本只在理论推导时露个脸真正上板子很少用它。我的习惯是控制器类先用后向差分快速跑通滤波器类或者对截止频率精度要求高的场合改用双线性并老老实实做预畸变。2. 转换的核心法则把z⁻¹当成一拍延迟2.1 从Z变换定义到位移定理要把Z变换方程翻成差分方程前提是真正理解z⁻¹的物理含义。Z变换的定义是 X(z) Σ x[n]·z⁻ⁿ这里的z在纯数学上是个复变量。但对于因果离散序列z⁻¹实际扮演的角色是一个延迟算子把它乘在某个序列上就等于把这个序列整体往后推一拍。这个结论来自位移定理也是整条转换链路的地基值得亲手推一遍。设序列x[n]的单边Z变换为X(z)考虑延迟k拍的序列y[n] x[n-k]它的Z变换是 Σ(n从k到∞) x[n-k]·z⁻ⁿ。做一次下标替换令m n-k则n mk求和变成 Σ(m从0到∞) x[m]·z⁻⁽ᵐ⁺ᵏ⁾ z⁻ᵏ·Σ x[m]·z⁻ᵐ z⁻ᵏ·X(z)。整个过程只用到指数运算的拆分没有别的花招。推完得到的结论非常直白在z域里乘以z⁻ᵏ等价于在时域里延迟k个采样周期。这条定理一旦吃透Z变换方程转差分方程就变成了一个机械的替换动作误差也主要出在别的地方身上。2.2 常用的Z变换对与延迟映射表实际动手时不需要每次都从定义推把常用的映射关系背下来、做成对照表效率最高。下面这张表列的是推导过程中最常打交道的几个连续环节和它们对应的z域表达、离散差分项左边是s域右边是可直接写进代码的z域形式和时域递推关系。s域环节含义z域替换后向差分时域离散形式s微分(1 - z⁻¹) / T(x[n] - x[n-1]) / T1/s积分T / (1 - z⁻¹)T · Σ x[i]e^(-sT)纯延迟z⁻¹x[n-1]1/(τs1)一阶惯性KT/(τT) / (1 - a·z⁻¹)见第3节推导a/(sa)一阶低通......这张表的用法很直接拿到一个由若干个基本环节串并联组成的连续传递函数逐项把s和1/s按表中规则换掉再通分整理成Y(z)/X(z) 多项式之比最后把z⁻¹对应成延迟拍数差分方程就出来了。这里要特别提醒一句逐项替换只对低阶、结构简单的系统有效高阶系统如果拆开逐项换再相乘误差会累积正确做法是先得到完整的G(s)再整体做替换然后统一代入。2.3 零初始条件这个默认前提不能丢几乎所有工程推导在写位移定理时都默认零初始条件也就是系统从静止状态开始所有历史值最初都是零。单边Z变换的完整位移定理其实带初值项比如 Z{x[n-1]} z⁻¹X(z) x[-1]那个x[-1]就是初值。工程上忽略它是因为大部分系统上电时状态确实是零忽略后式子干净很多。但这个前提不能忘一旦场景变了它还真的会咬人。什么时候初值项会跳出来捣乱比如控制器运行中途被人为复位、从某个非零稳态重新启动、或者你把一段已经跑了一段时间的数据截取出来单独处理这些情况下x[-1]、x[-2]不为零直接套简化公式就会得到错误结果。还有一个更隐蔽的坑有些代码在中断里每次都把历史数组清零表面上看是干净实际上破坏了差分方程的连续性会导致输出在复位瞬间出现跳变。正确的做法是把历史值作为控制器状态保留下来只在系统真正重新上电时才清零。3. 从传递函数到差分方程的完整推导3.1 一阶惯性环节三种方法的系数对比拿最容易上手的一阶惯性环节开刀连续传递函数 G(s) K/(τs 1)。用后向差分代入 s (1 - z⁻¹)/T得到 G(z) K / [τ(1 - z⁻¹)/T 1]。把分母通分整理分子分母同时乘以T得 G(z) KT / [τ(1 - z⁻¹) T] KT / (τ T - τz⁻¹)。再上下同除(τT)化成标准形式 G(z) [KT/(τT)] / [1 - (τ/(τT))·z⁻¹]。把这个结果翻译成差分方程就一目了然了y[n] a·y[n-1] b·x[n]其中 a τ/(τT)b KT/(τT)。代入具体数字验证一下取K1、τ0.05s、T0.001s则 a 0.05/0.051 ≈ 0.980392b 0.001/0.051 ≈ 0.019608。做个自检稳态增益应该等于b/(1-a) 0.019608/0.019608 1和连续系统的直流增益K1一致说明推导没错。这个b/(1-a)等于稳态增益的自检方法非常好用每次做完离散化都建议套一下能拦下相当一部分笔误。同一环节用双线性变换再走一遍做个对照。代入 s (2/T)(1 - z⁻¹)/(1 z⁻¹)整理后 G(z) K(1 z⁻¹) / [(1 2τ/T) (1 - 2τ/T)z⁻¹]。仍然取K1、τ0.05、T0.001则2τ/T 100于是 G(z) (1 z⁻¹)/(101 - 99z⁻¹) 0.009901·(1 z⁻¹)/(1 - 0.980198·z⁻¹)。对应的差分方程是 y[n] 0.980198·y[n-1] 0.009901·x[n] 0.009901·x[n-1]。稳态自检0.009901×2/(1-0.980198) 0.019802/0.019802 1同样没问题。3.2 二阶系统与一般形式的直接读数法一阶系统掌握了真正高频出现的是二阶以及更高阶系统。好在有个直接读数法能省掉大量笔墨任何因果线性离散系统的传递函数都可以写成 H(z) (b0 b1·z⁻¹ b2·z⁻² …) / (1 a1·z⁻¹ a2·z⁻² …)。这个形式一旦整理出来差分方程就是照着系数往下抄y[n] -a1·y[n-1] - a2·y[n-2] … b0·x[n] b1·x[n-1] b2·x[n-2] …。特别注意分母系数的符号和位置这是最容易出错的地方。传递函数分母写成 1 a1·z⁻¹ 时差分方程里对应的是 -a1·y[n-1]如果分母写成 1 - a1·z⁻¹那差分方程里就是 a1·y[n-1]。两种写法都没错但一定要前后统一混着用必然翻车。举一个实际例子某二阶低通滤波器离散化后得到 H(z) (0.0201 0.0402·z⁻¹ 0.0201·z⁻²)/(1 - 1.5610·z⁻¹ 0.6414·z⁻²)照抄即可得到 y[n] 1.5610·y[n-1] - 0.6414·y[n-2] 0.0201·x[n] 0.0402·x[n-1] 0.0201·x[n-2]。整个过程的机械程度很高只要前面传递函数整理正确这一步几乎不可能出错。3.3 双线性变换的预畸变怎么算双线性变换虽然好用但它有个必须正视的副作用频率轴被非线性压缩模拟频率ω_a和数字频率ω_d之间的关系是 ω_a (2/T)·tan(ω_d·T/2)。这意味着如果你在连续域把截止频率设在某个值直接做双线性变换后数字滤波器的实际截止频率会和预期对不上频率越高偏差越大。解决办法是预畸变先根据想要的数字频率反推一个模拟频率用这个修正后的频率去设计连续滤波器再做双线性变换正好把偏差抵消掉。举个数采样周期T0.001s想要数字截止频率对应50Hz即ω_d 2π·50 314.16 rad/s。算ω_d·T/2 314.16×0.001/2 0.15708tan(0.15708) ≈ 0.15844于是ω_a (2/0.001)×0.15844 ≈ 316.88 rad/s换算成频率约50.43Hz。可以看到修正前后只差0.43Hz采样率越高偏差越小。所以如果采样率是信号频率的几十倍以上预畸变带来的修正很小可以偷懒不做但只要截止频率接近采样率的十分之一甚至更高就必须老老实实预畸变否则滤波器实际特性和设计目标会明显偏离。3.4 采样周期T的取值与系数数量级采样周期T是整条链路里最需要凭经验拍板的参数它同时影响精度、稳定性和实时性取大了发散取小了浪费算力。一个常用的经验法则是T取系统时间常数的1/10到1/20。还是拿τ0.05s的系统说T在2.5ms到5ms之间比较合适。另一个更严格的准则是按闭环带宽或截止频率ω_c来判断一般要求采样角频率ω_s 2π/T至少是ω_c的20倍以上保险起见取30到50倍。T的选择还会直接影响系数的数量级这一点在定点数实现时尤其关键。看前面一阶系统的后向差分结果a τ/(τT)T越小a越接近1b KT/(τT)就越小。当T0.001s、τ0.05s时a≈0.98、b≈0.0196两者相差约50倍如果T再缩到0.0001sb就只剩0.002左右了。用16位定点数时这么小的b很容易被截断成零导致控制失效。所以定点实现往往要把整个差分方程两边同时放大一个系数做标幺化处理这个技巧在第5节会具体展开。4. 位置式PID的离散化热词背后到底在说什么4.1 位置式PID的Z变换推导位置式PID用差分方程这个说法最近被反复提到本质就是要把连续PID的积分和微分换成离散形式。连续PID的时域表达式是 u(t) Kp·e(t) Ki·∫e(t)dt Kd·de(t)/dt这里Ki、Kd分别对应积分、微分系数有时写成Kp/Ti和Kp·Td。对三项分别处理比例项不变积分项用累加代替积分即∫e dt ≈ T·Σe[i]对应z域是T/(1-z⁻¹)微分项用后向差分代替即de/dt ≈ (e[n]-e[n-1])/T对应z域是(1-z⁻¹)/T。把三项合起来得到位置式PID的z域表达式U(z) Kp·E(z) Ki·T·E(z)/(1-z⁻¹) (Kd/T)·(1-z⁻¹)·E(z)。为了把它统一成一个差分方程两边同乘(1-z⁻¹)整理就得到时域递推式。不过位置式PID在时域里其实更好理解直接写成 u[n] Kp·e[n] Ki·T·Σ(i0..n)e[i] (Kd/T)(e[n]-e[n-1])右边那个求和项就是积分累加每拍加上一个Ki·T·e[n]即可。这种写法计算量小但有个致命弱点积分项是一个不断累加的全局量一旦被限幅或者系统重启整个积分历史都要跟着处理工程上容易出问题。4.2 增量式PID是怎么从位置式推出来的为了解决位置式的积分管理难题工程里更常用的是增量式PID它其实是从位置式差分方程作差推出来的。把第n拍和第n-1拍的位置式输出相减得到 Δu[n] u[n] - u[n-1]代入展开积分项相减只剩 Ki·T·e[n]比例项变成 Kp·(e[n]-e[n-1])微分项变成 (Kd/T)·(e[n]-2e[n-1]e[n-2])。把同类项按误差归并得到增量式的标准形式Δu[n] A·e[n] - B·e[n-1] C·e[n-2]其中 A Kp(1 T/Ti Td/T)B Kp(1 2Td/T)C Kp·Td/T这里用的是Ti、Td形式。这三个系数一旦算好每拍只需要保存最近两个误差值计算量极小而且天然不需要单独维护一个大积分累加器。唯一要额外注意的是增量式输出的是这一步要加多少真正的控制量要自己累加u[n] u[n-1] Δu[n]所以那个u[n-1]本身要作为状态保存。4.3 两套C代码实现与积分饱和处理先把位置式PID的C代码写出来重点看抗积分饱和是怎么处理的。typedef struct { float Kp, Ki, Kd, T; float integral; // 积分累加器 float e_prev; // 上一次误差 float out_max, out_min; } PID_Pos_t; float pid_position(PID_Pos_t *p, float sp, float fb) { float e sp - fb; p-integral p-Ki * p-T * e; // 积分累加 float d p-Kd / p-T * (e - p-e_prev); // 微分 float u p-Kp * e p-integral d; // 抗积分饱和输出越限时把这一步的积分贡献回退 if (u p-out_max) { u p-out_max; p-integral - p-Ki * p-T * e; } else if (u p-out_min) { u p-out_min; p-integral - p-Ki * p-T * e; } p-e_prev e; return u; }这段代码里最关键的是输出越限时把这一步的积分贡献回退那两行。如果不做这个回退积分项会一直往上累加等到误差反向时控制器要花很长时间才能把积分器拉回来表现出来就是超调大、恢复慢这就是典型的积分饱和。回退的做法相当于告诉积分器这一步的累加没起作用别记了。再看增量式PIDtypedef struct { float A, B, C; float e1, e2; // e[n-1], e[n-2] float out; // 累积输出 float out_max, out_min; } PID_Inc_t; float pid_incremental(PID_Inc_t *p, float sp, float fb) { float e sp - fb; float du p-A * e - p-B * p-e1 p-C * p-e2; p-e2 p-e1; p-e1 e; p-out du; if (p-out p-out_max) p-out p-out_max; else if (p-out p-out_min) p-out p-out_min; return p-out; }增量式代码里积分管理被藏进了累加输出p-out里每拍只做三次乘加逻辑非常干净。切换手动/自动模式或者限幅时也不会像位置式那样留下积分残渣。它的代价是少了积分项的显式表达一旦系统需要做无扰切换bumpless transfer得额外记录输出值。两套代码里我都把输出限幅加上了这个看似不起眼的动作实际上是防止执行机构被顶到极限的关键。5. 实操验证用Python把推导结果验一遍5.1 用双线性变换自动生成系数手工推完系数最好别急着写进固件先用脚本自动生成一遍做交叉验证能省下大量在板子上试错的时间。Python的scipy里就有现成的连续转离散工具一行调用就把系数吐出来了。from scipy import signal import numpy as np # 连续一阶系统 K/(tau*s 1) K, tau, T 1.0, 0.05, 0.001 num_c [K] den_c [tau, 1.0] # 双线性离散化 num_d, den_d, dt signal.cont2discrete((num_c, den_c), T, methodbilinear) num_d num_d[0] # 去掉多余的维度 print(离散分子系数:, num_d) print(离散分母系数:, den_d) # 手动推导的对照值0.009901*(1 z^-1)/(1 - 0.980198*z^-1) print(手动对照分子:, [0.009901, 0.009901]) print(手动对照分母:, [1.0, -0.980198])跑完会发现脚本输出和手动推导完全一致这等于给推导上了双保险。把method换成backward_diff或bilinear就能对比不同离散化方法的结果差异非常直观。我个人的习惯是任何要写进固件的滤波或控制系数都必须先过一遍这个脚本参数对不上就往回查推导绝不靠上板子试试看来验证。5.2 时域仿真对比连续与离散响应系数验证完再跑一遍时域仿真对比连续系统和离散系统的阶跃响应看两者是否贴合。这一步能直观暴露采样周期是否合适、有没有混叠、稳态增益对不对。import numpy as np import matplotlib.pyplot as plt from scipy import signal K, tau, T 1.0, 0.05, 0.001 sys_c signal.TransferFunction([K], [tau, 1.0]) t_c, y_c signal.step(sys_c, Tnp.arange(0, 0.3, 0.0001)) num_d, den_d, _ signal.cont2discrete(([K], [tau, 1.0]), T, methodbilinear) sys_d signal.TransferFunction(num_d[0], den_d, dtT) t_d, y_d signal.dstep(sys_d, tnp.arange(0, 0.3, T)) plt.plot(t_c, y_c, b-, label连续) plt.step(t_d, np.squeeze(y_d), r--, wherepost, label离散) plt.legend(); plt.grid(True); plt.show()正常结果应该是两条曲线几乎重合离散曲线在采样点上顶着连续曲线走。如果离散曲线明显滞后或者振荡先怀疑采样周期太大如果稳态值不重合回头查系数归一化如果振荡发散大概率是离散化方法选错或系数符号写反了。这个对比图是我每次调完参数必看的一张图比盯着数字看高效得多。5.3 定点数实现时的系数放大与舍入前面提过一阶系统里a接近0.98、b只有0.0196低温定点实现时b很容易被截断成0这是定点移植最常见的翻车点之一。解决办法是给整个差分方程做标幺化把输入x先放大再进方程或者把系数统一乘以一个放大倍数最后输出的时候再缩回去。举例说明。用Q15格式取值范围[-1, 1)分辨率2⁻¹⁵ ≈ 3.05e-5实现 y[n] 0.980392·y[n-1] 0.019608·x[n]。直接把0.019608转成Q15是642转换本身没问题但如果输入x的幅度也接近最小值乘出来的结果会低于一个LSB累计之后就损失掉了。更稳妥的做法是先把系数放大2⁴16倍再用Q15表示即让整个系统以输入放大16倍的方式运行输出时右移4位。这样b的有效位数提高了4位小信号精度明显改善。当然放大之后要留意中间结果的位宽会不会溢出Q15乘Q15得Q30累加前最好先做饱和保护。6. 常见问题与排查速查表6.1 输出发散、超调大、稳态有差怎么排查实际调试里遇到的症状就那么几类但原因可能藏在推导、系数、实现的不同环节。我习惯按先看现象再缩范围最后定位的顺序来查。发散一般先怀疑采样周期过大导致离散模型失稳或者系数符号搞反也有可能是双线性变换的公式用错了幂次方向。稳态有差先看稳态增益b/(1-a)是否等于连续系统的直流增益这一步能立刻锁定是否归一化出错。高频振荡多半和微分项有关位置式PID里微分系数是Kd/TT取得很小时这个值会被放得很大只要有一点噪声就会被放大成剧烈的抖动这时候要么加低通滤波要么适当增大T。6.2 常见问题速查表把上面这些症状和对应的排查动作整理成一张表遇到问题直接对号入座比盲猜快得多。症状可能原因排查动作输出单调发散采样周期过大、系数符号反检查T是否超过时间常数核对分母系数符号高频剧烈振荡微分项增益过大、噪声混叠增大T、给微分加一阶低通稳态有静差稳态增益不为1、积分未起作用算b/(1-a)是否等于直流增益阶跃响应有过冲双线性零点、预畸变缺失试后向差分、检查预畸变复位后跳变历史状态被清零、初值项未处理保留历史数组作为状态定点实现失效小系数被截断成0整体放大后再定点化这张表我基本是贴在工位上的遇到新问题就往里加一行久了就成了一本自己的故障手册。要强调的是表里的可能原因永远只是起点真正定位还是要结合数据用仿真脚本复现问题而不是凭感觉改参数。6.3 几个我踩过的坑第一个坑是z⁻¹的方向搞反。早期我把 z⁻¹ 当成了未来值写出来的差分方程用了 x[n1]代码里当然找不到未来的数据直接编译逻辑就错了。其实 z⁻¹ 永远是延迟往后看绝不会往前看这一点想清楚就不会再错。第二个坑是混淆采样周期T和角标T。写公式时T既是采样周期又被当成时间常数符号自己把自己绕晕后来我统一把采样周期写成Ts把时间常数写成τ再没出过这类低级错误。第三个坑最隐蔽把离散化之后得到的差分方程直接原样丢进中断但中断周期和推导时用的T不一致。比如推导用T1ms代码定时器实际是2ms触发一次整个系统的动态特性直接翻倍。这个坑排查了很久才找到后来我在代码里加了一行断言把实际中断周期和设计值打印出来对比问题一眼就现形了。第四个坑是关于位置式PID从手动切自动的瞬间。位置式输出的积分项如果不在切换时做一次同步切入的瞬间输出会跳变执行机构跟着抖一下。解决办法是在切换前把积分器预置成当前输出减去比例和微分贡献的值这样切换前后输出连续也就是常说的无扰切换。这个细节很多入门资料不讲但实际项目里如果忽略轻则设备抖动重则触发保护停机。最后再分享一个我常用的检查习惯每次把Z变换方程整理成差分方程后先令输入为一个常数看看输出能否收敛到连续的直流增益再令输入为一个阶跃看前几拍数值是否符合直觉。这两步加起来不到五分钟却往往能在写代码之前就拦下错误。差分方程看着简单真正让它在硬件上稳定跑起来靠的从来不是运气而是推导时的细心和排查时的耐心。