1use secs_dsp::{
4 fft, interp, pairwise_sum,
5 smooth::{Fraction, smooth_power},
6};
7
8pub fn measurement_grid() -> Vec<f64> {
10 crate::design::geomspace(5.0, 22000.0, 280)
11}
12
13pub fn spl_grid() -> Vec<f64> {
15 crate::design::geomspace(5.0, 20000.0, 1390)
16}
17
18pub 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
33pub 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
78pub 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
91pub 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}