diff --git a/audio_filters/README.md b/audio_filters/README.md
index 4419bd8bdbf9..a28f2d6a9949 100644
--- a/audio_filters/README.md
+++ b/audio_filters/README.md
@@ -1,9 +1,58 @@
-# Audio Filter
+# Audio Filters
-Audio filters work on the frequency of an audio signal to attenuate unwanted frequency and amplify wanted ones.
-They are used within anything related to sound, whether it is radio communication or a hi-fi system.
+Audio filters work on the frequency of an audio signal to attenuate unwanted
+frequencies and amplify wanted ones. They are used within anything related to
+sound, whether it is radio communication or a hi-fi system. If you have ever
+turned up the bass or cut the treble on a stereo, tuned a radio to a station, or
+removed the background hum from a recording, you have used an audio filter.
+
+Curious to learn more? These are great starting points:
*
*
*
*
+*
+
+## What's in this directory
+
+| File | Description |
+| ---- | ----------- |
+| [`iir_filter.py`](iir_filter.py) | A generic N-order [Infinite Impulse Response (IIR)](https://en.wikipedia.org/wiki/Infinite_impulse_response) filter. This is the engine every filter below runs on: give it a set of coefficients and it processes a stream of samples one at a time. |
+| [`butterworth_filter.py`](butterworth_filter.py) | A collection of second-order [Butterworth](https://en.wikipedia.org/wiki/Butterworth_filter) / biquad filter designs from the RBJ Audio EQ Cookbook. Each function returns a ready-to-use `IIRFilter`. |
+| [`equal_loudness_filter.py`](equal_loudness_filter.py) | An [equal-loudness](https://en.wikipedia.org/wiki/Equal-loudness_contour) filter that compensates for the human ear's non-linear response to sound by cascading a Yule-Walker filter and a Butterworth high-pass filter. Includes a dependency-free `yulewalk` implementation. |
+| [`show_response.py`](show_response.py) | Helpers to plot the [magnitude and phase response](https://en.wikipedia.org/wiki/Frequency_response) of any filter so you can *see* what it does. |
+| [`loudness_curve.json`](loudness_curve.json) | The Robinson-Dadson equal-loudness contour data used by the equal-loudness filter. |
+
+## Filter designs in `butterworth_filter.py`
+
+| Function | Effect |
+| -------- | ------ |
+| `make_lowpass` | Passes frequencies below the cutoff, attenuates those above it. |
+| `make_highpass` | Passes frequencies above the cutoff, attenuates those below it. |
+| `make_bandpass` | Passes a band of frequencies around the center (constant skirt gain). |
+| `make_bandpass_peak` | Passes a band of frequencies around the center (constant 0 dB peak gain). |
+| `make_notch` | Rejects a narrow band around the center — great for removing mains hum. |
+| `make_allpass` | Passes all frequencies but changes their phase relationship. |
+| `make_peak` | Boosts or cuts a band around the center by a given gain (parametric EQ). |
+| `make_lowshelf` | Boosts or cuts everything below the cutoff. |
+| `make_highshelf` | Boosts or cuts everything above the cutoff. |
+
+## Try it out
+
+```python
+from audio_filters.butterworth_filter import make_lowpass
+from audio_filters.show_response import show_frequency_response
+
+# A 5 kHz low-pass filter for CD-quality audio (44.1 kHz sample rate)
+filt = make_lowpass(5000, 44100)
+
+# Process samples one at a time...
+filtered = [filt.process(sample) for sample in my_audio_samples]
+
+# ...or visualise what the filter does to the spectrum:
+show_frequency_response(make_lowpass(5000, 44100), 44100)
+```
+
+Every module has runnable doctests — read them for concrete, copy-pasteable
+examples of each filter in action.
diff --git a/audio_filters/butterworth_filter.py b/audio_filters/butterworth_filter.py
index 4e6ea1b18fb4..686084714b26 100644
--- a/audio_filters/butterworth_filter.py
+++ b/audio_filters/butterworth_filter.py
@@ -7,6 +7,16 @@
Code based on https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
Alternatively you can use scipy.signal.butter, which should yield the same results.
+
+https://en.wikipedia.org/wiki/Butterworth_filter
+
+Notation used throughout this module (from the RBJ Audio EQ Cookbook):
+ w0 -- normalised angular frequency, ``2 * pi * frequency / samplerate``
+ alpha -- bandwidth parameter, ``sin(w0) / (2 * q_factor)``
+ b0..b2 -- feed-forward (numerator) coefficients of the biquad
+ a0..a2 -- feed-back (denominator) coefficients of the biquad
+The a/b coefficient names match ``IIRFilter.set_coefficients`` and the standard
+biquad transfer function, so they are kept consistent across every filter here.
"""
@@ -232,3 +242,81 @@ def make_highshelf(
filt = IIRFilter(2)
filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
return filt
+
+
+def make_notch(
+ frequency: int,
+ samplerate: int,
+ q_factor: float = 1 / sqrt(2),
+) -> IIRFilter:
+ """
+ Creates a notch (band-reject) filter that strongly attenuates a narrow band
+ of frequencies around ``frequency`` while leaving the rest of the spectrum
+ unchanged. It is the complement of the band-pass filter and is commonly used
+ to remove a single tone such as 50/60 Hz mains hum.
+
+ https://en.wikipedia.org/wiki/Band-stop_filter
+
+ >>> filter = make_notch(1000, 48000)
+ >>> filter.a_coeffs + filter.b_coeffs # doctest: +NORMALIZE_WHITESPACE
+ [1.0922959556412573, -1.9828897227476208, 0.9077040443587427, 1.0,
+ -1.9828897227476208, 1.0]
+ """
+ w0 = tau * frequency / samplerate # centre frequency, in radians/sample
+ _sin = sin(w0)
+ _cos = cos(w0)
+ alpha = _sin / (2 * q_factor) # controls how narrow the rejected band is
+
+ # Feed-forward: a pair of zeros placed exactly on the notch frequency, so
+ # that frequency is fully cancelled while the rest of the spectrum passes.
+ b0 = 1.0
+ b1 = -2 * _cos
+ b2 = 1.0
+
+ # Feed-back: matching poles just inside the unit circle keep the notch
+ # narrow and the surrounding gain flat.
+ 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_bandpass_peak(
+ frequency: int,
+ samplerate: int,
+ q_factor: float = 1 / sqrt(2),
+) -> IIRFilter:
+ """
+ Creates a band-pass filter with constant 0 dB peak gain.
+
+ Unlike :func:`make_bandpass`, whose skirt (edge) gain is held constant so the
+ peak gain grows with ``q_factor``, this variant normalises the response so
+ the peak always reaches 0 dB regardless of the chosen ``q_factor``. Both
+ forms come from the RBJ Audio EQ Cookbook.
+
+ https://en.wikipedia.org/wiki/Band-pass_filter
+
+ >>> filter = make_bandpass_peak(1000, 48000)
+ >>> filter.a_coeffs + filter.b_coeffs # doctest: +NORMALIZE_WHITESPACE
+ [1.0922959556412573, -1.9828897227476208, 0.9077040443587427,
+ 0.09229595564125725, 0, -0.09229595564125725]
+ """
+ w0 = tau * frequency / samplerate
+ _sin = sin(w0)
+ _cos = cos(w0)
+ alpha = _sin / (2 * q_factor)
+
+ b0 = alpha
+ b1 = 0
+ b2 = -alpha
+
+ a0 = 1 + alpha
+ a1 = -2 * _cos
+ a2 = 1 - alpha
+
+ filt = IIRFilter(2)
+ filt.set_coefficients([a0, a1, a2], [b0, b1, b2])
+ return filt
diff --git a/audio_filters/equal_loudness_filter.py b/audio_filters/equal_loudness_filter.py
new file mode 100644
index 000000000000..387772aa1f50
--- /dev/null
+++ b/audio_filters/equal_loudness_filter.py
@@ -0,0 +1,198 @@
+from __future__ import annotations
+
+from json import loads
+from pathlib import Path
+
+import numpy as np
+from scipy.linalg import toeplitz
+from scipy.signal import lfilter, unit_impulse
+
+from audio_filters.butterworth_filter import make_highpass
+from audio_filters.iir_filter import IIRFilter
+
+data = loads((Path(__file__).resolve().parent / "loudness_curve.json").read_text())
+
+
+def _polystab(poly: np.ndarray) -> np.ndarray:
+ """
+ Stabilize a polynomial by reflecting any roots that lie outside the unit
+ circle back inside it. This keeps the resulting IIR filter stable without
+ changing its magnitude response.
+
+ https://en.wikipedia.org/wiki/Minimum_phase
+
+ >>> np.round(_polystab(np.array([1.0, 2.0, 1.0])), 6)
+ array([1., 2., 1.])
+ >>> np.round(_polystab(np.array([1.0, 2.0, 1.01])), 6)
+ array([1. , 1.980198, 0.990099])
+ """
+ if poly.size <= 1:
+ return poly
+ roots = np.roots(poly)
+ nonzero = np.where(roots != 0)[0]
+ outside = 0.5 * (np.sign(np.abs(roots[nonzero]) - 1) + 1)
+ roots[nonzero] = (1 - outside) * roots[nonzero] + outside / np.conj(roots[nonzero])
+ stabilized = np.poly(roots)
+ if not np.imag(poly).any():
+ stabilized = np.real(stabilized)
+ return stabilized
+
+
+def _numerator(
+ impulse_response: np.ndarray, denominator: np.ndarray, numerator_order: int
+) -> np.ndarray:
+ """
+ Least-squares estimate of the numerator polynomial of a transfer function
+ given its impulse response and (already known) denominator polynomial.
+
+ >>> num = _numerator(np.array([1.0, 0.0, 0.0]), np.array([1.0, 0.0, 0.0]), 1)
+ >>> np.round(num, 6)
+ array([1., 0.])
+ """
+ length = impulse_response.size
+ impulse = lfilter([1.0], denominator.ravel(), unit_impulse(length))
+ toep = toeplitz(impulse, unit_impulse(numerator_order + 1))
+ return np.linalg.lstsq(toep.conj(), impulse_response.ravel().conj(), rcond=None)[
+ 0
+ ].conj()
+
+
+def yulewalk(
+ order: int, frequencies: np.ndarray, magnitudes: np.ndarray, npt: int = 512
+) -> tuple[np.ndarray, np.ndarray]:
+ """
+ Design a recursive (IIR) digital filter that approximates an arbitrary
+ frequency response using the modified Yule-Walker method. This is a
+ dependency-free re-implementation of MATLAB/Octave's ``yulewalk`` so that
+ the equal-loudness filter below no longer relies on a third-party package.
+
+ https://en.wikipedia.org/wiki/Autoregressive_model#Yule%E2%80%93Walker_equations
+
+ :param order: order of the filter to design
+ :param frequencies: sample points on ``[0, 1]`` where 1 is the Nyquist
+ frequency, in increasing order and starting at 0
+ :param magnitudes: desired (linear) magnitude at each point in ``frequencies``
+ :param npt: number of points used to estimate the frequency response
+ :return: ``(a_coeffs, b_coeffs)``, the denominator and numerator polynomials
+
+ >>> a, b = yulewalk(4, np.array([0.0, 0.5, 1.0]), np.array([1.0, 0.5, 0.0]))
+ >>> len(a), len(b)
+ (5, 5)
+ >>> bool(np.all(np.abs(np.roots(a)) < 1)) # the designed filter is stable
+ True
+
+ Mismatched inputs and non-increasing frequencies are rejected:
+
+ >>> yulewalk(4, np.array([0.0, 1.0]), np.array([1.0]))
+ Traceback (most recent call last):
+ ...
+ ValueError: frequencies and magnitudes must have the same length
+ >>> yulewalk(4, np.array([0.0, 1.0, 0.5]), np.array([1.0, 0.5, 0.0]))
+ Traceback (most recent call last):
+ ...
+ ValueError: frequencies must be in increasing order
+ """
+ frequencies = np.asarray(frequencies, dtype=float).ravel()
+ magnitudes = np.asarray(magnitudes, dtype=float).ravel()
+ if frequencies.size != magnitudes.size:
+ msg = "frequencies and magnitudes must have the same length"
+ raise ValueError(msg)
+ if np.any(np.diff(frequencies) < 0):
+ msg = "frequencies must be in increasing order"
+ raise ValueError(msg)
+
+ npt = npt + 1
+ # Linearly interpolate the target response onto a dense grid, then mirror it
+ # to build the full (symmetric) magnitude spectrum.
+ response = np.interp(np.linspace(0, 1, npt), frequencies, magnitudes)
+ response = np.concatenate([response, response[-2:0:-1]])
+
+ total = response.size
+ half = (total + 1) // 2
+ window_len = 4 * order
+ index = np.arange(window_len)
+
+ # Autocorrelation from the power spectrum, tapered with a Hamming window.
+ correlation = np.real(np.fft.ifft(response * response))
+ correlation = correlation[:window_len] * (
+ 0.54 + 0.46 * np.cos(np.pi * index / (window_len - 1))
+ )
+ cepstral_window = np.concatenate([[0.5], np.ones(half - 1), np.zeros(total - half)])
+
+ # Solve the Yule-Walker normal equations for the denominator coefficients.
+ rmat = toeplitz(correlation[order : window_len - 1], correlation[order:0:-1])
+ rhs = -correlation[order + 1 : window_len]
+ denominator = np.concatenate([[1.0], np.linalg.lstsq(rmat, rhs, rcond=None)[0]])
+ denominator = _polystab(denominator)
+
+ half_correlation = correlation.copy()
+ half_correlation[0] = correlation[0] / 2
+ numerator = _numerator(half_correlation, denominator, order)
+
+ padded_num = np.zeros(total)
+ padded_num[: numerator.size] = numerator
+ padded_den = np.zeros(total)
+ padded_den[: denominator.size] = denominator
+
+ spectrum = 2 * np.real(np.fft.fft(padded_num) / np.fft.fft(padded_den))
+ complex_log = np.log(np.abs(spectrum)) + 1j * np.angle(spectrum)
+ cepstrum = np.fft.ifft(
+ np.exp(np.fft.fft(cepstral_window * np.fft.ifft(complex_log)))
+ )
+ numerator = np.real(_numerator(cepstrum[:window_len], denominator, order))
+ return denominator, numerator
+
+
+class EqualLoudnessFilter:
+ r"""
+ An equal-loudness filter which compensates for the human ear's non-linear
+ response to sound. This filter corrects this by cascading a Yule-Walker
+ filter and a Butterworth filter.
+
+ Designed for use with samplerate of 44.1kHz and above. If you're using a
+ lower samplerate, use with caution.
+
+ Code based on the matlab implementation at https://bit.ly/3eqh2HU
+ (url shortened for ruff)
+
+ Target curve: https://i.imgur.com/3g2VfaM.png
+ Yulewalk response: https://i.imgur.com/J9LnJ4C.png
+ Butterworth and overall response: https://i.imgur.com/3g2VfaM.png
+
+ Images and original matlab implementation by David Robinson, 2001
+
+ https://en.wikipedia.org/wiki/Equal-loudness_contour
+
+ >>> filt = EqualLoudnessFilter()
+ >>> isinstance(filt.yulewalk_filter, IIRFilter)
+ True
+ """
+
+ def __init__(self, samplerate: int = 44100) -> None:
+ self.yulewalk_filter = IIRFilter(10)
+ self.butterworth_filter = make_highpass(150, samplerate)
+
+ # pad the data to nyquist
+ curve_freqs = np.array(data["frequencies"] + [max(20000.0, samplerate / 2)])
+ curve_gains = np.array(data["gains"] + [140])
+
+ # Convert to angular frequency
+ freqs_normalized = curve_freqs / samplerate * 2
+ # Invert the curve and normalize to 0dB
+ gains_normalized = np.power(10, (np.min(curve_gains) - curve_gains) / 20)
+
+ # Compute the coefficients using a least-squares fit to the curve with
+ # the built-in ``yulewalk`` implementation above (no third-party deps).
+ ya, yb = yulewalk(10, freqs_normalized, gains_normalized)
+ self.yulewalk_filter.set_coefficients(ya.tolist(), yb.tolist())
+
+ def process(self, sample: float) -> float:
+ """
+ Process a single sample through both filters
+
+ >>> filt = EqualLoudnessFilter()
+ >>> filt.process(0.0)
+ 0.0
+ """
+ tmp = self.yulewalk_filter.process(sample)
+ return self.butterworth_filter.process(tmp)
diff --git a/audio_filters/equal_loudness_filter.py.broken.txt b/audio_filters/equal_loudness_filter.py.broken.txt
deleted file mode 100644
index 88cba8533cf7..000000000000
--- a/audio_filters/equal_loudness_filter.py.broken.txt
+++ /dev/null
@@ -1,61 +0,0 @@
-from json import loads
-from pathlib import Path
-
-import numpy as np
-from yulewalker import yulewalk
-
-from audio_filters.butterworth_filter import make_highpass
-from audio_filters.iir_filter import IIRFilter
-
-data = loads((Path(__file__).resolve().parent / "loudness_curve.json").read_text())
-
-
-class EqualLoudnessFilter:
- r"""
- An equal-loudness filter which compensates for the human ear's non-linear response
- to sound.
- This filter corrects this by cascading a yulewalk filter and a butterworth filter.
-
- Designed for use with samplerate of 44.1kHz and above. If you're using a lower
- samplerate, use with caution.
-
- Code based on matlab implementation at https://bit.ly/3eqh2HU
- (url shortened for ruff)
-
- Target curve: https://i.imgur.com/3g2VfaM.png
- Yulewalk response: https://i.imgur.com/J9LnJ4C.png
- Butterworth and overall response: https://i.imgur.com/3g2VfaM.png
-
- Images and original matlab implementation by David Robinson, 2001
- """
-
- def __init__(self, samplerate: int = 44100) -> None:
- self.yulewalk_filter = IIRFilter(10)
- self.butterworth_filter = make_highpass(150, samplerate)
-
- # pad the data to nyquist
- curve_freqs = np.array(data["frequencies"] + [max(20000.0, samplerate / 2)])
- curve_gains = np.array(data["gains"] + [140])
-
- # Convert to angular frequency
- freqs_normalized = curve_freqs / samplerate * 2
- # Invert the curve and normalize to 0dB
- gains_normalized = np.power(10, (np.min(curve_gains) - curve_gains) / 20)
-
- # Scipy's `yulewalk` function is a stub, so we're using the
- # `yulewalker` library instead.
- # This function computes the coefficients using a least-squares
- # fit to the specified curve.
- ya, yb = yulewalk(10, freqs_normalized, gains_normalized)
- self.yulewalk_filter.set_coefficients(ya, yb)
-
- def process(self, sample: float) -> float:
- """
- Process a single sample through both filters
-
- >>> filt = EqualLoudnessFilter()
- >>> filt.process(0.0)
- 0.0
- """
- tmp = self.yulewalk_filter.process(sample)
- return self.butterworth_filter.process(tmp)
diff --git a/audio_filters/iir_filter.py b/audio_filters/iir_filter.py
index fa3e6c54b33f..2144e6bd4c7a 100644
--- a/audio_filters/iir_filter.py
+++ b/audio_filters/iir_filter.py
@@ -51,8 +51,37 @@ def set_coefficients(self, a_coeffs: list[float], b_coeffs: list[float]) -> None
... fs=48000)
>>> filt = IIRFilter(2)
>>> filt.set_coefficients(a_coeffs, b_coeffs)
+
+ The leading :math:`a_0` coefficient may be omitted and defaults to 1.0:
+
+ >>> filt = IIRFilter(2)
+ >>> filt.set_coefficients([-1.9, 0.9], [1.0, -2.0, 1.0])
+ >>> filt.a_coeffs
+ [1.0, -1.9, 0.9]
+
+ Passing the wrong number of coefficients raises a ``ValueError``:
+
+ >>> IIRFilter(2).set_coefficients([1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0])
+ Traceback (most recent call last):
+ ...
+ ValueError: Expected a_coeffs to have 3 elements for 2-order filter, got 4
+ >>> IIRFilter(2).set_coefficients([1.0, 2.0, 3.0], [1.0, 2.0])
+ Traceback (most recent call last):
+ ...
+ ValueError: Expected b_coeffs to have 3 elements for 2-order filter, got 2
+
+ A genuinely too-short ``a_coeffs`` is reported with its real length
+ rather than after the optional ``a_0`` has been filled in:
+
+ >>> IIRFilter(2).set_coefficients([1.0], [1.0, 2.0, 3.0])
+ Traceback (most recent call last):
+ ...
+ ValueError: Expected a_coeffs to have 3 elements for 2-order filter, got 1
"""
- if len(a_coeffs) < self.order:
+ if len(a_coeffs) == self.order:
+ # The leading a_0 coefficient is optional; default it to 1.0.
+ # Only a single missing coefficient is filled in this way, so that
+ # genuinely too-short inputs are reported with their real length.
a_coeffs = [1.0, *a_coeffs]
if len(a_coeffs) != self.order + 1:
@@ -65,7 +94,7 @@ def set_coefficients(self, a_coeffs: list[float], b_coeffs: list[float]) -> None
if len(b_coeffs) != self.order + 1:
msg = (
f"Expected b_coeffs to have {self.order + 1} elements "
- f"for {self.order}-order filter, got {len(a_coeffs)}"
+ f"for {self.order}-order filter, got {len(b_coeffs)}"
)
raise ValueError(msg)
diff --git a/audio_filters/show_response.py b/audio_filters/show_response.py
index f9c9537c047c..5830db0a2a38 100644
--- a/audio_filters/show_response.py
+++ b/audio_filters/show_response.py
@@ -1,3 +1,9 @@
+"""
+Plot the magnitude and phase response of an audio filter.
+
+https://en.wikipedia.org/wiki/Frequency_response
+"""
+
from __future__ import annotations
from abc import abstractmethod