Python Detection Statistics API¶

Threshold eta for a given false-alarm probability.

Exact closed-form inversion of Pfa = exp(-eta^2/2):

eta = sqrt(-2 * ln(pfa))

The threshold is independent of dwell and SNR; it depends only on the desired Pfa.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_threshold
>>> round(det_threshold(pfa=1e-6), 4)
5.2565

Detection probability for given per-sample amplitude SNR and dwell.

Computes Pd = Q_1(a, eta) where a = sqrt(2 * dwell) * snr.

At snr = 0, det_pd returns Pfa (the false-alarm rate, as expected for a noise-only input). As snr or dwell increase, Pd approaches 1.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_pd, det_threshold
>>> thr = det_threshold(pfa=1e-6)
>>> round(det_pd(snr=1.613, dwell=8, threshold=thr), 2)  # Pd 0.9
0.9
>>> round(det_pd(snr=0.0, dwell=8, threshold=thr), 6)    # Pd = Pfa
1e-06

Minimum dwell such that Pd >= pd_min for the given SNR and Pfa.

Iterates dwell = 1, 2, ..., max_dwell, computing det_pd() at each step. Returns the first dwell that satisfies the Pd requirement, or -1 if none is found within max_dwell iterations.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_dwell
>>> det_dwell(snr=0.5, pd_min=0.9, pfa=1e-6, max_dwell=256)
84

Minimum per-sample amplitude SNR achieving Pd >= pd_min.

Binary search over SNR in [0, hi] where hi is doubled from 1.0 until det_pd(hi, dwell, threshold) >= pd_min. 64 bisection iterations yield ~1e-19 relative precision on the final interval.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_snr, det_pd, det_threshold
>>> snr = det_snr(dwell=8, pd_min=0.9, pfa=1e-6)
>>> round(snr, 3)
1.613
>>> pd = det_pd(snr=snr, dwell=8, threshold=det_threshold(pfa=1e-6))
>>> abs(pd - 0.9) < 1e-9   # det_snr inverts det_pd, to tolerance
True

CFAR threshold eta_nc for a non-coherent detector of n_noncoh looks.

Solves marcum_q(n_noncoh, 0, eta_nc) = pfa (the order-M central tail, monotone decreasing in eta_nc) by bisection. For n_noncoh = 1 this is the exact closed form sqrt(-2 ln pfa) (== det_threshold).

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_threshold_noncoherent
>>> from doppler.detection import det_threshold
>>> round(det_threshold_noncoherent(pfa=1e-3, n_noncoh=4), 3)
5.111
>>> det_threshold_noncoherent(pfa=1e-6, n_noncoh=1) == det_threshold(
...     pfa=1e-6)
True

Detection probability for n_noncoh non-coherent looks.

Computes Pd = Q_{n_noncoh}(a, threshold) with the non-centrality a = sqrt(2 * n_coh * n_noncoh) * snr. At n_noncoh = 1 this is exactly det_pd(snr, n_coh, threshold); at snr = 0 it returns the per-test Pfa.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_pd_noncoherent, det_pd
>>> from doppler.detection import det_threshold_noncoherent
>>> from doppler.detection import det_threshold
>>> eta = det_threshold(pfa=1e-6)
>>> det_pd_noncoherent(snr=0.5, n_coh=8, n_noncoh=1, threshold=eta) \
...     == det_pd(snr=0.5, dwell=8, threshold=eta)  # -> coherent
True
>>> eta4 = det_threshold_noncoherent(pfa=1e-3, n_noncoh=4)
>>> round(det_pd_noncoherent(
...     snr=0.3, n_coh=16, n_noncoh=4, threshold=eta4), 2)
0.19

Minimum non-coherent looks achieving Pd >= pd_min at fixed n_coh.

Iterates n_noncoh = 1, 2, ..., max_n_noncoh, recomputing the threshold (det_threshold_noncoherent, which grows with the look count) at each step. Returns the first look count that meets the Pd requirement, or -1 if none does within max_n_noncoh. Used by the acquisition engine's (M, N_nc) split.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_n_noncoh
>>> det_n_noncoh(
...     snr=2.0, n_coh=16, pd_min=0.9, pfa=1e-3, max_n_noncoh=64)
1

EMA coefficient for a target estimator SNR (DC level in noise).

Sizes a first-order EMA y = (1-alpha)*y + alpha*x that estimates a DC level from noisy i.i.d. measurements x. Per sample the estimator SNR (mean^2 / variance) is snr_in; the EMA improves it by its variance reduction (2-alpha)/alpha, so the output SNR is snr_out = snr_in * (2-alpha)/alpha. Solving for the coefficient:

alpha = 2 * snr_in / (snr_in + snr_out) (SNRs linear)

Returns 1.0 (no averaging) when snr_out_db <= snr_in_db. Typical inputs: a signal-free power reference |n|^2 is exponential (0 dB per sample); a lock signal at known C/N0 has per-look SNR from its coherent integration (minus squaring loss), and this picks the smoothing bandwidth that makes the lock decision variable meet a chosen decision SNR.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_ema_alpha
>>> det_ema_alpha(0.0, 0.0)      # no gain requested -> no averaging
1.0
>>> round(1 / det_ema_alpha(0.0, 20.0), 1)   # 20 dB gain ~ 50 looks
50.5
>>> round(1 / det_ema_alpha(10.0, 30.0), 1)  # same 20 dB gain, shifted
50.5

Verify count: consecutive looks needed to compound to a budget.

n consecutive independent looks at per-look probability p compound to p^n, so the smallest n with p_look^n <= p_target is ceil(ln p_target / ln p_look) (clamped to >= 1). One function serves both sides of a lock detector (lockdet_core.h): the declare count from (per-look pfa, false-declare budget) and the drop count from (per-look miss rate 1 - pd, false-drop budget). Degenerate inputs resolve naturally: a target already met by one look returns 1; p_look >= 1 can never compound below a smaller target and returns INT_MAX.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_verify_count
>>> det_verify_count(1e-3, 1e-6)   # two 1e-3 looks reach 1e-6
2
>>> det_verify_count(1e-3, 1e-9)
3
>>> det_verify_count(0.5, 1e-3)    # drop side: pd = 0.5 per look
10
>>> det_verify_count(1e-3, 0.5)    # budget already met -> 1
1

Expected looks until a run of n consecutive successes completes.

The mean waiting time of the consecutive-run process a lockdet verify counter implements: at per-look success probability p, the first run of n straight successes takes on average

E[T] = (1 - p^n) / (p^n * (1 - p)) looks,

which is the declare latency bought by a verify count of n (multiply by the look period for time). Limits are handled exactly: p = 1 gives n (the run completes immediately), p = 0 gives infinity.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_verify_delay
>>> det_verify_delay(1.0, 8)             # certain hits: exactly n
8.0
>>> round(det_verify_delay(0.5, 2), 6)   # 2 straight coin heads: 6
6.0
>>> round(det_verify_delay(0.9, 8), 1)
13.2

Upper quantile of F(n, n) — the exact H0 law for a ratio test whose noise reference is estimated from as many samples as the signal sum.

A chi-square threshold (det_threshold_noncoherent) prices a statistic normalised by a KNOWN noise power. When the noise power is instead estimated from n same-burst samples (the BurstDespreader lock test: sum Re^2 against sum Im^2), the ratio's tail fattens to F(n, n) and the chi-square gate realizes tens of times the priced pfa (41x at n = 16, pfa = 1e-3). This helper returns the exact gate: P(chi2_n / chi2_n > g) = I_{1/(1+g)}(n/2, n/2) = pfa, solved on the regularized incomplete beta — valid for every n >= 1, odd included. As n grows the estimate hardens and g approaches the known-noise value. Threshold a BurstDespreader as lock_stat > sqrt(stat_n * det_threshold_f(pfa, stat_n)).

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_threshold_f
>>> round(det_threshold_f(1e-3, 2), 6)  # exact: (1 - pfa)/pfa
999.0
>>> round(det_threshold_f(1e-3, 4), 4)
53.4358
>>> round(det_threshold_f(1e-3, 64), 4)  # hardens toward known-noise
2.1931

Power threshold p from Pfa for the power detector.

Exact closed-form: P(Exponential(1) > p) = exp(-p) = Pfa, so

p = -ln(Pfa)

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_threshold_power
>>> round(det_threshold_power(pfa=1e-6), 3)   # -ln(1e-6) = 6*ln(10)
13.816

Detection probability for the power detector.

Pd = Q_1(sqrt(2·dwell·snr_power), sqrt(2·power_threshold))

The result equals det_pd() at the equivalent amplitude SNR: power SNR s corresponds to amplitude SNR sqrt(s), and the Q_1 arguments match.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_pd_power, det_threshold_power
>>> thr = det_threshold_power(pfa=1e-6)
>>> round(det_pd_power(
...     snr_power=2.6017, dwell=8, power_threshold=thr), 2)
0.9

Minimum dwell such that Pd >= pd_min for the power detector.

Parameters:

Returns:

Examples:

>>> from doppler.detection import det_dwell_power
>>> det_dwell_power(
...     snr_power=0.25, pd_min=0.9, pfa=1e-6, max_dwell=256)
84

Minimum per-sample power SNR achieving Pd >= pd_min.

Parameters:

Returns:

Examples:

>>> from doppler.detection import (det_snr_power, det_pd_power,
...                                det_threshold_power)
>>> sp = det_snr_power(dwell=8, pd_min=0.9, pfa=1e-6)
>>> round(sp, 4)
2.6017
>>> pd = det_pd_power(snr_power=sp, dwell=8,
...                   power_threshold=det_threshold_power(pfa=1e-6))
>>> abs(pd - 0.9) < 1e-9   # det_snr_power inverts det_pd_power
True

Marcum Q function Q_M(a, b) for integer M >= 1.

Probability that a Rice(a, sigma=1) random variable exceeds b. For M=1: Q_1(a, b) = P(Rice(a,1) > b). General integer M relates to the noncentral chi-squared CDF with 2M degrees of freedom.

Computed via the Poisson-weighted chi-squared series (exact for M=1, converges in ~60 terms for practical a, b <= 15):

Q_M(a, b) = sum_{k=0}^inf w_k * Q_{M+k}(0, b)

where: w_k = exp(-u) * u^k/k! (u = a^2/2) Q_n(0,b) = exp(-v) * sum_{j=0}^{n-1} v^j/j! (v = b^2/2)

Each iteration advances both the Poisson weight and the chi-sum in O(1) using the recurrences w_{k+1} = w_k * u/(k+1) and Q_{n+1}(0,b) = Q_n(0,b) + exp(-v)*v^n/n!. Total cost: O(K) where K ~ max(u, M) + safety margin.

Special cases:

• a = 0: Q_M(0, b) = exp(-b^2/2) * sum_{j=0}^{M-1} (b^2/2)^j/j!
• b <= 0: Q_M(a, b) = 1.0

Parameters:

Returns:

Examples:

>>> from doppler.detection import marcum_q
>>> round(marcum_q(m=1, a=0.0, b=1.0), 5)  # P(Rayleigh>1) = exp(-.5)
0.60653
>>> round(marcum_q(m=1, a=0.0, b=2.0), 5)   # exp(-2)
0.13534
>>> round(marcum_q(m=2, a=0.0, b=2.0), 5)   # 3*exp(-2)
0.40601
>>> round(marcum_q(m=1, a=2.0, b=1.0), 5)   # signal present (a=2)
0.91811