//! Sample-rate conversion.
//!
//! Two implementations:
//!
//! - [`linear`] — straight linear interpolation. Cheap but introduces audible
//!   aliasing on downsamples; only suitable for non-critical paths or as a
//!   sanity baseline.
//! - [`SincResampler`] — windowed-sinc convolution with a Kaiser window.
//!   Higher quality at the cost of `O(N · taps)` work per output sample.
//!
//! Both treat any out-of-range input as silence (zero-pad). Mono `f32` only;
//! callers downmix multi-channel inputs upstream.

use alloc::vec;
use alloc::vec::Vec;

use libm::{sinf, sqrtf};

/// Linear resampler. Output length is `ceil(input.len() * to_sr / from_sr)`,
/// or `0` if the input is empty.
///
/// Pass-through when `from_sr == to_sr` (clones the input).
///
/// # Example
///
/// ```
/// use audiofp::dsp::resample::linear;
///
/// // Halve the rate: every other sample is kept.
/// let x = vec![0.0_f32, 1.0, 0.0, 1.0, 0.0, 1.0];
/// let y = linear(&x, 2, 1);
/// assert_eq!(y.len(), 3);
/// ```
#[must_use]
pub fn linear(input: &[f32], from_sr: u32, to_sr: u32) -> Vec<f32> {
    assert!(from_sr > 0 && to_sr > 0, "sample rates must be non-zero");
    if input.is_empty() {
        return Vec::new();
    }
    if from_sr == to_sr {
        return input.to_vec();
    }

    let n_in = input.len();
    let n_out = ((n_in as u64 * to_sr as u64).div_ceil(from_sr as u64)) as usize;
    let ratio = from_sr as f64 / to_sr as f64;

    let mut out = Vec::with_capacity(n_out);
    for n in 0..n_out {
        let pos = n as f64 * ratio;
        let i = pos.floor() as usize;
        let frac = (pos - pos.floor()) as f32;
        let a = input[i.min(n_in - 1)];
        let b = input[(i + 1).min(n_in - 1)];
        out.push(a * (1.0 - frac) + b * frac);
    }
    out
}

/// Quality knobs for [`SincResampler`].
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct SincQuality {
    /// Half the filter width in *input* samples. The filter spans
    /// `2 * half_taps` input samples around each output position.
    /// Typical: 16 (low-latency), 32 (good), 64 (excellent).
    pub half_taps: usize,
    /// Kaiser β parameter. Larger values trade transition-band sharpness
    /// for stopband attenuation. Typical: 8.6 (≈ -80 dB stopband), 12.0
    /// (≈ -120 dB).
    pub kaiser_beta: f32,
    /// Number of polyphase steps for the precomputed kernel table.
    /// Each step stores `2 * half_taps + 1` kernel coefficients for a
    /// different fractional offset. Default 256.
    pub polyphase_steps: u16,
}

impl Default for SincQuality {
    fn default() -> Self {
        Self {
            half_taps: 32,
            kaiser_beta: 8.6,
            polyphase_steps: 256,
        }
    }
}

/// Windowed-sinc Kaiser resampler.
///
/// Each output sample is computed by convolving the input with a sinc
/// kernel multiplied by a Kaiser window. The cutoff is automatically set
/// to `min(from_sr, to_sr) / 2` (in the input's sample-rate frame), which
/// suppresses aliasing on downsamples and limits image content on
/// upsamples.
///
/// # Example
///
/// ```
/// use audiofp::dsp::resample::SincResampler;
///
/// // 1.0 s of 1 kHz tone at 16 kHz → resample to 8 kHz.
/// let sr_in = 16_000u32;
/// let sr_out = 8_000u32;
/// let n = sr_in as usize;
/// let x: Vec<f32> = (0..n)
///     .map(|i| (2.0 * std::f32::consts::PI * 1000.0 * i as f32 / sr_in as f32).sin())
///     .collect();
///
/// let r = SincResampler::new(sr_in, sr_out);
/// let y = r.process(&x);
/// // Output length is roughly half of input.
/// assert!((y.len() as i64 - sr_out as i64).abs() < 16);
/// ```
pub struct SincResampler {
    from_sr: u32,
    to_sr: u32,
    quality: SincQuality,
    /// Precomputed reciprocal of the kernel's DC gain; multiplied per
    /// output sample to normalise constant signals to unity.
    inv_dc_gain: f32,
    /// Precomputed polyphase kernel table: `[steps][2*half_taps+1]`.
    /// Each row holds the sinc·Kaiser coefficients for one fractional
    /// offset, avoiding per-sample kernel evaluation in `process`.
    kernel_table: Vec<f32>,
}

impl SincResampler {
    /// Build a resampler with default quality (`half_taps=32`, `β=8.6`).
    #[must_use]
    pub fn new(from_sr: u32, to_sr: u32) -> Self {
        Self::with_quality(from_sr, to_sr, SincQuality::default())
    }

    /// Build a resampler with custom quality.
    ///
    /// # Panics
    ///
    /// Panics if `from_sr == 0`, `to_sr == 0`, or `quality.half_taps == 0`.
    #[must_use]
    pub fn with_quality(from_sr: u32, to_sr: u32, quality: SincQuality) -> Self {
        assert!(from_sr > 0 && to_sr > 0, "sample rates must be non-zero");
        assert!(quality.half_taps > 0, "half_taps must be > 0");
        assert!(quality.polyphase_steps > 0, "polyphase_steps must be > 0");

        let cutoff = from_sr.min(to_sr) as f32 / from_sr as f32 / 2.0;
        let inv_i0_beta = 1.0 / modified_bessel_i0(quality.kaiser_beta);
        let half = quality.half_taps;
        let steps = quality.polyphase_steps as usize;
        let taps = 2 * half + 1;

        // Precompute the polyphase kernel table.
        let mut kernel_table = vec![0.0_f32; steps * taps];
        let two_cutoff = 2.0 * cutoff;
        let half_f = half as f32;
        for s in 0..steps {
            let frac = s as f32 / steps as f32;
            let row = &mut kernel_table[s * taps..(s + 1) * taps];
            for (k, cell) in row.iter_mut().enumerate() {
                let x = (k as isize - half as isize) as f32 - frac;
                *cell = sinc(x * two_cutoff)
                    * kaiser_window_input(x, half_f, quality.kaiser_beta, inv_i0_beta)
                    * two_cutoff;
            }
        }

        // Compute DC gain from the centre step (frac = 0).
        let mut dc_gain = 0.0_f32;
        for &v in &kernel_table[..taps] {
            dc_gain += v;
        }
        if dc_gain.abs() < 1e-10 {
            dc_gain = 1.0;
        }
        let inv_dc_gain = 1.0 / dc_gain;

        Self {
            from_sr,
            to_sr,
            quality,
            inv_dc_gain,
            kernel_table,
        }
    }

    /// Borrow the quality knobs this resampler was built with.
    #[must_use]
    pub fn quality(&self) -> &SincQuality {
        &self.quality
    }

    /// Resample `input` into the configured target rate.
    ///
    /// Output length is `ceil(input.len() * to_sr / from_sr)`. Out-of-range
    /// neighbours are treated as zero (zero-pad boundary).
    #[must_use]
    pub fn process(&self, input: &[f32]) -> Vec<f32> {
        let mut out = Vec::new();
        self.process_into(input, &mut out);
        out
    }

    /// Resample `input` into a caller-owned buffer.
    ///
    /// Same semantics as [`process`](Self::process) but writes into `out`
    /// (cleared and resized to the correct length) instead of allocating
    /// a new buffer, so callers in a hot loop can reuse the allocation
    /// across chunks.
    pub fn process_into(&self, input: &[f32], out: &mut Vec<f32>) {
        if input.is_empty() {
            out.clear();
            return;
        }
        if self.from_sr == self.to_sr {
            out.clear();
            out.extend_from_slice(input);
            return;
        }

        let n_in = input.len();
        let n_out = ((n_in as u64 * self.to_sr as u64).div_ceil(self.from_sr as u64)) as usize;
        let ratio = self.from_sr as f64 / self.to_sr as f64;
        let half = self.quality.half_taps;
        let steps = self.quality.polyphase_steps as usize;
        let taps = 2 * half + 1;

        // Determine the output-sample range where the kernel is fully
        // inside the input buffer (no bounds check needed).
        let n_safe_start = ((half as f64) / ratio).ceil() as usize;
        let n_safe_end = (((n_in.saturating_sub(1 + half)) as f64) / ratio).floor() as usize + 1;
        let n_safe_end = n_safe_end.min(n_out);

        out.clear();
        out.reserve(n_out);

        // Left boundary: kernel may extend past the start of input.
        for n in 0..n_safe_start.min(n_out) {
            let pos = n as f64 * ratio;
            let i_centre = pos.floor() as isize;
            let frac = (pos - pos.floor()) as f32;
            let step = ((frac * steps as f32) as usize).min(steps - 1);
            let kernel = &self.kernel_table[step * taps..(step + 1) * taps];

            let mut acc = 0.0_f32;
            for (k, &coeff) in kernel.iter().enumerate() {
                let idx = i_centre + k as isize - half as isize;
                if idx < 0 || (idx as usize) >= n_in {
                    continue;
                }
                acc += input[idx as usize] * coeff;
            }
            out.push(acc * self.inv_dc_gain);
        }

        // Middle: kernel is fully inside input — no bounds check.
        for n in n_safe_start..n_safe_end {
            let pos = n as f64 * ratio;
            let i_centre = pos.floor() as usize;
            let frac = (pos - pos.floor()) as f32;
            let step = ((frac * steps as f32) as usize).min(steps - 1);
            let kernel = &self.kernel_table[step * taps..(step + 1) * taps];

            let base = i_centre.wrapping_sub(half);
            let acc: f32 = input[base..base + taps]
                .iter()
                .zip(kernel.iter())
                .map(|(&s, &k)| s * k)
                .sum();
            out.push(acc * self.inv_dc_gain);
        }

        // Right boundary: kernel may extend past the end of input.
        for n in n_safe_end..n_out {
            let pos = n as f64 * ratio;
            let i_centre = pos.floor() as isize;
            let frac = (pos - pos.floor()) as f32;
            let step = ((frac * steps as f32) as usize).min(steps - 1);
            let kernel = &self.kernel_table[step * taps..(step + 1) * taps];

            let mut acc = 0.0_f32;
            for (k, &coeff) in kernel.iter().enumerate() {
                let idx = i_centre + k as isize - half as isize;
                if idx < 0 || (idx as usize) >= n_in {
                    continue;
                }
                acc += input[idx as usize] * coeff;
            }
            out.push(acc * self.inv_dc_gain);
        }
    }
}

/// Normalised sinc: `sin(π·x) / (π·x)`, with `sinc(0) = 1`.
#[inline]
fn sinc(x: f32) -> f32 {
    if x.abs() < 1e-10 {
        return 1.0;
    }
    let pi_x = core::f32::consts::PI * x;
    sinf(pi_x) / pi_x
}

/// Kaiser window evaluated at `x` measured in input samples from the centre,
/// with width ±`half`. Returns 0 outside the window.
#[inline]
fn kaiser_window_input(x: f32, half: f32, beta: f32, inv_i0_beta: f32) -> f32 {
    let t = x / half;
    if t.abs() > 1.0 {
        return 0.0;
    }
    let arg = beta * sqrtf((1.0 - t * t).max(0.0));
    modified_bessel_i0(arg) * inv_i0_beta
}

/// Modified Bessel function of the first kind, order 0, via series
/// expansion. Converges fast for `|x| ≤ 30` (well past audio Kaiser β).
#[inline]
fn modified_bessel_i0(x: f32) -> f32 {
    let mut sum = 1.0_f64;
    let mut term = 1.0_f64;
    let half_x_sq = (x as f64 * x as f64) / 4.0;
    for k in 1..=30 {
        let kf = k as f64;
        term *= half_x_sq / (kf * kf);
        sum += term;
        if term / sum < 1e-12 {
            break;
        }
    }
    sum as f32
}

#[cfg(test)]
mod tests {
    use super::*;
    use alloc::vec;
    use approx::assert_relative_eq;
    use core::f32::consts::PI;

    #[test]
    fn linear_passthrough_when_rates_match() {
        let x: Vec<f32> = (0..32).map(|i| i as f32 * 0.1).collect();
        let y = linear(&x, 16_000, 16_000);
        assert_eq!(x, y);
    }

    #[test]
    fn linear_halves_length_on_2_to_1() {
        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let y = linear(&x, 2, 1);
        assert_eq!(y.len(), 3);
        // Linear at t=0,2,4 → exact samples 1, 3, 5.
        assert_relative_eq!(y[0], 1.0, max_relative = 1e-5);
        assert_relative_eq!(y[1], 3.0, max_relative = 1e-5);
        assert_relative_eq!(y[2], 5.0, max_relative = 1e-5);
    }

    #[test]
    fn linear_doubles_length_on_1_to_2() {
        let x = vec![0.0, 2.0, 4.0];
        let y = linear(&x, 1, 2);
        assert_eq!(y.len(), 6);
        // Halfway between 0 and 2 → 1.
        assert_relative_eq!(y[1], 1.0, max_relative = 1e-5);
    }

    #[test]
    fn sinc_resampler_passthrough_when_rates_match() {
        let r = SincResampler::new(16_000, 16_000);
        let x: Vec<f32> = (0..16).map(|i| i as f32 * 0.05).collect();
        assert_eq!(r.process(&x), x);
    }

    #[test]
    fn sinc_resampler_preserves_dc() {
        let r = SincResampler::new(44_100, 16_000);
        let x = vec![1.0_f32; 4096];
        let y = r.process(&x);
        // Skip the boundary region where zero-padding pulls the DC down.
        let mid_lo = y.len() / 4;
        let mid_hi = 3 * y.len() / 4;
        for &s in &y[mid_lo..mid_hi] {
            assert_relative_eq!(s, 1.0, epsilon = 1e-2);
        }
    }

    #[test]
    fn sinc_resampler_output_length_matches_ratio() {
        let r = SincResampler::new(48_000, 8_000);
        let x = vec![0.0_f32; 48_000];
        let y = r.process(&x);
        // 48k → 8k means 6× downsample.
        assert_eq!(y.len(), 8_000);
    }

    #[test]
    fn sinc_resampler_preserves_pure_tone() {
        // 1 kHz tone at 16 kHz, downsampled to 8 kHz. The tone is well below
        // 4 kHz Nyquist, so it should survive.
        let sr_in = 16_000u32;
        let sr_out = 8_000u32;
        let freq = 1000.0_f32;
        let n = (sr_in as usize) * 2;
        let x: Vec<f32> = (0..n)
            .map(|i| libm::sinf(2.0 * PI * freq * i as f32 / sr_in as f32))
            .collect();

        let r = SincResampler::new(sr_in, sr_out);
        let y = r.process(&x);

        // Compute peak frequency via crude DFT bin scan around 1000 Hz.
        // Use centre region to skip boundary roll-off.
        let lo = y.len() / 4;
        let hi = 3 * y.len() / 4;
        let seg = &y[lo..hi];
        let m = seg.len();

        // Probe a handful of frequencies and pick the strongest.
        let mut best = (0.0_f32, 0.0_f32);
        for f in (500..=1500).step_by(10) {
            let mut re = 0.0_f32;
            let mut im = 0.0_f32;
            for (k, &s) in seg.iter().enumerate() {
                let theta = 2.0 * PI * f as f32 * k as f32 / sr_out as f32;
                re += s * libm::cosf(theta);
                im -= s * libm::sinf(theta);
            }
            let mag = sqrtf(re * re + im * im) / m as f32;
            if mag > best.1 {
                best = (f as f32, mag);
            }
        }
        assert!(
            (best.0 - 1000.0).abs() <= 20.0,
            "peak at {} Hz, expected near 1000 Hz",
            best.0
        );
    }

    #[test]
    fn empty_input_produces_empty_output() {
        assert!(linear(&[], 44_100, 16_000).is_empty());
        let r = SincResampler::new(44_100, 16_000);
        assert!(r.process(&[]).is_empty());
    }

    #[test]
    fn linear_output_length_matches_ratio() {
        // Output length follows ceil(n_in * to / from).
        for (n_in, from, to, expect) in [
            (100, 8_000, 16_000, 200),
            (100, 16_000, 8_000, 50),
            (1000, 44_100, 22_050, 500),
            (3, 2, 1, 2), // ceil(3/2) = 2
        ] {
            let x = vec![1.0_f32; n_in];
            let y = linear(&x, from, to);
            assert_eq!(y.len(), expect, "n_in={n_in} {from}->{to}");
        }
    }

    #[test]
    fn sinc_resampler_with_higher_quality_has_lower_passband_droop() {
        // A 1 kHz tone at sr=16k stays well below the cutoff after a
        // 16k→8k downsample. Higher Kaiser β should preserve amplitude
        // closer to 1.0 at midband.
        let sr_in = 16_000u32;
        let sr_out = 8_000u32;
        let freq = 1000.0_f32;
        let n = sr_in as usize;
        let x: Vec<f32> = (0..n)
            .map(|i| libm::sinf(2.0 * PI * freq * i as f32 / sr_in as f32))
            .collect();

        let low_q = SincResampler::with_quality(
            sr_in,
            sr_out,
            SincQuality {
                half_taps: 8,
                kaiser_beta: 4.0,
                polyphase_steps: 256,
            },
        );
        let high_q = SincResampler::with_quality(
            sr_in,
            sr_out,
            SincQuality {
                half_taps: 64,
                kaiser_beta: 12.0,
                polyphase_steps: 256,
            },
        );

        let y_low = low_q.process(&x);
        let y_high = high_q.process(&x);

        // RMS ≈ 0.707 for a clean unit-amplitude sine. Use the centre
        // region to avoid boundary roll-off.
        let mid = |v: &[f32]| {
            let lo = v.len() / 4;
            let hi = 3 * v.len() / 4;
            let s: f32 = v[lo..hi].iter().map(|x| x * x).sum();
            sqrtf(s / (hi - lo) as f32)
        };
        let rms_low = mid(&y_low);
        let rms_high = mid(&y_high);

        // High quality should be closer to ideal RMS than low quality.
        let ideal = 1.0_f32 / sqrtf(2.0);
        assert!(
            (rms_high - ideal).abs() <= (rms_low - ideal).abs() + 1e-3,
            "low_q={rms_low} high_q={rms_high} ideal={ideal}",
        );
    }

    #[test]
    fn linear_doubling_then_halving_recovers_signal() {
        // Round-trip: 1× → 2× → 1× should recover the original up to
        // bilinear interpolation error.
        let n = 64;
        let x: Vec<f32> = (0..n)
            .map(|i| libm::sinf(2.0 * PI * i as f32 / 16.0))
            .collect();
        let up = linear(&x, 1, 2);
        let down = linear(&up, 2, 1);
        assert_eq!(down.len(), n);
        for (a, b) in x.iter().zip(down.iter()) {
            assert!((a - b).abs() < 0.05, "a={a} b={b}");
        }
    }

    #[test]
    fn modified_bessel_i0_matches_known_values() {
        // I0(0) = 1, I0(1) ≈ 1.2660658, I0(5) ≈ 27.2398718
        assert_relative_eq!(modified_bessel_i0(0.0), 1.0, max_relative = 1e-6);
        assert_relative_eq!(modified_bessel_i0(1.0), 1.266_065_8, max_relative = 1e-5);
        assert_relative_eq!(modified_bessel_i0(5.0), 27.239_872, max_relative = 1e-5);
    }

    #[test]
    fn quality_getter_returns_the_quality_passed_in() {
        let q = SincQuality {
            half_taps: 32,
            kaiser_beta: 8.6,
            polyphase_steps: 128,
        };
        let r = SincResampler::with_quality(8_000, 16_000, q);
        assert_eq!(*r.quality(), q);
    }

    // -----------------------------------------------------------------
    // Constructor panic coverage.
    //
    // `linear` and `SincResampler::with_quality` each have documented
    // panics on invalid input. Pin them with `should_panic` tests so
    // a refactor that loosens any assertion is caught.
    // -----------------------------------------------------------------

    #[test]
    #[should_panic(expected = "sample rates must be non-zero")]
    fn linear_panics_on_zero_from_sr() {
        let _ = linear(&[0.0; 16], 0, 8_000);
    }

    #[test]
    #[should_panic(expected = "sample rates must be non-zero")]
    fn linear_panics_on_zero_to_sr() {
        let _ = linear(&[0.0; 16], 8_000, 0);
    }

    #[test]
    #[should_panic(expected = "half_taps must be > 0")]
    fn sinc_resampler_panics_on_zero_half_taps() {
        let q = SincQuality {
            half_taps: 0,
            kaiser_beta: 8.6,
            polyphase_steps: 128,
        };
        let _ = SincResampler::with_quality(8_000, 16_000, q);
    }

    #[test]
    #[should_panic(expected = "polyphase_steps must be > 0")]
    fn sinc_resampler_panics_on_zero_polyphase_steps() {
        let q = SincQuality {
            half_taps: 32,
            kaiser_beta: 8.6,
            polyphase_steps: 0,
        };
        let _ = SincResampler::with_quality(8_000, 16_000, q);
    }
}
