sound_to_patch/estimate/envelope/
tail.rs

1//! The tail after a note-off (#70): a linear release plus a decaying tail.
2
3use super::POOL_NOTES_MIN;
4use super::adsr::NoteEnvelope;
5use super::track::HOP;
6use crate::audio::{Band, SampleRate, stereo_screen};
7use crate::cast;
8use crate::dsp::{StftParams, Summary, describe};
9use crate::estimate::pitch::fit::nelder_mead;
10use crate::notes::Note;
11use crate::units::{Decibels, Seconds};
12
13/// A tail RT60 counts only when the frames above the floor span this share of
14/// it: a fall of 15 dB.
15const RT60_SPAN_SHARE: f64 = 0.25;
16/// The fit reads the tail down to this far under the sustain level.
17pub const TAIL_FLOOR: Decibels = Decibels::from_const(-55.0);
18/// The fewest frames of a tail worth a fit.
19const TAIL_FRAMES_MIN: usize = 8;
20/// The starting releases, in seconds.
21const RELEASE_STARTS: [f64; 4] = [0.02, 0.08, 0.3, 1.0];
22/// The starting tail RT60s, in seconds.
23const RT60_STARTS: [f64; 3] = [0.4, 1.2, 3.0];
24/// The starting tail weight.
25const WEIGHT_START: f64 = 0.1;
26/// A tail reads every frame up to here.
27const DENSE_FRAMES: usize = 200;
28/// Past [`DENSE_FRAMES`] a tail reads every this many frames.
29const SPARSE_STEP: usize = 5;
30/// The simplex steps of one fit.
31const FIT_ITERATIONS: usize = 400;
32/// The releases a fit can read, in seconds. A level that never falls would
33/// otherwise run the fit to a release of 1e19 s and more.
34pub const TAIL_RELEASE_RANGE: (f64, f64) = (0.001, 10.0);
35/// The tail weights a fit can read.
36pub const TAIL_WEIGHT_RANGE: (f64, f64) = (1e-4, 10.0);
37/// The RT60s a fit can read, in seconds.
38pub const TAIL_RT60_RANGE: (f64, f64) = (0.01, 30.0);
39
40/// The tail model fitted to a note after its note-off.
41#[derive(Clone, Copy, Debug, PartialEq)]
42pub struct TailFit {
43    /// The linear release.
44    pub release: Seconds,
45    /// The weight of the tail at the note-off, over the sustain level.
46    pub weight: f64,
47    /// The RT60 of the tail.
48    pub rt60: Seconds,
49    /// The frames fitted span enough of the RT60 to read it.
50    pub rt60_observed: bool,
51    /// The RMS of the fit in dB, over the frames fitted.
52    pub residual: Decibels,
53    /// The shortest RT60 that the filter of the track can show.
54    pub rt60_floor: Seconds,
55    /// The median size of the second difference of the held level in dB, over
56    /// [`TAIL_HOLD_FRAMES`] before the knee, [`TAIL_HOLD_LAG`] frames each way.
57    pub hold_bend: Decibels,
58}
59
60impl TailFit {
61    /// The tail is heavy but fits too smoothly for a reverb, whose level wobbles:
62    /// the note's own exponential release, with no linear part to find.
63    #[must_use]
64    pub fn is_smooth(&self) -> bool {
65        self.weight > TAIL_WEIGHT_MAX && self.residual.get() <= TAIL_SMOOTH.get()
66    }
67
68    /// The fit shows an RT60 of a tail: one heavy enough, read long enough, past
69    /// what the filter shows by itself, and not a release of the note's own.
70    #[must_use]
71    pub fn reads_rt60(&self) -> bool {
72        self.weight >= TAIL_WEIGHT_MIN
73            && self.rt60_observed
74            && self.rt60.get() >= self.rt60_floor.get()
75            && !self.is_smooth()
76    }
77}
78
79/// The model in dB, `t` seconds after the note-off, no lower than the floor.
80fn model_db(t: f64, release: f64, weight: f64, rt60: f64) -> f64 {
81    let amplitude = (1.0 - t / release).max(0.0) + weight * 10f64.powf(-3.0 * t / rt60);
82    (20.0 * amplitude.max(1e-12).log10()).max(TAIL_FLOOR.get())
83}
84
85/// Fits the model to `db`, the level in dB over the sustain level on frames
86/// from the note-off. `None` when too few frames lie above the floor.
87#[must_use]
88pub fn fit_tail(db: &[f64]) -> Option<TailFit> {
89    let floor = TAIL_FLOOR.get();
90    let last = db.iter().rposition(|v| *v > floor)?;
91    let above = db.iter().filter(|v| **v > floor).count();
92    if above < TAIL_FRAMES_MIN {
93        return None;
94    }
95    // Past the release's first frames a tail moves slowly, so a sparser read
96    // costs less without moving the fit.
97    let data: Vec<(f64, f64, f64)> = db[..=last]
98        .iter()
99        .enumerate()
100        .filter(|(i, _)| *i < DENSE_FRAMES || i % SPARSE_STEP == 0)
101        .map(|(i, v)| {
102            let weight = if i < DENSE_FRAMES {
103                1.0
104            } else {
105                cast::usize_to_f64(SPARSE_STEP)
106            };
107            (cast::usize_to_f64(i) * HOP.get(), v.max(floor), weight)
108        })
109        .collect();
110    let frames: f64 = data.iter().map(|d| d.2).sum();
111    let bounded = |p: &[f64; 3]| {
112        (
113            p[0].exp().clamp(TAIL_RELEASE_RANGE.0, TAIL_RELEASE_RANGE.1),
114            p[1].exp().clamp(TAIL_WEIGHT_RANGE.0, TAIL_WEIGHT_RANGE.1),
115            p[2].exp().clamp(TAIL_RT60_RANGE.0, TAIL_RT60_RANGE.1),
116        )
117    };
118    let cost = |p: &[f64; 3]| -> f64 {
119        let (release, weight, rt60) = bounded(p);
120        data.iter()
121            .map(|&(t, v, w)| w * (v - model_db(t, release, weight, rt60)).powi(2))
122            .sum()
123    };
124    let mut best: Option<([f64; 3], f64)> = None;
125    for release in RELEASE_STARTS {
126        for rt60 in RT60_STARTS {
127            let start = [release.ln(), WEIGHT_START.ln(), rt60.ln()];
128            let (p, value) = nelder_mead(cost, start, [0.3, 0.5, 0.3], FIT_ITERATIONS, 1e-6);
129            if best.is_none_or(|(_, b)| value < b) {
130                best = Some((p, value));
131            }
132        }
133    }
134    let (p, value) = best?;
135    let (release, weight, rt60) = bounded(&p);
136    Some(TailFit {
137        release: Seconds::clamped(release),
138        weight,
139        rt60: Seconds::clamped(rt60),
140        rt60_observed: cast::usize_to_f64(above) * HOP.get() >= RT60_SPAN_SHARE * rt60,
141        residual: Decibels::clamped((value / frames).sqrt()),
142        rt60_floor: Seconds::clamped(0.0),
143        hold_bend: Decibels::clamped(0.0),
144    })
145}
146
147/// A tail lighter than this at the note-off has no RT60 to read: on dry notes
148/// the fit often leaves a little weight, and its RT60 is noise.
149pub const TAIL_WEIGHT_MIN: f64 = 0.05;
150/// Over this tail weight the release is not identifiable, and no release is
151/// written.
152pub const TAIL_WEIGHT_MAX: f64 = 0.5;
153/// Over this residual the release is not identifiable, and no release is
154/// written.
155pub const TAIL_RESIDUAL_MAX: Decibels = Decibels::from_const(3.0);
156/// A written release reads low when the median residual is over this and the
157/// median [`TailFit::hold_bend`] over [`TAIL_HOLD_BEND`]: a tremolo carries on
158/// through the release, which the fit cannot follow.
159pub const TAIL_WOBBLE_RESIDUAL: Decibels = Decibels::from_const(2.0);
160/// See [`TAIL_WOBBLE_RESIDUAL`].
161pub const TAIL_HOLD_BEND: Decibels = Decibels::from_const(0.3);
162/// The frames before the knee that [`TailFit::hold_bend`] reads.
163pub const TAIL_HOLD_FRAMES: usize = 300;
164/// The frames each way of [`TailFit::hold_bend`]. Slower than the knee's bend, to
165/// see a tremolo under 10 Hz.
166pub const TAIL_HOLD_LAG: usize = 25;
167/// A tail over [`TAIL_WEIGHT_MAX`] that fits within this RMS reads as the note's
168/// own exponential release: a reverb's level wobbles more.
169pub const TAIL_SMOOTH: Decibels = Decibels::from_const(1.1);
170/// The shortest RT60 a band can show is this many periods of its cutoff: a dry
171/// chord's release through a band reads 3 to 5 periods, whatever its length up
172/// to 50 ms.
173pub const TAIL_BAND_PERIODS: f64 = 6.0;
174/// A tail outlasts a gap when its level before the next note lies above this
175/// level under the sustain level.
176pub const GAP_LEVEL: Decibels = Decibels::from_const(-40.0);
177/// Side over mid grows by this much from the note body to the tail.
178pub const SIDE_GROWTH: Decibels = Decibels::from_const(3.0);
179/// The flatness rises by this much from the note body to the tail.
180pub const FLATNESS_RISE: f64 = 0.1;
181/// The tail rate holds over the notes when its interquartile range stays under
182/// this share of its median.
183pub const SAME_RATE_SPREAD: f64 = 0.25;
184/// The share of the bands from 1 kHz falls by this much from the first half of
185/// a tail to the second.
186pub const HIGH_BAND_FALL: Decibels = Decibels::from_const(3.0);
187/// The body of a note runs up to this long before the knee.
188const BODY_SPAN: Seconds = Seconds::from_const(0.2);
189/// The level before the next note reads this span.
190const GAP_SPAN: Seconds = Seconds::from_const(0.02);
191/// The frames of a vote's spectrum.
192const VOTE_FRAME: usize = 1024;
193
194/// The votes that a tail is reverb: `None` when nothing could be read.
195#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, serde::Serialize)]
196pub struct TailVotes {
197    /// The tail outlasts the silent gap before the next note.
198    pub outlasts_gap: Option<bool>,
199    /// Side over mid grows from the note body to the tail.
200    pub side_grows: Option<bool>,
201    /// Spectral flatness rises in the tail.
202    pub flatness_rises: Option<bool>,
203    /// The tail rate is the same for every note.
204    pub same_rate: Option<bool>,
205    /// The high band decays faster than the low band.
206    pub high_band_faster: Option<bool>,
207}
208
209impl TailVotes {
210    /// Returns the votes in favour.
211    #[must_use]
212    pub fn count(self) -> usize {
213        [
214            self.outlasts_gap,
215            self.side_grows,
216            self.flatness_rises,
217            self.same_rate,
218            self.high_band_faster,
219        ]
220        .into_iter()
221        .filter(|v| *v == Some(true))
222        .count()
223    }
224}
225
226/// The share of the votes reads as a majority; `None` without any.
227fn majority(values: &[bool]) -> Option<bool> {
228    (!values.is_empty()).then(|| 2 * values.iter().filter(|v| **v).count() > values.len())
229}
230
231fn rms_db(x: &[f64]) -> Option<f64> {
232    let energy = x.iter().map(|v| v * v).sum::<f64>();
233    (energy > 0.0).then(|| 10.0 * (energy / cast::usize_to_f64(x.len())).log10())
234}
235
236/// The sample at `t` seconds in a signal of `len` samples. A time past the
237/// last index maps to `len`.
238fn sample_index(t: f64, fs: f64, len: usize) -> usize {
239    cast::round_to_usize(t.max(0.0) * fs)
240        .unwrap_or(len)
241        .min(len)
242}
243
244/// Reads the votes over the pooled notes that have a tail fit.
245pub(super) fn read_votes(
246    mono: &[f64],
247    stereo: Option<(&[f64], &[f64])>,
248    rate: SampleRate,
249    notes: &[Note],
250    pool: &[&NoteEnvelope],
251) -> TailVotes {
252    let index = |t: f64| sample_index(t, rate.hz().get(), mono.len());
253    let mut gaps = Vec::new();
254    let mut sides = Vec::new();
255    let mut flats = Vec::new();
256    let mut highs = Vec::new();
257    let mut rates = Vec::new();
258    for note in pool {
259        let (Some(knee), Some(fit)) = (note.knee.map(Seconds::get), note.tail) else {
260            continue;
261        };
262        let next = notes.iter().map(|n| n.onset().get()).find(|t| *t > knee);
263        let stop = next.map_or(note.end.get(), |n| n.min(note.end.get()));
264        let body = index((knee - BODY_SPAN.get()).max(note.onset.get()))..index(knee);
265        let tail = index(knee + fit.release.get())..index(stop);
266        if fit.reads_rt60() {
267            rates.push(fit.rt60.get());
268        }
269        // Only a gap with no other note in it shows a tail. The fitted note is the
270        // first of its onset and pitch, so a second one counts as another note.
271        let own = notes
272            .iter()
273            .position(|n| n.onset() == note.onset && n.pitch() == note.pitch);
274        let silent_gap = |next: f64| {
275            !notes.iter().enumerate().any(|(i, n)| {
276                Some(i) != own && n.onset().get() < next && n.offset().get() > next - GAP_SPAN.get()
277            })
278        };
279        if let Some(next) = next
280            .filter(|n| *n - knee > fit.release.get() + GAP_SPAN.get())
281            .filter(|n| silent_gap(*n))
282        {
283            let before = index(next - GAP_SPAN.get())..index(next);
284            if let (Some(gap), Some(level)) = (rms_db(&mono[before]), rms_db(&mono[body.clone()])) {
285                gaps.push(gap - level > GAP_LEVEL.get());
286            }
287        }
288        if tail.is_empty() || body.is_empty() {
289            continue;
290        }
291        if let Some((left, right)) = stereo {
292            let sm = |r: std::ops::Range<usize>| {
293                stereo_screen(left.get(r.clone())?, right.get(r)?, rate, None)
294                    .ok()
295                    .and_then(|s| s.band(Band::Full).sm_db)
296                    .map(Decibels::get)
297            };
298            if let (Some(a), Some(b)) = (sm(body.clone()), sm(tail.clone())) {
299                sides.push(b - a >= SIDE_GROWTH.get());
300            }
301        }
302        let spectra = |r: std::ops::Range<usize>| {
303            describe(
304                &mono[r],
305                rate,
306                StftParams::new(VOTE_FRAME).ok()?,
307                None,
308                false,
309            )
310            .ok()
311        };
312        if let (Some(a), Some(b)) = (spectra(body), spectra(tail.clone()))
313            && let (Some(a), Some(b)) = (a.flatness.value(), b.flatness.value())
314        {
315            flats.push(b.median().get() - a.median().get() >= FLATNESS_RISE);
316        }
317        let middle = tail.start + tail.len() / 2;
318        if let (Some(a), Some(b)) = (spectra(tail.start..middle), spectra(middle..tail.end)) {
319            let high = |d: &crate::dsp::Descriptors| {
320                d.band_share.value().map(|bands| {
321                    bands
322                        .iter()
323                        .filter(|p| p.lo.get() >= 1000.0)
324                        .map(|p| p.share.get())
325                        .sum::<f64>()
326                })
327            };
328            if let (Some(a), Some(b)) = (high(&a), high(&b))
329                && a > 0.0
330                && b > 0.0
331            {
332                highs.push(10.0 * (b / a).log10() <= -HIGH_BAND_FALL.get());
333            }
334        }
335    }
336    let same_rate = (rates.len() >= POOL_NOTES_MIN)
337        .then(|| Summary::of(&rates).is_some_and(|s| s.iqr() <= SAME_RATE_SPREAD * s.median()));
338    TailVotes {
339        outlasts_gap: majority(&gaps),
340        side_grows: majority(&sides),
341        flatness_rises: majority(&flats),
342        same_rate,
343        high_band_faster: majority(&highs),
344    }
345}
346
347#[cfg(test)]
348mod tests {
349    use super::*;
350
351    fn curve(release: f64, weight: f64, rt60: f64) -> Vec<f64> {
352        (0..4000)
353            .map(|i| model_db(cast::usize_to_f64(i) * HOP.get(), release, weight, rt60))
354            .collect()
355    }
356
357    #[test]
358    fn the_fit_recovers_the_model_it_reads() {
359        for (release, weight, rt60) in [(0.05, 0.15, 0.6), (0.2, 0.4, 1.5), (0.02, 0.05, 3.0)] {
360            let fit = fit_tail(&curve(release, weight, rt60)).unwrap();
361            assert!((fit.release.get() / release - 1.0).abs() < 0.03, "{fit:?}");
362            assert!((fit.weight / weight - 1.0).abs() < 0.03, "{fit:?}");
363            assert!((fit.rt60.get() / rt60 - 1.0).abs() < 0.03, "{fit:?}");
364            assert!(fit.residual.get() < 0.1, "{fit:?}");
365        }
366    }
367
368    #[test]
369    fn a_level_that_never_falls_stays_inside_the_ranges() {
370        let fit = fit_tail(&[0.0; 200]).unwrap();
371        let inside = |v: f64, (lo, hi): (f64, f64)| (lo..=hi).contains(&v);
372        assert!(inside(fit.release.get(), TAIL_RELEASE_RANGE), "{fit:?}");
373        assert!(inside(fit.weight, TAIL_WEIGHT_RANGE), "{fit:?}");
374        assert!(inside(fit.rt60.get(), TAIL_RT60_RANGE), "{fit:?}");
375        assert!(fit.residual.get() < 0.1, "{fit:?}");
376        let loud = fit_tail(&[40.0; 200]).unwrap();
377        assert!(inside(loud.weight, TAIL_WEIGHT_RANGE), "{loud:?}");
378    }
379
380    #[test]
381    fn a_time_past_the_last_index_maps_to_the_end() {
382        assert_eq!(sample_index(1e30, 44_100.0, 100), 100);
383        assert_eq!(sample_index(1.0, 50.0, 100), 50);
384        assert_eq!(sample_index(-1.0, 50.0, 100), 0);
385    }
386
387    #[test]
388    fn a_tail_under_the_floor_or_too_short_has_no_fit() {
389        assert_eq!(fit_tail(&[-60.0; 100]), None);
390        let mut one_frame = [-60.0; 100];
391        one_frame[50] = -1.0;
392        assert_eq!(fit_tail(&one_frame), None);
393        assert_eq!(fit_tail(&[-1.0; 7]), None);
394        assert!(fit_tail(&[-1.0; 8]).is_some());
395    }
396
397    #[test]
398    fn an_rt60_counts_over_a_fall_of_15_db_of_it() {
399        let short: Vec<f64> = curve(0.05, 0.3, 3.0)[..500].to_vec();
400        assert!(!fit_tail(&short).unwrap().rt60_observed);
401        assert!(fit_tail(&curve(0.05, 0.3, 3.0)).unwrap().rt60_observed);
402        // 0.9 s reads a quarter of 3 s, and not a half.
403        assert!(
404            fit_tail(&curve(0.05, 0.3, 3.0)[..900])
405                .unwrap()
406                .rt60_observed
407        );
408        // One late frame above the floor adds no span.
409        let mut spike = vec![-5.0; 300];
410        spike.resize(1500, -60.0);
411        spike[1499] = -5.0;
412        assert!(!fit_tail(&spike).unwrap().rt60_observed);
413    }
414
415    #[test]
416    fn the_votes_count_only_those_in_favour() {
417        let votes = TailVotes {
418            outlasts_gap: Some(true),
419            side_grows: Some(false),
420            flatness_rises: None,
421            same_rate: Some(true),
422            high_band_faster: None,
423        };
424        assert_eq!(votes.count(), 2);
425        assert_eq!(TailVotes::default().count(), 0);
426    }
427
428    #[test]
429    fn a_majority_needs_more_than_half() {
430        assert_eq!(majority(&[]), None);
431        assert_eq!(majority(&[true, false]), Some(false));
432        assert_eq!(majority(&[true, true, false]), Some(true));
433    }
434}
sound_to_patch/estimate/envelope/tail.rs

sound_to_patch/estimate/envelope/
tail.rs
