orion-sdr 0.0.60

// Copyright (c) 2026 G & R Associates LLC
// SPDX-License-Identifier: MIT OR Apache-2.0

// src/sync/dvb_t_gi_sync.rs
//
// Guard-interval (cyclic-prefix) acquisition for a preamble-less OFDM frame —
// the DVB-T way. A conformant DVB-T signal carries NO Schmidl & Cox preamble;
// the receiver finds the symbol boundary and fractional carrier-frequency
// offset directly from the cyclic prefix (the van de Beek / Sandell / Börjesson
// ML estimator, "ML estimation of time and frequency offset in OFDM systems",
// IEEE Trans. Signal Processing, 1997).
//
// Each OFDM symbol is `cp_len + n_fft` samples, where the first `cp_len`
// samples (the guard interval) are a copy of the symbol's last `cp_len`
// samples. So for a correct symbol-start offset `d`, the two windows
// `r[d .. d+cp_len]` and `r[d+n_fft .. d+n_fft+cp_len]` are identical up to
// noise and a phase ramp `exp(j·2π·ε·n_fft)` from a residual CFO `ε` (in
// subcarrier-spacing units). Define the correlation and energy terms
//
//     γ(d) = Σ_{k=0}^{cp_len-1} r[d+k] · conj(r[d+n_fft+k])
//     Φ(d) = ½ Σ_{k=0}^{cp_len-1} ( |r[d+k]|² + |r[d+n_fft+k]|² ).
//
// The ML timing estimate maximizes the log-likelihood metric
//
//     Λ(d) = |γ(d)| − ρ · Φ(d),
//
// where ρ = SNR/(SNR+1) = |correlation coefficient| weights the energy term by
// the signal reliability (ρ→1 at high SNR fully subtracts Φ; ρ→0 at low SNR
// leaves a pure correlation peak-pick). Ranking by |γ| alone — dropping ρ·Φ —
// is only the correct ML rule in the ρ→0 limit and spuriously rewards
// high-energy offsets; the full metric is used here. The fractional CFO comes
// from the angle at the winning offset:
//
//     ε̂ = −∠γ / (2π)        (cycles per n_fft samples)
//     cfo_hz = ε̂ · fs / n_fft = −∠γ · fs / (2π · n_fft).
//
// The estimate is unambiguous only within ±½ a subcarrier spacing (±fs/2n_fft);
// integer-CFO and frame/super-frame lock come from the scattered pilots and TPS
// downstream. This provides the timing + fractional-CFO acquisition an external
// IQ capture needs before any symbol can be demapped.
//
// DIVERGENCE FROM THE PUBLISHED ESTIMATOR. Van de Beek's Λ is derived over a
// SINGLE cyclic prefix. Here γ and Φ may optionally be accumulated coherently
// over up to `max_symbols` consecutive symbols at the same candidate offset
// (like gr-dtv's moving-average correlator), which sharpens the lock for a
// known-length batch with a stable channel. This is coherent, so a large
// residual CFO rotates successive symbols' γ and partially cancels the sum —
// keep `max_symbols` small (a few symbols) when a meaningful CFO may be present,
// or set it to 1 for the strict single-symbol estimator. Fractional-CFO
// estimation itself is unaffected (it reads the accumulated angle, which the
// per-symbol phase ramp shares).

use num_complex::Complex32 as C32;

/// Tuning for [`dvb_t_gi_sync_with`]: the ML energy-term weight and the coherent
/// accumulation bound.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GiSyncConfig {
    /// The van de Beek energy-term weight `ρ = SNR/(SNR+1) ∈ [0, 1]` in the
    /// timing metric `|γ| − ρ·Φ`. The default (`0.95`) suits the moderate-to-high
    /// SNR of a decodable DVB-T signal; lower it toward 0 for very low SNR (where
    /// the pure correlation magnitude is the better statistic).
    pub rho: f32,
    /// Maximum number of consecutive symbols whose CP correlation is accumulated
    /// coherently at each candidate offset. `1` is the strict single-symbol van
    /// de Beek estimator; a small value (the default, `4`) sharpens the lock for
    /// a batch while staying robust to residual CFO. Larger values risk coherent
    /// cancellation under CFO (see the module header).
    pub max_symbols: usize,
    /// Unwrapping guard: the fraction of the peak's correlation score
    /// (`|γ|/Φ`) that the period boundary at or before the peak must itself
    /// reach before it is reported instead. `0.0` disables unwrapping entirely
    /// and restores a plain argmax.
    ///
    /// Λ is periodic in the symbol period, so a peak at phase `−δ` can only
    /// surface at `period − δ`, at the very top of the search range. Reporting
    /// that offset is the right phase but the *next* symbol — nearly a whole
    /// symbol late, which overruns a caller that sized its buffer for the frame.
    /// TX symbol windowing puts the peak exactly there: the taper attenuates
    /// each symbol's leading cyclic-prefix samples but not their unwindowed
    /// copies in the symbol's interior, biasing the estimate early by roughly a
    /// third of `roll_off`. A negative offset is not representable, so the only
    /// admissible answer is the origin.
    ///
    /// The asymmetry is what justifies preferring it: a consumer absorbs a
    /// slightly EARLY estimate for free — the guard sits ahead of the FFT window
    /// and the residual is a phase ramp the channel estimate divides out — but
    /// cannot absorb a symbol-late one.
    ///
    /// **Two conditions must both hold**, and neither suffices alone:
    ///
    /// 1. The peak sits within `cp_len/2` below a period boundary, so its phase
    ///    really is a small negative number. A genuine lead-in of 200 samples
    ///    peaks at 200, nowhere near a boundary.
    /// 2. That boundary's own correlation score reaches `origin_score_ratio` of
    ///    the peak's, so a symbol plausibly does start there.
    ///
    /// Condition 2 compares **scores** (`|γ|/Φ`), not the ML metric Λ, and that
    /// choice is load-bearing. Λ rewards low energy — an all-silence window has
    /// `Φ ≈ 0` and so `Λ ≈ 0`, which sits *near the peak* — so comparing Λ would
    /// accept the origin of a buffer that begins with silence. The normalized
    /// score has no such bias: silence-against-signal correlates at ~0.03.
    ///
    /// Measured across guard intervals and shaping, the score ratio separates
    /// the two cases with a wide gap — 0.61–1.00 when the frame really does
    /// start at the boundary, 0.03–0.41 when it does not — so the default `0.5`
    /// sits between with ~20% margin either side. (For reference, Λ's own
    /// deficit does *not* separate them: at G1/32 a true wrap sits 5.5% below
    /// the peak and a genuine 200-sample lead-in 6.3%.)
    ///
    /// Note this considers exactly **two** candidates, the argmax and that one
    /// boundary. Scanning for the earliest of all near-equal offsets would be
    /// wrong: Λ's peak sits on a plateau roughly `cp_len` wide, so any width able
    /// to jump the wrap also walks tens of samples down that plateau and
    /// re-introduces the timing error it was meant to remove.
    pub origin_score_ratio: f32,
}

impl Default for GiSyncConfig {
    fn default() -> Self {
        Self {
            rho: 0.95,
            max_symbols: 4,
            origin_score_ratio: 0.5,
        }
    }
}

/// One guard-interval acquisition result: a candidate symbol-start offset with
/// its correlation strength and fractional-CFO estimate. Mirrors
/// [`crate::sync::OfdmSyncResult`]'s role for the preamble path.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GiSyncResult {
    /// Sample offset of the first full OFDM symbol's start (its cyclic prefix).
    pub start_sample: usize,
    /// Fractional CFO estimate (Hz), unambiguous within ±`fs / (2·n_fft)`.
    pub cfo_hz: f32,
    /// Normalized correlation score in `[0, 1]`: `|γ(d)| / Φ(d)` at the winning
    /// offset, higher is a stronger cyclic-prefix match (1 for identical
    /// noiseless windows). Reported for quality/thresholding; the *selection*
    /// uses the ML metric `|γ| − ρ·Φ`, not this ratio.
    pub score: f32,
}

/// Searches `iq` for the best guard-interval-aligned symbol start in the offset
/// range `0..search_len`, using the cyclic prefix of a `(n_fft, cp_len)` OFDM
/// symbol and the [default](GiSyncConfig::default) tuning. Returns `None` if
/// `iq` is too short to hold a full symbol plus the search span.
///
/// `search_len` should span at least one full symbol period (`n_fft + cp_len`)
/// so the true peak is included; callers typically pass one symbol period to
/// lock onto the symbol grid. See [`dvb_t_gi_sync_with`] to tune the ML energy
/// weight and the coherent-accumulation bound.
pub fn dvb_t_gi_sync(
    iq: &[C32],
    n_fft: usize,
    cp_len: usize,
    fs: f32,
    search_len: usize,
) -> Option<GiSyncResult> {
    dvb_t_gi_sync_with(iq, n_fft, cp_len, fs, search_len, &GiSyncConfig::default())
}

/// Like [`dvb_t_gi_sync`], with an explicit [`GiSyncConfig`] (`ρ` weight, the
/// coherent-accumulation bound `max_symbols`, and the `origin_score_ratio` unwrap guard).
/// Ranks offsets by the van de Beek ML timing metric `|γ(d)| − ρ·Φ(d)` and
/// reports the preceding period boundary instead when the peak wrapped past it —
/// see [`GiSyncConfig::origin_score_ratio`] for why that is not simply the argmax.
pub fn dvb_t_gi_sync_with(
    iq: &[C32],
    n_fft: usize,
    cp_len: usize,
    fs: f32,
    search_len: usize,
    cfg: &GiSyncConfig,
) -> Option<GiSyncResult> {
    if cp_len == 0 || n_fft == 0 {
        return None;
    }
    // Need, for the largest tested offset d = search_len-1, the window
    // r[d + n_fft .. d + n_fft + cp_len].
    let need = search_len.saturating_sub(1) + n_fft + cp_len;
    if iq.len() < need || search_len == 0 {
        return None;
    }

    let period = n_fft + cp_len;
    let max_syms = cfg.max_symbols.max(1);
    // Score every offset, keeping γ and Φ so the winner needs no second pass over
    // the (search_len × cp_len × max_syms) correlation — which dominates the cost.
    let mut scored: Vec<(f32, C32, f32)> = Vec::with_capacity(search_len);
    let mut max_metric = f32::NEG_INFINITY;
    let mut min_metric = f32::INFINITY;
    for d in 0..search_len {
        // Accumulate γ and Φ over up to `max_syms` consecutive symbols whose
        // CP+tail windows fit within the buffer (coherent for γ).
        let mut gamma = C32::default();
        let mut phi = 0.0f32;
        let mut base = d;
        let mut used = 0usize;
        while used < max_syms && base + n_fft + cp_len <= iq.len() {
            for k in 0..cp_len {
                let a = iq[base + k];
                let b = iq[base + n_fft + k];
                gamma += a * b.conj();
                phi += a.norm_sqr() + b.norm_sqr();
            }
            base += period;
            used += 1;
        }
        phi *= 0.5; // Φ = ½ Σ(|a|² + |b|²)

        // van de Beek ML timing metric: |γ| − ρ·Φ.
        let metric = gamma.norm() - cfg.rho * phi;
        max_metric = max_metric.max(metric);
        min_metric = min_metric.min(metric);
        scored.push((metric, gamma, phi));
    }

    let argmax = scored
        .iter()
        .enumerate()
        .max_by(|(_, a), (_, b)| a.0.total_cmp(&b.0))
        .map(|(d, _)| d)
        .unwrap_or(0);

    // Normalized single-symbol correlation |γ|/Φ ∈ [0, 1] at one offset.
    //
    // Deliberately NOT accumulated over `max_symbols`, unlike the selection
    // metric. Accumulation is what makes the unwrap test hard to get right: with
    // four symbols summed, an offset whose *first* symbol is silence still scores
    // well on the strength of the three behind it, so an accumulated score cannot
    // tell "a symbol starts here" from "a symbol starts one period from here".
    // A single symbol answers exactly the question being asked.
    let single_score = |d: usize| {
        if d + n_fft + cp_len > iq.len() {
            return 0.0;
        }
        let mut gamma = C32::default();
        let mut phi = 0.0f32;
        for k in 0..cp_len {
            let a = iq[d + k];
            let b = iq[d + n_fft + k];
            gamma += a * b.conj();
            phi += a.norm_sqr() + b.norm_sqr();
        }
        phi *= 0.5;
        if phi > 0.0 {
            (gamma.norm() / phi).min(1.0)
        } else {
            0.0
        }
    };

    // Unwrap a peak that landed just below a period boundary: that is phase `−δ`
    // for small δ — the next symbol, not this one. Both conditions are required;
    // see `GiSyncConfig::origin_score_ratio`.
    let phase = argmax % period;
    let origin = argmax - phase; // the period boundary at or before the peak
    let best_d = if cfg.origin_score_ratio > 0.0
        && phase != 0
        && period - phase <= cp_len.div_ceil(2)
        && single_score(origin) >= cfg.origin_score_ratio.clamp(0.0, 1.0) * single_score(argmax)
    {
        origin
    } else {
        argmax
    };
    let (_, best_gamma, best_phi) = scored[best_d];

    // Reported score: |γ| / Φ at the winning offset, in [0, 1] (1 for identical
    // noiseless windows). Selection used the ML metric above, not this ratio.
    let score = if best_phi > 0.0 {
        (best_gamma.norm() / best_phi).min(1.0)
    } else {
        0.0
    };
    // Fractional CFO: negate the accumulated angle so a positive applied CFO
    // reads as a positive estimate (ε̂ = −∠γ / 2π).
    let cfo_hz = -best_gamma.im.atan2(best_gamma.re) * fs / (core::f32::consts::TAU * n_fft as f32);
    Some(GiSyncResult {
        start_sample: best_d,
        cfo_hz,
        score,
    })
}

/// Refines a coarse offset to the best cyclic-prefix ML metric in a small window
/// `±radius` around `coarse`, using the [default](GiSyncConfig::default) tuning.
/// A convenience for a receiver that already has an approximate symbol boundary
/// (e.g. from a prior frame) and wants a cheap local re-lock rather than a full-
/// period search.
pub fn dvb_t_gi_refine(
    iq: &[C32],
    n_fft: usize,
    cp_len: usize,
    fs: f32,
    coarse: usize,
    radius: usize,
) -> Option<GiSyncResult> {
    dvb_t_gi_refine_with(
        iq,
        n_fft,
        cp_len,
        fs,
        coarse,
        radius,
        &GiSyncConfig::default(),
    )
}

/// Like [`dvb_t_gi_refine`], with an explicit [`GiSyncConfig`].
///
/// `cfg.origin_score_ratio` is deliberately ignored here (a plain argmax is used). The
/// tie-break exists to resolve a peak that wrapped to the far end of a
/// full-period search; a refine already knows the phase to within `±radius`, so
/// there is nothing to unwrap — and applying it would just bias every refine
/// toward the low end of its own window.
pub fn dvb_t_gi_refine_with(
    iq: &[C32],
    n_fft: usize,
    cp_len: usize,
    fs: f32,
    coarse: usize,
    radius: usize,
    cfg: &GiSyncConfig,
) -> Option<GiSyncResult> {
    let start = coarse.saturating_sub(radius);
    let span = 2 * radius + 1;
    let sub = iq.get(start..)?;
    let local = GiSyncConfig {
        origin_score_ratio: 0.0,
        ..*cfg
    };
    let mut r = dvb_t_gi_sync_with(sub, n_fft, cp_len, fs, span.min(sub.len()), &local)?;
    r.start_sample += start;
    Some(r)
}

// ── Integer-CFO estimation (continual-pilot spectral correlation) ────────────
//
// The guard-interval estimator above resolves the CFO only within ±½ a
// subcarrier spacing; a real front end can be off by whole subcarriers, which
// slides the entire spectrum by that integer `k`. DVB-T's 45 continual pilots sit
// at FIXED carrier positions on every symbol (§4.5.4) and are boosted (16/9
// power), so they anchor the integer offset: after fractional correction and
// symbol alignment, FFT one symbol and, for each trial shift `k`, sum the energy
// landing at the continual-pilot bins shifted by `k`. The `k` that maximizes that
// pilot-position energy is the integer CFO — the boosted pilots dominate any
// coincidental data energy at the same 45 positions. This is the DVB-T-native
// counterpart to the OFDM preamble path's training-symbol integer-CFO recovery
// (`sync::ofdm_sync`), which a preamble-less frame cannot use.

use crate::waveform::dvb_t::continual_pilot_bins;

/// One integer-CFO estimate: the offset in whole subcarrier spacings and a
/// confidence ratio.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct IntegerCfoResult {
    /// Integer CFO in subcarrier spacings (the spectrum is shifted by this many
    /// bins). `cfo_hz = bins · fs / n_fft`; correct by rotating the time domain by
    /// `−cfo_hz`.
    pub bins: i32,
    /// Confidence: the winning shift's continual-pilot energy divided by the mean
    /// over all trial shifts (`1.0` = no discrimination; higher is a clearer
    /// lock). Useful for thresholding whether an integer offset is present.
    pub confidence: f32,
}

/// Estimates the integer carrier-frequency offset of a DVB-T symbol from its
/// **frequency-domain** bins `freq` (one symbol's FFT output, `n_fft` long, in
/// rustfft bin order — i.e. already CP-removed, FFT'd, and ideally fractional-CFO
/// corrected). Searches trial shifts `k ∈ [−max_bins, max_bins]` and returns the
/// one maximizing the energy at the 45 continual-pilot bins shifted by `k`.
///
/// Returns `None` if `freq.len() < n_fft` or `max_bins == 0`. `max_bins` bounds
/// the search; the continual pilots span the active band, so shifts larger than
/// the guard-band margin slide pilots out of the band and lose discrimination —
/// a few tens of subcarriers is a generous front-end range.
pub fn dvb_t_integer_cfo(freq: &[C32], n_fft: usize, max_bins: i32) -> Option<IntegerCfoResult> {
    if freq.len() < n_fft || n_fft == 0 || max_bins <= 0 {
        return None;
    }
    let pilot_bins = continual_pilot_bins();

    // Pilot-position energy for a trial shift `k` (bins wrap mod n_fft).
    let energy_at = |k: i32| -> f32 {
        pilot_bins
            .iter()
            .map(|&b| {
                let idx = (b as i32 + k).rem_euclid(n_fft as i32) as usize;
                freq[idx].norm_sqr()
            })
            .sum()
    };

    let mut best_k = 0i32;
    let mut best_energy = f32::NEG_INFINITY;
    let mut sum_energy = 0.0f32;
    let mut count = 0u32;
    for k in -max_bins..=max_bins {
        let e = energy_at(k);
        sum_energy += e;
        count += 1;
        if e > best_energy {
            best_energy = e;
            best_k = k;
        }
    }

    let mean = sum_energy / count as f32;
    let confidence = if mean > 0.0 { best_energy / mean } else { 0.0 };
    Some(IntegerCfoResult {
        bins: best_k,
        confidence,
    })
}