做数字信号处理、嵌入式控制这一类活儿的人迟早会跟Z变换正面撞上。课本上那个定义式 X(z) Σ x[n]·z^(-n)第一眼看过去就像把一串离散的数硬塞进一个无穷级数里看不出它跟滤波器系数、跟板子上跑的差分代码、跟示波器上那条发散或者振铃的曲线有什么关系。可只要你认真设计过一支数字滤波器或者用双线性变换把一个连续域控制器搬到单片机上就会明白绕不过去为什么改一个系数输出就自激了为什么相位延迟总是对不上为什么仿真完美上板就废答案几乎都藏在Z变换里。这篇内容我按自己这些年反复翻书、反复踩坑的顺序把Z变换从头梳理一遍从定义和收敛域到性质、系统函数、极零点、逆变换再到落地成可运行的代码和排查技巧。适合刚学完信号与系统还没把公式和工程连起来的学生也适合转做数字信号处理、需要快速把这套工具捡回来的工程师。1. 为什么离散系统里Z变换是绕不过去的门槛我真正意识到Z变换的分量是在一次调二阶IIR滤波器的时候。当时连续域那套东西我闭着眼都能写极点实部为负就稳定阻尼比一算超调量大概多少心里有数。结果换成离散域之后同样的思路全部失效稳定判据从左半平面变成了单位圆内频响从虚轴变成了圆上走一圈。那一刻我才明白Z变换不是拉普拉斯变换换了个字母它是一套完整的、自洽的语言专门用来描述离散时间系统。它最大的价值是把差分方程这种看起来只能一步步递推的东西变成了代数方程。差分方程之所以烦人是因为 y[n] 依赖 y[n-1]、y[n-2]求解必须从初值开始一格一格往前推想看清整体行为非常困难。而Z变换一出手延迟操作变成乘 z^(-1)卷积变成乘法整个递推关系直接变成两个多项式的比值极零点、稳定性、频率响应全部可以一次性看清。还有一点少有人一开始就讲透离散系统里频率这个概念本身就是有边界的。连续域的频率可以从0一直延伸到无穷而采样之后频率被折叠到 [0, fs/2] 这个区间Z变换的 z e^(jω) 正好把单位圆映射成这条频率轴ω 从 0 走到 π 就对应直流到奈奎斯特频率。这个几何图像是后面理解一切数字滤波器行为的根。1.1 一个让我重新翻书的真实场景给一个电流采样环节做数字低通我在连续域用二阶巴特沃斯设计好截止频率定在采样率的十分之一然后用双线性变换离散化浮点仿真曲线又平又顺一点问题没有。烧进定点DSP之后输出开始缓慢地自激幅度一点点涨上去最后顶到满量程。排查了两天问题出在系数上。我的极点在离散域里模长是 0.94非常靠近单位圆。理论上它还是稳定的但定点量化带来的系数误差让实际计算的极点稍微往外挪了一点点模长变成 1.0002于是就发散了。当时我对Z变换的理解只到会算没有到知道极点位置对量化有多敏感所以完全没想到这一层。这件事之后我养成了两个习惯。第一任何IIR设计完第一件事就是把极零点算出来看最大极点模离单位圆还有多远。第二不管数字多好看都要用定点系数再跑一遍对比。这两条后来救了我很多次。1.2 拉普拉斯、傅里叶和Z变换到底是什么关系很多人把这三个变换当成三门独立的课其实它们是一家人只是站在不同视角看同一个问题。傅里叶变换问的是这个信号里各个频率成分各占多少拉普拉斯变换在傅里叶的基础上加了一个衰减因子把不收敛的信号也拉进可分析的范畴而Z变换就是拉普拉斯变换在离散域的对应物。变换适用对象变换变量收敛区域主要用途傅里叶变换连续信号jω虚轴若不收敛取广义频谱分析拉普拉斯变换连续系统s σ jωs平面某条带或半平面连续系统求解与稳定性DTFT离散序列e^(jω)单位圆若非绝对可和取广义离散频谱分析Z变换离散序列z r·e^(jω)z平面某个环状区域差分方程、稳定性、极零点分析把它们串起来的关键是一句映射关系z e^(sT)T 是采样周期。s 平面的左半平面实部小于0连续系统稳定区正好映射到 z 平面的单位圆内s 平面的虚轴傅里叶变换那条线映射到单位圆本身s 平面实部为正的不稳定区映射到单位圆外。所以离散系统稳定性看极点是否在单位圆内这句话本质上是连续系统稳定性看极点是否在左半平面经过指数映射之后的结果。理解了这层映射很多东西就不用死记了。比如为什么离散系统的频率响应是周期性的因为 z e^(jω) 里 ω 加上 2π 之后 z 完全不变所以频响天然以 2π 为周期。再比如为什么采样会产生混叠因为映射 z e^(sT) 不是单射很多个不同的 s 会映到同一个 z 上。提示把 z e^(sT) 这个关系刻在脑子里后面双线性变换、预畸变、频率折叠这些内容都能自己推出来不需要背结论。1.3 Z变换真正解决的三个核心问题第一是求解。把差分方程两边同时做Z变换利用位移性质把 y[n-k] 变成 z^(-k)Y(z)整个递推关系直接退化成一个关于 Y(z) 的代数方程移项就得到输出。这跟用拉普拉斯变换解微分方程是完全一样的套路只是对象从微分换成了差分。第二是稳定性判定。一个离散系统的稳定性完全由它传递函数的极点位置决定极点全部在单位圆内就稳定有一个在外面就发散在圆上就是临界振荡。这个判据简单到可以口算前提是你能把极零点算出来。第三是频率响应的快速估算。系统的频率响应就是 H(z) 在单位圆上取值而 |H(e^(jω))| 可以写成所有零点到该点的距离之积除以所有极点到该点的距离之积。不需要真的去算每一点的数值光看极零点图就能估出通带、阻带和峰值的形状这在设计阶段极其省时间。2. 定义与收敛域把公式背后的约束讲清楚定义本身没什么好纠结的真正容易出错的是收敛域。我见过太多人只写 X(z) 的表达式不写收敛域然后两个人在同一个式子上得出完全相反的结论一个说是因果稳定系统一个说是反因果发散系统。根源在于Z变换是一个无穷级数它并不对所有的 z 都收敛。z 是复变量|z| 大一点或者小一点级数的敛散性会变所以任何一个Z变换都必须带上收敛域才有完整意义。表达式相同的两个Z变换如果收敛域不同对应的时域序列可以完全不同。2.1 双边定义和单边定义工程上到底用哪个双边Z变换的定义是X(z) Σ (n 从 -∞ 到 ∞) x[n]·z^(-n)单边Z变换的定义是X₁(z) Σ (n 从 0 到 ∞) x[n]·z^(-n)两者的区别只在求和起点。对于 n 0 时全部为零的因果序列两者完全等价。那为什么还要单独定义单边变换因为解带初始条件的差分方程时单边变换非常方便。举个直观的例子对 x[n-1] 做双边变换得到的是 z^(-1)X(z)干脆利落但对因果序列做单边变换时得到的是 z^(-1)X(z) x[-1]那个 x[-1] 就是初始条件项它会自动出现在方程里。很多教材讲解初始条件响应时用的就是单边变换就是为了让初值自然地冒出来不用额外处理。我个人的习惯是分析稳态特性和滤波器设计用双边求解带初值的瞬态响应用单边。两者不要混着用混用最容易漏掉初值项。2.2 收敛域为什么是一环而不是一片把 z 写成极坐标 z r·e^(jω)代入定义式|X(z)| ≤ Σ |x[n]|·r^(-n)你会发现右边这个级数的敛散性只跟 r也就是 |z|有关跟幅角 ω 一点关系都没有。所以收敛域天然就是一个以原点为中心的环状区域可以是 |z| r₁可以是 |z| r₂也可以是 r₁ |z| r₂。这个结论非常重要它意味着收敛域的边界必然由极点位置决定——因为让级数发散的地方只能来自表达式中的极点。具体规则可以总结成一张表序列类型时域特征收敛域形状边界由谁决定有限长序列只在有限个 n 上非零整个z平面可能排除 z0 或 z∞无极点右边序列n n₀ 时为零|z| r₁圆外区域最外侧的极点模长左边序列n n₀ 时为零|z| r₂圆内区域最内侧的极点模长双边序列两边都无限延伸r₁ |z| r₂环形区域内外两侧极点最实用的一条经验收敛域内绝对不能包含极点。所以你在画收敛域的时候边界一定是被某个极点撑出来的不会凭空出现。2.3 几个必须记住的变换对以及它们是怎么来的背表不如会推。指数序列 a^n·u[n] 的Z变换推导只用了等比数列求和X(z) Σ (n≥0) a^n·z^(-n) Σ (n≥0) (a·z^(-1))^n这是一个公比为 a·z^(-1) 的等比级数收敛条件是 |a·z^(-1)| 1也就是 |z| |a|。在收敛域内求和结果是X(z) 1 / (1 - a·z^(-1)) z / (z - a)注意那个收敛条件 |z| |a|它和极点位置 z a 完全吻合——极点被收敛域排除在外了这就是前面说的规律。把常见的变换对整理成表方便随时查时域序列 x[n]Z变换 X(z)收敛域δ[n]1整个z平面u[n]1 / (1 - z^(-1))|z| 1a^n·u[n]1 / (1 - a·z^(-1))|z| |a|-a^n·u[-n-1]1 / (1 - a·z^(-1))|z| |a|n·a^n·u[n]a·z^(-1) / (1 - a·z^(-1))²|z| |a|cos(ω₀n)·u[n](1 - z^(-1)cosω₀) / (1 - 2z^(-1)cosω₀ z^(-2))|z| 1sin(ω₀n)·u[n]z^(-1)sinω₀ / (1 - 2z^(-1)cosω₀ z^(-2))|z| 1这张表里第三行和第四行要特别小心它们的表达式一模一样只有收敛域不同。前者对应右边序列因果后者对应左边序列反因果。这就是为什么我说不写收敛域的Z变换是不完整的。注意余弦和正弦那一行的分母结构 1 - 2z^(-1)cosω₀ z^(-2) 是一对共轭极点极点在 z cosω₀ ± j·sinω₀ e^(±jω₀)模长正好为1。这解释了为什么理想正弦振荡器是临界稳定系统。3. 核心性质工程直觉比公式推导更重要性质这一块考试的时候大家都会背但真正落到工程上能用起来的其实就那么几条。我的建议是先理解每一条性质在物理上意味着什么再去记公式这样遇到没见过的问题也能自己推。3.1 位移性质延迟一拍就是乘以 z 的负一次方位移性质是Z变换里用得最多的一条没有之一。它的内容很简单若 x[n] 的Z变换是 X(z)则 x[n - n₀] 的Z变换是 z^(-n₀)·X(z)。单边变换里要补初值项比如 x[n-1] 对应 z^(-1)X(z) x[-1]。这条性质之所以重要是因为它把时间上延迟这个操作直接翻译成了乘一个 z^(-1)。数字系统里的每一个延迟单元、每一级流水线、每一个上一拍的采样值在Z域里都是一个 z^(-1)。滤波器实现结构图里的延迟框画的就是它。由此可以推出一个很实用的结论纯延迟环节 H(z) z^(-n₀)它的频率响应是 e^(-jωn₀)幅值恒为1相位是 -ωn₀也就是线性的相位延迟。这就是为什么线性相位FIR滤波器能靠对称结构实现——它的群延迟是常数不会对信号造成相位失真。3.2 卷积定理为什么串联系统可以直接相乘卷积定理说的是若 y[n] x[n] * h[n]则 Y(z) X(z)·H(z)。这条性质的工程价值极高。时域里的卷积是一个非常耗时的操作尤其是长序列但换到Z域就变成了简单的多项式相乘。更关键的是两个系统级联时整体传递函数就是各自传递函数的乘积这让复杂系统的分析变得极其简单——你可以把一长串处理模块拆成若干段逐个分析极零点最后再合并。在滤波器设计里这正是级联型cascade结构的理论基础。一个八阶滤波器直接实现时系数量化误差会严重影响极点位置但如果拆成四个二阶节biquad级联每一节的极点是独立配置的量化误差只影响局部整体鲁棒性大幅提升。我前面说的那次定点自激最后的解决方案就是改成二级级联。3.3 初值定理和终值定理稳态误差怎么算这两条定理在控制系统里特别实用因为很多时候你只关心最终稳定在多少不需要把整个响应算出来。初值定理x[0] lim (z→∞) X(z)。前提是 X(z) 表达式是 z 的真分式或者严格真分式。它的用处是快速验证逆变换结果对不对——如果算出来的 x[0] 跟直接看序列第一项对不上说明中间某一步出错了。终值定理lim (n→∞) x[n] lim (z→1) (1 - z^(-1))·X(z)。它的使用有严格前提X(z) 的极点必须全部在单位圆内唯一的例外是允许在 z 1 处有一个单阶极点对应一个恒定的直流分量。举个实际例子。系统输入为单位阶跃 U(z) 1/(1 - z^(-1))误差传递函数的稳态值想要求出来直接代入终值定理就行lim (z→1) (1 - z^(-1))·X(z)如果 X(z) 里含有 1/(1 - z^(-1)) 这个因子那么 (1 - z^(-1)) 正好把它约掉剩下的就是稳态值。这就是控制系统里计算稳态误差的标准方法。注意终值定理最容易被误用。如果系统有一个极点在 z -1 上对应 ω π 的振荡或者有一对共轭极点在单位圆上终值定理算出来的数字毫无意义因为系统根本不收敛。用之前先检查极点。3.4 性质速查表性质名称时域关系Z域关系工程用途线性a·x₁[n] b·x₂[n]a·X₁(z) b·X₂(z)叠加原理分解复杂信号位移x[n - n₀]z^(-n₀)·X(z)延迟单元、流水线建模z域尺度a^n·x[n]X(z/a)加窗、指数加权时域卷积x[n] * h[n]X(z)·H(z)级联系统、滤波器实现时域乘 nn·x[n]-z·dX(z)/dz斜坡响应分析初值x[0]lim (z→∞) X(z)结果验证终值lim (n→∞) x[n]lim (z→1) (1-z^(-1))X(z)稳态误差计算累加Σ (k≤n) x[k]X(z)/(1 - z^(-1))积分器建模最后一行累加其实是个隐藏福利数字积分器就是累加器它的Z变换就是在原信号上乘一个 1/(1 - z^(-1))也就是在 z 1 处加了一个极点。这就是为什么PI控制器里那个I会在原点附近制造一个极点带来无穷大的直流增益和零稳态误差同时也带来相位滞后。4. 从差分方程到系统函数极零点、稳定性与频率响应前面的都是工具这一节才是真正把工具组装起来用的地方。一个离散系统通常用差分方程描述而Z变换的任务就是把这个差分方程变成传递函数再从传递函数里读出所有我们关心的信息。4.1 传递函数 H(z) 是怎么导出来的一般的线性常系数差分方程写成Σ (k0 到 N) aₖ·y[n-k] Σ (k0 到 M) bₖ·x[n-k]注意 a₀ 通常归一化为1。两边同时做双边Z变换利用位移性质每一项 y[n-k] 变成 z^(-k)Y(z)x[n-k] 变成 z^(-k)X(z)Y(z)·Σ aₖz^(-k) X(z)·Σ bₖz^(-k)于是H(z) Y(z)/X(z) (Σ bₖz^(-k)) / (Σ aₖz^(-k))分子多项式的根是零点分母多项式的根是极点。整个系统的行为就由这两组根完全决定。这也顺便说明了一件事为什么数字滤波器的系数不能随便取因为系数的微小变化会直接改变根的位置进而影响稳定性。如果要写成 z 的正幂次形式两边同乘 z^NH(z) z^(N-M)·(Σ bₖz^(M-k)) / (Σ aₖz^(N-k))这时候要注意 z 0 处可能出现的极点或零点别在数根的时候漏掉。4.2 极零点图的正确读法极零点图就是一个复平面横轴是实部纵轴是虚部单位圆画在中间极点和零点分别用不同的记号标出来。图本身很简单难的是怎么读。我的经验是按频率沿着单位圆走一圈来读。单位圆上角度为 ω 的那个点对应的是数字频率 ω 处的响应。零点离单位圆越近它对应的频率附近幅频响应就越低形成凹口极点离单位圆越近对应频率附近幅频响应就越高形成尖峰。如果零点正好落在单位圆上那个频率的增益严格为零这就是我们常说的陷波。还有一条规律靠近原点的极零点对幅频形状几乎没影响只影响整体的幅度缩放和相位。真正决定通带阻带形状的永远是那些靠近单位圆的家伙。所以设计滤波器时重点盯住靠近单位圆的那几个极零点就够了。极零点位置对幅频的影响对相位的影响常见场景极点靠近单位圆内侧该频率处出现高增益尖峰该频率附近相位急剧变化谐振器、窄带带通零点落在单位圆上该频率增益严格为零相位发生跳变陷波器、去除特定干扰零点在单位圆内该频率附近增益下降相位超前高通、预加重极点在原点纯延迟幅值不变线性相位滞后延迟链极点或零点成共轭对幅频对称无虚部系数——实数系数的必然结果最后一行值得多说一句。因为实际系统的系数都是实数多项式有实系数所以复数根必然以共轭对的形式出现。这意味着任何实系数滤波器的幅频响应都是关于 ω 0 和 ω π 对称的。这不是巧合是数学上的必然。4.3 稳定性判据单位圆为什么是生死线对于一个因果系统冲激响应在 n 0 时为零稳定的充要条件是所有极点都在单位圆内也就是 |pⱼ| 1。这个结论的直觉解释是极点的模长 |p| 对应时域中 a^n 的底数如果 |p| 1那么 a^n 随着 n 增大衰减到零系统冲激响应绝对可和如果 |p| 1冲激响应指数增长系统发散。对于非因果系统判据要放宽成收敛域包含单位圆。因为只有收敛域包含单位圆频率响应才存在系统才是稳定这里指有界输入有界输出稳定的。实际工程中还有一类特殊情况极点在单位圆上。这时候系统是临界稳定或者叫临界振荡冲激响应既不衰减也不发散而是等幅振荡。数字积分器极点 z 1和数字振荡器极点 z e^(±jω₀)都属于这一类。理论上有界输入不一定有界输出但工程上如果输入是直流或者有限时长的信号它是可以用的只要你能接受那个等幅输出。提示判断稳定性时如果你手算的极点是 0.99、0.998 这种数字千万不要觉得小于1就安全。在定点实现里0.998 和 1.002 之间可能只差一个最低有效位。这类系统必须做定点仿真验证而且要预留足够的稳定裕度我个人经验是最大极点模最好别超过 0.98。4.4 z e^(jω) 的几何意义频率响应怎么一眼估出来把 H(z) 写成因式分解形式H(z) K·Π(z - zᵢ) / Π(z - pⱼ)令 z e^(jω)取模|H(e^(jω))| |K|·Π|e^(jω) - zᵢ| / Π|e^(jω) - pⱼ|这里每个 |e^(jω) - zᵢ| 的几何意义是单位圆上的点 e^(jω) 到零点 zᵢ 的直线距离。分母同理是到极点的距离。所以幅频响应就是到所有零点的距离之积除以到所有极点的距离之积再乘一个常数。这个几何解释非常直观。当 ω 走到某个极点附近时分母中的那一项距离变得很小整个式子的值就变得很大形成峰值当 ω 走到某个零点附近时分子中的那一项变小整个式子的值减小形成凹陷如果零点正好在单位圆上距离为零增益严格为零。用这个方法可以快速估出一支滤波器的中心频率和带宽中心频率就在极点的幅角位置带宽和极点离单位圆的远近成反比——越靠近圆峰越尖锐带宽越窄。这个估算在设计阶段非常有用可以让你在动手写代码之前就大致知道结果对不对。5. 逆Z变换的三种求法以及什么时候用哪个正变换好办直接套公式。逆变换才是真正需要技术的地方因为没有一个万能的公式得根据具体情况选方法。我常用的有三种各有各的适用场景下面逐个说清楚。5.1 部分分式展开法最通用也最容易出错部分分式展开的思路和求拉普拉斯逆变换一模一样把复杂的分式拆成若干个简单分式之和每个简单分式对应一个已知的变换对直接查表就行。但这里有一个非常关键的细节坑了无数人做部分分式之前应该先展开 X(z)/z而不是 X(z) 本身。原因是标准的部分分式公式针对的是分子次数低于分母次数的真分式而 X(z) 往往是 z 的有理函数分子分母次数相同直接展开会多出一个常数项最后反变换的时候就会丢掉一个 z 因子结果全错。完整流程是第一步求 X(z)/z第二步对 X(z)/z 做部分分式展开第三步两边同乘 z 得到 X(z)第四步逐项查表反变换。举个具体例子。设X(z) 1 / [(1 - 0.5z^(-1))(1 - 0.25z^(-1))]假设收敛域是 |z| 0.5因果序列。改写一下X(z)/z z / [(z - 0.5)(z - 0.25)] A/(z - 0.5) B/(z - 0.25)求系数A 0.5/(0.5 - 0.25) 2B 0.25/(0.25 - 0.5) -1。于是X(z) 2z/(z - 0.5) - z/(z - 0.25)查表得 x[n] [2·(0.5)^n - (0.25)^n]·u[n]。验证一下 n 0 时2 - 1 1跟原式分子常数项一致说明没错。如果直接在 z^(-1) 变量下做部分分式得到的结论其实是一样的但前提是分子在 z^(-1) 下的阶数低于分母。这两种写法都行我个人更习惯用 z^(-1) 变量因为最终结果直接对应延迟单元写代码的时候系数一一对应不容易抄错。5.2 长除法适合短序列和快速验证长除法就是把 X(z) 按 z^(-1) 的升幂展开成幂级数展开出来的系数序列就是 x[n]。原理很简单因为 Z变换的定义本身就是幂级数。还是用上例X(z) 1 / [(1 - 0.5z^(-1))(1 - 0.25z^(-1))] 1 / (1 - 0.75z^(-1) 0.125z^(-2))分子1除以分母得到1 0.75z^(-1) 0.4375z^(-2) ...所以 x[0] 1x[1] 0.75x[2] 0.4375。用上一节的闭式解验算x[1] 2×0.5 - 0.25 0.75x[2] 2×0.25 - 0.0625 0.4375。对上了。长除法的好处是不需要求根、不需要解方程纯机械操作特别适合在纸上验证前面几步的结果。缺点是当序列很长时这个幂级数不会收敛成简洁形式你只能得到一个数值序列。所以它适合确认前几项对不对不适合求通项。5.3 留数法理论推导和验证用留数法的公式是x[n] Σ Res[X(z)·z^(n-1)]对所有位于收敛域内的极点求和这里的留数按复变函数的规则算。对于单阶极点 p留数是 lim(z→p) (z-p)·X(z)·z^(n-1)。对于高阶极点需要用求导公式。留数法在理论上最完整适合处理复杂情况比如有多重极点或者需要严格证明某个结论。但实际工作里用得不多因为计算量大容易算错。我一般只在两种情况下用一是验证部分分式的结果二是在写论文或者做推导需要严谨表达的时候。5.4 用代码做数值验证手算容易出错最稳的验证方式是直接用代码算一遍。我用 SymPy 做符号计算几行就能确认部分分式展开对不对import sympy as sp z sp.symbols(z) Xz 1 / ((1 - 0.5/z) * (1 - 0.25/z)) print(X(z) , sp.simplify(Xz)) print(X(z)/z , sp.simplify(Xz / z)) print(部分分式展开) print(sp.apart(Xz / z, z))输出里会看到 X(z)/z 2/(z - 0.5) - 1/(z - 0.25)和手算完全一致。如果想拿到前若干项时域系数用长除法展开更直接import sympy as sp z sp.symbols(z) Xz 1 / ((1 - 0.5/z) * (1 - 0.25/z)) series sp.series(Xz, z, sp.oo, 8).removeO() print(series)至于数值验证最快的办法是把这个Z变换对应的差分方程写出来用 lfilter 跑一遍冲激响应看是不是和理论值吻合。这个方法我在下一节会详细展开。注意无论用哪种方法算逆变换一定要用初值定理先验一下 x[0]。如果 x[0] 对不上后面的都不用看了肯定是表达式或者收敛域写错了。这一步只花十秒钟能省掉半小时的排查。6. 动手实操从差分方程到可运行代码的完整链路前面全是理论这一节我拿一个真实的二阶系统走一遍完整流程写差分方程、推导传递函数、求极零点、判断稳定性、算频率响应、写代码验证。整个过程你可以直接抄去改参数用。6.1 选一个系统把传递函数推出来选一个二阶谐振器型的系统差分方程如下y[n] 0.0725·x[n] 0.145·x[n-1] 0.0725·x[n-2] 1.6·y[n-1] - 0.89·y[n-2]移项整理成标准形式y[n] - 1.6·y[n-1] 0.89·y[n-2] 0.0725·x[n] 0.145·x[n-1] 0.0725·x[n-2]对比一下系数a₀ 1a₁ -1.6a₂ 0.89b₀ 0.0725b₁ 0.145b₂ 0.0725。做Z变换得到H(z) (0.0725 0.145z^(-1) 0.0725z^(-2)) / (1 - 1.6z^(-1) 0.89z^(-2))分子的三个系数正好是 [1, 2, 1] 的缩放也就是说分子 0.0725·(1 z^(-1))²所以在 z -1 处有一个二阶零点。除以 z^(-2) 换成正幂次形式后还要注意 z 0 处有两个极点。分母的根解 z² - 1.6z 0.89 0判别式是 2.56 - 3.56 -1所以 z 0.8 ± 0.5j。两个共轭极点。项目数值含义零点z -1二阶在奈奎斯特频率处形成零点高频完全被压制极点 10.8 0.5j谐振点之一极点 20.8 - 0.5j共轭谐振点极点模长√(0.64 0.25) √0.89 ≈ 0.9434小于1因果系统稳定极点幅角atan(0.5/0.8) ≈ 0.5586 rad谐振频率位置系统的直流增益可以直接代入 z 1 算分子是 0.0725×4 0.29分母是 1 - 1.6 0.89 0.29所以直流增益正好是1。这说明我选的分子系数是经过归一化的通带增益为1方便观察。谐振频率换算成实际频率ω 0.5586 rad/sample如果采样率是 48 kHz那么谐振频率是 0.5586/(2π)×48000 ≈ 4268 Hz。6.2 Python 代码逐行验证第一步算极零点并检查稳定性import numpy as np from scipy import signal b np.array([0.0725, 0.145, 0.0725]) a np.array([1.0, -1.6, 0.89]) z, p, k signal.tf2zpk(b, a) print(零点, np.round(z, 4)) print(极点, np.round(p, 4)) print(极点模长, np.round(np.abs(p), 4)) print(最大极点模长, np.max(np.abs(p)))预期输出极点模长 0.9434小于1稳定。第二步跑冲激响应同时用手写递推验证一致n np.arange(0, 80) x np.zeros_like(n, dtypefloat) x[0] 1.0 # 方式一直接用 lfilter y_auto signal.lfilter(b, a, x) # 方式二手写递推逐项验证 y_manual np.zeros_like(x) for i in range(len(x)): s b[0] * x[i] if i 1: s b[1] * x[i-1] - a[1] * y_manual[i-1] if i 2: s b[2] * x[i-2] - a[2] * y_manual[i-2] y_manual[i] s print(两种实现是否一致, np.allclose(y_auto, y_manual)) print(前12点冲激响应, np.round(y_auto[:12], 4))手写递推这一段的写法值得注意因为 a₀ 归一化为1所以我把 y[n] 单独留在左边然后把 -a₁y[n-1] - a₂y[n-2] 移到右边变成加号。这个正负号是最容易写错的地方很多人的代码跑出来结果不对就是这里符号搞反了。第三步算频率响应并找出峰值位置w, h signal.freqz(b, a, worN4096) mag_db 20 * np.log10(np.abs(h) 1e-12) peak_idx np.argmax(mag_db) peak_w w[peak_idx] fs 48000.0 print(峰值数字频率 rad/sample, round(peak_w, 4)) print(峰值实际频率 Hz, round(peak_w / (2*np.pi) * fs, 1)) print(峰值增益 dB, round(mag_db[peak_idx], 2))理论谐振频率是 0.5586 rad/sample代码算出来的峰值位置应该在这个数附近可能会有一点点偏移因为零点也会对峰值位置产生牵引作用。这个偏差本身就是一个值得关注的现象——它说明极点幅角就是谐振频率这个说法只是近似零点存在的时候会偏移。第四步画极零点和频响图import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(12, 5)) # 极零点图 theta np.linspace(0, 2*np.pi, 400) axes[0].plot(np.cos(theta), np.sin(theta), k--, linewidth1) axes[0].plot(z.real, z.imag, o, markersize8, label零点) axes[0].plot(p.real, p.imag, x, markersize10, label极点) axes[0].set_aspect(equal) axes[0].grid(True, alpha0.3) axes[0].legend() axes[0].set_title(极零点分布) # 幅频响应 axes[1].plot(w / np.pi, mag_db) axes[1].grid(True, alpha0.3) axes[1].set_xlabel(归一化频率 (×π rad/sample)) axes[1].set_ylabel(幅度 (dB)) axes[1].set_title(幅频响应) plt.tight_layout() plt.show()图出来之后你应该能看到极点几乎贴着单位圆内壁零点都堆在 z -1 那个位置幅频响应在 0.56 rad/sample 附近有一个明显的谐振峰在奈奎斯特频率处掉得很深。这就是极零点图直接告诉你频响形状的直观演示。6.3 参数选择和实操中踩过的坑这块是文档里不会写的部分都是我自己撞出来的。关于采样率的选择。谐振频率定下来之后采样率不能随便取。如果采样率太低谐振频率靠近奈奎斯特频率极点会挤在单位圆右半边靠近 -1 的地方系数量化的敏感度会急剧上升。我的经验是让目标频率落在采样率的 5% 到 20% 之间这个区间里设计裕度和实现难度都比较合理。关于极点到单位圆的距离和量化精度。极点模长 0.9434距离单位圆还有 0.0566 的裕度看起来很安全。但如果用 Q15 定点表示系数的量化步长是 2^(-15) ≈ 3.05e-5。对于二阶系统极点位置对系数的敏感度大约和极点模长成正比系数误差 3e-5 会导致极点模长变化量级在 1e-5 到 1e-4 之间。0.0566 的裕度是够的但如果极点模长做到了 0.998裕度只有 0.002量化误差就可能把它推出去。所以判断能不能用定点不能只看模长小于1要看裕度够不够大。关于高阶级联结构。我前面反复强调过四阶以上的滤波器不要用直接型Direct Form实现一定要拆成二阶节级联。原因是高阶多项式的根对系数极其敏感一个八阶直接型滤波器系数量化之后极点可能完全跑偏。拆成二阶节之后每一节的极点只受本节的四个系数影响误差被隔离在局部。关于双线性变换的预畸变。如果你是先用连续域指标设计好再离散化双线性变换会引入频率畸变映射关系是 ω_digital 2·arctan(ω_analog·T/2)。也就是说数字域的截止频率和模拟域的截止频率不是线性对应的。正确的做法是先按目标数字频率做预畸变把模拟设计频率往前推让映射之后正好落在你要的位置。这一步漏掉截止频率会偏而且频率越高偏得越厉害。提示预畸变的完整流程是先算 ω_a (2/T)·tan(ω_d·T/2)用 ω_a 去做模拟滤波器设计再离散化。如果需要我可以单独开一篇把这个流程拆细因为里面的坑比这里写的还多。7. 常见问题与排查技巧实录最后这一节把我这些年遇到过的典型问题整理成速查表配上线上的排查思路。列出的都是真实遇到过的不是编的。7.1 极点在单位圆上到底算不算稳定这是被问得最多的问题。答案是严格意义上的有界输入有界输出稳定要求极点严格在单位圆内。极点在单位圆上是临界稳定输入有界时输出可能无界——典型例子就是累加器输入单位阶跃输出是斜坡显然是发散的。但工程上并不是临界稳定就完全不能用。数字积分器就是极点在 z 1 处的系统几乎所有控制器都在用。区别在于输入信号的特性决定了输出会不会累积发散。纯直流输入进积分器输出一直涨但如果你有饱和限幅就是可控的而振荡器极点在 z e^(±jω)输入持续的能量才会让它涨输入一停它就保持等幅。所以实际判断的时候我会分三步。先算极点模长看是否有超过1的有就直接判定不稳定。再看是否有等于1的有的话确认这是设计故意为之积分器或者振荡器还是意外。最后看那些小于1但非常接近1的比如0.999这些是准临界需要重点做定点仿真验证。极点模长稳定性分类时域表现工程处理方式明显小于1 0.95稳定裕度充足冲激响应快速衰减可直接实现0.95 到 0.99稳定裕度偏小衰减较慢振铃较长建议用级联结构做定点验证大于 0.99 小于1稳定但极度敏感衰减很慢接近等幅慎重使用必须浮点或高精度定点等于1临界稳定等幅或缓慢漂移仅用于积分器需配限幅大于1不稳定指数发散设计错误必须重新设计7.2 收敛域漏写会带来什么后果前面提过同一个表达式对应不同的收敛域就对应完全不同的序列。这里给一个具体例子也是一个经典陷阱。X(z) 1/(1 - 0.5z^(-1))这个表达式对应三个可能收敛域 |z| 0.5对应右边序列 x[n] 0.5^n·u[n]因果系统稳定。收敛域 |z| 0.5对应左边序列 x[n] -0.5^n·u[-n-1]反因果不稳定。收敛域空不可能或者包含单位圆的其他形式需要具体判断。如果你只写了那个分式然后去判断稳定性得到极点在 0.5在单位圆内所以稳定那只是碰巧对了其中一个分支。真正严谨的结论是只有当收敛域是 |z| 0.5 时这个因果系统才稳定。我在看别人的设计文档时只要看到只写传递函数不写收敛域就会多留个心眼因为后面很可能藏着因果性或者稳定性的误判。7.3 问题排查速查表现象可能原因排查动作冲激响应指数发散有极点模长大于1算极零点检查是否量化导致输出自激但浮点仿真正常定点系数量化让极点跑出单位圆对比浮点和定点系数下的极点位置频率响应峰值位置和理论对不上零点对峰值的牵引或双线性变换未预畸变用代码算实际峰值位置检查预畸变稳态值算错终值定理结果离谱系统有单位圆上极点终值定理前提不满足先检查极点位置再决定能不能用逆变换结果比理论少一个 z 因子部分分式时用了 X(z) 而不是 X(z)/z重新用 X(z)/z 展开手写递推和 lfilter 结果不一致差分方程移项时符号写反仔细核对 a 系数的正负号高阶滤波器定点后完全失效用了直接型高阶多项式根对系数敏感拆成二阶节级联截止频率整体偏移双线性变换频率畸变未预畸变用 ω_a (2/T)tan(ω_d·T/2) 重算幅频响应不对称系数量化破坏了共轭对称性检查是否强制共轭配对或者用了非实数系数这张表基本覆盖了我在实际项目里遇到的八成问题。剩下两成通常是多个原因叠加那就得一步步来先确认极点位置和稳定性再确认频率响应最后才是量化误差。排查顺序永远是从大到小先看结构性错误再看精度问题。我个人的习惯是任何数字滤波器设计完之后先跑一个自动检查脚本算极零点、算最大极点模长、算直流增益、跑一遍冲激响应看有没有发散。这四个检查加起来不到二十行代码但能在设计阶段就挡掉绝大部分问题比烧到板子上之后对着示波器发呆划算得多。极点模长这个数字我现在基本是设计完第一眼就要看的比任何曲线都直观。