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
6const MAX_FRAME: usize = 8192;
7const 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#[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, Debug, PartialEq)]
29pub struct Spectrum {
30    pub frame_size: usize,
31    pub frames: usize,
32    pub resolution_hz: f64,
33    pub rms: f64,
34    pub centroid_hz: f64,
35    pub rolloff85_hz: f64,
36    pub bands: Vec<Band>,
37    pub peaks: Vec<Peak>,
38}
39
40pub fn analyze(
41    planes: &[&[f64]],
42    sample_rate: f64,
43    max_peaks: usize,
44    frame_secs: Option<f64>,
45) -> Spectrum {
46    let Framed {
47        mags,
48        bin_hz,
49        frame_size,
50        frames,
51        spanned,
52    } = 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: crate::measure::peaks::resolved(&mags, bin_hz, spanned, 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. A Hann spans samples under a frame.
94pub fn magnitudes(
95    samples: &[f64],
96    sample_rate: f64,
97    frame_secs: Option<f64>,
98) -> (Vec<f64>, f64, usize, usize) {
99    let f = magnitudes_of(&[samples], sample_rate, frame_secs);
100    (f.mags, f.bin_hz, f.frame_size, f.frames)
101}
102
103/// `spanned`: the samples the Hann covers.
104pub struct Framed {
105    pub mags: Vec<f64>,
106    pub bin_hz: f64,
107    pub frame_size: usize,
108    pub frames: usize,
109    pub spanned: usize,
110}
111
112pub fn magnitudes_of(planes: &[&[f64]], sample_rate: f64, frame_secs: Option<f64>) -> Framed {
113    let len = crate::measure::held(planes);
114    let frame_size = match frame_secs {
115        Some(secs) => pinned_frame(secs, sample_rate).min(MAX_PINNED_FRAME),
116        None => len.next_power_of_two().clamp(MIN_FRAME, MAX_FRAME),
117    };
118    let hop = frame_size / 2;
119    let bins = frame_size / 2 + 1;
120    let spanned = len.clamp(1, frame_size);
121    let window = crate::stft::hann_periodic(spanned);
122
123    let mut acc = vec![0f64; bins];
124    let mut frames = 0usize;
125    let mut start = 0usize;
126    loop {
127        for samples in planes {
128            let mut re = vec![0f64; frame_size];
129            let mut im = vec![0f64; frame_size];
130            for (n, w) in window.iter().enumerate() {
131                re[n] = samples.get(start + n).copied().unwrap_or(0.0) * w;
132            }
133            fft(&mut re, &mut im);
134            for (k, a) in acc.iter_mut().enumerate() {
135                *a += re[k] * re[k] + im[k] * im[k];
136            }
137        }
138        frames += 1;
139        start += hop;
140        if start + frame_size > len {
141            break;
142        }
143    }
144
145    let scale = 4.0 / spanned as f64;
146    let averaged = (frames * planes.len()) as f64;
147    Framed {
148        mags: acc.iter().map(|p| (p / averaged).sqrt() * scale).collect(),
149        bin_hz: sample_rate / frame_size as f64,
150        frame_size,
151        frames,
152        spanned,
153    }
154}
155
156/// Third-octave from 20 Hz up, the resolution a listener's critical bands roughly follow.
157pub fn third_octave_edges(sample_rate: f64) -> Vec<(f64, f64)> {
158    let nyquist = sample_rate / 2.0;
159    let ratio = 2f64.powf(1.0 / 3.0);
160    let mut out = Vec::new();
161    let mut lo = 20.0;
162    while lo < nyquist {
163        let hi = (lo * ratio).min(nyquist);
164        out.push((lo, hi));
165        lo = hi;
166    }
167    out
168}
169
170fn bands(power: &[f64], bin_hz: f64, sample_rate: f64) -> Vec<Band> {
171    let mut out = Vec::new();
172    for (lo, hi) in third_octave_edges(sample_rate) {
173        let sum: f64 = power
174            .iter()
175            .enumerate()
176            .filter(|(k, _)| {
177                let hz = *k as f64 * bin_hz;
178                hz >= lo && hz < hi
179            })
180            .map(|(_, p)| p)
181            .sum();
182        out.push(Band {
183            lo_hz: lo,
184            hi_hz: hi,
185            db: db(sum.sqrt()),
186        });
187    }
188    out
189}
sva_samples/measure/spectrum.rs

sva_samples/measure/
spectrum.rs
