Monte Carlo vs Marcum Q Theory¶

import math

import numpy as np

from doppler.detection import (
    det_pd,
    det_pd_power,
    det_threshold,
    det_threshold_power,
    marcum_q,
)

N = 64  # frame length (complex samples)
DWELL = 4  # coherent integrations
SIGMA = 1.0  # noise std dev per real/imag component
PFA = 1e-2  # false-alarm rate (moderate for MC statistics)
N_TRIALS = 30_000  # trials per SNR point
SNR_HIST = 1.5  # SNR used for the survival-function panels

# Envelope detector thresholds.
ETA = det_threshold(PFA)  # Q₁ argument
THETA = ETA * math.sqrt(2.0 / math.pi)  # env_stat threshold

# Power detector thresholds.
P_THRESH = det_threshold_power(PFA)  # pow_stat threshold = -ln(Pfa)

RNG = np.random.default_rng(0)


def simulate(snr: float, n_trials: int = N_TRIALS) -> tuple:
    """Run Monte Carlo trials and return both test statistics.

    Samples the matched-filter output and noise reference directly from
    their theoretical distributions.  The FFT cross-correlation of a
    single-bin tone with AWGN produces R[τ] = W·exp(j2πτ/N), making
    |R[τ]| identical at every lag — so both test statistics would be 1
    regardless of SNR.  Direct sampling avoids that degeneracy and is
    the canonical way to verify the Marcum Q model.

    After M coherent integrations:
    - Signal channel: R₀ ~ CN(M·A, 2·M·σ²/N)
    - Noise reference: N−1 i.i.d. samples from CN(0, 2·M·σ²/N)

    Parameters
    ----------
    snr : float
        Post-correlation amplitude SNR: snr = A·√(N/2)/σ.  Zero → H0.
    n_trials : int
        Number of independent trials.

    Returns
    -------
    env_stat : np.ndarray, shape (n_trials,)
        |R₀| / mean(|noise_ref|).
    pow_stat : np.ndarray, shape (n_trials,)
        |R₀|² / mean(|noise_ref|²).
    """
    A = snr * SIGMA * math.sqrt(2.0 / N)
    sig_amp = DWELL * A  # coherently accumulated
    sigma_ac = SIGMA * math.sqrt(DWELL / N)  # per-component std dev

    # Matched-filter output at lag 0: CN(sig_amp, 2·sigma_ac²).
    re0 = sig_amp + RNG.standard_normal(n_trials) * sigma_ac
    im0 = RNG.standard_normal(n_trials) * sigma_ac
    mag0 = np.hypot(re0, im0)

    # N−1 independent noise reference samples from CN(0, 2·sigma_ac²).
    re_n = RNG.standard_normal((n_trials, N - 1)) * sigma_ac
    im_n = RNG.standard_normal((n_trials, N - 1)) * sigma_ac
    mag_n = np.hypot(re_n, im_n)

    env_stat = mag0 / mag_n.mean(axis=1)
    pow_stat = mag0**2 / (mag_n**2).mean(axis=1)
    return env_stat, pow_stat



print(f"Running H0 ({N_TRIALS:,} trials)…")
env_h0, pow_h0 = simulate(0.0)

print(f"Running H1 snr={SNR_HIST} ({N_TRIALS:,} trials)…")
env_h1, pow_h1 = simulate(SNR_HIST)

snr_sweep = np.linspace(0.0, 3.0, 25)
print(f"Running Pd sweep ({len(snr_sweep)} SNR × {N_TRIALS:,} trials)…")
pd_env_mc = np.empty(len(snr_sweep))
pd_pow_mc = np.empty(len(snr_sweep))
for i, s in enumerate(snr_sweep):
    ev, pw = simulate(float(s))
    pd_env_mc[i] = (ev > THETA).mean()
    pd_pow_mc[i] = (pw > P_THRESH).mean()

pd_env_th = np.array([det_pd(float(s), DWELL, ETA) for s in snr_sweep])
pd_pow_th = np.array(
    [det_pd_power(float(s) ** 2, DWELL, P_THRESH) for s in snr_sweep]
)

print("Done.")

# Validate MC against theory.
# CFAR design point: both detectors' H0 exceedance rates must sit at
# the Pfa their thresholds were derived for (30k trials → MC sigma
# ≈ 6e-4, so ±50% of Pfa = 1e-2 is a generous band).
pfa_env_mc = float((env_h0 > THETA).mean())
pfa_pow_mc = float((pow_h0 > P_THRESH).mean())
print(f"Pfa: env = {pfa_env_mc:.4f}, pow = {pfa_pow_mc:.4f} (design {PFA})")
assert abs(pfa_env_mc - PFA) < 0.5 * PFA, "envelope Pfa off design point"
assert abs(pfa_pow_mc - PFA) < 0.5 * PFA, "power Pfa off design point"

# Marcum-Q Pd model: the MC sweep must track theory at every SNR.
err_env = float(np.max(np.abs(pd_env_mc - pd_env_th)))
err_pow = float(np.max(np.abs(pd_pow_mc - pd_pow_th)))
print(f"max |MC - theory| Pd: env = {err_env:.4f}, pow = {err_pow:.4f}")
assert err_env < 0.02, "envelope Pd deviates from Marcum-Q theory"
assert err_pow < 0.02, "power Pd deviates from Marcum-Q theory"