sound_to_patch/dsp/
descriptors.rs

1//! Global spectral and level descriptors of a stem, each with one named
2//! convention.
3//!
4//! Shape descriptors use magnitude weights, band shares use power. The
5//! published tools disagree on weights, floors and rolloff arguments, and
6//! the conventions differ by up to 15x on the same sound, so a value is
7//! comparable only together with its convention.
8
9use crate::audio::SampleRate;
10use crate::cast;
11use crate::confidence::{ConfidenceError, Estimate, Level, Reason, ReasonCode, Signal, Spread};
12use crate::units::{Decibels, Hz, Normalized, Ratio};
13
14use super::frames::{ActiveFrames, MIN_ACTIVE_FRAMES, Stats, Summary};
15use super::{Stft, StftParams, Window};
16
17/// The convention of the shape descriptors.
18pub const SHAPE_CONVENTION: &str = "magnitude weights |X_k| divided by the window sum, \
19periodic Hann frames, statistics over active frames as median and interquartile range";
20/// The convention of the band shares.
21pub const BAND_CONVENTION: &str = "power |X_k|^2 summed over active frames per band, divided \
22by the power in the bands reported; bands start at 20, 60, 250, 1000, 4000 and 10000 Hz and end \
23at the next start or 20 kHz; a band that starts at or above the Nyquist frequency is left out, \
24and one that crosses it ends there";
25/// The convention of the flatness.
26pub const FLATNESS_CONVENTION: &str = "geometric over arithmetic mean of the magnitude, \
27floored at 1e-5 after the window-sum scale, bins from 100 Hz to 10 kHz";
28/// The convention of the levels.
29pub const LEVEL_CONVENTION: &str = "peak and RMS of the mono mix over the whole stem in dB \
30relative to full scale 1.0; a full-scale sine reads -3.01 dB RMS; crest is peak minus RMS";
31
32/// The band edges of the band shares in hertz. The edges are a project
33/// choice.
34pub const BAND_EDGES_HZ: [Hz; 7] = [
35    Hz::from_const(20.0),
36    Hz::from_const(60.0),
37    Hz::from_const(250.0),
38    Hz::from_const(1000.0),
39    Hz::from_const(4000.0),
40    Hz::from_const(10_000.0),
41    Hz::from_const(20_000.0),
42];
43
44/// The flatness floor, applied after the window-sum scale.
45pub const FLATNESS_FLOOR: Normalized = Normalized::from_const(1e-5);
46
47const FLATNESS_BAND_HZ: (f64, f64) = (100.0, 10_000.0);
48
49/// The descriptors of one frame.
50#[derive(Clone, Copy, Debug, PartialEq)]
51pub struct FrameDescriptors {
52    /// `sum f_k p_k` with `p_k = a_k / sum a_k`.
53    pub centroid: Hz,
54    /// `sqrt(sum (f_k - centroid)^2 p_k)`.
55    pub spread: Hz,
56    /// The first bin whose cumulative magnitude reaches 85 % of the total.
57    pub rolloff85: Hz,
58    /// The same at 95 %.
59    pub rolloff95: Hz,
60    /// Geometric over arithmetic mean of the floored magnitude in 100 Hz to
61    /// 10 kHz.
62    pub flatness: Normalized,
63}
64
65/// Computes the descriptors of one magnitude frame, already divided by the
66/// window sum, with bins `bin_hz` apart. Returns `None` for an all-zero
67/// frame, or for one with a negative, NaN or infinite magnitude.
68#[must_use]
69pub fn frame_descriptors(magnitudes: &[f64], bin_hz: Hz) -> Option<FrameDescriptors> {
70    if magnitudes.iter().any(|a| !a.is_finite() || *a < 0.0) {
71        return None;
72    }
73    let total: f64 = magnitudes.iter().sum();
74    if !total.is_finite() || total <= 0.0 {
75        return None;
76    }
77    let width = bin_hz.get();
78    let frequency = |k: usize| cast::usize_to_f64(k) * width;
79    let centroid = magnitudes
80        .iter()
81        .enumerate()
82        .map(|(k, a)| frequency(k) * a)
83        .sum::<f64>()
84        / total;
85    let spread = (magnitudes
86        .iter()
87        .enumerate()
88        .map(|(k, a)| (frequency(k) - centroid).powi(2) * a)
89        .sum::<f64>()
90        / total)
91        .sqrt();
92    let rolloff = |share: f64| {
93        let mut cumulative = 0.0;
94        let bin = magnitudes
95            .iter()
96            .position(|a| {
97                cumulative += a;
98                cumulative >= share * total
99            })
100            .unwrap_or(magnitudes.len() - 1);
101        Hz::clamped(frequency(bin))
102    };
103    let band: Vec<f64> = magnitudes
104        .iter()
105        .enumerate()
106        .filter(|&(k, _)| (FLATNESS_BAND_HZ.0..=FLATNESS_BAND_HZ.1).contains(&frequency(k)))
107        .map(|(_, a)| a.max(FLATNESS_FLOOR.get()))
108        .collect();
109    let flatness = if band.is_empty() {
110        1.0
111    } else {
112        let count = cast::usize_to_f64(band.len());
113        let log_mean = band.iter().map(|a| a.ln()).sum::<f64>() / count;
114        let mean = band.iter().sum::<f64>() / count;
115        log_mean.exp() / mean
116    };
117    Some(FrameDescriptors {
118        centroid: Hz::clamped(centroid),
119        spread: Hz::clamped(spread),
120        rolloff85: rolloff(0.85),
121        rolloff95: rolloff(0.95),
122        flatness: Normalized::clamped(flatness),
123    })
124}
125
126/// The share of power in one band.
127#[derive(Clone, Copy, Debug, PartialEq)]
128pub struct BandShare {
129    /// The lower edge, inclusive.
130    pub lo: Hz,
131    /// The upper edge, exclusive.
132    pub hi: Hz,
133    /// The share of the power in all bands.
134    pub share: Normalized,
135}
136
137/// The descriptors of a stem.
138#[derive(Clone, Debug, PartialEq)]
139pub struct Descriptors {
140    /// The frame length `N` of the analysis.
141    pub frame_len: usize,
142    /// The number of active frames.
143    pub active_frames: usize,
144    /// The spectral centroid, magnitude weights.
145    pub centroid: Estimate<Stats<Hz>, Hz>,
146    /// The centroid over the median f0 of the note list.
147    pub centroid_over_f0: Estimate<Stats<Ratio>, Ratio>,
148    /// The spectral spread around the centroid, magnitude weights.
149    pub spread: Estimate<Stats<Hz>, Hz>,
150    /// The 85 % magnitude rolloff.
151    pub rolloff85: Estimate<Stats<Hz>, Hz>,
152    /// The 95 % magnitude rolloff.
153    pub rolloff95: Estimate<Stats<Hz>, Hz>,
154    /// The spectral flatness, magnitude, 100 Hz to 10 kHz.
155    pub flatness: Estimate<Stats<Normalized>, Normalized>,
156    /// The power shares of the bands between [`BAND_EDGES_HZ`].
157    pub band_share: Estimate<Vec<BandShare>, Normalized>,
158    /// The peak level of the mono mix in dB relative to full scale.
159    pub peak_dbfs: Estimate<Decibels, Decibels>,
160    /// The RMS level of the mono mix in dB relative to full scale.
161    pub rms_dbfs: Estimate<Decibels, Decibels>,
162    /// Peak over RMS in dB.
163    pub crest_db: Estimate<Decibels, Decibels>,
164}
165
166/// Computes the descriptors of a mono signal.
167///
168/// Frames are periodic Hann of `params.frame_len()` samples with hop
169/// `frame_len / 4`. `f0` is the median fundamental of the note list, if one
170/// was given. `hard_limited` is the flat-top flag of the stem; the levels
171/// carry it as a reason.
172///
173/// # Errors
174///
175/// Returns a [`ConfidenceError`] only if a reason cannot be built, which
176/// the counts here never cause.
177pub fn describe(
178    samples: &[f64],
179    rate: SampleRate,
180    params: StftParams,
181    f0: Option<Hz>,
182    hard_limited: bool,
183) -> Result<Descriptors, ConfidenceError> {
184    let params = params.with_window(Window::Hann);
185    let active = ActiveFrames::of(samples, params);
186    let count = u64::try_from(active.count()).unwrap_or(u64::MAX);
187    let evidence = Evidence {
188        frames: Reason::new(ReasonCode::ActiveFrames, Signal::Count(count))?,
189        silent: Reason::new(ReasonCode::SilentStem, Signal::Flag(true))?,
190        short: (active.count() > 0 && active.count() < MIN_ACTIVE_FRAMES)
191            .then(too_few_frames)
192            .transpose()?,
193    };
194    let (per_frame, band_power) = frame_values(samples, rate, params, &active);
195    let pick = |f: fn(&FrameDescriptors) -> f64| per_frame.iter().map(f).collect::<Vec<_>>();
196    let ratios = f0.filter(|f| f.get() > 0.0).map(|f0| {
197        pick(|d| d.centroid.get())
198            .iter()
199            .map(|c| c / f0.get())
200            .collect::<Vec<_>>()
201    });
202    let centroid_over_f0 = match ratios {
203        Some(ratios) => evidence.stats(&ratios, Level::Medium, Ratio::clamped)?,
204        None => Estimate::none(vec![Reason::new(
205            ReasonCode::NoNoteList,
206            Signal::Flag(true),
207        )?])?,
208    };
209    let peak = samples.iter().fold(0.0_f64, |max, x| max.max(x.abs()));
210    let rms = if samples.is_empty() {
211        0.0
212    } else {
213        (samples.iter().map(|x| x * x).sum::<f64>() / cast::usize_to_f64(samples.len())).sqrt()
214    };
215    let peak_db = Decibels::from_amplitude_gain(peak).ok();
216    let rms_db = Decibels::from_amplitude_gain(rms).ok();
217    let crest = peak_db
218        .zip(rms_db)
219        .map(|(p, r)| Decibels::clamped(p.get() - r.get()));
220    Ok(Descriptors {
221        frame_len: params.frame_len(),
222        active_frames: active.count(),
223        centroid: evidence.stats(&pick(|d| d.centroid.get()), Level::High, Hz::clamped)?,
224        centroid_over_f0,
225        spread: evidence.stats(&pick(|d| d.spread.get()), Level::Medium, Hz::clamped)?,
226        rolloff85: evidence.stats(&pick(|d| d.rolloff85.get()), Level::Medium, Hz::clamped)?,
227        rolloff95: evidence.stats(&pick(|d| d.rolloff95.get()), Level::Medium, Hz::clamped)?,
228        flatness: evidence.stats(
229            &pick(|d| d.flatness.get()),
230            Level::Medium,
231            Normalized::clamped,
232        )?,
233        band_share: evidence.band_share(&band_power, rate.nyquist())?,
234        peak_dbfs: evidence.level(peak_db, hard_limited)?,
235        rms_dbfs: evidence.level(rms_db, hard_limited)?,
236        crest_db: evidence.level(crest, hard_limited)?,
237    })
238}
239
240/// The descriptors of each active frame and the power per band summed over
241/// them.
242fn frame_values(
243    samples: &[f64],
244    rate: SampleRate,
245    params: StftParams,
246    active: &ActiveFrames,
247) -> (Vec<FrameDescriptors>, Vec<f64>) {
248    let window_sum: f64 = Window::Hann.samples(params.frame_len()).iter().sum();
249    let bin_hz = params.bin_width(rate);
250    let mut stft = Stft::new(params);
251    let mut per_frame = Vec::with_capacity(active.count());
252    let mut band_power = vec![0.0; BAND_EDGES_HZ.len() - 1];
253    for frame in active.indices() {
254        let magnitudes: Vec<f64> = stft
255            .frame(samples, frame)
256            .iter()
257            .map(|value| value.norm() / window_sum)
258            .collect();
259        for (k, a) in magnitudes.iter().enumerate() {
260            let f = cast::usize_to_f64(k) * bin_hz.get();
261            if let Some(band) = BAND_EDGES_HZ
262                .windows(2)
263                .position(|edge| (edge[0].get()..edge[1].get()).contains(&f))
264            {
265                band_power[band] += a * a;
266            }
267        }
268        if let Some(descriptors) = frame_descriptors(&magnitudes, bin_hz) {
269            per_frame.push(descriptors);
270        }
271    }
272    (per_frame, band_power)
273}
274
275/// The reason of a statistic with fewer than [`MIN_ACTIVE_FRAMES`] frames.
276pub(crate) fn too_few_frames() -> Result<Reason, ConfidenceError> {
277    Reason::new(
278        ReasonCode::TooFewFrames,
279        Signal::Count(u64::try_from(MIN_ACTIVE_FRAMES).unwrap_or(u64::MAX)),
280    )
281}
282
283/// The reasons every descriptor cites.
284struct Evidence {
285    frames: Reason,
286    silent: Reason,
287    /// Set when the active frames are too few for a statistic.
288    short: Option<Reason>,
289}
290
291impl Evidence {
292    /// The median and IQR of per-frame values, the IQR as the spread; no
293    /// estimate without a frame or with fewer than [`MIN_ACTIVE_FRAMES`].
294    fn stats<T: Spread>(
295        &self,
296        values: &[f64],
297        level: Level,
298        unit: fn(f64) -> T,
299    ) -> Result<Estimate<Stats<T>, T>, ConfidenceError> {
300        if let Some(short) = self.short {
301            return Estimate::none(vec![self.frames, short]);
302        }
303        match Summary::of(values) {
304            Some(summary) => Estimate::new(
305                summary.to_stats(unit),
306                level,
307                Some(unit(summary.iqr())),
308                vec![self.frames],
309            ),
310            None => Estimate::none(vec![self.silent, self.frames]),
311        }
312    }
313
314    /// The share of each band. A band that starts at or above the Nyquist
315    /// frequency holds no bin and is left out; a band that crosses it ends
316    /// there.
317    fn band_share(
318        &self,
319        band_power: &[f64],
320        nyquist: Hz,
321    ) -> Result<Estimate<Vec<BandShare>, Normalized>, ConfidenceError> {
322        if let Some(short) = self.short {
323            return Estimate::none(vec![self.frames, short]);
324        }
325        let kept: Vec<(&[Hz], f64)> = BAND_EDGES_HZ
326            .windows(2)
327            .zip(band_power.iter().copied())
328            .filter(|(edge, _)| edge[0] < nyquist)
329            .collect();
330        // Over the kept bands, so the shares sum to 1; a band that starts
331        // exactly at the Nyquist frequency held only its bin.
332        let total: f64 = kept.iter().map(|(_, power)| power).sum();
333        if total <= 0.0 {
334            return Estimate::none(vec![self.silent, self.frames]);
335        }
336        let shares = kept
337            .into_iter()
338            .map(|(edge, power)| BandShare {
339                lo: edge[0],
340                hi: if edge[1] < nyquist { edge[1] } else { nyquist },
341                share: Normalized::clamped(power / total),
342            })
343            .collect();
344        Estimate::new(shares, Level::High, None, vec![self.frames])
345    }
346
347    fn level(
348        &self,
349        value: Option<Decibels>,
350        hard_limited: bool,
351    ) -> Result<Estimate<Decibels, Decibels>, ConfidenceError> {
352        match value {
353            Some(db) => Estimate::new(
354                db,
355                Level::High,
356                None,
357                vec![Reason::new(
358                    ReasonCode::HardLimited,
359                    Signal::Flag(hard_limited),
360                )?],
361            ),
362            None => Estimate::none(vec![self.silent]),
363        }
364    }
365}
366
367#[cfg(test)]
368mod tests {
369    use super::*;
370
371    #[test]
372    fn the_constants_pass_the_unit_check() {
373        let edges =
374            [20.0, 60.0, 250.0, 1000.0, 4000.0, 10_000.0, 20_000.0].map(|v| Hz::new(v).unwrap());
375        assert_eq!(BAND_EDGES_HZ, edges);
376        assert_eq!(FLATNESS_FLOOR, Normalized::new(1e-5).unwrap());
377    }
378
379    #[test]
380    fn a_frame_with_a_bad_magnitude_has_no_descriptors() {
381        let bin = Hz::new(10.0).unwrap();
382        for bad in [f64::NAN, f64::INFINITY, -1.0] {
383            assert_eq!(frame_descriptors(&[1.0, bad, 1.0], bin), None, "{bad}");
384        }
385        assert_eq!(frame_descriptors(&[f64::MAX, f64::MAX], bin), None);
386    }
387
388    #[test]
389    fn a_frame_with_a_zero_bin_has_a_finite_flatness() {
390        let mut magnitudes = vec![0.01; 1025];
391        magnitudes[300] = 0.0;
392        let d = frame_descriptors(&magnitudes, Hz::new(21.533).unwrap()).unwrap();
393        assert!(d.flatness.get() > 0.0 && d.flatness.get().is_finite());
394    }
395
396    #[test]
397    fn a_single_bin_puts_centroid_and_rolloffs_on_it() {
398        let mut magnitudes = vec![0.0; 100];
399        magnitudes[10] = 1.0;
400        let d = frame_descriptors(&magnitudes, Hz::new(10.0).unwrap()).unwrap();
401        assert_eq!(d.centroid.get(), 100.0);
402        assert_eq!(d.spread.get(), 0.0);
403        assert_eq!((d.rolloff85.get(), d.rolloff95.get()), (100.0, 100.0));
404        assert!(frame_descriptors(&[0.0; 4], Hz::new(10.0).unwrap()).is_none());
405    }
406
407    #[test]
408    fn silence_has_no_estimate_with_a_reason() {
409        let rate = SampleRate::new(44_100).unwrap();
410        let d = describe(
411            &vec![0.0; 44_100],
412            rate,
413            StftParams::new(4096).unwrap(),
414            None,
415            false,
416        )
417        .unwrap();
418        assert_eq!(d.active_frames, 0);
419        assert!(d.centroid.value().is_none());
420        assert!(d.peak_dbfs.value().is_none());
421        assert!(d.crest_db.value().is_none());
422        assert_eq!(
423            d.centroid_over_f0.confidence().reasons()[0].code(),
424            ReasonCode::NoNoteList
425        );
426    }
427}
sound_to_patch/dsp/descriptors.rs

sound_to_patch/dsp/
descriptors.rs
