滤波器设计（二）模拟到数字

系列文章目录

• 【音频处理】如何“认识”一个滤波器？
• 【音频处理】IIR滤波器设计（一）
• 【音频处理】IIR滤波器设计（二）模拟到数字

前言

在 【音频处理】IIR滤波器设计（一） 中，我们介绍了多种滤波器，并给出它们的差分方程、变换方程等。针对每种滤波器，我们都举了一个具体的实例来说明。同时，还讨论了零点和极点对频响的影响，已经如何用平面几何的方法计算频响。最后引出了双二阶滤波器（Biquad），并给出了它的 C++ 实现代码。

本文将理论付诸于实践，开始制作一些音频滤波器。我们知道 Biquad 的系数决定了其频谱响应和其他性质，那么如何确定这些系数呢？本文主要介绍一种从模拟信号滤波器到数字信号滤波器的转换方法。

一、模拟信号滤波器到数字信号滤波器的转换

在数字信号普及以前，科学家们在模拟信号领域做了很多研究，设计非常多优秀的滤波器。我们只需要找到一种能够将模拟信号滤波器转换为数字信号滤波器的方法，那么就能够高效、便捷地进行数字信号滤波器的设计。当然，滤波器设计还有很多其他的办法，在一些其他书籍中（例如 Audio Effects: Theory, Implementation and Application）你会看到其他设计方法。

虽然模拟信号滤波器设计并不在我们的讨论范围内，但这两个设计世界中有许多相似之处。例如模拟和数字滤波器设计中都设计到传递函数，你可以通过操作产生零点和极点，它们也都使用一个变换将时域转换到频域。一个根本的区别是，在模拟世界中，没有 Nyquist 限制，它包含 $-\infty$ 到 $+\infty$ 多有频率。另外，模拟滤波器使用电容器、电感器等元件来完成积分或者微分功能，数字滤波器则使用 delay 模块来产生偏移。下表总结了这些异同点。

数字滤波器设计	模拟滤波器设计
使用传递函数来关联输入与输出	使用传递函数来关联输入与输出
使用 delay 模块产生相位偏移和延迟	电容、电感器等元件的集成
使用 z 变换（离散时间到频域的转换）	使用 Laplace 变换（连续时间到频域的转换）
零点和极点位于 z 域中	零点和极点位于 s 域中
Nyquist 限制	所有频率
极点必须在单位圆内，以保证信号稳定	极点必须在s平面的左手部分，以保证信号稳定

二、S域与Z域

s 域与 z 域都是复数域，s 域表示为:
$$s=\sigma +j\omega$$
$\sigma \in R,\omega \in R$，它的图像长这样：

关于 z 域和 s 域的关系可以总结为：

$z=e^{s{T}_s}$ ，z 域是 s 域的进一步映射，z 变换和拉普拉斯变换通过某种映射连接在一起

由 $s=\sigma +j\omega$，当 $\sigma =0$时，$s=j\omega$，对应的是 s 域的虚轴，而此时 $x=e^{j\omega T}$ 对应的是单位圆，也就是说 z 变换将 s 域的虚轴映射成 z 域的单位圆。

当 $\sigma >0$时，$s=\sigma +j\omega$，对应的是 s 域的正半轴，而此时 $z=e^{\sigma }{e}^{j\omega T}$，由于 $e^{\sigma }>1$ ，也就是说此时 z 变换将 s 域正半轴映射到了z域的单位圆外部。

当 $\sigma <0$ 时，$s=\sigma +j\omega$，对应的是 s 域的负半轴，而此时 $z=e^{\sigma }{e}^{j\omega T}$ ，由于 $e^{\sigma }<1$，也就是说此时 z 变换将 s 域负半轴映射到了 z 域的单位圆内部。

三、模拟信号滤波器

我们的目标是将模拟信号滤波器转换为数字滤波器，因此首先需要对模拟信号滤波器进行介绍。先看一个最基本的滤波器，其传递函数为：
$$H(s)=\frac{{\omega }_c}{s+{\omega }_c}=\frac{1}{\frac{1}{{\omega }_c}s+1}$$
现在我们只考虑稳态时系统的响应，也就是让 $s=j\omega$，于是传递函数变化为：
$$H(\omega )=\frac{1}{j\frac{\omega }{{\omega }_c}+1}$$
进一步分析这个传递函数，显然：
当 $\omega \to 0$ 时，$|H(\omega )|\to 1$
当 $\omega \to {\omega }_c$ 时，$|H(\omega )|\to \frac{1}{\sqrt{2}}$
当 $\omega \to \infty$ 时，$|H(\omega )|\to 0$
这说明，$H(\omega )$在低频中响应最大，在 ${\omega }_c$ 处响应衰减至 $\frac{1}{\sqrt{2}}$，当频率增加到无限大时，响应几乎是 0。因此这个是以低通滤波器。从数学角度开看，该滤波器幅度响应为：
$$|H(\omega )|=\frac{1}{\sqrt{{\left(\omega /{\omega }_c\right)}^2+1}}=\frac{{\omega }_c}{\sqrt{{\omega }^2+{\omega }_c^2}}$$
相位响应为：
$$\angle H(\omega )=-\mathrm{atan}\left(\frac{\omega }{{\omega }_c}\right)$$
画图长这样：

现在有了一个低通滤波器，那么高通、带通、带阻等滤波器要怎么设计呢？具体详细的步骤大家可以参考 如何快速设计一个IIR滤波器。我这边就不再赘述，直接摘取总结：

我们有了一个标准的低通模拟滤波器 $H(s)=\frac{1}{s+1}$，通过适当的变形，我们可以得到相应的高通、带通及带阻模拟滤波器。但是我们现在是要设计数字滤波器啊，这不是跑题了吗？——也不见得，按照我们的直觉，模拟滤波器和数字滤波器应该存在某种联系，如果我们找到了这个联系，就可以通过模拟滤波器来得到数字滤波器，问题就解决了。那这个联系是什么呢？——双线性变换。

四、双线性变换

现在我们已经知道了 z 域是 s 域的扩展，它们之间通过某种神奇的映射联系在一起。在 如何理解离散傅里叶变换及Z变换
详细介绍了 z 与 s 之间的关系，具体来说：
$$z=e^{s{T}_s}$$
其中 $T_s=1/f_s$, fs 为采样率。我们对这个公式进行变换：
$$\begin{array}{l}\ln (z)=\ln \left(e^{sT}\right)\\ \ln (z)=sT\\ \text{or}\\ sT=\ln (z)\end{array}$$
其中 $\ln (z)$ 的泰勒展开为：
$$sT=2\left[\frac{z-1}{z+1}+\frac{1}{3}{\left(\frac{z-1}{z+1}\right)}^3-\frac{1}{5}{\left(\frac{z-1}{z+1}\right)}^5+\frac{1}{7}{\left(\frac{z-1}{z+1}\right)}^7-\dots \right]$$
我们只取第一项，对结果进行近似，得到：
$$s=\frac{2}{T}\frac{z-1}{z+1}$$
这就是双线性变化了。同时可以变换得到：
$$z=\frac{1+sT/2}{1-sT/2}$$
当 $s=0$ 时，即 0 频率，此时 $z=1$
当 $s\to \infty$ 时，即无穷大频率，此时 $z=-1$，对应频率 为 $\frac{f_s}{2}$
也就是说，双线性变换的本质是将 s 域中无穷大的频率映射到 z 域中的 $\frac{f_s}{2}$。画个图，长这样：

也就是说，在双线性变换下，s 域中的频率与 z 域中频率不同，发生了一定的弯曲。这也意味着截止频率在 s 域和 z 域是不一样的，我们需要找到这种截止频率的对应关系：
$$\begin{array}{rl}z& =e^{sT}=e^{j{\omega }_{d}T}\end{array}$$
$$\begin{array}{rl}s& =\frac{2}{T}\frac{z-1}{z+1}\\ & =\frac{2}{T}\frac{e^{j{\omega }_{d}T}-1}{e^{j{\omega }_{d}T}+1}\\ \text{令}{\omega }_{d}T=\theta & \\ s& =\frac{2}{T}\frac{e^{j\theta }-1}{e^{j\theta }+1}\\ & =\frac{2}{T}\frac{e^{j\theta /2}(e^{j\theta /2}-e^{-j\theta /2})}{e^{j\theta /2}(e^{j\theta /2}+e^{-j\theta /2})}\\ & =\frac{2}{T}\frac{(e^{j\theta /2}-e^{-j\theta /2})}{(e^{j\theta /2}+e^{-j\theta /2})}\\ \text{利用欧拉公式}& e^{jx}=\cos x+j\sin x\\ s& =\frac{2}{T}\frac{j\sin \theta /2}{\cos \theta /2}\\ & =\frac{2}{T}j\tan \theta /2\end{array}$$
又因为 $s=j{\omega }_a$(只考虑稳态时)，因此有：
$$\begin{array}{rl}j{\omega }_a& =\frac{2}{T}j\tan \theta /2\\ & =\frac{2}{T}j\tan (\frac{{\omega }_{d}T}{2})\end{array}$$
可得，s 域与 z 域截止频率之间的关系为：
$${\omega }_a=\frac{2}{T}\tan (\frac{{\omega }_{d}T}{2})$$
其中 ${\omega }_a$ 表示模拟频率，${\omega }_d$ 表示映射在 z 域中的频率，$T=\frac{1}{f_s}$

五、滤波器设计的步骤

有了前面的铺垫，接下来我们就能够按照步骤进行滤波器设计了

1. 选定归一化的低通滤波器，即$H(s)=\frac{1}{s+1}$
2. 确定数字滤波器的截止频率 ${\omega }_d$
3. 通过 ${\omega }_a=\frac{2}{T}\tan (\frac{{\omega }_{d}T}{2})$ 计算模拟频率
4. 采取合适的变换得到 $H^{\prime }(s)$，例如低通滤波器变换为 $s=\frac{s}{{\omega }_a}$
5. 在 $H^{\prime }(s)$ 中采取 $s=\frac{2}{T}\frac{z-1}{z+1}$ 进行变化，得到数字滤波器 $H(z)$。
6. 用代数方法处理数字传递函数H(z)，使其变成熟悉的形式，以便你能够 识别系数，这往往是最困难的部分。

举个例子，假设我们要设计一个滤波器，截止频率为$f_c=1000hz$，采样率为 $f_s=44100$，步骤为：

1. 选定归一化的低通滤波器，即$H(s)=\frac{1}{s+1}$
2. 截止频率 ${\omega }_d=2\pi \ast 1000=6283.2$rad/sec
3. ${\omega }_a=\frac{2}{T}\tan (\frac{{\omega }_{d}T}{2})=6293.85$ rad/sec
4. 对 $H(x)$ 中 s 进行替换，$s\to \frac{s}{{\omega }_a}$，得到 $H^{\prime }(s)=\frac{{\omega }_a}{s+{\omega }_a}$
5. 用 $s=\frac{2}{T}\frac{z-1}{z+1}$ 对 $H^{\prime }(s)$ 进行替换，得到 $H^{\prime }(s)=\frac{0.0667(z+1)}{z-0.8667}$
6. 整理成熟悉的形式 $H(z)=\frac{0.0667+0.0667z^{-1}}{1-0.8667z^{-1}}=\frac{a_0+a_1{z}^{-1}}{1+b_1{z}^{-1}}$

根据最终的结果，我们可以确认该滤波器在 Biquad 中的参数，即
$$\begin{gathered} a_0=0.0667 \\ {a}_1=0.0667 \\ {b}_1=-0.8667 \end{gathered}$$
该滤波器的差分方程为：
$$y(n)=0.0667x(n)+0.0667x(n-1)+0.8667y(n-1)$$
频谱响应长这样（通过 Filter playground绘制得到）

总结

本文介绍了一种将模拟信号滤波器转换为数字信号滤波器的方法，首先说明了 s 域与 z 域之间的关系，利用双线性变换来近似地转换 s 域与 z 域，然后推导出了截止频率在 s 域与 z 域之间的关系，最后给出了数字滤波器设计的具体步骤。

参考

如何理解离散傅里叶变换及Z变换

如何快速设计一个IIR滤波器