harmos_signal/spectrum/

indicator.rs

//! What one band of a spectrum measures: how much energy is in it, where that
//! energy sits, and how evenly it is spread.
//!
//! Every reading here takes a frequency axis and a *linear* magnitude column.
//! A caller holding decibels converts first with [`Decibel::to_linear`]; a
//! reading that silently converted would be guessing which reference it was
//! handed.
//!
//! [`Decibel::to_linear`]: super::Decibel::to_linear

use std::ops::Range;

/// The energy of one band, by Parseval: `Σ |X|² · df`.
///
/// The band is inclusive at both ends. A band holding no bin has no energy.
#[must_use]
pub fn band_energy(frequency: &[f64], magnitude: &[f64], low: f64, high: f64) -> f64 {
    let band = band(frequency, low, high);
    if band.is_empty() {
        return 0.0;
    }

    let step = step_of(frequency, &band);
    magnitude[band]
        .iter()
        .map(|value| value * value * step)
        .sum()
}

/// The RMS amplitude of one band: `sqrt(energy / bandwidth)`.
///
/// The effective amplitude every component in the band adds up to, in the same
/// unit the magnitudes are in.
#[must_use]
pub fn band_rms(frequency: &[f64], magnitude: &[f64], low: f64, high: f64) -> f64 {
    let bandwidth = high - low;
    if bandwidth <= 0.0 {
        return 0.0;
    }

    (band_energy(frequency, magnitude, low, high) / bandwidth).sqrt()
}

/// The spectral centroid of one band: its centre of energy, in Hz.
///
/// `Σ f · |X|² / Σ |X|²` — where the band would balance if energy were weight.
/// A band with no bins, or no energy in them, has no centroid.
#[must_use]
pub fn centroid(frequency: &[f64], magnitude: &[f64], low: f64, high: f64) -> Option<f64> {
    let band = band(frequency, low, high);
    if band.is_empty() {
        return None;
    }

    let (weighted, total) = frequency[band.clone()].iter().zip(&magnitude[band]).fold(
        (0.0, 0.0),
        |(weighted, total), (at, value)| {
            let power = value * value;
            (weighted + at * power, total + power)
        },
    );

    (total > 0.0).then_some(weighted / total)
}

/// The peak-to-average power ratio of one band, as a linear ratio.
///
/// High when a few components dominate, near one when the band is flat. In
/// decibels it is `10 · log10` of this.
#[must_use]
pub fn papr(frequency: &[f64], magnitude: &[f64], low: f64, high: f64) -> Option<f64> {
    let band = &magnitude[band(frequency, low, high)];
    if band.is_empty() {
        return None;
    }

    let peak = band.iter().copied().fold(0.0_f64, f64::max);
    let mean = band.iter().sum::<f64>() / band.len() as f64;

    (mean > 0.0).then(|| (peak * peak) / (mean * mean))
}

/// The frequency beneath which `fraction` of the band's energy lies, in Hz.
///
/// The roll-off point: `0.95` answers where ninety-five percent of the energy
/// has accumulated, walking upwards from the bottom of the band.
#[must_use]
pub fn rolloff(
    frequency: &[f64],
    magnitude: &[f64],
    low: f64,
    high: f64,
    fraction: f64,
) -> Option<f64> {
    let band = band(frequency, low, high);
    if band.is_empty() {
        return None;
    }

    let step = step_of(frequency, &band);
    let total: f64 = magnitude[band.clone()]
        .iter()
        .map(|value| value * value * step)
        .sum();
    if total <= 0.0 {
        return None;
    }

    let target = total * fraction.clamp(0.0, 1.0);
    let mut accumulated = 0.0;

    for (at, value) in frequency[band.clone()].iter().zip(&magnitude[band.clone()]) {
        accumulated += value * value * step;
        if accumulated >= target {
            return Some(*at);
        }
    }

    frequency.get(band.end - 1).copied()
}

/// The narrowest span of the band holding `fraction` of its energy.
///
/// Answers `(lower, upper, bandwidth)` in Hz, where the bandwidth includes the
/// width of the bins at both ends. Unlike [`rolloff`], which always starts at
/// the bottom, this slides a window: a signal sitting high in the band is
/// measured where it is.
#[must_use]
pub fn occupied_bandwidth(
    frequency: &[f64],
    magnitude: &[f64],
    low: f64,
    high: f64,
    fraction: f64,
) -> Option<(f64, f64, f64)> {
    let band = band(frequency, low, high);
    if band.is_empty() {
        return None;
    }

    let step = step_of(frequency, &band);
    let energies: Vec<f64> = magnitude[band.clone()]
        .iter()
        .map(|value| value * value * step)
        .collect();
    let total: f64 = energies.iter().sum();
    if total <= 0.0 {
        return None;
    }

    let target = total * fraction.clamp(0.0, 1.0);
    // Prefix sums, so the energy of any window is one subtraction.
    let mut prefix = vec![0.0; energies.len() + 1];
    for (index, energy) in energies.iter().enumerate() {
        prefix[index + 1] = prefix[index] + energy;
    }

    let (mut lower, mut upper, mut width) = (0, energies.len() - 1, energies.len());
    let mut left = 0;

    for right in 0..energies.len() {
        // Shrink from the left while the window still clears the target, so
        // what is checked below is the narrowest window ending here.
        while left < right && prefix[right + 1] - prefix[left + 1] >= target {
            left += 1;
        }

        if prefix[right + 1] - prefix[left] >= target && right - left + 1 < width {
            width = right - left + 1;
            lower = left;
            upper = right;
        }
    }

    let lower = frequency[band.start + lower];
    let upper = frequency[band.start + upper];

    Some((lower, upper, upper - lower + step))
}

/// The bins of a sorted frequency axis inside `low..=high`.
fn band(frequency: &[f64], low: f64, high: f64) -> Range<usize> {
    let start = frequency.partition_point(|at| *at < low);
    let end = frequency.partition_point(|at| *at <= high);

    start..end.max(start)
}

/// The bin width to integrate a band with.
///
/// A uniform axis states it outright; a band of one bin borrows a neighbour's
/// spacing; anything else averages the band's own.
fn step_of(frequency: &[f64], band: &Range<usize>) -> f64 {
    match band.len() {
        0 => 1.0,
        1 => match (band.start.checked_sub(1), frequency.get(band.end)) {
            (Some(before), _) => frequency[band.start] - frequency[before],
            (None, Some(after)) => after - frequency[band.start],
            (None, None) => 1.0,
        },
        length => (frequency[band.end - 1] - frequency[band.start]) / (length - 1) as f64,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// A uniform axis of `count` bins `step` Hz apart, starting at `start`.
    fn axis(start: f64, step: f64, count: usize) -> Vec<f64> {
        (0..count).map(|bin| start + bin as f64 * step).collect()
    }

    #[test]
    fn one_bin_of_energy_is_its_square_times_the_bin_width() {
        let frequency = axis(0.0, 10.0, 3);

        assert!((band_energy(&frequency, &[5.0, 0.0, 0.0], 0.0, 5.0) - 250.0).abs() < 1e-9);
    }

    #[test]
    fn band_energy_adds_the_bins_it_covers() {
        let frequency = axis(0.0, 10.0, 3);

        assert!((band_energy(&frequency, &[1.0, 2.0, 3.0], 0.0, 30.0) - 140.0).abs() < 1e-9);
    }

    #[test]
    fn a_band_with_no_bins_has_no_energy() {
        let frequency = axis(100.0, 10.0, 3);

        assert_eq!(band_energy(&frequency, &[1.0, 2.0, 3.0], 0.0, 50.0), 0.0);
    }

    #[test]
    fn a_flat_band_rms_is_the_flat_amplitude() {
        let frequency = axis(0.0, 10.0, 5);

        assert!((band_rms(&frequency, &[5.0; 5], 0.0, 50.0) - 5.0).abs() < 1e-9);
        assert_eq!(band_rms(&frequency, &[5.0; 5], 50.0, 50.0), 0.0);
    }

    #[test]
    fn a_symmetric_band_balances_at_its_middle() {
        let frequency = vec![50.0, 100.0, 150.0];

        let centre = centroid(&frequency, &[1.0, 2.0, 1.0], 0.0, 200.0);
        assert!(centre.is_some_and(|centre| (centre - 100.0).abs() < 1e-9));
    }

    #[test]
    fn energy_pulls_the_centroid_towards_itself() {
        let frequency = vec![100.0, 200.0];

        // (100·1 + 200·9) / 10 = 190.
        let centre = centroid(&frequency, &[1.0, 3.0], 0.0, 300.0);
        assert!(centre.is_some_and(|centre| (centre - 190.0).abs() < 1e-9));
    }

    #[test]
    fn a_band_with_no_energy_has_no_centroid() {
        assert!(centroid(&[100.0, 200.0], &[0.0, 0.0], 0.0, 300.0).is_none());
        assert!(centroid(&[100.0], &[1.0], 200.0, 300.0).is_none());
    }

    #[test]
    fn a_flat_band_peaks_at_its_own_average() {
        let frequency = vec![50.0, 100.0, 150.0];

        assert!(
            papr(&frequency, &[5.0; 3], 0.0, 200.0).is_some_and(|papr| (papr - 1.0).abs() < 1e-9)
        );
        assert!(papr(&frequency, &[0.1, 10.0, 0.1], 0.0, 200.0).is_some_and(|papr| papr > 5.0));
        assert!(papr(&frequency, &[5.0; 3], 300.0, 400.0).is_none());
    }

    #[test]
    fn rolloff_walks_up_until_the_share_is_met() {
        let frequency = axis(10.0, 10.0, 5);

        let half = rolloff(&frequency, &[1.0; 5], 0.0, 60.0, 0.5);
        assert!(half.is_some_and(|at| (20.0..=30.0).contains(&at)));

        let all = rolloff(&frequency, &[1.0; 5], 0.0, 60.0, 1.0);
        assert!(all.is_some_and(|at| (at - 50.0).abs() < 1e-9));
    }

    #[test]
    fn energy_concentrated_low_rolls_off_immediately() {
        let frequency = axis(10.0, 10.0, 5);

        let at = rolloff(&frequency, &[10.0, 1.0, 1.0, 1.0, 1.0], 0.0, 60.0, 0.95);
        assert!(at.is_some_and(|at| (at - 10.0).abs() < 1e-9));
    }

    #[test]
    fn occupied_bandwidth_finds_the_peak_rather_than_starting_low() {
        let frequency = axis(10.0, 10.0, 5);

        let (lower, upper, width) =
            occupied_bandwidth(&frequency, &[0.1, 0.1, 10.0, 0.1, 0.1], 0.0, 60.0, 0.95)
                .expect("a band with energy has an occupied width");

        assert!((lower - 30.0).abs() < 1e-9);
        assert!((upper - 30.0).abs() < 1e-9);
        assert!((width - 10.0).abs() < 1e-9);
    }

    #[test]
    fn spread_energy_occupies_the_whole_band() {
        let frequency = axis(10.0, 10.0, 5);

        let (_, _, width) = occupied_bandwidth(&frequency, &[1.0; 5], 0.0, 60.0, 0.99)
            .expect("a band with energy has an occupied width");

        assert!(width >= 40.0);
    }

    #[test]
    fn a_silent_band_occupies_nothing() {
        assert!(occupied_bandwidth(&[10.0, 20.0], &[0.0, 0.0], 0.0, 30.0, 0.99).is_none());
        assert!(occupied_bandwidth(&[10.0, 20.0], &[1.0, 1.0], 100.0, 200.0, 0.99).is_none());
    }

    #[test]
    fn a_single_bin_band_borrows_a_neighbours_spacing() {
        let frequency = axis(0.0, 10.0, 3);

        // Only bin 1 is inside, and its width is the axis step either way.
        assert!((band_energy(&frequency, &[1.0, 2.0, 3.0], 6.0, 14.0) - 40.0).abs() < 1e-9);
        assert!((band_energy(&frequency, &[1.0, 2.0, 3.0], 0.0, 4.0) - 10.0).abs() < 1e-9);
    }
}