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 and nearest-edge Gaussian
34/// smoothing at 1/12 octave, before choosing a graph level reference.
35pub fn smoothed_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    (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}
68
69/// Mean of a measurement-grid curve over 200–2000 Hz.
70pub fn midrange_reference(curve: &[f64]) -> f64 {
71    let mid: Vec<_> = measurement_grid()
72        .into_iter()
73        .zip(curve)
74        .filter(|(f, _)| (200.0..=2000.0).contains(f))
75        .map(|(_, v)| *v)
76        .collect();
77    pairwise_sum(&mid) / mid.len() as f64
78}
79
80/// Smoothed measurement response with a zero mean over 200–2000 Hz, exactly
81/// as upstream step 5 displays each curve independently.
82pub fn magnitude_db(samples: &[f64], sample_rate: f64) -> Vec<f64> {
83    let mut smoothed = smoothed_db(samples, sample_rate);
84    let reference = midrange_reference(&smoothed);
85    smoothed.iter_mut().for_each(|v| *v -= reference);
86    smoothed
87}
88
89/// Pre-calibration position-one SPL, using the measurement's native sample rate
90/// and a minimum 65536-sample FFT, with 1/24-octave power smoothing.
91pub fn measurement_spl(samples: &[f64], sample_rate: f64) -> Vec<f64> {
92    let (f, h) = spectrum(samples, sample_rate, 65536);
93    let power: Vec<_> = h.iter().map(|v| v.re.hypot(v.im).powi(2)).collect();
94    let db: Vec<_> = smooth_power(&f, &power, Fraction::Scalar(24.0))
95        .expect("FFT grid")
96        .iter()
97        .map(|p| 10.0 * (p + 1e-30).log10())
98        .collect();
99    spl_grid().iter().map(|&g| interp(g, &f, &db)).collect()
100}
101
102/// NumPy's linear percentile, ignoring nonfinite values for graph axis ranges.
103pub fn percentile(values: &[f64], percent: f64) -> f64 {
104    let mut values: Vec<_> = values.iter().copied().filter(|v| v.is_finite()).collect();
105    if values.is_empty() {
106        return 0.0;
107    }
108    values.sort_by(f64::total_cmp);
109    let index = (values.len() - 1) as f64 * percent / 100.0;
110    let lo = index.floor() as usize;
111    let hi = index.ceil() as usize;
112    values[lo] + (values[hi] - values[lo]) * index.fract()
113}
114
115#[cfg(test)]
116mod tests {
117    use super::*;
118
119    #[test]
120    fn split_measurement_transform_preserves_upstream_normalization() {
121        let samples: Vec<_> = (0..4096)
122            .map(|index| {
123                let time = index as f64 / 48_000.0;
124                (2.0 * std::f64::consts::PI * 997.0 * time).sin() * (-30.0 * time).exp()
125            })
126            .collect();
127        let mut split = smoothed_db(&samples, 48_000.0);
128        let reference = midrange_reference(&split);
129        split.iter_mut().for_each(|value| *value -= reference);
130        assert_eq!(split, magnitude_db(&samples, 48_000.0));
131    }
132}
secs_engine/display.rs

secs_engine/
display.rs
