Python实现所有算法-音频过滤器.下（巴特沃斯）
上节简单的写了一下音频滤波器的定义和作用。而这篇文章将主要集中精力在巴特沃斯过滤器上，在末尾将会给出：使用 Butterworth 设计的二阶 IIR 滤波器。
另外，因为微信这个垃圾的公式排版，我也使用了：
来进行一个排版
$H(z)=\frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}+...+b_{k}z^{-k}}{a_{0}+a_{1}z^{-1}+a_{2}z^{-2}+...+a_{k}z^{-k}}$在构建N阶IIR滤波器的时候，使用的传递函数是这样的。
这个
然后会重写成这样：
y[n]={\frac{1}{a_{0}}}\left(\left(b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]+...+b_{k}x[n-k]\right)-\left(a_{1}y[n-1]+a_{2}y[n-2]+...+a_{k}y[n-k]\right)\right)
因为这个重写的公式太长了，截图就是这样的
接下来就是大家不爱看的东西了，我尽量选通俗易懂的话。先解释频繁出现的IIR滤波器是什么？
IIR（Infinite Impulse Response）——无限脉冲响应滤波器，是数字信号处理的基本组成部分。
1、IIR滤波器是什么
IIR滤波器是用于数字信号处理(DSP)应用的两种主要数字滤波器之一(另一种是FIR)。“IIR”的意思是“无限脉冲响应”。
2、 IIR为什么脉冲响应是“无限的”?
脉冲响应是“无限”的，因为滤波器中有反馈;如果你输入一个脉冲(一个“1”样本后面跟着多个“0”样本)，理论上就会输出无限个非零值。
3、 IIR过滤器的替代方案是什么?
DSP滤波器也可以是“有限脉冲响应”(FIR)。FIR滤波器不使用反馈，所以对于N个系数的FIR滤波器，输入N个脉冲响应的样本后输出总是为零。
4、IIR滤波器(与FIR滤波器相比)的优点是什么?
与类似的FIR滤波器相比，IIR滤波器可以用更少的内存和计算来实现给定的滤波特性。
5、IIR滤波器(与FIR滤波器相比)的缺点是什么?
（1）它们更容易受到有限长度算法问题的影响，比如计算产生的噪声和极限环。(这是反馈造成的直接结果:当输出没有被完美地计算出来并得到反馈时，不完美可能会加剧。）
（2）它们使用定点算法更难(更慢)实现.
（3）对于多速率(抽取和插值)应用，它们没有FIR滤波器的计算优势。
再说一个，什么叫线性时不变系统？
一、线性
通信系统中的线性不再是数学中坐标轴上的直线，也不是所有的直线都符合线性特征，通信系统中的线性要满足一个条件。
为什么线性系统怎么重要。如果已知一个系统是线性系统，一旦有了任何特定输入的输出，就可以在不进行实际测试的情况下计算出所有其他输出。这就是线性系统的优点。也就是系统是可控的
二、时不变系统
一句话概括时不变系统：如果你的系统是时不变系统，如果将相同的输入输入到系统中，无论何时将输入输入到系统中，都将获得相同的输出。
那么我们还假设你的系统是时不变系统，你有一个确定的输入，你只需在知道某个时间下这个输入的所有特征。如果你明天还输入相同的输入，那么你的不用重新掌握它的特征。
相对的如果你的系统是时变系统：
投入了一定的投入，花了大量的时间和精力调查了投入的所有特征。但是无法重用结果，因为当重用它时，时间不同，结果也会不同。就像。。。你有一台电脑。每次运行完全相同的程序时，都会得到不同的结果。
唯一可以重用结果的情况是进一步努力找出规则，以显示输出是如何随时间变化的。不幸的是，你不能保证你能找到规则。即使你很幸运找到了规则，你的系统模型也会非常复杂。
三、总结
线性系统，你知道一个确定的结果，你就知道其他所有的结果，
时不变系统，你的特定输入，其输出的结果不会随着时间而变化。
这个IIR二阶的传递函数是从wiki里面拿到的
因为这个IIR是给下一级使用的，这里设计成一个类
类的初始化
接下来是使用一个现成的包来验证一下结果
返回的是分子和分母的多项式
计算a，b
设置系数用于 IIR 滤波器。这些都应该是 size order + 1，a_0 可以省略，它将使用 1.0 作为默认值。
这就是这样，下面的都是对参数的约束
a，b就是系数
计算f(n)，对着公式编程序
这个就是两项计算的公式
按照这个写的
关于滤波器还有若干的问题需要说明：一阶滤波器，思路就是把一个连续的滤波器形式，通过离散化的方式，转换成差分方程。从wiki或者文章里面拿到的滤波器公式，通常是用传递函数表达的，这是S域下的表达形式，是连续的，这种我们称之为模拟滤波器。模拟滤波器传递函数，目的是用来设计滤波电路，针对的是连续时间的模拟信号，组成元器件是电阻，电容，电感。而数字滤波器实现方法是把滤波器所要完成的运算编成程序并让计算机执行,也就是采用在代码的形式。它面对的是离散时间的数字信号，是把输入序列通过一定的运算变换成输出序列。有没有办法能把连续的模拟滤波器变成离散的数字滤波器？
其中最常使用的一种叫做双线性变换：
把这个公式带入传递函数就可以得到一个z域的差分方程了。
以上变换这段参考:
后面截至频率什么的没有写
但是知道接下来应该看的是：自动控制原理。
我们在BW滤波器里面将要实现这些算法
所有滤波器传递函数均源自模拟原型，并已使用双线性变换 (BLT) 进行数字化。BLT 频率扭曲已被考虑到重要的频率重定位（这是使用 BLT 时必需的正常“预扭曲”）和带宽重新调整（因为使用 BLT 从模拟映射到数字时带宽被压缩）。这个就是上面文章里面说的：
出现这个问题，需要使用数学工具矫正
模拟截止频率与数字截止频率的关系
证明自己看吧，我也看不懂了，越看越觉得自己不该念这个书。。。这是TM的什么人间疾苦。
不过看了下。。。有欧拉公式就行
先定义双线性的函数，在上面写了
经过一次变换
最后才写成这个
省去一段推导，给出结果：
低通滤波器的结果
参数是，频率，采样率，Q值
第一个的滤波器的算法设计
在最后给出一个在双线性变换下使用补偿频率扭曲的补偿算法推导：
首先是归一化的操作
使用这个三角函数
还有这个
对公式进行重构，构个头就是重新带进去
三个
还有因子也重新写
对分子和分母中的所有项都是通用的，可以因式分解，因此在上面的替换中被忽略，所以有：
又重写了
此外，所有项，分子和分母，都可以乘以一个共同的：
好大个系数
最后就是简化出来的双二阶系数公式
啧，万物之源是数学，还不滚去学习。
在最终的演示中
这个错误
说是3.8才有
pip install typing_extensions
from typing_extensions import Protocol
可以写成这样
话说有朝一日我还能拍出这种照片
是不是有点古城的味道了嗷！
考虑到这个算法的复杂性，这里将代码写在了一起，可以直接调用：
from __future__ import annotations
# from audio_filters.iir_filter import IIRFilter
import numpy as np
import matplotlib.pyplot as plt
from typing_extensions import Protocol
from math import pi
from math import cos, sin, sqrt, tau
class IIRFilter:
def __init__(self, order: int) -> None:
self.order = order
# a_{0} ... a_{k}
self.a_coeffs = [1.0] + [0.0] * order
# b_{0} ... b_{k}
self.b_coeffs = [1.0] + [0.0] * order
# x[n-1] ... x[n-k]
self.input_history = [0.0] * self.order
# y[n-1] ... y[n-k]
self.output_history = [0.0] * self.order
def set_coefficients(self, a_coeffs: list[float], b_coeffs: list[float]) -> None:
if len(a_coeffs) < self.order:
a_coeffs = [1.0] + a_coeffs
if len(a_coeffs) != self.order + 1:
raise ValueError(
f"预期 a_coeffs to 有 {self.order + 1} elements for {self.order}"
f"-order filter, got {len(a_coeffs)}"
)
if len(b_coeffs) != self.order + 1:
raise ValueError(
f"Expected b_coeffs to have {self.order + 1} elements for {self.order}"
f"-order filter, got {len(a_coeffs)}"
)
self.a_coeffs = a_coeffs
self.b_coeffs = b_coeffs
def process(self, sample: float) -> float:
result = 0.0
# 从索引 1 开始，最后执行索引 0
for i in range(1, self.order + 1):
result += (self.b_coeffs[i] * self.input_history[i-1]-self.a_coeffs[i]*self.output_history[i-1]
)
result = (result + self.b_coeffs[0] * sample) / self.a_coeffs[0]
self.input_history[1:] = self.input_history[:-1]
self.output_history[1:] = self.output_history[:-1]
self.input_history[0] = sample
self.output_history[0] = result
return result
def make_lowpass(
frequency: int, samplerate: int, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
b0 = (1 - _cos) / 2
b1 = 1 - _cos
a0 = 1 + alpha
a1 = -2 * _cos
a2 = 1 - alpha
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b0])
return filt
def make_highpass(
frequency: int, samplerate: int, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
b0 = (1 + _cos) / 2
b1 = -1 - _cos
a0 = 1 + alpha
a1 = -2 * _cos
a2 = 1 - alpha
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b0])
return filt
def make_bandpass(
frequency: int, samplerate: int, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
"""
创建带通滤波器
>>> filter = make_bandpass(1000, 48000)
>>> filter.a_coeffs + filter.b_coeffs # doctest: +NORMALIZE_WHITESPACE
[1.0922959556412573, -1.9828897227476208, 0.9077040443587427, 0.06526309611002579,
0, -0.06526309611002579]
"""
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
b0 = _sin / 2
b1 = 0
b2 = -b0
a0 = 1 + alpha
a1 = -2 * _cos
a2 = 1 - alpha
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
return filt
def make_allpass(
frequency: int, samplerate: int, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
b0 = 1 - alpha
b1 = -2 * _cos
b2 = 1 + alpha
filt = IIRFilter(2)
filt.set_coefficients([b2, b1, b0], [b0, b1, b2])
return filt
def make_peak(
frequency: int, samplerate: int, gain_db: float, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
big_a = 10 ** (gain_db / 40)
b0 = 1 + alpha * big_a
b1 = -2 * _cos
b2 = 1 - alpha * big_a
a0 = 1 + alpha / big_a
a1 = -2 * _cos
a2 = 1 - alpha / big_a
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
return filt
def make_lowshelf(
frequency: int, samplerate: int, gain_db: float, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
big_a = 10 ** (gain_db / 40)
pmc = (big_a + 1) - (big_a - 1) * _cos
ppmc = (big_a + 1) + (big_a - 1) * _cos
mpc = (big_a - 1) - (big_a + 1) * _cos
pmpc = (big_a - 1) + (big_a + 1) * _cos
aa2 = 2 * sqrt(big_a) * alpha
b0 = big_a * (pmc + aa2)
b1 = 2 * big_a * mpc
b2 = big_a * (pmc - aa2)
a0 = ppmc + aa2
a1 = -2 * pmpc
a2 = ppmc - aa2
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
return filt
def make_highshelf(
frequency: int, samplerate: int, gain_db: float, q_factor: float = 1 / sqrt(2)
) -> IIRFilter:
w0 = tau * frequency / samplerate
_sin = sin(w0)
_cos = cos(w0)
alpha = _sin / (2 * q_factor)
big_a = 10 ** (gain_db / 40)
pmc = (big_a + 1) - (big_a - 1) * _cos
ppmc = (big_a + 1) + (big_a - 1) * _cos
mpc = (big_a - 1) - (big_a + 1) * _cos
pmpc = (big_a - 1) + (big_a + 1) * _cos
aa2 = 2 * sqrt(big_a) * alpha
b0 = big_a * (ppmc + aa2)
b1 = -2 * big_a * pmpc
b2 = big_a * (ppmc - aa2)
a0 = pmc + aa2
a1 = 2 * mpc
a2 = pmc - aa2
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
return filt
class FilterType(Protocol):
def process(self, sample: float) -> float:
return 0.0
def get_bounds(
fft_results: np.ndarray, samplerate: int
) -> tuple[int | float, int | float]:
lowest = min([-20, np.min(fft_results[1: samplerate // 2 - 1])])
highest = max([20, np.max(fft_results[1: samplerate // 2 - 1])])
return lowest, highest
def show_frequency_response(filter: FilterType, samplerate: int) -> None:
size = 512
inputs = [1] + [0] * (size - 1)
outputs = [filter.process(item) for item in inputs]
filler = [0] * (samplerate - size) # zero-padding
outputs += filler
fft_out = np.abs(np.fft.fft(outputs))
fft_db = 20 * np.log10(fft_out)
# Frequencies on log scale from 24 to nyquist frequency
plt.xlim(24, samplerate / 2 - 1)
plt.xlabel("Frequency (Hz)")
plt.xscale("log")
# Display within reasonable bounds
bounds = get_bounds(fft_db, samplerate)
plt.ylim(max([-80, bounds[0]]), min([80, bounds[1]]))
plt.ylabel("Gain (dB)")
plt.plot(fft_db)
plt.show()
def show_phase_response(filter: FilterType, samplerate: int) -> None:
size = 512
inputs = [1] + [0] * (size - 1)
outputs = [filter.process(item) for item in inputs]
filler = [0] * (samplerate - size)
outputs += filler
fft_out = np.angle(np.fft.fft(outputs))
plt.xlim(24, samplerate / 2 - 1)
plt.xlabel("Frequency (Hz)")
plt.xscale("log")
plt.ylim(-2 * pi, 2 * pi)
plt.ylabel("Phase shift (Radians)")
plt.plot(np.unwrap(fft_out, -2 * pi))
plt.show()
filt = IIRFilter(4)
show_phase_response(filt, 48000)
应该是可以直接出现这样的图
https://stackoverflow.com/questions/54274630/can-not-import-protocol-from-typinghttps://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.htmlhttps://blog.csdn.net/abc123mma/article/details/120251384https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.htmlhttps://blog.csdn.net/SHC2012377/article/details/120910496https://en.wikipedia.org/wiki/Digital_biquad_filterhttps://zhuanlan.zhihu.com/p/357619650https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html腾讯云开发者
