sva_samples/measure/
spectrum.rs

1// Concern: resolves every channel's samples into magnitudes, third-octave bands and peaks | Non-concern: naming a peak as a note (pitch.rs), choosing the window | IO: (planes, sample rate) -> Spectrum
2
3use crate::fft::fft;
4
5const MAX_FRAME: usize = 8192;
6const MIN_FRAME: usize = 64;
7
8/// One FFT's memory, which is all that bounds a caller's own frame.
9pub const MAX_PINNED_FRAME: usize = 1 << 20;
10
11/// A caller checks this against [`MAX_PINNED_FRAME`]: a shrunk frame is not a pinned one.
12pub fn pinned_frame(secs: f64, sample_rate: f64) -> usize {
13    let asked = (secs * sample_rate).round().max(0.0);
14    if !asked.is_finite() || asked > MAX_PINNED_FRAME as f64 {
15        return MAX_PINNED_FRAME * 2;
16    }
17    (asked as usize).next_power_of_two().max(MIN_FRAME)
18}
19const FLOOR_DB: f64 = -160.0;
20
21#[derive(Clone, Copy, Debug, PartialEq)]
22pub struct Band {
23    pub lo_hz: f64,
24    pub hi_hz: f64,
25    pub db: f64,
26}
27
28#[derive(Clone, Copy, Debug, PartialEq)]
29pub struct Peak {
30    pub hz: f64,
31    pub db: f64,
32}
33
34#[derive(Clone, Debug, PartialEq)]
35pub struct Spectrum {
36    pub frame_size: usize,
37    pub frames: usize,
38    pub resolution_hz: f64,
39    pub rms: f64,
40    pub centroid_hz: f64,
41    pub rolloff85_hz: f64,
42    pub bands: Vec<Band>,
43    pub peaks: Vec<Peak>,
44}
45
46pub fn analyze(
47    planes: &[&[f64]],
48    sample_rate: f64,
49    max_peaks: usize,
50    frame_secs: Option<f64>,
51) -> Spectrum {
52    let (mags, bin_hz, frame_size, frames) = magnitudes_of(planes, sample_rate, frame_secs);
53    let power: Vec<f64> = mags.iter().map(|m| m * m).collect();
54    let total: f64 = power.iter().sum();
55
56    let centroid_hz = if total > 0.0 {
57        power
58            .iter()
59            .enumerate()
60            .map(|(k, p)| k as f64 * bin_hz * p)
61            .sum::<f64>()
62            / total
63    } else {
64        0.0
65    };
66
67    let mut running = 0.0;
68    let mut rolloff85_hz = 0.0;
69    for (k, p) in power.iter().enumerate() {
70        running += p;
71        if total > 0.0 && running >= 0.85 * total {
72            rolloff85_hz = k as f64 * bin_hz;
73            break;
74        }
75    }
76
77    let rms = crate::measure::envelope::level(planes);
78
79    Spectrum {
80        frame_size,
81        frames,
82        resolution_hz: bin_hz,
83        rms,
84        centroid_hz,
85        rolloff85_hz,
86        bands: bands(&power, bin_hz, sample_rate),
87        peaks: peaks(&mags, bin_hz, max_peaks),
88    }
89}
90
91/// Welch-averaged, Hann-windowed magnitudes normalized so a full-scale sine reads 1.0 at its
92/// own bin. Without `frame_secs` the size follows the WINDOW, so one node's two windows land
93/// on two bin grids and cannot be compared.
94pub fn magnitudes(
95    samples: &[f64],
96    sample_rate: f64,
97    frame_secs: Option<f64>,
98) -> (Vec<f64>, f64, usize, usize) {
99    magnitudes_of(&[samples], sample_rate, frame_secs)
100}
101
102pub fn magnitudes_of(
103    planes: &[&[f64]],
104    sample_rate: f64,
105    frame_secs: Option<f64>,
106) -> (Vec<f64>, f64, usize, usize) {
107    let len = crate::measure::held(planes);
108    let frame_size = match frame_secs {
109        Some(secs) => pinned_frame(secs, sample_rate).min(MAX_PINNED_FRAME),
110        None => len.next_power_of_two().clamp(MIN_FRAME, MAX_FRAME),
111    };
112    let hop = frame_size / 2;
113    let bins = frame_size / 2 + 1;
114    let window = crate::stft::hann_periodic(frame_size);
115
116    let mut acc = vec![0f64; bins];
117    let mut frames = 0usize;
118    let mut start = 0usize;
119    loop {
120        for samples in planes {
121            let mut re = vec![0f64; frame_size];
122            let mut im = vec![0f64; frame_size];
123            for n in 0..frame_size {
124                let s = samples.get(start + n).copied().unwrap_or(0.0);
125                re[n] = s * window[n];
126            }
127            fft(&mut re, &mut im);
128            for (k, a) in acc.iter_mut().enumerate() {
129                *a += re[k] * re[k] + im[k] * im[k];
130            }
131        }
132        frames += 1;
133        start += hop;
134        if start + frame_size > len {
135            break;
136        }
137    }
138
139    let scale = 4.0 / frame_size as f64;
140    let averaged = (frames * planes.len()) as f64;
141    let mags = acc.iter().map(|p| (p / averaged).sqrt() * scale).collect();
142    (mags, sample_rate / frame_size as f64, frame_size, frames)
143}
144
145pub fn db(amplitude: f64) -> f64 {
146    if amplitude <= 0.0 {
147        FLOOR_DB
148    } else {
149        (20.0 * amplitude.log10()).max(FLOOR_DB)
150    }
151}
152
153/// Third-octave from 20 Hz up, the resolution a listener's critical bands roughly follow.
154pub fn third_octave_edges(sample_rate: f64) -> Vec<(f64, f64)> {
155    let nyquist = sample_rate / 2.0;
156    let ratio = 2f64.powf(1.0 / 3.0);
157    let mut out = Vec::new();
158    let mut lo = 20.0;
159    while lo < nyquist {
160        let hi = (lo * ratio).min(nyquist);
161        out.push((lo, hi));
162        lo = hi;
163    }
164    out
165}
166
167fn bands(power: &[f64], bin_hz: f64, sample_rate: f64) -> Vec<Band> {
168    let mut out = Vec::new();
169    for (lo, hi) in third_octave_edges(sample_rate) {
170        let sum: f64 = power
171            .iter()
172            .enumerate()
173            .filter(|(k, _)| {
174                let hz = *k as f64 * bin_hz;
175                hz >= lo && hz < hi
176            })
177            .map(|(_, p)| p)
178            .sum();
179        out.push(Band {
180            lo_hz: lo,
181            hi_hz: hi,
182            db: db(sum.sqrt()),
183        });
184    }
185    out
186}
187
188/// Local maxima, parabolically interpolated so a peak between two bins reports its real
189/// frequency instead of the nearest bin's.
190pub fn peaks(mags: &[f64], bin_hz: f64, max_peaks: usize) -> Vec<Peak> {
191    let mut found: Vec<Peak> = Vec::new();
192    for k in 1..mags.len().saturating_sub(1) {
193        if mags[k] <= mags[k - 1] || mags[k] < mags[k + 1] || mags[k] <= 0.0 {
194            continue;
195        }
196        let (a, b, c) = (mags[k - 1], mags[k], mags[k + 1]);
197        let denom = a - 2.0 * b + c;
198        let delta = if denom.abs() > f64::EPSILON {
199            (0.5 * (a - c) / denom).clamp(-0.5, 0.5)
200        } else {
201            0.0
202        };
203        found.push(Peak {
204            hz: (k as f64 + delta) * bin_hz,
205            db: db(b - 0.25 * (a - c) * delta),
206        });
207    }
208    found.sort_by(|x, y| y.db.total_cmp(&x.db));
209    found.truncate(max_peaks);
210    found
211}
sva_samples/measure/spectrum.rs

sva_samples/measure/
spectrum.rs
