系列文章目录

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

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

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

文章目录

系列文章目录
前言
一、模拟信号滤波器到数字信号滤波器的转换
二、S域与Z域
三、模拟信号滤波器
四、双线性变换
五、滤波器设计的步骤
总结
参考

前言

在【音频处理】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=σ+jωs = \sigma + j\omega s=σ+jω
σ∈R,ω∈R\sigma \in R, \omega \in Rσ∈R,ω∈R，它的图像长这样：

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

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

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

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

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

三、模拟信号滤波器

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

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

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

四、双线性变换

现在我们已经知道了 z 域是 s 域的扩展，它们之间通过某种神奇的映射联系在一起。在如何理解离散傅里叶变换及Z变换
详细介绍了 z 与 s 之间的关系，具体来说：
z=esTsz = e^{sT_s} z=esTs
其中 Ts=1/fsT_s = 1/f_sTs=1/fs, fs 为采样率。我们对这个公式进行变换：
ln⁡(z)=ln⁡(esT)ln⁡(z)=sTorsT=ln⁡(z)\begin{array}{l} \ln (z)=\ln \left(e^{s T}\right) \\ \ln (z)=s T \\ \text{or} \\ s T=\ln (z) \\ \end{array} ln(z)=ln(esT)ln(z)=sTorsT=ln(z)
其中 ln⁡(z)\ln(z)ln(z) 的泰勒展开为：
sT=2[z−1z+1+13(z−1z+1)3−15(z−1z+1)5+17(z−1z+1)7−…]s T=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}-\ldots\right] sT=2[z+1z−1+31(z+1z−1)3−51(z+1z−1)5+71(z+1z−1)7−…]
我们只取第一项，对结果进行近似，得到：
s=2Tz−1z+1s=\frac{2}{T} \frac{z-1}{z+1} s=T2z+1z−1
这就是双线性变化了。同时可以变换得到：
z=1+sT/21−sT/2z = \frac{1+sT/2}{1-sT/2} z=1−sT/21+sT/2
当 s=0s = 0s=0 时，即 0 频率，此时 z=1z = 1z=1
当 s→∞s \rightarrow \inftys→∞ 时，即无穷大频率，此时 z=−1z = -1z=−1，对应频率为 1/fs1/f_s1/fs
也就是说，双线性变换的本质是将 s 域中无穷大的频率映射到 z 域中的 1/fs1/f_s1/fs。画个图，长这样：

也就是说，在双线性变换下，s 域中的频率与 z 域中频率不同，发生了一定的弯曲。这也意味着截止频率在 s 域和 z 域是不一样的，我们需要找到这种截止频率的对应关系：
z=esT=ejωdT\begin{aligned} z &= e^{sT} = e^{j\omega_d T} \\ \end{aligned} z=esT=ejωdT
s=2Tz−1z+1=2TejωdT−1ejωdT+1令ωdT=θs=2Tejθ−1ejθ+1=2Tejθ/2(ejθ/2−e−jθ/2)ejθ/2(ejθ/2+e−jθ/2)=2T(ejθ/2−e−jθ/2)(ejθ/2+e−jθ/2)利用欧拉公式ejx=cos⁡x+jsin⁡xs=2Tjsin⁡θ/2cos⁡θ/2=2Tjtan⁡θ/2\begin{aligned} s &=\frac{2}{T} \frac{z-1}{z+1} \\ &= \frac{2}{T} \frac{e^{j\omega_dT} - 1}{e^{j\omega_dT} + 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{aligned} s令ωdT=θs利用欧拉公式s=T2z+1z−1=T2ejωdT+1ejωdT−1=T2ejθ+1ejθ−1=T2ejθ/2(ejθ/2+e−jθ/2)ejθ/2(ejθ/2−e−jθ/2)=T2(ejθ/2+e−jθ/2)(ejθ/2−e−jθ/2)ejx=cosx+jsinx=T2cosθ/2jsinθ/2=T2jtanθ/2
又因为 s=jωas=j\omega_as=jωa(只考虑稳态时)，因此有：
jωa=2Tjtan⁡θ/2=2Tjtan⁡(ωdT2)\begin{aligned} j\omega_a &= \frac{2}{T} j\tan{\theta/2} \\ &= \frac{2}{T} j\tan(\frac{\omega_d T}{2}) \end{aligned} jωa=T2jtanθ/2=T2jtan(2ωdT)
可得，s 域与 z 域截止频率之间的关系为：
ωa=2Ttan⁡(ωdT2)\omega_a = \frac{2}{T} \tan(\frac{\omega_d T}{2}) ωa=T2tan(2ωdT)
其中 ωa\omega_aωa 表示模拟频率，ωd\omega_dωd 表示映射在 z 域中的频率，T=1fsT = \frac{1}{f_s}T=fs1

五、滤波器设计的步骤

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

选定归一化的低通滤波器，即H(s)=1s+1H(s)=\frac{1}{s+1}H(s)=s+11
确定数字滤波器的截止频率 ωd\omega_dωd
通过 ωa=2Ttan⁡(ωdT2)\omega_a = \frac{2}{T} \tan(\frac{\omega_d T}{2})ωa=T2tan(2ωdT) 计算模拟频率
采取合适的变换得到 H′(s)H'(s)H′(s)，例如低通滤波器变换为 s=sωas=\frac{s}{\omega_a}s=ωas
在 H′(s)H'(s)H′(s) 中采取 s=2Tz−1z+1s=\frac{2}{T} \frac{z-1}{z+1}s=T2z+1z−1 进行变化，得到数字滤波器 H(z)H(z)H(z)。
用代数方法处理数字传递函数H(z)，使其变成熟悉的形式，以便你能够识别系数，这往往是最困难的部分。

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

选定归一化的低通滤波器，即H(s)=1s+1H(s)=\frac{1}{s+1}H(s)=s+11
截止频率 ωd=2π∗1000=6283.2\omega_d=2\pi*1000 = 6283.2ωd=2π∗1000=6283.2rad/sec
ωa=2Ttan⁡(ωdT2)=6293.85\omega_a = \frac{2}{T} \tan(\frac{\omega_d T}{2}) = 6293.85ωa=T2tan(2ωdT)=6293.85 rad/sec
对 H(x)H(x)H(x) 中 s 进行替换，s→sωas \rightarrow \frac{s}{\omega_a}s→ωas，得到 H′(s)=ωas+ωaH'(s) = \frac{\omega_a}{s + \omega_a}H′(s)=s+ωaωa
用 s=2Tz−1z+1s=\frac{2}{T} \frac{z-1}{z+1}s=T2z+1z−1 对 H′(s)H'(s)H′(s) 进行替换，得到 H′(s)=0.0667(z+1)z−0.8667H'(s) = \frac{0.0667(z + 1)}{z - 0.8667}H′(s)=z−0.86670.0667(z+1)
整理成熟悉的形式 H(z)=0.0667+0.0667z−11−0.8667z−1=a0+a1z−11+b1z−1H(z)=\frac{0.0667+0.0667 z^{-1}}{1-0.8667 z^{-1}}=\frac{a_{0}+a_{1} z^{-1}}{1+b_{1} z^{-1}}H(z)=1−0.8667z−10.0667+0.0667z−1=1+b1z−1a0+a1z−1

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

总结

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

参考

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

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

滤波器设计（二）模拟到数字相关推荐

双线性变换 matlab,matlab和双线性变换的滤波器设计.doc
matlab和双线性变换的滤波器设计.doc 武汉理工大学MATLAB课程设计报告书题目MATLAB课程设计基于MATLAB和双线性变换的滤波器设计初始条件MATLAB仿真软件数字信号处理与图像处理基 ...
【数字信号处理2】IIR 滤波器设计
一.实验目的 1.掌握冲激响应法和双线性变换法设计IIR滤波器的原理及具体设计方法,熟悉用双线性设计法设计低通.带通和高通IIR数字滤波器的计算机程序: 2.熟悉模拟Butterworth滤波器的设计 ...
如何使用双线性变换法将模拟电路滤波器设计成为数字滤波器？
信号处理(数字模拟信号) 1.1 双线性变换设计递归滤波器设计滤波器1 试着写出双线性变换法设计IIR数字高通滤波器的主要步骤将数字高通的频率指标转换为模拟高通的频率指标(其中将高通截止频率通过预 ...
数字信号处理6：IIR滤波器设计
IIR滤波器设计文章目录 IIR滤波器设计 1. 简介 2. 设计步骤简明 3. 拉普拉斯变换和Z变换 3.1 拉普拉斯变换 3.2 Z变换 4. 双线性变换法 4.1 模拟域与数字域的映射 4.2 ...
FIR数字滤波器的FPGA实现(二)-串行FIR滤波器设计(1)
(二)FIR数字滤波器的FPGA实现-串行FIR滤波器设计文章目录 (二)FIR数字滤波器的FPGA实现-串行FIR滤波器设计 0 串行FIR滤波器基本原理 1 基于移位寄存器的串行 FIR 滤波器 ...
如何利用CIC滤波器、CIC补偿滤波器和半带滤波器设计一个高频数字抽取滤波器
设计了采样频率为640 MHz.过采样率为64的高频数字抽取滤波器.该数字抽取滤波器由CIC(Cascaded Integrator Comb)滤波器(降16倍).CIC补偿滤波器(降2倍)和半带滤波 ...
数字信号处理5：FIR滤波器设计
文章目录 1. 滤波器初识 2. 最直观的滤波方式:频域滤波 3. 傅里叶变换中的加窗 4. FIR滤波器设计 5. 总结之前的一系列博客中,详细分解了从卷积到FFT的相关知识,不过那些属于理论,是 ...
FIR数字滤波器的FPGA实现(二)-串行FIR滤波器设计(2)
(二)FIR数字滤波器的FPGA实现-串行FIR滤波器设计文章目录 (二)FIR数字滤波器的FPGA实现-串行FIR滤波器设计 0 串行FIR滤波器基本原理 1 基于移位寄存器的串行 FIR 滤波器 ...
音频D类功放LC滤波器设计（二）
上一节我们分析了D类功放的频谱,这一节就来具体看看滤波器该怎么设计. 截止频率的确定首先,要设计滤波器,自然需要知道截止频率设计到多少比较合适. 我们上一节分析了频谱,可以知道,频谱里面除了包含音频 ...

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