harmos_signal/spectrum/

taper.rs

//! The two shaping vocabularies: what a taper does to samples before the
//! transform, and what a band does to bins before the reconstruction.

use std::f64::consts::PI;

/// A window function applied to samples before the forward transform.
///
/// Every taper trades main-lobe width against sidelobe height. A rectangular
/// taper is no taper at all: it keeps amplitude exact for a tone that lands on
/// a bin and leaks badly for one that does not.
///
/// The forward transform divides by this taper's coherent gain, so a tone's
/// recovered amplitude is the tone's amplitude whichever taper measured it.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Taper {
    /// No shaping: exact on-bin amplitude, worst leakage off it.
    Rectangular,
    /// Hamming — a lower first sidelobe than Hann at a higher noise floor.
    Hamming,
    /// Hann, once spelled Hanning — the general-purpose choice, zero at both ends.
    Hann,
    /// Blackman — the lowest sidelobes here, and the widest main lobe.
    Blackman,
}

impl Taper {
    /// This taper's coefficient for sample `index` of a `length`-sample segment.
    #[must_use]
    pub fn coefficient(&self, index: usize, length: usize) -> f64 {
        let length = length as f64;
        let index = index as f64;

        match self {
            Self::Rectangular => 1.0,
            Self::Hamming => 0.54 - 0.46 * (2.0 * PI * index / length).cos(),
            Self::Hann => 0.5 * (1.0 - (2.0 * PI * index / length).cos()),
            Self::Blackman => {
                0.42 - 0.5 * (2.0 * PI * index / length).cos()
                    + 0.08 * (4.0 * PI * index / length).cos()
            }
        }
    }

    /// The coefficients of a whole `length`-sample segment.
    #[must_use]
    pub fn coefficients(&self, length: usize) -> Vec<f64> {
        (0..length)
            .map(|index| self.coefficient(index, length))
            .collect()
    }
}

/// A frequency-domain window: which bins survive a reconstruction, and how
/// sharply the surviving band ends.
///
/// Every shape is zero outside `start..=end`. They differ in what happens
/// inside it, which is the difference between a brick-wall filter and one whose
/// edges do not ring.
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Band {
    /// A brick wall: full pass inside the band, nothing outside it.
    Rectangular {
        /// The lowest frequency the band keeps, in Hz.
        start: f64,
        /// The highest frequency the band keeps, in Hz.
        end: f64,
    },
    /// Raised-cosine edges of a stated width, flat in between.
    RaisedCosine {
        /// The lowest frequency the band keeps, in Hz.
        start: f64,
        /// The highest frequency the band keeps, in Hz.
        end: f64,
        /// The width of each cosine edge, in Hz.
        transition: f64,
    },
    /// A Kaiser shape, whose `beta` trades edge width against sidelobe height.
    Kaiser {
        /// The lowest frequency the band keeps, in Hz.
        start: f64,
        /// The highest frequency the band keeps, in Hz.
        end: f64,
        /// The shape parameter: larger is smoother and wider-edged.
        beta: f64,
    },
    /// A Gaussian shape of a stated standard deviation, centred in the band.
    Gaussian {
        /// The lowest frequency the band keeps, in Hz.
        start: f64,
        /// The highest frequency the band keeps, in Hz.
        end: f64,
        /// The standard deviation, in Hz.
        sigma: f64,
    },
}

impl Band {
    /// What fraction of one bin at `frequency` this band lets through, `0.0..=1.0`.
    #[must_use]
    pub fn coefficient(&self, frequency: f64) -> f64 {
        let (start, end) = self.bounds();
        if frequency < start || frequency > end {
            return 0.0;
        }

        match self {
            Self::Rectangular { .. } => 1.0,
            Self::RaisedCosine { transition, .. } => {
                if frequency < start + transition {
                    0.5 * (1.0 - (PI * (frequency - start) / transition).cos())
                } else if frequency > end - transition {
                    0.5 * (1.0 - (PI * (end - frequency) / transition).cos())
                } else {
                    1.0
                }
            }
            Self::Kaiser { beta, .. } => {
                let half = (end - start) / 2.0;
                let offset = (frequency - (start + end) / 2.0) / half;
                bessel_i0(beta * (1.0 - offset * offset).sqrt()) / bessel_i0(*beta)
            }
            Self::Gaussian { sigma, .. } => {
                let offset = (frequency - (start + end) / 2.0) / sigma;
                (-0.5 * offset * offset).exp()
            }
        }
    }

    /// The band's frequency bounds, in Hz.
    #[must_use]
    pub fn bounds(&self) -> (f64, f64) {
        match self {
            Self::Rectangular { start, end }
            | Self::RaisedCosine { start, end, .. }
            | Self::Kaiser { start, end, .. }
            | Self::Gaussian { start, end, .. } => (*start, *end),
        }
    }
}

/// The modified Bessel function of the first kind, order zero.
///
/// The Abramowitz and Stegun polynomial approximation, which is what a Kaiser
/// shape is normally computed from and is accurate to about seven digits.
fn bessel_i0(argument: f64) -> f64 {
    let magnitude = argument.abs();

    if magnitude < 3.75 {
        let square = (argument / 3.75).powi(2);
        1.0 + square
            * (3.515_622_9
                + square
                    * (3.089_942_4
                        + square
                            * (1.206_749_2
                                + square
                                    * (0.265_973_2
                                        + square * (0.036_076_8 + square * 0.004_581_3)))))
    } else {
        let ratio = 3.75 / magnitude;
        (magnitude.exp() / magnitude.sqrt())
            * (0.398_942_28
                + ratio
                    * (0.013_285_92
                        + ratio
                            * (0.002_253_19
                                + ratio
                                    * (-0.001_575_65
                                        + ratio
                                            * (0.009_162_81
                                                + ratio
                                                    * (-0.020_577_06
                                                        + ratio
                                                            * (0.026_355_37
                                                                + ratio
                                                                    * (-0.016_476_33
                                                                        + ratio
                                                                            * 0.003_923_77))))))))
    }
}

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

    #[test]
    fn a_rectangular_taper_is_no_taper() {
        for index in 0..64 {
            assert_eq!(Taper::Rectangular.coefficient(index, 64), 1.0);
        }
    }

    #[test]
    fn hann_starts_at_zero_and_hamming_at_its_pedestal() {
        assert!(Taper::Hann.coefficient(0, 128) < 1e-12);
        assert!((Taper::Hamming.coefficient(0, 128) - 0.08).abs() < 1e-12);
    }

    #[test]
    fn every_taper_is_symmetric_about_its_centre() {
        let length = 100;
        for taper in [Taper::Hamming, Taper::Hann, Taper::Blackman] {
            let left = taper.coefficient(length / 2 - 10, length);
            let right = taper.coefficient(length / 2 + 10, length);
            assert!((left - right).abs() < 1e-12, "{taper:?} is not symmetric");
        }
    }

    #[test]
    fn coefficients_are_the_whole_segment() {
        let coefficients = Taper::Blackman.coefficients(32);

        assert_eq!(coefficients.len(), 32);
        assert_eq!(coefficients[7], Taper::Blackman.coefficient(7, 32));
    }

    #[test]
    fn a_rectangular_band_is_a_brick_wall() {
        let band = Band::Rectangular {
            start: 100.0,
            end: 200.0,
        };

        assert_eq!(band.coefficient(99.0), 0.0);
        assert_eq!(band.coefficient(100.0), 1.0);
        assert_eq!(band.coefficient(150.0), 1.0);
        assert_eq!(band.coefficient(200.0), 1.0);
        assert_eq!(band.coefficient(201.0), 0.0);
    }

    #[test]
    fn a_raised_cosine_band_rises_and_falls_inside_its_edges() {
        let band = Band::RaisedCosine {
            start: 100.0,
            end: 200.0,
            transition: 20.0,
        };

        assert_eq!(band.coefficient(50.0), 0.0);
        assert_eq!(band.coefficient(150.0), 1.0);
        assert!((0.0..1.0).contains(&band.coefficient(110.0)));
        assert!((0.0..1.0).contains(&band.coefficient(190.0)));
    }

    #[test]
    fn a_kaiser_band_peaks_at_its_centre() {
        let band = Band::Kaiser {
            start: 0.0,
            end: 100.0,
            beta: 6.0,
        };

        assert!((band.coefficient(50.0) - 1.0).abs() < 1e-9);
        assert!(band.coefficient(10.0) < band.coefficient(40.0));
        assert_eq!(band.coefficient(101.0), 0.0);
    }

    #[test]
    fn a_gaussian_band_falls_off_by_its_sigma() {
        let band = Band::Gaussian {
            start: 0.0,
            end: 100.0,
            sigma: 25.0,
        };

        assert!((band.coefficient(50.0) - 1.0).abs() < 1e-12);
        // One sigma from the centre is exp(-1/2) of the peak.
        assert!((band.coefficient(75.0) - (-0.5f64).exp()).abs() < 1e-12);
    }

    #[test]
    fn every_band_reports_the_bounds_it_was_built_with() {
        assert_eq!(
            Band::Rectangular {
                start: 1.0,
                end: 2.0
            }
            .bounds(),
            (1.0, 2.0)
        );
        assert_eq!(
            Band::RaisedCosine {
                start: 1.0,
                end: 2.0,
                transition: 0.1
            }
            .bounds(),
            (1.0, 2.0)
        );
        assert_eq!(
            Band::Kaiser {
                start: 1.0,
                end: 2.0,
                beta: 6.0
            }
            .bounds(),
            (1.0, 2.0)
        );
        assert_eq!(
            Band::Gaussian {
                start: 1.0,
                end: 2.0,
                sigma: 0.5
            }
            .bounds(),
            (1.0, 2.0)
        );
    }
}