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: (planes, sample rate) -> Crest
2
3use sva_formula::filter::Shape;
4
5use crate::biquad::{State, design};
6use crate::measure::envelope::level;
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(planes: &[&[f64]], sample_rate: f64) -> Crest {
40    crate::measure::held(planes);
41    let mut bands: Vec<BandCrest> = third_octave_edges(sample_rate)
42        .into_iter()
43        .map(|(lo, hi)| band(planes, sample_rate, lo, hi))
44        .collect();
45
46    let loudest = bands.iter().map(|b| b.rms).fold(0.0, f64::max);
47    let floor = loudest * 10f64.powf(-COUNTED_UNDER_DB / 20.0);
48    for b in &mut bands {
49        b.counted = b.rms > floor && b.rms > 0.0;
50    }
51
52    let counted: Vec<&BandCrest> = bands.iter().filter(|b| b.counted).collect();
53    let widest = counted
54        .iter()
55        .max_by(|a, b| a.crest_db.total_cmp(&b.crest_db));
56    let tightest = counted
57        .iter()
58        .min_by(|a, b| a.crest_db.total_cmp(&b.crest_db));
59
60    Crest {
61        broadband_crest_db: crest_db(
62            planes
63                .iter()
64                .flat_map(|p| p.iter())
65                .fold(0.0, |a, x| a.max(x.abs())),
66            level(planes),
67        ),
68        spread_db: widest.zip(tightest).map(|(w, t)| w.crest_db - t.crest_db),
69        widest_band_hz: widest.map(|b| b.centre_hz),
70        tightest_band_hz: tightest.map(|b| b.centre_hz),
71        counted_under_db: COUNTED_UNDER_DB,
72        bands,
73    }
74}
75
76/// A crest is a RATIO within one band, so the cascade's gain closed form cancels out of it.
77fn band(planes: &[&[f64]], sample_rate: f64, lo: f64, hi: f64) -> BandCrest {
78    let centre = (lo * hi).sqrt();
79    let q = centre / (hi - lo);
80    let coeffs = design(Shape::Bandpass, centre, q, 0.0, sample_rate);
81    let settle = settle_secs(centre, q) * sample_rate;
82    let len = planes[0].len();
83    let skip = (settle.round() as usize).min(len);
84    let mut peak = 0.0f64;
85    let powers: Vec<f64> = planes
86        .iter()
87        .map(|samples| {
88            let (mut first, mut second) = (State::default(), State::default());
89            let mut power = 0.0f64;
90            for (i, x) in samples.iter().enumerate() {
91                let y = second.step(&coeffs, first.step(&coeffs, *x));
92                if i < skip {
93                    continue;
94                }
95                peak = peak.max(y.abs());
96                power += y * y;
97            }
98            power
99        })
100        .collect();
101    let power = powers[1..].iter().fold(powers[0], |sum, p| sum + p);
102    let rms = match (len - skip) * planes.len() {
103        0 => 0.0,
104        n => (power / n as f64).sqrt(),
105    };
106    BandCrest {
107        lo_hz: lo,
108        hi_hz: hi,
109        centre_hz: centre,
110        peak,
111        rms,
112        crest_db: crest_db(peak, rms),
113        counted: false,
114    }
115}
116
117fn settle_secs(centre_hz: f64, q: f64) -> f64 {
118    SETTLE_TAUS * CASCADE_Q * q / (std::f64::consts::PI * centre_hz)
119}
120
121fn crest_db(peak: f64, rms: f64) -> f64 {
122    if rms <= 0.0 { db(0.0) } else { db(peak / rms) }
123}
sva_samples/measure/crest.rs

sva_samples/measure/
crest.rs
