系列文章目录

【音频处理】如何“认识”一个滤波器？

【音频处理】IIR滤波器设计（一）Biquad 滤波器

【音频处理】IIR滤波器设计（二）模拟到数字

前言

在开始学习 IIR 滤波器之前，你还记得 z变换吗？它可以简化数字滤波器的分析与设计，接下来内容将建立在 z变换的基础上进行推导。如果你忘记了，那么可以在【音频处理】如何“认识”一个滤波器？复习下，除了z变换外，这篇博客中也提到的零点、极点等概念，相信这些会对你有帮助的。

接下来，将从最简单的滤波器开始介绍，并最终引出今天的主角双二阶滤波器（Biquad Filter），要介绍的滤波器包括：

一阶前馈滤波器（1st order feed-forward filter）
一阶反馈滤波器（1st order feedback filter)
二阶前馈滤波器（2nd order feed-forward filter）
二阶反馈滤波器（2nd order feedback filter）
搁架式均衡器（Shelving filter）
双二阶滤波器（Biquad Filter）

一阶前馈滤波器（1st order feed-forward filter）

作为最简单的滤波器，一阶前馈滤波器的差分方程非常易懂:
y(n)=a0x(n)+a1x(n−1)y(n) = a_0x(n) + a_1x(n-1) y(n)=a0x(n)+a1x(n−1)

z变换后：
H(z)=Y(z)X(z)=a0+a1z−1H(z) = \frac{Y(z)}{X(z)} = a_0 + a_1z^{-1} H(z)=X(z)Y(z)=a0+a1z−1

其中 z=ejωz = e^{j\omega}z=ejω

计算频谱响应

我们假设 a0=a1=0.5a_0 = a_1 = 0.5a0=a1=0.5
H(z)=Y(z)X(z)=0.5+0.5z−1=0.5+0.5zH(z) = \frac{Y(z)}{X(z)} = 0.5 + 0.5z^{-1} = 0.5 + \frac{0.5}{z} H(z)=X(z)Y(z)=0.5+0.5z−1=0.5+z0.5

现在来计算它的频谱响应，我们取四个点：

DC: 0
Nyquist: π\piπ
12\frac{1}{2}21 Nyquist: π/2\pi/2π/2
14\frac{1}{4}41 Nyquist: π/4\pi/4π/4

DC(0 Hz)

ω=0e−jω=e0=1H(ω)=0.5+0.5e−jω=1∣H(ω)∣=1\begin{aligned} &\omega = 0 \\ &e^{-j\omega} = e^{0} = 1 \\ &H(\omega) = 0.5 + 0.5e^{-j\omega} = 1 \\ \vert &H(\omega) \vert = 1 \end{aligned} ∣ω=0e−jω=e0=1H(ω)=0.5+0.5e−jω=1H(ω)∣=1

Nyquist: π\piπ

ω=πe−jω=e−jπ=cos⁡(π)−jsin⁡(π)=−1H(ω)=0.5+0.5e−jω=0∣H(ω)∣=0\begin{aligned} &\omega = \pi \\ &e^{-j\omega} = e^{-j\pi} = \cos(\pi) - j\sin(\pi) = -1 \\ &H(\omega) = 0.5 + 0.5e^{-j\omega} = 0 \\ \vert &H(\omega) \vert = 0 \end{aligned} ∣ω=πe−jω=e−jπ=cos(π)−jsin(π)=−1H(ω)=0.5+0.5e−jω=0H(ω)∣=0

12\frac{1}{2}21 Nyquist: π/2\pi/2π/2

ω=π/2e−jω=e−jπ/2=cos⁡(π/2)−jsin⁡(π/2)=−j1H(ω)=0.5+0.5e−jω=0.5−j0.5∣H(ω)∣=0.7071\begin{aligned} &\omega = \pi/2 \\ &e^{-j\omega} = e^{-j\pi/2} = \cos(\pi/2) - j\sin(\pi/2) = -j1\\ &H(\omega) = 0.5 + 0.5e^{-j\omega} = 0.5 - j0.5 \\ \vert &H(\omega) \vert = 0.7071 \end{aligned} ∣ω=π/2e−jω=e−jπ/2=cos(π/2)−jsin(π/2)=−j1H(ω)=0.5+0.5e−jω=0.5−j0.5H(ω)∣=0.7071

14\frac{1}{4}41 Nyquist: π/4\pi/4π/4

ω=π/4e−jω=e−jπ/4=cos⁡(π/4)−jsin⁡(π/4)=0.707−j0.707H(ω)=0.5+0.5e−jω=0.5+0.5(0.707−j0.707)=0.835−j0.353∣H(ω)∣=0.9065\begin{aligned} &\omega = \pi/4 \\ &e^{-j\omega} = e^{-j\pi/4} = \cos(\pi/4) - j\sin(\pi/4) = 0.707 - j0.707\\ &H(\omega) = 0.5 + 0.5e^{-j\omega} = 0.5 + 0.5(0.707 - j0.707) = 0.835 - j0.353 \\ \vert &H(\omega) \vert = 0.9065 \end{aligned} ∣ω=π/4e−jω=e−jπ/4=cos(π/4)−jsin(π/4)=0.707−j0.707H(ω)=0.5+0.5e−jω=0.5+0.5(0.707−j0.707)=0.835−j0.353H(ω)∣=0.9065

将上述 4 个点连起来后，可以得当 a0=a1=0.5a_0 = a_1 = 0.5a0=a1=0.5 时滤波器的频响曲线

零点与频响的几何关系

当 a0=0.5a_0 = 0.5a0=0.5 a1=0.5a_1 = 0.5a1=0.5 时
H(z)=Y(z)X(z)=0.5+0.5z−1=0.5+0.5zH(z) = \frac{Y(z)}{X(z)} = 0.5 + 0.5z^{-1} = 0.5 + \frac{0.5}{z} H(z)=X(z)Y(z)=0.5+0.5z−1=0.5+z0.5

滤波器存在一个零点 z=−1+j0z= -1 + j0z=−1+j0

接下来，我们用几何的方式来确定频谱响应与零点直接的关系。对于只有一个零点的一阶前馈滤波器，可以这样确定频响：

给定一个频率，在复平面坐标的单位圆外圈定位它
画一条线，连接该点与零点
线的长度将是该频率的增幅

以 000，π\piπ，π/2\pi/2π/2, π/4\pi/4π/4 为例，用上述方法进行计算，可以得到下图

你可能会注意到，用几何关系得到的值与直接计算得到结果是不同的。在 0Hz 中其增益为 2.0，但实际上只有它的一半。事实上，几何方法得到的增益是实际幅度的两倍。为了得到正确的结果，我们需要做的最后一件事是从H(z)H(z)H(z)中去掉整体增益系数，这样整体的增益（或者衰减）就可以只由一个变量控制。这实现起来相当简单，提取一个 a0a_0a0 作为公因子即可：

H(z)=a0+a1z−1=a0(1+a1a0z−1)Let α1=a1a0H(z)=a0(1+α1z−1)\begin{aligned} H(z)&= a_0 + a_1z^{-1} \\ &= a_0(1 + \frac{a_1}{a_0}z^{-1}) \\ \text{Let } &\alpha_1 = \frac{a_1}{a_0} \\ H(z)&= a_0(1 + \alpha_1z^{-1}) \end{aligned} H(z)Let H(z)=a0+a1z−1=a0(1+a0a1z−1)α1=a0a1=a0(1+α1z−1)

因此，在复平面得到幅度后，需要做的最后一步是用 a0a_0a0 来缩放它，这样结果就相互匹配了。

一阶反馈滤波器（1st order feedback filter)

接下来介绍 一阶反馈滤波器，它与 一阶前馈滤波器 有很多相似之处，但也会有几个关键区别。

它的差分方程为：
y(n)=a0x(n)−b1y(n−1)y(n) = a_0x(n) - b_1y(n-1) y(n)=a0x(n)−b1y(n−1)

进行z变换：
Y(z)=a0X(z)−b1Y(z)z−1Y(z) = a_0X(z) - b_1Y(z)z^{-1} Y(z)=a0X(z)−b1Y(z)z−1
得到转换方程：
H(z)=Y(z)X(z)=a01+b1z−1H(z) = \frac{Y(z)}{X(z)} = \frac{a_0}{1+b_1z^{-1}} H(z)=X(z)Y(z)=1+b1z−1a0
接着，提取 a0a_0a0 作为增益控制系数：
H(z)=a01+b1z−1=a011+b1z−1\begin{aligned} H(z) &= \frac{a_0}{1+b_1z^{-1}} \\ &= a_0\frac{1}{1+b_1z^{-1}} \end{aligned} H(z)=1+b1z−1a0=a01+b1z−11

我们再做一些变化，方便计算零点和极点：
H(z)=a01+b1z−1=a011+b1z=a0zz+b1\begin{aligned} H(z) &= \frac{a_0}{1+b_1z^{-1}} \\ &= a_0\frac{1}{1+\frac{b_1}{z}} \\ &= a_0\frac{z}{z + b_1} \end{aligned} H(z)=1+b1z−1a0=a01+zb11=a0z+b1z

可以看到，当 z=−b1z=-b_1z=−b1 时分母为零，也就是说极点在 z=−b1z=-b_1z=−b1 位置。你可能还会注意到传递函数还有一个零点 z=0z=0z=0，这个零点被称为 “trivial zero”，它对频谱响应没有任何影响。

在z=0处的极点或零点是平凡的（trivial），可以为了分析而忽略，因为它对频率响应没有影响。

我们假设 a0=1.0a_0 = 1.0a0=1.0 b1=0.9b_1 = 0.9b1=0.9
H(z)=Y(z)X(z)=a01+b1z−1=11+0.9z−1H(z) = \frac{Y(z)}{X(z)} = \frac{a_0}{1+b_1z^{-1}} = \frac{1}{1 + 0.9z^{-1}} H(z)=X(z)Y(z)=1+b1z−1a0=1+0.9z−11
接下来开始计算该滤波器的频响曲线。

极点与频响的几何关系

当 a0=1.0a_0 = 1.0a0=1.0 b1=0.9b_1 = 0.9b1=0.9 时，滤波器存在一个极点 z=−0.9+j0z = -0.9 + j0z=−0.9+j0，对于只有一个零点的一阶反馈滤波器，可以这样确定频响：

给定一个频率，在复平面坐标的单位圆外圈定位它
画一条线，连接该点与极点
线段长度的倒数将是频率的增幅

以 000，π\piπ，π/2\pi/2π/2, π/4\pi/4π/4 为例，用上述方法进行计算，可以得到下图