sva_samples/measure/
peaks.rs

1// Concern: picks a Hann spectrum's peaks, those louder peaks' sidelobes cannot explain | Non-concern: computing the spectrum (spectrum.rs), naming a peak (pitch.rs) | IO: (magnitudes, span) -> peaks
2
3pub const FLOOR_DB: f64 = -160.0;
4
5#[derive(Clone, Copy, Debug, PartialEq)]
6pub struct Peak {
7    pub hz: f64,
8    pub db: f64,
9}
10
11pub fn db(amplitude: f64) -> f64 {
12    if amplitude <= 0.0 {
13        FLOOR_DB
14    } else {
15        (20.0 * amplitude.log10()).max(FLOOR_DB)
16    }
17}
18
19/// Local maxima, parabolically interpolated so a peak between two bins reports its real
20/// frequency instead of the nearest bin's.
21pub fn peaks(mags: &[f64], bin_hz: f64, max_peaks: usize) -> Vec<Peak> {
22    let mut found: Vec<Peak> = located(mags, bin_hz).into_iter().map(|(_, p)| p).collect();
23    found.truncate(max_peaks);
24    found
25}
26
27/// Loudest first by amplitude: `db` floors.
28fn located(mags: &[f64], bin_hz: f64) -> Vec<(usize, Peak)> {
29    let mut found: Vec<(f64, usize, Peak)> = Vec::new();
30    for k in 1..mags.len().saturating_sub(1) {
31        if mags[k] <= mags[k - 1] || mags[k] < mags[k + 1] || mags[k] <= 0.0 {
32            continue;
33        }
34        let (a, b, c) = (mags[k - 1], mags[k], mags[k + 1]);
35        let denom = a - 2.0 * b + c;
36        let delta = if denom.abs() > f64::EPSILON {
37            (0.5 * (a - c) / denom).clamp(-0.5, 0.5)
38        } else {
39            0.0
40        };
41        let amplitude = b - 0.25 * (a - c) * delta;
42        found.push((
43            amplitude,
44            k,
45            Peak {
46                hz: (k as f64 + delta) * bin_hz,
47                db: db(amplitude),
48            },
49        ));
50    }
51    found.sort_by(|x, y| y.0.total_cmp(&x.0));
52    found.into_iter().map(|(_, k, peak)| (k, peak)).collect()
53}
54
55/// Loudest first, peaks no kept peak's sidelobes explain, nor a DC's at bin 0.
56pub fn resolved(mags: &[f64], bin_hz: f64, spanned: usize, max_peaks: usize) -> Vec<Peak> {
57    resolved_counting(mags, bin_hz, spanned, max_peaks, |_| true)
58}
59
60/// As [`resolved`], until `max_peaks` kept pass `counts`.
61pub fn resolved_counting(
62    mags: &[f64],
63    bin_hz: f64,
64    spanned: usize,
65    max_peaks: usize,
66    counts: impl Fn(&Peak) -> bool,
67) -> Vec<Peak> {
68    let r = spanned as f64 / (2 * (mags.len() - 1)) as f64;
69    let half = r / 2.0;
70    let dc = match mags {
71        [zero, one, ..] if zero >= one => *zero,
72        _ => 0.0,
73    };
74    let mut kept: Vec<(usize, Peak)> = Vec::new();
75    let mut counted = 0;
76    for (k, peak) in located(mags, bin_hz) {
77        if counted == max_peaks {
78            break;
79        }
80        let leaked: f64 = kept
81            .iter()
82            .map(|&(q, _)| {
83                let near = (k.abs_diff(q) as f64 * r - half).max(0.0);
84                let image = (k + q) as f64 * r - half;
85                mags[q] / kernel(half) * (leakage(near) + leakage(image))
86            })
87            .sum::<f64>()
88            + dc / kernel(half) * leakage(k as f64 * r - half);
89        if mags[k] > leaked {
90            counted += usize::from(counts(&peak));
91            kept.push((k, peak));
92        }
93    }
94    kept.into_iter().map(|(_, peak)| peak).collect()
95}
96
97pub(crate) fn leakage(m: f64) -> f64 {
98    match m > 1.0 {
99        true => (1.0 / (std::f64::consts::PI * m * (m * m - 1.0))).min(1.0),
100        false => 1.0,
101    }
102}
103
104pub(crate) fn kernel(d: f64) -> f64 {
105    let d = d.abs();
106    match d {
107        0.0 => 1.0,
108        _ if d >= 1.0 => 0.5,
109        _ => {
110            let x = std::f64::consts::PI * d;
111            x.sin() / (x * (1.0 - d * d))
112        }
113    }
114}
sva_samples/measure/peaks.rs

sva_samples/measure/
peaks.rs
