//! Continuous, block-boundary-independent 4x true-peak interpolation.

use std::collections::VecDeque;

const OVERSAMPLE: usize = 4;
const RADIUS: usize = 8;
const BUFFER_LEN: usize = RADIUS * 2 + 1;

#[derive(Clone, Debug)]
pub struct TruePeakDetector {
    samples: VecDeque<f32>,
}

impl Default for TruePeakDetector {
    fn default() -> Self {
        Self {
            samples: VecDeque::with_capacity(BUFFER_LEN + 1),
        }
    }
}

impl TruePeakDetector {
    pub fn process(&mut self, sample: f32) -> f32 {
        let mut peak = sample.abs();
        self.samples.push_back(sample);
        if self.samples.len() < BUFFER_LEN {
            return peak;
        }

        let mut contiguous = [0.0f32; BUFFER_LEN];
        for (target, source) in contiguous.iter_mut().zip(self.samples.iter()) {
            *target = *source;
        }
        for phase in 1..OVERSAMPLE {
            let position = RADIUS as f32 + phase as f32 / OVERSAMPLE as f32;
            peak = peak.max(interpolate(&contiguous, position).abs());
        }
        self.samples.pop_front();
        peak
    }
}

fn sinc(value: f32) -> f32 {
    if value.abs() < 1.0e-6 {
        1.0
    } else {
        let x = std::f32::consts::PI * value;
        x.sin() / x
    }
}

fn blackman(value: f32) -> f32 {
    let span = (RADIUS * 2) as f32;
    let phase = 2.0 * std::f32::consts::PI * (value + RADIUS as f32) / span;
    0.42 - 0.5 * phase.cos() + 0.08 * (2.0 * phase).cos()
}

fn interpolate(samples: &[f32], position: f32) -> f32 {
    let base = position.floor() as isize;
    let mut sum = 0.0;
    let mut norm = 0.0;
    for index in (base - RADIUS as isize + 1)..=(base + RADIUS as isize) {
        if !(0..samples.len() as isize).contains(&index) {
            continue;
        }
        let distance = position - index as f32;
        let weight = sinc(distance) * blackman(distance);
        sum += samples[index as usize] * weight;
        norm += weight;
    }
    if norm.abs() > 1.0e-6 {
        sum / norm
    } else {
        0.0
    }
}

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

    fn run_in_blocks(samples: &[f32], block_size: usize) -> f32 {
        let mut detector = TruePeakDetector::default();
        let mut peak = 0.0f32;
        for block in samples.chunks(block_size) {
            for &sample in block {
                peak = peak.max(detector.process(sample));
            }
        }
        for _ in 0..BUFFER_LEN {
            peak = peak.max(detector.process(0.0));
        }
        peak
    }

    #[test]
    fn result_is_independent_of_capture_block_boundaries() {
        let samples: Vec<f32> = (0..4_800)
            .map(|index| {
                (2.0 * std::f32::consts::PI * 11_025.0 * index as f32 / 48_000.0 + 0.31).sin() * 0.9
            })
            .collect();
        let reference = run_in_blocks(&samples, 1);
        for size in [64, 127, 128, 192, 511] {
            assert!((run_in_blocks(&samples, size) - reference).abs() < 1.0e-7);
        }
    }

    #[test]
    fn detects_an_intersample_peak_above_sample_peak() {
        let samples: Vec<f32> = (0..4_800)
            .map(|index| {
                (2.0 * std::f32::consts::PI * 11_025.0 * index as f32 / 48_000.0 + 0.31).sin() * 0.9
            })
            .collect();
        let sample_peak = samples
            .iter()
            .fold(0.0f32, |peak, value| peak.max(value.abs()));
        let detected = run_in_blocks(&samples, 128);
        assert!(detected > sample_peak + 0.001);
        assert!(
            detected <= 0.91,
            "unexpected interpolation overshoot: {detected}"
        );
    }
}