sva_samples/measure/
crest.rs1use sva_formula::filter::Shape;
4
5use crate::biquad::{State, design};
6use crate::measure::envelope::rms;
7use crate::measure::spectrum::{db, third_octave_edges};
8
9const COUNTED_UNDER_DB: f64 = 60.0;
11
12const 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#[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
72fn 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}