Python實現所有算法-音頻過濾器.下（巴特沃斯）

$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}}$

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)

pip install typing_extensions
from typing_extensions import Protocol

from __future__ import annotations# from audio_filters.iir_filter import IIRFilterimport numpy as npimport matplotlib.pyplot as pltfrom typing_extensions import Protocolfrom math import pifrom 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-typing

https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html

https://blog.csdn.net/abc123mma/article/details/120251384

https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html

https://blog.csdn.net/SHC2012377/article/details/120910496

https://en.wikipedia.org/wiki/Digital_biquad_filter

https://zhuanlan.zhihu.com/p/357619650

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html