Skip to main content

secs_engine/
display.rs

1//! Display-only transforms from `step5_measure.magnitude_db` and the pre-run SPL preview.
2
3use secs_dsp::{
4    fft, interp, pairwise_sum,
5    smooth::{Fraction, smooth_power},
6};
7
8/// Measurement plot's logarithmic 280-point grid, in Hz.
9pub fn measurement_grid() -> Vec<f64> {
10    crate::design::geomspace(5.0, 22000.0, 280)
11}
12
13/// Calibration SPL plot's 1390-point grid, in Hz.
14pub fn spl_grid() -> Vec<f64> {
15    crate::design::geomspace(5.0, 20000.0, 1390)
16}
17
18/// Calibration GD plot's 1070-point grid, in Hz.
19pub fn gd_grid() -> Vec<f64> {
20    crate::design::geomspace(5.0, 5000.0, 1070)
21}
22
23fn spectrum(samples: &[f64], sample_rate: f64, minimum: usize) -> (Vec<f64>, Vec<crate::C>) {
24    let n = samples.len().max(minimum).next_power_of_two();
25    let mut padded = vec![0.0; n];
26    padded[..samples.len()].copy_from_slice(samples);
27    let frequencies = (0..=n / 2)
28        .map(|i| i as f64 * sample_rate / n as f64)
29        .collect();
30    (frequencies, fft::rfft(&padded))
31}
32
33/// Unwindowed measurement FFT, log interpolation, nearest-edge Gaussian smoothing
34/// at 1/12 octave and a zero mean over 200–2000 Hz, exactly as step 5 displays it.
35pub fn magnitude_db(samples: &[f64], sample_rate: f64) -> Vec<f64> {
36    let grid = measurement_grid();
37    if samples.is_empty() {
38        return vec![0.0; grid.len()];
39    }
40    let (frequencies, h) = spectrum(samples, sample_rate, 8192);
41    let log_f: Vec<_> = frequencies.iter().map(|f| f.max(1e-6).log10()).collect();
42    let db: Vec<_> = h
43        .iter()
44        .map(|v| 20.0 * v.re.hypot(v.im).max(1e-12).log10())
45        .collect();
46    let sampled: Vec<_> = grid
47        .iter()
48        .map(|f| interp(f.log10(), &log_f, &db))
49        .collect();
50    let sigma = (grid.len() - 1) as f64 / (grid[grid.len() - 1] / grid[0]).log2() / (12.0 * 2.3548);
51    let radius = (4.0 * sigma + 0.5) as isize;
52    let mut weights: Vec<_> = (-radius..=radius)
53        .map(|k| (-0.5 * (k as f64 / sigma).powi(2)).exp())
54        .collect();
55    let total = pairwise_sum(&weights);
56    weights.iter_mut().for_each(|w| *w /= total);
57    let mut smoothed: Vec<_> = (0..grid.len())
58        .map(|i| {
59            (-radius..=radius)
60                .zip(&weights)
61                .map(|(k, w)| {
62                    sampled[(i as isize + k).clamp(0, sampled.len() as isize - 1) as usize] * w
63                })
64                .sum::<f64>()
65        })
66        .collect();
67    let mid: Vec<_> = grid
68        .iter()
69        .zip(&smoothed)
70        .filter(|(f, _)| (200.0..=2000.0).contains(*f))
71        .map(|(_, v)| *v)
72        .collect();
73    let reference = pairwise_sum(&mid) / mid.len() as f64;
74    smoothed.iter_mut().for_each(|v| *v -= reference);
75    smoothed
76}
77
78/// Pre-calibration position-one SPL, using the measurement's native sample rate
79/// and a minimum 65536-sample FFT, with 1/24-octave power smoothing.
80pub fn measurement_spl(samples: &[f64], sample_rate: f64) -> Vec<f64> {
81    let (f, h) = spectrum(samples, sample_rate, 65536);
82    let power: Vec<_> = h.iter().map(|v| v.re.hypot(v.im).powi(2)).collect();
83    let db: Vec<_> = smooth_power(&f, &power, Fraction::Scalar(24.0))
84        .expect("FFT grid")
85        .iter()
86        .map(|p| 10.0 * (p + 1e-30).log10())
87        .collect();
88    spl_grid().iter().map(|&g| interp(g, &f, &db)).collect()
89}
90
91/// NumPy's linear percentile, ignoring nonfinite values for graph axis ranges.
92pub fn percentile(values: &[f64], percent: f64) -> f64 {
93    let mut values: Vec<_> = values.iter().copied().filter(|v| v.is_finite()).collect();
94    if values.is_empty() {
95        return 0.0;
96    }
97    values.sort_by(f64::total_cmp);
98    let index = (values.len() - 1) as f64 * percent / 100.0;
99    let lo = index.floor() as usize;
100    let hi = index.ceil() as usize;
101    values[lo] + (values[hi] - values[lo]) * index.fract()
102}