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;
4pub use crate::measure::peaks::{Peak, db, peaks};
5
6pub const MAX_FRAME: usize = 8192;
7pub const MIN_FRAME: usize = 64;
8
9/// One FFT's memory, which is all that bounds a caller's own frame.
10pub const MAX_PINNED_FRAME: usize = 1 << 20;
11
12/// A caller checks this against [`MAX_PINNED_FRAME`]: a shrunk frame is not a pinned one.
13pub fn pinned_frame(secs: f64, sample_rate: f64) -> usize {
14    let asked = (secs * sample_rate).round().max(0.0);
15    if !asked.is_finite() || asked > MAX_PINNED_FRAME as f64 {
16        return MAX_PINNED_FRAME * 2;
17    }
18    (asked as usize).next_power_of_two().max(MIN_FRAME)
19}
20
21/// An unpinned spectrum's frame over `len` samples.
22pub fn fitted_frame(len: usize) -> usize {
23    len.next_power_of_two().clamp(MIN_FRAME, MAX_FRAME)
24}
25
26#[derive(Clone, Copy, Debug, PartialEq)]
27pub struct Band {
28    pub lo_hz: f64,
29    pub hi_hz: f64,
30    pub db: f64,
31}
32
33#[derive(Clone, Debug, PartialEq)]
34pub struct Spectrum {
35    pub frame_size: usize,
36    pub frames: usize,
37    pub resolution_hz: f64,
38    pub rms: f64,
39    pub centroid_hz: f64,
40    pub rolloff85_hz: f64,
41    pub bands: Vec<Band>,
42    pub peaks: Vec<Peak>,
43}
44
45pub fn analyze(
46    planes: &[&[f64]],
47    sample_rate: f64,
48    max_peaks: usize,
49    frame_secs: Option<f64>,
50) -> Spectrum {
51    let Framed {
52        mags,
53        bin_hz,
54        frame_size,
55        frames,
56        spanned,
57    } = magnitudes_of(planes, sample_rate, frame_secs);
58    let power = one_sided(&mags);
59    let total: f64 = power.iter().sum();
60
61    let centroid_hz = if total > 0.0 {
62        power
63            .iter()
64            .enumerate()
65            .map(|(k, p)| k as f64 * bin_hz * p)
66            .sum::<f64>()
67            / total
68    } else {
69        0.0
70    };
71
72    let mut running = 0.0;
73    let mut rolloff85_hz = 0.0;
74    for (k, p) in power.iter().enumerate() {
75        running += p;
76        if total > 0.0 && running >= 0.85 * total {
77            rolloff85_hz = k as f64 * bin_hz;
78            break;
79        }
80    }
81
82    let rms = crate::measure::envelope::level(planes);
83
84    Spectrum {
85        frame_size,
86        frames,
87        resolution_hz: bin_hz,
88        rms,
89        centroid_hz,
90        rolloff85_hz,
91        bands: bands(
92            &power,
93            1.5 * frame_size as f64 / spanned as f64,
94            bin_hz,
95            sample_rate,
96        ),
97        peaks: crate::measure::peaks::resolved(&mags, bin_hz, spanned, max_peaks),
98    }
99}
100
101/// 0 and Nyquist alone have no mirror folded in.
102fn one_sided(mags: &[f64]) -> Vec<f64> {
103    let last = mags.len() - 1;
104    mags.iter()
105        .enumerate()
106        .map(|(k, m)| match k == 0 || k == last {
107            true => m * m / 2.0,
108            false => m * m,
109        })
110        .collect()
111}
112
113/// A full-scale sine reads 1.0 at its own bin. Without `frame_secs` the size follows the
114/// range, so two ranges land on two bin grids.
115pub fn magnitudes(
116    samples: &[f64],
117    sample_rate: f64,
118    frame_secs: Option<f64>,
119) -> (Vec<f64>, f64, usize, usize) {
120    let f = magnitudes_of(&[samples], sample_rate, frame_secs);
121    (f.mags, f.bin_hz, f.frame_size, f.frames)
122}
123
124pub struct Framed {
125    pub mags: Vec<f64>,
126    pub bin_hz: f64,
127    pub frame_size: usize,
128    pub frames: usize,
129    /// The samples the Hann spans, fewer than `frame_size` where the range is.
130    pub spanned: usize,
131}
132
133pub fn magnitudes_of(planes: &[&[f64]], sample_rate: f64, frame_secs: Option<f64>) -> Framed {
134    let len = crate::measure::held(planes);
135    let frame_size = match frame_secs {
136        Some(secs) => pinned_frame(secs, sample_rate).min(MAX_PINNED_FRAME),
137        None => fitted_frame(len),
138    };
139    let bins = frame_size / 2 + 1;
140    let spanned = len.clamp(1, frame_size);
141    let window = crate::stft::hann_periodic(spanned);
142
143    let starts = starts(len, frame_size);
144    let mut acc = vec![0f64; bins];
145    for start in &starts {
146        for samples in planes {
147            let mut re = vec![0f64; frame_size];
148            let mut im = vec![0f64; frame_size];
149            for (n, w) in window.iter().enumerate() {
150                re[n] = samples.get(start + n).copied().unwrap_or(0.0) * w;
151            }
152            fft(&mut re, &mut im);
153            for (k, a) in acc.iter_mut().enumerate() {
154                *a += re[k] * re[k] + im[k] * im[k];
155            }
156        }
157    }
158
159    let scale = 4.0 / spanned as f64;
160    let averaged = (starts.len() * planes.len()) as f64;
161    Framed {
162        mags: acc.iter().map(|p| (p / averaged).sqrt() * scale).collect(),
163        bin_hz: sample_rate / frame_size as f64,
164        frame_size,
165        frames: starts.len(),
166        spanned,
167    }
168}
169
170/// Whole frames, first sample to last, evenly at most half a frame apart.
171pub(crate) fn starts(len: usize, frame: usize) -> Vec<usize> {
172    let Some(past) = len.checked_sub(frame).filter(|past| *past > 0) else {
173        return vec![0];
174    };
175    let gaps = past.div_ceil(frame / 2);
176    (0..=gaps)
177        .map(|i| (2 * i * past + gaps) / (2 * gaps))
178        .collect()
179}
180
181/// Third-octave from 20 Hz up.
182pub fn third_octave_edges(sample_rate: f64) -> Vec<(f64, f64)> {
183    let nyquist = sample_rate / 2.0;
184    let ratio = libm::cbrt(2.0);
185    let mut out = Vec::new();
186    let mut lo = 20.0;
187    while lo < nyquist {
188        let hi = (lo * ratio).min(nyquist);
189        out.push((lo, hi));
190        lo = hi;
191    }
192    out
193}
194
195/// Over `enbw`, the noise bandwidth in bins, a band reads as the one sine carrying its power.
196fn bands(power: &[f64], enbw: f64, bin_hz: f64, sample_rate: f64) -> Vec<Band> {
197    let mut out = Vec::new();
198    for (lo, hi) in third_octave_edges(sample_rate) {
199        let sum: f64 = power
200            .iter()
201            .enumerate()
202            .filter(|(k, _)| {
203                let hz = *k as f64 * bin_hz;
204                hz >= lo && hz < hi
205            })
206            .map(|(_, p)| p)
207            .sum();
208        out.push(Band {
209            lo_hz: lo,
210            hi_hz: hi,
211            db: db((sum / enbw).sqrt()),
212        });
213    }
214    out
215}
sva_samples/measure/spectrum.rs

sva_samples/measure/
spectrum.rs
