Clean up the audio_filters directory

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:
 
 * <https://www.masteringbox.com/filter-types/>
 * <http://ethanwiner.com/filters.html>
 * <https://en.wikipedia.org/wiki/Audio_filter>
 * <https://en.wikipedia.org/wiki/Electronic_filter>
+* <https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html>
+
+## 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.

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

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)