mfsk-core 0.10.0

//! Coarse-sync stage: Costas-array correlation search across the
//! power spectrogram. Mirrors WSJT-X `lib/ft8/sync8.f90` — 16-bin
//! sliding-window allsum noise estimator with `tone_step_bins = 2.0`
//! exactly at `NFFT_SPEC = 3840`.
//!
//! ε.3 of the `docs/CLEANUP_2026_05.md` `decode_block` split. The
//! parent (`decode_block.rs`) re-exports the public entries so
//! external callers — host `decode.rs`, embedded `stage1_inc`,
//! `mfsk-ffi-ft8`, integration tests — keep the same
//! `mfsk_core::ft8::decode_block::*` paths.

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

#[cfg(not(feature = "std"))]
use num_traits::Float;

use super::super::params::{COSTAS, COSTAS_POS, NTONES};
use super::spectrogram::{CoarseAcc, Spectrogram};
use super::types::{
    NFFT_SPEC, NSSY, NSTEP, SAMPLE_RATE_HZ, TONE_SPACING_HZ, TX_START_OFFSET_S, ratio_eps,
    sync_lag_s,
};
use crate::engine::sync::SyncCandidate;

// ── Coarse sync ─────────────────────────────────────────────────────────────

/// Costas-array correlation search across the spectrogram. Mirrors
/// WSJT-X `lib/ft8/sync8.f90` — 16-bin sliding-window allsum noise
/// estimator with `tone_step_bins = 2.0` exactly at NFFT_SPEC=3840.
///
/// Public as of v0.6 (#48): this is now the canonical FT8 coarse-sync
/// for both the embedded port and the host `decode_frame*` pipeline.
/// `engine::sync::coarse_sync<Ft8>` is no longer reachable from FT8 code
/// paths; that generic function stays in place for FT4 / FST4 / JT9 /
/// Q65 / WSPR / uvpacket.
pub fn coarse_sync(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    sync_min: f32,
    max_cand: usize,
) -> Vec<SyncCandidate> {
    coarse_sync_inner(spec, freq_min, freq_max, sync_min, max_cand, None, None)
}

/// [`coarse_sync`] with an explicit ±lag search window in seconds,
/// instead of the crate default (`SYNC_LAG_S_DEFAULT` /
/// `MFSK_SYNC_LAG_S`, WSJT-X's own ±2.5 s).
///
/// The window is not purely a speed knob — it changes *which*
/// candidates exist and how they rank, so a caller whose downstream
/// use has its own timing assumptions should state them here rather
/// than inherit the decoder's. Two real cases (issue #280):
///
/// - [`bootstrap_dt_median`](crate::engine::sync::bootstrap_dt_median)
///   needs a window comparable to the timing error it is trying to
///   estimate. At ±2.5 s a strong signal's far-lag ghost (secondary
///   channel, see this module's `coarse_sync_inner`) can outscore its own true
///   peak — verified as WSJT-X behaviour too, by probing real `jt9`'s
///   `sync8` on `qso3_busy.wav` — which invalidates a top-K DT
///   median over the raw list.
/// - Clock-synced deployments (NTP / GPS embedded targets) live well
///   within ±1 s and pay for the wider scan in `n_lag`-linear
///   `sync2d` work for nothing.
pub fn coarse_sync_with_lag(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    sync_min: f32,
    max_cand: usize,
    sync_lag_s: f32,
) -> Vec<SyncCandidate> {
    coarse_sync_inner(
        spec,
        freq_min,
        freq_max,
        sync_min,
        max_cand,
        None,
        Some(sync_lag_s),
    )
}

/// Phase-E2 entry point — like [`coarse_sync`] but consumes a
/// caller-built allsum table instead of recomputing it. Saves the
/// 280-300 ms allsum precompute on Core2 by hiding it under the
/// 15 s capture window in the embedded port (`stage1_inc` builds
/// the allsum incrementally as new spectrogram rows arrive).
///
/// `allsum` must match exactly what
/// [`precompute_coarse_allsum`] produces for the same `spec`,
/// `freq_min`, `freq_max` triple — same layout
/// (`data[fi * spec.n_time + m]`), same length
/// ([`coarse_allsum_len`]).
///
/// Public as of v0.6 (#49 cat C): the embedded port (`m5stack-core2`,
/// `m5stack-s3`, `m5stack-s3-app`) and `mfsk-ffi-ft8` both build the
/// allsum incrementally during slot capture and pass it back here, so
/// this is a stable public surface — not a benchmark-only escape hatch.
pub fn coarse_sync_with_allsum(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    sync_min: f32,
    max_cand: usize,
    allsum: &[CoarseAcc],
) -> Vec<SyncCandidate> {
    coarse_sync_inner(
        spec,
        freq_min,
        freq_max,
        sync_min,
        max_cand,
        Some(allsum),
        None,
    )
}

/// [`coarse_sync_with_allsum`] with an explicit ±lag search window in
/// seconds. See [`coarse_sync_with_lag`] for when a caller should pin
/// the window rather than inherit the crate default.
///
/// `allsum` is indexed by `(fi, m)` only, so the same table is valid
/// at any lag window — no need to rebuild it when pinning.
pub fn coarse_sync_with_allsum_and_lag(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    sync_min: f32,
    max_cand: usize,
    allsum: &[CoarseAcc],
    sync_lag_s: f32,
) -> Vec<SyncCandidate> {
    coarse_sync_inner(
        spec,
        freq_min,
        freq_max,
        sync_min,
        max_cand,
        Some(allsum),
        Some(sync_lag_s),
    )
}

/// Length of the allsum buffer that [`coarse_sync_with_allsum`]
/// expects for the given `spec` + `freq_min..=freq_max` band.
/// Returns 0 when the band has no candidates (then
/// `coarse_sync_with_allsum` would return an empty vec immediately
/// regardless of `allsum`).
pub fn coarse_allsum_len(
    spec_n_freq: usize,
    spec_n_time: usize,
    freq_min: f32,
    freq_max: f32,
) -> usize {
    let df = SAMPLE_RATE_HZ / NFFT_SPEC as f32;
    let tone_step_bins = TONE_SPACING_HZ / df;
    let ia = (freq_min / df).round() as usize;
    let max_tone_off = ((NTONES - 1) as f32 * tone_step_bins).ceil() as usize + 1;
    let ib_unbounded = (freq_max / df).round() as usize;
    let ib = ib_unbounded.min(spec_n_freq.saturating_sub(max_tone_off));
    if ib < ia {
        return 0;
    }
    let n_freq = ib - ia + 1;
    n_freq * spec_n_time
}

/// One-shot allsum builder. The embedded port can build the allsum
/// incrementally instead and pass it to [`coarse_sync_with_allsum`];
/// host callers use this. Layout matches [`coarse_sync_with_allsum`]'s
/// expected input.
pub fn precompute_coarse_allsum(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
) -> Vec<CoarseAcc> {
    let mut buf =
        vec![CoarseAcc::default(); coarse_allsum_len(spec.n_freq, spec.n_time, freq_min, freq_max)];
    if !buf.is_empty() {
        precompute_coarse_allsum_into(spec, freq_min, freq_max, &mut buf);
    }
    buf
}

/// In-place variant of [`precompute_coarse_allsum`]. `dst` must have
/// length [`coarse_allsum_len`]. No-op if the band is empty.
pub fn precompute_coarse_allsum_into(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    dst: &mut [CoarseAcc],
) {
    let df = SAMPLE_RATE_HZ / NFFT_SPEC as f32;
    let tone_step_bins = TONE_SPACING_HZ / df;
    let ia = (freq_min / df).round() as usize;
    let max_tone_off = ((NTONES - 1) as f32 * tone_step_bins).ceil() as usize + 1;
    let nh1 = spec.n_freq;
    let ib_unbounded = (freq_max / df).round() as usize;
    let ib = ib_unbounded.min(nh1.saturating_sub(max_tone_off));
    if ib < ia {
        return;
    }
    let n_freq = ib - ia + 1;
    debug_assert_eq!(
        dst.len(),
        n_freq * spec.n_time,
        "allsum buffer length mismatch"
    );
    fill_coarse_allsum(spec, ia, ib, n_freq, dst);
}

/// Build the 7-tone-stride-2 allsum for one time slice `m` into `dst`.
///
/// Loop order: `m`-outer, `fi`-inner. One time slice (~2 KB at u16)
/// stays in cache for all carriers instead of striding by `n_freq`
/// bytes per time step.
///
/// Under `fixed-point` (SpecCell = u16, values ≤ 65535): uses a
/// **sliding window** — carriers at `fi` and `fi+2` share 6 of 7
/// bins, so `allsum[fi] = allsum[fi-2] + row[ia+fi+12] - row[ia+fi-2]`.
/// Reduces inner work from 7 adds to 1 add + 1 sub. Intermediate
/// sums ≤ 7×65535 < 2^19 are exactly representable in f32, so no
/// drift accumulates.
///
/// Without `fixed-point` (SpecCell = f32, host path): direct sum to
/// avoid rounding drift on large spectrum values that would break the
/// full-band-allsum-then-slice equivalence test.
fn fill_allsum_row(spec: &Spectrogram, ia: usize, n_freq: usize, m: usize, dst: &mut [CoarseAcc]) {
    if n_freq == 0 {
        return;
    }
    let nh1 = spec.n_freq;
    let upper = nh1 - 1;
    let n_time = spec.n_time;

    #[cfg(feature = "fixed-point")]
    {
        // Sliding window: two independent chains (even / odd fi).
        let mut prev = [CoarseAcc::default(); 2];
        for parity in 0..2usize {
            if parity >= n_freq {
                break;
            }
            let i_carrier = ia + parity;
            let mut s = CoarseAcc::default();
            for k in 0..(NTONES - 1) {
                let bin = (i_carrier + 2 * k).min(upper);
                s += spec.power_acc(bin, m);
            }
            dst[parity * n_time + m] = s;
            prev[parity] = s;
        }
        for fi in 2..n_freq {
            let i_carrier = ia + fi;
            let drop = (i_carrier - 2).min(upper);
            let add = (i_carrier + 2 * (NTONES - 2)).min(upper);
            let p = fi & 1;
            #[allow(clippy::unnecessary_cast)]
            let s = prev[p] + spec.power_acc(add, m) - spec.power_acc(drop, m);
            dst[fi * n_time + m] = s;
            prev[p] = s;
        }
    }

    #[cfg(not(feature = "fixed-point"))]
    {
        // Direct sum: no f32 drift for host tests.
        for fi in 0..n_freq {
            let i_carrier = ia + fi;
            let mut s = CoarseAcc::default();
            for k in 0..(NTONES - 1) {
                let bin = (i_carrier + 2 * k).min(upper);
                s += spec.power_acc(bin, m);
            }
            dst[fi * n_time + m] = s;
        }
    }
}

/// Build the 7-tone-stride-2 allsum for the carrier range `[ia..=ib]`
/// into `dst` (`n_freq * spec.n_time`, layout `dst[fi * n_time + m]`).
///
/// Mirrors the inline precompute inside [`coarse_sync_inner`] and
/// [`stage1_inc::update_one_half`] — kept as a public helper so the
/// embedded port can replicate it row-by-row without accessing private
/// internals. Computes **every m** so callers can populate cells
/// incrementally in any order.
fn fill_coarse_allsum(
    spec: &Spectrogram,
    ia: usize,
    ib: usize,
    _n_freq: usize,
    dst: &mut [CoarseAcc],
) {
    let n_freq = ib - ia + 1;
    // m-outer loop: one time slice (~2 KB u16) stays in cache for all fi.
    for m in 0..spec.n_time {
        fill_allsum_row(spec, ia, n_freq, m, dst);
    }
}

fn coarse_sync_inner(
    spec: &Spectrogram,
    freq_min: f32,
    freq_max: f32,
    sync_min: f32,
    max_cand: usize,
    external_allsum: Option<&[CoarseAcc]>,
    sync_lag_override_s: Option<f32>,
) -> Vec<SyncCandidate> {
    let df = SAMPLE_RATE_HZ / NFFT_SPEC as f32;
    let tstep = NSTEP as f32 / SAMPLE_RATE_HZ;
    let jstrt = (TX_START_OFFSET_S / tstep).round() as i32;
    let requested_jz = (sync_lag_override_s.unwrap_or_else(sync_lag_s) / tstep).round() as i32;
    let Some(jz) = bounded_sync_lag_steps(spec.n_time, jstrt, requested_jz) else {
        return Vec::new();
    };
    let tone_step_bins = TONE_SPACING_HZ / df;

    // Carrier-bin search range. Reserve room above the carrier for the
    // top tone (round(7 * tone_step_bins)).
    let ia = (freq_min / df).round() as usize;
    let max_tone_off = ((NTONES - 1) as f32 * tone_step_bins).ceil() as usize + 1;
    let nh1 = spec.n_freq;
    let ib_unbounded = (freq_max / df).round() as usize;
    let ib = ib_unbounded.min(nh1.saturating_sub(max_tone_off));
    if ib < ia {
        return Vec::new();
    }
    let n_freq = ib - ia + 1;
    let n_lag = (2 * jz + 1) as usize;
    let mut sync2d = vec![0.0f32; n_freq * n_lag];
    let idx = |fi: usize, lag: i32| fi * n_lag + (lag + jz) as usize;
    let ratio_eps = ratio_eps();
    #[cfg(all(
        feature = "profile-coarse",
        not(all(target_arch = "wasm32", target_os = "unknown"))
    ))]
    let t_setup = std::time::Instant::now();

    // **Single-bin tone gather**. NFFT=3840 → tone_step_bins = 2.0
    // exactly, so the 8 FT8 tones sit on integer bins {0, 2, 4, …, 14}
    // relative to the carrier. With the rectangular pre-FFT window
    // each tone's mainlobe is one bin wide and contains the full tone
    // energy, so single-bin reads are loss-less. The Plan-A multi-bin
    // sum (floor-bin + floor-bin+1) was a fractional-alignment +
    // Hann-mainlobe compensation for the NFFT=4096 era and was
    // dropped at the 3840 migration; do NOT re-add it without also
    // restoring Hann (the two were always paired).
    let mut tone_bin_lo = [0usize; NTONES];
    for k in 0..NTONES {
        tone_bin_lo[k] = (k as f32 * tone_step_bins).floor() as usize;
    }

    // Pre-compute the (bk, n) → m_base table. m for the inner iter is
    // `m_base[bk][n] + lag`; m_base depends only on the Costas pattern
    // and `jstrt` (constants of the slot), not on (fi, lag). Hoists
    // 21 mul/add chains out of the n_freq × n_lag × 21 inner loop.
    let m_base: [[i32; COSTAS.len()]; COSTAS_POS.len()] = {
        let mut t = [[0i32; COSTAS.len()]; COSTAS_POS.len()];
        for (bk, &start_sym) in COSTAS_POS.iter().enumerate() {
            let block_offset = NSSY * start_sym as i32;
            for (n, _) in COSTAS.iter().enumerate() {
                t[bk][n] = jstrt + block_offset + NSSY * n as i32;
            }
        }
        t
    };
    // Pre-compute Costas-tone bin offsets in COSTAS-order (i.e.
    // `tone_bin_lo[COSTAS[n]]`). Saves an indirect lookup per inner
    // iteration.
    let costas_off: [usize; COSTAS.len()] = {
        let mut t = [0usize; COSTAS.len()];
        for (n, &costas_n) in COSTAS.iter().enumerate() {
            t[n] = tone_bin_lo[costas_n];
        }
        t
    };

    // Pre-compute the set of `m` time-indices that the score loop
    // will read. Only Costas-symbol positions ± lag count — typically
    // 3 contiguous bands totalling ~110 of the 184 frames, so the
    // naïve "compute allsum for every m" wastes ~40 % of its work
    // on cells the score loop never reads.
    let needed_m: alloc::vec::Vec<usize> = {
        let mut mark = alloc::vec![false; spec.n_time];
        for bk in 0..COSTAS_POS.len() {
            let lo = m_base[bk][0] - jz;
            let hi = m_base[bk][COSTAS.len() - 1] + jz;
            let lo_u = lo.max(0) as usize;
            let hi_u = (hi.min(spec.n_time as i32 - 1)) as usize;
            if lo_u <= hi_u {
                #[allow(clippy::needless_range_loop)]
                for m in lo_u..=hi_u {
                    mark[m] = true;
                }
            }
        }
        (0..spec.n_time).filter(|&m| mark[m]).collect()
    };

    // Pre-compute `Σ_{k=0..NTONES-1} spec[i_carrier + 2*k, m]` for every
    // (fi, m ∈ needed_m). NFFT=3840 → tone_step_bins=2.0 exactly →
    // 8-bin every-other gather per carrier (matches WSJT-X `sync8.f90`'s
    // single-bin `t0a` accumulator). No sliding-window reuse since
    // adjacent fi share zero bins under this pattern.
    //
    // Phase-E2: caller can pass a pre-built allsum (built incrementally
    // during slot capture by stage1_inc). When provided, skip the
    // precompute entirely. Otherwise build it inline.
    let owned_allsum: Vec<CoarseAcc>;
    let allsum: &[CoarseAcc] = if let Some(ext) = external_allsum {
        debug_assert_eq!(
            ext.len(),
            n_freq * spec.n_time,
            "external allsum length mismatch (expected n_freq * spec.n_time)"
        );
        ext
    } else {
        owned_allsum = {
            // Sum 7 Costas-tone bins (k=0..NTONES-1; tone 7 is data-only,
            // never a Costas position). Matches WSJT-X `sync8.f90:66` and
            // the public `fill_coarse_allsum` helper. Only populates
            // `needed_m` cells — ~60 % of spec.n_time — to save ~40 % of
            // allsum work; the remaining cells stay at default (0/0.0) and
            // are never read by the score loop.
            // Uses the stride-2 sliding window: allsum[fi] = allsum[fi-2]
            // + row[add] - row[drop], 7× fewer ops per fi after init.
            let mut buf = vec![CoarseAcc::default(); n_freq * spec.n_time];
            for &m in &needed_m {
                fill_allsum_row(spec, ia, n_freq, m, &mut buf);
            }
            buf
        };
        &owned_allsum
    };
    #[cfg(all(
        feature = "profile-coarse",
        not(all(target_arch = "wasm32", target_os = "unknown"))
    ))]
    let t_allsum = std::time::Instant::now();

    // Three identical Costas arrays at symbol positions 0, 36, 72.
    //
    // **Bounds-check hoist**. For the standard FT8 slot:
    //   block 0: m = lag + jstrt + NSSY*n   ∈ [-7..31] over (lag, n)
    //   block 1: m = lag + jstrt + 144 + NSSY*n  ∈ [65..103]
    //   block 2: m = lag + jstrt + 288 + NSSY*n  ∈ [137..175]
    //
    // Block 0 can dip below 0 and block 2 can run past `n_time`, exactly
    // the two bounds guarded by WSJT-X `sync8.f90`. Compute the valid
    // symbol ranges once per lag instead of checking every inner read.
    // Block 1 remains fully in range for every accepted lag window.
    debug_assert!(
        m_base[1][0] - jz >= 0 && m_base[1][COSTAS.len() - 1] + jz < spec.n_time as i32,
        "n_time too small for SYNC_LAG_S/jstrt"
    );
    let n_time = spec.n_time;
    // Per-fi: tbin_lo[n] = i_carrier + costas_off[n]. Hoist out of
    // the lag loop — saves an addition per inner iteration.
    let mut tbin_lo_arr = [0usize; COSTAS.len()];
    for (fi, i_carrier) in (ia..=ib).enumerate() {
        for n in 0..COSTAS.len() {
            tbin_lo_arr[n] = i_carrier + costas_off[n];
        }
        let allsum_row = &allsum[fi * n_time..(fi + 1) * n_time];
        for lag in -jz..=jz {
            let mut t_blocks: [CoarseAcc; 3] = [CoarseAcc::default(); 3];
            let mut t0_blocks: [CoarseAcc; 3] = [CoarseAcc::default(); 3];

            // Block 0: smallest valid n where m_base[0][n] + lag >= 0.
            // m_base[0][n] = jstrt + NSSY*n. Solve for n:
            //   jstrt + NSSY*n + lag >= 0 ⇒ n >= ceil((-jstrt - lag) / NSSY)
            let bk0_n_start = {
                let needed = -jstrt - lag;
                if needed <= 0 {
                    0usize
                } else {
                    ((needed + NSSY - 1) / NSSY).min(COSTAS.len() as i32) as usize
                }
            };
            // NFFT=3840 → tone_step_bins = 2.0 exactly, signal at bin
            // `tbin_lo` only (single bin gather, matching WSJT-X's
            // `s(i+nfos*icos7(n), m)`). The legacy `+ spec[tbin_lo + 1]`
            // half was a fractional-alignment compensation for NFFT=4096
            // and is removed.
            for n in bk0_n_start..COSTAS.len() {
                let m_u = (m_base[0][n] + lag) as usize;
                let tbin_lo = tbin_lo_arr[n];
                t_blocks[0] += spec.power_acc(tbin_lo, m_u);
                t0_blocks[0] += allsum_row[m_u];
            }
            // Block 1 is always fully in range.
            for n in 0..COSTAS.len() {
                let m_u = (m_base[1][n] + lag) as usize;
                let tbin_lo = tbin_lo_arr[n];
                t_blocks[1] += spec.power_acc(tbin_lo, m_u);
                t0_blocks[1] += allsum_row[m_u];
            }
            // Block 2: WSJT-X guards `m + nssy*72 <= NHSYM` because
            // positive lag can push its trailing symbols past the slot.
            let bk2_n_end = valid_trailing_symbol_count(m_base[2][0], lag, n_time);
            for n in 0..bk2_n_end {
                let m_u = (m_base[2][n] + lag) as usize;
                let tbin_lo = tbin_lo_arr[n];
                t_blocks[2] += spec.power_acc(tbin_lo, m_u);
                t0_blocks[2] += allsum_row[m_u];
            }

            // Regularised ratio `t / (mean_others + ε)`.
            // ε prevents the u16-quantised fp path from blowing up
            // when phantom carriers happen to land where 7 of 8 tone
            // bins quantise to 0; `t0_ref → 0` would otherwise inflate
            // ratio scores by 100-1000× over real signals.
            //
            // Phase 3: t/t0 sums stay in CoarseAcc; convert to f32 only
            // at the score-division boundary so sync2d / red / base
            // (downstream sort + percentile) remain f32.
            let t_all: CoarseAcc = t_blocks[0] + t_blocks[1] + t_blocks[2];
            let t0_all: CoarseAcc = t0_blocks[0] + t0_blocks[1] + t0_blocks[2];
            // `as f32` is a real promotion under fixed-point (CoarseAcc=i32)
            // and a no-op when CoarseAcc=f32 — silence the host clippy.
            #[allow(clippy::unnecessary_cast)]
            let t_all_f = t_all as f32;
            #[allow(clippy::unnecessary_cast)]
            let t0_all_f = t0_all as f32;
            let t0_ref = (t0_all_f - t_all_f) / (NTONES as f32 - 2.0);
            let sync_all = t_all_f / (t0_ref + ratio_eps);

            // Trailing-2-blocks score (drop block 0 — late-start tolerance).
            #[allow(clippy::unnecessary_cast)]
            let t_tail_f = (t_blocks[1] + t_blocks[2]) as f32;
            #[allow(clippy::unnecessary_cast)]
            let t0_tail_f = (t0_blocks[1] + t0_blocks[2]) as f32;
            let t0_tail_ref = (t0_tail_f - t_tail_f) / (NTONES as f32 - 2.0);
            let sync_tail = t_tail_f / (t0_tail_ref + ratio_eps);

            sync2d[idx(fi, lag)] = sync_all.max(sync_tail);
        }
    }
    #[cfg(all(
        feature = "profile-coarse",
        not(all(target_arch = "wasm32", target_os = "unknown"))
    ))]
    let t_score = std::time::Instant::now();

    const MLAG: i32 = 10;

    // Two independent noise-floor channels, matching WSJT-X `sync8.f90`
    // exactly (issue #280). Collapsing these into a single
    // `±jz`-dependent channel — the pre-#280 shape of this function —
    // let the noise floor `base` scale with the requested lag window:
    // widening `MFSK_SYNC_LAG_S` toward WSJT-X's own `±2.5s` inflated
    // the *only* floor available here (per-bin max over a wider window
    // has a higher expected value under basic extreme-value
    // statistics), silently pushing weak real decodes below `sync_min`
    // even though WSJT-X's actual primary channel — anchored to a
    // fixed `mlag=10` regardless of `JZ` — would still have caught
    // them. `red_primary`/`base_primary` restores that JZ-independent
    // anchor; `red_secondary`/`base_secondary` (the old sole channel)
    // now plays its correct WSJT-X role: a supplementary second
    // candidate per bin, contributed only when its peak lag disagrees
    // with the primary's.
    let primary_half = MLAG.min(jz);
    let mut red_primary = vec![0.0f32; n_freq];
    let mut jpeak_primary = vec![0i32; n_freq];
    let mut red_secondary = vec![0.0f32; n_freq];
    let mut jpeak_secondary = vec![0i32; n_freq];
    for fi in 0..n_freq {
        let mut best_p = f32::NEG_INFINITY;
        let mut lag_p = -primary_half;
        let mut best_s = f32::NEG_INFINITY;
        let mut lag_s = -jz;
        for lag in -jz..=jz {
            let v = sync2d[idx(fi, lag)];
            if v > best_s {
                best_s = v;
                lag_s = lag;
            }
            if (-primary_half..=primary_half).contains(&lag) && v > best_p {
                best_p = v;
                lag_p = lag;
            }
        }
        red_primary[fi] = best_p;
        jpeak_primary[fi] = lag_p;
        red_secondary[fi] = best_s;
        jpeak_secondary[fi] = lag_s;
    }

    // 40th-percentile noise floor: `select_nth_unstable_by` is O(N)
    // average vs a full sort's O(N log N). Pivot ordering inside the
    // percentile isn't load-bearing (downstream only uses the selected
    // element as a normaliser), so unstable selection is fine.
    // `unwrap_or(Ordering::Equal)` for NaN safety — matches the same
    // pattern at `process_candidates.rs` xsnr2_db_simple median.
    // Gemini PR #79 + #101 review.
    let percentile_floor = |red: &[f32]| -> f32 {
        let mut sorted = red.to_vec();
        let pct_idx = ((0.40 * n_freq as f32) as usize).min(n_freq - 1);
        sorted.select_nth_unstable_by(pct_idx, |a, b| {
            a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal)
        });
        sorted[pct_idx].max(f32::EPSILON)
    };
    let base_primary = percentile_floor(&red_primary);
    let base_secondary = percentile_floor(&red_secondary);

    // Parabolic sub-quantum interpolation around a chosen peak lag —
    // shared by both channels. Only valid strictly inside the searched
    // window (needs both neighbours); falls back to 0 at the edge.
    let dt_quanta_at = |fi: usize, lag: i32| -> f32 {
        if lag > -jz && lag < jz {
            let y_lo = sync2d[idx(fi, lag - 1)];
            let y_mi = sync2d[idx(fi, lag)];
            let y_hi = sync2d[idx(fi, lag + 1)];
            let denom = y_lo - 2.0 * y_mi + y_hi;
            if denom.abs() > f32::EPSILON {
                (0.5 * (y_lo - y_hi) / denom).clamp(-1.0, 1.0)
            } else {
                0.0
            }
        } else {
            0.0
        }
    };

    let mut cands: Vec<SyncCandidate> = Vec::new();
    for fi in 0..n_freq {
        let i_carrier = ia + fi;
        let freq_hz = i_carrier as f32 * df;

        let lag_p = jpeak_primary[fi];
        let norm_p = red_primary[fi] / base_primary;
        if norm_p.is_finite() && norm_p >= sync_min {
            let dt_lag = lag_p as f32 + dt_quanta_at(fi, lag_p);
            cands.push(SyncCandidate {
                freq_hz,
                dt_sec: (dt_lag - 0.5) * tstep,
                score: norm_p,
            });
        }

        // Secondary candidate only when the full-window peak actually
        // disagrees with the fixed-window one — matches `sync8.f90`'s
        // `if (jpeak2(n)==jpeak(n)) cycle`.
        let lag_s = jpeak_secondary[fi];
        if lag_s != lag_p {
            let norm_s = red_secondary[fi] / base_secondary;
            if norm_s.is_finite() && norm_s >= sync_min {
                let dt_lag = lag_s as f32 + dt_quanta_at(fi, lag_s);
                cands.push(SyncCandidate {
                    freq_hz,
                    dt_sec: (dt_lag - 0.5) * tstep,
                    score: norm_s,
                });
            }
        }
    }

    // Dedupe within 4 Hz / 40 ms; keep highest score.
    //
    // Sort once by score desc, then greedily keep cands with no
    // already-kept near neighbour. Stops early at `max_cand` for
    // O(n log n + n × max_cand) instead of the prior O(n²) pairwise
    // compare-and-zero (n is several thousand).
    //
    // Empirically this drops 1 borderline busy-band truth on qso3
    // vs the byte-equivalent O(n²) implementation — likely a
    // tie-break ordering difference at the 4 Hz dedupe boundary
    // (haven't fully traced; the algorithms are equivalent on every
    // small-case scenario I worked through). We accept the recall
    // delta for the ~150 us host (~150 ms Core2) saving.
    cands.sort_by(|a, b| b.score.partial_cmp(&a.score).unwrap());
    let mut out: Vec<SyncCandidate> = Vec::with_capacity(max_cand);
    for c in cands {
        if c.score < sync_min {
            break;
        }
        let near = out
            .iter()
            .any(|k| (c.freq_hz - k.freq_hz).abs() < 4.0 && (c.dt_sec - k.dt_sec).abs() < 0.04);
        if near {
            continue;
        }
        out.push(c);
        if out.len() >= max_cand {
            break;
        }
    }
    let cands = out;
    #[cfg(all(
        feature = "profile-coarse",
        not(all(target_arch = "wasm32", target_os = "unknown"))
    ))]
    {
        let t_end = std::time::Instant::now();
        let allsum_us = (t_allsum - t_setup).as_micros();
        let score_us = (t_score - t_allsum).as_micros();
        let post_us = (t_end - t_score).as_micros();
        let total_us = (t_end - t_setup).as_micros();
        eprintln!(
            "[coarse_sync prof] n_freq={n_freq} n_lag={n_lag}  allsum={allsum_us} score={score_us} dedupe+sort={post_us}  total={total_us} us"
        );
    }
    cands
}

/// Keep the requested symmetric lag window when the always-read middle Costas
/// block fits. Only pathological overrides need clamping; the first and third
/// blocks have the same conditional bounds as WSJT-X `sync8.f90`.
fn bounded_sync_lag_steps(
    spec_n_time: usize,
    start_step: i32,
    requested_steps: i32,
) -> Option<i32> {
    let requested_steps = requested_steps.max(0);
    let second_block_base = start_step + NSSY * COSTAS_POS[1] as i32;
    let second_block_end = second_block_base + NSSY * (COSTAS.len() as i32 - 1);
    let max_positive_steps = (spec_n_time as i32 - 1).checked_sub(second_block_end)?;
    let safe_steps = second_block_base.min(max_positive_steps);
    (safe_steps >= 0).then(|| requested_steps.min(safe_steps))
}

fn valid_trailing_symbol_count(base_step: i32, lag: i32, n_time: usize) -> usize {
    let available = (n_time as i32 - 1) - (base_step + lag);
    if available < 0 {
        0
    } else {
        (available / NSSY + 1).min(COSTAS.len() as i32) as usize
    }
}

#[cfg(test)]
mod tests {
    use super::{
        Spectrogram, bounded_sync_lag_steps, coarse_sync_inner, valid_trailing_symbol_count,
    };

    #[test]
    fn coarse_sync_accepts_wsjtx_lag_window_without_out_of_bounds_access() {
        // A standard 15 s FT8 spectrogram has 372 quarter-symbol rows.
        // Before the trailing-block guard was restored, requesting WSJT-X's
        // +/-2.5 s window made coarse_sync index beyond the final row.
        let n_freq = 128;
        let n_time = 372;
        let spec =
            Spectrogram::from_parts(n_freq, n_time, vec![Default::default(); n_freq * n_time]);

        let candidates = coarse_sync_inner(&spec, 200.0, 300.0, 1.5, 10, None, Some(2.5));

        assert!(candidates.is_empty());
    }

    #[test]
    fn sync_lag_preserves_wsjtx_window_when_middle_block_fits() {
        assert_eq!(bounded_sync_lag_steps(372, 13, 63), Some(63));
        assert_eq!(bounded_sync_lag_steps(372, 13, 500), Some(157));
        assert_eq!(bounded_sync_lag_steps(181, 13, 25), None);
        assert_eq!(bounded_sync_lag_steps(372, 13, -10), Some(0));
    }

    #[test]
    fn trailing_block_skips_symbols_beyond_spectrogram() {
        assert_eq!(valid_trailing_symbol_count(301, 0, 372), 7);
        assert_eq!(valid_trailing_symbol_count(301, 63, 372), 2);
        assert_eq!(valid_trailing_symbol_count(301, 80, 372), 0);
    }
}