sound_to_patch/estimate/envelope/
track.rs

1//! Level tracks on a 1 ms hop: the centred RMS of the full band and the
2//! zero-phase magnitude of single harmonics.
3
4use crate::audio::SampleRate;
5use crate::cast;
6use crate::dsp::{FilterError, band_magnitude};
7use crate::units::{Hz, Seconds};
8
9/// The hop of every envelope track.
10pub const HOP: Seconds = Seconds::from_const(0.001);
11/// The RMS window of a note spans this many periods of its fundamental.
12///
13/// A whole number of periods of a periodic wave leaves no ripple, and one
14/// period follows an attack best. On generated saws a window of two
15/// periods read the 10-90 % rise of a 20 ms attack at 55 Hz as 34 ms and
16/// one period as 19.4 ms; at 220 Hz both read 18 ms for 16 ms.
17pub const RMS_PERIODS: f64 = 1.0;
18/// The RMS window of the pump spans this many periods of the lowest note.
19///
20/// A chord is not periodic in one period of its lowest note, and one
21/// period lets its beating through: on the pump's saw chord the floor rose
22/// from 1.8 to 3.7 dB.
23pub const PUMP_RMS_PERIODS: f64 = 2.0;
24/// The RMS window spans at least this long, in whole periods.
25///
26/// A window of whole periods leaves no ripple, and a short one follows an
27/// attack best. On generated saws a floor of 5 ms read the 10-90 % rise of
28/// a 10 ms attack at 440 Hz as 9.8 ms for 8 ms, and 2 ms read 7.6 ms; at
29/// 220 Hz, where 5 ms is 1.1 periods, a 20 ms attack read 17.7 and 14.9
30/// ms for 16 ms.
31pub const RMS_WINDOW_MIN: Seconds = Seconds::from_const(0.002);
32/// The low-pass of a harmonic track sits at most this far from the
33/// harmonic. A wider band follows an attack faster.
34pub const HARMONIC_CUTOFF_MAX: Hz = Hz::from_const(40.0);
35/// The low-pass of a harmonic track sits at least this far from the
36/// harmonic. A 12 Hz low-pass smears a rise by about 30 ms.
37pub const HARMONIC_CUTOFF_MIN: Hz = Hz::from_const(12.0);
38
39/// The level floor of a track, in dB: digital silence reads this.
40const FLOOR_DB: f64 = -200.0;
41
42/// A level track in dB, one value per [`HOP`] from `start`.
43#[derive(Clone, Debug, PartialEq)]
44pub struct LevelTrack {
45    start: Seconds,
46    db: Vec<f64>,
47    window: Seconds,
48}
49
50impl LevelTrack {
51    /// Pairs dB values with the time of the first.
52    #[must_use]
53    pub fn new(start: Seconds, db: Vec<f64>) -> Self {
54        Self {
55            start,
56            db,
57            window: Seconds::clamped(0.0),
58        }
59    }
60
61    /// The centred RMS of `samples` in dB, from `from` to `to`, over the
62    /// window of [`rms_window`]. `samples` starts at sample `first` of the
63    /// stem, so a note reads only its own span. The window shrinks at the
64    /// ends of the samples.
65    #[must_use]
66    pub fn rms(
67        samples: &[f64],
68        first: usize,
69        rate: SampleRate,
70        (from, to): (Seconds, Seconds),
71        f0: Hz,
72        periods: f64,
73    ) -> Self {
74        let fs = rate.hz().get();
75        let window = rms_window(f0, periods);
76        let half = cast::round_to_usize(0.5 * window.get() * fs)
77            .unwrap_or(0)
78            .max(1);
79        let mut sums = Vec::with_capacity(samples.len() + 1);
80        sums.push(0.0);
81        let mut total = 0.0;
82        for x in samples {
83            total += x * x;
84            sums.push(total);
85        }
86        let db = frame_times(from, to)
87            .map(|t| {
88                let centre = cast::round_to_usize(t * fs)
89                    .unwrap_or(0)
90                    .saturating_sub(first);
91                let lo = centre.saturating_sub(half).min(samples.len());
92                let hi = (centre + half).min(samples.len());
93                if hi <= lo {
94                    return FLOOR_DB;
95                }
96                power_db((sums[hi] - sums[lo]) / cast::usize_to_f64(hi - lo))
97            })
98            .collect();
99        Self {
100            start: from,
101            db,
102            window,
103        }
104    }
105
106    /// Returns the RMS window of the track, 0 s when it has none.
107    #[must_use]
108    pub fn window(&self) -> Seconds {
109        self.window
110    }
111
112    /// The summed power of single partials in dB, from `from` to `to`.
113    ///
114    /// Each partial is a `(frequency, cutoff)` pair: the zero-phase band
115    /// magnitude of [`band_magnitude`] over `frequency ± cutoff`, so a
116    /// track has no lag. `samples` starts at `offset` in the stem.
117    ///
118    /// # Errors
119    ///
120    /// Returns a [`FilterError`] when a band reaches 0 Hz or the Nyquist
121    /// frequency.
122    pub fn partials(
123        samples: &[f64],
124        offset: Seconds,
125        rate: SampleRate,
126        (from, to): (Seconds, Seconds),
127        partials: &[(Hz, Hz)],
128    ) -> Result<Self, FilterError> {
129        let fs = rate.hz().get();
130        let mut power = vec![0.0; frame_times(from, to).count()];
131        for &(frequency, cutoff) in partials {
132            let magnitude = band_magnitude(
133                samples,
134                rate,
135                Hz::clamped(frequency.get() - cutoff.get()),
136                Hz::clamped(frequency.get() + cutoff.get()),
137            )?;
138            for (p, t) in power.iter_mut().zip(frame_times(from, to)) {
139                let i = cast::round_to_usize((t - offset.get()) * fs).unwrap_or(0);
140                if let Some(m) = magnitude.get(i) {
141                    // A partial of amplitude A reads A; its power is A^2 / 2,
142                    // as the RMS track reads it.
143                    *p += 0.5 * m * m;
144                }
145            }
146        }
147        Ok(Self {
148            start: from,
149            db: power.into_iter().map(power_db).collect(),
150            window: Seconds::clamped(0.0),
151        })
152    }
153
154    /// The magnitude of one harmonic in dB, from `from` to `to`: a sliding
155    /// DFT at `frequency` over one period of `f0`, centred on each frame.
156    ///
157    /// The window spans a whole period of the fundamental, so every other
158    /// harmonic of the note falls on a null, and its step response rises
159    /// without overshoot, unlike a Butterworth low-pass. A partial of
160    /// amplitude A reads A, its power A^2 / 2 as the RMS track reads it.
161    /// `samples` starts at sample `first` of the stem. At the start of the
162    /// samples the window keeps its whole period and runs from the first
163    /// sample, so the nulls hold; at their end it shrinks.
164    #[must_use]
165    pub fn harmonic(
166        samples: &[f64],
167        first: usize,
168        rate: SampleRate,
169        (from, to): (Seconds, Seconds),
170        frequency: Hz,
171        f0: Hz,
172    ) -> Self {
173        let fs = rate.hz().get();
174        let len = cast::round_to_usize(fs / f0.get().max(f64::MIN_POSITIVE))
175            .unwrap_or(1)
176            .max(1);
177        let step = std::f64::consts::TAU * frequency.get() / fs;
178        let (mut re, mut im) = (vec![0.0; samples.len() + 1], vec![0.0; samples.len() + 1]);
179        for (n, x) in samples.iter().enumerate() {
180            // Reduce the phase per sample so it stays exact on long stems.
181            let (sin, cos) =
182                ((step * cast::usize_to_f64(first + n)) % std::f64::consts::TAU).sin_cos();
183            re[n + 1] = re[n] + x * cos;
184            im[n + 1] = im[n] - x * sin;
185        }
186        let db = frame_times(from, to)
187            .map(|t| {
188                let centre = cast::round_to_usize(t * fs)
189                    .unwrap_or(0)
190                    .saturating_sub(first);
191                let lo = centre.saturating_sub(len / 2).min(samples.len());
192                let hi = (lo + len).min(samples.len());
193                if hi <= lo {
194                    return FLOOR_DB;
195                }
196                let n = cast::usize_to_f64(hi - lo);
197                let amplitude = 2.0 * (re[hi] - re[lo]).hypot(im[hi] - im[lo]) / n;
198                power_db(0.5 * amplitude * amplitude)
199            })
200            .collect();
201        Self {
202            start: from,
203            db,
204            window: Seconds::clamped(0.0),
205        }
206    }
207
208    /// Returns the time of the first value.
209    #[must_use]
210    pub fn start(&self) -> Seconds {
211        self.start
212    }
213
214    /// Returns the dB values.
215    #[must_use]
216    pub fn db(&self) -> &[f64] {
217        &self.db
218    }
219
220    /// Returns the dB values to change in place.
221    pub fn db_mut(&mut self) -> &mut [f64] {
222        &mut self.db
223    }
224
225    /// Returns the time of value `i`.
226    #[must_use]
227    pub fn time(&self, i: usize) -> f64 {
228        self.start.get() + cast::usize_to_f64(i) * HOP.get()
229    }
230
231    /// Returns the index of the value at or after time `t`, clamped to the
232    /// track.
233    #[must_use]
234    pub fn index(&self, t: f64) -> usize {
235        cast::ceil_to_usize((t - self.start.get()) / HOP.get() - 1e-9)
236            .unwrap_or(0)
237            .min(self.db.len())
238    }
239}
240
241/// The RMS window over `periods` periods of `f0`, or the fewest whole
242/// periods that span [`RMS_WINDOW_MIN`] when that is more.
243#[must_use]
244pub fn rms_window(f0: Hz, periods: f64) -> Seconds {
245    let f0 = f0.get().max(f64::MIN_POSITIVE);
246    let whole = (RMS_WINDOW_MIN.get() * f0).ceil();
247    Seconds::clamped(periods.max(whole) / f0)
248}
249
250/// The frame times from `from` up to `to`, one per [`HOP`].
251fn frame_times(from: Seconds, to: Seconds) -> impl Iterator<Item = f64> {
252    let n = cast::floor_to_usize((to.get() - from.get()) / HOP.get() + 1e-9).map_or(0, |n| n + 1);
253    (0..n).map(move |i| from.get() + cast::usize_to_f64(i) * HOP.get())
254}
255
256/// A mean power in dB, floored at digital silence.
257fn power_db(power: f64) -> f64 {
258    if power > 0.0 {
259        (10.0 * power.log10()).max(FLOOR_DB)
260    } else {
261        FLOOR_DB
262    }
263}
264
265#[cfg(test)]
266mod tests {
267    use super::*;
268    use std::f64::consts::TAU;
269
270    #[test]
271    fn the_constants_pass_the_unit_check() {
272        assert_eq!(HOP, Seconds::new(0.001).unwrap());
273        assert_eq!(RMS_WINDOW_MIN, Seconds::new(0.002).unwrap());
274        assert_eq!(HARMONIC_CUTOFF_MAX, Hz::new(40.0).unwrap());
275        assert_eq!(HARMONIC_CUTOFF_MIN, Hz::new(12.0).unwrap());
276    }
277
278    fn rate() -> SampleRate {
279        SampleRate::new(44_100).unwrap()
280    }
281
282    fn sine(f: f64, amplitude: f64, seconds: f64) -> Vec<f64> {
283        (0..cast::round_to_usize(seconds * 44_100.0).unwrap())
284            .map(|n| amplitude * (TAU * f * cast::usize_to_f64(n) / 44_100.0).sin())
285            .collect()
286    }
287
288    #[test]
289    fn a_steady_sine_reads_its_rms_with_no_ripple() {
290        let x = sine(110.0, 0.5, 0.5);
291        let track = LevelTrack::rms(
292            &x,
293            0,
294            rate(),
295            (Seconds::new(0.1).unwrap(), Seconds::new(0.4).unwrap()),
296            Hz::new(110.0).unwrap(),
297            RMS_PERIODS,
298        );
299        assert_eq!(track.db().len(), 301);
300        let expected = 10.0 * (0.125_f64).log10();
301        for db in track.db() {
302            assert!((db - expected).abs() < 0.05, "{db}");
303        }
304    }
305
306    #[test]
307    fn a_partial_track_reads_one_harmonic_of_two() {
308        let x: Vec<f64> = sine(110.0, 0.5, 0.6)
309            .iter()
310            .zip(sine(220.0, 0.25, 0.6))
311            .map(|(a, b)| a + b)
312            .collect();
313        let track = LevelTrack::partials(
314            &x,
315            Seconds::new(0.0).unwrap(),
316            rate(),
317            (Seconds::new(0.2).unwrap(), Seconds::new(0.4).unwrap()),
318            &[(Hz::new(220.0).unwrap(), Hz::new(40.0).unwrap())],
319        )
320        .unwrap();
321        let expected = 10.0 * (0.5 * 0.25 * 0.25_f64).log10();
322        for db in track.db() {
323            assert!((db - expected).abs() < 0.1, "{db}");
324        }
325    }
326
327    #[test]
328    fn a_harmonic_track_reads_one_harmonic_with_no_leak_from_the_others() {
329        let x: Vec<f64> = sine(110.0, 0.5, 0.6)
330            .iter()
331            .zip(sine(220.0, 0.25, 0.6))
332            .zip(sine(330.0, 0.4, 0.6))
333            .map(|((a, b), c)| a + b + c)
334            .collect();
335        let track = LevelTrack::harmonic(
336            &x,
337            0,
338            rate(),
339            (Seconds::new(0.2).unwrap(), Seconds::new(0.4).unwrap()),
340            Hz::new(220.0).unwrap(),
341            Hz::new(110.0).unwrap(),
342        );
343        let expected = 10.0 * (0.5 * 0.25 * 0.25_f64).log10();
344        for db in track.db() {
345            assert!((db - expected).abs() < 0.05, "{db}");
346        }
347    }
348
349    #[test]
350    fn a_harmonic_track_reads_a_whole_period_at_the_start_of_the_samples() {
351        // A 220 Hz sine from the first sample, with its second harmonic:
352        // the frame at 0 s reads the fundamental alone, at amplitude 1.
353        let fs = 22_050.0;
354        let x: Vec<f64> = (0..2205)
355            .map(|n| {
356                let t = cast::usize_to_f64(n) / fs;
357                (TAU * 220.0 * t).sin() + 0.5 * (TAU * 440.0 * t).sin()
358            })
359            .collect();
360        let track = LevelTrack::harmonic(
361            &x,
362            0,
363            SampleRate::new(22_050).unwrap(),
364            (Seconds::clamped(0.0), Seconds::clamped(0.01)),
365            Hz::new(220.0).unwrap(),
366            Hz::new(220.0).unwrap(),
367        );
368        let want = 10.0 * 0.5_f64.log10();
369        assert!((track.db()[0] - want).abs() < 0.1, "{}", track.db()[0]);
370    }
371
372    #[test]
373    fn a_harmonic_track_is_centred_on_its_frame() {
374        // A 110 Hz sine switched on at 0.2 s: the frame there sees half a
375        // window of it, so it reads half the amplitude, 6 dB down.
376        let mut x = sine(110.0, 0.5, 0.6);
377        for v in &mut x[..8820] {
378            *v = 0.0;
379        }
380        let track = LevelTrack::harmonic(
381            &x,
382            0,
383            rate(),
384            (Seconds::new(0.2).unwrap(), Seconds::new(0.4).unwrap()),
385            Hz::new(110.0).unwrap(),
386            Hz::new(110.0).unwrap(),
387        );
388        let full = 10.0 * (0.5 * 0.5 * 0.5_f64).log10();
389        assert!(
390            (track.db()[0] - (full - 6.02)).abs() < 0.3,
391            "{}",
392            track.db()[0]
393        );
394        assert!((track.db()[20] - full).abs() < 0.05, "{}", track.db()[20]);
395    }
396
397    #[test]
398    fn the_window_is_whole_periods_of_at_least_2_ms() {
399        let window = |f0: f64, periods: f64| rms_window(Hz::new(f0).unwrap(), periods).get();
400        assert!((window(110.0, 1.0) - 1.0 / 110.0).abs() < 1e-12);
401        assert!((window(110.0, 2.0) - 2.0 / 110.0).abs() < 1e-12);
402        // 1.14 ms periods at 880 Hz: two of them span 2.27 ms.
403        assert!((window(880.0, 1.0) - 2.0 / 880.0).abs() < 1e-12);
404    }
405
406    #[test]
407    fn silence_reads_the_floor_and_times_step_by_the_hop() {
408        let track = LevelTrack::rms(
409            &vec![0.0; 4410],
410            0,
411            rate(),
412            (Seconds::new(0.01).unwrap(), Seconds::new(0.02).unwrap()),
413            Hz::new(100.0).unwrap(),
414            RMS_PERIODS,
415        );
416        assert!(track.db().iter().all(|&db| db == FLOOR_DB));
417        assert!((track.time(10) - 0.02).abs() < 1e-12);
418        assert_eq!(track.index(0.0155), 6);
419        assert_eq!(track.index(0.0150), 5);
420    }
421}
sound_to_patch/estimate/envelope/track.rs

sound_to_patch/estimate/envelope/
track.rs
