2-D Acquisition Grid¶

import math

import numpy as np

from doppler.detection import det_dwell, det_pd, det_threshold
from doppler.spectral import Corr2D, CorrDetector2D

N_DOPPLER = 16  # Doppler search bins  (rows)
N_CODE_PHASE = 16  # code-phase search bins (columns, = PRN length in chips)
N = N_DOPPLER * N_CODE_PHASE  # total cells

# True signal location in the search grid.
# Flat index 0 → noise reference uses all other N-1=255 cells, maximising
# the reference population and minimising CFAR scalloping loss.
DOPPLER_BIN_TRUE = 5  # true Doppler bin
CODE_PHASE_BIN_TRUE = 11  # true code-phase offset (chips)

SNR_DB = 3.0  # per-sample amplitude SNR (dB)
SIGMA = 1.0  # noise std dev per real/imag component

PFA = 1e-3  # false-alarm probability (per cell per dwell)
PD_MIN = 0.90  # minimum detection probability requirement
MAX_DWELL = 64  # upper search limit for det_dwell()

N_TRIALS = 3_000  # MC trials
RNG = np.random.default_rng(0)

# Theory: minimum dwell and CFAR threshold.
snr_amp = 10.0 ** (SNR_DB / 20.0)

# Bonferroni correction: N cells tested per dwell → per-cell Pfa must be
# much tighter so that the system (max-of-N) Pfa meets the target.
pfa_cell = 1.0 - (1.0 - PFA) ** (1.0 / N)

# det_threshold() returns η (Marcum Q argument).
# The gate on test_stat = peak_mag / noise_est is θ = η · √(2/π).
eta = det_threshold(pfa_cell)
theta = eta * math.sqrt(2.0 / math.pi)

M = det_dwell(snr_amp, PD_MIN, pfa_cell, MAX_DWELL)
if M <= 0:
    raise RuntimeError(
        f"SNR = {SNR_DB} dB is not achievable at Pd ≥ {PD_MIN}, "
        f"Pfa_sys = {PFA:.0e} within {MAX_DWELL} dwells."
    )

pd_cell = det_pd(snr_amp, M, eta)
pd_theory = 1.0 - (1.0 - pd_cell) * (1.0 - pfa_cell) ** (N - 1)

print(
    f"Grid            : {N_DOPPLER} Doppler × {N_CODE_PHASE} code-phase  "
    f"(N = {N} cells)"
)
print(f"Amplitude SNR   : {SNR_DB:+.1f} dB  (linear {snr_amp:.3f})")
print(f"Pfa_sys target  : {PFA:.0e}  →  pfa_cell = {pfa_cell:.2e}")
print(f"                        →  η = {eta:.4f},  θ = {theta:.4f}")
print(f"Pd target       : {PD_MIN:.2f}   →  required dwell M = {M}")
print(f"Theory Pd @ M   : {pd_theory:.4f}  (cell: {pd_cell:.4f})")

# Build a CAZAC-style (flat-spectrum) reference by assigning a random phase
# to each 2D FFT bin and synthesizing via IFFT.  Flat spectrum → exactly
# zero circular autocorrelation outside lag (0,0).  This prevents the
# coherent ±1 BPSK sidelobe contamination that would otherwise inflate the
# CFAR noise estimate over M dwells and cause a systematic ~3% Pd loss.
#
# In a GPS/CDMA receiver the 2D ref template encodes each (Doppler, code-
# phase) hypothesis; using a flat-power-spectrum template is equivalent to
# an acquisition grid with orthogonal hypotheses (like a CAZAC / Zadoff-Chu).

phases_spec = RNG.uniform(0, 2.0 * np.pi, (N_DOPPLER, N_CODE_PHASE)).astype(
    np.float32
)
# Unit-amplitude spectrum → IFFT synthesizes the flat-autocorrelation template.
# numpy ifft2 normalises by 1/N, giving ||ref2d||² = 1; we scale by sqrt(N) so
# that the matched-filter peak (= ||ref2d||² × N) equals N, matching the SNR
# formula A = snr · σ / √N (signal per frame at true cell = A·N).
ref2d = (np.sqrt(N) * np.fft.ifft2(np.exp(1j * phases_spec))).astype(
    np.complex64
)

# Correlator normalizes by N per frame, so the per-frame amplitude at the
# true cell equals A exactly.  The per-component noise at the output is
# σ / √(2N) per frame.  The amplitude SNR seen by det_pd is therefore:
#     snr = A · √N / σ   →   A = snr · σ / √N
A = snr_amp * SIGMA / math.sqrt(N)
ns = np.float32(SIGMA / math.sqrt(2.0))  # per-component noise std


def signal_frame() -> np.ndarray:
    """One CF32 acquisition frame with signal at
    (DOPPLER_BIN_TRUE, CODE_PHASE_BIN_TRUE).
    """
    sig = np.roll(
        np.roll(ref2d * A, DOPPLER_BIN_TRUE, axis=0),
        CODE_PHASE_BIN_TRUE,
        axis=1,
    )
    noise = (
        RNG.standard_normal((N_DOPPLER, N_CODE_PHASE))
        + 1j * RNG.standard_normal((N_DOPPLER, N_CODE_PHASE))
    ).astype(np.complex64) * ns
    return sig + noise

signal_block = np.concatenate([signal_frame().ravel() for _ in range(M)])

det = CorrDetector2D(
    ref2d, dwell=M, noise_lo=1, noise_hi=N - 1, threshold=0.0
)
for *_, stat in det.push(signal_block):
    detected = stat > theta