Function bodies 463 total

    fn gap_in_middle() {
        let contour = vec![
            frame(0.0, Some(100.0)),
            frame(0.01, Some(101.0)),
            frame(0.02, None), // gap
            frame(0.03, None),
            frame(0.04, Some(99.0)),
            frame(0.05, Some(100.0)),
        ];
        let runs = voiced_runs(&contour);
        assert_eq!(runs, vec![(0, 1), (4, 5)]);
    }

    fn all_unvoiced() {
        let contour = vec![frame(0.0, None), frame(0.01, None)];
        let runs = voiced_runs(&contour);
        assert!(runs.is_empty());
    }

    fn empty_contour() {
        let runs = voiced_runs(&[]);
        assert!(runs.is_empty());
    }

    fn run_duration() {
        // 10 frames at 10ms hop = 100ms = 0.1s
        assert!((run_duration_secs(0, 9, 10.0) - 0.1).abs() < 0.001);
    }

    fn percentile_basic() {
        let data: Vec<f32> = (0..100).map(|i| i as f32).collect();
        assert!((percentile(&data, 0.0) - 0.0).abs() < 0.5);
        assert!((percentile(&data, 0.5) - 50.0).abs() < 1.0);
        assert!((percentile(&data, 1.0) - 99.0).abs() < 0.5);
    }

    fn merge_close_runs_empty() {
        assert_eq!(merge_close_runs(&[], 3), Vec::<(usize, usize)>::new());
    }

    fn merge_close_runs_single() {
        assert_eq!(merge_close_runs(&[(0, 10)], 3), vec![(0, 10)]);
    }

    fn merge_close_runs_bridges_small_gap() {
        // Two runs separated by 2 unvoiced frames → merged with max_gap=3
        let runs = vec![(0, 5), (8, 15)]; // gap: 8 - 5 - 1 = 2
        assert_eq!(merge_close_runs(&runs, 3), vec![(0, 15)]);
    }

    fn merge_close_runs_keeps_large_gap() {
        // Two runs separated by 10 unvoiced frames → NOT merged with max_gap=3
        let runs = vec![(0, 5), (16, 25)]; // gap: 16 - 5 - 1 = 10
        assert_eq!(merge_close_runs(&runs, 3), vec![(0, 5), (16, 25)]);
    }

    fn merge_close_runs_chains_multiple() {
        // Three runs with small gaps between each → all merge into one
        let runs = vec![(0, 10), (12, 20), (23, 30)]; // gaps: 1, 2
        assert_eq!(merge_close_runs(&runs, 3), vec![(0, 30)]);
    }

    fn merge_close_runs_mixed() {
        // Small gap then large gap → first two merge, third stays separate
        let runs = vec![(0, 10), (12, 20), (30, 40)]; // gaps: 1, 9
        assert_eq!(merge_close_runs(&runs, 3), vec![(0, 20), (30, 40)]);
    }

    fn default() -> Self {
        Self {
            frame_size_ms: 40.0,
            hop_size_ms: 10.0,
            quefrency_min_ms: 2.5,
            quefrency_max_ms: 16.7,
            energy_gate_relative_db: 30.0,
            smooth_time_frames: 5,    // 5 frames each side = 11 total (~110ms at 10ms hop)
            smooth_quefrency_bins: 5, // 5 bins each side = 11 total
        }
    }

pub fn compute_cpps(samples: &[f32], sample_rate: u32, config: &CppsConfig) -> Option<f32> {
    let sr = sample_rate as f32;
    let frame_size = (config.frame_size_ms / 1000.0 * sr) as usize;
    let hop_size = (config.hop_size_ms / 1000.0 * sr) as usize;

    if frame_size == 0 || hop_size == 0 || samples.len() < frame_size {
        return None;
    }

    let fft_size = frame_size.next_power_of_two();

    let mut planner = FftPlanner::new();
    let fft_forward = planner.plan_fft_forward(fft_size);
    let fft_inverse = planner.plan_fft_inverse(fft_size);

    // Quefrency range in samples
    let q_min = (config.quefrency_min_ms / 1000.0 * sr) as usize;
    let q_max = (config.quefrency_max_ms / 1000.0 * sr).ceil() as usize;
    let q_max = q_max.min(fft_size / 2 - 1);

    if q_min >= q_max {
        return None;
    }

    // Pass 1: compute per-frame RMS to find max energy for relative gate
    let mut frame_rms_values = Vec::new();
    let mut pos = 0;
    while pos + frame_

fn smooth_2d(
    matrix: &[Vec<f32>],
    num_frames: usize,
    num_quefrency: usize,
    half_time: usize,
    half_quefrency: usize,
) -> Vec<Vec<f32>> {
    // Build 2D prefix sum for O(1) region averages
    // prefix[t+1][q+1] = sum of matrix[0..=t][0..=q]
    let mut prefix = vec![vec![0.0f64; num_quefrency + 1]; num_frames + 1];

    for t in 0..num_frames {
        for q in 0..num_quefrency {
            prefix[t + 1][q + 1] = matrix[t][q] as f64
                + prefix[t][q + 1]
                + prefix[t + 1][q]
                - prefix[t][q];
        }
    }

    let mut result = vec![vec![0.0f32; num_quefrency]; num_frames];

    for t in 0..num_frames {
        let t0 = t.saturating_sub(half_time);
        let t1 = (t + half_time).min(num_frames - 1);

        for q in 0..num_quefrency {
            let q0 = q.saturating_sub(half_quefrency);
            let q1 = (q + half_quefrency).min(num_quefrency - 1);

            let sum = prefix[t1 + 1][q1 + 1]
                - 

fn cepstral_peak_prominence(row: &[f32], len: usize) -> Option<f32> {
    if len < 2 {
        return None;
    }

    // Find peak in dB domain (input is already in dB)
    let mut peak_val = f32::NEG_INFINITY;
    let mut peak_idx = 0;

    for (i, &v) in row.iter().take(len).enumerate() {
        if v > peak_val {
            peak_val = v;
            peak_idx = i;
        }
    }

    // Linear regression across the quefrency range
    let n = len as f32;
    let mut sum_x = 0.0f32;
    let mut sum_y = 0.0f32;
    let mut sum_xy = 0.0f32;
    let mut sum_xx = 0.0f32;

    for (i, &v) in row.iter().take(len).enumerate() {
        let x = i as f32;
        sum_x += x;
        sum_y += v;
        sum_xy += x * v;
        sum_xx += x * x;
    }

    let denom = n * sum_xx - sum_x * sum_x;
    if denom.abs() < 1e-10 {
        return Some(0.0);
    }

    let slope = (n * sum_xy - sum_x * sum_y) / denom;
    let intercept = (sum_y - slope * sum_x) / n;

    let regression_at_peak = slope

fn frame_rms(samples: &[f32]) -> f32 {
    if samples.is_empty() {
        return 0.0;
    }
    let sum_sq: f32 = samples.iter().map(|&s| s * s).sum();
    (sum_sq / samples.len() as f32).sqrt()
}

    fn sine_wave(freq_hz: f32, sample_rate: u32, duration_secs: f32) -> Vec<f32> {
        let num_samples = (sample_rate as f32 * duration_secs) as usize;
        (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                0.5 * (2.0 * PI * freq_hz * t).sin()
            })
            .collect()
    }

    fn white_noise(sample_rate: u32, duration_secs: f32, seed: u32) -> Vec<f32> {
        let n = (sample_rate as f32 * duration_secs) as usize;
        let mut rng_state: u32 = seed;
        (0..n)
            .map(|_| {
                rng_state = rng_state.wrapping_mul(1103515245).wrapping_add(12345);
                0.5 * ((rng_state as f32 / u32::MAX as f32) * 2.0 - 1.0)
            })
            .collect()
    }

    fn pure_tone_high_cpps() {
        let samples = sine_wave(100.0, 44100, 1.0);
        let cpps = compute_cpps(&samples, 44100, &CppsConfig::default());
        assert!(cpps.is_some(), "Should compute CPPS for a sine wave");
        let val = cpps.unwrap();
        assert!(
            val > 2.0,
            "Pure tone should have CPPS > 2 dB with smoothing, got {val:.2}"
        );
    }

    fn noise_low_cpps() {
        let samples = white_noise(44100, 1.0, 42);
        let cpps = compute_cpps(&samples, 44100, &CppsConfig::default());
        assert!(cpps.is_some(), "Should compute CPPS for noise");
        let val = cpps.unwrap();
        assert!(
            val < 5.0,
            "White noise should have low CPPS, got {val:.2}"
        );
    }

    fn silence_returns_none() {
        let samples = vec![0.0; 44100];
        let cpps = compute_cpps(&samples, 44100, &CppsConfig::default());
        assert!(cpps.is_none(), "Silence should return None (all frames gated)");
    }

    fn tone_higher_than_noise() {
        let tone = sine_wave(100.0, 44100, 1.0);
        let tone_cpps = compute_cpps(&tone, 44100, &CppsConfig::default()).unwrap();

        let noise = white_noise(44100, 1.0, 42);
        let noise_cpps = compute_cpps(&noise, 44100, &CppsConfig::default()).unwrap();

        assert!(
            tone_cpps > noise_cpps,
            "Tone CPPS ({tone_cpps:.2}) should exceed noise CPPS ({noise_cpps:.2})"
        );
    }

    fn smooth_2d_identity_with_zero_radius() {
        let matrix = vec![
            vec![1.0, 2.0, 3.0],
            vec![4.0, 5.0, 6.0],
            vec![7.0, 8.0, 9.0],
        ];
        let smoothed = smooth_2d(&matrix, 3, 3, 0, 0);
        for t in 0..3 {
            for q in 0..3 {
                assert!(
                    (smoothed[t][q] - matrix[t][q]).abs() < 1e-5,
                    "Zero-radius smoothing should be identity"
                );
            }
        }
    }

    fn smooth_2d_full_average() {
        // With radius large enough to cover the whole matrix, every cell = global mean
        let matrix = vec![
            vec![1.0, 2.0, 3.0],
            vec![4.0, 5.0, 6.0],
            vec![7.0, 8.0, 9.0],
        ];
        let smoothed = smooth_2d(&matrix, 3, 3, 10, 10);
        let global_mean = 45.0 / 9.0; // 5.0
        for t in 0..3 {
            for q in 0..3 {
                assert!(
                    (smoothed[t][q] - global_mean).abs() < 1e-4,
                    "Full-window smoothing should equal global mean, got {:.4} at ({t},{q})",
                    smoothed[t][q]
                );
            }
        }
    }

    fn relative_gate_adapts_to_level() {
        // Quiet signal (0.01 amplitude) should still work — gate is relative
        let samples: Vec<f32> = sine_wave(100.0, 44100, 1.0)
            .iter()
            .map(|&s| s * 0.01)
            .collect();
        let cpps = compute_cpps(&samples, 44100, &CppsConfig::default());
        assert!(cpps.is_some(), "Quiet signal should still compute CPPS with relative gate");
    }

    fn cpps_values_in_clinical_range_for_clean_tone() {
        // A clean periodic signal should produce CPPS roughly in the clinical
        // range. Pure sine is idealized, but values should be well above 3 dB.
        let samples = sine_wave(150.0, 44100, 2.0);
        let cpps = compute_cpps(&samples, 44100, &CppsConfig::default()).unwrap();
        assert!(
            cpps > 3.0,
            "Clean periodic signal should have CPPS > 3 dB (severe dysphonia cutoff), got {cpps:.2}"
        );
    }

pub fn compute_hnr_db(
    samples: &[f32],
    sample_rate: u32,
    contour: &[PitchFrame],
    hop_size_ms: f32,
) -> Option<f32> {
    let sr = sample_rate as f32;
    let hop_samples = (hop_size_ms / 1000.0 * sr) as usize;

    let mut hnr_values: Vec<f32> = Vec::new();

    for (i, frame) in contour.iter().enumerate() {
        let Some(f0) = frame.frequency else {
            continue;
        };

        // Pitch period in samples — the lag at which we compute autocorrelation
        let period_samples = (sr / f0).round() as usize;
        if period_samples == 0 {
            continue;
        }

        // We need at least 2 periods of audio to compute autocorrelation at this lag.
        // Center the window around the frame position.
        let center = i * hop_samples;
        let window_size = period_samples * 3; // 3 periods gives a stable estimate
        let start = center.saturating_sub(window_size / 2);
        let end = (start + window_size).min(samples.len());

   

pub(crate) fn normalized_autocorrelation(signal: &[f32], lag: usize) -> f32 {
    if lag >= signal.len() {
        return 0.0;
    }

    let n = signal.len() - lag;
    if n == 0 {
        return 0.0;
    }

    // Numerator: dot product of signal with shifted copy
    let mut cross_sum = 0.0_f64;
    // Denominators: energy of each segment
    let mut energy_a = 0.0_f64;
    let mut energy_b = 0.0_f64;

    for i in 0..n {
        let a = signal[i] as f64;
        let b = signal[i + lag] as f64;
        cross_sum += a * b;
        energy_a += a * a;
        energy_b += b * b;
    }

    let denom = (energy_a * energy_b).sqrt();
    if denom == 0.0 {
        return 0.0;
    }

    (cross_sum / denom) as f32
}

    fn frame(time: f32, freq: Option<f32>) -> PitchFrame {
        PitchFrame {
            time,
            frequency: freq,
        }
    }

    fn sine_wave(freq: f32, sample_rate: u32, duration: f32) -> Vec<f32> {
        let n = (sample_rate as f32 * duration) as usize;
        (0..n)
            .map(|i| (2.0 * PI * freq * i as f32 / sample_rate as f32).sin())
            .collect()
    }

    fn pure_tone_high_hnr() {
        // A pure sine wave has no noise → HNR should be very high
        let sr = 44100;
        let samples = sine_wave(100.0, sr, 1.0);

        let hop_ms = 10.0;
        let num_frames = 80; // ~0.8s worth
        let contour: Vec<_> = (0..num_frames)
            .map(|i| frame(i as f32 * hop_ms / 1000.0, Some(100.0)))
            .collect();

        let hnr = compute_hnr_db(&samples, sr, &contour, hop_ms).unwrap();
        assert!(
            hnr > 20.0,
            "Pure tone should have HNR > 20 dB, got {hnr:.1} dB"
        );
    }

    fn noisy_signal_low_hnr() {
        // Mix a sine wave with white noise at ~1:1 ratio → low HNR
        let sr = 44100;
        let n = (sr as f32 * 1.0) as usize;

        // Simple deterministic "noise" using a linear congruential generator.
        // Real white noise would be random, but we want reproducible tests.
        let mut rng_state: u32 = 42;
        let samples: Vec<f32> = (0..n)
            .map(|i| {
                let signal = (2.0 * PI * 100.0 * i as f32 / sr as f32).sin();
                // LCG pseudo-random noise in [-1, 1]
                rng_state = rng_state.wrapping_mul(1103515245).wrapping_add(12345);
                let noise = (rng_state as f32 / u32::MAX as f32) * 2.0 - 1.0;
                // Equal parts signal and noise
                0.5 * signal + 0.5 * noise
            })
            .collect();

        let hop_ms = 10.0;
        let num_frames = 80;
        let contour: Vec<_> = (0..num_frames)
            .map(|i| frame(i as f32 * hop_ms / 10

    fn autocorrelation_self() {
        // Autocorrelation at lag 0 should be 1.0 (signal correlates perfectly with itself)
        let signal: Vec<f32> = (0..100).map(|i| (i as f32 * 0.1).sin()).collect();
        let r = normalized_autocorrelation(&signal, 0);
        assert!(
            (r - 1.0).abs() < 0.001,
            "Self-correlation should be 1.0, got {r}"
        );
    }

    fn no_voiced_frames() {
        let contour = vec![frame(0.0, None), frame(0.01, None)];
        assert!(compute_hnr_db(&[0.0; 44100], 44100, &contour, 10.0).is_none());
    }

pub fn local_jitter_percent(contour: &[PitchFrame]) -> Option<f32> {
    // Collect consecutive voiced frequencies.
    // We need pairs of adjacent voiced frames — if there's an unvoiced frame
    // in between, the pair is broken (we can't measure cycle-to-cycle variation
    // across a gap).
    let mut periods: Vec<f32> = Vec::new();
    let mut perturbations: Vec<f32> = Vec::new();

    // prev_period tracks the period of the previous voiced frame.
    // It's reset to None whenever we hit an unvoiced frame.
    let mut prev_period: Option<f32> = None;

    for frame in contour {
        match frame.frequency {
            Some(f0) => {
                let period = 1.0 / f0;
                periods.push(period);

                if let Some(prev) = prev_period {
                    perturbations.push((period - prev).abs());
                }

                prev_period = Some(period);
            }
            None => {
                // Gap in voicing — reset the chain
       

pub fn local_jitter_percent_gated(
    contour: &[PitchFrame],
    frame_tiers: &[u8],
    hop_size_ms: f32,
) -> Option<f32> {
    if contour.len() != frame_tiers.len() {
        return None;
    }

    // Build a filtered contour: only keep tier 1/2 voiced frames,
    // treat tier 3 or unvoiced as gaps.
    let mut periods: Vec<f32> = Vec::new();
    let mut perturbations: Vec<f32> = Vec::new();
    let mut prev_period: Option<f32> = None;
    let mut consecutive = 0_usize;
    let mut max_consecutive = 0_usize;
    let mut total_gated_frames = 0_usize;

    for (frame, &tier) in contour.iter().zip(frame_tiers.iter()) {
        match frame.frequency {
            Some(f0) if tier <= 2 => {
                let period = 1.0 / f0;
                periods.push(period);
                total_gated_frames += 1;
                consecutive += 1;
                max_consecutive = max_consecutive.max(consecutive);

                if let Some(prev) = prev_period {
                    perturb

    fn frame(time: f32, freq: Option<f32>) -> PitchFrame {
        PitchFrame {
            time,
            frequency: freq,
        }
    }

    fn perfect_signal_zero_jitter() {
        // All frames at exactly 100 Hz — zero perturbation
        let contour: Vec<_> = (0..50)
            .map(|i| frame(i as f32 * 0.01, Some(100.0)))
            .collect();

        let jitter = local_jitter_percent(&contour).unwrap();
        assert!(
            jitter < 0.001,
            "Perfect signal should have ~0% jitter, got {jitter:.4}%"
        );
    }

    fn known_jitter() {
        // Alternating between 100 Hz and 110 Hz.
        // Periods: 0.01s and 0.00909s
        // |diff| = 0.000909... each time
        // Mean period = (0.01 + 0.00909) / 2 ≈ 0.009545
        // Jitter = 0.000909 / 0.009545 ≈ 9.5%
        let contour: Vec<_> = (0..20)
            .map(|i| {
                let f0 = if i % 2 == 0 { 100.0 } else { 110.0 };
                frame(i as f32 * 0.01, Some(f0))
            })
            .collect();

        let jitter = local_jitter_percent(&contour).unwrap();
        assert!(
            (jitter - 9.5).abs() < 1.0,
            "Expected ~9.5% jitter, got {jitter:.2}%"
        );
    }

    fn gap_breaks_chain() {
        // Two voiced segments separated by a gap — the gap should NOT
        // contribute a perturbation measurement
        let contour = vec![
            frame(0.0, Some(100.0)),
            frame(0.01, Some(100.0)),
            frame(0.02, None), // gap
            frame(0.03, Some(200.0)), // different pitch, but gap breaks chain
            frame(0.04, Some(200.0)),
        ];

        let jitter = local_jitter_percent(&contour).unwrap();
        // Within each segment, pitch is constant → 0% jitter
        assert!(
            jitter < 0.001,
            "Should be ~0% jitter (gap breaks the chain), got {jitter:.4}%"
        );
    }

    fn too_few_frames() {
        let contour = vec![frame(0.0, Some(100.0))];
        assert!(local_jitter_percent(&contour).is_none());
    }

    fn all_unvoiced() {
        let contour = vec![frame(0.0, None), frame(0.01, None)];
        assert!(local_jitter_percent(&contour).is_none());
    }

    fn gated_sufficient_tier1_data() {
        // 200 tier-1 frames at 100 Hz = 2s at 10ms hop
        let contour: Vec<_> = (0..200)
            .map(|i| frame(i as f32 * 0.01, Some(100.0)))
            .collect();
        let tiers = vec![1u8; 200];

        let jitter = local_jitter_percent_gated(&contour, &tiers, 10.0).unwrap();
        assert!(jitter < 0.001, "Perfect signal should have ~0% gated jitter, got {jitter:.4}%");
    }

    fn gated_rejects_tier3_frames() {
        // 200 frames all tier 3 — should fail
        let contour: Vec<_> = (0..200)
            .map(|i| frame(i as f32 * 0.01, Some(100.0)))
            .collect();
        let tiers = vec![3u8; 200];

        assert!(local_jitter_percent_gated(&contour, &tiers, 10.0).is_none());
    }

    fn gated_insufficient_consecutive() {
        // 10 tier-1, then 10 tier-3, then 10 tier-1... never reaching 15 consecutive
        let mut contour = Vec::new();
        let mut tiers = Vec::new();
        for chunk in 0..20 {
            for i in 0..10 {
                let idx = chunk * 10 + i;
                contour.push(frame(idx as f32 * 0.01, Some(100.0)));
                tiers.push(if chunk % 2 == 0 { 1 } else { 3 });
            }
        }

        assert!(local_jitter_percent_gated(&contour, &tiers, 10.0).is_none());
    }

pub fn max_phonation_time_secs(contour: &[PitchFrame], hop_size_ms: f32, max_bridge_ms: f32) -> f32 {
    let runs = contour::voiced_runs(contour);
    // Bridge gaps up to max_bridge_ms caused by pitch detector dropouts.
    //
    // MPT measures sustained phonation — the patient holds a vowel as long as
    // possible. Short gaps are detector artifacts (breathy signal loses
    // periodicity momentarily). Longer gaps indicate the patient actually
    // stopped (e.g., took a breath), which should end the measurement.
    let max_gap_frames = (max_bridge_ms / hop_size_ms) as usize;
    let runs = contour::merge_close_runs(&runs, max_gap_frames);

    runs.iter()
        .map(|&(start, end)| contour::run_duration_secs(start, end, hop_size_ms))
        .fold(0.0_f32, f32::max)
}

    fn frame(time: f32, freq: Option<f32>) -> PitchFrame {
        PitchFrame {
            time,
            frequency: freq,
        }
    }

    fn continuous_voicing() {
        // 100 voiced frames at 10ms hop = 1.0 seconds
        let contour: Vec<_> = (0..100)
            .map(|i| frame(i as f32 * 0.01, Some(100.0)))
            .collect();

        let mpt = max_phonation_time_secs(&contour, 10.0, 250.0);
        assert!(
            (mpt - 1.0).abs() < 0.02,
            "Expected ~1.0s MPT, got {mpt:.3}s"
        );
    }

    fn longest_run_wins() {
        // Two phonation attempts separated by a breathing pause (>500ms).
        // MPT should reflect the longer one.
        let mut contour: Vec<PitchFrame> = Vec::new();

        // First attempt: 200 frames = 2.0s
        for i in 0..200 {
            contour.push(frame(i as f32 * 0.01, Some(100.0)));
        }
        // Breathing pause: 60 frames = 600ms (> 500ms, not bridged)
        for i in 200..260 {
            contour.push(frame(i as f32 * 0.01, None));
        }
        // Second attempt: 500 frames = 5.0s
        for i in 260..760 {
            contour.push(frame(i as f32 * 0.01, Some(100.0)));
        }

        let mpt = max_phonation_time_secs(&contour, 10.0, 250.0);
        assert!(
            (mpt - 5.0).abs() < 0.1,
            "Expected ~5.0s (500 frames × 10ms), got {mpt:.3}s"
        );
    }

    fn bridges_detector_dropouts() {
        // Simulate pitch detector dropping frames in a 20s sustained vowel.
        // Gaps of 100ms and 200ms should be bridged (< 500ms threshold).
        let mut contour: Vec<PitchFrame> = Vec::new();
        // 10s voiced
        for i in 0..1000 {
            contour.push(frame(i as f32 * 0.01, Some(100.0)));
        }
        // 100ms gap (10 frames) — detector dropout
        for i in 1000..1010 {
            contour.push(frame(i as f32 * 0.01, None));
        }
        // 5s voiced
        for i in 1010..1510 {
            contour.push(frame(i as f32 * 0.01, Some(100.0)));
        }
        // 200ms gap (20 frames) — still a detector dropout
        for i in 1510..1530 {
            contour.push(frame(i as f32 * 0.01, None));
        }
        // 5s voiced
        for i in 1530..2030 {
            contour.push(frame(i as f32 * 0.01, Some(100.0)));
        }

        let mpt = max_phonation_time_secs(&contour, 10.0, 250.0);
        assert!(