nvrl
d77b88c7a9
KWeighting is one channel of the standard's pre-filter and RLB high-pass, computed from the analog prototypes so every sample rate is exact. lufs() turns a mean square summed over channels into LUFS. A full-scale 997 Hz sine reads -3.01 LUFS within 0.03 LU from 44.1 to 192 kHz.
269 lines
7.6 KiB
Rust
269 lines
7.6 KiB
Rust
//! Shared allocation-free DSP building blocks: filters, pitch detection, notes, WAV files.
|
||
pub mod loudness;
|
||
pub mod meter;
|
||
pub mod oversample;
|
||
pub mod svf;
|
||
pub mod tracker;
|
||
pub mod wav;
|
||
|
||
use std::f32::consts::TAU;
|
||
|
||
/// Transposed direct form II biquad (RBJ cookbook coefficients).
|
||
#[derive(Clone, Copy, Debug, Default)]
|
||
pub struct Biquad {
|
||
b0: f32,
|
||
b1: f32,
|
||
b2: f32,
|
||
a1: f32,
|
||
a2: f32,
|
||
z1: f32,
|
||
z2: f32,
|
||
}
|
||
impl Biquad {
|
||
fn from(b: [f32; 3], a: [f32; 3]) -> Self {
|
||
Self {
|
||
b0: b[0] / a[0],
|
||
b1: b[1] / a[0],
|
||
b2: b[2] / a[0],
|
||
a1: a[1] / a[0],
|
||
a2: a[2] / a[0],
|
||
z1: 0.,
|
||
z2: 0.,
|
||
}
|
||
}
|
||
pub fn lowpass(rate: f32, hz: f32, q: f32) -> Self {
|
||
let w = TAU * (hz / rate).clamp(1e-5, 0.49);
|
||
let (s, c) = w.sin_cos();
|
||
let alpha = s / (2. * q);
|
||
Self::from(
|
||
[(1. - c) / 2., 1. - c, (1. - c) / 2.],
|
||
[1. + alpha, -2. * c, 1. - alpha],
|
||
)
|
||
}
|
||
pub fn highpass(rate: f32, hz: f32, q: f32) -> Self {
|
||
let w = TAU * (hz / rate).clamp(1e-5, 0.49);
|
||
let (s, c) = w.sin_cos();
|
||
let alpha = s / (2. * q);
|
||
Self::from(
|
||
[(1. + c) / 2., -(1. + c), (1. + c) / 2.],
|
||
[1. + alpha, -2. * c, 1. - alpha],
|
||
)
|
||
}
|
||
/// Peaking EQ; used by tests to model resonant material.
|
||
pub fn peak(rate: f32, hz: f32, q: f32, gain_db: f32) -> Self {
|
||
let w = TAU * (hz / rate).clamp(1e-5, 0.49);
|
||
let (s, c) = w.sin_cos();
|
||
let a = 10f32.powf(gain_db / 40.);
|
||
let alpha = s / (2. * q);
|
||
Self::from(
|
||
[1. + alpha * a, -2. * c, 1. - alpha * a],
|
||
[1. + alpha / a, -2. * c, 1. - alpha / a],
|
||
)
|
||
}
|
||
/// Replace the coefficients and keep the filter state, for smooth parameter changes.
|
||
pub fn retune(&mut self, other: Self) {
|
||
let (z1, z2) = (self.z1, self.z2);
|
||
*self = other;
|
||
self.z1 = z1;
|
||
self.z2 = z2;
|
||
}
|
||
#[inline]
|
||
pub fn process(&mut self, x: f32) -> f32 {
|
||
let y = self.b0 * x + self.z1;
|
||
self.z1 = self.b1 * x - self.a1 * y + self.z2;
|
||
self.z2 = self.b2 * x - self.a2 * y;
|
||
y
|
||
}
|
||
pub fn reset(&mut self) {
|
||
self.z1 = 0.;
|
||
self.z2 = 0.;
|
||
}
|
||
}
|
||
|
||
/// Fourth-order Linkwitz–Riley crossover. Low plus high is an allpass: flat magnitude.
|
||
#[derive(Clone, Copy, Debug, Default)]
|
||
pub struct Crossover {
|
||
low: [Biquad; 2],
|
||
high: [Biquad; 2],
|
||
}
|
||
impl Crossover {
|
||
pub fn new(rate: f32, hz: f32) -> Self {
|
||
let q = std::f32::consts::FRAC_1_SQRT_2;
|
||
Self {
|
||
low: [Biquad::lowpass(rate, hz, q); 2],
|
||
high: [Biquad::highpass(rate, hz, q); 2],
|
||
}
|
||
}
|
||
pub fn retune(&mut self, rate: f32, hz: f32) {
|
||
let fresh = Self::new(rate, hz);
|
||
for i in 0..2 {
|
||
self.low[i].retune(fresh.low[i]);
|
||
self.high[i].retune(fresh.high[i]);
|
||
}
|
||
}
|
||
/// Returns (low, high).
|
||
#[inline]
|
||
pub fn process(&mut self, x: f32) -> (f32, f32) {
|
||
let low = self.low[0].process(x);
|
||
let low = self.low[1].process(low);
|
||
let high = self.high[0].process(x);
|
||
let high = self.high[1].process(high);
|
||
// LR4 bands are in phase at the crossover point, so they sum to an allpass.
|
||
(low, high)
|
||
}
|
||
pub fn reset(&mut self) {
|
||
for f in self.low.iter_mut().chain(self.high.iter_mut()) {
|
||
f.reset();
|
||
}
|
||
}
|
||
}
|
||
|
||
/// YIN fundamental-frequency estimator with preallocated scratch.
|
||
pub struct Yin {
|
||
diff: Vec<f32>,
|
||
pub threshold: f32,
|
||
}
|
||
impl Yin {
|
||
/// `max_lag` is the longest period in samples that will ever be requested.
|
||
pub fn new(max_lag: usize) -> Self {
|
||
Self {
|
||
diff: vec![0.; max_lag + 2],
|
||
threshold: 0.15,
|
||
}
|
||
}
|
||
pub fn max_lag(&self) -> usize {
|
||
self.diff.len() - 2
|
||
}
|
||
/// Estimates the fundamental of `x` in Hz with its clarity (1 = perfectly periodic).
|
||
/// `x` must hold at least twice the longest period. Never allocates.
|
||
pub fn estimate(
|
||
&mut self,
|
||
x: &[f32],
|
||
rate: f32,
|
||
min_hz: f32,
|
||
max_hz: f32,
|
||
) -> Option<(f32, f32)> {
|
||
let hi = ((rate / min_hz).ceil() as usize).min(self.max_lag());
|
||
let lo = ((rate / max_hz).floor() as usize).max(2);
|
||
if hi <= lo + 2 || x.len() < hi + 2 {
|
||
return None;
|
||
}
|
||
let window = x.len() - hi - 1;
|
||
if window < hi {
|
||
return None;
|
||
}
|
||
let energy: f32 = x[..window].iter().map(|v| v * v).sum();
|
||
if energy < 1e-9 * window as f32 {
|
||
return None;
|
||
}
|
||
let d = &mut self.diff;
|
||
d[0] = 1.;
|
||
let mut running = 0.;
|
||
for tau in 1..=hi + 1 {
|
||
let mut s = 0.;
|
||
for j in 0..window {
|
||
let e = x[j] - x[j + tau];
|
||
s += e * e;
|
||
}
|
||
running += s;
|
||
// Cumulative-mean normalisation: suppresses the zero-lag dip.
|
||
d[tau] = if running > 0. {
|
||
s * tau as f32 / running
|
||
} else {
|
||
1.
|
||
};
|
||
}
|
||
let mut best = None;
|
||
let mut tau = lo;
|
||
while tau < hi {
|
||
if d[tau] < self.threshold {
|
||
while tau + 1 < hi && d[tau + 1] < d[tau] {
|
||
tau += 1;
|
||
}
|
||
best = Some(tau);
|
||
break;
|
||
}
|
||
tau += 1;
|
||
}
|
||
let tau = match best {
|
||
Some(t) => t,
|
||
None => {
|
||
let (t, v) = (lo..hi)
|
||
.map(|t| (t, d[t]))
|
||
.min_by(|a, b| a.1.total_cmp(&b.1))?;
|
||
if v > 0.35 {
|
||
return None;
|
||
}
|
||
t
|
||
}
|
||
};
|
||
let (a, b, c) = (d[tau - 1], d[tau], d[tau + 1]);
|
||
let denom = a - 2. * b + c;
|
||
let shift = if denom.abs() > 1e-12 {
|
||
(0.5 * (a - c) / denom).clamp(-1., 1.)
|
||
} else {
|
||
0.
|
||
};
|
||
Some((rate / (tau as f32 + shift), (1. - b).clamp(0., 1.)))
|
||
}
|
||
}
|
||
|
||
pub fn midi_to_hz(note: f32) -> f32 {
|
||
440. * 2f32.powf((note - 69.) / 12.)
|
||
}
|
||
pub fn hz_to_midi(hz: f32) -> f32 {
|
||
69. + 12. * (hz.max(1e-3) / 440.).log2()
|
||
}
|
||
/// Scientific note name: MIDI 60 = C4, 43.65 Hz = F1. FL Studio labels MIDI 60 as C5.
|
||
pub fn note_name(note: u8) -> String {
|
||
const NAMES: [&str; 12] = [
|
||
"C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B",
|
||
];
|
||
format!("{}{}", NAMES[note as usize % 12], note as i32 / 12 - 1)
|
||
}
|
||
pub fn db_to_gain(db: f32) -> f32 {
|
||
10f32.powf(db / 20.)
|
||
}
|
||
pub fn gain_to_db(g: f32) -> f32 {
|
||
20. * g.abs().max(1e-9).log10()
|
||
}
|
||
/// One-pole smoothing coefficient for a time constant in milliseconds.
|
||
pub fn coefficient(rate: f32, ms: f32) -> f32 {
|
||
if ms <= 0. {
|
||
1.
|
||
} else {
|
||
1. - (-1. / (rate * ms * 0.001)).exp()
|
||
}
|
||
}
|
||
/// PolyBLEP sawtooth oscillator: cheap and close to band-limited, for tests and previews.
|
||
pub struct Saw {
|
||
phase: f32,
|
||
inc: f32,
|
||
}
|
||
impl Saw {
|
||
pub fn new(rate: f32, hz: f32) -> Self {
|
||
Self {
|
||
phase: 0.,
|
||
inc: hz / rate,
|
||
}
|
||
}
|
||
pub fn sample(&mut self) -> f32 {
|
||
let t = self.phase;
|
||
let dt = self.inc;
|
||
let blep = if t < dt {
|
||
let x = t / dt;
|
||
x + x - x * x - 1.
|
||
} else if t > 1. - dt {
|
||
let x = (t - 1.) / dt;
|
||
x * x + x + x + 1.
|
||
} else {
|
||
0.
|
||
};
|
||
self.phase += dt;
|
||
if self.phase >= 1. {
|
||
self.phase -= 1.;
|
||
}
|
||
2. * t - 1. - blep
|
||
}
|
||
}
|