diff --git a/audio_filters/equal_loudness_filter.py b/audio_filters/equal_loudness_filter.py new file mode 100644 index 000000000000..dbfd51370a9d --- /dev/null +++ b/audio_filters/equal_loudness_filter.py @@ -0,0 +1,160 @@ +from __future__ import annotations + +from audio_filters.iir_filter import IIRFilter + +# Coefficients from the "Original ReplayGain specification" (Equal Loudness Filter) +# - Yulewalk: 10th-order IIR +# - Butterworth: 2nd-order high-pass at 150 Hz +# Source tables / coefficient file: +# - https://wiki.hydrogenaudio.org/index.php?title=Original_ReplayGain_specification +# - https://replaygain.hydrogenaudio.org/equal_loud_coef.txt +# (We embed the coefficients to avoid external dependencies and file I/O.) +_REPLAYGAIN_COEFFS: dict[int, dict[str, list[float]]] = { + 44100: { + "yule_a": [ + 1.0, + -3.47845948550071, + 6.36317777566148, + -8.54751527471874, + 9.47693607801280, + -8.81498681370155, + 6.85401540936998, + -4.39470996079559, + 2.19611684890774, + -0.75104302451432, + 0.13149317958808, + ], + "yule_b": [ + 0.05418656406430, + -0.02911007808948, + -0.00848709379851, + -0.00851165645469, + -0.00834990904936, + 0.02245293253339, + -0.02596338512915, + 0.01624864962975, + -0.00240879051584, + 0.00674613682247, + -0.00187763777362, + ], + "butter_a": [1.0, -1.96977855582618, 0.97022847566350], + "butter_b": [0.98500175787242, -1.97000351574484, 0.98500175787242], + }, + 48000: { + "yule_a": [ + 1.0, + -3.84664617118067, + 7.81501653005538, + -11.34170355132042, + 13.05504219327545, + -12.28759895145294, + 9.48293806319790, + -5.87257861775999, + 2.75465861874613, + -0.86984376593551, + 0.13919314567432, + ], + "yule_b": [ + 0.03857599435200, + -0.02160367184185, + -0.00123395316851, + -0.00009291677959, + -0.01655260341619, + 0.02161526843274, + -0.02074045215285, + 0.00594298065125, + 0.00306428023191, + 0.00012025322027, + 0.00288463683916, + ], + "butter_a": [1.0, -1.97223372919527, 0.97261396931306], + "butter_b": [0.98621192462708, -1.97242384925416, 0.98621192462708], + }, + 32000: { + "yule_a": [ + 1.0, + -2.37898834973084, + 2.84868151156327, + -2.64577170229825, + 2.23697657451713, + -1.67148153367602, + 1.00595954808547, + -0.45953458054983, + 0.16378164858596, + -0.05032077717131, + 0.02347897407020, + ], + "yule_b": [ + 0.15457299681924, + -0.09331049056315, + -0.06247880153653, + 0.02163541888798, + -0.05588393329856, + 0.04781476674921, + 0.00222312597743, + 0.03174092540049, + -0.01390589421898, + 0.00651420667831, + -0.00881362733839, + ], + "butter_a": [1.0, -1.95835380975398, 0.95920349965459], + "butter_b": [0.97938932735214, -1.95877865470428, 0.97938932735214], + }, +} + + +class EqualLoudnessFilter: + r""" + Equal-loudness compensation filter (ReplayGain-style). + + This is a cascade of: + - 10th-order "yulewalk" IIR filter + - 2nd-order Butterworth high-pass filter at 150 Hz + + Coefficients are embedded for a few common sample rates, matching the + Original ReplayGain specification. :contentReference[oaicite:1]{index=1} + + >>> filt = EqualLoudnessFilter(44100) + >>> filt.process(0.0) + 0.0 + + >>> EqualLoudnessFilter(12345) + Traceback (most recent call last): + ... + ValueError: Unsupported samplerate 12345. Supported samplerates: 32000, 44100, 48000 + """ + + def __init__(self, samplerate: int = 44100) -> None: + if samplerate not in _REPLAYGAIN_COEFFS: + supported = ", ".join(str(sr) for sr in sorted(_REPLAYGAIN_COEFFS)) + msg = ( + f"Unsupported samplerate {samplerate}. " + f"Supported samplerates: {supported}" + ) + raise ValueError(msg) + + coeffs = _REPLAYGAIN_COEFFS[samplerate] + + self.yulewalk_filter = IIRFilter(10) + self.yulewalk_filter.set_coefficients(coeffs["yule_a"], coeffs["yule_b"]) + + self.butterworth_filter = IIRFilter(2) + self.butterworth_filter.set_coefficients(coeffs["butter_a"], coeffs["butter_b"]) + + def reset(self) -> None: + """Reset the internal filter histories to zero.""" + self.yulewalk_filter.input_history = [0.0] * self.yulewalk_filter.order + self.yulewalk_filter.output_history = [0.0] * self.yulewalk_filter.order + self.butterworth_filter.input_history = [0.0] * self.butterworth_filter.order + self.butterworth_filter.output_history = [0.0] * self.butterworth_filter.order + + 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..cbe92d204d6d 100644 --- a/audio_filters/iir_filter.py +++ b/audio_filters/iir_filter.py @@ -19,8 +19,8 @@ class IIRFilter: we can rewrite this to .. math:: 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) + \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) """ def __init__(self, order: int) -> None: @@ -44,15 +44,29 @@ def set_coefficients(self, a_coeffs: list[float], b_coeffs: list[float]) -> None This method works well with scipy's filter design functions - >>> # Make a 2nd-order 1000Hz butterworth lowpass filter - >>> import scipy.signal - >>> b_coeffs, a_coeffs = scipy.signal.butter(2, 1000, - ... btype='lowpass', - ... fs=48000) + >>> # Make a 2nd-order 1000Hz butterworth lowpass filter (SciPy optional) + >>> import scipy.signal # doctest: +SKIP + >>> b_coeffs, a_coeffs = scipy.signal.butter( # doctest: +SKIP + ... 2, 1000, btype="lowpass", fs=48000 + ... ) + >>> filt = IIRFilter(2) # doctest: +SKIP + >>> filt.set_coefficients(a_coeffs, b_coeffs) # doctest: +SKIP + + >>> # a0 can be omitted (defaults to 1.0) + >>> filt = IIRFilter(2) + >>> filt.set_coefficients([0.5, 0.25], [0.1, 0.2, 0.3]) + >>> filt.a_coeffs + [1.0, 0.5, 0.25] + + >>> # b_coeffs length check should report len(b_coeffs) >>> filt = IIRFilter(2) - >>> filt.set_coefficients(a_coeffs, b_coeffs) + >>> filt.set_coefficients([1.0, 0.5, 0.25], [0.1, 0.2]) + Traceback (most recent call last): + ... + ValueError: Expected b_coeffs to have 3 elements for 2-order filter, got 2 """ - if len(a_coeffs) < self.order: + # Allow omitting a0 (use 1.0 as default) when only a1..ak are provided + if len(a_coeffs) == self.order: a_coeffs = [1.0, *a_coeffs] if len(a_coeffs) != self.order + 1: @@ -65,7 +79,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)