Mirror/Python
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Clean up the audio_filters directory (#15087)

Browse Source 
* Clean up the audio_filters directory

- Restore equal_loudness_filter.py: replace the third-party 'yulewalker'
  dependency with a dependency-free 'yulewalk' implementation (numpy + scipy
  only, both already project dependencies) so the filter works again.
- Fix two bugs in iir_filter.set_coefficients: the a0-omitted branch never
  triggered, and the b_coeffs length error reported len(a_coeffs).
- Add doctests covering valid and erroneous inputs (coefficient-length errors,
  yulewalk input validation).
- Add Wikipedia URLs to the docstrings that were missing them.
- Add two new RBJ Audio EQ Cookbook filters: make_notch and make_bandpass_peak.
- Rewrite audio_filters/README.md to document every file and filter.

* audio_filters: clarify make_notch, document biquad notation, tighten set_coefficients guard

- Add a shared notation guide to butterworth_filter so the single-letter
  RBJ-cookbook names (w0/alpha/a0..a2/b0..b2) are self-explanatory across
  every filter instead of diverging in one function.
- Comment make_notch's feed-forward/feed-back sections in plain English.
- set_coefficients only fills the optional a_0 when exactly one coefficient
  is missing, so genuinely too-short inputs raise with their real length
  (new regression doctest).
This commit is contained in:
priya-sundaram-dev
2026-08-27 13:13:49 +02:00
committed by GitHub
parent 25a312ded5
commit 70b9e2fb35
6 changed files with 375 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files
audio_filters/README.md
+52 -3
View File
@@ -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
+88
View File
@@ -7,6 +7,16 @@ Create 2nd-order IIR filters with Butterworth design.
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
+198
View File
@@ -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)
audio_filters/equal_loudness_filter.py.broken.txt
-61
View File
@@ -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)
audio_filters/iir_filter.py
+31 -2
View File
@@ -51,8 +51,37 @@ class IIRFilter:
... 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 @@ class IIRFilter:
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)
audio_filters/show_response.py
+6
View File
@@ -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