Skip to main content

sva_samples/measure/
crest.rs

1// Concern: peak minus RMS in each third-octave band, and the spread across them | Non-concern: the discrete spectrum (spectrum.rs), the ERB envelope bank (bands.rs) | IO: (&[f64], sample rate) -> Crest
2
3use sva_formula::filter::Shape;
4
5use crate::biquad::{State, design};
6use crate::measure::envelope::rms;
7use crate::measure::spectrum::{db, third_octave_edges};
8
9/// Below this a band holds leakage, which would set the spread.
10const COUNTED_UNDER_DB: f64 = 60.0;
11
12/// Five time constants of overshoot, and `sqrt(sqrt(2)-1)` is what a second identical
13/// stage leaves of one stage's -3 dB width.
14const SETTLE_TAUS: f64 = 5.0;
15const CASCADE_Q: f64 = 1.55;
16
17#[derive(Clone, Copy, Debug, PartialEq)]
18pub struct BandCrest {
19    pub lo_hz: f64,
20    pub hi_hz: f64,
21    pub centre_hz: f64,
22    pub peak: f64,
23    pub rms: f64,
24    pub crest_db: f64,
25    pub counted: bool,
26}
27
28/// The spread, not the level: the worst-rated master had the most low-band range.
29#[derive(Clone, Debug, PartialEq)]
30pub struct Crest {
31    pub broadband_crest_db: f64,
32    pub spread_db: Option<f64>,
33    pub widest_band_hz: Option<f64>,
34    pub tightest_band_hz: Option<f64>,
35    pub counted_under_db: f64,
36    pub bands: Vec<BandCrest>,
37}
38
39pub fn analyze(samples: &[f64], sample_rate: f64) -> Crest {
40    let mut bands: Vec<BandCrest> = third_octave_edges(sample_rate)
41        .into_iter()
42        .map(|(lo, hi)| band(samples, sample_rate, lo, hi))
43        .collect();
44
45    let loudest = bands.iter().map(|b| b.rms).fold(0.0, f64::max);
46    let floor = loudest * 10f64.powf(-COUNTED_UNDER_DB / 20.0);
47    for b in &mut bands {
48        b.counted = b.rms > floor && b.rms > 0.0;
49    }
50
51    let counted: Vec<&BandCrest> = bands.iter().filter(|b| b.counted).collect();
52    let widest = counted
53        .iter()
54        .max_by(|a, b| a.crest_db.total_cmp(&b.crest_db));
55    let tightest = counted
56        .iter()
57        .min_by(|a, b| a.crest_db.total_cmp(&b.crest_db));
58
59    Crest {
60        broadband_crest_db: crest_db(
61            samples.iter().fold(0.0, |a, x| a.max(x.abs())),
62            rms(samples),
63        ),
64        spread_db: widest.zip(tightest).map(|(w, t)| w.crest_db - t.crest_db),
65        widest_band_hz: widest.map(|b| b.centre_hz),
66        tightest_band_hz: tightest.map(|b| b.centre_hz),
67        counted_under_db: COUNTED_UNDER_DB,
68        bands,
69    }
70}
71
72/// A crest is a RATIO within one band, so the cascade's gain closed form cancels out of it.
73fn band(samples: &[f64], sample_rate: f64, lo: f64, hi: f64) -> BandCrest {
74    let centre = (lo * hi).sqrt();
75    let q = centre / (hi - lo);
76    let coeffs = design(Shape::Bandpass, centre, q, 0.0, sample_rate);
77    let settle = settle_secs(centre, q) * sample_rate;
78    let skip = (settle.round() as usize).min(samples.len());
79    let (mut first, mut second) = (State::default(), State::default());
80    let mut peak = 0.0f64;
81    let mut power = 0.0f64;
82    for (i, x) in samples.iter().enumerate() {
83        let y = second.step(&coeffs, first.step(&coeffs, *x));
84        if i < skip {
85            continue;
86        }
87        peak = peak.max(y.abs());
88        power += y * y;
89    }
90    let measured = samples.len() - skip;
91    let rms = match measured {
92        0 => 0.0,
93        n => (power / n as f64).sqrt(),
94    };
95    BandCrest {
96        lo_hz: lo,
97        hi_hz: hi,
98        centre_hz: centre,
99        peak,
100        rms,
101        crest_db: crest_db(peak, rms),
102        counted: false,
103    }
104}
105
106fn settle_secs(centre_hz: f64, q: f64) -> f64 {
107    SETTLE_TAUS * CASCADE_Q * q / (std::f64::consts::PI * centre_hz)
108}
109
110fn crest_db(peak: f64, rms: f64) -> f64 {
111    if rms <= 0.0 { db(0.0) } else { db(peak / rms) }
112}