做数字信号处理、嵌入式控制这一类活儿的人，迟早会跟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变换，就是拉普拉斯变换在离散域的对应物。

把它们串起来的关键是一句映射关系：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₂。这个结论非常重要，它意味着收敛域的边界必然由极点位置决定——因为让级数发散的地方只能来自表达式中的极点。

具体规则可以总结成一张表：

最实用的一条经验：收敛域内绝对不能包含极点。所以你在画收敛域的时候，边界一定是被某个极点"撑"出来的，不会凭空出现。

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 完全吻合——极点被收敛域排除在外了，这就是前面说的规律。

把常见的变换对整理成表，方便随时查：

这张表里第三行和第四行要特别小心：它们的表达式一模一样，只有收敛域不同。前者对应右边序列（因果），后者对应左边序列（反因果）。这就是为什么我说不写收敛域的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变换就是在原信号上乘一个 1/(1 - z^(-1))，也就是在 z = 1 处加了一个极点。这就是为什么PI控制器里那个I，会在原点附近制造一个极点，带来无穷大的直流增益和零稳态误差，同时也带来相位滞后。

4. 从差分方程到系统函数：极零点、稳定性与频率响应

前面的都是工具，这一节才是真正把工具组装起来用的地方。一个离散系统通常用差分方程描述，而Z变换的任务就是把这个差分方程变成传递函数，再从传递函数里读出所有我们关心的信息。

4.1 传递函数 H(z) 是怎么导出来的

一般的线性常系数差分方程写成：

Σ (k=0 到 N) aₖ·y[n-k] = Σ (k=0 到 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^N：

H(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) = 2，B = 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] = 1，x[1] = 0.75，x[2] = 0.4375。用上一节的闭式解验算：x[1] = 2×0.5 - 0.25 = 0.75，x[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₀ = 1，a₁ = -1.6，a₂ = 0.89，b₀ = 0.0725，b₁ = 0.145，b₂ = 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 算：分子是 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, dtype=float) 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, worN=4096) 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--', linewidth=1) axes[0].plot(z.real, z.imag, 'o', markersize=8, label='零点') axes[0].plot(p.real, p.imag, 'x', markersize=10, label='极点') axes[0].set_aspect('equal') axes[0].grid(True, alpha=0.3) axes[0].legend() axes[0].set_title('极零点分布') # 幅频响应 axes[1].plot(w / np.pi, mag_db) axes[1].grid(True, alpha=0.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），这些是"准临界"，需要重点做定点仿真验证。

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 问题排查速查表

这张表基本覆盖了我在实际项目里遇到的八成问题。剩下两成通常是多个原因叠加，那就得一步步来，先确认极点位置和稳定性，再确认频率响应，最后才是量化误差。排查顺序永远是从大到小，先看结构性错误，再看精度问题。

我个人的习惯是，任何数字滤波器设计完之后，先跑一个自动检查脚本：算极零点、算最大极点模长、算直流增益、跑一遍冲激响应看有没有发散。这四个检查加起来不到二十行代码，但能在设计阶段就挡掉绝大部分问题，比烧到板子上之后对着示波器发呆划算得多。极点模长这个数字，我现在基本是设计完第一眼就要看的，比任何曲线都直观。