feat(common-dsp): add ITU-R BS.1770 K-weighting

KWeighting is one channel of the standard's pre-filter and RLB high-pass,
computed from the analog prototypes so every sample rate is exact. lufs()
turns a mean square summed over channels into LUFS. A full-scale 997 Hz sine
reads -3.01 LUFS within 0.03 LU from 44.1 to 192 kHz.

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;

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
+}

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.);
+}