timestretch 0.15.0

//! Phase vocoder time stretching with identity phase locking and sub-bass phase locking.

use crate::core::fft::{COMPLEX_ZERO, WINDOW_SUM_EPSILON, WINDOW_SUM_FLOOR_RATIO};
use crate::core::window::{WindowType, generate_window};
use crate::error::StretchError;
use crate::stretch::phase_locking::{PhaseLockingMode, apply_phase_locking_realtime};
use rustfft::{FftPlanner, num_complex::Complex};
use std::sync::Arc;

const TWO_PI_F64: f64 = 2.0 * std::f64::consts::PI;
/// Fraction of bins to pre-allocate for spectral peak detection (1/4 of bins).
const PEAKS_CAPACITY_DIVISOR: usize = 4;
/// Blend factor for phase gradient integration (soft vertical coherence).
const PHASE_GRADIENT_BLEND: f64 = 0.20;
/// Minimum magnitude to consider a bin as a spectral peak (avoids noise peaks).
use crate::stretch::phase_locking::MIN_PEAK_MAGNITUDE;
/// Treat values this close to integers as integral synthesis positions.
const SYNTH_POS_EPSILON: f64 = 1e-9;
/// Floor used in adaptive locking feature extraction.
const ADAPTIVE_FEATURE_EPS: f64 = 1e-12;
/// Flatness above this is considered noise-like (prefer ROI).
const ADAPTIVE_NOISY_FLATNESS: f64 = 0.72;
/// Crest below this is considered weakly-structured/noisy (prefer ROI).
const ADAPTIVE_NOISY_CREST: f64 = 2.5;
/// Harmonic confidence above this prefers identity locking near unity.
const ADAPTIVE_HARMONIC_CONFIDENCE: f64 = 4.5;
/// Harmonic confidence above this prefers selective locking on moderate ratios.
const ADAPTIVE_SELECTIVE_HARMONIC_CONFIDENCE: f64 = 2.8;
/// Max ratio-distance from unity where identity is preferred on harmonic frames.
const ADAPTIVE_IDENTITY_RATIO_DISTANCE_MAX: f64 = 0.35;
/// Max ratio-distance from unity where selective locking is preferred.
const ADAPTIVE_SELECTIVE_RATIO_DISTANCE_MAX: f64 = 0.65;
/// Ratio-distance above this always prefers ROI for stability.
const ADAPTIVE_FORCE_ROI_RATIO_DISTANCE: f64 = 0.75;
/// Number of synthesis frames to keep transient-focused locking active.
const TRANSIENT_FOCUS_FRAMES: usize = 3;
/// Phase vocoder state for time stretching.
pub struct PhaseVocoder {
    fft_size: usize,
    hop_analysis: usize,
    stretch_ratio: f64,
    /// Absolute synthesis position (in samples) of the next frame start.
    synthesis_pos: f64,
    /// Number of synthesized samples already emitted by the streaming path.
    synthesis_emitted: usize,
    window: Vec<f32>,
    /// Phase accumulator for resynthesis (f64 for precision over long signals).
    phase_accum: Vec<f64>,
    /// Previous analysis phase (f64 to match accumulator precision).
    prev_phase: Vec<f64>,
    /// Marks bins that must be seeded from fresh analysis phase on next frame.
    ///
    /// Used by selective band resets so transient-triggered partial resets do
    /// not compute an invalid instantaneous-frequency jump from zeroed phase
    /// history.
    phase_seed_pending: Vec<bool>,
    /// Pre-planned forward FFT for analysis frames.
    fft_forward: Arc<dyn rustfft::Fft<f32>>,
    /// Pre-planned inverse FFT for synthesis frames.
    fft_inverse: Arc<dyn rustfft::Fft<f32>>,
    /// Scratch buffer for forward FFT execution.
    fft_forward_scratch: Vec<Complex<f32>>,
    /// Scratch buffer for inverse FFT execution.
    fft_inverse_scratch: Vec<Complex<f32>>,
    /// Pre-computed expected phase advance per bin (f64 for precision).
    expected_phase_advance: Vec<f64>,
    /// Reusable FFT buffer.
    fft_buffer: Vec<Complex<f32>>,
    /// Reusable magnitude buffer.
    magnitudes: Vec<f32>,
    /// Reusable phase buffer.
    new_phases: Vec<f32>,
    /// Reusable peaks buffer for identity phase locking.
    peaks: Vec<usize>,
    /// Reusable trough buffer for phase-lock influence-region discovery.
    phase_lock_troughs: Vec<usize>,
    /// Reusable backup of pre-lock phases for ROI clamping.
    phase_lock_pv_phases: Vec<f32>,
    /// Current frame's analysis phases (for identity phase locking).
    analysis_phases: Vec<f32>,
    /// Bin index at or below which sub-bass phase locking is applied.
    sub_bass_bin: usize,
    /// Phase locking algorithm to use.
    phase_locking_mode: PhaseLockingMode,
    /// Enables confidence-driven adaptive phase-lock mode switching.
    adaptive_phase_locking: bool,
    /// When positive, holds the phase-gradient coherence blend at this
    /// strength at wide ratios instead of tapering it out (batch/QA
    /// experiment knob, see ROADMAP Stage 11; the live engine never sets
    /// this). 0.0 = shipped taper behavior.
    wide_coherence_blend: f64,
    /// Short-lived state that tightens phase behavior right after a transient reset.
    transient_focus_frames: usize,
    /// Precomputed overlap-add gain (`synthesis_window / fft_size`).
    ola_gain: Vec<f32>,
    /// f64 view of `ola_gain` to avoid per-sample casts in fractional OLA path.
    ola_gain_f64: Vec<f64>,
    /// Precomputed window product (`analysis_window * synthesis_window`) for OLA normalization.
    ola_window_product: Vec<f32>,
    /// f64 view of `ola_window_product` to avoid per-sample casts in fractional OLA path.
    ola_window_product_f64: Vec<f64>,
    /// Backup of IF-estimated phases before phase locking overwrites them.
    /// Used to blend IF estimates with locked phases for non-peak bins.
    if_phases_backup: Vec<f32>,
    /// Reusable output buffer (avoids allocation per process() call).
    output_buf: Vec<f32>,
    /// Reusable window sum buffer (avoids allocation per process() call).
    window_sum_buf: Vec<f32>,
    /// Unnormalized overlap-add tail carried between streaming calls.
    streaming_tail: Vec<f32>,
    /// Window-sum tail matching `streaming_tail`.
    streaming_tail_window_sum: Vec<f32>,
    /// Most conservative ratio contributing to the carried streaming tail.
    ///
    /// When a tail spans a ratio change, keep the larger expansion ratio so
    /// the next normalization pass does not suddenly tighten the overlap
    /// floor at the chunk boundary.
    streaming_tail_ratio: f64,
    /// Reusable accumulation buffer for streaming overlap-add.
    streaming_accum_output: Vec<f32>,
    /// Reusable window-sum accumulation buffer for streaming overlap-add.
    streaming_accum_window_sum: Vec<f32>,
}

#[inline]
fn streaming_tail_normalize_ratio(carried_tail_ratio: f64, current_ratio: f64) -> f64 {
    carried_tail_ratio.max(current_ratio)
}

impl PhaseVocoder {
    /// Creates a new phase vocoder.
    pub fn new(
        fft_size: usize,
        hop_analysis: usize,
        stretch_ratio: f64,
        sample_rate: u32,
        sub_bass_cutoff: f32,
    ) -> Self {
        Self::with_window(
            fft_size,
            hop_analysis,
            stretch_ratio,
            sample_rate,
            sub_bass_cutoff,
            WindowType::BlackmanHarris,
        )
    }

    /// Creates a new phase vocoder with a specific window function.
    pub fn with_window(
        fft_size: usize,
        hop_analysis: usize,
        stretch_ratio: f64,
        sample_rate: u32,
        sub_bass_cutoff: f32,
        window_type: WindowType,
    ) -> Self {
        Self::with_options(
            fft_size,
            hop_analysis,
            stretch_ratio,
            sample_rate,
            sub_bass_cutoff,
            window_type,
            PhaseLockingMode::RegionOfInfluence,
        )
    }

    /// Creates a new phase vocoder with full configuration options.
    pub fn with_options(
        fft_size: usize,
        hop_analysis: usize,
        stretch_ratio: f64,
        sample_rate: u32,
        sub_bass_cutoff: f32,
        window_type: WindowType,
        phase_locking_mode: PhaseLockingMode,
    ) -> Self {
        let window = generate_window(window_type, fft_size);
        // Match synthesis window to analysis window type for a proper window product.
        // Using the same window type ensures the overlap-add normalization works
        // correctly and avoids spectral distortion from mismatched window shapes.
        //
        // Exception: BlackmanHarris analysis uses Hann for synthesis because BH^2
        // has poor COLA (constant overlap-add) properties at standard 75% overlap
        // (hop = fft_size/4). The BH*Hann product provides better overlap-add
        // flatness while BH still provides excellent sidelobe suppression for
        // the analysis stage.
        let synthesis_window_type = match window_type {
            WindowType::BlackmanHarris => WindowType::Hann,
            other => other,
        };
        let synthesis_window = generate_window(synthesis_window_type, fft_size);
        let inv_fft = 1.0 / fft_size as f32;
        let ola_gain: Vec<f32> = synthesis_window.iter().map(|&w| w * inv_fft).collect();
        let ola_gain_f64: Vec<f64> = ola_gain.iter().map(|&w| w as f64).collect();
        let ola_window_product: Vec<f32> = window
            .iter()
            .zip(synthesis_window.iter())
            .map(|(&a, &b)| a * b)
            .collect();
        let ola_window_product_f64: Vec<f64> =
            ola_window_product.iter().map(|&w| w as f64).collect();
        let num_bins = fft_size / 2 + 1;
        let mut planner = FftPlanner::new();
        let fft_forward = planner.plan_fft_forward(fft_size);
        let fft_inverse = planner.plan_fft_inverse(fft_size);
        let fft_forward_scratch_len = fft_forward.get_inplace_scratch_len();
        let fft_inverse_scratch_len = fft_inverse.get_inplace_scratch_len();

        let expected_phase_advance: Vec<f64> = (0..num_bins)
            .map(|bin| TWO_PI_F64 * bin as f64 * hop_analysis as f64 / fft_size as f64)
            .collect();

        // Compute the bin index for the sub-bass cutoff frequency.
        // Bins at or below this index get rigid phase locking to prevent
        // phase cancellation in the critical sub-bass region.
        let sub_bass_bin =
            (sub_bass_cutoff * fft_size as f32 / sample_rate as f32).round() as usize;
        let sub_bass_bin = sub_bass_bin.min(num_bins);

        Self {
            fft_size,
            hop_analysis,
            stretch_ratio,
            synthesis_pos: 0.0,
            synthesis_emitted: 0,
            window,
            phase_accum: vec![0.0f64; num_bins],
            prev_phase: vec![0.0f64; num_bins],
            phase_seed_pending: vec![true; num_bins],
            fft_forward,
            fft_inverse,
            fft_forward_scratch: vec![COMPLEX_ZERO; fft_forward_scratch_len],
            fft_inverse_scratch: vec![COMPLEX_ZERO; fft_inverse_scratch_len],
            expected_phase_advance,
            fft_buffer: vec![COMPLEX_ZERO; fft_size],
            magnitudes: vec![0.0; num_bins],
            new_phases: vec![0.0; num_bins],
            peaks: Vec::with_capacity(num_bins / PEAKS_CAPACITY_DIVISOR),
            phase_lock_troughs: Vec::with_capacity(num_bins / 2),
            phase_lock_pv_phases: Vec::with_capacity(num_bins),
            analysis_phases: vec![0.0; num_bins],
            sub_bass_bin,
            phase_locking_mode,
            adaptive_phase_locking: false,
            wide_coherence_blend: 0.0,
            transient_focus_frames: 0,
            ola_gain,
            ola_gain_f64,
            ola_window_product,
            ola_window_product_f64,
            if_phases_backup: vec![0.0; num_bins],
            output_buf: Vec::new(),
            window_sum_buf: Vec::new(),
            streaming_tail: Vec::new(),
            streaming_tail_window_sum: Vec::new(),
            streaming_tail_ratio: stretch_ratio,
            streaming_accum_output: Vec::new(),
            streaming_accum_window_sum: Vec::new(),
        }
    }

    /// Returns the FFT size.
    #[inline]
    pub fn fft_size(&self) -> usize {
        self.fft_size
    }

    /// Returns the analysis hop size.
    #[inline]
    pub fn hop_analysis(&self) -> usize {
        self.hop_analysis
    }

    /// Returns the sub-bass bin cutoff index.
    #[inline]
    pub fn sub_bass_bin(&self) -> usize {
        self.sub_bass_bin
    }

    /// Updates the stretch ratio without resetting phase state.
    ///
    /// Phase advance follows the new ratio from the next synthesized frame
    /// while all accumulated phase information is preserved. Production
    /// callers update the ratio as a smooth, small-step stream (the wide
    /// keylock stage slews it per block), so no seam masking is needed.
    /// Any overlap tail carried across a streaming emit boundary keeps its
    /// own normalization ratio (see `streaming_tail_ratio`), so output
    /// rendered at the previous ratio is still normalized at that ratio.
    #[inline]
    pub fn set_stretch_ratio(&mut self, stretch_ratio: f64) {
        self.stretch_ratio = stretch_ratio;
    }

    /// Resets the phase accumulator and previous-phase buffers.
    ///
    /// Call this at transient boundaries so that stale phase state from a
    /// previous tonal segment does not contaminate the next one. The PV will
    /// re-derive phases from the first analysis frame after the reset.
    #[inline]
    pub fn reset_phase_state(&mut self) {
        self.phase_accum.fill(0.0);
        self.prev_phase.fill(0.0);
        self.phase_seed_pending.fill(true);
        self.transient_focus_frames = 0;
    }

    /// Clears ALL per-stream state — phase history, the held streaming
    /// overlap tail, and synthesis position — without deallocating any
    /// buffer.
    ///
    /// After this call the vocoder behaves like a freshly constructed one
    /// for streaming purposes, which makes it the allocation-free
    /// equivalent of dropping and rebuilding the instance. Used by
    /// warm-start seek to re-prime from new material.
    pub fn reset_streaming_state(&mut self) {
        self.reset_phase_state();
        self.streaming_tail.clear();
        self.streaming_tail_window_sum.clear();
        self.streaming_tail_ratio = self.stretch_ratio;
        self.synthesis_pos = 0.0;
        self.synthesis_emitted = 0;
    }

    /// Enables or disables confidence-driven adaptive phase-lock switching.
    #[inline]
    pub fn set_adaptive_phase_locking(&mut self, enabled: bool) {
        self.adaptive_phase_locking = enabled;
    }

    /// Returns whether adaptive phase-lock switching is enabled.
    #[inline]
    pub fn adaptive_phase_locking(&self) -> bool {
        self.adaptive_phase_locking
    }

    /// Holds the phase-gradient coherence blend at `blend` strength at
    /// wide ratios instead of tapering it out past |ratio - 1| (Stage 11
    /// falsification knob: the taper leaves expansions ≥ 2x with no
    /// vertical-coherence help, the confirmed cause of robotic
    /// slowdowns). `0.0` restores the shipped taper; the shipped
    /// near-unity strength is 0.20. Clamped to `[0.0, 1.0]`.
    #[inline]
    pub fn set_wide_ratio_coherence_blend(&mut self, blend: f64) {
        self.wide_coherence_blend = blend.clamp(0.0, 1.0);
    }

    /// Selectively resets phase state for specific frequency bands.
    ///
    /// Only zeros `phase_accum` and `prev_phase` for bins within the bands
    /// indicated by `reset_mask`: `[sub_bass, low, mid, high]`.
    /// Band boundaries: sub-bass <100Hz, low 100-500Hz, mid 500-4000Hz, high >4000Hz.
    ///
    /// This avoids disrupting phase tracking in bands where no transient occurred
    /// (e.g., a hi-hat hit shouldn't reset the sustained bass phase).
    pub fn reset_phase_state_bands(&mut self, reset_mask: [bool; 4], sample_rate: u32) {
        let num_bins = self.fft_size / 2 + 1;
        let bin_freq = sample_rate as f32 / self.fft_size as f32;

        for bin in 0..num_bins {
            let freq = bin as f32 * bin_freq;
            let band_idx = if freq < 100.0 {
                0
            } else if freq < 500.0 {
                1
            } else if freq < 4000.0 {
                2
            } else {
                3
            };
            if reset_mask[band_idx] {
                self.phase_accum[bin] = 0.0;
                self.prev_phase[bin] = 0.0;
                self.phase_seed_pending[bin] = true;
            }
        }

        // Engage a short transient-focus window for audible bands so the
        // first few post-reset frames favor tighter attack coherence.
        if reset_mask[1] || reset_mask[2] || reset_mask[3] {
            self.transient_focus_frames = self.transient_focus_frames.max(TRANSIENT_FOCUS_FRAMES);
        }
    }

    /// Stretches a mono audio signal using phase vocoder with identity phase locking.
    pub fn process(&mut self, input: &[f32]) -> Result<Vec<f32>, StretchError> {
        // Batch calls are independent; clear any prior streaming overlap state.
        self.streaming_tail.clear();
        self.streaming_tail_window_sum.clear();
        self.streaming_tail_ratio = self.stretch_ratio;

        // Mirror-pad input to stabilize edge normalization.
        //
        // The first/last analysis frames have incomplete window overlap,
        // producing an unstable window-sum profile that distorts the
        // normalized output at signal edges. Padding with reflected
        // samples lets the overlap-add window sum reach its steady-state
        // value before processing actual content, eliminating edge
        // artifacts that degrade spectral metrics.
        //
        // Asymmetric padding: the start uses a shorter mirror so that
        // onsets at t=0 remain distinct — longer mirrors pre-condition
        // the phase state with identical spectral content and mask the
        // amplitude transition, causing onset detectors to miss it.
        // The end keeps a longer mirror (hop*8) for full spectral quality.
        //
        // The start padding scales with the ratio distance from unity:
        // at extreme ratios (>0.3 from 1.0) the onset time-scaling
        // amplifies small phase-state artefacts enough to shift onsets
        // beyond the scoring tolerance, so we use a shorter mirror.
        // Near unity the phase state is well-behaved and the longer
        // mirror gives better spectral metrics.
        // Expansion ratios produce longer output, so the PV processes
        // further into the tail padding region.  Extra end padding gives
        // the window-sum normalization more frames to stabilise, reducing
        // amplitude droop at the output tail that degrades LSD.
        let end_pad_mult = if self.stretch_ratio > 1.1 { 10 } else { 8 };
        let end_pad = (self.hop_analysis * end_pad_mult).min(input.len());
        // Graduated start padding: shorter mirrors preserve onset sharpness
        // (critical for TP scoring) while longer mirrors give better spectral
        // metrics (SC/LSD).  Ratios 0.2-0.3 from unity get an intermediate
        // 6-hop padding that balances both concerns.
        let ratio_dist = (self.stretch_ratio - 1.0).abs();
        let start_pad_mult = if ratio_dist > 0.3 {
            4
        } else if ratio_dist > 0.15 {
            6
        } else {
            8
        };
        let start_pad = (self.hop_analysis * start_pad_mult).min(input.len());
        if start_pad > 0 && input.len() >= self.fft_size {
            let padded_len = input.len() + start_pad + end_pad;
            let mut padded = vec![0.0f32; padded_len];
            // Reflect start (shorter — preserves onset sharpness)
            for i in 0..start_pad {
                padded[i] = input[start_pad - 1 - i];
            }
            // Copy original
            padded[start_pad..start_pad + input.len()].copy_from_slice(input);
            // Reflect end with cosine taper: the reflected samples fade
            // smoothly to zero so the PV sees a gradually decaying
            // continuation instead of a cusp.  This reduces phase
            // interference at the boundary and eliminates the sharp
            // amplitude droop in the last few output frames that causes
            // false onset detection at the signal end.
            for i in 0..end_pad {
                let t = (i + 1) as f32 / end_pad as f32;
                let fade = 0.5 * (1.0 + (std::f32::consts::PI * t).cos());
                padded[start_pad + input.len() + i] = input[input.len() - 1 - i] * fade;
            }

            let (_num_frames, output_len) = self.process_core(&padded, true)?;
            let mut output = self.output_buf[..output_len].to_vec();
            Self::normalize_output(
                &mut output,
                &self.window_sum_buf[..output_len],
                self.stretch_ratio,
            );

            // Trim output to remove padding artifacts.
            let trim_start = (start_pad as f64 * self.stretch_ratio).round() as usize;
            let expected_len = (input.len() as f64 * self.stretch_ratio).round() as usize;
            let trim_end = (trim_start + expected_len).min(output.len());
            if trim_start < output.len() {
                let mut result = output[trim_start..trim_end].to_vec();

                // Edge amplitude correction: mirror padding causes phase
                // interference near the signal boundary, reducing amplitude
                // in the first ~fft_size*ratio output samples.  Measure
                // this droop and apply a smooth gain ramp to match the
                // steady-state level, preventing false onset detection.
                let correction_len =
                    (self.fft_size as f64 * self.stretch_ratio.abs().max(1.0)).ceil() as usize;
                let correction_len = correction_len.min(result.len() / 4);
                if correction_len >= 128 && result.len() >= correction_len * 4 {
                    let ss_start = correction_len;
                    let ss_end = (correction_len * 3).min(result.len());
                    let ss_len = ss_end - ss_start;
                    let ss_energy: f64 = result[ss_start..ss_end]
                        .iter()
                        .map(|&s| (s as f64) * (s as f64))
                        .sum();
                    let ss_rms = (ss_energy / ss_len as f64).sqrt();
                    if ss_rms > 1e-6 {
                        let edge_energy: f64 = result[..correction_len]
                            .iter()
                            .map(|&s| (s as f64) * (s as f64))
                            .sum();
                        let edge_rms = (edge_energy / correction_len as f64).sqrt();
                        // Only correct if the edge is significantly quieter
                        if edge_rms < ss_rms * 0.7 && edge_rms > 1e-8 {
                            let gain = (ss_rms / edge_rms).min(4.0) as f32;
                            for (i, sample) in result.iter_mut().enumerate().take(correction_len) {
                                // Cubic Hermite ramp: zero derivative at both
                                // endpoints for smoother spectral transition.
                                let t = i as f32 / correction_len as f32;
                                let g = 1.0 + (gain - 1.0) * (1.0 - 3.0 * t * t + 2.0 * t * t * t);
                                *sample *= g;
                            }
                        }
                    }
                }

                // Same correction for the end edge (mirror padding droop).
                if correction_len >= 128 && result.len() >= correction_len * 4 {
                    let ss_start = result.len() / 4;
                    let ss_end = result.len() * 3 / 4;
                    let ss_len = ss_end - ss_start;
                    let ss_energy: f64 = result[ss_start..ss_end]
                        .iter()
                        .map(|&s| (s as f64) * (s as f64))
                        .sum();
                    let ss_rms = (ss_energy / ss_len as f64).sqrt();
                    if ss_rms > 1e-6 {
                        let end_start = result.len() - correction_len;
                        let end_energy: f64 = result[end_start..]
                            .iter()
                            .map(|&s| (s as f64) * (s as f64))
                            .sum();
                        let end_rms = (end_energy / correction_len as f64).sqrt();
                        if end_rms < ss_rms * 0.7 && end_rms > 1e-8 {
                            let gain = (ss_rms / end_rms).min(4.0) as f32;
                            let rlen = result.len();
                            for i in 0..correction_len {
                                // Cubic Hermite ramp: zero derivative at both
                                // endpoints (matches start correction shape).
                                let t = i as f32 / correction_len as f32;
                                let g = 1.0 + (gain - 1.0) * (3.0 * t * t - 2.0 * t * t * t);
                                result[rlen - correction_len + i] *= g;
                            }
                        }
                    }
                }

                return Ok(result);
            }
        }

        // Fallback: process without padding (short inputs or edge cases).
        let (_num_frames, output_len) = self.process_core(input, true)?;
        let mut output = self.output_buf[..output_len].to_vec();
        Self::normalize_output(
            &mut output,
            &self.window_sum_buf[..output_len],
            self.stretch_ratio,
        );
        Ok(output)
    }

    /// Pre-allocates all streaming-path buffers for a maximum input window.
    ///
    /// Call once at build time with the caller's rolling-window capacity so
    /// [`Self::process_streaming_into`] performs no allocations afterwards,
    /// even as the stretch ratio moves (bounded by `max_ratio`) and the
    /// render window size varies between calls.
    pub fn reserve_streaming_capacity(&mut self, max_window_frames: usize, max_ratio: f64) {
        fn reserve_to(buf: &mut Vec<f32>, capacity: usize) {
            if buf.capacity() < capacity {
                buf.reserve(capacity - buf.len());
            }
        }
        let ratio_mult = (max_ratio.max(1.0).ceil() as usize).saturating_add(1);
        let out_bound = max_window_frames
            .saturating_mul(ratio_mult)
            .saturating_add(self.fft_size.saturating_mul(2));
        reserve_to(&mut self.output_buf, out_bound);
        reserve_to(&mut self.window_sum_buf, out_bound);
        reserve_to(&mut self.streaming_accum_output, out_bound);
        reserve_to(&mut self.streaming_accum_window_sum, out_bound);
        reserve_to(&mut self.streaming_tail, out_bound);
        reserve_to(&mut self.streaming_tail_window_sum, out_bound);
        let num_bins = self.fft_size / 2 + 1;
        if self.peaks.capacity() < num_bins {
            self.peaks.reserve(num_bins - self.peaks.len());
        }
    }

    /// Streaming phase-vocoder pass that preserves phase across calls.
    ///
    /// `input` should include any required analysis overlap context from the
    /// caller (typically managed by a higher-level stream processor). This
    /// method keeps synthesis overlap/window tails internally and emits only
    /// hop-aligned samples that are final for this call.
    pub fn process_streaming(&mut self, input: &[f32]) -> Result<Vec<f32>, StretchError> {
        let mut output = Vec::with_capacity(
            ((input.len() as f64 * self.stretch_ratio).ceil() as usize)
                .saturating_add(self.fft_size),
        );
        self.process_streaming_into(input, &mut output)?;
        Ok(output)
    }

    /// Streaming phase-vocoder pass writing directly into `output`.
    ///
    /// This avoids temporary output allocations in real-time paths.
    pub fn process_streaming_into(
        &mut self,
        input: &[f32],
        output: &mut Vec<f32>,
    ) -> Result<(), StretchError> {
        if input.len() < self.fft_size {
            output.clear();
            return Ok(());
        }

        let (emit_len, output_len) = self.process_core(input, false)?;
        if output.capacity() < emit_len {
            return Err(StretchError::BufferOverflow {
                buffer: "phase_vocoder_stream_output",
                requested: emit_len,
                available: output.capacity(),
            });
        }

        let work_len = output_len
            .max(emit_len)
            .max(self.streaming_tail.len())
            .max(self.streaming_tail_window_sum.len());

        self.streaming_accum_output.resize(work_len, 0.0);
        self.streaming_accum_output.fill(0.0);
        self.streaming_accum_window_sum.resize(work_len, 0.0);
        self.streaming_accum_window_sum.fill(0.0);

        self.streaming_accum_output[..output_len].copy_from_slice(&self.output_buf[..output_len]);
        self.streaming_accum_window_sum[..output_len]
            .copy_from_slice(&self.window_sum_buf[..output_len]);

        let carried_tail_ratio = self.streaming_tail_ratio;
        let tail_len = self
            .streaming_tail
            .len()
            .min(self.streaming_tail_window_sum.len());
        for i in 0..tail_len {
            self.streaming_accum_output[i] += self.streaming_tail[i];
            self.streaming_accum_window_sum[i] += self.streaming_tail_window_sum[i];
        }

        // Keep the unresolved overlap region for the next chunk.
        self.streaming_tail.clear();
        self.streaming_tail_window_sum.clear();
        if emit_len < work_len {
            self.streaming_tail
                .extend_from_slice(&self.streaming_accum_output[emit_len..work_len]);
            self.streaming_tail_window_sum
                .extend_from_slice(&self.streaming_accum_window_sum[emit_len..work_len]);
        }

        output.resize(emit_len, 0.0);
        output[..emit_len].copy_from_slice(&self.streaming_accum_output[..emit_len]);
        let emitted_window_sum = &self.streaming_accum_window_sum[..emit_len];
        let emitted_max_window_sum = emitted_window_sum.iter().copied().fold(0.0f32, f32::max);
        let seam_len = tail_len.min(emit_len);
        if seam_len == 0 {
            Self::normalize_output_with_window_floor(
                output,
                emitted_window_sum,
                self.stretch_ratio,
                emitted_max_window_sum,
            );
        } else if seam_len == emit_len {
            Self::normalize_output_with_window_floor(
                output,
                emitted_window_sum,
                streaming_tail_normalize_ratio(carried_tail_ratio, self.stretch_ratio),
                emitted_max_window_sum,
            );
        } else {
            // Once the carried seam has fully drained inside this callback, switch
            // the remainder back to the current ratio's normalization floor so a
            // previous expansion ratio does not keep loosening the whole chunk.
            Self::normalize_output_with_window_floor(
                &mut output[..seam_len],
                &emitted_window_sum[..seam_len],
                streaming_tail_normalize_ratio(carried_tail_ratio, self.stretch_ratio),
                emitted_max_window_sum,
            );
            Self::normalize_output_with_window_floor(
                &mut output[seam_len..emit_len],
                &emitted_window_sum[seam_len..emit_len],
                self.stretch_ratio,
                emitted_max_window_sum,
            );
        }
        self.synthesis_emitted = self.synthesis_emitted.saturating_add(emit_len);
        self.streaming_tail_ratio = if self.streaming_tail.is_empty() {
            self.stretch_ratio
        } else if emit_len < tail_len {
            // Preserve the more conservative seam ratio only while unresolved
            // overlap from the previous chunk is still present in the carried tail.
            streaming_tail_normalize_ratio(carried_tail_ratio, self.stretch_ratio)
        } else {
            self.stretch_ratio
        };
        Ok(())
    }

    /// Flushes remaining streaming overlap/window tail at end of stream.
    pub fn flush_streaming(&mut self) -> Result<Vec<f32>, StretchError> {
        let mut output = Vec::with_capacity(self.streaming_tail.len());
        self.flush_streaming_into(&mut output)?;
        Ok(output)
    }

    /// Flushes remaining streaming overlap/window tail into `output`.
    pub fn flush_streaming_into(&mut self, output: &mut Vec<f32>) -> Result<(), StretchError> {
        if self.streaming_tail.is_empty() || self.streaming_tail_window_sum.is_empty() {
            self.streaming_tail.clear();
            self.streaming_tail_window_sum.clear();
            self.streaming_tail_ratio = self.stretch_ratio;
            self.synthesis_pos = 0.0;
            self.synthesis_emitted = 0;
            output.clear();
            return Ok(());
        }

        let len = self
            .streaming_tail
            .len()
            .min(self.streaming_tail_window_sum.len());
        if output.capacity() < len {
            return Err(StretchError::BufferOverflow {
                buffer: "phase_vocoder_flush_output",
                requested: len,
                available: output.capacity(),
            });
        }
        output.resize(len, 0.0);
        output.copy_from_slice(&self.streaming_tail[..len]);
        Self::normalize_output(
            output,
            &self.streaming_tail_window_sum[..len],
            streaming_tail_normalize_ratio(self.streaming_tail_ratio, self.stretch_ratio),
        );
        self.streaming_tail.clear();
        self.streaming_tail_window_sum.clear();
        self.streaming_tail_ratio = self.stretch_ratio;
        self.synthesis_pos = 0.0;
        self.synthesis_emitted = 0;
        Ok(())
    }

    /// Shared PV core used by both batch and streaming paths.
    ///
    /// Returns `(emit_len, output_len)` where `output_len` samples are
    /// accumulated (unnormalized) into `self.output_buf` and
    /// `self.window_sum_buf`, and `emit_len` is the number of samples
    /// finalized for streaming emission (`floor(next_synthesis_pos)`).
    fn process_core(
        &mut self,
        input: &[f32],
        reset_phase_state: bool,
    ) -> Result<(usize, usize), StretchError> {
        if input.len() < self.fft_size {
            return Err(StretchError::InputTooShort {
                provided: input.len(),
                minimum: self.fft_size,
            });
        }

        let num_bins = self.fft_size / 2 + 1;
        let num_frames = (input.len() - self.fft_size) / self.hop_analysis + 1;

        if reset_phase_state {
            self.phase_accum.fill(0.0);
            self.prev_phase.fill(0.0);
            self.phase_seed_pending.fill(true);
            self.transient_focus_frames = 0;
            self.synthesis_pos = 0.0;
            self.synthesis_emitted = 0;
        }

        let hop_ratio = self.stretch_ratio;
        let frame_advance = self.hop_analysis as f64 * hop_ratio;
        let fft_forward = Arc::clone(&self.fft_forward);
        let fft_inverse = Arc::clone(&self.fft_inverse);

        // Local synthesis timeline starts at the current emission cursor.
        let start_synthesis_pos =
            snap_near_integer((self.synthesis_pos - self.synthesis_emitted as f64).max(0.0));

        // Pre-compute required accumulation length with fractional placement.
        let mut max_write_idx = 0usize;
        let mut synthesis_scan_pos = start_synthesis_pos;
        for _ in 0..num_frames {
            let synthesis_pos = snap_near_integer(synthesis_scan_pos);
            let synthesis_floor = synthesis_pos.floor() as usize;
            let frac = synthesis_pos - synthesis_floor as f64;
            let frame_end = synthesis_floor.saturating_add(
                self.fft_size
                    .saturating_sub(1)
                    .saturating_add(usize::from(frac > SYNTH_POS_EPSILON)),
            );
            max_write_idx = max_write_idx.max(frame_end);
            synthesis_scan_pos = synthesis_pos + frame_advance;
        }
        let output_len = max_write_idx.saturating_add(1);
        let emit_len = snap_near_integer(synthesis_scan_pos).floor() as usize;

        // Reuse pre-allocated buffers, growing if needed (never shrinks).
        self.output_buf.resize(output_len, 0.0);
        self.output_buf.fill(0.0);
        self.window_sum_buf.resize(output_len, 0.0);
        self.window_sum_buf.fill(0.0);

        let mut synthesis_frame_pos = start_synthesis_pos;
        for frame_idx in 0..num_frames {
            let analysis_pos = frame_idx * self.hop_analysis;
            let synthesis_pos = snap_near_integer(synthesis_frame_pos);
            let synthesis_floor = synthesis_pos.floor() as usize;
            let frac = synthesis_pos - synthesis_floor as f64;

            self.analyze_frame(
                &input[analysis_pos..analysis_pos + self.fft_size],
                &fft_forward,
            );
            self.advance_phases(num_bins, hop_ratio);

            // Save IF-estimated phases before phase locking overwrites them.
            self.if_phases_backup[..num_bins].copy_from_slice(&self.new_phases[..num_bins]);

            // Phase locking: lock non-peak bins to their nearest peak using
            // the analysis phase relationship. Only applies above the sub-bass region.
            let locking_mode = self.select_phase_locking_mode(num_bins);
            apply_phase_locking_realtime(
                locking_mode,
                &self.magnitudes,
                &self.analysis_phases,
                &mut self.new_phases,
                num_bins,
                self.sub_bass_bin,
                &mut self.peaks,
                &mut self.phase_lock_troughs,
                &mut self.phase_lock_pv_phases,
            );

            // Blend IF estimates with locked phases for non-peak bins above sub-bass.
            // At ratio near 1.0, phase locking is very accurate so we trust it fully.
            // As the ratio increases, IF estimates become more valuable for preserving
            // frequency accuracy, so we blend in up to 10% IF (reduced from 30% to improve coherence).
            let transient_focus = self.transient_focus_active();
            let if_blend = if transient_focus {
                0.0
            } else {
                (0.06 * ((hop_ratio - 1.0).abs() / 0.5).min(1.0)).min(0.06)
            };
            if if_blend > 1e-6 {
                for bin in self.sub_bass_bin..num_bins {
                    if self.peaks.binary_search(&bin).is_ok() {
                        continue; // Peak bins keep their locked phase
                    }
                    // Distance-adaptive IF: bins far from peaks get more IF blend
                    // since they benefit less from phase locking and more from
                    // independent frequency tracking.
                    let nearest_dist = match self.peaks.binary_search(&bin) {
                        Ok(_) => 0,
                        Err(idx) => {
                            let lower = if idx > 0 {
                                bin - self.peaks[idx - 1]
                            } else {
                                usize::MAX
                            };
                            let upper = if idx < self.peaks.len() {
                                self.peaks[idx] - bin
                            } else {
                                usize::MAX
                            };
                            lower.min(upper)
                        }
                    };
                    let dist_scale = if nearest_dist > 4 {
                        1.0 + ((nearest_dist - 4) as f64 / 10.0).min(3.0)
                    } else {
                        1.0
                    };
                    let bin_if_blend = (if_blend * dist_scale).min(0.20);
                    let locked = self.new_phases[bin] as f64;
                    let if_est = self.if_phases_backup[bin] as f64;
                    // Wrapped-difference blend (see the gradient blend above).
                    self.new_phases[bin] =
                        (locked + bin_if_blend * wrap_phase_f64(if_est - locked)) as f32;
                }
            }

            self.reconstruct_spectrum(num_bins);
            fft_inverse.process_with_scratch(&mut self.fft_buffer, &mut self.fft_inverse_scratch);

            // Fractional overlap-add: when synthesis frame starts between samples,
            // distribute each sample between nearest output samples via linear interpolation.
            if frac <= SYNTH_POS_EPSILON {
                for i in 0..self.fft_size {
                    let idx = synthesis_floor + i;
                    self.output_buf[idx] += self.fft_buffer[i].re * self.ola_gain[i];
                    self.window_sum_buf[idx] += self.ola_window_product[i];
                }
            } else {
                let w0 = 1.0 - frac;
                let w1 = frac;
                for i in 0..self.fft_size {
                    let idx = synthesis_floor + i;
                    let sample = self.fft_buffer[i].re as f64 * self.ola_gain_f64[i];
                    let window_weight = self.ola_window_product_f64[i];

                    self.output_buf[idx] += (sample * w0) as f32;
                    self.output_buf[idx + 1] += (sample * w1) as f32;
                    self.window_sum_buf[idx] += (window_weight * w0) as f32;
                    self.window_sum_buf[idx + 1] += (window_weight * w1) as f32;
                }
            }

            self.decay_transient_focus();

            synthesis_frame_pos = synthesis_pos + frame_advance;
        }

        let next_local_synthesis_pos = snap_near_integer(synthesis_frame_pos);
        self.synthesis_pos = self.synthesis_emitted as f64 + next_local_synthesis_pos;

        Ok((emit_len, output_len))
    }

    /// Windows the input frame and transforms to frequency domain.
    #[inline]
    fn analyze_frame(
        &mut self,
        input_frame: &[f32],
        fft_forward: &std::sync::Arc<dyn rustfft::Fft<f32>>,
    ) {
        let len = input_frame.len().min(self.fft_buffer.len());
        for (i, (&sample, &w)) in input_frame
            .iter()
            .zip(self.window.iter())
            .enumerate()
            .take(len)
        {
            self.fft_buffer[i] = Complex::new(sample * w, 0.0);
        }
        fft_forward.process_with_scratch(&mut self.fft_buffer, &mut self.fft_forward_scratch);
    }

    /// Extracts magnitudes and advances phase accumulators for each bin.
    ///
    /// Uses a multi-pass approach for improved frequency tracking and phase coherence:
    /// 1. Compute magnitudes and raw analysis phases for all bins.
    /// 2. Detect spectral peaks and compute refined instantaneous frequencies via
    ///    parabolic interpolation of the log-magnitude spectrum.
    /// 3. Advance phases using instantaneous frequency (IF) estimation: compute the
    ///    true per-hop phase advance of each bin from the phase difference, then
    ///    scale it by the f64 stretch ratio to resynthesize at the correct rate.
    ///    This handles the stretch ratio directly and eliminates cumulative
    ///    phase drift.
    /// 4. Apply soft phase gradient integration to propagate coherent phase from
    ///    peaks to nearby non-peak bins.
    ///
    /// Sub-bass bins (below `sub_bass_bin`) use rigid phase propagation to prevent
    /// phase cancellation in the critical sub-bass region and are excluded from
    /// peak-based refinements.
    ///
    /// Phase accumulation uses f64 precision to prevent cumulative rounding errors
    /// over long signals. The final phases are converted back to f32 for the
    /// spectrum reconstruction step.
    #[inline]
    fn advance_phases(&mut self, num_bins: usize, hop_ratio: f64) {
        let hop_a = self.hop_analysis as f64;
        let fft = self.fft_size as f64;

        // --- Pass 1: Extract magnitudes and analysis phases ---
        for bin in 0..num_bins {
            let c = self.fft_buffer[bin];
            self.magnitudes[bin] = c.norm();
            self.analysis_phases[bin] = c.arg();
        }

        // --- Pass 2: Detect peaks for IF refinement + phase gradient ---
        let search_start = self.sub_bass_bin.max(1);
        self.peaks.clear();
        if num_bins >= 3 && search_start < num_bins.saturating_sub(1) {
            for bin in search_start..num_bins - 1 {
                if self.magnitudes[bin] > MIN_PEAK_MAGNITUDE
                    && self.magnitudes[bin] > self.magnitudes[bin - 1]
                    && self.magnitudes[bin] > self.magnitudes[bin + 1]
                {
                    self.peaks.push(bin);
                }
            }
        }

        // --- Pass 3: Advance phases using instantaneous frequency (IF) estimation ---
        //
        // For each bin we compute the true per-hop phase advance from the phase
        // difference between consecutive frames, then advance the synthesis phase
        // accumulator by that advance multiplied by the f64 stretch ratio. Doing
        // the scaling as a single f64 multiplication eliminates the cumulative
        // drift that a separate div+mul through an integer synthesis hop causes.
        //
        // For spectral peak bins, parabolic interpolation of the log-magnitude
        // spectrum refines the frequency estimate to sub-bin precision (~1 Hz
        // accuracy vs ~5 Hz for integer-bin estimation).
        for bin in 0..num_bins {
            let phase = self.analysis_phases[bin] as f64;

            // Seed bins that were explicitly reset (full or selective reset)
            // from the current analysis phase to avoid a bogus first IF jump.
            if self.phase_seed_pending[bin] {
                self.phase_accum[bin] = phase;
                self.new_phases[bin] = phase as f32;
                self.prev_phase[bin] = phase;
                self.phase_seed_pending[bin] = false;
                continue;
            }

            // DC and Nyquist of a real signal are real-valued: hold the
            // analysis phase (0 or PI) instead of accumulating a synthetic
            // phase that would break Hermitian symmetry at reconstruction.
            if bin == 0 || bin == num_bins - 1 {
                self.phase_accum[bin] = phase;
                self.new_phases[bin] = phase as f32;
                self.prev_phase[bin] = phase;
                continue;
            }

            if bin < self.sub_bass_bin {
                // Sub-bass IF estimation: same instantaneous-frequency approach
                // as standard bins, but without parabolic interpolation (sub-bass
                // bins are narrow enough that integer-bin IF is sufficient).
                // The identity phase locking in phase_locking.rs handles inter-bin
                // coherence for sub-bass via trough-based regions.
                let expected_diff = self.expected_phase_advance[bin];
                let phase_diff = phase - self.prev_phase[bin];
                let deviation = wrap_phase_f64(phase_diff - expected_diff);
                // Keep the accumulator wrapped: it is only ever consumed via
                // `from_polar` and phase differences, and an unwrapped value
                // loses precision to the per-frame f32 downcast as it grows.
                self.phase_accum[bin] =
                    wrap_phase_f64(self.phase_accum[bin] + (expected_diff + deviation) * hop_ratio);
            } else {
                // Standard IF estimation:
                //   phase_diff = current_phase - prev_phase
                //   expected_diff = 2*pi * bin * hop_analysis / fft_size
                //   deviation = wrap(phase_diff - expected_diff)
                //   phase_accum += (expected_diff + deviation) * stretch_ratio
                let expected_diff = self.expected_phase_advance[bin]; // 2*pi*bin*hop_a/fft
                let phase_diff = phase - self.prev_phase[bin];
                let deviation = wrap_phase_f64(phase_diff - expected_diff);

                // For peak bins, use parabolic interpolation of log-magnitude
                // to refine the frequency estimate to sub-bin precision.
                //
                // The synthesis phase advance is `(expected + deviation)` scaled by
                // the f64 stretch ratio in a single multiplication. This minimizes
                // floating-point roundoff (especially at ratio 1.0, where the
                // advance must be bit-exact with the analysis advance).
                let is_peak = self.peaks.binary_search(&bin).is_ok();
                let phase_advance = if is_peak
                    && bin >= 1
                    && bin + 1 < num_bins
                    && self.magnitudes[bin] > MIN_PEAK_MAGNITUDE
                {
                    // Parabolic interpolation on log-magnitudes for sub-bin accuracy:
                    //   alpha = log(M[k-1]), beta = log(M[k]), gamma = log(M[k+1])
                    //   p = 0.5 * (alpha - gamma) / (alpha - 2*beta + gamma)
                    //   refined_freq_bin = k + p
                    //
                    // Log-magnitude interpolation gives better accuracy for Gaussian
                    // spectral peaks (which approximate windowed sinusoids) compared
                    // to linear interpolation.
                    let m_prev = (self.magnitudes[bin - 1] as f64).max(1e-30);
                    let m_curr = (self.magnitudes[bin] as f64).max(1e-30);
                    let m_next = (self.magnitudes[bin + 1] as f64).max(1e-30);
                    let alpha = m_prev.ln();
                    let beta = m_curr.ln();
                    let gamma = m_next.ln();
                    let denom = alpha - 2.0 * beta + gamma;
                    if denom.abs() > 1e-12 {
                        let p = 0.5 * (alpha - gamma) / denom;
                        // Refined expected phase advance based on interpolated bin position
                        let refined_expected = TWO_PI_F64 * (bin as f64 + p) * hop_a / fft;
                        let refined_deviation = wrap_phase_f64(phase_diff - refined_expected);
                        (refined_expected + refined_deviation) * hop_ratio
                    } else {
                        (expected_diff + deviation) * hop_ratio
                    }
                } else {
                    (expected_diff + deviation) * hop_ratio
                };

                self.phase_accum[bin] = wrap_phase_f64(self.phase_accum[bin] + phase_advance);
            }

            self.new_phases[bin] = self.phase_accum[bin] as f32;
            self.prev_phase[bin] = phase;
        }

        // --- Phase gradient integration (soft vertical coherence) ---
        // Propagate phase from peaks to nearby non-peak bins using the analysis
        // phase gradient, blended with the independently-advanced phase.
        // Apply up to 2.5x ratio with a tapering blend around unity:
        // full strength near ratio≈1.0, gradually reduced as we move away on
        // either side (compression or expansion) to avoid over-locking artifacts.
        if !self.transient_focus_active()
            && !self.peaks.is_empty()
            && (hop_ratio < 2.5 || self.wide_coherence_blend > 0.0)
        {
            let gradient_blend = if self.wide_coherence_blend > 0.0 {
                // Stage 11 experiment: hold this blend at wide ratios.
                self.wide_coherence_blend
            } else {
                let ratio_distance = (hop_ratio - 1.0).abs();
                // For ratios >1.0, taper gradient locking faster to avoid
                // over-locking smearing at larger slowdowns.
                let taper_span = if hop_ratio > 1.0 { 1.2 } else { 1.5 };
                PHASE_GRADIENT_BLEND * (1.0 - (ratio_distance / taper_span).clamp(0.0, 1.0))
            };
            for bin in self.sub_bass_bin..num_bins {
                if self.peaks.binary_search(&bin).is_ok() {
                    continue; // Peak bins keep their phase (they are the anchors)
                }

                // Find the nearest peak via binary search
                let nearest_peak = match self.peaks.binary_search(&bin) {
                    Ok(_) => unreachable!(),
                    Err(idx) => {
                        let lower = if idx > 0 {
                            Some(self.peaks[idx - 1])
                        } else {
                            None
                        };
                        let upper = if idx < self.peaks.len() {
                            Some(self.peaks[idx])
                        } else {
                            None
                        };
                        match (lower, upper) {
                            (Some(l), Some(u)) => {
                                if bin - l <= u - bin {
                                    l
                                } else {
                                    u
                                }
                            }
                            (Some(l), None) => l,
                            (None, Some(u)) => u,
                            (None, None) => continue,
                        }
                    }
                };

                // Attenuate gradient blend for bins far from their anchor peak.
                // Distant bins have less acoustic relationship to the peak,
                // so independent phase evolution is more appropriate.
                let peak_distance = bin.abs_diff(nearest_peak);
                let distance_fade = if peak_distance > 8 {
                    (1.0 - ((peak_distance - 8) as f64 / 24.0).min(1.0)).max(0.0)
                } else {
                    1.0
                };
                let effective_blend = gradient_blend * distance_fade;

                let gradient =
                    self.analysis_phases[bin] as f64 - self.analysis_phases[nearest_peak] as f64;
                let propagated = self.new_phases[nearest_peak] as f64 + gradient;
                let independent = self.new_phases[bin] as f64;
                // Blend along the wrapped difference: phases are angles, and
                // a linear mix of raw values that sit different multiples of
                // 2*PI apart lands on a meaningless intermediate angle.
                self.new_phases[bin] = (independent
                    + effective_blend * wrap_phase_f64(propagated - independent))
                    as f32;
            }
        }
    }

    /// Reconstructs the complex spectrum from magnitudes and phases,
    /// then mirrors negative frequencies for inverse FFT.
    #[inline]
    fn reconstruct_spectrum(&mut self, num_bins: usize) {
        for i in 0..num_bins {
            self.fft_buffer[i] = Complex::from_polar(self.magnitudes[i], self.new_phases[i]);
        }
        // DC and Nyquist of a real signal are real: project onto the real
        // axis so the mirrored spectrum below is exactly Hermitian and the
        // inverse FFT is real by construction (their held analysis phases
        // are 0 or PI, so this preserves signed magnitude).
        self.fft_buffer[0] = Complex::new(self.fft_buffer[0].re, 0.0);
        self.fft_buffer[num_bins - 1] = Complex::new(self.fft_buffer[num_bins - 1].re, 0.0);
        for bin in 1..num_bins - 1 {
            self.fft_buffer[self.fft_size - bin] = self.fft_buffer[bin].conj();
        }
    }

    /// Normalizes output by window sum, clamping to prevent amplification in
    /// low-overlap regions (occurs when synthesis hop > analysis hop).
    #[inline]
    fn normalize_output(output: &mut [f32], window_sum: &[f32], stretch_ratio: f64) {
        let max_window_sum = window_sum.iter().copied().fold(0.0f32, f32::max);
        Self::normalize_output_with_window_floor(output, window_sum, stretch_ratio, max_window_sum);
    }

    #[inline]
    fn normalize_output_with_window_floor(
        output: &mut [f32],
        window_sum: &[f32],
        stretch_ratio: f64,
        max_window_sum: f32,
    ) {
        // For stretches >1.0, synthesis frames are farther apart and aggressive
        // flooring can over-attenuate low-overlap regions. Relax the floor more
        // strongly with ratio while keeping a safety minimum against blow-ups.
        let floor_ratio = if stretch_ratio > 1.0 {
            (WINDOW_SUM_FLOOR_RATIO / (stretch_ratio as f32 * stretch_ratio as f32))
                .clamp(0.005, WINDOW_SUM_FLOOR_RATIO)
        } else {
            WINDOW_SUM_FLOOR_RATIO
        };
        let min_window_sum = (max_window_sum * floor_ratio).max(WINDOW_SUM_EPSILON);
        let len = output.len().min(window_sum.len());
        for i in 0..len {
            output[i] /= window_sum[i].max(min_window_sum);
        }
    }

    /// Chooses phase-locking mode for the current frame.
    ///
    /// When adaptive switching is enabled, this computes simple confidence
    /// features from the current magnitude spectrum:
    /// - spectral flatness (noise-likeness),
    /// - crest factor (peak prominence),
    /// - stretch-ratio distance from unity.
    ///
    /// Heuristic:
    /// - noisy/weak frames or large ratio offsets -> ROI
    /// - harmonic/peaky frames near unity ratio -> Identity
    /// - harmonic frames at moderate ratio offsets -> Selective
    /// - otherwise -> configured `phase_locking_mode`
    #[inline]
    fn select_phase_locking_mode(&self, num_bins: usize) -> PhaseLockingMode {
        if self.transient_focus_active() {
            return PhaseLockingMode::Identity;
        }

        if !self.adaptive_phase_locking {
            return self.phase_locking_mode;
        }
        if num_bins == 0 || self.sub_bass_bin >= num_bins {
            return self.phase_locking_mode;
        }

        let start = self.sub_bass_bin;
        let mags = &self.magnitudes[start..num_bins];
        if mags.is_empty() {
            return self.phase_locking_mode;
        }

        let mut max_mag = 0.0f64;
        let mut sum = 0.0f64;
        let mut log_sum = 0.0f64;
        for &m in mags {
            let mf = (m as f64).max(ADAPTIVE_FEATURE_EPS);
            max_mag = max_mag.max(mf);
            sum += mf;
            log_sum += mf.ln();
        }

        let n = mags.len() as f64;
        let mean = sum / n;
        let geometric = (log_sum / n).exp();
        let flatness = geometric / mean.max(ADAPTIVE_FEATURE_EPS);
        let crest = max_mag / mean.max(ADAPTIVE_FEATURE_EPS);
        let ratio_distance = (self.stretch_ratio - 1.0).abs();
        let harmonic_confidence = crest / (1.0 + 4.0 * flatness);

        if ratio_distance >= ADAPTIVE_FORCE_ROI_RATIO_DISTANCE
            || flatness >= ADAPTIVE_NOISY_FLATNESS
            || crest <= ADAPTIVE_NOISY_CREST
        {
            return PhaseLockingMode::RegionOfInfluence;
        }

        if ratio_distance <= ADAPTIVE_IDENTITY_RATIO_DISTANCE_MAX
            && harmonic_confidence >= ADAPTIVE_HARMONIC_CONFIDENCE
        {
            return PhaseLockingMode::Identity;
        }

        if ratio_distance <= ADAPTIVE_SELECTIVE_RATIO_DISTANCE_MAX
            && harmonic_confidence >= ADAPTIVE_SELECTIVE_HARMONIC_CONFIDENCE
        {
            return PhaseLockingMode::Selective;
        }

        self.phase_locking_mode
    }

    #[inline]
    fn transient_focus_active(&self) -> bool {
        self.transient_focus_frames > 0
    }

    #[inline]
    fn decay_transient_focus(&mut self) {
        self.transient_focus_frames = self.transient_focus_frames.saturating_sub(1);
    }
}

impl std::fmt::Debug for PhaseVocoder {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("PhaseVocoder")
            .field("fft_size", &self.fft_size)
            .field("hop_analysis", &self.hop_analysis)
            .field("stretch_ratio", &self.stretch_ratio)
            .field("synthesis_pos", &self.synthesis_pos)
            .field("sub_bass_bin", &self.sub_bass_bin)
            .field("phase_locking_mode", &self.phase_locking_mode)
            .field("adaptive_phase_locking", &self.adaptive_phase_locking)
            .field("streaming_tail_len", &self.streaming_tail.len())
            .finish()
    }
}

/// Snaps values extremely close to integer grid points to the exact integer.
#[inline]
fn snap_near_integer(value: f64) -> f64 {
    let rounded = value.round();
    if (value - rounded).abs() <= SYNTH_POS_EPSILON {
        rounded
    } else {
        value
    }
}

/// Wraps a phase value to [-PI, PI] using f64 precision.
#[inline]
fn wrap_phase_f64(phase: f64) -> f64 {
    let pi = std::f64::consts::PI;
    let p = phase + pi;
    p - (p / TWO_PI_F64).floor() * TWO_PI_F64 - pi
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::stretch::phase_locking::apply_phase_locking;
    use std::f32::consts::PI;

    const TWO_PI: f32 = 2.0 * PI;

    /// Wraps a phase value to [-PI, PI] using efficient modulo arithmetic (f32).
    fn wrap_phase(phase: f32) -> f32 {
        let p = phase + PI;
        p - (p / TWO_PI).floor() * TWO_PI - PI
    }

    #[test]
    fn test_wrap_phase() {
        assert!((wrap_phase(0.0) - 0.0).abs() < 1e-6);
        assert!((wrap_phase(PI + 0.1) - (-PI + 0.1)).abs() < 1e-5);
        assert!((wrap_phase(-PI - 0.1) - (PI - 0.1)).abs() < 1e-5);
        // Test larger values
        assert!((wrap_phase(10.0 * PI + 0.5) - wrap_phase(0.5)).abs() < 1e-4);
        assert!((wrap_phase(-10.0 * PI - 0.5) - wrap_phase(-0.5)).abs() < 1e-4);
    }

    #[test]
    fn test_normalize_output_with_window_floor_scopes_expansion_floor_to_seam_prefix() {
        let reference_max_window_sum = 1.0f32;
        let window_sum = [1.0f32, 0.01, 0.01, 0.01];
        let mut split_output = vec![1.0f32; window_sum.len()];
        let mut all_expansion_output = vec![1.0f32; window_sum.len()];

        PhaseVocoder::normalize_output_with_window_floor(
            &mut split_output[..2],
            &window_sum[..2],
            1.18,
            reference_max_window_sum,
        );
        PhaseVocoder::normalize_output_with_window_floor(
            &mut split_output[2..],
            &window_sum[2..],
            0.82,
            reference_max_window_sum,
        );
        PhaseVocoder::normalize_output_with_window_floor(
            &mut all_expansion_output,
            &window_sum,
            1.18,
            reference_max_window_sum,
        );

        assert!(
            (split_output[0] - all_expansion_output[0]).abs() < 1e-6
                && (split_output[1] - all_expansion_output[1]).abs() < 1e-6,
            "the unresolved seam prefix should keep the prior expansion ratio floor"
        );
        assert!(
            split_output[2] < all_expansion_output[2] && split_output[3] < all_expansion_output[3],
            "once the seam prefix drains, the remainder should normalize against the current ratio instead of the prior expansion floor"
        );
    }

    #[test]
    fn test_adaptive_phase_locking_prefers_roi_for_flat_spectrum() {
        let mut pv = PhaseVocoder::with_options(
            1024,
            256,
            1.0,
            44100,
            0.0,
            WindowType::Hann,
            PhaseLockingMode::Identity,
        );
        pv.set_adaptive_phase_locking(true);
        pv.magnitudes.fill(1.0);
        let mode = pv.select_phase_locking_mode(1024 / 2 + 1);
        assert_eq!(mode, PhaseLockingMode::RegionOfInfluence);
    }

    #[test]
    fn test_adaptive_phase_locking_prefers_identity_for_harmonic_frame() {
        let mut pv = PhaseVocoder::with_options(
            1024,
            256,
            1.05,
            44100,
            0.0,
            WindowType::Hann,
            PhaseLockingMode::RegionOfInfluence,
        );
        pv.set_adaptive_phase_locking(true);
        pv.magnitudes.fill(0.001);
        pv.magnitudes[20] = 1.0;
        pv.magnitudes[60] = 0.8;
        pv.magnitudes[120] = 0.6;
        let mode = pv.select_phase_locking_mode(1024 / 2 + 1);
        assert_eq!(mode, PhaseLockingMode::Identity);
    }

    #[test]
    fn test_adaptive_phase_locking_prefers_selective_for_moderate_ratio_harmonic_frame() {
        let mut pv = PhaseVocoder::with_options(
            1024,
            256,
            1.45,
            44100,
            0.0,
            WindowType::Hann,
            PhaseLockingMode::RegionOfInfluence,
        );
        pv.set_adaptive_phase_locking(true);
        pv.magnitudes.fill(0.001);
        pv.magnitudes[20] = 1.0;
        pv.magnitudes[60] = 0.8;
        pv.magnitudes[120] = 0.6;
        let mode = pv.select_phase_locking_mode(1024 / 2 + 1);
        assert_eq!(mode, PhaseLockingMode::Selective);
    }

    #[test]
    fn test_adaptive_phase_locking_disabled_uses_configured_mode() {
        let mut pv = PhaseVocoder::with_options(
            1024,
            256,
            1.0,
            44100,
            0.0,
            WindowType::Hann,
            PhaseLockingMode::Identity,
        );
        pv.set_adaptive_phase_locking(false);
        pv.magnitudes.fill(1.0);
        let mode = pv.select_phase_locking_mode(1024 / 2 + 1);
        assert_eq!(mode, PhaseLockingMode::Identity);
    }

    #[test]
    fn test_reset_phase_state_bands_enables_transient_focus_for_audible_bands() {
        let sample_rate = 44_100u32;
        let mut pv = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 120.0);
        assert_eq!(pv.transient_focus_frames, 0);

        pv.reset_phase_state_bands([false, false, true, false], sample_rate);
        assert_eq!(
            pv.transient_focus_frames, TRANSIENT_FOCUS_FRAMES,
            "mid-band reset should engage transient focus"
        );
    }

    #[test]
    fn test_reset_phase_state_bands_sub_only_does_not_enable_transient_focus() {
        let sample_rate = 44_100u32;
        let mut pv = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 120.0);
        assert_eq!(pv.transient_focus_frames, 0);

        pv.reset_phase_state_bands([true, false, false, false], sample_rate);
        assert_eq!(
            pv.transient_focus_frames, 0,
            "sub-only reset should not engage transient focus"
        );
    }

    #[test]
    fn test_select_phase_locking_mode_forces_identity_during_transient_focus() {
        let mut pv = PhaseVocoder::with_options(
            1024,
            256,
            1.4,
            44_100,
            120.0,
            WindowType::Hann,
            PhaseLockingMode::RegionOfInfluence,
        );
        pv.transient_focus_frames = 1;
        pv.set_adaptive_phase_locking(true);
        pv.magnitudes.fill(1.0);

        let mode = pv.select_phase_locking_mode(1024 / 2 + 1);
        assert_eq!(
            mode,
            PhaseLockingMode::Identity,
            "transient focus should force identity locking"
        );
    }

    #[test]
    fn test_reset_phase_state_bands_marks_only_target_band_for_seeding() {
        let sample_rate = 44_100u32;
        let mut pv = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 120.0);
        pv.phase_accum.fill(1.0);
        pv.prev_phase.fill(1.0);
        pv.phase_seed_pending.fill(false);

        // Reset only the low band [100, 500) Hz.
        pv.reset_phase_state_bands([false, true, false, false], sample_rate);

        let bin_hz = sample_rate as f32 / pv.fft_size as f32;
        for bin in 0..pv.phase_accum.len() {
            let freq = bin as f32 * bin_hz;
            let in_low_band = (100.0..500.0).contains(&freq);
            if in_low_band {
                assert_eq!(pv.phase_accum[bin], 0.0, "low-band phase_accum not reset");
                assert_eq!(pv.prev_phase[bin], 0.0, "low-band prev_phase not reset");
                assert!(
                    pv.phase_seed_pending[bin],
                    "low-band seed flag should be set"
                );
            } else {
                assert_eq!(pv.phase_accum[bin], 1.0, "non-low-band phase_accum changed");
                assert_eq!(pv.prev_phase[bin], 1.0, "non-low-band prev_phase changed");
                assert!(
                    !pv.phase_seed_pending[bin],
                    "non-low-band seed flag should stay clear"
                );
            }
        }
    }

    #[test]
    fn test_advance_phases_seeds_pending_bins_without_if_jump() {
        let sample_rate = 44_100u32;
        let mut pv = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 120.0);
        let num_bins = pv.fft_size / 2 + 1;

        // Simulate steady state where no bins are pending.
        pv.phase_seed_pending.fill(false);
        pv.phase_accum.fill(0.5);
        pv.prev_phase.fill(0.5);

        let target_bin = 64usize;
        pv.phase_accum[target_bin] = 0.0;
        pv.prev_phase[target_bin] = 0.0;
        pv.phase_seed_pending[target_bin] = true;

        for bin in 0..num_bins {
            let phase = 0.2 + bin as f32 * 0.001;
            pv.fft_buffer[bin] = Complex::from_polar(1.0, phase);
        }

        pv.advance_phases(num_bins, 1.0);

        let seeded_phase = pv.analysis_phases[target_bin] as f64;
        assert!(
            (pv.phase_accum[target_bin] - seeded_phase).abs() < 1e-9,
            "pending bin should seed directly from analysis phase"
        );
        assert!(
            (pv.prev_phase[target_bin] - seeded_phase).abs() < 1e-9,
            "pending bin prev phase should match seeded analysis phase"
        );
        assert!(
            !pv.phase_seed_pending[target_bin],
            "pending flag should clear after first seeded frame"
        );
    }

    #[test]
    fn test_advance_phases_skips_phase_gradient_during_transient_focus() {
        let sample_rate = 44_100u32;
        let num_bins = 1024 / 2 + 1;
        let peak_bin = 20usize;
        let target_bin = 21usize;

        let configure = |pv: &mut PhaseVocoder| {
            pv.phase_seed_pending.fill(false);
            pv.phase_accum.fill(0.0);

            for bin in 0..num_bins {
                pv.prev_phase[bin] = -pv.expected_phase_advance[bin];
                let magnitude = if bin == peak_bin {
                    2.0
                } else if bin == peak_bin - 1 || bin == peak_bin + 1 {
                    0.5
                } else {
                    0.1
                };
                pv.fft_buffer[bin] = Complex::from_polar(magnitude, 0.0);
            }
        };

        let mut focused = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 0.0);
        configure(&mut focused);
        focused.transient_focus_frames = 1;
        focused.advance_phases(num_bins, 1.0);

        let mut unfocused = PhaseVocoder::new(1024, 256, 1.0, sample_rate, 0.0);
        configure(&mut unfocused);
        unfocused.advance_phases(num_bins, 1.0);

        // The accumulator stays wrapped, so phases compare modulo 2*PI.
        let independent_phase = wrap_phase_f64(focused.expected_phase_advance[target_bin]) as f32;
        assert!(
            wrap_phase_f64((focused.new_phases[target_bin] - independent_phase) as f64).abs()
                < 1e-6,
            "transient focus should keep non-peak bins on their independently advanced phase instead of reapplying the gradient field"
        );
        assert!(
            wrap_phase_f64((unfocused.new_phases[target_bin] - independent_phase) as f64).abs()
                > 1e-3,
            "without transient focus the same non-peak bin should still be pulled by phase-gradient integration"
        );
    }

    #[test]
    fn test_phase_vocoder_identity() {
        // Stretch ratio 1.0 should approximately preserve the signal
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;

        // Generate a 440 Hz sine wave
        let num_samples = fft_size * 4;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);
        let output = pv.process(&input).unwrap();

        // Output length should be approximately the same
        let len_ratio = output.len() as f64 / input.len() as f64;
        assert!(
            (len_ratio - 1.0).abs() < 0.1,
            "Length ratio {} too far from 1.0",
            len_ratio
        );

        // Check that the output contains a similar frequency
        // (RMS should be similar)
        let input_rms: f32 = (input.iter().map(|x| x * x).sum::<f32>() / input.len() as f32).sqrt();
        let output_rms: f32 =
            (output.iter().map(|x| x * x).sum::<f32>() / output.len() as f32).sqrt();

        assert!(
            (output_rms - input_rms).abs() < input_rms * 0.5,
            "RMS mismatch: input={}, output={}",
            input_rms,
            output_rms
        );
    }

    #[test]
    fn test_phase_vocoder_stretch() {
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;

        // Use a longer signal for more accurate length ratio
        let num_samples = fft_size * 8;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let stretch_ratio = 2.0;
        let mut pv = PhaseVocoder::new(fft_size, hop, stretch_ratio, sample_rate, 120.0);
        let output = pv.process(&input).unwrap();

        // Output should be approximately 2x longer (with tolerance for edge effects)
        let len_ratio = output.len() as f64 / input.len() as f64;
        assert!(
            (len_ratio - stretch_ratio).abs() < 0.35,
            "Length ratio {} too far from {}",
            len_ratio,
            stretch_ratio
        );
    }

    #[test]
    fn test_phase_vocoder_compress() {
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;

        let num_samples = fft_size * 4;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let stretch_ratio = 0.5;
        let mut pv = PhaseVocoder::new(fft_size, hop, stretch_ratio, sample_rate, 120.0);
        let output = pv.process(&input).unwrap();

        let len_ratio = output.len() as f64 / input.len() as f64;
        assert!(
            (len_ratio - stretch_ratio).abs() < 0.2,
            "Length ratio {} too far from {}",
            len_ratio,
            stretch_ratio
        );
    }

    #[test]
    fn test_phase_vocoder_input_too_short() {
        let mut pv = PhaseVocoder::new(4096, 1024, 1.0, 44100, 120.0);
        let result = pv.process(&[0.0; 100]);
        assert!(result.is_err());
    }

    #[test]
    fn test_sub_bass_bin_calculation() {
        // 120 Hz cutoff at 44100 Hz with FFT size 4096
        // Expected bin: 120 * 4096 / 44100 ≈ 11.15 → 11
        let pv = PhaseVocoder::new(4096, 1024, 1.0, 44100, 120.0);
        assert_eq!(pv.sub_bass_bin, 11);

        // 0 Hz cutoff should give bin 0 (no sub-bass locking)
        let pv = PhaseVocoder::new(4096, 1024, 1.0, 44100, 0.0);
        assert_eq!(pv.sub_bass_bin, 0);

        // High cutoff at 48000 Hz
        let pv = PhaseVocoder::new(4096, 1024, 1.0, 48000, 200.0);
        let expected = (200.0f32 * 4096.0 / 48000.0).round() as usize;
        assert_eq!(pv.sub_bass_bin, expected);
    }

    #[test]
    fn test_sub_bass_phase_locking_preserves_low_freq() {
        // A 60 Hz sine should be handled by sub-bass rigid phase locking.
        // Compare output quality with sub-bass locking (120 Hz cutoff)
        // vs without (0 Hz cutoff).
        let sample_rate = 44100u32;
        let fft_size = 4096;
        let hop = fft_size / 4;
        let num_samples = fft_size * 8;
        let freq = 60.0f32; // Well below 120 Hz cutoff

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * freq * i as f32 / sample_rate as f32).sin())
            .collect();

        // Process with sub-bass locking enabled (120 Hz cutoff)
        let mut pv_locked = PhaseVocoder::new(fft_size, hop, 1.5, sample_rate, 120.0);
        let output_locked = pv_locked.process(&input).unwrap();

        // Process without sub-bass locking (0 Hz cutoff)
        let mut pv_unlocked = PhaseVocoder::new(fft_size, hop, 1.5, sample_rate, 0.0);
        let output_unlocked = pv_unlocked.process(&input).unwrap();

        // Both should produce output
        assert!(!output_locked.is_empty());
        assert!(!output_unlocked.is_empty());

        // Both should have similar RMS (we aren't destroying energy)
        let rms_locked =
            (output_locked.iter().map(|x| x * x).sum::<f32>() / output_locked.len() as f32).sqrt();
        let rms_unlocked = (output_unlocked.iter().map(|x| x * x).sum::<f32>()
            / output_unlocked.len() as f32)
            .sqrt();

        assert!(
            rms_locked > 0.1,
            "Sub-bass locked output should have significant energy, got RMS={}",
            rms_locked
        );
        assert!(
            rms_unlocked > 0.1,
            "Unlocked output should have significant energy, got RMS={}",
            rms_unlocked
        );
    }

    #[test]
    fn test_sub_bass_locking_does_not_affect_high_freq() {
        // A 1000 Hz sine should NOT be affected by sub-bass phase locking
        // (it's above the 120 Hz cutoff).
        let sample_rate = 44100u32;
        let fft_size = 4096;
        let hop = fft_size / 4;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 1000.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv_with = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);
        let output_with = pv_with.process(&input).unwrap();

        let mut pv_without = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 0.0);
        let output_without = pv_without.process(&input).unwrap();

        // Output lengths should be the same
        assert_eq!(output_with.len(), output_without.len());

        // RMS should be very similar since 1000 Hz is above the cutoff
        let rms_with =
            (output_with.iter().map(|x| x * x).sum::<f32>() / output_with.len() as f32).sqrt();
        let rms_without = (output_without.iter().map(|x| x * x).sum::<f32>()
            / output_without.len() as f32)
            .sqrt();

        assert!(
            (rms_with - rms_without).abs() < rms_with * 0.3,
            "1000 Hz signal should be similar with/without sub-bass locking: {} vs {}",
            rms_with,
            rms_without
        );
    }

    #[test]
    fn test_phase_vocoder_with_blackman_harris() {
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::with_window(
            fft_size,
            hop,
            1.5,
            sample_rate,
            120.0,
            WindowType::BlackmanHarris,
        );
        let output = pv.process(&input).unwrap();

        // Should produce valid stretched output
        assert!(!output.is_empty());
        let len_ratio = output.len() as f64 / input.len() as f64;
        assert!(
            (len_ratio - 1.5).abs() < 0.3,
            "BH window length ratio {} too far from 1.5",
            len_ratio
        );
    }

    #[test]
    fn test_phase_vocoder_with_kaiser() {
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::with_window(
            fft_size,
            hop,
            1.5,
            sample_rate,
            120.0,
            WindowType::Kaiser(800),
        );
        let output = pv.process(&input).unwrap();

        assert!(!output.is_empty());
        let len_ratio = output.len() as f64 / input.len() as f64;
        assert!(
            (len_ratio - 1.5).abs() < 0.3,
            "Kaiser window length ratio {} too far from 1.5",
            len_ratio
        );
    }

    #[test]
    fn test_phase_vocoder_different_windows_produce_different_output() {
        let sample_rate = 44100;
        let fft_size = 4096;
        let hop = fft_size / 4;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv_hann =
            PhaseVocoder::with_window(fft_size, hop, 1.5, sample_rate, 120.0, WindowType::Hann);
        let output_hann = pv_hann.process(&input).unwrap();

        let mut pv_bh = PhaseVocoder::with_window(
            fft_size,
            hop,
            1.5,
            sample_rate,
            120.0,
            WindowType::BlackmanHarris,
        );
        let output_bh = pv_bh.process(&input).unwrap();

        // Both should produce valid output of similar length
        assert!(!output_hann.is_empty());
        assert!(!output_bh.is_empty());

        // Outputs should differ (different windows produce different spectral characteristics)
        let min_len = output_hann.len().min(output_bh.len());
        let diff: f32 = output_hann[..min_len]
            .iter()
            .zip(&output_bh[..min_len])
            .map(|(a, b)| (a - b).abs())
            .sum::<f32>()
            / min_len as f32;
        assert!(
            diff > 1e-6,
            "Different windows should produce different output, avg diff = {}",
            diff
        );
    }

    // --- phase locking integration (detailed tests in phase_locking module) ---

    #[test]
    fn test_phase_lock_identity_no_peaks() {
        // Flat magnitude spectrum: no local maxima → no peaks → phases unchanged
        let num_bins = 16;
        let magnitudes = vec![1.0f32; num_bins]; // all equal, no peaks
        let analysis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.1).collect();
        let mut synthesis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.2).collect();
        let original_phases = synthesis_phases.clone();
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            num_bins,
            0,
            &mut peaks,
        );

        // With no peaks found, phases should remain unchanged
        assert_eq!(synthesis_phases, original_phases);
    }

    #[test]
    fn test_phase_lock_identity_single_peak() {
        // Single peak at bin 5 with a realistic spectral lobe shape.
        // Trough-bounded identity locking propagates the peak's phase
        // rotation to all bins within its influence region (between troughs).
        let num_bins = 16;
        // Create a Gaussian-like lobe centered at bin 5, with troughs at 0 and 15
        let magnitudes: Vec<f32> = (0..num_bins)
            .map(|i| {
                let dist = (i as f32 - 5.0).abs();
                0.01 + 0.99 * (-dist * dist / 8.0).exp()
            })
            .collect();
        let analysis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.3).collect();
        let mut synthesis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.5).collect();
        let peak_synth = synthesis_phases[5];
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            num_bins,
            0, // start_bin = 0
            &mut peaks,
        );

        // Peak at bin 5 should keep its phase
        assert!(
            (synthesis_phases[5] - peak_synth).abs() < 1e-6,
            "Peak bin should keep its phase"
        );

        // The phase rotation from the peak
        let phase_rotation = peak_synth - analysis_phases[5];

        // Bins in the peak's influence region should have:
        // synth[bin] = analysis[bin] + phase_rotation
        // With a single Gaussian lobe and no other peaks, all bins should
        // be in the peak's influence region.
        for bin in 1..num_bins - 1 {
            if bin == 5 {
                continue;
            }
            let expected = analysis_phases[bin] + phase_rotation;
            assert!(
                (synthesis_phases[bin] - expected).abs() < 1e-5,
                "Bin {} should be locked to peak: got {}, expected {}",
                bin,
                synthesis_phases[bin],
                expected
            );
        }
    }

    #[test]
    fn test_phase_lock_start_bin_above_num_bins() {
        // start_bin >= num_bins: early return, no changes
        let num_bins = 8;
        let magnitudes = vec![0.0f32; num_bins];
        let analysis_phases = vec![0.0f32; num_bins];
        let mut synthesis_phases = vec![1.0f32; num_bins];
        let original = synthesis_phases.clone();
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            num_bins,
            num_bins, // start_bin == num_bins → early return
            &mut peaks,
        );

        assert_eq!(synthesis_phases, original);
    }

    #[test]
    fn test_phase_lock_num_bins_less_than_3() {
        // num_bins < 3: early return
        let magnitudes = vec![1.0f32; 2];
        let analysis_phases = vec![0.0f32; 2];
        let mut synthesis_phases = vec![0.5f32; 2];
        let original = synthesis_phases.clone();
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            2,
            0,
            &mut peaks,
        );

        assert_eq!(synthesis_phases, original);
    }

    #[test]
    fn test_phase_lock_sub_bass_region_skipped() {
        // Peaks exist only below start_bin → no peaks found above sub-bass
        let num_bins = 16;
        let mut magnitudes = vec![0.1f32; num_bins];
        magnitudes[2] = 1.0; // peak below start_bin=5
        let analysis_phases = vec![0.0f32; num_bins];
        let mut synthesis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32).collect();
        let original = synthesis_phases.clone();
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            num_bins,
            5, // start_bin=5, peak at bin 2 is below
            &mut peaks,
        );

        // No peaks above start_bin → no changes
        assert_eq!(synthesis_phases, original);
    }

    #[test]
    fn test_phase_lock_multiple_peaks() {
        // Two peaks with realistic spectral lobe shapes.
        // Trough-bounded identity locking assigns each bin to the peak
        // whose influence region (bounded by troughs) contains it.
        let num_bins = 16;
        // Create two Gaussian lobes: peak at bin 3, peak at bin 10
        // with a clear trough between them (around bin 7)
        let magnitudes: Vec<f32> = (0..num_bins)
            .map(|i| {
                let d3 = (i as f32 - 3.0).abs();
                let d10 = (i as f32 - 10.0).abs();
                let lobe3 = 1.0 * (-d3 * d3 / 4.0).exp();
                let lobe10 = 0.8 * (-d10 * d10 / 4.0).exp();
                0.001 + lobe3.max(lobe10) // ensure non-zero floor
            })
            .collect();
        let analysis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.1).collect();
        let mut synthesis_phases: Vec<f32> = (0..num_bins).map(|i| i as f32 * 0.2).collect();
        let synth_peak3 = synthesis_phases[3];
        let synth_peak10 = synthesis_phases[10];
        let mut peaks = Vec::new();

        apply_phase_locking(
            PhaseLockingMode::Identity,
            &magnitudes,
            &analysis_phases,
            &mut synthesis_phases,
            num_bins,
            1, // start_bin=1
            &mut peaks,
        );

        // Verify both peaks are found
        assert!(peaks.contains(&3), "Should find peak at bin 3");
        assert!(peaks.contains(&10), "Should find peak at bin 10");

        // Phase rotation for each peak
        let rotation_3 = synth_peak3 - analysis_phases[3];
        let rotation_10 = synth_peak10 - analysis_phases[10];

        // Bin 2 is in peak 3's influence region (between start_bin boundary and trough)
        let expected_2 = analysis_phases[2] + rotation_3;
        assert!(
            (synthesis_phases[2] - expected_2).abs() < 1e-5,
            "Bin 2 should lock to peak 3: got {}, expected {}",
            synthesis_phases[2],
            expected_2
        );

        // Bin 12 is in peak 10's influence region
        let expected_12 = analysis_phases[12] + rotation_10;
        assert!(
            (synthesis_phases[12] - expected_12).abs() < 1e-5,
            "Bin 12 should lock to peak 10: got {}, expected {}",
            synthesis_phases[12],
            expected_12
        );
    }

    // --- normalize_output internals ---

    #[test]
    fn test_normalize_output_uniform_window_sum() {
        // When window_sum is uniform, output should be divided by that value
        let mut output = vec![2.0f32; 10];
        let window_sum = vec![2.0f32; 10];
        PhaseVocoder::normalize_output(&mut output, &window_sum, 1.0);
        for &s in &output {
            assert!((s - 1.0).abs() < 1e-6, "Expected 1.0, got {}", s);
        }
    }

    #[test]
    fn test_normalize_output_low_window_sum_clamped() {
        // Very small window sums should be clamped to min_window_sum
        // to prevent amplification
        let mut output = vec![1.0f32; 10];
        let mut window_sum = vec![1.0f32; 10];
        // One sample has near-zero window sum (low-overlap region)
        window_sum[5] = 1e-10;
        PhaseVocoder::normalize_output(&mut output, &window_sum, 1.0);

        // The clamped sample should NOT be amplified wildly
        // min_window_sum = max(1.0) * WINDOW_SUM_FLOOR_RATIO = 0.1
        // So output[5] = 1.0 / 0.1 = 10.0
        assert!(
            output[5] <= 11.0,
            "Low window sum should be clamped, got {}",
            output[5]
        );
        // Normal samples should be ~1.0
        assert!((output[0] - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_normalize_output_all_zero_window_sum() {
        // All-zero window sum: should use WINDOW_SUM_EPSILON floor
        let mut output = vec![1.0f32; 5];
        let window_sum = vec![0.0f32; 5];
        PhaseVocoder::normalize_output(&mut output, &window_sum, 1.0);
        // Each sample = 1.0 / WINDOW_SUM_EPSILON
        for &s in &output {
            assert!(s.is_finite(), "Output should be finite, got {}", s);
        }
    }

    // --- wrap_phase edge cases ---

    #[test]
    fn test_wrap_phase_exact_boundaries() {
        // Exactly PI should wrap to -PI (or very close)
        let result = wrap_phase(PI);
        assert!(
            (result - (-PI)).abs() < 1e-5 || (result - PI).abs() < 1e-5,
            "wrap_phase(PI) = {} should be near ±PI",
            result
        );

        // Exactly -PI
        let result = wrap_phase(-PI);
        assert!(
            (result - (-PI)).abs() < 1e-5 || (result - PI).abs() < 1e-5,
            "wrap_phase(-PI) = {} should be near ±PI",
            result
        );

        // Exactly 0
        assert!((wrap_phase(0.0)).abs() < 1e-6);
    }

    #[test]
    fn test_wrap_phase_very_large_values() {
        // Very large positive and negative values
        let result = wrap_phase(1000.0 * PI);
        assert!(
            (-PI..=PI).contains(&result),
            "wrap_phase(1000*PI) = {} should be in [-PI, PI]",
            result
        );

        let result = wrap_phase(-999.0 * PI);
        assert!(
            (-PI..=PI).contains(&result),
            "wrap_phase(-999*PI) = {} should be in [-PI, PI]",
            result
        );
    }

    // --- set_stretch_ratio ---

    #[test]
    fn test_set_stretch_ratio_preserves_phase_state() {
        // Process some audio, then change ratio and process more.
        // Phase should be continuous (no reset).
        let fft_size = 4096;
        let hop = 1024;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);
        let output1 = pv.process(&input).unwrap();
        assert!(!output1.is_empty());

        // Change ratio and process again — should work without error
        pv.set_stretch_ratio(1.5);
        let output2 = pv.process(&input).unwrap();
        assert!(!output2.is_empty());
        assert!(output2.len() > output1.len()); // 1.5x should be longer
    }

    #[test]
    fn test_process_streaming_and_flush_produce_finite_output() {
        let fft_size = 2048;
        let hop = 512;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 10;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                (2.0 * PI * 220.0 * t).sin() * 0.6 + (2.0 * PI * 880.0 * t).sin() * 0.25
            })
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.1, sample_rate, 120.0);
        let mut total_output = Vec::new();
        let mut analysis_buffer = Vec::new();

        for chunk in input.chunks(700) {
            analysis_buffer.extend_from_slice(chunk);
            if analysis_buffer.len() < fft_size {
                continue;
            }

            let out = pv.process_streaming(&analysis_buffer).unwrap();
            total_output.extend_from_slice(&out);

            let num_frames = (analysis_buffer.len() - fft_size) / hop + 1;
            let consumed = num_frames * hop;
            analysis_buffer.drain(..consumed);
        }

        if !analysis_buffer.is_empty() {
            analysis_buffer.resize(fft_size, 0.0);
            let out = pv.process_streaming(&analysis_buffer).unwrap();
            total_output.extend_from_slice(&out);
        }

        let tail = pv.flush_streaming().unwrap();
        total_output.extend_from_slice(&tail);

        assert!(
            !total_output.is_empty(),
            "Streaming path should produce output"
        );
        assert!(total_output.iter().all(|s| s.is_finite()));
    }

    #[test]
    fn test_streaming_tail_ratio_preserves_overlap_history_for_large_ratio_change() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 5;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                (2.0 * PI * 220.0 * t).sin() * 0.55 + (2.0 * PI * 660.0 * t).sin() * 0.20
            })
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.18, sample_rate, 120.0);
        let first_chunk_len = fft_size * 2;
        let first = pv.process_streaming(&input[..first_chunk_len]).unwrap();
        assert!(!first.is_empty(), "first streaming chunk should emit audio");
        assert!(
            !pv.streaming_tail.is_empty(),
            "first streaming chunk should retain overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.18).abs() < 1e-12,
            "tail ratio should match the ratio that generated the carried overlap"
        );

        let first_frames = (first_chunk_len - fft_size) / hop + 1;
        let consumed = first_frames * hop;
        pv.set_stretch_ratio(0.82);
        let second = pv
            .process_streaming(&input[consumed..consumed + first_chunk_len])
            .unwrap();
        assert!(
            !second.is_empty(),
            "second streaming chunk should emit audio"
        );
        assert!(
            !pv.streaming_tail.is_empty(),
            "second streaming chunk should continue carrying overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 0.82).abs() < 1e-12,
            "once the previous overlap has fully emitted, the carried tail should re-arm to the current ratio"
        );

        let _ = pv.flush_streaming().unwrap();
        assert!(
            (pv.streaming_tail_ratio - 0.82).abs() < 1e-12,
            "flush should clear carried overlap history and re-arm the current ratio"
        );
    }

    #[test]
    fn test_streaming_tail_ratio_preserves_carried_overlap_for_small_cross_unity_modulation() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 5;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                (2.0 * PI * 220.0 * t).sin() * 0.55 + (2.0 * PI * 660.0 * t).sin() * 0.20
            })
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.04, sample_rate, 120.0);
        let first_chunk_len = fft_size * 2;
        let first = pv.process_streaming(&input[..first_chunk_len]).unwrap();
        assert!(!first.is_empty(), "first streaming chunk should emit audio");
        assert!(
            !pv.streaming_tail.is_empty(),
            "first streaming chunk should retain overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.04).abs() < 1e-12,
            "tail ratio should match the ratio that generated the carried overlap"
        );

        let first_frames = (first_chunk_len - fft_size) / hop + 1;
        let consumed = first_frames * hop;
        pv.set_stretch_ratio(0.96);
        let second = pv
            .process_streaming(&input[consumed..consumed + first_chunk_len])
            .unwrap();
        assert!(
            !second.is_empty(),
            "second streaming chunk should emit audio"
        );
        assert!(
            !pv.streaming_tail.is_empty(),
            "second streaming chunk should continue carrying overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 0.96).abs() < 1e-12,
            "once the previous overlap has fully emitted, small cross-unity modulation should re-arm to the current ratio"
        );

        let _ = pv.flush_streaming().unwrap();
        assert!(
            (pv.streaming_tail_ratio - 0.96).abs() < 1e-12,
            "flush should clear carried overlap history and re-arm the current ratio"
        );
    }

    #[test]
    fn test_streaming_tail_ratio_only_preserves_prior_ratio_while_overlap_remains_unresolved() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 5;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                (2.0 * PI * 220.0 * t).sin() * 0.55 + (2.0 * PI * 660.0 * t).sin() * 0.20
            })
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.18, sample_rate, 120.0);
        let first_chunk_len = fft_size * 2;
        let first = pv.process_streaming(&input[..first_chunk_len]).unwrap();
        assert!(!first.is_empty(), "first streaming chunk should emit audio");
        assert!(
            !pv.streaming_tail.is_empty(),
            "first streaming chunk should retain overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.18).abs() < 1e-12,
            "tail ratio should match the ratio that generated the carried overlap"
        );

        let first_frames = (first_chunk_len - fft_size) / hop + 1;
        let consumed = first_frames * hop;
        let second_chunk_len = fft_size + hop;
        pv.set_stretch_ratio(0.82);
        let second = pv
            .process_streaming(&input[consumed..consumed + second_chunk_len])
            .unwrap();
        assert!(
            !second.is_empty(),
            "short second streaming chunk should still emit audio"
        );
        assert!(
            !pv.streaming_tail.is_empty(),
            "short second streaming chunk should continue carrying overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.18).abs() < 1e-12,
            "the prior expansion ratio should persist only while unresolved overlap from that chunk remains in the carried tail"
        );
    }

    #[test]
    fn test_streaming_tail_ratio_holds_prior_seam_across_repeated_short_interval_modulation() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 8;
        let input: Vec<f32> = (0..num_samples)
            .map(|i| {
                let t = i as f32 / sample_rate as f32;
                (2.0 * PI * 220.0 * t).sin() * 0.55 + (2.0 * PI * 660.0 * t).sin() * 0.20
            })
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.04, sample_rate, 120.0);
        let first_chunk_len = fft_size * 2;
        let first = pv.process_streaming(&input[..first_chunk_len]).unwrap();
        assert!(!first.is_empty(), "first streaming chunk should emit audio");
        assert!(
            !pv.streaming_tail.is_empty(),
            "first streaming chunk should retain overlap tail"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.04).abs() < 1e-12,
            "tail ratio should match the ratio that generated the carried overlap"
        );

        let mut consumed = ((first_chunk_len - fft_size) / hop + 1) * hop;
        let short_chunk_len = fft_size + hop;

        pv.set_stretch_ratio(0.96);
        let second = pv
            .process_streaming(&input[consumed..consumed + short_chunk_len])
            .unwrap();
        assert!(
            !second.is_empty(),
            "second short streaming chunk should still emit audio"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.04).abs() < 1e-12,
            "the unresolved prior seam should keep its expansion ratio through the first cross-unity modulation step"
        );

        consumed += ((short_chunk_len - fft_size) / hop + 1) * hop;
        pv.set_stretch_ratio(1.02);
        let third = pv
            .process_streaming(&input[consumed..consumed + short_chunk_len])
            .unwrap();
        assert!(
            !third.is_empty(),
            "third short streaming chunk should still emit audio"
        );
        assert!(
            (pv.streaming_tail_ratio - 1.04).abs() < 1e-12,
            "repeated short-interval modulation should keep the unresolved prior seam ratio instead of re-arming on each toggle"
        );

        consumed += ((short_chunk_len - fft_size) / hop + 1) * hop;
        pv.set_stretch_ratio(0.98);
        let fourth = pv
            .process_streaming(&input[consumed..consumed + first_chunk_len])
            .unwrap();
        assert!(
            !fourth.is_empty(),
            "a longer follow-up chunk should still emit audio"
        );
        assert!(
            (pv.streaming_tail_ratio - 0.98).abs() < 1e-12,
            "once the older overlap has drained, the carried tail should re-arm to the current modulation ratio"
        );

        let _ = pv.flush_streaming().unwrap();
        assert!(
            (pv.streaming_tail_ratio - 0.98).abs() < 1e-12,
            "flush should preserve the re-armed current ratio after repeated short-interval modulation"
        );
    }

    #[test]
    fn test_flush_streaming_is_idempotent() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;
        let input: Vec<f32> = (0..fft_size * 3)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);
        let _ = pv.process_streaming(&input).unwrap();

        let _first = pv.flush_streaming().unwrap();
        let second = pv.flush_streaming().unwrap();
        assert!(
            second.is_empty(),
            "Second flush_streaming() call should be empty"
        );
    }

    // --- sub_bass_bin edge cases ---

    #[test]
    fn test_sub_bass_bin_clamped_to_num_bins() {
        // Very high cutoff: sub_bass_bin should be clamped to num_bins
        let pv = PhaseVocoder::new(256, 64, 1.0, 44100, 30000.0);
        let num_bins = 256 / 2 + 1;
        assert!(
            pv.sub_bass_bin() <= num_bins,
            "sub_bass_bin {} should be <= num_bins {}",
            pv.sub_bass_bin(),
            num_bins
        );
    }

    #[test]
    fn test_sub_bass_all_bins_rigid() {
        // With cutoff >= Nyquist, all bins should use rigid locking.
        // This should still produce valid output (no crash).
        let fft_size = 512;
        let hop = 128;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        // Cutoff at Nyquist: all bins are "sub-bass" → all rigid locking
        let mut pv = PhaseVocoder::new(fft_size, hop, 1.5, sample_rate, 22050.0);
        let output = pv.process(&input).unwrap();
        assert!(!output.is_empty());
        assert!(output.iter().all(|s| s.is_finite()));
    }

    // --- reconstruct_spectrum conjugate symmetry ---

    #[test]
    fn test_reconstruct_spectrum_produces_real_output() {
        // After reconstruct_spectrum + inverse FFT, output should be real-valued
        // (imaginary parts near zero). This verifies conjugate symmetry is correct.
        let fft_size = 256;
        let hop = 64;
        let sample_rate = 44100u32;
        let num_samples = fft_size * 4;

        let input: Vec<f32> = (0..num_samples)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);
        let output = pv.process(&input).unwrap();

        // If conjugate symmetry is wrong, we'd get complex residues causing
        // large imaginary parts. The output being finite and reasonable is evidence.
        assert!(output.iter().all(|s| s.is_finite()));
        let rms = (output.iter().map(|x| x * x).sum::<f32>() / output.len() as f32).sqrt();
        assert!(
            rms > 0.01,
            "Output should have significant energy, got RMS={}",
            rms
        );
    }

    // --- PV reuse (buffers grow but don't shrink) ---

    #[test]
    fn test_phase_vocoder_reuse_across_different_lengths() {
        let fft_size = 1024;
        let hop = 256;
        let sample_rate = 44100u32;

        let mut pv = PhaseVocoder::new(fft_size, hop, 1.0, sample_rate, 120.0);

        // Process a long signal
        let long_input: Vec<f32> = (0..fft_size * 8)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();
        let output1 = pv.process(&long_input).unwrap();
        assert!(!output1.is_empty());

        // Process a shorter signal — buffers should still work (they don't shrink)
        let short_input: Vec<f32> = (0..fft_size * 2)
            .map(|i| (2.0 * PI * 440.0 * i as f32 / sample_rate as f32).sin())
            .collect();
        let output2 = pv.process(&short_input).unwrap();
        assert!(!output2.is_empty());
        assert!(output2.len() < output1.len());
    }
}