diff --git a/Cargo.lock b/Cargo.lock index da34901e..09095698 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1308,6 +1308,10 @@ dependencies = [ "thiserror 2.0.21", ] +[[package]] +name = "filmcraft-audio-dsp" +version = "0.1.0" + [[package]] name = "filmcraft-automation" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index a809ee9e..0093b367 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -23,6 +23,7 @@ filmcraft-hevc = { path = "crates/hevc" } filmcraft-prores = { path = "crates/prores" } filmcraft-aac = { path = "crates/aac" } filmcraft-h264enc = { path = "crates/h264enc" } +filmcraft-audio-dsp = { path = "crates/audio-dsp" } filmcraft-frame = { path = "crates/frame" } filmcraft-media = { path = "crates/media" } filmcraft-codecs = { path = "crates/codecs" } diff --git a/crates/audio-dsp/Cargo.toml b/crates/audio-dsp/Cargo.toml new file mode 100644 index 00000000..bc1a0483 --- /dev/null +++ b/crates/audio-dsp/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "filmcraft-audio-dsp" +description = "Real-time audio DSP for FilmCraft: BS.1770 / EBU R128 loudness metering and audio effects" +version.workspace = true +edition.workspace = true +license.workspace = true +rust-version.workspace = true + +[lints] +workspace = true + +# Intentionally dependency-free (layer L1): planar `&mut [f32]` slices in, planar slices out. +[dependencies] diff --git a/crates/audio-dsp/README.md b/crates/audio-dsp/README.md new file mode 100644 index 00000000..99c85a0f --- /dev/null +++ b/crates/audio-dsp/README.md @@ -0,0 +1,57 @@ +# filmcraft-audio-dsp + +Real-time audio DSP for FilmCraft. Layer **L1**, **no dependencies**, compiles for +`wasm32-unknown-unknown`, no `unsafe`. Clean-room: written from the public standards +(ITU-R BS.1770-4, EBU R128 / Tech 3341 / Tech 3342) and textbook DSP (RBJ cookbook biquads, +Householder FDN, phase vocoder). + +## Conventions + +- Planar `f32` audio: `&mut [&mut [f32]]`, one slice per channel, equal lengths. +- Effects are built (allocating) with `(sample_rate, channels)`; `process`, `set_param` and + `reset` never allocate. +- Parameters change between blocks and are smoothed per sample; output does not depend on how + the stream is cut into blocks (tested for every effect with 0/1/random-size blocks). +- Recursive state is flushed below ~1e-30, so no denormals without CPU FTZ flags. + +## Loudness (`loudness`) + +`LoudnessMeter::new(sample_rate, channels)` → `process(&[&[f32]])` / `process_interleaved`, then +`momentary()`, `short_term()`, `integrated()`, `loudness_range()`, `sample_peak()`, +`true_peak()` / `true_peak_dbtp()`, `summary()`. + +- K-weighting derived for any sample rate by bilinear transform (matches the 48 kHz BS.1770 table). +- Gated integrated loudness (−70 LUFS / −10 LU) and LRA (−70 / −20 LU, 10th–95th percentile) use a + 0.01 LU histogram with exact energy sums: fixed memory for programmes of any length. +- True peak: 4× (2× at 88.2/96 kHz) polyphase Kaiser-sinc interpolation. +- `normalize_gain_db(measured, target)` and `normalize_gain_db_peak_limited(...)` for Normalize / + Essential Sound loudness auto-match. + +## Effects (`effects`) + +All implement `AudioEffect` (`id`, `params`, `set_param`, `param`, `reset`, `process`, `latency`). +`effects()` lists `EffectInfo { id, name, category, params, create }`; each `ParamSpec` has +id, name, range, default, `Unit`, log-scale hint and choice labels so UIs can be generated. + +| id | effect | +|---|---| +| `parametric_eq` | 8 bands, each peaking / shelf / LP / HP / notch / band-pass, RBJ biquads in TDF-II | +| `simple_eq` | low shelf, mid peak, high shelf | +| `compressor` | peak/RMS detector, soft knee, attack/release, make-up; stereo-linked | +| `gate` | downward expander / gate with range and hold | +| `limiter` | 5 ms look-ahead, true-peak-aware, ceiling guaranteed (latency reported) | +| `delay` | feedback delay, ms or tempo-synced divisions, mix | +| `reverb` | 8×8 Householder FDN, pre-delay, RT60 decay, HF damping, size, mix | +| `amplify`, `channel_volume`, `balance` | gain, per-channel gain, balance / constant-power pan | +| `invert`, `swap_channels`, `fill_left`, `fill_right` | polarity and routing | +| `dehum` | 50/60 Hz notch comb with up to 10 harmonics | +| `denoise` | STFT spectral gating, minimum-statistics noise floor (latency = FFT size) | +| `pitch_shifter` | ±12 semitones, phase vocoder with peak shifting + phase locking (latency = FFT size) | + +## Known limits + +- Effects assume ≤ the channel count they were built with; extra channels pass through + (channel-generic gain effects process all of them). +- DeNoise adapts continuously; a perfectly stationary tone is (by design) treated as noise after + ~1.5 s. There is no explicit "learn noise print" mode yet. +- Filter type / choice changes switch immediately (continuous parameters are smoothed). diff --git a/crates/audio-dsp/src/biquad.rs b/crates/audio-dsp/src/biquad.rs new file mode 100644 index 00000000..b213d7b9 --- /dev/null +++ b/crates/audio-dsp/src/biquad.rs @@ -0,0 +1,219 @@ +//! Second-order IIR sections: RBJ "Audio EQ Cookbook" designs, transposed direct form II. + +use crate::flush; +use std::f64::consts::PI; + +/// Filter shape. +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub enum FilterType { + Peaking, + LowShelf, + HighShelf, + LowPass, + HighPass, + Notch, + BandPass, + AllPass, +} + +impl FilterType { + /// Order used by the parametric EQ's `type` choice parameter. + pub const ALL: [FilterType; 7] = + [FilterType::Peaking, FilterType::LowShelf, FilterType::HighShelf, FilterType::LowPass, FilterType::HighPass, FilterType::Notch, FilterType::BandPass]; + pub fn from_index(i: usize) -> FilterType { + Self::ALL[i.min(Self::ALL.len() - 1)] + } +} + +/// Normalised biquad coefficients (`a0 == 1`). +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Coeffs { + pub b0: f64, + pub b1: f64, + pub b2: f64, + pub a1: f64, + pub a2: f64, +} + +impl Coeffs { + /// Pass-through. + pub const IDENTITY: Coeffs = Coeffs { b0: 1.0, b1: 0.0, b2: 0.0, a1: 0.0, a2: 0.0 }; + + fn norm(b0: f64, b1: f64, b2: f64, a0: f64, a1: f64, a2: f64) -> Coeffs { + Coeffs { b0: b0 / a0, b1: b1 / a0, b2: b2 / a0, a1: a1 / a0, a2: a2 / a0 } + } + + /// RBJ cookbook design. `freq` in Hz (clamped below Nyquist), `q` > 0, `gain_db` used by + /// peaking and shelving types (shelves use `q` as the shelf "Q", S=1 ≙ Q=1/√2). + pub fn design(kind: FilterType, sample_rate: f64, freq: f64, q: f64, gain_db: f64) -> Coeffs { + let f = freq.clamp(1.0, sample_rate * 0.49); + let q = q.max(1e-3); + let w0 = 2.0 * PI * f / sample_rate; + let (sn, cs) = w0.sin_cos(); + let alpha = sn / (2.0 * q); + let a = 10f64.powf(gain_db / 40.0); + match kind { + FilterType::LowPass => Self::norm((1.0 - cs) / 2.0, 1.0 - cs, (1.0 - cs) / 2.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha), + FilterType::HighPass => Self::norm((1.0 + cs) / 2.0, -(1.0 + cs), (1.0 + cs) / 2.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha), + FilterType::BandPass => Self::norm(alpha, 0.0, -alpha, 1.0 + alpha, -2.0 * cs, 1.0 - alpha), + FilterType::Notch => Self::norm(1.0, -2.0 * cs, 1.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha), + FilterType::AllPass => Self::norm(1.0 - alpha, -2.0 * cs, 1.0 + alpha, 1.0 + alpha, -2.0 * cs, 1.0 - alpha), + FilterType::Peaking => Self::norm(1.0 + alpha * a, -2.0 * cs, 1.0 - alpha * a, 1.0 + alpha / a, -2.0 * cs, 1.0 - alpha / a), + FilterType::LowShelf => { + let k = 2.0 * a.sqrt() * alpha; + Self::norm( + a * ((a + 1.0) - (a - 1.0) * cs + k), + 2.0 * a * ((a - 1.0) - (a + 1.0) * cs), + a * ((a + 1.0) - (a - 1.0) * cs - k), + (a + 1.0) + (a - 1.0) * cs + k, + -2.0 * ((a - 1.0) + (a + 1.0) * cs), + (a + 1.0) + (a - 1.0) * cs - k, + ) + } + FilterType::HighShelf => { + let k = 2.0 * a.sqrt() * alpha; + Self::norm( + a * ((a + 1.0) + (a - 1.0) * cs + k), + -2.0 * a * ((a - 1.0) + (a + 1.0) * cs), + a * ((a + 1.0) + (a - 1.0) * cs - k), + (a + 1.0) - (a - 1.0) * cs + k, + 2.0 * ((a - 1.0) - (a + 1.0) * cs), + (a + 1.0) - (a - 1.0) * cs - k, + ) + } + } + } + + /// |H(e^{jω})| at `freq` Hz. + pub fn magnitude(&self, freq: f64, sample_rate: f64) -> f64 { + let w = 2.0 * PI * freq / sample_rate; + let (s1, c1) = w.sin_cos(); + let (s2, c2) = (2.0 * w).sin_cos(); + // H = (b0 + b1 z^-1 + b2 z^-2) / (1 + a1 z^-1 + a2 z^-2), z^-1 = e^{-jω}. + let nr = self.b0 + self.b1 * c1 + self.b2 * c2; + let ni = -(self.b1 * s1 + self.b2 * s2); + let dr = 1.0 + self.a1 * c1 + self.a2 * c2; + let di = -(self.a1 * s1 + self.a2 * s2); + ((nr * nr + ni * ni) / (dr * dr + di * di)).sqrt() + } + + /// Magnitude in dB. + pub fn magnitude_db(&self, freq: f64, sample_rate: f64) -> f64 { + 20.0 * self.magnitude(freq, sample_rate).log10() + } +} + +/// One biquad section's state (TDF-II), f64 internally. +#[derive(Clone, Copy, Debug, Default)] +pub struct Biquad { + s1: f64, + s2: f64, +} + +impl Biquad { + pub const fn new() -> Self { + Biquad { s1: 0.0, s2: 0.0 } + } + #[inline(always)] + pub fn tick(&mut self, c: &Coeffs, x: f64) -> f64 { + let y = c.b0 * x + self.s1; + self.s1 = flush(c.b1 * x - c.a1 * y + self.s2); + self.s2 = flush(c.b2 * x - c.a2 * y); + y + } + pub fn reset(&mut self) { + *self = Biquad::new(); + } + /// Filter a buffer in place. + pub fn process(&mut self, c: &Coeffs, buf: &mut [f32]) { + for v in buf { + *v = self.tick(c, *v as f64) as f32; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + const SR: f64 = 48000.0; + + fn measured_db(c: &Coeffs, freq: f64) -> f64 { + let n = SR as usize; + let mut x = sine(freq, 0.5, SR, n, 0.0); + let mut bq = Biquad::new(); + bq.process(c, &mut x); + // Skip the transient; measure over an integer number of periods-ish long tail. + let tail = &x[n / 2..]; + db(tone_amplitude(tail, freq, SR) / 0.5) + } + + #[test] + fn peaking_hits_gain_at_centre() { + let c = Coeffs::design(FilterType::Peaking, SR, 1000.0, 1.0, 6.0); + assert!((c.magnitude_db(1000.0, SR) - 6.0).abs() < 1e-9); + assert!(c.magnitude_db(20.0, SR).abs() < 0.05); + assert!((measured_db(&c, 1000.0) - 6.0).abs() < 0.05); + } + + #[test] + fn lowpass_highpass_minus_3db_at_cutoff() { + let q = std::f64::consts::FRAC_1_SQRT_2; + let lp = Coeffs::design(FilterType::LowPass, SR, 2000.0, q, 0.0); + let hp = Coeffs::design(FilterType::HighPass, SR, 2000.0, q, 0.0); + assert!((lp.magnitude_db(2000.0, SR) + 3.0103).abs() < 0.01); + assert!((hp.magnitude_db(2000.0, SR) + 3.0103).abs() < 0.01); + assert!(lp.magnitude_db(100.0, SR).abs() < 0.01); + assert!(hp.magnitude_db(15000.0, SR).abs() < 0.1); + // 12 dB/oct slope well past cutoff. + assert!(lp.magnitude_db(8000.0, SR) < -22.0); + assert!((measured_db(&lp, 2000.0) + 3.0103).abs() < 0.05); + assert!((measured_db(&hp, 500.0) - hp.magnitude_db(500.0, SR)).abs() < 0.05); + } + + #[test] + fn notch_is_deep_at_centre() { + let c = Coeffs::design(FilterType::Notch, SR, 60.0, 10.0, 0.0); + assert!(c.magnitude(60.0, SR) < 1e-6); + assert!(c.magnitude_db(1000.0, SR).abs() < 0.01); + assert!(measured_db(&c, 60.0) < -40.0); + } + + #[test] + fn shelves_reach_their_gain() { + let q = std::f64::consts::FRAC_1_SQRT_2; + let ls = Coeffs::design(FilterType::LowShelf, SR, 200.0, q, 9.0); + let hs = Coeffs::design(FilterType::HighShelf, SR, 5000.0, q, -9.0); + assert!((ls.magnitude_db(10.0, SR) - 9.0).abs() < 0.05); + assert!((ls.magnitude_db(200.0, SR) - 4.5).abs() < 0.01); + assert!(ls.magnitude_db(10000.0, SR).abs() < 0.05); + assert!((hs.magnitude_db(23000.0, SR) + 9.0).abs() < 0.1); + assert!((hs.magnitude_db(5000.0, SR) + 4.5).abs() < 0.01); + assert!((measured_db(&ls, 40.0) - ls.magnitude_db(40.0, SR)).abs() < 0.05); + assert!((measured_db(&hs, 12000.0) - hs.magnitude_db(12000.0, SR)).abs() < 0.05); + } + + #[test] + fn bandpass_and_allpass() { + let bp = Coeffs::design(FilterType::BandPass, SR, 1000.0, 2.0, 0.0); + assert!(bp.magnitude_db(1000.0, SR).abs() < 1e-9); + let ap = Coeffs::design(FilterType::AllPass, SR, 1000.0, 2.0, 0.0); + for f in [50.0, 1000.0, 9000.0] { + assert!(ap.magnitude_db(f, SR).abs() < 1e-9); + } + } + + #[test] + fn state_flushes_to_zero_no_denormals() { + let c = Coeffs::design(FilterType::LowPass, SR, 100.0, 0.7, 0.0); + let mut bq = Biquad::new(); + bq.tick(&c, 1.0); + let mut last = 1.0; + for _ in 0..2_000_000 { + last = bq.tick(&c, 0.0); + assert!(!(last as f32).is_subnormal()); + } + assert_eq!(last, 0.0); + } +} diff --git a/crates/audio-dsp/src/effects/channel.rs b/crates/audio-dsp/src/effects/channel.rs new file mode 100644 index 00000000..9fe275d8 --- /dev/null +++ b/crates/audio-dsp/src/effects/channel.rs @@ -0,0 +1,340 @@ +//! Gain and channel utilities: amplify, channel volume, balance/pan, invert, swap, fill. + +use crate::{AudioEffect, ParamSpec, ParamValues, Smoothed, Unit, block_len, db_to_gain, param_plumbing}; +use std::f32::consts::FRAC_PI_2; + +const RAMP_MS: f32 = 20.0; + +/// Amplify: one smoothed gain for all channels. +pub struct Amplify { + pv: ParamValues, + gain: Smoothed, +} + +impl Amplify { + pub const PARAMS: &'static [ParamSpec] = &[ParamSpec::new("gain", "Gain", -96.0, 24.0, 0.0, Unit::Decibels)]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let mut s = Amplify { pv: ParamValues::new(Self::PARAMS), gain: Smoothed::with_ms(1.0, sample_rate, RAMP_MS) }; + s.apply_params(true); + s + } + fn apply_params(&mut self, snap: bool) { + let db = self.pv.v("gain"); + self.gain.set(if db <= -96.0 { 0.0 } else { db_to_gain(db) }); + if snap { + self.gain.snap(); + } + } +} + +impl AudioEffect for Amplify { + param_plumbing!("amplify"); + fn reset(&mut self) { + self.gain.snap(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + for i in 0..n { + let g = self.gain.tick(); + for ch in channels.iter_mut() { + ch[i] *= g; + } + } + } +} + +/// Per-channel volume (L R C LFE Ls Rs order). +pub struct ChannelVolume { + pv: ParamValues, + gains: [Smoothed; 6], +} + +const CHANNEL_IDS: [&str; 6] = ["left", "right", "center", "lfe", "left_surround", "right_surround"]; + +impl ChannelVolume { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("left", "Left", -96.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("right", "Right", -96.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("center", "Center", -96.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("lfe", "LFE", -96.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("left_surround", "Left Surround", -96.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("right_surround", "Right Surround", -96.0, 24.0, 0.0, Unit::Decibels), + ]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let mut s = ChannelVolume { pv: ParamValues::new(Self::PARAMS), gains: [Smoothed::with_ms(1.0, sample_rate, RAMP_MS); 6] }; + s.apply_params(true); + s + } + fn apply_params(&mut self, snap: bool) { + for (g, id) in self.gains.iter_mut().zip(CHANNEL_IDS) { + let db = self.pv.v(id); + g.set(if db <= -96.0 { 0.0 } else { db_to_gain(db) }); + if snap { + g.snap(); + } + } + } +} + +impl AudioEffect for ChannelVolume { + param_plumbing!("channel_volume"); + fn reset(&mut self) { + self.gains.iter_mut().for_each(Smoothed::snap); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + for i in 0..n { + for (c, g) in self.gains.iter_mut().enumerate() { + let g = g.tick(); + if let Some(ch) = channels.get_mut(c) { + ch[i] *= g; + } + } + } + } +} + +/// Constant-power gains `(left, right)` for a pan position in −1 … +1 (centre = −3 dB each). +pub fn constant_power_pan(p: f32) -> (f32, f32) { + let theta = (p.clamp(-1.0, 1.0) + 1.0) * 0.25 * std::f32::consts::PI; + (theta.cos(), theta.sin()) +} + +/// Balance gains for a stereo source: the side opposite the direction is attenuated with a +/// cosine (constant-power) law, the other side stays at unity. Centre = unity both sides. +pub fn balance_gains(p: f32) -> (f32, f32) { + let p = p.clamp(-1.0, 1.0); + ((p.max(0.0) * FRAC_PI_2).cos(), ((-p).max(0.0) * FRAC_PI_2).cos()) +} + +/// Stereo balance / pan. "Balance" attenuates one side of a stereo source; "Pan" sums to mono +/// and positions it with the constant-power law (normalised so the centre is unity for a +/// correlated source). +pub struct Balance { + pv: ParamValues, + gl: Smoothed, + gr: Smoothed, + pan_mode: bool, + mode_mix: Smoothed, +} + +impl Balance { + pub const PARAMS: &'static [ParamSpec] = + &[ParamSpec::new("balance", "Balance", -100.0, 100.0, 0.0, Unit::Pan), ParamSpec::choice("mode", "Mode", &["Balance", "Pan (constant power)"], 0)]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let mut s = Balance { + pv: ParamValues::new(Self::PARAMS), + gl: Smoothed::with_ms(1.0, sample_rate, RAMP_MS), + gr: Smoothed::with_ms(1.0, sample_rate, RAMP_MS), + pan_mode: false, + mode_mix: Smoothed::with_ms(0.0, sample_rate, RAMP_MS), + }; + s.apply_params(true); + s + } + fn apply_params(&mut self, snap: bool) { + let p = self.pv.v("balance") / 100.0; + self.pan_mode = self.pv.idx("mode") == 1; + let (l, r) = if self.pan_mode { + let (l, r) = constant_power_pan(p); + (l * std::f32::consts::SQRT_2, r * std::f32::consts::SQRT_2) + } else { + balance_gains(p) + }; + self.gl.set(l); + self.gr.set(r); + self.mode_mix.set(if self.pan_mode { 1.0 } else { 0.0 }); + if snap { + self.gl.snap(); + self.gr.snap(); + self.mode_mix.snap(); + } + } +} + +impl AudioEffect for Balance { + param_plumbing!("balance"); + fn reset(&mut self) { + self.gl.snap(); + self.gr.snap(); + self.mode_mix.snap(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + if channels.len() < 2 { + return; + } + let (a, b) = channels.split_at_mut(1); + let (l, r) = (&mut a[0], &mut b[0]); + for i in 0..n { + let gl = self.gl.tick(); + let gr = self.gr.tick(); + let mm = self.mode_mix.tick(); + let mono = 0.5 * (l[i] + r[i]); + let il = l[i] + mm * (mono - l[i]); + let ir = r[i] + mm * (mono - r[i]); + l[i] = il * gl; + r[i] = ir * gr; + } + } +} + +/// Polarity inversion. +pub struct Invert { + pv: ParamValues, +} + +impl Invert { + pub const PARAMS: &'static [ParamSpec] = &[ParamSpec::choice("channels", "Channels", &["All", "Left", "Right"], 0)]; + pub fn new(_sample_rate: f32, _channels: usize) -> Self { + Invert { pv: ParamValues::new(Self::PARAMS) } + } + fn apply_params(&mut self, _snap: bool) {} +} + +impl AudioEffect for Invert { + param_plumbing!("invert"); + fn reset(&mut self) {} + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + let which = self.pv.idx("channels"); + for (c, ch) in channels.iter_mut().enumerate() { + if which == 0 || which == c + 1 { + ch[..n].iter_mut().for_each(|v| *v = -*v); + } + } + } +} + +/// Parameterless channel routing utilities. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Routing { + /// Swap left and right. + Swap, + /// Fill left with right. + FillLeft, + /// Fill right with left. + FillRight, +} + +/// Swap / Fill Left with Right / Fill Right with Left. +pub struct Route { + pv: ParamValues, + kind: Routing, +} + +impl Route { + pub const PARAMS: &'static [ParamSpec] = &[]; + pub fn new(kind: Routing) -> Self { + Route { pv: ParamValues::new(Self::PARAMS), kind } + } + fn apply_params(&mut self, _snap: bool) {} +} + +impl AudioEffect for Route { + fn id(&self) -> &'static str { + match self.kind { + Routing::Swap => "swap_channels", + Routing::FillLeft => "fill_left", + Routing::FillRight => "fill_right", + } + } + fn params(&self) -> &'static [ParamSpec] { + Self::PARAMS + } + fn set_param(&mut self, id: &str, value: f32) -> bool { + let ok = self.pv.set(id, value); + self.apply_params(false); + ok + } + fn param(&self, id: &str) -> Option { + self.pv.get(id) + } + fn reset(&mut self) {} + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + if channels.len() < 2 { + return; + } + let (a, b) = channels.split_at_mut(1); + let (l, r) = (&mut a[0][..n], &mut b[0][..n]); + match self.kind { + Routing::Swap => l.swap_with_slice(r), + Routing::FillLeft => l.copy_from_slice(r), + Routing::FillRight => r.copy_from_slice(l), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn amplify_reaches_exact_gain() { + let mut a = Amplify::new(48000.0, 2); + a.set_param("gain", -6.0); + let mut l = vec![1.0f32; 2000]; + let mut r = vec![0.5f32; 2000]; + a.process(&mut [&mut l, &mut r]); + assert!((l[1999] - db_to_gain(-6.0)).abs() < 1e-6); + assert!((r[1999] - 0.5 * db_to_gain(-6.0)).abs() < 1e-6); + // Smoothed, not stepped. + assert!(l[0] > 0.99); + a.set_param("gain", -96.0); + a.process(&mut [&mut l, &mut r]); + assert_eq!(l[1999], 0.0); + } + + #[test] + fn channel_volume_per_channel() { + let mut c = ChannelVolume::new(48000.0, 2); + c.set_param("right", -20.0); + c.reset(); + let mut l = vec![1.0f32; 10]; + let mut r = vec![1.0f32; 10]; + c.process(&mut [&mut l, &mut r]); + assert_eq!(l[5], 1.0); + assert!((r[5] - 0.1).abs() < 1e-6); + } + + #[test] + fn pan_is_constant_power() { + for i in -10..=10 { + let (l, r) = constant_power_pan(i as f32 / 10.0); + assert!((l * l + r * r - 1.0).abs() < 1e-6); + } + assert_eq!(balance_gains(0.0), (1.0, 1.0)); + let (l, r) = balance_gains(1.0); + assert!(l.abs() < 1e-6 && r == 1.0); + let mut b = Balance::new(48000.0, 2); + b.set_param("mode", 1.0); + b.set_param("balance", 100.0); + b.reset(); + let mut l = vec![0.5f32; 4]; + let mut r = vec![0.5f32; 4]; + b.process(&mut [&mut l, &mut r]); + assert!(l[3].abs() < 1e-6); + assert!((r[3] - 0.5 * std::f32::consts::SQRT_2).abs() < 1e-6); + } + + #[test] + fn invert_swap_fill() { + let mut l = vec![1.0f32, 2.0]; + let mut r = vec![3.0f32, 4.0]; + Route::new(Routing::Swap).process(&mut [&mut l, &mut r]); + assert_eq!((l.clone(), r.clone()), (vec![3.0, 4.0], vec![1.0, 2.0])); + Route::new(Routing::FillLeft).process(&mut [&mut l, &mut r]); + assert_eq!(l, vec![1.0, 2.0]); + l[0] = 9.0; + Route::new(Routing::FillRight).process(&mut [&mut l, &mut r]); + assert_eq!(r, vec![9.0, 2.0]); + let mut inv = Invert::new(48000.0, 2); + inv.set_param("channels", 2.0); + inv.process(&mut [&mut l, &mut r]); + assert_eq!((l, r), (vec![9.0, 2.0], vec![-9.0, -2.0])); + } +} diff --git a/crates/audio-dsp/src/effects/dynamics.rs b/crates/audio-dsp/src/effects/dynamics.rs new file mode 100644 index 00000000..14f89fff --- /dev/null +++ b/crates/audio-dsp/src/effects/dynamics.rs @@ -0,0 +1,552 @@ +//! Dynamics: compressor, expander/gate, look-ahead true-peak limiter. +//! +//! Detection is stereo-linked (the loudest channel drives one common gain) so the image does +//! not shift. Gain computers work in dB; attack/release are one-pole smoothers on the gain. + +use crate::oversample::{Interpolator, PolyphaseBank, true_peak_factor}; +use crate::{AudioEffect, ParamSpec, ParamValues, Smoothed, Unit, block_len, db_to_gain, param_plumbing}; +use std::collections::VecDeque; + +/// One-pole coefficient for a time constant in ms. +fn coef(ms: f32, sr: f32) -> f32 { + if ms <= 0.0 { 0.0 } else { (-1.0 / (ms * 0.001 * sr)).exp() } +} + +#[inline] +fn lin_to_db(x: f32) -> f32 { + 20.0 * (x.max(1e-10)).log10() +} + +/// Static compressor curve with a quadratic soft knee (output level in dB for input level `x`). +pub fn compressor_curve(x: f32, threshold: f32, ratio: f32, knee: f32) -> f32 { + let d = x - threshold; + if knee > 0.0 && 2.0 * d.abs() <= knee { + x + (1.0 / ratio - 1.0) * (d + knee / 2.0).powi(2) / (2.0 * knee) + } else if d > 0.0 { + threshold + d / ratio + } else { + x + } +} + +/// Static downward-expander curve (output level in dB), attenuation limited to `range` dB. +pub fn expander_curve(x: f32, threshold: f32, ratio: f32, range: f32) -> f32 { + if x >= threshold { x } else { x + ((x - threshold) * (ratio - 1.0)).max(-range) } +} + +// --------------------------------------------------------------------------------------------- + +/// Feed-forward compressor (peak or RMS detector, soft knee, attack/release, make-up gain). +pub struct Compressor { + pv: ParamValues, + sr: f32, + threshold: f32, + ratio: f32, + knee: f32, + rms: bool, + att: f32, + rel: f32, + rms_coef: f32, + makeup: Smoothed, + /// Detector mean square (RMS mode). + ms: f32, + /// Smoothed gain reduction (dB, ≤ 0). + env: f32, + gr_meter: f32, +} + +impl Compressor { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("threshold", "Threshold", -60.0, 0.0, -20.0, Unit::Decibels), + ParamSpec::new("ratio", "Ratio", 1.0, 30.0, 4.0, Unit::Ratio), + ParamSpec::new("knee", "Knee", 0.0, 24.0, 6.0, Unit::Decibels), + ParamSpec::log("attack", "Attack", 0.01, 300.0, 10.0, Unit::Milliseconds), + ParamSpec::log("release", "Release", 1.0, 3000.0, 100.0, Unit::Milliseconds), + ParamSpec::new("makeup", "Make-up Gain", 0.0, 40.0, 0.0, Unit::Decibels), + ParamSpec::choice("detector", "Detector", &["Peak", "RMS"], 0), + ]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let mut s = Compressor { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + threshold: 0.0, + ratio: 1.0, + knee: 0.0, + rms: false, + att: 0.0, + rel: 0.0, + rms_coef: coef(10.0, sample_rate), + makeup: Smoothed::with_ms(1.0, sample_rate, 20.0), + ms: 0.0, + env: 0.0, + gr_meter: 0.0, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + self.threshold = self.pv.v("threshold"); + self.ratio = self.pv.v("ratio"); + self.knee = self.pv.v("knee"); + self.att = coef(self.pv.v("attack"), self.sr); + self.rel = coef(self.pv.v("release"), self.sr); + self.rms = self.pv.idx("detector") == 1; + self.makeup.set(db_to_gain(self.pv.v("makeup"))); + if snap { + self.makeup.snap(); + } + } + + /// Static curve at the current settings. + pub fn curve(&self, x_db: f32) -> f32 { + compressor_curve(x_db, self.threshold, self.ratio, self.knee) + } + + /// Current gain reduction in dB (≤ 0), for metering. + pub fn gain_reduction_db(&self) -> f32 { + self.gr_meter + } +} + +impl AudioEffect for Compressor { + param_plumbing!("compressor"); + fn reset(&mut self) { + self.makeup.snap(); + self.ms = 0.0; + self.env = 0.0; + self.gr_meter = 0.0; + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + for i in 0..n { + let mut peak = 0.0f32; + let mut sq = 0.0f32; + for ch in channels.iter() { + let x = ch[i]; + peak = peak.max(x.abs()); + sq = sq.max(x * x); + } + let level_db = if self.rms { + self.ms = self.rms_coef * self.ms + (1.0 - self.rms_coef) * sq; + if self.ms < 1e-30 { + self.ms = 0.0; + } + 10.0 * self.ms.max(1e-20).log10() + } else { + lin_to_db(peak) + }; + let target = self.curve(level_db) - level_db; + let c = if target < self.env { self.att } else { self.rel }; + self.env = c * self.env + (1.0 - c) * target; + if self.env.abs() < 1e-9 { + self.env = 0.0; + } + let g = db_to_gain(self.env) * self.makeup.tick(); + for ch in channels.iter_mut() { + ch[i] *= g; + } + } + self.gr_meter = self.env; + } +} + +// --------------------------------------------------------------------------------------------- + +/// Downward expander / noise gate with hold (ratio ≥ 10 behaves as a gate). +pub struct Gate { + pv: ParamValues, + sr: f32, + threshold: f32, + ratio: f32, + range: f32, + att: f32, + rel: f32, + hold: u32, + hold_left: u32, + env: f32, +} + +impl Gate { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("threshold", "Threshold", -90.0, 0.0, -40.0, Unit::Decibels), + ParamSpec::log("ratio", "Ratio", 1.0, 100.0, 10.0, Unit::Ratio), + ParamSpec::new("range", "Range", 0.0, 90.0, 60.0, Unit::Decibels), + ParamSpec::log("attack", "Attack", 0.01, 100.0, 1.0, Unit::Milliseconds), + ParamSpec::new("hold", "Hold", 0.0, 2000.0, 50.0, Unit::Milliseconds), + ParamSpec::log("release", "Release", 1.0, 3000.0, 150.0, Unit::Milliseconds), + ]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let mut s = Gate { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + threshold: 0.0, + ratio: 1.0, + range: 0.0, + att: 0.0, + rel: 0.0, + hold: 0, + hold_left: 0, + env: 0.0, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, _snap: bool) { + self.threshold = self.pv.v("threshold"); + self.ratio = self.pv.v("ratio"); + self.range = self.pv.v("range"); + self.att = coef(self.pv.v("attack"), self.sr); + self.rel = coef(self.pv.v("release"), self.sr); + self.hold = (self.pv.v("hold") * 0.001 * self.sr) as u32; + } + + pub fn curve(&self, x_db: f32) -> f32 { + expander_curve(x_db, self.threshold, self.ratio, self.range) + } +} + +impl AudioEffect for Gate { + param_plumbing!("gate"); + fn reset(&mut self) { + self.env = 0.0; + self.hold_left = 0; + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + for i in 0..n { + let peak = channels.iter().fold(0.0f32, |a, c| a.max(c[i].abs())); + let level = lin_to_db(peak); + let target = self.curve(level) - level; + if target >= self.env { + // Opening (or staying open): attack. + self.env = self.att * self.env + (1.0 - self.att) * target; + if target >= 0.0 { + self.hold_left = self.hold; + } + } else if self.hold_left > 0 { + self.hold_left -= 1; + } else { + self.env = self.rel * self.env + (1.0 - self.rel) * target; + } + let g = db_to_gain(self.env); + for ch in channels.iter_mut() { + ch[i] *= g; + } + } + } +} + +// --------------------------------------------------------------------------------------------- + +/// Look-ahead brick-wall limiter. +/// +/// The required gain for each (optionally true-peak) detected sample is held for the look-ahead +/// window with a sliding minimum, released with a one-pole, and then smoothed with a moving +/// average exactly as long as the look-ahead, so the gain has reached its target when the peak +/// leaves the delay line. A final clip at the ceiling guarantees the sample peak never exceeds it. +pub struct Limiter { + pv: ParamValues, + sr: f32, + in_gain: Smoothed, + ceiling: f32, + rel: f32, + true_peak: bool, + /// Look-ahead (samples) and detector alignment delay. + look: usize, + det_delay: usize, + bank: PolyphaseBank, + interp: Vec, + /// Per-channel delay line (ring), length `look + det_delay + 1`. + delay: Vec>, + wpos: usize, + /// Monotonic deque for the sliding minimum: (sample index, gain). + minq: VecDeque<(u64, f32)>, + /// Box filter ring and running sum. + box_ring: Vec, + box_pos: usize, + box_sum: f64, + rel_state: f32, + t: u64, + gr_meter: f32, +} + +impl Limiter { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("input_gain", "Input Boost", 0.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::new("ceiling", "Ceiling", -24.0, 0.0, -1.0, Unit::Decibels), + ParamSpec::log("release", "Release", 1.0, 1000.0, 60.0, Unit::Milliseconds), + ParamSpec::toggle("true_peak", "True Peak", true), + ]; + /// Fixed look-ahead time. + pub const LOOKAHEAD_MS: f32 = 5.0; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let channels = channels.max(1); + let bank = PolyphaseBank::new(true_peak_factor(sample_rate as f64).max(2)); + let look = ((Self::LOOKAHEAD_MS * 0.001 * sample_rate).round() as usize).max(1); + let det_delay = bank.delay(); + let dlen = look + det_delay + 1; + let mut s = Limiter { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + in_gain: Smoothed::with_ms(1.0, sample_rate, 20.0), + ceiling: 1.0, + rel: 0.0, + true_peak: true, + look, + det_delay, + bank, + interp: vec![Interpolator::default(); channels], + delay: vec![vec![0.0; dlen]; channels], + wpos: 0, + minq: VecDeque::with_capacity(look + 2), + box_ring: vec![1.0; look + 1], + box_pos: 0, + box_sum: (look + 1) as f64, + rel_state: 1.0, + t: 0, + gr_meter: 0.0, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + self.in_gain.set(db_to_gain(self.pv.v("input_gain"))); + self.ceiling = db_to_gain(self.pv.v("ceiling")); + self.rel = coef(self.pv.v("release"), self.sr); + self.true_peak = self.pv.on("true_peak"); + if snap { + self.in_gain.snap(); + } + } + + /// Current gain reduction in dB (≤ 0). + pub fn gain_reduction_db(&self) -> f32 { + self.gr_meter + } +} + +impl AudioEffect for Limiter { + param_plumbing!("limiter"); + fn latency(&self) -> usize { + self.look + self.det_delay + } + fn reset(&mut self) { + self.in_gain.snap(); + self.interp.iter_mut().for_each(Interpolator::reset); + self.delay.iter_mut().for_each(|d| d.iter_mut().for_each(|v| *v = 0.0)); + self.wpos = 0; + self.minq.clear(); + self.box_ring.iter_mut().for_each(|v| *v = 1.0); + self.box_pos = 0; + self.box_sum = self.box_ring.len() as f64; + self.rel_state = 1.0; + self.t = 0; + self.gr_meter = 0.0; + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + let nch = channels.len().min(self.delay.len()); + let dlen = self.delay[0].len(); + let window = (self.look + 1) as u64; + for i in 0..n { + let g_in = self.in_gain.tick(); + let det_pos = (self.wpos + dlen - self.det_delay) % dlen; + let out_pos = (self.wpos + 1) % dlen; // oldest = look + det_delay samples ago + let mut det = 0.0f32; + for ch in 0..nch { + let x = channels[ch][i] * g_in; + self.delay[ch][self.wpos] = x; + let aligned = self.delay[ch][det_pos].abs(); + let tp = if self.true_peak { self.interp[ch].push_peak(&self.bank, x) } else { 0.0 }; + det = det.max(aligned.max(tp)); + } + let greq = if det > self.ceiling { self.ceiling / det } else { 1.0 }; + // Sliding minimum over the last `window` values. + while let Some(&(_, v)) = self.minq.back() { + if v >= greq { + self.minq.pop_back(); + } else { + break; + } + } + self.minq.push_back((self.t, greq)); + while let Some(&(idx, _)) = self.minq.front() { + if idx + window <= self.t { + self.minq.pop_front(); + } else { + break; + } + } + let m = self.minq.front().map_or(1.0, |&(_, v)| v); + // Release. + self.rel_state = if m < self.rel_state { m } else { self.rel_state + (m - self.rel_state) * (1.0 - self.rel) }; + // Moving average over look+1 samples. + self.box_sum += self.rel_state as f64 - self.box_ring[self.box_pos] as f64; + self.box_ring[self.box_pos] = self.rel_state; + self.box_pos = (self.box_pos + 1) % self.box_ring.len(); + let g = (self.box_sum / self.box_ring.len() as f64) as f32; + for ch in 0..nch { + let y = self.delay[ch][out_pos] * g; + channels[ch][i] = y.clamp(-self.ceiling, self.ceiling); + } + self.wpos = (self.wpos + 1) % dlen; + self.t += 1; + // Periodically re-sum to cancel floating-point drift. + if self.t.is_multiple_of(65536) { + self.box_sum = self.box_ring.iter().map(|&v| v as f64).sum(); + } + self.gr_meter = lin_to_db(g); + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::loudness::LoudnessMeter; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + #[test] + fn compressor_static_curve() { + let c = |x| compressor_curve(x, -20.0, 4.0, 0.0); + assert_eq!(c(-30.0), -30.0); + assert_eq!(c(-20.0), -20.0); + assert_eq!(c(0.0), -15.0); + // Soft knee: continuous at both knee edges, midpoint below the hard curve. + let k = |x| compressor_curve(x, -20.0, 4.0, 10.0); + assert!((k(-25.0) + 25.0).abs() < 1e-5); + assert!((k(-15.0) - c(-15.0)).abs() < 1e-5); + assert!(k(-20.0) < -20.0); + assert!((k(-20.0) - (-20.0 + (0.25 - 1.0) * 25.0 / 20.0)).abs() < 1e-5); + // Monotonic. + let mut prev = f32::MIN; + for i in 0..800 { + let y = k(-80.0 + i as f32 * 0.1); + assert!(y >= prev); + prev = y; + } + } + + #[test] + fn compressor_steady_state_matches_curve() { + for detector in [0.0, 1.0] { + let mut c = Compressor::new(SR, 2); + c.set_param("threshold", -20.0); + c.set_param("ratio", 4.0); + c.set_param("knee", 0.0); + c.set_param("makeup", 3.0); + c.set_param("detector", detector); + // Constant (DC) level: peak and RMS detectors agree. + for level_db in [-40.0f32, -10.0, -2.0] { + c.reset(); + let v = db_to_gain(level_db); + let mut l = vec![v; 48000]; + let mut r = vec![-v; 48000]; + c.process(&mut [&mut l, &mut r]); + let out_db = lin_to_db(l[47999]); + let expect = compressor_curve(level_db, -20.0, 4.0, 0.0) + 3.0; + assert!((out_db - expect).abs() < 0.05, "det {detector} in {level_db}: {out_db} vs {expect}"); + assert_eq!(l[47999], -r[47999]); + } + } + } + + #[test] + fn compressor_attack_time() { + let mut c = Compressor::new(SR, 1); + c.set_param("threshold", -20.0); + c.set_param("ratio", 20.0); + c.set_param("knee", 0.0); + c.set_param("attack", 10.0); + let mut x = vec![1.0f32; 4800]; + c.process(&mut [&mut x]); + // After one time constant ~63% of the final reduction (in dB) is applied. + let final_db = compressor_curve(0.0, -20.0, 20.0, 0.0); + let at_tau = lin_to_db(x[480]); + assert!((at_tau / final_db - 0.632).abs() < 0.02, "{at_tau} / {final_db}"); + } + + #[test] + fn gate_attenuates_below_threshold_and_holds() { + let mut g = Gate::new(SR, 1); + g.set_param("threshold", -30.0); + g.set_param("ratio", 100.0); + g.set_param("range", 40.0); + g.set_param("hold", 100.0); + g.set_param("release", 10.0); + let loud = db_to_gain(-10.0); + let quiet = db_to_gain(-50.0); + let mut x: Vec = (0..48000).map(|i| if i < 24000 { loud } else { quiet }).collect(); + g.process(&mut [&mut x]); + assert!((lin_to_db(x[23999]) + 10.0).abs() < 0.01); + // During hold (100 ms) the gate stays open. + assert!((lin_to_db(x[24000 + 2400]) + 50.0).abs() < 0.01); + // Well after hold + release: attenuated by the full range. + assert!((lin_to_db(x[47999]) + 90.0).abs() < 0.1, "{}", lin_to_db(x[47999])); + assert!((expander_curve(-40.0, -30.0, 2.0, 60.0) + 50.0).abs() < 1e-6); + } + + #[test] + fn limiter_never_exceeds_ceiling() { + let mut rng = Rng::new(7); + for tp in [0.0, 1.0] { + let mut lim = Limiter::new(SR, 2); + lim.set_param("ceiling", -1.0); + lim.set_param("input_gain", 12.0); + lim.set_param("true_peak", tp); + let ceiling = db_to_gain(-1.0); + let mut all_l = Vec::new(); + let mut all_r = Vec::new(); + for block in 0..200 { + let len = 1 + (rng.next_u64() % 700) as usize; + let mut l: Vec = (0..len).map(|_| rng.uniform() * 0.8).collect(); + let mut r: Vec = (0..len).map(|_| rng.uniform() * 0.3).collect(); + if block % 7 == 0 { + l[len / 2] = 4.0; // impulse + r[0] = -3.0; + } + lim.process(&mut [&mut l, &mut r]); + assert!(l.iter().chain(&r).all(|v| v.abs() <= ceiling + 1e-7)); + all_l.extend_from_slice(&l); + all_r.extend_from_slice(&r); + } + if tp == 1.0 { + // With true-peak detection the reconstructed peak stays close to the ceiling. + let mut m = LoudnessMeter::new(SR as f64, 2); + m.process(&[&all_l, &all_r]); + assert!(m.true_peak_dbtp() < -1.0 + 0.6, "TP {}", m.true_peak_dbtp()); + } + } + } + + #[test] + fn limiter_is_transparent_below_ceiling_with_reported_latency() { + let mut lim = Limiter::new(SR, 1); + let lat = lim.latency(); + let x = sine(1000.0, 0.3, SR as f64, 4800, 0.0); + let mut y = x.clone(); + lim.process(&mut [&mut y]); + for i in lat..4800 { + assert!((y[i] - x[i - lat]).abs() < 1e-6); + } + } + + #[test] + fn limiter_reduces_loud_sine_to_ceiling_smoothly() { + let mut lim = Limiter::new(SR, 1); + lim.set_param("ceiling", -6.0); + let mut x = sine(200.0, 1.0, SR as f64, 48000, 0.0); + lim.process(&mut [&mut x]); + let tail = &x[24000..]; + let peak = tail.iter().fold(0.0f32, |a, v| a.max(v.abs())); + assert!((lin_to_db(peak) + 6.0).abs() < 0.3, "{}", lin_to_db(peak)); + } +} diff --git a/crates/audio-dsp/src/effects/eq.rs b/crates/audio-dsp/src/effects/eq.rs new file mode 100644 index 00000000..0fef89b7 --- /dev/null +++ b/crates/audio-dsp/src/effects/eq.rs @@ -0,0 +1,369 @@ +//! Parametric (8-band) and simple 3-band equalisers. + +use crate::biquad::{Biquad, Coeffs, FilterType}; +use crate::{AudioEffect, ParamSpec, ParamValues, Smoothed, Unit, block_len, param_plumbing}; + +/// Names of [`FilterType::ALL`] for the `type` choice parameters. +pub const FILTER_TYPE_NAMES: &[&str] = &["Peaking", "Low Shelf", "High Shelf", "Low Pass", "High Pass", "Notch", "Band Pass"]; + +/// Coefficients are recomputed at most every this many samples while parameters glide. +const COEFF_INTERVAL: u32 = 16; +const SMOOTH_MS: f32 = 30.0; + +#[derive(Clone, Debug)] +struct Band { + kind: FilterType, + /// log2(Hz). + freq: Smoothed, + gain: Smoothed, + q: Smoothed, + /// Band enable cross-fade (0 = bypassed, 1 = active). + mix: Smoothed, + coeffs: Coeffs, + dirty: bool, +} + +/// A bank of biquad bands shared by both EQs. +#[derive(Clone, Debug)] +pub(crate) struct BandBank { + sample_rate: f64, + bands: Vec, + /// `state[ch][band]`. + state: Vec>, + counter: u32, +} + +impl BandBank { + pub(crate) fn new(sample_rate: f32, channels: usize, nbands: usize) -> Self { + let band = Band { + kind: FilterType::Peaking, + freq: Smoothed::with_ms(1000f32.log2(), sample_rate, SMOOTH_MS), + gain: Smoothed::with_ms(0.0, sample_rate, SMOOTH_MS), + q: Smoothed::with_ms(1.0, sample_rate, SMOOTH_MS), + mix: Smoothed::with_ms(0.0, sample_rate, SMOOTH_MS), + coeffs: Coeffs::IDENTITY, + dirty: true, + }; + BandBank { sample_rate: sample_rate as f64, bands: vec![band; nbands], state: vec![vec![Biquad::new(); nbands]; channels.max(1)], counter: 0 } + } + + pub(crate) fn set_band(&mut self, i: usize, on: bool, kind: FilterType, freq: f32, gain: f32, q: f32, snap: bool) { + let b = &mut self.bands[i]; + if b.kind != kind { + b.kind = kind; + b.dirty = true; + } + b.freq.set(freq.max(1.0).log2()); + b.gain.set(gain); + b.q.set(q); + b.mix.set(if on { 1.0 } else { 0.0 }); + if snap { + b.freq.snap(); + b.gain.snap(); + b.q.snap(); + b.mix.snap(); + } + b.dirty = true; + } + + fn update_coeffs(&mut self, force: bool) { + let sr = self.sample_rate; + for b in &mut self.bands { + if b.dirty || force { + b.coeffs = Coeffs::design(b.kind, sr, 2f64.powf(b.freq.value() as f64), b.q.value() as f64, b.gain.value() as f64); + b.dirty = b.freq.is_smoothing() || b.gain.is_smoothing() || b.q.is_smoothing(); + } + } + } + + pub(crate) fn reset(&mut self) { + for b in &mut self.bands { + b.freq.snap(); + b.gain.snap(); + b.q.snap(); + b.mix.snap(); + b.dirty = true; + } + for ch in &mut self.state { + ch.iter_mut().for_each(Biquad::reset); + } + self.counter = 0; + self.update_coeffs(true); + } + + /// Analytic magnitude (dB) of the whole bank at its current (target-reached) settings. + pub(crate) fn response_db(&self, freq: f64) -> f64 { + let mut db = 0.0; + for b in &self.bands { + if b.mix.target() >= 0.5 { + let c = Coeffs::design(b.kind, self.sample_rate, 2f64.powf(b.freq.target() as f64), b.q.target() as f64, b.gain.target() as f64); + db += c.magnitude_db(freq, self.sample_rate); + } + } + db + } + + /// Process; `out_gain` is applied per sample. + pub(crate) fn process(&mut self, channels: &mut [&mut [f32]], out_gain: &mut Smoothed) { + let n = block_len(channels); + let nch = channels.len().min(self.state.len()); + for i in 0..n { + if self.counter == 0 { + self.update_coeffs(false); + } + self.counter = (self.counter + 1) % COEFF_INTERVAL; + for b in &mut self.bands { + b.freq.tick(); + b.gain.tick(); + b.q.tick(); + if b.freq.is_smoothing() || b.gain.is_smoothing() || b.q.is_smoothing() { + b.dirty = true; + } + } + let g = out_gain.tick(); + // Tick the band mixes once per sample (shared by all channels). + let mut mixes = [0.0f32; 16]; + for (m, b) in mixes.iter_mut().zip(self.bands.iter_mut()) { + *m = b.mix.tick(); + } + for ch in 0..nch { + let mut x = channels[ch][i] as f64; + for (k, b) in self.bands.iter().enumerate() { + let m = mixes[k] as f64; + let st = &mut self.state[ch][k]; + if m > 0.0 { + let y = st.tick(&b.coeffs, x); + x += m * (y - x); + } else { + st.reset(); + } + } + channels[ch][i] = (x * g as f64) as f32; + } + } + } +} + +macro_rules! eq_params { + ($(($n:literal, $on:expr, $ty:expr, $f:expr, $q:expr)),* $(,)?) => { + &[ + $( + ParamSpec::toggle(concat!("b", $n, ".on"), concat!("Band ", $n, " On"), $on), + ParamSpec::choice(concat!("b", $n, ".type"), concat!("Band ", $n, " Type"), FILTER_TYPE_NAMES, $ty), + ParamSpec::log(concat!("b", $n, ".freq"), concat!("Band ", $n, " Frequency"), 20.0, 20000.0, $f, Unit::Hertz), + ParamSpec::new(concat!("b", $n, ".gain"), concat!("Band ", $n, " Gain"), -30.0, 30.0, 0.0, Unit::Decibels), + ParamSpec::log(concat!("b", $n, ".q"), concat!("Band ", $n, " Q"), 0.1, 30.0, $q, Unit::Q), + )* + ParamSpec::new("output", "Output Gain", -30.0, 30.0, 0.0, Unit::Decibels), + ] + }; +} + +const BAND_IDS: [[&str; 5]; 8] = [ + ["b1.on", "b1.type", "b1.freq", "b1.gain", "b1.q"], + ["b2.on", "b2.type", "b2.freq", "b2.gain", "b2.q"], + ["b3.on", "b3.type", "b3.freq", "b3.gain", "b3.q"], + ["b4.on", "b4.type", "b4.freq", "b4.gain", "b4.q"], + ["b5.on", "b5.type", "b5.freq", "b5.gain", "b5.q"], + ["b6.on", "b6.type", "b6.freq", "b6.gain", "b6.q"], + ["b7.on", "b7.type", "b7.freq", "b7.gain", "b7.q"], + ["b8.on", "b8.type", "b8.freq", "b8.gain", "b8.q"], +]; + +/// 8-band parametric EQ. Every band can be any [`FilterType`]; defaults are +/// HP · low shelf · 4 × peaking · high shelf · LP with the pass filters disabled. +pub struct ParametricEq { + pv: ParamValues, + bank: BandBank, + out: Smoothed, +} + +impl ParametricEq { + pub const PARAMS: &'static [ParamSpec] = eq_params![ + ("1", false, 4, 30.0, 0.707), + ("2", true, 1, 100.0, 0.707), + ("3", true, 0, 250.0, 1.0), + ("4", true, 0, 700.0, 1.0), + ("5", true, 0, 2000.0, 1.0), + ("6", true, 0, 5000.0, 1.0), + ("7", true, 2, 10000.0, 0.707), + ("8", false, 3, 18000.0, 0.707), + ]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let mut s = ParametricEq { + pv: ParamValues::new(Self::PARAMS), + bank: BandBank::new(sample_rate, channels, 8), + out: Smoothed::with_ms(1.0, sample_rate, SMOOTH_MS), + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + for (i, ids) in BAND_IDS.iter().enumerate() { + self.bank.set_band( + i, + self.pv.on(ids[0]), + FilterType::from_index(self.pv.idx(ids[1])), + self.pv.v(ids[2]), + self.pv.v(ids[3]), + self.pv.v(ids[4]), + snap, + ); + } + self.out.set(crate::db_to_gain(self.pv.v("output"))); + if snap { + self.out.snap(); + self.bank.update_coeffs(true); + } + } + + /// Analytic magnitude response (dB) at the target settings, including output gain — for + /// drawing the EQ curve. + pub fn response_db(&self, freq: f64) -> f64 { + self.bank.response_db(freq) + self.pv.v("output") as f64 + } +} + +impl AudioEffect for ParametricEq { + param_plumbing!("parametric_eq"); + fn reset(&mut self) { + self.out.snap(); + self.bank.reset(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + self.bank.process(channels, &mut self.out); + } +} + +/// Simple 3-band EQ: low shelf, mid peak, high shelf. +pub struct SimpleEq { + pv: ParamValues, + bank: BandBank, + out: Smoothed, +} + +impl SimpleEq { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("low", "Low", -24.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::log("low_freq", "Low Frequency", 20.0, 1000.0, 200.0, Unit::Hertz), + ParamSpec::new("mid", "Mid", -24.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::log("mid_freq", "Mid Frequency", 100.0, 10000.0, 1000.0, Unit::Hertz), + ParamSpec::log("mid_q", "Mid Q", 0.1, 10.0, 0.7, Unit::Q), + ParamSpec::new("high", "High", -24.0, 24.0, 0.0, Unit::Decibels), + ParamSpec::log("high_freq", "High Frequency", 1000.0, 20000.0, 5000.0, Unit::Hertz), + ParamSpec::new("output", "Output Gain", -24.0, 24.0, 0.0, Unit::Decibels), + ]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let mut s = + SimpleEq { pv: ParamValues::new(Self::PARAMS), bank: BandBank::new(sample_rate, channels, 3), out: Smoothed::with_ms(1.0, sample_rate, SMOOTH_MS) }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + let sq = std::f32::consts::FRAC_1_SQRT_2; + let p = &self.pv; + let (lg, lf, mg, mf, mq, hg, hf) = (p.v("low"), p.v("low_freq"), p.v("mid"), p.v("mid_freq"), p.v("mid_q"), p.v("high"), p.v("high_freq")); + self.bank.set_band(0, true, FilterType::LowShelf, lf, lg, sq, snap); + self.bank.set_band(1, true, FilterType::Peaking, mf, mg, mq, snap); + self.bank.set_band(2, true, FilterType::HighShelf, hf, hg, sq, snap); + self.out.set(crate::db_to_gain(self.pv.v("output"))); + if snap { + self.out.snap(); + self.bank.update_coeffs(true); + } + } + + pub fn response_db(&self, freq: f64) -> f64 { + self.bank.response_db(freq) + self.pv.v("output") as f64 + } +} + +impl AudioEffect for SimpleEq { + param_plumbing!("simple_eq"); + fn reset(&mut self) { + self.out.snap(); + self.bank.reset(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + self.bank.process(channels, &mut self.out); + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + fn measure(fx: &mut dyn AudioEffect, freq: f64) -> f64 { + fx.reset(); + let n = SR as usize; + let mut l = sine(freq, 0.25, SR as f64, n, 0.0); + let mut r = l.clone(); + fx.process(&mut [&mut l, &mut r]); + let a = tone_amplitude(&l[n / 2..], freq, SR as f64); + let b = tone_amplitude(&r[n / 2..], freq, SR as f64); + assert!((a - b).abs() < 1e-6); + db(a / 0.25) + } + + #[test] + fn default_parametric_is_flat() { + let mut eq = ParametricEq::new(SR, 2); + for f in [30.0, 250.0, 1000.0, 8000.0, 16000.0] { + assert!(eq.response_db(f).abs() < 1e-9); + assert!(measure(&mut eq, f).abs() < 0.01); + } + } + + #[test] + fn parametric_bands_match_analytic_response() { + let mut eq = ParametricEq::new(SR, 2); + eq.set_param("b4.gain", 9.0); + eq.set_param("b4.freq", 1000.0); + eq.set_param("b4.q", 2.0); + eq.set_param("b2.gain", -6.0); + eq.set_param("b1.on", 1.0); + eq.set_param("b1.freq", 40.0); + eq.set_param("b8.on", 1.0); + eq.set_param("b8.type", 5.0); // notch + eq.set_param("b8.freq", 12000.0); + eq.set_param("b8.q", 5.0); + eq.set_param("output", 1.5); + assert!((eq.response_db(1000.0) - (9.0 + 1.5 + eq.bank.response_db(1000.0) - 9.0)).abs() < 1e-9); + for f in [30.0, 80.0, 400.0, 1000.0, 3000.0, 12000.0 * 1.2] { + let m = measure(&mut eq, f); + assert!((m - eq.response_db(f)).abs() < 0.05, "f {f}: measured {m} analytic {}", eq.response_db(f)); + } + assert!(measure(&mut eq, 12000.0) < -30.0); + } + + #[test] + fn simple_eq_gains() { + let mut eq = SimpleEq::new(SR, 2); + eq.set_param("low", 6.0); + eq.set_param("high", -6.0); + eq.set_param("mid", 3.0); + assert!((measure(&mut eq, 30.0) - eq.response_db(30.0)).abs() < 0.05); + assert!((eq.response_db(20.0) - 6.0).abs() < 0.3); + assert!((eq.response_db(20000.0) + 6.0).abs() < 0.5); + assert!((measure(&mut eq, 1000.0) - eq.response_db(1000.0)).abs() < 0.05); + } + + #[test] + fn gain_change_glides_without_jump() { + let mut eq = ParametricEq::new(SR, 1); + let mut x = vec![0.5f32; 4800]; + eq.process(&mut [&mut x[..2400]]); + eq.set_param("output", 12.0); + eq.process(&mut [&mut x[2400..]]); + let max_step = x.windows(2).map(|w| (w[1] - w[0]).abs()).fold(0.0, f32::max); + assert!(max_step < 0.01, "step {max_step}"); + assert!((x[4799] - 0.5 * crate::db_to_gain(12.0)).abs() < 1e-4); + } +} diff --git a/crates/audio-dsp/src/effects/mod.rs b/crates/audio-dsp/src/effects/mod.rs new file mode 100644 index 00000000..1be4f57b --- /dev/null +++ b/crates/audio-dsp/src/effects/mod.rs @@ -0,0 +1,267 @@ +//! Audio effects and their registry. + +pub mod channel; +pub mod dynamics; +pub mod eq; +pub mod pitch; +pub mod restoration; +mod stft; +pub mod time; + +pub use channel::{Amplify, Balance, ChannelVolume, Invert, Route, Routing}; +pub use dynamics::{Compressor, Gate, Limiter}; +pub use eq::{ParametricEq, SimpleEq}; +pub use pitch::PitchShifter; +pub use restoration::{DeHum, DeNoise}; +pub use time::{Delay, Reverb}; + +use crate::{AudioEffect, ParamSpec}; + +/// Effect category (for menus / the Effects panel's Audio Effects bin). +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub enum Category { + Filter, + Dynamics, + Time, + Channel, + Restoration, + Pitch, +} + +/// Registry entry describing one effect type. +#[derive(Clone, Copy)] +pub struct EffectInfo { + pub id: &'static str, + pub name: &'static str, + pub category: Category, + pub params: &'static [ParamSpec], + /// Construct an instance for `(sample_rate, channels)`. + pub create: fn(f32, usize) -> Box, +} + +impl std::fmt::Debug for EffectInfo { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("EffectInfo") + .field("id", &self.id) + .field("name", &self.name) + .field("category", &self.category) + .field("params", &self.params.len()) + .finish() + } +} + +macro_rules! entry { + ($id:literal, $name:literal, $cat:ident, $ty:ty) => { + EffectInfo { id: $id, name: $name, category: Category::$cat, params: <$ty>::PARAMS, create: |sr, ch| Box::new(<$ty>::new(sr, ch)) } + }; +} + +static REGISTRY: &[EffectInfo] = &[ + entry!("parametric_eq", "Parametric Equalizer", Filter, ParametricEq), + entry!("simple_eq", "Simple 3-Band EQ", Filter, SimpleEq), + entry!("compressor", "Compressor", Dynamics, Compressor), + entry!("gate", "Expander / Gate", Dynamics, Gate), + entry!("limiter", "Hard Limiter", Dynamics, Limiter), + entry!("delay", "Delay", Time, Delay), + entry!("reverb", "Reverb", Time, Reverb), + entry!("amplify", "Amplify", Channel, Amplify), + entry!("channel_volume", "Channel Volume", Channel, ChannelVolume), + entry!("balance", "Balance / Pan", Channel, Balance), + entry!("invert", "Invert", Channel, Invert), + EffectInfo { + id: "swap_channels", + name: "Swap Channels", + category: Category::Channel, + params: Route::PARAMS, + create: |_, _| Box::new(Route::new(Routing::Swap)), + }, + EffectInfo { + id: "fill_left", + name: "Fill Left with Right", + category: Category::Channel, + params: Route::PARAMS, + create: |_, _| Box::new(Route::new(Routing::FillLeft)), + }, + EffectInfo { + id: "fill_right", + name: "Fill Right with Left", + category: Category::Channel, + params: Route::PARAMS, + create: |_, _| Box::new(Route::new(Routing::FillRight)), + }, + entry!("dehum", "DeHum", Restoration, DeHum), + entry!("denoise", "DeNoise", Restoration, DeNoise), + entry!("pitch_shifter", "Pitch Shifter", Pitch, PitchShifter), +]; + +/// All registered effects. +pub fn effects() -> &'static [EffectInfo] { + REGISTRY +} + +/// Look up an effect by id. +pub fn effect_info(id: &str) -> Option<&'static EffectInfo> { + REGISTRY.iter().find(|e| e.id == id) +} + +/// Instantiate an effect by id. +pub fn create_effect(id: &str, sample_rate: f32, channels: usize) -> Option> { + effect_info(id).map(|e| (e.create)(sample_rate, channels)) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Unit; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + fn noise(seed: u64, n: usize, amp: f32) -> Vec { + let mut r = Rng::new(seed); + (0..n).map(|_| r.uniform() * amp).collect() + } + + /// A non-default setting for every parameter (to exercise real processing paths). + fn tweak(fx: &mut dyn AudioEffect) { + for p in fx.params() { + let v = match p.unit { + Unit::Toggle => 1.0, + Unit::Choice => p.max, + _ => p.min + (p.max - p.min) * 0.37, + }; + assert!(fx.set_param(p.id, v)); + } + } + + #[test] + fn registry_is_consistent() { + let mut ids = std::collections::HashSet::new(); + for e in effects() { + assert!(ids.insert(e.id), "duplicate id {}", e.id); + let fx = (e.create)(SR, 2); + assert_eq!(fx.id(), e.id); + assert_eq!(fx.params(), e.params); + let mut pids = std::collections::HashSet::new(); + for p in e.params { + assert!(pids.insert(p.id), "{}: duplicate param {}", e.id, p.id); + assert!(p.min <= p.default && p.default <= p.max, "{}.{}", e.id, p.id); + assert_eq!(fx.param(p.id), Some(p.default)); + if p.unit == Unit::Choice { + assert_eq!(p.choices.len() as f32, p.max + 1.0); + } + } + } + assert!(effect_info("reverb").is_some()); + assert!(create_effect("nope", SR, 2).is_none()); + let mut c = create_effect("compressor", SR, 2).unwrap(); + assert!(!c.set_param("nope", 1.0)); + assert!(c.set_param("ratio", 1000.0)); + assert_eq!(c.param("ratio"), Some(30.0)); + assert!(c.set_param("ratio", f32::NAN)); + assert_eq!(c.param("ratio"), Some(4.0)); + } + + #[test] + fn block_size_independence_with_param_changes() { + let n = 30000; + let src_l = noise(1, n, 0.7); + let src_r = noise(2, n, 0.7); + let change_at = 11111; + for e in effects() { + let run = |sizes: &mut dyn FnMut() -> usize| { + let mut fx = (e.create)(SR, 2); + let (mut l, mut r) = (src_l.clone(), src_r.clone()); + let mut pos = 0; + let mut changed = false; + while pos < n { + if !changed && pos == change_at { + tweak(fx.as_mut()); + changed = true; + } + let mut len = sizes().min(n - pos); + if !changed && pos + len > change_at { + len = change_at - pos; + } + fx.process(&mut [&mut l[pos..pos + len], &mut r[pos..pos + len]]); + pos += len; + } + (l, r) + }; + let big = run(&mut || usize::MAX); + let mut rng = Rng::new(5); + let small = run(&mut || match rng.next_u64() % 4 { + 0 => 0, + 1 => 1, + 2 => (rng.next_u64() % 64) as usize, + _ => (rng.next_u64() % 3000) as usize, + }); + for (a, b) in big.0.iter().chain(&big.1).zip(small.0.iter().chain(&small.1)) { + assert!((a - b).abs() <= 1e-6, "{}: block-size dependent output", e.id); + } + assert!(big.0.iter().chain(&big.1).all(|v| v.is_finite()), "{}", e.id); + } + } + + #[test] + fn zero_length_and_single_sample_blocks() { + for e in effects() { + let mut fx = (e.create)(SR, 2); + fx.process(&mut []); + fx.process(&mut [&mut [], &mut []]); + let (mut a, mut b) = ([0.25f32], [-0.25f32]); + fx.process(&mut [&mut a, &mut b]); + // Mono and more channels than configured don't panic. + fx.process(&mut [&mut [0.1f32; 3]]); + fx.process(&mut [&mut [0.1f32; 3], &mut [0.1; 3], &mut [0.1; 3]]); + assert!(a[0].is_finite() && b[0].is_finite(), "{}", e.id); + } + } + + #[test] + fn extreme_settings_stay_finite_and_denormal_free() { + for e in effects() { + for setting in [0.0f32, 1.0] { + let mut fx = (e.create)(SR, 2); + for p in e.params { + fx.set_param(p.id, if setting == 0.0 { p.min } else { p.max }); + } + let mut l = noise(3, 24000, 4.0); + let mut r = noise(4, 24000, 4.0); + l[100] = 1e6; + r[200] = -1e6; + fx.process(&mut [&mut l, &mut r]); + assert!(l.iter().chain(&r).all(|v| v.is_finite()), "{} setting {setting}", e.id); + // Long silence afterwards: output decays without denormals. + fx.reset(); + let mut l = noise(3, 4800, 0.5); + let mut r = noise(4, 4800, 0.5); + fx.process(&mut [&mut l, &mut r]); + for _ in 0..40 { + let mut l = vec![0.0f32; 48000]; + let mut r = vec![0.0f32; 48000]; + fx.process(&mut [&mut l, &mut r]); + assert!(l.iter().chain(&r).all(|v| v.is_finite() && !v.is_subnormal()), "{} setting {setting}: denormal/NaN in tail", e.id); + } + } + } + } + + #[test] + fn reset_restores_initial_behaviour() { + for e in effects() { + let mut fx = (e.create)(SR, 2); + let x = noise(8, 6000, 0.5); + let (mut a1, mut a2) = (x.clone(), x.clone()); + fx.process(&mut [&mut a1, &mut a2]); + let (mut j1, mut j2) = (noise(9, 3000, 0.9), noise(10, 3000, 0.9)); + fx.process(&mut [&mut j1, &mut j2]); + fx.reset(); + let (mut b1, mut b2) = (x.clone(), x.clone()); + fx.process(&mut [&mut b1, &mut b2]); + for (p, q) in a1.iter().zip(&b1) { + assert!((p - q).abs() < 1e-6, "{}: reset didn't clear state", e.id); + } + } + } +} diff --git a/crates/audio-dsp/src/effects/pitch.rs b/crates/audio-dsp/src/effects/pitch.rs new file mode 100644 index 00000000..fdaafd29 --- /dev/null +++ b/crates/audio-dsp/src/effects/pitch.rs @@ -0,0 +1,252 @@ +//! Pitch shifter: phase vocoder with peak-region shifting and identity phase locking. +//! +//! Every STFT frame (8× overlap) is analysed for spectral peaks; each peak's region of influence +//! (down to the magnitude minima between peaks) is moved as a block so that the peak lands on +//! `k·ratio`, preserving the window's main-lobe shape. All bins of a region receive the same phase +//! rotation, accumulated per frame from the peak's instantaneous frequency so the partial's +//! frequency is scaled exactly by `ratio` while the within-lobe phase relations stay coherent +//! (low phasiness). Duration is unchanged. At ratio 1 the output is the input delayed by the +//! frame size. + +use super::stft::{Spectral, StftChannel, fft_size}; +use crate::{AudioEffect, ParamSpec, ParamValues, Unit, block_len, param_plumbing}; +use std::f32::consts::PI; + +const OVERLAP: usize = 8; + +#[derive(Clone, Debug)] +struct PvState { + last_phase: Vec, + /// Accumulated phase rotation per output bin (from the previous frame). + rot: Vec, + new_rot: Vec, +} + +/// Frame scratch shared by all channels. +#[derive(Clone, Debug)] +struct Scratch { + mag: Vec, + phase: Vec, + out_re: Vec, + out_im: Vec, + peaks: Vec, +} + +/// Semitone pitch shifter (±12 st, fractional allowed). +pub struct PitchShifter { + pv: ParamValues, + n: usize, + hop: usize, + spec: Spectral, + stft: Vec, + state: Vec, + scratch: Scratch, + target_ratio: f32, + /// Ratio used by the current frame (glides per frame towards the target). + ratio: f32, +} + +/// Wrap a phase to (−π, π]. +#[inline] +fn wrap(p: f32) -> f32 { + p - 2.0 * PI * ((p + PI) / (2.0 * PI)).floor() +} + +impl PitchShifter { + pub const PARAMS: &'static [ParamSpec] = &[ParamSpec::new("semitones", "Pitch", -12.0, 12.0, 0.0, Unit::Semitones)]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let n = fft_size(sample_rate); + let bins = n / 2 + 1; + let ch = channels.max(1); + let mut s = PitchShifter { + pv: ParamValues::new(Self::PARAMS), + n, + hop: n / OVERLAP, + spec: Spectral::new(n, false), + stft: vec![StftChannel::new(n); ch], + state: vec![PvState { last_phase: vec![0.0; bins], rot: vec![0.0; bins], new_rot: vec![0.0; bins] }; ch], + scratch: Scratch { + mag: vec![0.0; bins], + phase: vec![0.0; bins], + out_re: vec![0.0; bins], + out_im: vec![0.0; bins], + peaks: Vec::with_capacity(bins), + }, + target_ratio: 1.0, + ratio: 1.0, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + self.target_ratio = 2f32.powf(self.pv.v("semitones") / 12.0); + if snap { + self.ratio = self.target_ratio; + } + } +} + +/// Process one analysed frame (spectrum in `spec.re/im`) into the shifted spectrum. +fn shift_frame(spec: &mut Spectral, st: &mut PvState, sc: &mut Scratch, ratio: f32, expect: f32) { + let bins = sc.mag.len(); + let mut max_mag = 0.0f32; + for k in 0..bins { + let (re, im) = (spec.re[k], spec.im[k]); + sc.mag[k] = (re * re + im * im).sqrt(); + sc.phase[k] = im.atan2(re); + max_mag = max_mag.max(sc.mag[k]); + sc.out_re[k] = 0.0; + sc.out_im[k] = 0.0; + st.new_rot[k] = 0.0; + } + // Peaks: local maxima over ±2 bins, above −100 dB of the frame maximum. + sc.peaks.clear(); + let floor = max_mag * 1e-5; + for k in 0..bins { + let m = sc.mag[k]; + if m <= floor { + continue; + } + let lo = k.saturating_sub(2); + let hi = (k + 2).min(bins - 1); + if (lo..k).all(|j| sc.mag[j] < m) && (k + 1..=hi).all(|j| sc.mag[j] <= m) { + sc.peaks.push(k); + } + } + let np = sc.peaks.len(); + let mut start = 0usize; + for i in 0..np { + let p = sc.peaks[i]; + // Region end: magnitude minimum between this peak and the next. + let end = if i + 1 < np { + let q = sc.peaks[i + 1]; + let mut mi = p + 1; + for j in p + 1..q { + if sc.mag[j] < sc.mag[mi] { + mi = j; + } + } + mi + } else { + bins + }; + // Instantaneous frequency of the peak (bins). + let d = wrap(sc.phase[p] - st.last_phase[p] - p as f32 * expect); + let f = p as f32 + d / expect; + let target = (p as f32 * ratio).round() as isize; + let delta = target - p as isize; + let prev_rot = if (0..bins as isize).contains(&target) { st.rot[target as usize] } else { 0.0 }; + let rot = wrap(prev_rot + (ratio - 1.0) * f * expect); + for k in start..end { + let t = k as isize + delta; + if t < 0 || t >= bins as isize { + continue; + } + let t = t as usize; + let (s, c) = (sc.phase[k] + rot).sin_cos(); + sc.out_re[t] += sc.mag[k] * c; + sc.out_im[t] += sc.mag[k] * s; + st.new_rot[t] = rot; + } + start = end; + } + for k in 0..bins { + st.last_phase[k] = sc.phase[k]; + spec.re[k] = sc.out_re[k]; + spec.im[k] = sc.out_im[k]; + } + std::mem::swap(&mut st.rot, &mut st.new_rot); +} + +impl AudioEffect for PitchShifter { + param_plumbing!("pitch_shifter"); + fn latency(&self) -> usize { + self.n + } + fn reset(&mut self) { + self.stft.iter_mut().for_each(StftChannel::reset); + for s in &mut self.state { + s.last_phase.iter_mut().for_each(|v| *v = 0.0); + s.rot.iter_mut().for_each(|v| *v = 0.0); + } + self.ratio = self.target_ratio; + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let len = block_len(channels); + let nch = channels.len().min(self.stft.len()); + let (n, hop) = (self.n, self.hop); + // Hann² summed over 8× overlap = 3/8 · 8 = 3. + let scale = 1.0 / (n as f32 * 3.0 / 8.0 * OVERLAP as f32); + let expect = 2.0 * PI * hop as f32 / n as f32; + // Sample-major so every channel hits its frame boundary on the same sample and the ratio + // glides identically for all channels (once per frame). + for i in 0..len { + let mut frame_done = false; + for ch in 0..nch { + let st = &mut self.state[ch]; + let spec = &mut self.spec; + let sc = &mut self.scratch; + let ratio = self.ratio; + let mut frame = |input: &[f32], out: &mut [f32]| { + spec.analyse(input); + shift_frame(spec, st, sc, ratio, expect); + spec.synthesise_add(out, scale); + frame_done = true; + }; + channels[ch][i] = self.stft[ch].tick(channels[ch][i], hop, &mut frame); + } + if frame_done { + self.ratio += (self.target_ratio - self.ratio) * 0.5; + if (self.ratio - self.target_ratio).abs() < 1e-5 { + self.ratio = self.target_ratio; + } + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + fn shifted(semis: f32, f_in: f64) -> (Vec, usize) { + let mut p = PitchShifter::new(SR, 2); + p.set_param("semitones", semis); + p.reset(); + let n = SR as usize * 2; + let mut l = sine(f_in, 0.5, SR as f64, n, 0.0); + let mut r = l.clone(); + p.process(&mut [&mut l, &mut r]); + assert_eq!(l, r); + (l, p.latency()) + } + + #[test] + fn octave_up_and_down() { + for (semis, f_in, f_out) in [(12.0, 440.0, 880.0), (-12.0, 440.0, 220.0), (7.0, 300.0, 300.0 * 1.498_307)] { + let (y, lat) = shifted(semis, f_in); + let tail = &y[lat + 4096..]; + let a_out = tone_amplitude(tail, f_out, SR as f64); + let a_in = tone_amplitude(tail, f_in, SR as f64); + assert!(db(a_out / 0.5).abs() < 3.0, "{semis}: level {} dB", db(a_out / 0.5)); + assert!(db(a_in / a_out) < -30.0, "{semis}: residual {} dB", db(a_in / a_out)); + // Output is a clean tone: almost all energy at the target frequency. + let r = rms(tail); + assert!((db(a_out / std::f64::consts::SQRT_2 / r)).abs() < 1.0); + } + } + + #[test] + fn zero_shift_is_delayed_identity() { + let (y, lat) = shifted(0.0, 440.0); + let x = sine(440.0, 0.5, SR as f64, y.len(), 0.0); + for i in lat + 2048..y.len() { + assert!((y[i] - x[i - lat]).abs() < 2e-3, "i {i}"); + } + } +} diff --git a/crates/audio-dsp/src/effects/restoration.rs b/crates/audio-dsp/src/effects/restoration.rs new file mode 100644 index 00000000..fb81006d --- /dev/null +++ b/crates/audio-dsp/src/effects/restoration.rs @@ -0,0 +1,334 @@ +//! Restoration: DeHum (harmonic notch comb) and DeNoise (STFT spectral gating). + +use super::stft::{Spectral, StftChannel, fft_size}; +use crate::biquad::{Biquad, Coeffs, FilterType}; +use crate::{AudioEffect, ParamSpec, ParamValues, Smoothed, Unit, block_len, db_to_gain, param_plumbing}; + +const MAX_HARMONICS: usize = 10; + +/// Removes mains hum: a cascade of notches at the fundamental (50/60 Hz) and its harmonics. +pub struct DeHum { + pv: ParamValues, + sr: f32, + coeffs: [Coeffs; MAX_HARMONICS], + /// Per-stage enable cross-fade (so changing the harmonic count is click-free). + stage_mix: [Smoothed; MAX_HARMONICS], + state: Vec<[Biquad; MAX_HARMONICS]>, + amount: Smoothed, +} + +impl DeHum { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::choice("frequency", "Frequency", &["50 Hz", "60 Hz"], 0), + ParamSpec::new("harmonics", "Harmonics", 1.0, MAX_HARMONICS as f32, 5.0, Unit::None), + ParamSpec::log("q", "Notch Q", 1.0, 100.0, 30.0, Unit::Q), + ParamSpec::new("amount", "Amount", 0.0, 100.0, 100.0, Unit::Percent), + ]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let mut s = DeHum { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + coeffs: [Coeffs::IDENTITY; MAX_HARMONICS], + stage_mix: [Smoothed::with_ms(0.0, sample_rate, 20.0); MAX_HARMONICS], + state: vec![[Biquad::new(); MAX_HARMONICS]; channels.max(1)], + amount: Smoothed::with_ms(1.0, sample_rate, 20.0), + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + let f0 = if self.pv.idx("frequency") == 1 { 60.0 } else { 50.0 }; + let h = self.pv.v("harmonics").round() as usize; + let q = self.pv.v("q") as f64; + for k in 0..MAX_HARMONICS { + let f = f0 * (k + 1) as f64; + let active = k < h && f < self.sr as f64 * 0.45; + self.coeffs[k] = Coeffs::design(FilterType::Notch, self.sr as f64, f, q, 0.0); + self.stage_mix[k].set(if active { 1.0 } else { 0.0 }); + if snap { + self.stage_mix[k].snap(); + } + } + self.amount.set(self.pv.v("amount") / 100.0); + if snap { + self.amount.snap(); + } + } +} + +impl AudioEffect for DeHum { + param_plumbing!("dehum"); + fn reset(&mut self) { + self.state.iter_mut().for_each(|s| s.iter_mut().for_each(Biquad::reset)); + self.stage_mix.iter_mut().for_each(Smoothed::snap); + self.amount.snap(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + let nch = channels.len().min(self.state.len()); + for i in 0..n { + let mut mix = [0.0f64; MAX_HARMONICS]; + for (m, s) in mix.iter_mut().zip(self.stage_mix.iter_mut()) { + *m = s.tick() as f64; + } + let amt = self.amount.tick() as f64; + for ch in 0..nch { + let dry = channels[ch][i] as f64; + let mut x = dry; + for k in 0..MAX_HARMONICS { + let st = &mut self.state[ch][k]; + if mix[k] > 0.0 { + let y = st.tick(&self.coeffs[k], x); + x += mix[k] * (y - x); + } else { + st.reset(); + } + } + channels[ch][i] = (dry + amt * (x - dry)) as f32; + } + } + } +} + +// --------------------------------------------------------------------------------------------- + +/// Minimum-statistics sub-windows. +const MS_SUBWINDOWS: usize = 8; +/// Length of the minimum-statistics search window (seconds). +const MS_WINDOW_S: f32 = 1.5; +/// Bias compensation of the minimum of smoothed periodograms (≈ +4 dB). +const MS_BIAS: f32 = 2.5; +const PSD_SMOOTH: f32 = 0.8; + +#[derive(Clone, Debug)] +struct NoiseState { + psd: Vec, + cur_min: Vec, + ring: Vec, + ring_pos: usize, + frames_in_sub: usize, + gains: Vec, + started: bool, +} + +impl NoiseState { + fn new(bins: usize) -> Self { + NoiseState { + psd: vec![0.0; bins], + cur_min: vec![f32::INFINITY; bins], + ring: vec![f32::INFINITY; bins * MS_SUBWINDOWS], + ring_pos: 0, + frames_in_sub: 0, + gains: vec![1.0; bins], + started: false, + } + } + fn reset(&mut self) { + self.psd.iter_mut().for_each(|v| *v = 0.0); + self.cur_min.iter_mut().for_each(|v| *v = f32::INFINITY); + self.ring.iter_mut().for_each(|v| *v = f32::INFINITY); + self.gains.iter_mut().for_each(|v| *v = 1.0); + self.ring_pos = 0; + self.frames_in_sub = 0; + self.started = false; + } +} + +/// Broadband noise reduction by spectral gating: the noise floor of every STFT bin is tracked +/// continuously with minimum statistics (so no "learn" step is needed), and bins that do not +/// rise above `floor × threshold` are attenuated by `reduction`, with per-bin attack/release +/// and light frequency smoothing to avoid musical noise. +pub struct DeNoise { + pv: ParamValues, + n: usize, + hop: usize, + sub_frames: usize, + spec: Spectral, + stft: Vec, + noise: Vec, + smoothed_gains: Vec, + floor: f32, + thr: f32, + release: f32, +} + +impl DeNoise { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("reduction", "Reduction", 0.0, 40.0, 12.0, Unit::Decibels), + ParamSpec::new("threshold", "Threshold", 0.0, 20.0, 6.0, Unit::Decibels), + ParamSpec::new("smoothing", "Smoothing", 0.0, 100.0, 50.0, Unit::Percent), + ]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let n = fft_size(sample_rate); + let hop = n / 4; + let bins = n / 2 + 1; + let sub_frames = ((MS_WINDOW_S * sample_rate / hop as f32 / MS_SUBWINDOWS as f32).round() as usize).max(1); + let ch = channels.max(1); + let mut s = DeNoise { + pv: ParamValues::new(Self::PARAMS), + n, + hop, + sub_frames, + spec: Spectral::new(n, true), + stft: vec![StftChannel::new(n); ch], + noise: vec![NoiseState::new(bins); ch], + smoothed_gains: vec![1.0; bins], + floor: 1.0, + thr: 1.0, + release: 0.5, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, _snap: bool) { + self.floor = db_to_gain(-self.pv.v("reduction")); + self.thr = 10f32.powf(self.pv.v("threshold") / 10.0); + // 0 % → fast (0.6 per frame), 100 % → slow (0.03 per frame). + self.release = 0.6 - 0.57 * self.pv.v("smoothing") / 100.0; + } +} + +impl AudioEffect for DeNoise { + param_plumbing!("denoise"); + fn latency(&self) -> usize { + self.n + } + fn reset(&mut self) { + self.stft.iter_mut().for_each(StftChannel::reset); + self.noise.iter_mut().for_each(NoiseState::reset); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let len = block_len(channels); + let nch = channels.len().min(self.stft.len()); + let (n, hop, sub_frames) = (self.n, self.hop, self.sub_frames); + let (floor, thr, release) = (self.floor, self.thr, self.release); + let scale = 1.0 / (n as f32 * 2.0); + for ch in 0..nch { + let st = &mut self.noise[ch]; + let spec = &mut self.spec; + let sg = &mut self.smoothed_gains; + let mut frame = |input: &[f32], out: &mut [f32]| { + spec.analyse(input); + let bins = n / 2 + 1; + for k in 0..bins { + let p = spec.re[k] * spec.re[k] + spec.im[k] * spec.im[k]; + st.psd[k] = if st.started { PSD_SMOOTH * st.psd[k] + (1.0 - PSD_SMOOTH) * p } else { p }; + st.cur_min[k] = st.cur_min[k].min(st.psd[k]); + } + st.started = true; + st.frames_in_sub += 1; + let rotate = st.frames_in_sub == sub_frames; + for k in 0..bins { + let mut m = st.cur_min[k]; + for u in 0..MS_SUBWINDOWS { + m = m.min(st.ring[u * bins + k]); + } + let noise = m * MS_BIAS; + let target = if st.psd[k] > noise * thr { 1.0 } else { floor }; + let g = &mut st.gains[k]; + let c = if target > *g { 0.7 } else { release }; + *g += (target - *g) * c; + if rotate { + st.ring[st.ring_pos * bins + k] = st.cur_min[k]; + st.cur_min[k] = f32::INFINITY; + } + } + if rotate { + st.ring_pos = (st.ring_pos + 1) % MS_SUBWINDOWS; + st.frames_in_sub = 0; + } + // Frequency smoothing (3-tap) to reduce musical noise. + for k in 0..bins { + let a = st.gains[k.saturating_sub(1)]; + let b = st.gains[k]; + let c = st.gains[(k + 1).min(bins - 1)]; + sg[k] = 0.25 * a + 0.5 * b + 0.25 * c; + } + for k in 0..bins { + spec.re[k] *= sg[k]; + spec.im[k] *= sg[k]; + } + spec.synthesise_add(out, scale); + }; + for v in channels[ch][..len].iter_mut() { + *v = self.stft[ch].tick(*v, hop, &mut frame); + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + #[test] + fn dehum_removes_hum_keeps_program() { + let mut d = DeHum::new(SR, 1); + d.set_param("frequency", 0.0); + d.set_param("harmonics", 4.0); + let n = SR as usize * 2; + let hum = |f: f64| sine(f, 0.2, SR as f64, n, 0.1); + let (h1, h3, p) = (hum(50.0), hum(150.0), sine(1000.0, 0.2, SR as f64, n, 0.0)); + let mut x: Vec = (0..n).map(|i| h1[i] + h3[i] + p[i]).collect(); + d.process(&mut [&mut x]); + let tail = &x[n / 2..]; + assert!(db(tone_amplitude(tail, 50.0, SR as f64) / 0.2) < -30.0); + assert!(db(tone_amplitude(tail, 150.0, SR as f64) / 0.2) < -30.0); + assert!(db(tone_amplitude(tail, 1000.0, SR as f64) / 0.2).abs() < 0.2); + // 60 Hz mode leaves 50 Hz mostly alone. + let mut d60 = DeHum::new(SR, 1); + d60.set_param("frequency", 1.0); + d60.reset(); + let mut y = hum(50.0); + d60.process(&mut [&mut y]); + assert!(db(tone_amplitude(&y[n / 2..], 50.0, SR as f64) / 0.2) > -3.0); + } + + #[test] + fn denoise_zero_reduction_is_delayed_identity() { + let mut d = DeNoise::new(SR, 1); + d.set_param("reduction", 0.0); + let lat = d.latency(); + let mut rng = Rng::new(9); + let x: Vec = (0..20000).map(|_| rng.uniform() * 0.5).collect(); + let mut y = x.clone(); + d.process(&mut [&mut y]); + for i in lat..x.len() { + assert!((y[i] - x[i - lat]).abs() < 1e-4, "i {i}: {} vs {}", y[i], x[i - lat]); + } + } + + #[test] + fn denoise_reduces_noise_and_keeps_bursts() { + let mut d = DeNoise::new(SR, 1); + d.set_param("reduction", 20.0); + let lat = d.latency(); + let n = SR as usize * 8; + let mut rng = Rng::new(21); + let noise: Vec = (0..n).map(|_| rng.uniform() * 0.02).collect(); + let tone = sine(1000.0, 0.3, SR as f64, n, 0.0); + // Tone bursts: on for 0.5 s, off for 0.5 s. + let on = |i: usize| (i / (SR as usize / 2)) % 2 == 1; + let x: Vec = (0..n).map(|i| noise[i] + if on(i) { tone[i] } else { 0.0 }).collect(); + let mut y = x.clone(); + d.process(&mut [&mut y]); + let half = SR as usize / 2; + // Evaluate the last seconds (after warm-up), aligned for latency, away from edges. + let seg = |k: usize| { + let s = k * half + lat + half / 4; + s..s + half / 2 + }; + let (off, onn) = (seg(14), seg(15)); + let red = db(rms(&y[off.clone()]) / rms(&x[off.start - lat..off.end - lat])); + assert!(red < -10.0, "noise reduction only {red} dB"); + let kept = db(tone_amplitude(&y[onn.clone()], 1000.0, SR as f64) / 0.3); + assert!(kept.abs() < 1.0, "tone changed by {kept} dB"); + } +} diff --git a/crates/audio-dsp/src/effects/stft.rs b/crates/audio-dsp/src/effects/stft.rs new file mode 100644 index 00000000..b5e35e61 --- /dev/null +++ b/crates/audio-dsp/src/effects/stft.rs @@ -0,0 +1,99 @@ +//! Streaming STFT overlap-add framework shared by the spectral effects. + +use crate::fft::Fft; +use std::f64::consts::PI; + +/// Per-channel streaming buffers. Latency is exactly `n` samples. +#[derive(Clone, Debug)] +pub(crate) struct StftChannel { + input: Vec, + output: Vec, + pos: usize, +} + +impl StftChannel { + pub(crate) fn new(n: usize) -> Self { + StftChannel { input: vec![0.0; n], output: vec![0.0; n], pos: 0 } + } + pub(crate) fn reset(&mut self) { + self.input.iter_mut().for_each(|v| *v = 0.0); + self.output.iter_mut().for_each(|v| *v = 0.0); + self.pos = 0; + } + /// Push one sample and pop one output sample. Every `hop` samples `frame` is called with the + /// last `n` input samples and must write an `n`-sample (windowed) frame that is overlap-added. + #[inline] + pub(crate) fn tick(&mut self, x: f32, hop: usize, frame: &mut impl FnMut(&[f32], &mut [f32])) -> f32 { + let n = self.input.len(); + self.input[n - hop + self.pos] = x; + let y = self.output[self.pos]; + self.pos += 1; + if self.pos == hop { + self.pos = 0; + self.output.copy_within(hop.., 0); + self.output[n - hop..].iter_mut().for_each(|v| *v = 0.0); + frame(&self.input, &mut self.output); + self.input.copy_within(hop.., 0); + } + y + } +} + +/// Periodic Hann window. +pub(crate) fn hann(n: usize) -> Vec { + (0..n).map(|i| (0.5 - 0.5 * (2.0 * PI * i as f64 / n as f64).cos()) as f32).collect() +} + +/// FFT size for spectral effects at a sample rate (~43 ms at 48 kHz). +pub(crate) fn fft_size(sample_rate: f32) -> usize { + if sample_rate <= 50_000.0 { + 2048 + } else if sample_rate <= 100_000.0 { + 4096 + } else { + 8192 + } +} + +/// Shared spectral scratch. +#[derive(Clone, Debug)] +pub(crate) struct Spectral { + pub fft: Fft, + pub window: Vec, + pub re: Vec, + pub im: Vec, +} + +impl Spectral { + /// `sqrt_window`: use √Hann for analysis and synthesis (sums to a constant at 75 % overlap); + /// otherwise Hann for both. + pub(crate) fn new(n: usize, sqrt_window: bool) -> Self { + let mut window = hann(n); + if sqrt_window { + window.iter_mut().for_each(|w| *w = w.sqrt()); + } + Spectral { fft: Fft::new(n), window, re: vec![0.0; n], im: vec![0.0; n] } + } + /// Window `input` into re/im and run the forward FFT. + pub(crate) fn analyse(&mut self, input: &[f32]) { + for i in 0..input.len() { + self.re[i] = input[i] * self.window[i]; + self.im[i] = 0.0; + } + self.fft.forward(&mut self.re, &mut self.im); + } + /// Enforce Hermitian symmetry from bins 0..=n/2, inverse FFT, window, scale and add to `out`. + pub(crate) fn synthesise_add(&mut self, out: &mut [f32], scale: f32) { + let n = self.re.len(); + self.im[0] = 0.0; + self.im[n / 2] = 0.0; + for k in 1..n / 2 { + self.re[n - k] = self.re[k]; + self.im[n - k] = -self.im[k]; + } + self.fft.inverse(&mut self.re, &mut self.im); + for i in 0..n { + out[i] += self.re[i] * self.window[i] * scale; + } + } +} diff --git a/crates/audio-dsp/src/effects/time.rs b/crates/audio-dsp/src/effects/time.rs new file mode 100644 index 00000000..d10db23d --- /dev/null +++ b/crates/audio-dsp/src/effects/time.rs @@ -0,0 +1,372 @@ +//! Time-based effects: feedback delay and an algorithmic (FDN) reverb. + +use crate::{AudioEffect, ParamSpec, ParamValues, Smoothed, Unit, block_len, flush32, param_plumbing}; + +/// Fractional read from a circular buffer: `delay` samples behind write position `w`. +#[inline] +fn read_frac(buf: &[f32], w: usize, delay: f32) -> f32 { + let len = buf.len(); + let d = delay.clamp(1.0, (len - 2) as f32); + let di = d.floor(); + let frac = d - di; + let i0 = (w + len - di as usize) % len; + let i1 = (i0 + len - 1) % len; + buf[i0] + (buf[i1] - buf[i0]) * frac +} + +/// Note divisions for tempo-synced delay, in quarter notes. +const DIVISIONS: &[&str] = &["1/1", "1/2", "1/4", "1/8", "1/16", "1/4 dotted", "1/8 dotted", "1/4 triplet", "1/8 triplet"]; +const DIVISION_QUARTERS: [f32; 9] = [4.0, 2.0, 1.0, 0.5, 0.25, 1.5, 0.75, 2.0 / 3.0, 1.0 / 3.0]; + +/// Delay time in ms for a tempo and division index. +pub fn tempo_delay_ms(bpm: f32, division: usize) -> f32 { + 60_000.0 / bpm.max(1.0) * DIVISION_QUARTERS[division.min(DIVISION_QUARTERS.len() - 1)] +} + +/// Feedback delay (per-channel), time in ms or tempo-synced. +pub struct Delay { + pv: ParamValues, + sr: f32, + bufs: Vec>, + w: usize, + /// Delay in samples (smoothed: glides like a tape delay instead of clicking). + time: Smoothed, + feedback: Smoothed, + mix: Smoothed, +} + +impl Delay { + pub const MAX_MS: f32 = 4000.0; + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::choice("mode", "Mode", &["Time", "Tempo"], 0), + ParamSpec::log("time", "Delay Time", 1.0, 2000.0, 250.0, Unit::Milliseconds), + ParamSpec::new("bpm", "Tempo", 20.0, 300.0, 120.0, Unit::Bpm), + ParamSpec::choice("division", "Division", DIVISIONS, 3), + ParamSpec::new("feedback", "Feedback", 0.0, 95.0, 30.0, Unit::Percent), + ParamSpec::new("mix", "Mix", 0.0, 100.0, 30.0, Unit::Percent), + ]; + + pub fn new(sample_rate: f32, channels: usize) -> Self { + let len = (Self::MAX_MS * 0.001 * sample_rate) as usize + 4; + let mut s = Delay { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + bufs: vec![vec![0.0; len]; channels.max(1)], + w: 0, + time: Smoothed::with_ms(1.0, sample_rate, 100.0), + feedback: Smoothed::with_ms(0.0, sample_rate, 20.0), + mix: Smoothed::with_ms(0.0, sample_rate, 20.0), + }; + s.apply_params(true); + s + } + + /// Effective delay time (ms) at the current settings. + pub fn delay_ms(&self) -> f32 { + let ms = if self.pv.idx("mode") == 1 { tempo_delay_ms(self.pv.v("bpm"), self.pv.idx("division")) } else { self.pv.v("time") }; + ms.min(Self::MAX_MS) + } + + fn apply_params(&mut self, snap: bool) { + self.time.set(self.delay_ms() * 0.001 * self.sr); + self.feedback.set(self.pv.v("feedback") / 100.0); + self.mix.set(self.pv.v("mix") / 100.0); + if snap { + self.time.snap(); + self.feedback.snap(); + self.mix.snap(); + } + } +} + +impl AudioEffect for Delay { + param_plumbing!("delay"); + fn reset(&mut self) { + self.bufs.iter_mut().for_each(|b| b.iter_mut().for_each(|v| *v = 0.0)); + self.w = 0; + self.time.snap(); + self.feedback.snap(); + self.mix.snap(); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + let nch = channels.len().min(self.bufs.len()); + let len = self.bufs[0].len(); + for i in 0..n { + let d = self.time.tick(); + let fb = self.feedback.tick(); + let mix = self.mix.tick(); + for ch in 0..nch { + let x = channels[ch][i]; + let buf = &mut self.bufs[ch]; + let yd = read_frac(buf, self.w, d); + buf[self.w] = flush32(x + fb * yd); + channels[ch][i] = x + mix * (yd - x); + } + self.w = (self.w + 1) % len; + } + } +} + +// --------------------------------------------------------------------------------------------- + +const FDN: usize = 8; +/// Base loop delays (ms) at size = 50 %; mutually prime-ish so modes don't pile up. +const BASE_MS: [f32; FDN] = [29.7, 37.1, 41.1, 43.7, 53.9, 59.3, 67.1, 73.3]; +const MAX_SIZE_SCALE: f32 = 2.0; +/// Output tap signs (left, right) for decorrelated stereo. +const SIGN_L: [f32; FDN] = [1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0]; +const SIGN_R: [f32; FDN] = [1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0]; + +fn size_scale(size_pct: f32) -> f32 { + 0.35 + (MAX_SIZE_SCALE - 0.35) * size_pct / 100.0 +} + +/// 8×8 feedback-delay-network reverb with a Householder feedback matrix, per-line RT60 gains, +/// one-pole HF damping in each loop, pre-delay, size and dry/wet mix. +pub struct Reverb { + pv: ParamValues, + sr: f32, + lines: Vec>, + w: usize, + len: [Smoothed; FDN], + lp: [f32; FDN], + gains: [f32; FDN], + damp: Smoothed, + rt60: Smoothed, + mix: Smoothed, + pre: Vec, + pre_w: usize, + pre_len: Smoothed, + counter: u32, +} + +impl Reverb { + pub const PARAMS: &'static [ParamSpec] = &[ + ParamSpec::new("predelay", "Pre-delay", 0.0, 200.0, 20.0, Unit::Milliseconds), + ParamSpec::log("decay", "Decay (RT60)", 0.1, 20.0, 2.0, Unit::Seconds), + ParamSpec::new("damping", "HF Damping", 0.0, 100.0, 40.0, Unit::Percent), + ParamSpec::new("size", "Room Size", 0.0, 100.0, 50.0, Unit::Percent), + ParamSpec::new("mix", "Mix", 0.0, 100.0, 25.0, Unit::Percent), + ]; + + pub fn new(sample_rate: f32, _channels: usize) -> Self { + let max_line = (BASE_MS[FDN - 1] * MAX_SIZE_SCALE * 0.001 * sample_rate) as usize + 4; + let pre_len = (0.2 * sample_rate) as usize + 4; + let sm = |v| Smoothed::with_ms(v, sample_rate, 50.0); + let mut s = Reverb { + pv: ParamValues::new(Self::PARAMS), + sr: sample_rate, + lines: vec![vec![0.0; max_line]; FDN], + w: 0, + len: [sm(1.0); FDN], + lp: [0.0; FDN], + gains: [0.0; FDN], + damp: sm(0.0), + rt60: sm(1.0), + mix: sm(0.0), + pre: vec![0.0; pre_len], + pre_w: 0, + pre_len: Smoothed::with_ms(1.0, sample_rate, 100.0), + counter: 0, + }; + s.apply_params(true); + s + } + + fn apply_params(&mut self, snap: bool) { + let scale = size_scale(self.pv.v("size")); + for (k, l) in self.len.iter_mut().enumerate() { + l.set(BASE_MS[k] * scale * 0.001 * self.sr); + } + self.pre_len.set((self.pv.v("predelay") * 0.001 * self.sr).max(1.0)); + self.damp.set(self.pv.v("damping") / 100.0 * 0.85); + self.rt60.set(self.pv.v("decay")); + self.mix.set(self.pv.v("mix") / 100.0); + if snap { + self.len.iter_mut().for_each(Smoothed::snap); + self.pre_len.snap(); + self.damp.snap(); + self.rt60.snap(); + self.mix.snap(); + self.update_gains(); + } + } + + fn update_gains(&mut self) { + let rt = self.rt60.value().max(0.01); + for k in 0..FDN { + // Attenuation per pass so the loop decays 60 dB in rt seconds. + self.gains[k] = 10f32.powf(-3.0 * self.len[k].value() / (rt * self.sr)); + } + } +} + +impl AudioEffect for Reverb { + param_plumbing!("reverb"); + fn reset(&mut self) { + self.lines.iter_mut().for_each(|l| l.iter_mut().for_each(|v| *v = 0.0)); + self.pre.iter_mut().for_each(|v| *v = 0.0); + self.lp = [0.0; FDN]; + self.w = 0; + self.pre_w = 0; + self.counter = 0; + self.apply_params(true); + } + fn process(&mut self, channels: &mut [&mut [f32]]) { + let n = block_len(channels); + let nch = channels.len(); + if nch == 0 { + return; + } + let line_len = self.lines[0].len(); + let pre_len = self.pre.len(); + let hh = 2.0 / FDN as f32; + for i in 0..n { + if self.counter == 0 { + self.update_gains(); + } + self.counter = (self.counter + 1) % 32; + for l in &mut self.len { + l.tick(); + } + self.rt60.tick(); + let damp = self.damp.tick(); + let mix = self.mix.tick(); + let pd = self.pre_len.tick(); + + let mut mono = 0.0; + for ch in channels.iter() { + mono += ch[i]; + } + mono /= nch as f32; + self.pre[self.pre_w] = mono; + let input = read_frac(&self.pre, self.pre_w, pd); + self.pre_w = (self.pre_w + 1) % pre_len; + + let mut v = [0.0f32; FDN]; + let mut sum = 0.0; + for k in 0..FDN { + let d = read_frac(&self.lines[k], self.w, self.len[k].value()); + self.lp[k] = flush32(d + (self.lp[k] - d) * damp); + v[k] = self.lp[k] * self.gains[k]; + sum += v[k]; + } + let (mut wl, mut wr) = (0.0, 0.0); + for k in 0..FDN { + let fbk = v[k] - hh * sum; + self.lines[k][self.w] = flush32(fbk + input * 0.5); + wl += SIGN_L[k] * v[k]; + wr += SIGN_R[k] * v[k]; + } + self.w = (self.w + 1) % line_len; + let (wl, wr) = (wl * 0.35, wr * 0.35); + for (c, ch) in channels.iter_mut().enumerate() { + let x = ch[i]; + let wet = if nch == 1 { + 0.5 * (wl + wr) + } else if c % 2 == 0 { + wl + } else { + wr + }; + ch[i] = x + mix * (wet - x); + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + const SR: f32 = 48000.0; + + #[test] + fn delay_impulse_and_feedback() { + let mut d = Delay::new(SR, 1); + d.set_param("time", 100.0); + d.set_param("feedback", 50.0); + d.set_param("mix", 100.0); + d.reset(); + let mut x = vec![0.0f32; 24000]; + x[0] = 1.0; + d.process(&mut [&mut x]); + assert!((x[4800] - 1.0).abs() < 1e-6); + assert!((x[9600] - 0.5).abs() < 1e-6); + assert!((x[14400] - 0.25).abs() < 1e-6); + assert_eq!(x[0], 0.0); + } + + #[test] + fn tempo_sync() { + assert_eq!(tempo_delay_ms(120.0, 2), 500.0); + assert_eq!(tempo_delay_ms(120.0, 3), 250.0); + assert!((tempo_delay_ms(120.0, 5) - 750.0).abs() < 1e-3); + let mut d = Delay::new(SR, 1); + d.set_param("mode", 1.0); + d.set_param("bpm", 90.0); + d.set_param("division", 2.0); + assert!((d.delay_ms() - 666.666).abs() < 0.01); + d.set_param("bpm", 20.0); + d.set_param("division", 0.0); + assert_eq!(d.delay_ms(), Delay::MAX_MS); + } + + #[test] + fn reverb_decays_and_is_bounded() { + let mut r = Reverb::new(SR, 2); + r.set_param("mix", 100.0); + r.set_param("decay", 1.0); + r.set_param("damping", 0.0); + r.reset(); + let n = (SR * 3.0) as usize; + let mut l = vec![0.0f32; n]; + let mut rr = vec![0.0f32; n]; + l[0] = 1.0; + rr[0] = 1.0; + r.process(&mut [&mut l, &mut rr]); + assert!(l.iter().chain(&rr).all(|v| v.is_finite() && v.abs() < 1.0)); + let win = (SR * 0.1) as usize; + let e = |s: usize| rms(&l[s..s + win]); + let e1 = e((SR * 0.2) as usize); + let e2 = e((SR * 1.2) as usize); + // RT60 = 1 s → ~60 dB drop per second (allow generous tolerance). + let drop = db(e1 / e2); + assert!((45.0..75.0).contains(&drop), "decay over 1 s = {drop} dB"); + // Stereo outputs are decorrelated but both carry energy. + assert!(rms(&rr[..win * 5]) > 0.0); + assert!(l.iter().zip(&rr).any(|(a, b)| (a - b).abs() > 1e-4)); + } + + #[test] + fn reverb_stable_with_longest_decay_and_noise() { + let mut r = Reverb::new(SR, 2); + r.set_param("decay", 20.0); + r.set_param("size", 100.0); + r.set_param("damping", 0.0); + r.set_param("mix", 100.0); + let mut rng = Rng::new(3); + let mut peak = 0.0f32; + for _ in 0..100 { + let mut l: Vec = (0..4800).map(|_| rng.uniform()).collect(); + let mut rr = l.clone(); + r.process(&mut [&mut l, &mut rr]); + peak = l.iter().chain(&rr).fold(peak, |a, v| a.max(v.abs())); + } + assert!(peak.is_finite() && peak < 20.0, "peak {peak}"); + // Then silence: energy decays. + let mut last = 0.0; + for k in 0..30 { + let mut l = vec![0.0f32; 48000]; + let mut rr = vec![0.0f32; 48000]; + r.process(&mut [&mut l, &mut rr]); + let e = rms(&l); + if k > 0 { + assert!(e < last, "energy must decay"); + } + last = e; + } + } +} diff --git a/crates/audio-dsp/src/fft.rs b/crates/audio-dsp/src/fft.rs new file mode 100644 index 00000000..5f7510b6 --- /dev/null +++ b/crates/audio-dsp/src/fft.rs @@ -0,0 +1,115 @@ +//! Minimal in-place radix-2 complex FFT (power-of-two sizes), allocation-free after planning. +//! +//! Kept in-crate (rather than depending on an FFT library) so the crate stays dependency-free; +//! the STFT effects only need a few fixed power-of-two sizes. + +use std::f64::consts::PI; + +/// A planned FFT of fixed power-of-two size. +#[derive(Clone, Debug)] +pub struct Fft { + n: usize, + /// e^{-2πik/n} for k in 0..n/2. + tw_re: Vec, + tw_im: Vec, + rev: Vec, +} + +impl Fft { + /// Plan an FFT of size `n` (must be a power of two ≥ 2). + pub fn new(n: usize) -> Fft { + assert!(n >= 2 && n.is_power_of_two(), "FFT size must be a power of two"); + let bits = n.trailing_zeros(); + let rev = (0..n as u32).map(|i| i.reverse_bits() >> (32 - bits)).collect(); + let (tw_re, tw_im) = (0..n / 2) + .map(|k| { + let a = -2.0 * PI * k as f64 / n as f64; + (a.cos() as f32, a.sin() as f32) + }) + .unzip(); + Fft { n, tw_re, tw_im, rev } + } + + pub fn len(&self) -> usize { + self.n + } + + pub fn is_empty(&self) -> bool { + false + } + + /// Forward transform (e^{-iωt}), unscaled. + pub fn forward(&self, re: &mut [f32], im: &mut [f32]) { + self.run(re, im, false); + } + + /// Inverse transform (e^{+iωt}), unscaled (divide by `n` yourself). + pub fn inverse(&self, re: &mut [f32], im: &mut [f32]) { + self.run(re, im, true); + } + + fn run(&self, re: &mut [f32], im: &mut [f32], inverse: bool) { + let n = self.n; + assert!(re.len() >= n && im.len() >= n); + for i in 0..n { + let j = self.rev[i] as usize; + if j > i { + re.swap(i, j); + im.swap(i, j); + } + } + let sign = if inverse { -1.0 } else { 1.0 }; + let mut len = 2; + while len <= n { + let half = len / 2; + let stride = n / len; + let mut start = 0; + while start < n { + for k in 0..half { + let wr = self.tw_re[k * stride]; + let wi = sign * self.tw_im[k * stride]; + let a = start + k; + let b = a + half; + let xr = re[b] * wr - im[b] * wi; + let xi = re[b] * wi + im[b] * wr; + re[b] = re[a] - xr; + im[b] = im[a] - xi; + re[a] += xr; + im[a] += xi; + } + start += len; + } + len <<= 1; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn matches_naive_dft_and_round_trips() { + let n = 64; + let fft = Fft::new(n); + let x: Vec = (0..n).map(|i| ((i * 7 % 13) as f32 - 6.0) / 6.0).collect(); + let mut re = x.clone(); + let mut im = vec![0.0; n]; + fft.forward(&mut re, &mut im); + for k in 0..n { + let (mut r, mut i) = (0.0f64, 0.0f64); + for (t, &v) in x.iter().enumerate() { + let a = -2.0 * PI * (k * t) as f64 / n as f64; + r += v as f64 * a.cos(); + i += v as f64 * a.sin(); + } + assert!((re[k] as f64 - r).abs() < 1e-4, "re {k}"); + assert!((im[k] as f64 - i).abs() < 1e-4, "im {k}"); + } + fft.inverse(&mut re, &mut im); + for t in 0..n { + assert!((re[t] / n as f32 - x[t]).abs() < 1e-5); + assert!((im[t] / n as f32).abs() < 1e-5); + } + } +} diff --git a/crates/audio-dsp/src/lib.rs b/crates/audio-dsp/src/lib.rs new file mode 100644 index 00000000..15edf5b0 --- /dev/null +++ b/crates/audio-dsp/src/lib.rs @@ -0,0 +1,272 @@ +//! Real-time audio DSP for FilmCraft (layer L1, dependency-free, `wasm32` clean). +//! +//! * [`loudness`] — ITU-R BS.1770-4 / EBU R128 loudness meter (momentary, short-term, +//! integrated, loudness range, sample peak, true peak) and the loudness-normalisation helper. +//! * [`effects`] — clip/track audio effects behind the common [`AudioEffect`] trait, with a +//! registry ([`effects()`]) describing every effect's parameters (range, unit, default) so the +//! engine and UI can build controls generically. +//! * [`biquad`] — RBJ-cookbook biquads (TDF-II) with analytic magnitude response. +//! +//! Conventions: audio is **planar f32** (`&mut [&mut [f32]]`, one slice per channel, all the same +//! length). `process` never allocates; all buffers are sized at construction. Parameters may be +//! changed between blocks and are smoothed per sample, so the output is independent of how the +//! stream is cut into blocks. Recursive state is flushed to zero when it decays below ~1e-30 so +//! denormals never appear (no FTZ/DAZ CPU flags needed, which keeps the crate `unsafe`-free and +//! portable to wasm). + +pub mod biquad; +pub mod effects; +pub mod fft; +pub mod loudness; +pub mod oversample; +mod smooth; + +pub use effects::{Category, EffectInfo, create_effect, effect_info, effects}; +pub use loudness::{LoudnessMeter, normalize_gain_db}; +pub use smooth::Smoothed; + +/// Unit of a parameter value (for display and control construction). +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub enum Unit { + /// Unitless scalar. + None, + Decibels, + Hertz, + Milliseconds, + Seconds, + Percent, + /// Filter quality factor. + Q, + /// Compression ratio (`x:1`). + Ratio, + Semitones, + Bpm, + /// Pan / balance position, −100 (left) … +100 (right). + Pan, + /// Boolean toggle: 0 = off, 1 = on. + Toggle, + /// Index into [`ParamSpec::choices`]. + Choice, +} + +/// Static description of one effect parameter. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ParamSpec { + /// Stable id used by `set_param` / automation (e.g. `"threshold"`, `"b3.freq"`). + pub id: &'static str, + /// Human-readable label. + pub name: &'static str, + pub min: f32, + pub max: f32, + pub default: f32, + pub unit: Unit, + /// Hint that the control should be logarithmic (frequencies, times). + pub log_scale: bool, + /// Option labels when `unit == Unit::Choice` (value is the index). + pub choices: &'static [&'static str], +} + +impl ParamSpec { + /// A continuous (linear) parameter. + pub const fn new(id: &'static str, name: &'static str, min: f32, max: f32, default: f32, unit: Unit) -> Self { + ParamSpec { id, name, min, max, default, unit, log_scale: false, choices: &[] } + } + /// A continuous parameter best shown on a logarithmic control. + pub const fn log(id: &'static str, name: &'static str, min: f32, max: f32, default: f32, unit: Unit) -> Self { + ParamSpec { id, name, min, max, default, unit, log_scale: true, choices: &[] } + } + /// An on/off toggle. + pub const fn toggle(id: &'static str, name: &'static str, default: bool) -> Self { + ParamSpec { id, name, min: 0.0, max: 1.0, default: if default { 1.0 } else { 0.0 }, unit: Unit::Toggle, log_scale: false, choices: &[] } + } + /// A discrete choice; the value is the index into `choices`. + pub const fn choice(id: &'static str, name: &'static str, choices: &'static [&'static str], default: usize) -> Self { + ParamSpec { id, name, min: 0.0, max: (choices.len() - 1) as f32, default: default as f32, unit: Unit::Choice, log_scale: false, choices } + } + /// Clamp (and for choices/toggles, round) a value into this parameter's domain. + /// NaN maps to the default. + pub fn sanitize(&self, v: f32) -> f32 { + if v.is_nan() { + return self.default; + } + let v = v.clamp(self.min, self.max); + match self.unit { + Unit::Choice | Unit::Toggle => v.round(), + _ => v, + } + } +} + +/// A real-time audio effect operating in place on planar channels. +/// +/// Implementations are constructed (allocating) with a sample rate and channel count; after that +/// `process`, `set_param` and `reset` never allocate. If `process` receives more channels than the +/// effect was built for, the extra channels are passed through untouched (channel-generic effects +/// such as gain process all of them). +pub trait AudioEffect: Send { + /// Registry id (e.g. `"compressor"`). + fn id(&self) -> &'static str; + /// Parameter descriptions (same as the registry entry's). + fn params(&self) -> &'static [ParamSpec]; + /// Set a parameter by id; the value is clamped to its range. Returns `false` for unknown ids. + /// Takes effect smoothly from the next processed sample. + fn set_param(&mut self, id: &str, value: f32) -> bool; + /// Current (target) value of a parameter. + fn param(&self, id: &str) -> Option; + /// Clear all internal state (delay lines, envelopes, filter memories); smoothed parameters + /// jump to their targets. + fn reset(&mut self); + /// Process one block in place. All channel slices must have the same length (the shortest + /// length is used otherwise). Zero-length blocks are allowed. + fn process(&mut self, channels: &mut [&mut [f32]]); + /// Processing latency in samples (look-ahead / STFT delay) for delay compensation. + fn latency(&self) -> usize { + 0 + } +} + +/// Parameter value storage shared by the effect implementations. +#[derive(Clone, Debug)] +pub(crate) struct ParamValues { + specs: &'static [ParamSpec], + values: Vec, +} + +impl ParamValues { + pub(crate) fn new(specs: &'static [ParamSpec]) -> Self { + ParamValues { specs, values: specs.iter().map(|s| s.default).collect() } + } + pub(crate) fn index_of(&self, id: &str) -> Option { + self.specs.iter().position(|s| s.id == id) + } + /// Returns `true` if the id exists. + pub(crate) fn set(&mut self, id: &str, v: f32) -> bool { + match self.index_of(id) { + Some(i) => { + self.values[i] = self.specs[i].sanitize(v); + true + } + None => false, + } + } + pub(crate) fn get(&self, id: &str) -> Option { + self.index_of(id).map(|i| self.values[i]) + } + /// Value of a parameter that is known to exist (panics otherwise — programming error). + pub(crate) fn v(&self, id: &str) -> f32 { + match self.get(id) { + Some(v) => v, + None => panic!("unknown parameter {id}"), + } + } + pub(crate) fn on(&self, id: &str) -> bool { + self.v(id) >= 0.5 + } + pub(crate) fn idx(&self, id: &str) -> usize { + self.v(id).max(0.0) as usize + } +} + +/// Generates the id/params/set_param/param methods of [`AudioEffect`] for a struct with a +/// `pv: ParamValues` field, a `PARAMS` const and an `apply_params(&mut self, snap: bool)` method. +macro_rules! param_plumbing { + ($id:expr) => { + fn id(&self) -> &'static str { + $id + } + fn params(&self) -> &'static [$crate::ParamSpec] { + Self::PARAMS + } + fn set_param(&mut self, id: &str, value: f32) -> bool { + if self.pv.set(id, value) { + self.apply_params(false); + true + } else { + false + } + } + fn param(&self, id: &str) -> Option { + self.pv.get(id) + } + }; +} +pub(crate) use param_plumbing; + +/// Decibels → linear amplitude. +#[inline] +pub fn db_to_gain(db: f32) -> f32 { + 10f32.powf(db / 20.0) +} + +/// Linear amplitude → decibels (−inf for 0). +#[inline] +pub fn gain_to_db(g: f32) -> f32 { + 20.0 * g.abs().log10() +} + +/// Flush tiny values to zero (denormal protection for recursive state). +#[inline(always)] +pub(crate) fn flush(x: f64) -> f64 { + if x.abs() < 1e-30 { 0.0 } else { x } +} + +/// f32 variant of [`flush`] for f32 delay lines. +#[inline(always)] +pub(crate) fn flush32(x: f32) -> f32 { + if x.abs() < 1e-25 { 0.0 } else { x } +} + +/// Common block length of planar channels. +#[inline] +pub(crate) fn block_len(channels: &[&mut [f32]]) -> usize { + channels.iter().map(|c| c.len()).min().unwrap_or(0) +} + +#[cfg(test)] +pub(crate) mod testutil { + /// Tiny deterministic PRNG (xorshift64*) for tests. + pub struct Rng(u64); + impl Rng { + pub fn new(seed: u64) -> Self { + Rng(seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1) + } + pub fn next_u64(&mut self) -> u64 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + x.wrapping_mul(0x2545_F491_4F6C_DD1D) + } + /// Uniform in [-1, 1). + pub fn uniform(&mut self) -> f32 { + ((self.next_u64() >> 40) as f32 / (1u64 << 24) as f32) * 2.0 - 1.0 + } + } + + pub fn sine(freq: f64, amp: f64, sr: f64, n: usize, phase: f64) -> Vec { + (0..n).map(|i| (amp * (2.0 * std::f64::consts::PI * freq * i as f64 / sr + phase).sin()) as f32).collect() + } + + pub fn rms(x: &[f32]) -> f64 { + if x.is_empty() { + return 0.0; + } + (x.iter().map(|&v| (v as f64) * (v as f64)).sum::() / x.len() as f64).sqrt() + } + + pub fn db(x: f64) -> f64 { + 20.0 * x.log10() + } + + /// Amplitude of the `freq` component of `x` (single-bin DFT / Goertzel-style correlation). + pub fn tone_amplitude(x: &[f32], freq: f64, sr: f64) -> f64 { + let (mut re, mut im) = (0.0f64, 0.0f64); + for (i, &v) in x.iter().enumerate() { + let w = 2.0 * std::f64::consts::PI * freq * i as f64 / sr; + re += v as f64 * w.cos(); + im += v as f64 * w.sin(); + } + 2.0 * (re * re + im * im).sqrt() / x.len() as f64 + } +} diff --git a/crates/audio-dsp/src/loudness.rs b/crates/audio-dsp/src/loudness.rs new file mode 100644 index 00000000..1fbe8711 --- /dev/null +++ b/crates/audio-dsp/src/loudness.rs @@ -0,0 +1,693 @@ +//! ITU-R BS.1770-4 / EBU R128 loudness metering. +//! +//! * K-weighting: the two BS.1770 stages (high-shelf "pre-filter" and the RLB high-pass) are +//! derived from their analog prototypes (centre frequency, Q, gain) through the bilinear +//! transform, so any sample rate is supported; at 48 kHz they reproduce the coefficient table in +//! the recommendation. +//! * Momentary (400 ms), short-term (3 s) windows are evaluated every 100 ms (75 % overlap for the +//! 400 ms gating blocks, as specified). +//! * Integrated loudness uses the two-stage gate (−70 LUFS absolute, −10 LU relative). Gating +//! blocks are accumulated in a fixed 0.01 LU histogram that also stores exact energy sums, so the +//! meter can run for arbitrarily long programmes without allocating and the result is exact up +//! to the placement of blocks within 0.01 LU of the relative threshold. +//! * Loudness range (EBU Tech 3342): short-term values, −70 LUFS absolute and −20 LU relative +//! gate, LRA = 95th − 10th percentile. +//! * Sample peak and true peak (polyphase oversampling to ≥ 176.4 kHz, see [`crate::oversample`]). + +use crate::biquad::{Biquad, Coeffs}; +use crate::oversample::{Interpolator, PolyphaseBank, true_peak_factor}; +use std::f64::consts::PI; + +const HIST_MIN: f64 = -70.0; +const HIST_MAX: f64 = 20.0; +const HIST_STEP: f64 = 0.01; +const HIST_BINS: usize = ((HIST_MAX - HIST_MIN) / HIST_STEP) as usize; +/// Short-term window in 100 ms sub-blocks. +const ST_SUBBLOCKS: usize = 30; +/// Momentary window in 100 ms sub-blocks. +const M_SUBBLOCKS: usize = 4; + +/// Loudness of a mean-square energy (already channel-weighted), in LUFS. +#[inline] +pub fn energy_to_lufs(e: f64) -> f64 { + if e <= 0.0 { f64::NEG_INFINITY } else { -0.691 + 10.0 * e.log10() } +} + +#[cfg(test)] +fn lufs_to_energy(l: f64) -> f64 { + 10f64.powf((l + 0.691) / 10.0) +} + +/// Stage 1 of K-weighting: the high-shelf ("head") pre-filter for `sample_rate`. +pub fn k_prefilter(sample_rate: f64) -> Coeffs { + // Analog prototype parameters of the BS.1770 shelving stage. + let f0 = 1_681.974_450_955_533; + let gain_db = 3.999_843_853_973_347; + let q = 0.707_175_236_955_419_6; + let k = (PI * f0 / sample_rate).tan(); + let vh = 10f64.powf(gain_db / 20.0); + let vb = vh.powf(0.499_666_774_154_541_6); + let a0 = 1.0 + k / q + k * k; + Coeffs { + b0: (vh + vb * k / q + k * k) / a0, + b1: 2.0 * (k * k - vh) / a0, + b2: (vh - vb * k / q + k * k) / a0, + a1: 2.0 * (k * k - 1.0) / a0, + a2: (1.0 - k / q + k * k) / a0, + } +} + +/// Stage 2 of K-weighting: the RLB high-pass for `sample_rate`. +pub fn k_rlb(sample_rate: f64) -> Coeffs { + let f0 = 38.135_470_876_024_44; + let q = 0.500_327_037_323_877_3; + let k = (PI * f0 / sample_rate).tan(); + let a0 = 1.0 + k / q + k * k; + Coeffs { b0: 1.0, b1: -2.0, b2: 1.0, a1: 2.0 * (k * k - 1.0) / a0, a2: (1.0 - k / q + k * k) / a0 } +} + +/// Default BS.1770 channel weights for a channel count (L R C LFE Ls Rs [Lb Rb] order for 5.1 / +/// 7.1; LFE is excluded; surrounds get 1.41 ≈ +1.5 dB). +pub fn default_channel_weights(channels: usize) -> Vec { + match channels { + 5 => vec![1.0, 1.0, 1.0, 1.41, 1.41], + 6 => vec![1.0, 1.0, 1.0, 0.0, 1.41, 1.41], + 8 => vec![1.0, 1.0, 1.0, 0.0, 1.41, 1.41, 1.41, 1.41], + n => vec![1.0; n], + } +} + +/// Histogram of block loudness values with exact energy sums per bin. +#[derive(Clone, Debug)] +struct Histogram { + count: Vec, + energy: Vec, +} + +impl Histogram { + fn new() -> Self { + Histogram { count: vec![0; HIST_BINS], energy: vec![0.0; HIST_BINS] } + } + fn clear(&mut self) { + self.count.iter_mut().for_each(|c| *c = 0); + self.energy.iter_mut().for_each(|e| *e = 0.0); + } + fn bin_of(l: f64) -> usize { + (((l - HIST_MIN) / HIST_STEP).floor().max(0.0) as usize).min(HIST_BINS - 1) + } + fn bin_centre(i: usize) -> f64 { + HIST_MIN + (i as f64 + 0.5) * HIST_STEP + } + /// Add a block if it passes the absolute gate. + fn add(&mut self, energy: f64) { + let l = energy_to_lufs(energy); + if l >= HIST_MIN { + let b = Self::bin_of(l); + self.count[b] += 1; + self.energy[b] += energy; + } + } + /// Power mean (energy) of bins from `from` upward; None if empty. + fn mean_energy_from(&self, from: usize) -> Option<(f64, u64)> { + let (mut e, mut n) = (0.0, 0u64); + for i in from..HIST_BINS { + e += self.energy[i]; + n += self.count[i]; + } + (n > 0).then(|| (e / n as f64, n)) + } + /// Relative-gate start bin for a gate of `rel` LU below the absolute-gated mean. + fn relative_gate_bin(&self, rel: f64) -> Option { + let (e, _) = self.mean_energy_from(0)?; + let thr = energy_to_lufs(e) + rel; + Some(if thr < HIST_MIN { 0 } else { Self::bin_of(thr) }) + } +} + +/// Streaming BS.1770-4 / EBU R128 loudness meter for planar or interleaved f32 audio. +#[derive(Clone, Debug)] +pub struct LoudnessMeter { + sample_rate: f64, + channels: usize, + weights: Vec, + pre: Coeffs, + rlb: Coeffs, + filt: Vec<(Biquad, Biquad)>, + /// Per-channel sum of squares in the current 100 ms sub-block. + acc: Vec, + sub_len: usize, + sub_pos: usize, + /// Ring of (weighted energy sum, sample count) for the last 30 sub-blocks. + ring: [(f64, usize); ST_SUBBLOCKS], + ring_pos: usize, + subblocks_seen: u64, + momentary: f64, + short_term: f64, + max_momentary: f64, + max_short_term: f64, + integrated_hist: Histogram, + lra_hist: Histogram, + sample_peak: Vec, + true_peak: Vec, + bank: PolyphaseBank, + interp: Vec, +} + +impl LoudnessMeter { + /// A meter for `channels` channels at `sample_rate` Hz with default BS.1770 weights. + pub fn new(sample_rate: f64, channels: usize) -> Self { + let channels = channels.max(1); + let sub_len = ((sample_rate * 0.1).round() as usize).max(1); + LoudnessMeter { + sample_rate, + channels, + weights: default_channel_weights(channels), + pre: k_prefilter(sample_rate), + rlb: k_rlb(sample_rate), + filt: vec![(Biquad::new(), Biquad::new()); channels], + acc: vec![0.0; channels], + sub_len, + sub_pos: 0, + ring: [(0.0, 0); ST_SUBBLOCKS], + ring_pos: 0, + subblocks_seen: 0, + momentary: f64::NEG_INFINITY, + short_term: f64::NEG_INFINITY, + max_momentary: f64::NEG_INFINITY, + max_short_term: f64::NEG_INFINITY, + integrated_hist: Histogram::new(), + lra_hist: Histogram::new(), + sample_peak: vec![0.0; channels], + true_peak: vec![0.0; channels], + bank: PolyphaseBank::new(true_peak_factor(sample_rate)), + interp: vec![Interpolator::default(); channels], + } + } + + pub fn sample_rate(&self) -> f64 { + self.sample_rate + } + pub fn channels(&self) -> usize { + self.channels + } + + /// Override the per-channel weights (missing entries are 1.0, extra ignored). + pub fn set_channel_weights(&mut self, weights: &[f64]) { + for (i, w) in self.weights.iter_mut().enumerate() { + *w = weights.get(i).copied().unwrap_or(1.0); + } + } + + /// Clear all measurements and filter state. + pub fn reset(&mut self) { + self.integrated_hist.clear(); + self.lra_hist.clear(); + self.filt.iter_mut().for_each(|(a, b)| { + a.reset(); + b.reset() + }); + self.acc.iter_mut().for_each(|a| *a = 0.0); + self.sample_peak.iter_mut().for_each(|p| *p = 0.0); + self.true_peak.iter_mut().for_each(|p| *p = 0.0); + self.interp.iter_mut().for_each(Interpolator::reset); + self.sub_pos = 0; + self.ring = [(0.0, 0); ST_SUBBLOCKS]; + self.ring_pos = 0; + self.subblocks_seen = 0; + self.momentary = f64::NEG_INFINITY; + self.short_term = f64::NEG_INFINITY; + self.max_momentary = f64::NEG_INFINITY; + self.max_short_term = f64::NEG_INFINITY; + } + + #[inline(always)] + fn feed(&mut self, ch: usize, x: f32) { + let (pre, rlb) = &mut self.filt[ch]; + let y = rlb.tick(&self.rlb, pre.tick(&self.pre, x as f64)); + self.acc[ch] += y * y; + let a = x.abs(); + if a > self.sample_peak[ch] { + self.sample_peak[ch] = a; + } + let tp = self.interp[ch].push_peak(&self.bank, x).max(a); + if tp > self.true_peak[ch] { + self.true_peak[ch] = tp; + } + } + + /// Feed a block of planar audio (one slice per channel; extra channels are ignored, missing + /// channels count as silence). + pub fn process(&mut self, channels: &[&[f32]]) { + let n = channels.iter().take(self.channels).map(|c| c.len()).min().unwrap_or(0); + let mut pos = 0; + while pos < n { + let take = (n - pos).min(self.sub_len - self.sub_pos); + for (ch, data) in channels.iter().take(self.channels).enumerate() { + for &x in &data[pos..pos + take] { + self.feed(ch, x); + } + } + for ch in channels.len()..self.channels { + for _ in 0..take { + self.feed(ch, 0.0); + } + } + pos += take; + self.sub_pos += take; + if self.sub_pos == self.sub_len { + self.finish_subblock(); + } + } + } + + /// Feed a block of interleaved audio with `self.channels()` channels per frame. + pub fn process_interleaved(&mut self, data: &[f32]) { + let nch = self.channels; + let frames = data.len() / nch; + let mut f = 0; + while f < frames { + let take = (frames - f).min(self.sub_len - self.sub_pos); + for i in f..f + take { + for ch in 0..nch { + self.feed(ch, data[i * nch + ch]); + } + } + f += take; + self.sub_pos += take; + if self.sub_pos == self.sub_len { + self.finish_subblock(); + } + } + } + + fn finish_subblock(&mut self) { + let mut e = 0.0; + for ch in 0..self.channels { + e += self.weights[ch] * self.acc[ch]; + self.acc[ch] = 0.0; + } + self.ring[self.ring_pos] = (e, self.sub_pos); + self.ring_pos = (self.ring_pos + 1) % ST_SUBBLOCKS; + self.sub_pos = 0; + self.subblocks_seen += 1; + + let window = |k: usize, ring: &[(f64, usize); ST_SUBBLOCKS], pos: usize, sub_len: usize| { + let mut s = 0.0; + for i in 0..k { + s += ring[(pos + ST_SUBBLOCKS - 1 - i) % ST_SUBBLOCKS].0; + } + s / (k * sub_len) as f64 + }; + let em = window(M_SUBBLOCKS, &self.ring, self.ring_pos, self.sub_len); + let es = window(ST_SUBBLOCKS, &self.ring, self.ring_pos, self.sub_len); + self.momentary = energy_to_lufs(em); + self.short_term = energy_to_lufs(es); + if self.subblocks_seen >= M_SUBBLOCKS as u64 { + self.integrated_hist.add(em); + self.max_momentary = self.max_momentary.max(self.momentary); + } + if self.subblocks_seen >= ST_SUBBLOCKS as u64 { + self.lra_hist.add(es); + self.max_short_term = self.max_short_term.max(self.short_term); + } + } + + /// Momentary loudness (last 400 ms), LUFS; updated every 100 ms. −∞ for silence. + pub fn momentary(&self) -> f64 { + self.momentary + } + /// Short-term loudness (last 3 s), LUFS; updated every 100 ms. + pub fn short_term(&self) -> f64 { + self.short_term + } + /// Maximum momentary loudness so far. + pub fn max_momentary(&self) -> f64 { + self.max_momentary + } + /// Maximum short-term loudness so far. + pub fn max_short_term(&self) -> f64 { + self.max_short_term + } + + /// Gated integrated loudness (BS.1770-4), LUFS. −∞ if nothing passed the gates. + pub fn integrated(&self) -> f64 { + let h = &self.integrated_hist; + match h.relative_gate_bin(-10.0).and_then(|b| h.mean_energy_from(b)) { + Some((e, _)) => energy_to_lufs(e), + None => f64::NEG_INFINITY, + } + } + + /// Relative gate threshold used for the integrated loudness (LUFS). + pub fn relative_threshold(&self) -> f64 { + match self.integrated_hist.mean_energy_from(0) { + Some((e, _)) => energy_to_lufs(e) - 10.0, + None => f64::NEG_INFINITY, + } + } + + /// Loudness range (EBU Tech 3342), LU. 0 when too little signal. + pub fn loudness_range(&self) -> f64 { + let h = &self.lra_hist; + let Some(start) = h.relative_gate_bin(-20.0) else { return 0.0 }; + let n: u64 = h.count[start..].iter().sum(); + if n == 0 { + return 0.0; + } + let percentile = |p: f64| { + // Nearest-rank on the sorted gated values. + let rank = ((p * (n - 1) as f64).round() as u64).min(n - 1); + let mut c = 0u64; + for i in start..HIST_BINS { + c += h.count[i]; + if c > rank { + return Histogram::bin_centre(i); + } + } + Histogram::bin_centre(HIST_BINS - 1) + }; + (percentile(0.95) - percentile(0.10)).max(0.0) + } + + /// Maximum absolute sample value over all channels (linear). + pub fn sample_peak(&self) -> f64 { + self.sample_peak.iter().fold(0.0f32, |a, &b| a.max(b)) as f64 + } + /// Per-channel sample peak (linear). + pub fn sample_peak_channel(&self, ch: usize) -> f64 { + self.sample_peak.get(ch).copied().unwrap_or(0.0) as f64 + } + /// Maximum true peak over all channels (linear). + pub fn true_peak(&self) -> f64 { + self.true_peak.iter().fold(0.0f32, |a, &b| a.max(b)) as f64 + } + /// Per-channel true peak (linear). + pub fn true_peak_channel(&self, ch: usize) -> f64 { + self.true_peak.get(ch).copied().unwrap_or(0.0) as f64 + } + /// Maximum true peak in dBTP. + pub fn true_peak_dbtp(&self) -> f64 { + 20.0 * self.true_peak().log10() + } + + /// Snapshot of all programme measurements. + pub fn summary(&self) -> LoudnessSummary { + LoudnessSummary { + integrated_lufs: self.integrated(), + loudness_range_lu: self.loudness_range(), + max_momentary_lufs: self.max_momentary, + max_short_term_lufs: self.max_short_term, + sample_peak_dbfs: 20.0 * self.sample_peak().log10(), + true_peak_dbtp: self.true_peak_dbtp(), + } + } +} + +/// Programme loudness measurements (EBU R128 set). +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct LoudnessSummary { + pub integrated_lufs: f64, + pub loudness_range_lu: f64, + pub max_momentary_lufs: f64, + pub max_short_term_lufs: f64, + pub sample_peak_dbfs: f64, + pub true_peak_dbtp: f64, +} + +/// Measure a whole planar programme in one call. +pub fn measure(channels: &[&[f32]], sample_rate: f64) -> LoudnessSummary { + let mut m = LoudnessMeter::new(sample_rate, channels.len()); + m.process(channels); + m.summary() +} + +/// Gain (dB) that moves a programme measured at `measured_lufs` to `target_lufs`. +/// Returns 0 for silent / unmeasurable input. +pub fn normalize_gain_db(measured_lufs: f64, target_lufs: f64) -> f64 { + if measured_lufs.is_finite() && target_lufs.is_finite() { target_lufs - measured_lufs } else { 0.0 } +} + +/// Like [`normalize_gain_db`] but limited so that the resulting true peak does not exceed +/// `max_true_peak_dbtp` (a linear gain can't fix peaks, so the loudness target is sacrificed — +/// the behaviour of "Normalize" / Essential Sound auto-match without a limiter). +pub fn normalize_gain_db_peak_limited(measured_lufs: f64, target_lufs: f64, true_peak_dbtp: f64, max_true_peak_dbtp: f64) -> f64 { + let g = normalize_gain_db(measured_lufs, target_lufs); + if true_peak_dbtp.is_finite() { g.min(max_true_peak_dbtp - true_peak_dbtp) } else { g } +} + +/// Expected loudness energy weighting of a pure tone at `freq` through K-weighting (|H|², for +/// tests and analytical tooling). +pub fn k_weighting_power(freq: f64, sample_rate: f64) -> f64 { + let a = k_prefilter(sample_rate).magnitude(freq, sample_rate); + let b = k_rlb(sample_rate).magnitude(freq, sample_rate); + (a * b).powi(2) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::testutil::*; + + fn dbfs(db: f64) -> f64 { + 10f64.powf(db / 20.0) + } + + /// Stereo programme made of (seconds, dBFS) sine segments. + fn segments(freq: f64, sr: f64, segs: &[(f64, f64)]) -> Vec { + let mut out = Vec::new(); + let mut phase_idx = 0usize; + for &(secs, level) in segs { + let n = (secs * sr).round() as usize; + let a = dbfs(level); + for _ in 0..n { + out.push((a * (2.0 * PI * freq * phase_idx as f64 / sr).sin()) as f32); + phase_idx += 1; + } + } + out + } + + fn stereo_measure(x: &[f32], sr: f64) -> LoudnessMeter { + let mut m = LoudnessMeter::new(sr, 2); + // Feed in odd-sized blocks to exercise the sub-block splitting. + for chunk in x.chunks(1237) { + m.process(&[chunk, chunk]); + } + m + } + + #[test] + fn k_weighting_matches_bs1770_table_at_48k() { + let p = k_prefilter(48000.0); + assert!((p.b0 - 1.535_124_859_586_97).abs() < 1e-8); + assert!((p.b1 + 2.691_696_189_406_38).abs() < 1e-8); + assert!((p.b2 - 1.198_392_810_852_85).abs() < 1e-8); + assert!((p.a1 + 1.690_659_293_182_41).abs() < 1e-8); + assert!((p.a2 - 0.732_480_774_215_85).abs() < 1e-8); + let r = k_rlb(48000.0); + assert!((r.a1 + 1.990_047_454_833_98).abs() < 1e-8); + assert!((r.a2 - 0.990_072_250_366_21).abs() < 1e-8); + } + + #[test] + fn k_weighting_shape_is_sample_rate_independent() { + for sr in [32000.0, 44100.0, 48000.0, 96000.0, 192000.0] { + let db = |f| 10.0 * k_weighting_power(f, sr).log10(); + // ~+0.69 dB at 1 kHz, ~+4 dB high shelf, strong low-cut. + assert!((db(1000.0) - 0.691).abs() < 0.03, "sr {sr}: {}", db(1000.0)); + assert!((db(10000.0) - 4.0).abs() < 0.3, "sr {sr}"); + assert!(db(20.0) < -10.0, "sr {sr}"); + } + } + + #[test] + fn tech3341_case1_and_2_stereo_sine() { + for (sr, freq) in [(48000.0, 1000.0), (44100.0, 997.0), (96000.0, 1000.0)] { + for level in [-23.0, -33.0] { + let x = segments(freq, sr, &[(20.0, level)]); + let m = stereo_measure(&x, sr); + assert!((m.integrated() - level).abs() < 0.1, "sr {sr} I {}", m.integrated()); + assert!((m.momentary() - level).abs() < 0.1, "M {}", m.momentary()); + assert!((m.short_term() - level).abs() < 0.1, "S {}", m.short_term()); + assert!((m.max_momentary() - level).abs() < 0.1); + assert!(m.loudness_range() < 0.1); + } + } + } + + #[test] + fn interleaved_matches_planar() { + let sr = 48000.0; + let x = segments(1000.0, sr, &[(5.0, -20.0)]); + let il: Vec = x.iter().flat_map(|&v| [v, v]).collect(); + let mut m = LoudnessMeter::new(sr, 2); + m.process_interleaved(&il); + let p = stereo_measure(&x, sr); + assert!((m.integrated() - p.integrated()).abs() < 1e-9); + } + + #[test] + fn tech3341_case3_4_gating() { + let sr = 48000.0; + // Case 3: -36 / -23 / -36 dBFS for 10 / 60 / 10 s → -23.0 LUFS. + let x = segments(1000.0, sr, &[(10.0, -36.0), (60.0, -23.0), (10.0, -36.0)]); + let m = stereo_measure(&x, sr); + assert!((m.integrated() + 23.0).abs() < 0.1, "{}", m.integrated()); + // Case 4: -72 / -36 / -23 / -36 / -72 dBFS for 10/10/60/10/10 s → -23.0 LUFS. + let x = segments(1000.0, sr, &[(10.0, -72.0), (10.0, -36.0), (60.0, -23.0), (10.0, -36.0), (10.0, -72.0)]); + let m = stereo_measure(&x, sr); + assert!((m.integrated() + 23.0).abs() < 0.1, "{}", m.integrated()); + } + + #[test] + fn tech3341_case5_relative_gate() { + let sr = 48000.0; + let x = segments(1000.0, sr, &[(20.0, -26.0), (20.1, -20.0), (20.0, -26.0)]); + let m = stereo_measure(&x, sr); + assert!((m.integrated() + 23.0).abs() < 0.1, "{}", m.integrated()); + } + + #[test] + fn silence_and_below_gate() { + let sr = 48000.0; + let mut m = LoudnessMeter::new(sr, 2); + let z = vec![0.0f32; 48000 * 2]; + m.process(&[&z, &z]); + assert_eq!(m.integrated(), f64::NEG_INFINITY); + assert_eq!(m.loudness_range(), 0.0); + let x = segments(1000.0, sr, &[(5.0, -80.0)]); + let m = stereo_measure(&x, sr); + assert_eq!(m.integrated(), f64::NEG_INFINITY); + } + + #[test] + fn tech3342_loudness_range() { + let sr = 48000.0; + let cases: [(&[(f64, f64)], f64); 4] = [ + (&[(20.0, -20.0), (20.0, -30.0)], 10.0), + (&[(20.0, -20.0), (20.0, -15.0)], 5.0), + (&[(20.0, -40.0), (20.0, -20.0)], 20.0), + (&[(20.0, -50.0), (20.0, -35.0), (20.0, -20.0), (20.0, -35.0), (20.0, -50.0)], 15.0), + ]; + for (segs, expect) in cases { + let x = segments(1000.0, sr, segs); + let m = stereo_measure(&x, sr); + assert!((m.loudness_range() - expect).abs() < 1.0, "LRA {} expected {expect}", m.loudness_range()); + } + } + + #[test] + fn log_sweep_matches_analytic_k_weighting() { + let sr = 48000.0; + let secs = 20.0; + let n = (sr * secs) as usize; + let (f1, f2) = (20.0f64, 20000.0f64); + let k = (f2 / f1).ln(); + let amp = dbfs(-20.0); + let mut phase = 0.0f64; + let mut x = Vec::with_capacity(n); + let mut expected_e = 0.0; + for i in 0..n { + let f = f1 * (k * i as f64 / n as f64).exp(); + x.push((amp * phase.sin()) as f32); + phase += 2.0 * PI * f / sr; + // Stereo: 2 channels × A²/2 × |H(f)|². + expected_e += amp * amp * k_weighting_power(f, sr); + } + expected_e /= n as f64; + // Ungated average of the same signal: compare with the mean of momentary energies, i.e. + // the integrated value is gated, so build the expectation the same way using short blocks. + let m = stereo_measure(&x, sr); + let expect = energy_to_lufs(expected_e); + // Gating removes the quiet low-frequency start, so integrated ≥ ungated expectation; the + // difference must be small for a sweep that is mostly well above the gates. + let i = m.integrated(); + assert!(i >= expect - 0.1 && i - expect < 0.6, "I {i} vs ungated analytic {expect}"); + } + + #[test] + fn surround_weights() { + let sr = 48000.0; + let x = segments(1000.0, sr, &[(10.0, -30.0)]); + let z = vec![0.0f32; x.len()]; + // Signal only in Ls of a 5.1 layout → +1.5 dB weighting vs. front channel. + let mut front = LoudnessMeter::new(sr, 6); + front.process(&[&x, &z, &z, &z, &z, &z]); + let mut surr = LoudnessMeter::new(sr, 6); + surr.process(&[&z, &z, &z, &z, &x, &z]); + let mut lfe = LoudnessMeter::new(sr, 6); + lfe.process(&[&z, &z, &z, &x, &z, &z]); + assert!((surr.integrated() - front.integrated() - 10.0 * 1.41f64.log10()).abs() < 0.01); + assert_eq!(lfe.integrated(), f64::NEG_INFINITY); + } + + #[test] + fn true_peak_of_intersample_peak_signal() { + for sr in [48000.0, 44100.0] { + // fs/4 sine with 45° phase: every sample is at ±0.5·√2·A, true peak A (0 dBFS). + let n = sr as usize; + let x: Vec = (0..n).map(|i| (PI / 2.0 * i as f64 + PI / 4.0).sin() as f32).collect(); + let m = stereo_measure(&x, sr); + assert!((20.0 * m.sample_peak().log10() + 3.01).abs() < 0.01); + // EBU Tech 3341 true-peak tolerance: +0.2 / −0.4 dB. + let tp = m.true_peak_dbtp(); + assert!((-0.4..=0.2).contains(&tp), "sr {sr} TP {tp}"); + } + // A tone at an off-grid frequency with a 1 kHz-ish content: TP ≈ amplitude. + let x = sine(997.0, dbfs(-6.0), 48000.0, 48000, 0.3); + let m = stereo_measure(&x, 48000.0); + assert!((m.true_peak_dbtp() + 6.0).abs() < 0.2); + // Tech 3341 case 15-ish: 0 dBFS 1 kHz-ish tone at 48 kHz sampled off-peak. + let x = sine(11025.0, 1.0, 44100.0, 44100, 0.0); + let m = stereo_measure(&x, 44100.0); + assert!((-0.4..=0.2).contains(&m.true_peak_dbtp()), "{}", m.true_peak_dbtp()); + } + + #[test] + fn reset_clears_everything() { + let sr = 48000.0; + let x = segments(1000.0, sr, &[(5.0, -20.0)]); + let mut m = stereo_measure(&x, sr); + m.reset(); + assert_eq!(m.integrated(), f64::NEG_INFINITY); + assert_eq!(m.true_peak(), 0.0); + let x = segments(1000.0, sr, &[(5.0, -30.0)]); + m.process(&[&x, &x]); + assert!((m.integrated() + 30.0).abs() < 0.1); + } + + #[test] + fn normalize_helpers() { + assert_eq!(normalize_gain_db(-30.0, -23.0), 7.0); + assert_eq!(normalize_gain_db(f64::NEG_INFINITY, -23.0), 0.0); + assert_eq!(normalize_gain_db_peak_limited(-30.0, -23.0, -3.0, -1.0), 2.0); + assert_eq!(normalize_gain_db_peak_limited(-30.0, -23.0, -10.0, -1.0), 7.0); + // Applying the gain lands on target. + let sr = 48000.0; + let x = segments(1000.0, sr, &[(10.0, -31.3)]); + let m = stereo_measure(&x, sr); + let g = 10f64.powf(normalize_gain_db(m.integrated(), -16.0) / 20.0) as f32; + let y: Vec = x.iter().map(|v| v * g).collect(); + let m2 = stereo_measure(&y, sr); + assert!((m2.integrated() + 16.0).abs() < 0.02); + assert!((lufs_to_energy(energy_to_lufs(0.3)) - 0.3).abs() < 1e-12); + } + + #[test] + fn tiny_blocks_equal_big_blocks() { + let sr = 48000.0; + let x = segments(440.0, sr, &[(4.0, -18.0)]); + let mut a = LoudnessMeter::new(sr, 1); + a.process(&[&x]); + let mut b = LoudnessMeter::new(sr, 1); + b.process(&[&[]]); + for s in x.chunks(1) { + b.process(&[s]); + } + assert_eq!(a.integrated(), b.integrated()); + assert_eq!(a.true_peak(), b.true_peak()); + // Mono is summed without the dual-mono +3 dB (BS.1770). + let expect = -0.691 + 20.0 * rms(&x).log10() + 10.0 * k_weighting_power(440.0, sr).log10(); + assert!((a.integrated() - expect).abs() < 0.1, "{} vs {expect}", a.integrated()); + } +} diff --git a/crates/audio-dsp/src/oversample.rs b/crates/audio-dsp/src/oversample.rs new file mode 100644 index 00000000..6bdfa4c3 --- /dev/null +++ b/crates/audio-dsp/src/oversample.rs @@ -0,0 +1,144 @@ +//! Polyphase FIR interpolator used for true-peak estimation (BS.1770-4 Annex 2 style). +//! +//! The prototype is a Kaiser-windowed sinc low-pass (cut-off at the original Nyquist) split into +//! `factor` phases of [`TAPS_PER_PHASE`] taps; each phase is normalised to unity DC gain. + +use std::f64::consts::PI; + +/// Taps per polyphase branch. +pub const TAPS_PER_PHASE: usize = 16; + +/// Oversampling factor giving an effective rate of at least 176.4 kHz (4× at 44.1/48 kHz, +/// 2× at 88.2/96 kHz, 1× above), as recommended by BS.1770-4 for true-peak metering. +pub fn true_peak_factor(sample_rate: f64) -> usize { + if sample_rate < 88_200.0 { + 4 + } else if sample_rate < 176_400.0 { + 2 + } else { + 1 + } +} + +fn bessel_i0(x: f64) -> f64 { + let mut sum = 1.0; + let mut term = 1.0; + let q = x * x / 4.0; + for k in 1..50 { + term *= q / (k * k) as f64; + sum += term; + if term < sum * 1e-17 { + break; + } + } + sum +} + +/// Polyphase coefficient bank shared by all channels. +#[derive(Clone, Debug)] +pub struct PolyphaseBank { + factor: usize, + /// `coeffs[p * TAPS_PER_PHASE + j]` multiplies `x[n - j]` for output phase `p`. + coeffs: Vec, +} + +impl PolyphaseBank { + pub fn new(factor: usize) -> Self { + let factor = factor.max(1); + let len = factor * TAPS_PER_PHASE; + let centre = (len - 1) as f64 / 2.0; + let beta = 8.0; + let i0b = bessel_i0(beta); + let proto: Vec = (0..len) + .map(|m| { + let t = (m as f64 - centre) / factor as f64; + let sinc = if t.abs() < 1e-12 { 1.0 } else { (PI * t).sin() / (PI * t) }; + let r = (m as f64 - centre) / (centre + 1.0); + let w = bessel_i0(beta * (1.0 - r * r).max(0.0).sqrt()) / i0b; + sinc * w + }) + .collect(); + let mut coeffs = vec![0.0f32; len]; + for p in 0..factor { + let sum: f64 = (0..TAPS_PER_PHASE).map(|j| proto[p + factor * j]).sum(); + for j in 0..TAPS_PER_PHASE { + coeffs[p * TAPS_PER_PHASE + j] = (proto[p + factor * j] / sum) as f32; + } + } + PolyphaseBank { factor, coeffs } + } + pub fn factor(&self) -> usize { + self.factor + } + /// Approximate delay (in input samples) between an input sample and the interpolated + /// values around it. + pub fn delay(&self) -> usize { + TAPS_PER_PHASE / 2 + } +} + +/// Per-channel interpolator history. +#[derive(Clone, Debug)] +pub struct Interpolator { + /// Double-written ring so a contiguous window of TAPS_PER_PHASE is always available. + hist: [f32; 2 * TAPS_PER_PHASE], + pos: usize, +} + +impl Default for Interpolator { + fn default() -> Self { + Interpolator { hist: [0.0; 2 * TAPS_PER_PHASE], pos: 0 } + } +} + +impl Interpolator { + pub fn reset(&mut self) { + *self = Self::default(); + } + + /// Push one input sample; returns the maximum absolute value of the `factor` interpolated + /// output samples produced for it. + #[inline] + pub fn push_peak(&mut self, bank: &PolyphaseBank, x: f32) -> f32 { + if bank.factor == 1 { + return x.abs(); + } + // Newest sample at window[TAPS-1-0]... store so that window[j] = x[n - j]. + self.pos = if self.pos == 0 { TAPS_PER_PHASE - 1 } else { self.pos - 1 }; + self.hist[self.pos] = x; + self.hist[self.pos + TAPS_PER_PHASE] = x; + let win = &self.hist[self.pos..self.pos + TAPS_PER_PHASE]; + let mut peak = 0.0f32; + for p in 0..bank.factor { + let c = &bank.coeffs[p * TAPS_PER_PHASE..(p + 1) * TAPS_PER_PHASE]; + let y: f32 = c.iter().zip(win).map(|(a, b)| a * b).sum(); + peak = peak.max(y.abs()); + } + peak + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn phases_have_unity_dc_and_interpolate_sines() { + let bank = PolyphaseBank::new(4); + for p in 0..4 { + let s: f32 = bank.coeffs[p * TAPS_PER_PHASE..(p + 1) * TAPS_PER_PHASE].iter().sum(); + assert!((s - 1.0).abs() < 1e-5); + } + // Sine at fs/4 with 45° phase: samples at ±0.707, true peak 1.0. + let mut ip = Interpolator::default(); + let mut peak = 0.0f32; + for n in 0..400 { + let x = (std::f64::consts::FRAC_PI_2 * n as f64 + std::f64::consts::FRAC_PI_4).sin() as f32; + let p = ip.push_peak(&bank, x); + if n > 40 { + peak = peak.max(p); + } + } + assert!((peak - 1.0).abs() < 0.02, "peak {peak}"); + } +} diff --git a/crates/audio-dsp/src/smooth.rs b/crates/audio-dsp/src/smooth.rs new file mode 100644 index 00000000..396cdfba --- /dev/null +++ b/crates/audio-dsp/src/smooth.rs @@ -0,0 +1,105 @@ +//! Per-sample parameter smoothing (linear ramps) to avoid zipper noise. + +/// A value that ramps linearly to its target over a fixed number of samples. +#[derive(Clone, Copy, Debug)] +pub struct Smoothed { + current: f64, + target: f64, + step: f64, + remaining: u32, + ramp: u32, +} + +impl Smoothed { + /// `ramp_samples` = length of a full transition (0 = no smoothing). + pub fn new(value: f32, ramp_samples: u32) -> Self { + Smoothed { current: value as f64, target: value as f64, step: 0.0, remaining: 0, ramp: ramp_samples } + } + /// Convenience: ramp length from a time in milliseconds. + pub fn with_ms(value: f32, sample_rate: f32, ms: f32) -> Self { + Self::new(value, (sample_rate * ms / 1000.0).round().max(0.0) as u32) + } + /// Start ramping towards `target`. + pub fn set(&mut self, target: f32) { + let t = target as f64; + if t == self.target { + return; + } + self.target = t; + if self.ramp == 0 { + self.current = t; + self.remaining = 0; + } else { + self.remaining = self.ramp; + self.step = (t - self.current) / self.ramp as f64; + } + } + /// Jump straight to the target. + pub fn snap(&mut self) { + self.current = self.target; + self.remaining = 0; + } + /// Set and jump. + pub fn set_immediate(&mut self, v: f32) { + self.target = v as f64; + self.snap(); + } + /// Advance one sample and return the new value. + #[inline] + pub fn tick(&mut self) -> f32 { + if self.remaining > 0 { + self.remaining -= 1; + if self.remaining == 0 { + self.current = self.target; + } else { + self.current += self.step; + } + } + self.current as f32 + } + /// Current value without advancing. + #[inline] + pub fn value(&self) -> f32 { + self.current as f32 + } + pub fn target(&self) -> f32 { + self.target as f32 + } + pub fn is_smoothing(&self) -> bool { + self.remaining > 0 + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn ramps_linearly_and_lands_exactly() { + let mut s = Smoothed::new(0.0, 4); + s.set(1.0); + let v: Vec = (0..6).map(|_| s.tick()).collect(); + assert_eq!(v, vec![0.25, 0.5, 0.75, 1.0, 1.0, 1.0]); + assert!(!s.is_smoothing()); + } + + #[test] + fn retarget_mid_ramp_is_continuous() { + let mut s = Smoothed::new(0.0, 10); + s.set(1.0); + for _ in 0..5 { + s.tick(); + } + let mid = s.value(); + s.set(0.0); + let n = s.tick(); + assert!((n - mid).abs() <= 0.1 + 1e-6); + } + + #[test] + fn zero_ramp_is_immediate() { + let mut s = Smoothed::new(0.0, 0); + s.set(3.0); + assert_eq!(s.value(), 3.0); + } +}