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

系列文章目录

• 【音频处理】如何“认识”一个滤波器？
• 【音频处理】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)=a_{0}x(n)+a_{1}x(n-1)$$

z变换后：
$$H(z)=\frac{Y(z)}{X(z)}=a_0+a_1{z}^{-1}$$

其中 $z=e^{j\omega }$

计算频谱响应

我们假设 $a_0=a_1=0.5$
$$H(z)=\frac{Y(z)}{X(z)}=0.5+0.5z^{-1}=0.5+\frac{0.5}{z}$$

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

• DC: 0
• Nyquist: $\pi$
• $\frac{1}{2}$ Nyquist: $\pi /2$
• $\frac{1}{4}$ Nyquist: $\pi /4$

DC(0 Hz)

$$\begin{array}{rl}& \omega =0\\ & e^{-j\omega }=e^0=1\\ & H(\omega )=0.5+0.5e^{-j\omega }=1\\ |& H(\omega )|=1\end{array}$$

Nyquist: $\pi$

$$\begin{array}{rl}& \omega =\pi \\ & e^{-j\omega }=e^{-j\pi }=\cos (\pi )-j\sin (\pi )=-1\\ & H(\omega )=0.5+0.5e^{-j\omega }=0\\ |& H(\omega )|=0\end{array}$$

$\frac{1}{2}$ Nyquist: $\pi /2$

$$\begin{array}{rl}& \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\\ |& H(\omega )|=0.7071\end{array}$$

$\frac{1}{4}$ Nyquist: $\pi /4$

$$\begin{array}{rl}& \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\\ |& H(\omega )|=0.9065\end{array}$$

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

零点与频响的几何关系

当 $a_0=0.5$ $a_1=0.5$ 时
$$H(z)=\frac{Y(z)}{X(z)}=0.5+0.5z^{-1}=0.5+\frac{0.5}{z}$$

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

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

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

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

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

$$\begin{array}{rl}H(z)& =a_0+a_1{z}^{-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 }_1{z}^{-1})\end{array}$$

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

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

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

它的差分方程为：
$$y(n)=a_{0}x(n)-b_{1}y(n-1)$$

进行z变换：
$$Y(z)=a_{0}X(z)-b_{1}Y(z)z^{-1}$$
得到转换方程：
$$H(z)=\frac{Y(z)}{X(z)}=\frac{a_0}{1+b_1{z}^{-1}}$$
接着，提取 $a_0$ 作为增益控制系数：
$$\begin{array}{rl}H(z)& =\frac{a_0}{1+b_1{z}^{-1}}\\ & =a_0\frac{1}{1+b_1{z}^{-1}}\end{array}$$

我们再做一些变化，方便计算零点和极点：
$$\begin{array}{rl}H(z)& =\frac{a_0}{1+b_1{z}^{-1}}\\ & =a_0\frac{1}{1+\frac{b_1}{z}}\\ & =a_0\frac{z}{z+b_1}\end{array}$$

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

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

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

极点与频响的几何关系

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

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

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

计算频谱响应

我们假设 $a_0=1.0$ $b_1=0.9$
$$H(z)=\frac{Y(z)}{X(z)}=\frac{a_0}{1+b_1{z}^{-1}}=\frac{1}{1+0.9z^{-1}}$$

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

• DC: 0
• Nyquist: $\pi$
• $\frac{1}{2}$ Nyquist: $\pi /2$
• $\frac{1}{4}$ Nyquist: $\pi /4$

DC(0 Hz)

$$\begin{array}{rl}& \omega =0\\ & e^{-j\omega }=e^0=1\\ & H(\omega )=\frac{1}{1+0.9e^{-j\omega }}=\frac{1}{1.9}=0.5263\\ |& H(\omega )|=0.5263\end{array}$$

Nyquist: $\pi$

$$\begin{array}{rl}& \omega =\pi \\ & e^{-j\omega }=e^{-j\pi }=\cos (\pi )-j\sin (\pi )=-1\\ & H(\omega )=\frac{1}{1+0.9e^{-j\omega }}=\frac{1}{0.1}=10\\ |& H(\omega )|=10\end{array}$$

$\frac{1}{2}$ Nyquist: $\pi /2$

$$\begin{array}{rl}& \omega =\pi /2\\ & e^{-j\omega }=e^{-j\pi /2}=\cos (\pi /2)-j\sin (\pi /2)=-j1\\ & H(\omega )=\frac{1}{1+0.9e^{-j\omega }}=\frac{1}{1-j0.9}\\ |& H(\omega )|=\frac{|1|}{|1-j0.9|}=0.743\end{array}$$

$\frac{1}{4}$ Nyquist: $\pi /4$

$$\begin{array}{rl}& \omega =\pi /4\\ & e^{-j\omega }=e^{-j\pi /4}=\cos (\pi /4)-j\sin (\pi /4)=0.707-j0.707\\ & H(\omega )=\frac{1}{1+0.9e^{-j\omega }}=\frac{1}{1+0.9(0.707-j0.707)}=\frac{1}{1.636-j0.636}\\ |& H(\omega )|=\frac{|1|}{|1.636-j0.636|}=0.57\end{array}$$

将上述 4 个点连起来后，可以得当 $a_0=1.0$ $b_1=0.9$ 时滤波器的频响曲线

当然，也可以转换成 dB 值。通过下图可知，该滤波器抑制低频信号，提升高频信息，在 0hz 处将会有 -6dB 的抑制，而在 22.5hz 则有 20dB 的增益。

二阶前馈滤波器（2nd order feed-forward filter）

二阶前馈滤波器的分析与一阶滤波器相似，但我们要处理的数学问题更多。它的差分方程为：
$$y(n)=a_{0}x(n)+a_{1}x(n-1)+a_{2}x(n-2)$$

进行z变换
$$Y(z)=a_{0}X(z)+a_{1}X(z)z^{-1}+a_{2}X(z)z^{-2}$$

得到转换方程
$$H(z)=\frac{Y(z)}{X(z)}=a_0+a_1{z}^{-1}+a_2{z}^{-2}$$

提取 $a0$
$$\begin{array}{rl}\text{Let }{\alpha }_1& =\frac{a_1}{a_0},{\alpha }_2=\frac{a_2}{a_0}\\ H(z)& =\frac{Y(z)}{X(z)}=a_0(1+{\alpha }_1{z}^{-1}+{\alpha }_2{z}^{-2})\end{array}$$

此时 $H(z)$ 是个二次方程，对它进行因式分解，可以有：
$$\begin{array}{rl}H(z)& =\frac{Y(z)}{X(z)}=a_0(1-Z_1{z}^{-1})(1-Z_2{z}^{-1})\\ & =a_0(1+(-Z_1-Z_2)z^{-1}+Z_1{Z}_2{z}^{-2})\end{array}$$

从上述式子可以得到几个有用的信息：

1. 零点在 $z=Z_1$ 和 $z=Z_2$ 上，它们在复平面上，所以应该是一个复数
2. ${\alpha }_1$ 和 ${\alpha }_2$ 是实数，因此 $(-Z_1-Z_2)$ 和 $Z_1{Z}_2$ 也应该是实数。

据此我们可以推断出 $Z_1$ 和 $Z_2$ 为共轭关系，才能满足上述两个条件。因此：
$$\begin{array}{r}{Z}_1=R{e}^{j\theta }=a+jb\\ Z_2=R{e}^{-j\theta }=a-jb\end{array}$$

将其带入转换方程，得到：
$$\begin{array}{rl}H(z)& =a_0(1+(-Z_1-Z_2)z^{-1}+Z_1{Z}_2{z}^{-2})\\ & =a_0(1-R{e}^{j\theta }{z}^{-1}-R{e}^{-j\theta }{z}^{-1}+R^2(e^{j\theta }{e}^{-j\theta })z^{-2})\\ \text{noting that }(e^{j\theta }{e}^{-j\theta })& =1\\ H(z)& =a_0(1-(R{e}^{j\theta }+R{e}^{-j\theta })z^{-1}+R^2{z}^{-2})\\ & =a_0(1-R(\cos (\Theta )+j\sin (\Theta )+\cos (\Theta )-j\sin (\Theta ))z^{-1}+R^2{z}^{-2})\\ & =a_0(1-2R\cos (\Theta )z^{-1}+R^2{z}^{-2})\end{array}$$
对比最初的转换方程 $H(z)=a_0(1+{\alpha }_1{z}^{-1}+{\alpha }_2{z}^{-2})$，可以得到：
$$\begin{array}{l}{\alpha }_1=-2R\cos (\Theta )\\ {\alpha }_2=R^2\end{array}$$

根据上述等式，给定 ${\alpha }_1$ 和 ${\alpha }_2$ 就能够求出 $R$ 和 $\Theta$，从而得到 $Z_1$ 和 $Z_2$，也就得到了零点。

我们以 $a_0=1.0$，$a_1=-1.27$，$a_2=0.81$ 为例。这时 ${\alpha }_1=-1.27$，${\alpha }_2=0.81$，那么可以得到 $R$：
$$\begin{array}{c}{R}^2={\alpha }_2=0.81\\ R=\sqrt{0.81}=0.9\end{array}$$
接着：
$$\begin{array}{rl}-2\mathrm{Rcos}(\Theta )& =-1.27\\ 2(0.9)\cos (\Theta )& =1.27\\ \cos (\Theta )& =\frac{1.27}{2(0.9)}\\ \Theta & =\arccos (0.705)\\ \Theta & =45^{\circ }\end{array}$$

至此，得到两个零点在复平面上的位置在 $0.9e^{j0.25\pi }$ 和 $0.9e^{j-0.25\pi }$ 上，也就是复平面上，半径为 0.9 ，角度为 $\pm 45^{\circ }$ 的位置，如下图：

估计频响

对于两个零点的情况，估计频响的步骤与之前类似，但多一个步骤。当估计一个以上零点的频响时，需要这么做：

1. 给定一个频率，在复平面坐标的单位圆外圈定位它
2. 画一条线，连接该点与所有零点
3. 将所有线段的长度相乘，得到的积将是频率的增幅

最后一条规则用公式来说明的话，应该是这样：
$${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=a_0\prod _{i=1}^N{U}_i$$
其中 $N$ 为滤波器的阶数（也就是零点的个数），$U_i$为与零点的距离长度。

以 $0$，$\pi$，$\pi /2$, $\pi /4$ 为例，用上述方法进行计算频响：

当 $\omega =0$ 时，与两个零点的距离为 0.7、0.7，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=0.7\ast 0.7=0.49$

当 $\omega =\pi /4$ 时，与两个零点的距离为 0.1、 1.4，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=0.1\ast 1.4=0.14$

当 $\omega =\pi /2$ 时，与两个零点的距离为 0.7、1.8，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=0.7\ast 1.8=1.26$

当 $\omega =\pi$ 时，与两个零点的距离为 1.75、1.75，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1.75\ast 1.75=3.1$

连接上述 4 个点，可以得到大致的频响曲线：

计算频谱响应

当 $a_0=1.0$，$a_1=-1.27$，$a_2=0.81$ 时
$$\begin{array}{rl}H(z)& =1-1.27z^{-1}+0.81z^{-2}\\ & =1-1.27e^{-j1\omega }+0.81e^{-j2\omega }\\ & =1-1.27[\cos (\omega )-j\sin (\omega )]+0.81[\cos (2\omega )-j\sin (2\omega )]\end{array}$$

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

• DC: 0
• Nyquist: $\pi$
• $\frac{1}{2}$ Nyquist: $\pi /2$
• $\frac{1}{4}$ Nyquist: $\pi /4$

DC: 0

H ( ω ) ∣ ω = 0 = 1 − 1.27 [ cos ⁡ ( ω ) − j sin ⁡ ( ω ) ] + 0.81 [ cos ⁡ ( 2 ω ) − j sin ⁡ ( 2 ω ) ] = 1 − 1.27 [ cos ⁡ ( 0 ) − j sin ⁡ ( 0 ) ] + 0.81 [ cos ⁡ ( 2 ∗ 0 ) − j sin ⁡ ( 2 ∗ 0 ) ] = 0.54 + j 0 ∣ H ( ω ) ∣ = a 2 + b 2 = 0.5 4 2 + 0 2 = 0.54 \begin{aligned} H(\omega)|_ {\omega=0}&= 1-1.27[\cos (\omega)-j \sin (\omega)]+0.81[\cos (2 \omega)-j \sin (2 \omega)] \\ & = 1-1.27[\cos (0)-j \sin (0)]+0.81[\cos (2 * 0)-j \sin (2 * 0)] \\ &= 0.54+j 0 \\ |H(\omega)| &=\sqrt{a^{2}+b^{2}} \\ &=\sqrt{0.54^{2}+0^{2}}=0.54 \end{aligned} H(ω)∣ω=0​∣H(ω)∣​=1−1.27[cos(ω)−jsin(ω)]+0.81[cos(2ω)−jsin(2ω)]=1−1.27[cos(0)−jsin(0)]+0.81[cos(2∗0)−jsin(2∗0)]=0.54+j0=a2+b2 ​=0.542+02 ​=0.54​

Nyquist: $\pi$

H ( ω ) ∣ ω = π = 1 − 1.27 [ cos ⁡ ( ω ) − j sin ⁡ ( ω ) ] + 0.81 [ cos ⁡ ( 2 ω ) − j sin ⁡ ( 2 ω ) ] = 1 − 1.27 [ cos ⁡ ( π ) − j sin ⁡ ( π ) ] + 0.81 [ cos ⁡ ( 2 π ) − j sin ⁡ ( 2 π ) ] = 3.08 + j 0 ∣ H ( ω ) ∣ = a 2 + b 2 = 3.0 8 2 + 0 2 = 3.08 \begin{aligned} \left.H(\omega)\right|_{\omega=\pi} &=1-1.27[\cos (\omega)-j \sin (\omega)]+0.81[\cos (2 \omega)-j \sin (2 \omega)] \\ &=1-1.27[\cos (\pi)-j \sin (\pi)]+0.81[\cos (2 \pi)-j \sin (2 \pi)] \\ &=3.08+j 0 \\ |H(\omega)| &=\sqrt{a^{2}+b^{2}} \\ &=\sqrt{3.08^{2}+0^{2}}=3.08 \end{aligned} H(ω)∣ω=π​∣H(ω)∣​=1−1.27[cos(ω)−jsin(ω)]+0.81[cos(2ω)−jsin(2ω)]=1−1.27[cos(π)−jsin(π)]+0.81[cos(2π)−jsin(2π)]=3.08+j0=a2+b2 ​=3.082+02 ​=3.08​

$\frac{1}{2}$ Nyquist: $\pi /2$

$$\begin{array}{rl}{H(\omega )|}_{\omega =\pi /2}& =1-1.27[\cos (\omega )-j\sin (\omega )]+0.81[\cos (2\omega )-j\sin (2\omega )]\\ & =1-1.27[\cos (\pi /2)-j\sin (\pi /2)]+0.81[\cos (\pi )-j\sin (\pi )]\\ & =0.19+j1.27\\ |H(\omega )|& =\sqrt{a^2+b^2}\\ & =\sqrt{0.19^2+1.27^2}=1.28\end{array}$$

$\frac{1}{4}$ Nyquist: $\pi /4$

$$\begin{array}{rl}{H(\omega )|}_{\omega =\pi /2}& =1-1.27[\cos (\omega )-j\sin (\omega )]+0.81[\cos (2\omega )-j\sin (2\omega )]\\ & =1-1.27[\cos (\pi /4)-j\sin (\pi /4)]+0.81[\cos (\pi /2)-j\sin (\pi /2)]\\ & =0.11+j0.08\\ |H(\omega )|& =\sqrt{a^2+b^2}\\ & =\sqrt{0.11^2+0.08^2}=0.136\end{array}$$

可以看到直接计算的结果，与通过几何关系得到结果非常接近。

二阶反馈滤波器（2nd order feedback filter）

二阶反馈滤波器的差分方程为：
$$y(n)=a_{0}x(n)-b_{1}y(n-1)-b_{2}y(n-2)$$
进行 z变换：
$$Y(z)=a_{0}X(z)-b_{1}Y(z)z^{-1}-b_{2}Y(z)z^{-2}$$
得到转换方程：
$$H(z)=a_0\frac{1}{1+b_1{z}^{-1}+b_2{z}^{-2}}$$

此时 $H(z)$ 的分母是个二次方程，对它进行因式分解，可以有：
$$H(z)=a_0\frac{1}{\left(1-P_1{z}^{-1}\right)\left(1-P_2{z}^{-1}\right)}$$

与之前的分析类似，可以知道 $P_1$ 和 $P_2$ 为共轭复数：
$$\begin{array}{l}{P}_1={\mathrm{R}\mathrm{e}}^{j\Theta }\\ P_2={\mathrm{R}\mathrm{e}}^{-j\Theta }\end{array}$$

将其带入式子中，得到：
$$\begin{array}{l}\left(1-{\mathrm{Re}}^{j\Theta }{z}^{-1}\right)\left(1-{\mathrm{Re}}^{-j\Theta }{z}^{-1}\right)\\ =1-{\mathrm{Re}}^{j\Theta }{z}^{-1}-{\mathrm{Re}}^{-j\Theta }{z}^{-1}+R^2\left(e^{j\Theta }{e}^{-j\Theta }\right)z^{-2}\\ =1-2R\cos (\Theta )z^{-1}+R^2{z}^{-2}\end{array}$$
$$\begin{array}{rl}H(z)& =a_0\frac{1}{1+b_1{z}^{-1}+b_2{z}^{-2}}\\ & =a_0\frac{1}{\left(1-2R\cos (\Theta )z^{-1}+R^2{z}^{-2}\right)}\end{array}$$

通过对比可知：
$$\begin{array}{l}{b}_1=-2R\cos (\Theta )\\ b_2=R^2\end{array}$$

我们以 $a_0=1.0$，$b_1=-1.34$，$a_2=0.902$ 为例，可以得到：
$$\begin{array}{rl}{R}^2& =b_2=0.902\\ R& =\sqrt{0.902}=0.95\end{array}$$
$$\begin{array}{rl}-2\mathrm{Rcos}(\Theta )& =-1.34\\ 2(0.95)\cos (\Theta )& =1.34\\ \cos (\Theta )& =\frac{1.34}{2(0.95)}\\ \Theta & =\arccos (0.707)\\ \Theta & =45^{\circ }\end{array}$$

这样就得到了两个极点，它们位于复平面 半径为 0.95 ，角度为 $\pm 45^{\circ }$ 的位置，如下图：

估计频响

当估计一个以上零点的频响时，需要这么做：

1. 给定一个频率，在复平面坐标的单位圆外圈定位它
2. 画一条线，连接该点与所有零点
3. 将所有线段的长度的倒数相乘，得到的积将是频率的增幅

最后一条规则用公式来说明的话，应该是这样：
$${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=a_0\frac{1}{{\prod }_{i=1}^N{V}_i}$$
其中 $N$ 为滤波器的阶数（也就是零点的个数），$V_i$为与极点的距离长度。

以 $0$，$\pi$，$\pi /2$, $\pi /4$ 为例，用上述方法进行计算频响：

当 $\omega =0$ 时，与两个零点的距离为 0.71、0.71，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1/(0.71\ast 0.71)=1.98$

当 $\omega =\pi /4$ 时，与两个零点的距离为 0.05、1.41，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1/(0.05\ast 1.41)=14.2$

当 $\omega =\pi /2$ 时，与两个零点的距离为 0.71、1.84，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1/(0.71\ast 1.84)=0.76$

当 $\omega =\pi$ 时，与两个零点的距离为 1.8、1.8，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1/(1.8\ast 1.8)=0.31$

计算频谱响应

当 $a_0=1.0$，$b_1=-1.34$，$b_2=0.902$ 时
$$\begin{array}{rl}H(\mathrm{z})& =a_0\frac{1}{1+b_1{z}^{-1}+b_2{z}^{-2}}\\ & =\frac{1}{1-1.34z^{-1}+0.902z^{-2}}\end{array}$$

现在来计算它的频谱响应，公式太难敲了，这里给一个 0Hz 的结果，剩下的大家可以当做练习计算一下。正确的结果在下面的表格。

DC(0Hz)

$$\begin{array}{rl}H(z)& =\frac{1}{1-1.34+0.902}=\frac{1}{0.562+j0}\\ |H(\omega )|& =\frac{\mid \text{ numerator }\mid }{\mid \text{ denominator }\mid }=\frac{|1|}{|0.562+j0|}\\ & =\frac{1}{\sqrt{a^2+b^2}}\\ & =\frac{1}{\sqrt{0.562^2}}\\ & =1.78=5.00dB\end{array}$$

Frequency($\omega$)	$|H(\omega )|$
Nyquist($\pi$)	-10.2 dB
1/2 Nyquist($\pi /2$)	-2.56 dB
1/4 Nyquist($\pi /4$)	23.15 dB

搁架式滤波器（Shelving filter）

搁架式均衡器包含一个零点和一个极点，所以也被称为 " 1st order pole/zero filter "，它的差分方程为：
$$y(n)=a_{0}x(n)+{\mathrm{a}}_{1}x(n-1)-{\mathrm{b}}_{1}y(n-1)$$

进行 z变换：
$$\begin{array}{rl}Y(z)+b_{1}Y(z)z^{-1}& =a_{0}X(z)+a_{1}X(z)z^{-1}\\ Y(z)\left[1-b_1{z}^{-1}\right]& =X(z)\left[a_0+a_1{z}^{-1}\right]\end{array}$$

得到转换方程：
$$H(z)=\frac{Y(z)}{X(z)}=a_0\frac{1+{\alpha }_1{z}^{-1}}{1+b_1{z}^{-1}}$$
其中 ${\alpha }_1=\frac{a_1}{a_0}$

估计频响

对转换方程做一些变化：
$$\begin{array}{rl}H(\mathrm{z})& =a_0\frac{1+{\alpha }_1{z}^{-1}}{1+b_1{z}^{-1}}\\ & =a_0\frac{1+\frac{{\alpha }_1}{z}}{1+\frac{b_1}{z}}\end{array}$$
可以清楚的看到，它有一个零点$z=-\alpha$，和一个极点 $z=-b_1$

举个例子，当 $a_0=1.0$，$a_1=-0.92$，$b_1=-0.71$ 时，零点和极点的位置入下图：

当同时有零点和极点时，估计频谱响应和以前步骤没差别，只是需要同时处理极点和零点两种线段：

1. 给定一个频率，在复平面坐标的单位圆外圈定位它
2. 画线，连接该点与所有零点
3. 画线，连接该点与所有极点
4. 将所有与零点相连线段的长度相乘，得到零点幅度
5. 将所有与极点相连线段的长度的倒数相乘，得到极点幅度
6. 零点幅度除以极点幅度，就是频响

最后一条规则用公式来说明的话，应该是这样：
$${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=\frac{a_0{\prod }_{i=1}^N{U}_i}{{\prod }_{i=1}^N{V}_i}$$

以 $0$，$\pi$，$\pi /2$, $\pi /4$ 为例，用上述方法进行计算频响：

当 $\omega =0$ 时，与零点和极点的距离分别为 0.08、0.29，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=0.08\ast 1/0.29=0.27$

当 $\omega =\pi /4$ 时，与零点和极点的距离分别为 1.75、1.70，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1.75\ast 1.0/1.70=0.17$

当 $\omega =\pi /2$ 时，与零点和极点的距离分别为 1.40、1.30，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1.40\ast 1.0/1.30=1.07$

当 $\omega =\pi$ 时，与零点和极点的距离分别为 1.92、1.71，因此在 ${\left|H\left(e^{j\omega }\right)\right|}_{\omega }=1.92\ast 1.0/1.71=1.12$

计算频谱响应

现在来计算它的频谱响应，公式太难敲了，这里给一个 0Hz 的结果，剩下的大家可以当做练习计算一下。正确的结果在下面的表格。

DC(0Hz)

$$\begin{array}{rl}H(\omega )& =\frac{1-0.92[\cos (\omega )-j\sin (\omega )]}{1-0.71[\cos (\omega )-j\sin (\omega )]}\\ & =\frac{1-0.92[\cos (0)-j\sin (0)]}{1-0.71[\cos (0)-j\sin (0)]}=\frac{1-0.92[1-j0]}{1-0.71[1-j0]}\\ & =\frac{1-0.92}{1-0.71}=\frac{0.08+j0}{0.29+j0}\end{array}$$
$$\begin{array}{rl}|H(\omega )|& =\frac{|0.08+j0|}{|0.29+j0|}\\ & =\frac{\sqrt{a_{num}^2+b_{num}^2}}{\sqrt{a_{\text{denom}}^2+b_{\text{denom}}^2}}\\ & =\frac{\sqrt{0.08^2}}{\sqrt{0.29^2}}=0.276=-11.2dB\end{array}$$

Frequency($\omega$)	$|H(\omega )|$
Nyquist($\pi$)	1.00 dB
1/2 Nyquist($\pi /2$)	0.82 dB
1/4 Nyquist($\pi /4$)	0.375 dB

双二阶滤波器（Biquad Filter）

双二阶滤波器最多有两个零点和两个极点，它的差分方程为：
$$y(n)=a_{0}x(n)+a_{1}x(n-1)+a_{2}x(n-2)-b_{1}y(n-1)-b_{2}y(n-2)$$

进行 z变换：
$$\begin{array}{rl}Y(z)+b_{1}Y(z)z^{-1}+b_{2}Y(z)z^{-2}& =a_{0}X(z)+a_{1}X(z)z^{-1}+a_{2}X(z)z^{-2}\\ Y(z)\left[1+b_1{z}^{-1}+b_2{z}^{-2}\right]& =X(z)\left[a_0+a_1{z}^{-1}+a_2{z}^{-2}\right]\end{array}$$

得到转换方程为：
$$\begin{array}{rl}H(z)=\frac{Y(z)}{X(z)}& =\frac{a_0+a_1{z}^{-1}+a_2{z}^{-2}}{1+b_1{z}^{-1}+b_2{z}^{-2}}\\ & =a_0\frac{1+{\alpha }_1{z}^{-1}+{\alpha }_2{z}^{-2}}{1+b_1{z}^{-1}+b_2{z}^{-2}}\end{array}$$
其中${\alpha }_1=\frac{a_1}{a_0}{\alpha }_2=\frac{a_2}{a_0}$

估计频响

双二阶滤波器最多能从分子和分母中分别产生一对共轭零点和共轭极点。当然，我们也可以将某些系数设置为 0.0，例如当 $a_1=0.0$ 时，就会得到一阶前馈和二阶反馈滤波器的结合。

通过前面的分析，可以将双二阶滤波器的转换方程转换下面这种简单的形式：
$$H(z)=a_0\frac{1-2R_z\cos (\theta )z^{-1}+R_z^2{z}^{-2}}{1-2R_p\cos (\varphi )z^{-1}+R_p^2{z}^{-2}}$$

举个例子来说明如何估计频响，我们假设
$$\begin{array}{l}{a}_0=1.0\\ a_1=0.73\\ a_2=1.00\\ b_1=-0.78\\ b_2=0.88\end{array}$$

计算零点位置：
$$\begin{array}{rl}{a}_2& =R_z^2=1.00\\ R_z& =\sqrt{1.00}=1.00\\ & \\ a_1& =-2R\cos (\Theta )=0.73\\ 2(1.00)\cos (\Theta )& =-0.73\\ \cos (\Theta )& =-0.365\\ \Theta & =\arccos (-0.365)\\ \Theta & =111.1^o\end{array}$$

计算极点位置：
$$\begin{array}{rl}{b}_2& =R_p^2=0.88\\ R_p& =\sqrt{0.88}=0.94\\ & \\ 2(0.94)\cos (\varphi )& =0.78\\ \cos (\varphi )& =\frac{0.78}{2(0.94)}\\ \varphi & =\arccos (0.414)\\ \varphi & =65.5^{\circ }\end{array}$$

接下来步骤与 搁架式滤波器 中估计频响步骤一致，不再赘述。

双二阶滤波器的实现

双二阶滤波器有多种形式，其中
$$y(n)=a_{0}x(n)+a_{1}x(n-1)+a_{2}x(n-2)-b_{1}y(n-1)-b_{2}y(n-2)$$
被称为 “Direct Form”，它实现起来非常简单：

yn = a0*xn + a1*x(n−1) + a2*x(n−2)—b1*y(n−1)—b2*y(n−2)

// --- update state registers
x(n−2) = x(n−1)
x(n−1) = xn
y(n−2) = y(n−1) y(n−1) = yn

此外还有三种常用的形式，如下图：

所有这四种形式都会产生相同的差分方程。许多DSP工程师根据其计算效率、内存存储要求以及对系数中舍入误差的敏感性来选择需要的形式。

在 libaa 中你可以找到四种形式的所有实现。具体细节就不再赘述了。

总结

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