sound_to_patch/dsp/
decay.rs

1//! Decay times from backward integration of the energy after a stop.
2//!
3//! A stem has no impulse response, so the Schroeder integral runs on the
4//! squared signal of a tail after an abrupt stop, such as the last note-off
5//! of a phrase. The integral smooths the random fluctuation of a noisy
6//! tail. Bleed raises the floor, so the caller takes the stop from the note
7//! list, never from any silence.
8
9use std::error::Error;
10use std::fmt;
11
12use crate::audio::SampleRate;
13use crate::cast;
14use crate::units::{Decibels, Seconds};
15
16/// Length of the blocks that estimate the level and the noise floor.
17const BLOCK_SECONDS: f64 = 0.01;
18/// Share of the tail at its end that estimates the noise floor.
19const FLOOR_SHARE: f64 = 0.1;
20/// The floor must lie this far below the bottom of a fit range.
21const FLOOR_MARGIN_DB: f64 = 10.0;
22
23/// The level range of a decay fit, in decibels below the start of the
24/// energy decay curve.
25#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
26pub enum FitRange {
27    /// From -5 dB to -25 dB.
28    T20,
29    /// From -5 dB to -35 dB.
30    T30,
31}
32
33impl FitRange {
34    /// Returns the top and the bottom of the range.
35    #[must_use]
36    pub fn limits(self) -> (Decibels, Decibels) {
37        let (top, bottom) = match self {
38            Self::T20 => (-5.0, -25.0),
39            Self::T30 => (-5.0, -35.0),
40        };
41        (Decibels::from_const(top), Decibels::from_const(bottom))
42    }
43}
44
45/// Why a decay fit has no value.
46#[derive(Clone, Copy, Debug, PartialEq)]
47pub enum DecayFailure {
48    /// The tail is digital silence.
49    Silent,
50    /// The tail is too short to fit a line.
51    TooShort,
52    /// The noise floor sits less than 10 dB below the bottom of the range.
53    FloorTooHigh {
54        /// The floor relative to the level at the stop.
55        floor: Decibels,
56        /// The highest floor the range accepts.
57        limit: Decibels,
58    },
59    /// The decay curve does not fall to the bottom of the range before the
60    /// tail meets the floor.
61    RangeNotReached,
62}
63
64impl fmt::Display for DecayFailure {
65    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
66        match self {
67            Self::Silent => write!(f, "the tail is silent"),
68            Self::TooShort => write!(f, "the tail is too short for a decay fit"),
69            Self::FloorTooHigh { floor, limit } => write!(
70                f,
71                "the noise floor at {floor} is above the {limit} that the fit range needs"
72            ),
73            Self::RangeNotReached => {
74                write!(f, "the decay does not reach the bottom of the fit range")
75            }
76        }
77    }
78}
79
80impl Error for DecayFailure {}
81
82/// The decay of one tail: the noise floor and a decay time per fit range.
83#[derive(Clone, Debug, PartialEq)]
84pub struct Decay {
85    floor: Option<Decibels>,
86    truncation: Seconds,
87    t20: Result<Seconds, DecayFailure>,
88    t30: Result<Seconds, DecayFailure>,
89}
90
91impl Decay {
92    /// Analyses the samples of a tail that starts at a stop.
93    ///
94    /// The analysis starts at the loudest 10 ms block, which is the first
95    /// block unless the tail opens with silence or a pre-delay. That block
96    /// is the reference level.
97    ///
98    /// 1. The noise floor is the mean power of the last 10 % of the tail,
99    ///    relative to the power of the reference block.
100    /// 2. A line fitted to the 10 ms block levels down to 10 dB above the
101    ///    floor gives the time where the decay meets the floor.
102    /// 3. The Schroeder integral runs from each sample to that time, plus
103    ///    the energy the decay would add after it at the fitted rate.
104    /// 4. A line fitted to the decay curve between the range limits gives
105    ///    the decay time, `-60 / slope`: the time a 60 dB drop takes.
106    ///
107    /// A fit fails when the floor lies less than 10 dB below the bottom of
108    /// its range or the curve never reaches that bottom.
109    #[must_use]
110    pub fn of(tail: &[f64], rate: SampleRate) -> Self {
111        match Self::analyse(tail, rate) {
112            Ok(decay) => decay,
113            Err(failure) => Self {
114                floor: None,
115                truncation: rate.seconds(tail.len()),
116                t20: Err(failure),
117                t30: Err(failure),
118            },
119        }
120    }
121
122    fn analyse(tail: &[f64], rate: SampleRate) -> Result<Self, DecayFailure> {
123        let fs = rate.hz().get();
124        let block = cast::round_to_usize(BLOCK_SECONDS * fs).unwrap_or(1).max(1);
125        if tail.iter().all(|&x| x == 0.0) {
126            return Err(DecayFailure::Silent);
127        }
128        let loudest = tail
129            .chunks_exact(block)
130            .map(|chunk| chunk.iter().map(|x| x * x).sum::<f64>())
131            .enumerate()
132            .fold(
133                (0, 0.0),
134                |best, (index, sum)| if sum > best.1 { (index, sum) } else { best },
135            )
136            .0;
137        let offset = loudest * block;
138        Self::analyse_from(&tail[offset..], rate, block, offset)
139    }
140
141    /// Analyses a tail whose first block is the reference and that starts
142    /// `offset` samples after the stop.
143    fn analyse_from(
144        tail: &[f64],
145        rate: SampleRate,
146        block: usize,
147        offset: usize,
148    ) -> Result<Self, DecayFailure> {
149        let fs = rate.hz().get();
150        let energy: Vec<f64> = tail.iter().map(|x| x * x).collect();
151        if energy.len() < 3 * block {
152            return Err(DecayFailure::TooShort);
153        }
154        let levels: Vec<f64> = energy
155            .chunks_exact(block)
156            .map(|chunk| chunk.iter().sum::<f64>() / cast::usize_to_f64(block))
157            .collect();
158        let reference = levels[0];
159        if reference <= 0.0 {
160            return Err(DecayFailure::Silent);
161        }
162        let floor_start = energy.len()
163            - cast::ceil_to_usize(FLOOR_SHARE * cast::usize_to_f64(energy.len()))
164                .unwrap_or(block)
165                .max(block);
166        let noise = energy[floor_start..].iter().sum::<f64>()
167            / cast::usize_to_f64(energy.len() - floor_start);
168        let floor = (noise > 0.0).then(|| power_db(noise / reference));
169
170        // Block levels from the start down to 10 dB above the floor.
171        let floor_db = floor.unwrap_or(f64::NEG_INFINITY);
172        let points: Vec<(f64, f64)> = levels
173            .iter()
174            .enumerate()
175            .map(|(index, &level)| {
176                let centre = (cast::usize_to_f64(index) + 0.5) * cast::usize_to_f64(block) / fs;
177                (
178                    centre,
179                    if level > 0.0 {
180                        power_db(level / reference)
181                    } else {
182                        f64::NEG_INFINITY
183                    },
184                )
185            })
186            .take_while(|&(_, level)| level >= floor_db + FLOOR_MARGIN_DB)
187            .collect();
188        if points.len() < 2 {
189            return Err(DecayFailure::FloorTooHigh {
190                floor: Decibels::clamped(floor_db),
191                limit: Decibels::clamped(-FLOOR_MARGIN_DB),
192            });
193        }
194        let (slope, intercept) = fit_line(&points).ok_or(DecayFailure::TooShort)?;
195        let end = cast::usize_to_f64(energy.len()) / fs;
196        let (cut, compensation) = if slope < 0.0 && floor.is_some() {
197            let crossing = ((floor_db - intercept) / slope).clamp(0.0, end);
198            // Power decays as exp(-t / tau) with tau = 10 / (ln 10 |slope|).
199            let tau = 10.0 / (std::f64::consts::LN_10 * slope.abs());
200            (crossing, noise * tau * fs)
201        } else {
202            (end, 0.0)
203        };
204        let cut_index = cast::floor_to_usize(cut * fs)
205            .unwrap_or(energy.len())
206            .clamp(1, energy.len());
207        let mut curve = backward_integral(&energy[..cut_index]);
208        for value in &mut curve {
209            *value += compensation;
210        }
211        Ok(Self {
212            floor: floor.map(Decibels::clamped),
213            truncation: rate.seconds(offset + cut_index),
214            t20: fit_range(&curve, fs, floor, FitRange::T20),
215            t30: fit_range(&curve, fs, floor, FitRange::T30),
216        })
217    }
218
219    /// Returns the noise floor relative to the level at the stop, or `None`
220    /// when the end of the tail is digital silence.
221    #[must_use]
222    pub fn floor(&self) -> Option<Decibels> {
223        self.floor
224    }
225
226    /// Returns the time after the stop where the integral was truncated.
227    #[must_use]
228    pub fn truncation(&self) -> Seconds {
229        self.truncation
230    }
231
232    /// Returns the decay time from the -5 to -25 dB fit.
233    ///
234    /// # Errors
235    ///
236    /// Returns the [`DecayFailure`] of the fit.
237    pub fn t20(&self) -> Result<Seconds, DecayFailure> {
238        self.t20
239    }
240
241    /// Returns the decay time from the -5 to -35 dB fit.
242    ///
243    /// # Errors
244    ///
245    /// Returns the [`DecayFailure`] of the fit.
246    pub fn t30(&self) -> Result<Seconds, DecayFailure> {
247        self.t30
248    }
249}
250
251fn power_db(ratio: f64) -> f64 {
252    10.0 * ratio.log10()
253}
254
255/// The Schroeder integral of an energy sequence: element `n` is the sum
256/// from `n` to the end. Summed backward, so the small late terms add first.
257#[must_use]
258pub fn backward_integral(energy: &[f64]) -> Vec<f64> {
259    let mut curve = vec![0.0; energy.len()];
260    let mut sum = 0.0;
261    for (slot, &value) in curve.iter_mut().zip(energy).rev() {
262        sum += value;
263        *slot = sum;
264    }
265    curve
266}
267
268fn fit_range(
269    curve: &[f64],
270    fs: f64,
271    floor: Option<f64>,
272    range: FitRange,
273) -> Result<Seconds, DecayFailure> {
274    let (top, bottom) = range.limits();
275    let (top, bottom) = (top.get(), bottom.get());
276    if let Some(floor) = floor
277        && floor > bottom - FLOOR_MARGIN_DB
278    {
279        return Err(DecayFailure::FloorTooHigh {
280            floor: Decibels::clamped(floor),
281            limit: Decibels::clamped(bottom - FLOOR_MARGIN_DB),
282        });
283    }
284    let start = curve
285        .first()
286        .copied()
287        .filter(|&v| v > 0.0)
288        .ok_or(DecayFailure::Silent)?;
289    let levels: Vec<f64> = curve.iter().map(|&v| power_db(v / start)).collect();
290    if levels.last().is_none_or(|&last| last > bottom) {
291        return Err(DecayFailure::RangeNotReached);
292    }
293    let points: Vec<(f64, f64)> = levels
294        .iter()
295        .enumerate()
296        .filter(|&(_, &level)| level <= top && level >= bottom)
297        .map(|(n, &level)| (cast::usize_to_f64(n) / fs, level))
298        .collect();
299    let (slope, _) = fit_line(&points).ok_or(DecayFailure::TooShort)?;
300    if slope >= 0.0 {
301        return Err(DecayFailure::RangeNotReached);
302    }
303    Ok(Seconds::clamped(-60.0 / slope))
304}
305
306/// Least-squares line through `(x, y)` points: `(slope, intercept)`.
307pub(crate) fn fit_line(points: &[(f64, f64)]) -> Option<(f64, f64)> {
308    if points.len() < 2 {
309        return None;
310    }
311    let n = cast::usize_to_f64(points.len());
312    let mean_x = points.iter().map(|p| p.0).sum::<f64>() / n;
313    let mean_y = points.iter().map(|p| p.1).sum::<f64>() / n;
314    let sxx: f64 = points.iter().map(|p| (p.0 - mean_x).powi(2)).sum();
315    let sxy: f64 = points.iter().map(|p| (p.0 - mean_x) * (p.1 - mean_y)).sum();
316    (sxx > 0.0).then(|| {
317        let slope = sxy / sxx;
318        (slope, mean_y - slope * mean_x)
319    })
320}
321
322#[cfg(test)]
323mod tests {
324    use super::*;
325
326    #[test]
327    fn the_fit_limits_pass_the_unit_check() {
328        let db = |v| Decibels::new(v).unwrap();
329        assert_eq!(FitRange::T20.limits(), (db(-5.0), db(-25.0)));
330        assert_eq!(FitRange::T30.limits(), (db(-5.0), db(-35.0)));
331    }
332
333    fn rate() -> SampleRate {
334        SampleRate::new(16_000).unwrap()
335    }
336
337    fn exponential(rt60: f64, seconds: f64) -> Vec<f64> {
338        let fs = 16_000.0;
339        (0..cast::floor_to_usize(seconds * fs).unwrap())
340            .map(|n| (-6.907_755_278_982_137 * cast::usize_to_f64(n) / fs / rt60).exp())
341            .collect()
342    }
343
344    #[test]
345    fn the_backward_integral_sums_to_the_end() {
346        assert_eq!(backward_integral(&[1.0, 2.0, 3.0]), vec![6.0, 5.0, 3.0]);
347    }
348
349    #[test]
350    fn a_pure_exponential_decays_60_db_per_rt60() {
351        // Amplitude exp(-ln(1000) t / RT60) is -60 dB at RT60 exactly.
352        let decay = Decay::of(&exponential(0.4, 2.0), rate());
353        for time in [decay.t20().unwrap(), decay.t30().unwrap()] {
354            assert!((time.get() - 0.4).abs() < 0.4e-3, "{time:?}");
355        }
356    }
357
358    #[test]
359    fn silence_and_short_tails_fail() {
360        assert_eq!(
361            Decay::of(&vec![0.0; 16_000], rate()).t20(),
362            Err(DecayFailure::Silent)
363        );
364        assert_eq!(
365            Decay::of(&[1.0; 100], rate()).t30(),
366            Err(DecayFailure::TooShort)
367        );
368    }
369
370    #[test]
371    fn a_flat_tail_never_reaches_the_range() {
372        let decay = Decay::of(&vec![0.5; 16_000], rate());
373        assert!(
374            matches!(decay.t20(), Err(DecayFailure::FloorTooHigh { .. })),
375            "{decay:?}"
376        );
377    }
378}
sound_to_patch/dsp/decay.rs

sound_to_patch/dsp/
decay.rs
