IIR 滤波器设计方法对比

脉冲响应不变法 vs. 双线性变换法

直观理解 · 深入对比 · 交互体验

💡 为什么要从模拟滤波器出发？

数字滤波器的设计其实是"站在巨人的肩膀上"！在数字信号处理出现之前，模拟滤波器（使用电阻、电容、电感等元件）已经发展得非常成熟，有很多经典的设计方法和成熟的理论基础，比如巴特沃斯滤波器、切比雪夫滤波器等。

IIR数字滤波器设计的核心思想：既然模拟滤波器已经设计得很好，为什么不直接"借用"它们的优秀特性呢？我们要做的就是找到一种方法，把模拟滤波器的设计"转换"成数字滤波器。这就像把一张优秀的纸质地图（模拟）转化成电子地图（数字）——关键是找到好的转换方法！

本文对比的两种方法：

🎯 用生活场景理解两种方法

🧵 脉冲响应不变法 = 折叠长布

想象：一条很长的彩色布料（模拟信号）需要装进一个固定大小的盒子（数字系统）。
做法：将布料折叠多次，使其适应盒子大小。
结果：不同位置的图案会叠加在一起（混叠），原始图案部分丢失。

🗺️ 双线性变换法 = 地图投影

想象：将球形的地球（无限大的模拟频率）投影到平面地图上（有限的数字频率）。
做法：使用特殊的投影方法（双线性变换），将整个球面压缩映射到平面。
结果：虽然形状扭曲变形（频率畸变），但每个点都有唯一对应，没有重叠。

脉冲响应不变法 (IIM)

📝 核心原理与步骤

IIM 不是简单的变量代换，而是基于部分分式展开的变换。

Step 1
对 $H_a(s)$ 进行部分分式展开：
$$H_a(s) = \sum_{k=1}^{N} \frac{A_k}{s - p_k}$$

其中 $p_k$ 是单极点

Step 2
利用变换对进行映射：
$$\frac{1}{s - p_k} \xrightarrow{IIM} \frac{1}{1 - e^{p_k T}z^{-1}}$$

*注: 有时分子会乘以 $T$ 保持增益一致

Step 3
得到最终 H(z)：
$$H(z) = \sum_{k=1}^{N} \frac{A_k}{1 - e^{p_k T}z^{-1}}$$

1. S平面到Z平面的映射关系（动画演示）

🔍 映射原理解析：为什么只显示 [-π/T, π/T] 到单位圆？

关键理解：

  • 时域采样导致频域周期化：由于时域采样 $h[n] = h(nT)$，数字频率响应 $H(e^{j\omega})$ 是模拟频率响应 $H_a(j\Omega)$ 的周期延拓，周期为 $2\pi/T$
  • 主频段 [-π/T, π/T]：这是模拟频率的"基本周期"，完整映射到一圈单位圆（$\omega \in [-\pi, \pi]$）
  • 其他频段的命运：由于周期性，$[\pi/T, 3\pi/T]$、$[-3\pi/T, -\pi/T]$ 等所有频段都会重复映射到同一个单位圆上
  • 与原理的关系：映射公式 $z = e^{sT}$ 在虚轴上变为 $z = e^{j\Omega T}$，这是一个周期函数，周期为 $2\pi/T$
⚠️ 混叠产生的根本原因：
  • S平面虚轴上相距 $2\pi/T$ 的所有点（如 $j\Omega$、$j(\Omega + 2\pi/T)$、$j(\Omega - 2\pi/T)$）都映射到Z平面的同一点
  • 这就像把无限长的纸条按固定长度反复折叠，不同段会叠在一起
  • 从图中可以看到：多个彩色条带（代表不同频段）的箭头都指向同一个单位圆，产生"多对一"映射

2. 设计过程与详细示例

例题：已知一阶模拟低通滤波器 $H_a(s) = \frac{2}{s+2}$，采样周期 $T=1$ 秒，使用IIM方法设计数字滤波器。

步骤1：部分分式展开 Ha(s)

理论：将系统函数展开为简单极点形式

$$H_a(s) = \sum_{k=1}^{N} \frac{A_k}{s - p_k}$$

其中 $p_k$ 是极点，$A_k$ 是留数

本例：$H_a(s) = \frac{2}{s+2} = \frac{2}{s-(-2)}$ 已是标准形式

• 极点：$p_1 = -2$（位于左半平面，稳定）

• 留数：$A_1 = 2$

• 系统阶数：$N = 1$
步骤2：应用IIM变换公式

变换规则：

$$\frac{A_k}{s-p_k} \rightarrow \frac{A_k T}{1-e^{p_k T}z^{-1}}$$

推导过程：

1) 逆变换：$h(t) = A_k e^{p_k t}$

2) 采样：$h[n] = A_k e^{p_k nT}$

3) z变换：求和得上式

代入参数：$A_1=2$, $p_1=-2$, $T=1$

$$H(z) = \frac{2 \cdot 1}{1 - e^{-2 \cdot 1}z^{-1}}$$

$$= \frac{2}{1 - e^{-2}z^{-1}}$$

计算过程详解：
• 分子：$A_1 T = 2 \times 1 = 2$
• 分母中的指数：$e^{p_1 T} = e^{(-2) \times 1} = e^{-2}$
步骤3：数值计算与验证

数值计算：计算指数项

验证要点：

  • 极点位置（稳定性）
  • 零频增益

$e^{-2} \approx 0.1353352832...$

取4位有效数字：$e^{-2} \approx 0.1353$

$$H(z) = \frac{2}{1 - 0.1353z^{-1}}$$

稳定性：极点 $z_1 = 0.1353$，$|z_1| < 1$ ✓

零频增益：

• 模拟：$H_a(0) = 2/2 = 1$

• 数字：$H(1) = 2/(1-0.1353) \approx 2.31$

⚠️ 零频增益不匹配！IIM保持脉冲响应采样点一致，但不保证频率响应完全相同。

3. 优缺点总结

  • 时域脉冲响应在采样点上完全匹配
  • 频率映射为线性关系：ω = ΩT
  • 物理意义清晰，易于理解
  • 固有的频率混叠问题
  • 仅适用于带限滤波器（低通、带通）
  • 不适合高通和带阻滤波器

双线性变换法 (BLT)

📝 核心原理与步骤

BLT 是纯粹的代数代换，基于梯形积分公式。

Step 1
写出变换公式：
$$s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}$$

这是 $z = \frac{1+sT/2}{1-sT/2}$ 的逆变换

Step 2
频率预畸变 (Pre-warping)：
由于频率映射是非线性的，需修正模拟截止频率 $\Omega$：
$$\Omega_{pre} = \frac{2}{T}\tan\left(\frac{\omega_c}{2}\right)$$
Step 3
直接代入替换：
$$H(z) = H_a(s)\Big|_{s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}$$

1. S平面到Z平面的映射关系（动画演示）

🔍 映射原理解析：为什么是三步过程？

完整的映射链条：

  • 第一步（S→S₁）：频率预畸变
    • 将S平面整个虚轴（$j\Omega, -\infty < \Omega < \infty$）通过反正切函数压缩到S₁平面的有限范围 $[-\pi/T, \pi/T]$
    • 映射关系：$\Omega_1 = \frac{2}{T}\tan(\frac{\Omega T}{2})$（反向预畸变公式）
    • 这一步是非线性压缩，将无限频率范围映射到有限范围
  • 第二步（S₁→Z）：标准线性映射
    • S₁平面的 $[-\pi/T, \pi/T]$ 区间线性映射到Z平面单位圆 $[-\pi, \pi]$
    • 这一步类似IIM，但因为S₁已经是有限范围，所以不会有混叠
  • 与原理的关系：BLT的代数公式 $s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$ 实际上隐含了这两步变换，通过巧妙的数学变换一次性完成
✓ 无混叠的秘密：
  • BLT先把整个S平面虚轴（无限长）压缩到S₁平面的有限区间 $[-\pi/T, \pi/T]$
  • 然后这个有限区间一对一映射到Z平面单位圆，每个频率都有唯一对应
  • 代价是频率非线性畸变：高频被严重压缩（地图投影的"极地变形"）
  • 从图中看：S平面整个左半平面（阴影区域）→ 中间压缩的S₁平面 → Z平面单位圆内，是完整的"一对一"映射

2. 设计过程与详细示例

例题：同样从 $H_a(s) = \frac{2}{s+2}$，$T=1$ 出发，使用BLT方法设计。

步骤1：确定BLT变换公式

标准形式：

$$s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$$

或：$$s = \frac{2}{T}\frac{z-1}{z+1}$$

来源：梯形积分法

代入 T=1：

$$s = 2\frac{1-z^{-1}}{1+z^{-1}}$$

系数 2/T 保证了s平面虚轴映射到z平面单位圆
步骤2：代入并展开

操作：将s的表达式代入 $H_a(s)$

$$H(z) = H_a(s)\Big|_{s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}$$

$$H(z) = \frac{2}{2\frac{1-z^{-1}}{1+z^{-1}}+2}$$

分母展开：

$$= 2\frac{1-z^{-1}}{1+z^{-1}}+2$$

$$= 2\left(\frac{1-z^{-1}+(1+z^{-1})}{1+z^{-1}}\right)$$

$$= 2\cdot\frac{2}{1+z^{-1}} = \frac{4}{1+z^{-1}}$$
步骤3：化简为标准形式

目标：化为

$$H(z) = \frac{b_0 + b_1z^{-1} + ...}{1 + a_1z^{-1} + ...}$$

$$H(z) = \frac{2}{\frac{4}{1+z^{-1}}} = \frac{2(1+z^{-1})}{4}$$

$$= \frac{1+z^{-1}}{2}$$

$$H(z) = 0.5(1 + z^{-1})$$

分析：

• 零点：$z = -1$（单位圆上）

• 极点：无（FIR滤波器！）

• 零频增益：$H(1) = 0.5 \times 2 = 1$ ✓

💡 有趣现象：IIR模拟滤波器变成了FIR数字滤波器！这是因为只有一个极点且T较大。

3. 优缺点总结

  • 完全消除混叠（一对一映射）
  • 适用于所有滤波器类型
  • 稳定性保持（稳定→稳定）
  • 设计简单（代数替换）
  • 频率畸变（非线性映射）
  • 需要频率预畸变校正
  • 脉冲响应与原型不同

🎮 交互式频率变换对比

拖动滑块改变模拟频率 Ω，在同一张图上实时观察两种方法的映射曲线。
蓝色直线 = IIM（线性，会混叠） | 红色曲线 = BLT（非线性，无混叠）

IIM：ω = Ω·T

0.000

BLT：ω = 2·arctan(Ω·T/2)

0.000

✓ 始终在 (-π, π) 内

🎬 深入理解：什么是混叠（Aliasing）？

🌟 通俗理解：车轮倒转之谜

你有没有看过电影或视频中，汽车加速时，轮子看起来像是在倒着转？这就是混叠的经典例子！

原理：摄像机每秒拍摄24帧（采样率），如果轮子转得太快，每次拍摄时轮辐的位置变化超过一定角度，大脑就会误以为轮子在倒转。实际上轮子一直在正向旋转，只是我们的"采样"（拍摄）不够快，无法正确捕捉真实的运动！

核心问题：当我们用有限的"快照"（采样点）去描述连续的运动时，如果运动太快（高频成分），就可能被错误地理解成慢速运动（低频成分）——这就是高频混叠成低频！

📐 技术原理：采样定理与频率折叠

奈奎斯特采样定理（Nyquist-Shannon Sampling Theorem）：

要完整恢复一个频率为 $f$ 的信号，采样频率 $f_s$ 必须满足： $$f_s > 2f$$ 其中 $f_s/2$ 被称为奈奎斯特频率，是采样系统能够正确表示的最高频率。

混叠的数学表达：

当模拟信号包含频率 $f > f_s/2$ 的成分时，这些高频成分会在采样后被"折叠"到低频区域： $$f_{aliased} = |f \bmod f_s - f_s/2|$$ 例如：如果采样率是 30 Hz（每秒30帧），奈奎斯特频率是 15 Hz。一个 20 Hz 的真实频率会被混叠成 10 Hz 的假频率！

在IIR滤波器设计中：

脉冲响应不变法通过时域采样实现，导致频域上出现周期延拓： $$H(e^{j\omega}) = \frac{1}{T}\sum_{k=-\infty}^{\infty} H_a\left(j\frac{\omega - 2\pi k}{T}\right)$$ 这意味着所有相距 $2\pi/T$ 的模拟频率成分都会叠加在同一个数字频率上，产生混叠失真！

📹 混叠现象视频演示

下面三个视频展示了经典的混叠现象：旋转物体在不同采样率下的视觉效果。注意观察物体的旋转方向和速度如何随着采样率变化而"失真"。

🚁 直升机旋翼

🎡 玩具风车

🛞 货车车轮

💡 观察要点：