diff --git a/plugins/common/dsp/src/lib.rs b/plugins/common/dsp/src/lib.rs index b6c8f2a..9c73ddc 100644 --- a/plugins/common/dsp/src/lib.rs +++ b/plugins/common/dsp/src/lib.rs @@ -1,4 +1,5 @@ //! Shared allocation-free DSP building blocks: filters, pitch detection, notes, WAV files. +pub mod loudness; pub mod meter; pub mod oversample; pub mod svf; diff --git a/plugins/common/dsp/src/loudness.rs b/plugins/common/dsp/src/loudness.rs new file mode 100644 index 0000000..c5fbb67 --- /dev/null +++ b/plugins/common/dsp/src/loudness.rs @@ -0,0 +1,71 @@ +//! ITU-R BS.1770 K-weighting: the pre-filter (a high shelf for the head) and the RLB high-pass. +//! Mean square of the weighted signal, summed over channels, gives the loudness in LUFS as +//! −0.691 + 10·log10(sum). Coefficients come from the standard's analog prototypes, so any +//! sample rate is exact, not only 48 kHz. + +/// Biquad in transposed direct form II, in f64: the RLB pole sits close to 1 at high rates. +#[derive(Clone, Copy, Debug, Default)] +struct Stage { + b: [f64; 3], + a: [f64; 2], + z: [f64; 2], +} +impl Stage { + #[inline] + fn process(&mut self, x: f64) -> f64 { + let y = self.b[0] * x + self.z[0]; + self.z[0] = self.b[1] * x - self.a[0] * y + self.z[1]; + self.z[1] = self.b[2] * x - self.a[1] * y; + y + } +} + +/// One channel of K-weighting. Allocation-free. +#[derive(Clone, Copy, Debug)] +pub struct KWeighting { + shelf: Stage, + high: Stage, +} +impl KWeighting { + pub fn new(rate: f32) -> Self { + let rate = rate as f64; + // Pre-filter: high shelf of about +4 dB above 1.7 kHz. + let (f0, g, q) = (1681.974450955533, 3.999843853973347, 0.7071752369554196); + let k = (std::f64::consts::PI * f0 / rate).tan(); + let vh = 10f64.powf(g / 20.); + let vb = vh.powf(0.4996667741545416); + let a0 = 1. + k / q + k * k; + let shelf = Stage { + b: [ + (vh + vb * k / q + k * k) / a0, + 2. * (k * k - vh) / a0, + (vh - vb * k / q + k * k) / a0, + ], + a: [2. * (k * k - 1.) / a0, (1. - k / q + k * k) / a0], + z: [0.; 2], + }; + // RLB: second-order high-pass at 38 Hz. + let (f0, q) = (38.13547087602444, 0.5003270373238773); + let k = (std::f64::consts::PI * f0 / rate).tan(); + let a0 = 1. + k / q + k * k; + let high = Stage { + b: [1., -2., 1.], + a: [2. * (k * k - 1.) / a0, (1. - k / q + k * k) / a0], + z: [0.; 2], + }; + Self { shelf, high } + } + #[inline] + pub fn process(&mut self, x: f32) -> f32 { + self.high.process(self.shelf.process(x as f64)) as f32 + } + pub fn reset(&mut self) { + self.shelf.z = [0.; 2]; + self.high.z = [0.; 2]; + } +} + +/// Loudness in LUFS of a mean square summed over channels. +pub fn lufs(mean_square: f64) -> f32 { + (-0.691 + 10. * mean_square.max(1e-20).log10()) as f32 +} diff --git a/plugins/common/dsp/tests/dsp.rs b/plugins/common/dsp/tests/dsp.rs index 94fa182..c606985 100644 --- a/plugins/common/dsp/tests/dsp.rs +++ b/plugins/common/dsp/tests/dsp.rs @@ -80,3 +80,29 @@ fn wave_round_trip_keeps_root_and_wavetable_chunks() { assert_eq!(wav::decode(&bytes).unwrap(), w); assert!(wav::decode(&bytes[..40]).is_err()); } + +#[test] +fn k_weighting_reads_a_full_scale_997_hz_sine_at_minus_3_lufs() { + use common_dsp::loudness::{lufs, KWeighting}; + let level = |rate: f32, hz: f64| { + let mut k = KWeighting::new(rate); + let n = rate as usize; + let mut sum = 0f64; + for i in 0..2 * n { + let x = (std::f64::consts::TAU * hz * i as f64 / rate as f64).sin() as f32; + let y = k.process(x) as f64; + if i >= n { + sum += y * y; + } + } + lufs(sum / n as f64) + }; + // BS.1770-4: a 0 dBFS 997 Hz sine in one channel reads −3.01 LKFS, at every rate. + for rate in [44100., 48000., 96000., 192000.] { + let got = level(rate, 997.); + // BS.1770 allows meters 0.1 LU; the prototypes land within 0.03 up to 192 kHz. + assert!((got + 3.01).abs() < 0.05, "{rate}: {got}"); + } + // The RLB high-pass counts a 30 Hz sub much less than the mids. + assert!(level(48000., 30.) < level(48000., 997.) - 5.); +}