harmos_signal/spectrum/

harmonic.rs

//! Peaks, the fundamental beneath a taper's DC leakage, and the distortion the
//! rest of the spectrum accounts for.

/// The bins that are local maxima, strongest first.
///
/// A bin is a peak when it stands above both neighbours; the first and last
/// bins qualify on their one neighbour alone. Equal peaks keep bin order, so
/// the answer is stable.
///
/// Indices rather than values, like every reducer in this crate: the caller
/// reads its own columns at the bins it is handed.
#[must_use]
pub fn peaks(magnitude: &[f64]) -> Vec<usize> {
    let mut found: Vec<usize> = Vec::new();

    if magnitude.len() >= 2 {
        if magnitude[0] > magnitude[1] {
            found.push(0);
        }
        found.extend(
            magnitude
                .windows(3)
                .enumerate()
                .filter(|(_, window)| window[1] > window[0] && window[1] > window[2])
                .map(|(index, _)| index + 1),
        );
        let last = magnitude.len() - 1;
        if magnitude[last] > magnitude[last - 1] {
            found.push(last);
        }
    }

    found.sort_by(|left, right| magnitude[*right].total_cmp(&magnitude[*left]));
    found
}

/// The peaks at or above `fraction` of the strongest one, in bin order.
///
/// `fraction` is a share of the largest magnitude in the spectrum, so `0.1`
/// keeps everything within 20 dB of the peak.
#[must_use]
pub fn harmonics(magnitude: &[f64], fraction: f64) -> Vec<usize> {
    let strongest = magnitude.iter().copied().fold(0.0_f64, f64::max);
    let floor = fraction * strongest;

    let mut found: Vec<usize> = peaks(magnitude)
        .into_iter()
        .filter(|bin| magnitude[*bin] >= floor)
        .collect();

    found.sort_unstable();
    found
}

/// The bin of the fundamental: the strongest peak that is not DC.
///
/// Peaks, not raw magnitudes, and that is the whole point. A taper leaks DC
/// into the bins beside it hard enough that the leakage can out-measure the
/// signal — but leakage decays away from DC and so is never a local maximum,
/// while the signal is. Anything within half a bin of DC is skipped outright.
#[must_use]
pub fn fundamental(frequency: &[f64], magnitude: &[f64]) -> Option<usize> {
    let guard = step_of(frequency) * 0.5;

    peaks(magnitude)
        .into_iter()
        .find(|bin| frequency.get(*bin).is_some_and(|at| *at > guard))
}

/// Total harmonic distortion as a percentage of the fundamental.
///
/// Every bin that is neither DC nor the fundamental counts as distortion:
/// `100 · sqrt(Σ m²) / m_fundamental`. A spectrum with no fundamental, or one
/// whose fundamental is zero, has no answer.
#[must_use]
pub fn thd(frequency: &[f64], magnitude: &[f64]) -> Option<f64> {
    let bin = fundamental(frequency, magnitude)?;
    let strongest = magnitude[bin];
    if strongest < f64::EPSILON {
        return None;
    }

    let tolerance = step_of(frequency) * 0.5;
    let distortion: f64 = magnitude
        .iter()
        .zip(frequency)
        .filter(|(_, at)| **at > tolerance && (**at - frequency[bin]).abs() > tolerance)
        .map(|(value, _)| value * value)
        .sum();

    Some(distortion.sqrt() / strongest * 100.0)
}

/// The bin step an axis implies, or one hertz when it implies nothing.
fn step_of(frequency: &[f64]) -> f64 {
    match frequency {
        [first, second, ..] => second - first,
        _ => 1.0,
    }
}

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

    /// A spectrum of `length` empty bins with the named ones filled.
    fn spectrum(length: usize, filled: &[(usize, f64)]) -> Vec<f64> {
        let mut magnitude = vec![0.0; length];
        for (bin, value) in filled {
            magnitude[*bin] = *value;
        }
        magnitude
    }

    #[test]
    fn peaks_come_back_strongest_first() {
        let magnitude = spectrum(101, &[(10, 5.0), (30, 10.0), (50, 7.0), (70, 2.0)]);

        assert_eq!(peaks(&magnitude), vec![30, 50, 10, 70]);
    }

    #[test]
    fn a_run_with_no_local_maximum_has_no_peaks() {
        assert!(peaks(&[1.0, 1.0, 1.0, 1.0]).is_empty());
        assert!(peaks(&[]).is_empty());
        assert!(peaks(&[1.0]).is_empty());
    }

    #[test]
    fn the_edges_qualify_on_their_one_neighbour() {
        assert_eq!(peaks(&[5.0, 1.0, 0.0]), vec![0]);
        assert_eq!(peaks(&[0.0, 1.0, 5.0]), vec![2]);
    }

    #[test]
    fn harmonics_keep_what_clears_the_fraction_in_bin_order() {
        let magnitude = spectrum(201, &[(0, 5.0), (50, 10.0), (150, 3.0)]);

        assert_eq!(harmonics(&magnitude, 0.1), vec![0, 50, 150]);
        assert_eq!(harmonics(&magnitude, 0.4), vec![0, 50]);
        assert_eq!(harmonics(&magnitude, 0.9), vec![50]);
    }

    #[test]
    fn the_fundamental_is_the_strongest_peak_off_dc() {
        let frequency = vec![0.0, 50.0, 100.0, 150.0];
        let magnitude = vec![5.0, 20.0, 10.0, 5.0];

        assert_eq!(fundamental(&frequency, &magnitude), Some(1));
    }

    #[test]
    fn dc_leakage_as_strong_as_dc_is_still_not_the_fundamental() {
        // A Hann taper can leak DC into bin 1 at DC's own apparent amplitude.
        // Bin 1 is a local maximum there, so only the DC guard rejects it.
        let frequency: Vec<f64> = (0..5001).map(|bin| f64::from(bin) * 1000.0).collect();
        let magnitude = spectrum(5001, &[(0, 75.0), (1, 75.0), (50, 5.0)]);

        assert_eq!(fundamental(&frequency, &magnitude), Some(50));
    }

    #[test]
    fn dc_leakage_that_merely_decays_is_never_a_peak() {
        let frequency: Vec<f64> = (0..100).map(|bin| f64::from(bin) * 100.0).collect();
        let magnitude = spectrum(100, &[(0, 50.0), (1, 30.0), (2, 10.0), (40, 8.0)]);

        assert_eq!(fundamental(&frequency, &magnitude), Some(40));
    }

    #[test]
    fn a_pure_tone_distorts_by_nothing() {
        let frequency: Vec<f64> = (0..100).map(|bin| f64::from(bin) * 10.0).collect();
        let magnitude = spectrum(100, &[(10, 4.0)]);

        assert!(thd(&frequency, &magnitude).is_some_and(|thd| thd.abs() < 1e-12));
    }

    #[test]
    fn one_harmonic_at_half_the_fundamental_is_fifty_percent() {
        let frequency: Vec<f64> = (0..100).map(|bin| f64::from(bin) * 10.0).collect();
        let magnitude = spectrum(100, &[(10, 4.0), (20, 2.0)]);

        assert!(thd(&frequency, &magnitude).is_some_and(|thd| (thd - 50.0).abs() < 1e-9));
    }

    #[test]
    fn a_spectrum_with_no_peak_has_no_distortion_to_report() {
        assert!(thd(&[0.0, 10.0, 20.0], &[0.0, 0.0, 0.0]).is_none());
    }
}