File detection_core.h¶

FileList > detection > detection_core.h

Go to the source code of this file

Detection-theory utilities for the amplitude-ratio test statistic. More...

Public Functions¶

Detailed Description¶

The doppler detector forms the test statistic:

test_stat = peak_mag / noise_est

With M-point coherent integration (dwell = M) and per-sample amplitude SNR snr (signal amplitude / noise amplitude, linear):

Under H0 (noise only): test_stat ~ Rayleigh(1) Under H1 (signal+noise): test_stat ~ Rice(a, 1), a = sqrt(2*M) * snr

False-alarm probability (threshold-only, M-independent):

Pfa = exp(-eta^2/2) => eta = sqrt(-2 ln Pfa) (exact)

Detection probability:

Pd = Q_1(a, eta) (Marcum Q function, order 1)

All functions are stateless and thread-safe.

Public Functions Documentation¶

function det_dwell¶

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

int det_dwell (
    double snr,
    double pd_min,
    double pfa,
    int max_dwell
) 

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:

• snr Per-sample amplitude SNR (linear).
• pd_min Required detection probability, e.g. 0.9.
• pfa False-alarm probability; used to derive eta.
• max_dwell Search upper bound; prevents infinite loops for low SNR.

Returns:

Minimum dwell >= 1, or -1 if not achievable.

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

function det_dwell_gauss¶

Looks a Gaussian statistic must average to separate H1 from H0.

int det_dwell_gauss (
    double mean,
    double var,
    double pd,
    double pfa
) 

The classic sizing: with a per-look H0 variance var and an H1 mean mean (H0 mean zero), block-averaging n looks shrinks the H0 spread as 1/n, and the smallest n whose H0 and H1 tails clear both budgets is

n = var * ((Q_inv(pfa) - Q_inv(pd)) / mean)^2

Q_inv(pd) is negative for pd > 0.5, so the difference is the total separation both tails must fit inside.

Parameters:

• mean H1 mean of one look, > 0 (H0 mean is taken as zero).
• var H0 variance of one look, > 0.
• pd Required detection probability, in (0, 1).
• pfa Allowed false-alarm probability, in (0, 1) and below pd.

Returns:

Looks needed, rounded up and clamped to >= 1; -1 on invalid input.

>>> from doppler.detection import det_dwell_gauss
>>> det_dwell_gauss(mean=0.4, var=0.5, pd=0.99, pfa=1e-5)
136
>>> det_dwell_gauss(mean=0.8, var=0.5, pd=0.99, pfa=1e-5)   # 2x mean
34
>>> det_dwell_gauss(mean=0.0, var=0.5, pd=0.99, pfa=1e-5)   # no signal
-1

function det_dwell_power¶

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

int det_dwell_power (
    double snr_power,
    double pd_min,
    double pfa,
    int max_dwell
) 

Parameters:

• snr_power Per-sample power SNR (linear).
• pd_min Required detection probability.
• pfa False-alarm probability; used to derive p.
• max_dwell Search upper bound.

Returns:

Minimum dwell >= 1, or -1 if not achievable.

>>> 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

function det_ema_alpha¶

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

double det_ema_alpha (
    double snr_in_db,
    double snr_out_db
) 

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:

• snr_in_db Per-sample estimator SNR, dB (mean^2 / variance).
• snr_out_db Desired EMA-output estimator SNR, dB.

Returns:

EMA coefficient alpha in (0, 1].

>>> 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

function det_n_noncoh¶

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

int det_n_noncoh (
    double snr,
    int n_coh,
    double pd_min,
    double pfa,
    int max_n_noncoh
) 

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:

• snr Per-sample amplitude SNR (linear).
• n_coh Coherent integration length in samples (dwell * N).
• pd_min Required detection probability, e.g. 0.9.
• pfa Per-test false-alarm probability.
• max_n_noncoh Search upper bound on the look count.

Returns:

Minimum n_noncoh >= 1, or -1 if not achievable.

>>> 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

function det_pd¶

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

double det_pd (
    double snr,
    int dwell,
    double threshold
) 

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:

• snr Per-sample amplitude SNR (signal / noise amplitude, linear). snr = 0 gives Pd = Pfa.
• dwell Coherent integration depth; must be >= 1.
• threshold Test-stat threshold eta, e.g. from det_threshold().

Returns:

Detection probability in [0, 1].

>>> 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

function det_pd_noncoherent¶

Detection probability for n_noncoh non-coherent looks.

double det_pd_noncoherent (
    double snr,
    int n_coh,
    int n_noncoh,
    double threshold
) 

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:

• snr Per-sample amplitude SNR (signal / noise amplitude).
• n_coh Coherent integration length in samples (dwell * N).
• n_noncoh Number of non-coherent looks; must be >= 1.
• threshold Threshold eta_nc, e.g. from det_threshold_noncoherent().

Returns:

Detection probability in [0, 1].

>>> 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

function det_pd_power¶

Detection probability for the power detector.

double det_pd_power (
    double snr_power,
    int dwell,
    double power_threshold
) 

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

Parameters:

• snr_power Per-sample power SNR (signal power / noise power at the correlator output, linear). 0 gives Pd = Pfa.
• dwell Coherent integration depth; must be >= 1.
• power_threshold Threshold p, e.g. from det_threshold_power().

Returns:

Detection probability in [0, 1].

>>> 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

function det_q_inv¶

Upper-tail quantile of the standard normal: the eta with Q(eta) = p.

double det_q_inv (
    double p
) 

Q(eta) = 0.5*erfc(eta/sqrt(2)), so this is sqrt(2)*erfcinv(2p). Everything below is expressed in it, and a caller thresholding its own zero-mean Gaussian statistic wants det_q_inv(pfa) * sd_H0.

Signed, and that matters. Above the median the quantile is negative, which is exactly why det_dwell_gauss()'s Q_inv(pfa) - Q_inv(pd) is a sum of two tails rather than a difference: every caller's pd is above 0.5. Clamping it to zero there halves the dwell without failing anything.

Parameters:

• p Tail probability in (0, 1).

Returns:

Quantile in H0 sigmas positive below the median, exactly 0 at it, negative above. Fails closed (0.0) for p outside (0, 1).

>>> from doppler.detection import det_q_inv, det_threshold
>>> round(det_q_inv(p=5e-6), 4)     # the carrier lock metric's 4.42 sigma
4.4172
>>> round(det_q_inv(p=0.5), 4)      # the median
0.0
>>> round(det_q_inv(p=0.99), 4)     # above it: NEGATIVE, by design
-2.3263
>>> round(det_threshold(pfa=5e-6), 4)   # the OTHER law -- not this one
4.9409

function det_snr¶

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

double det_snr (
    int dwell,
    double pd_min,
    double pfa
) 

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:

• dwell Coherent integration depth; must be >= 1.
• pd_min Required detection probability.
• pfa False-alarm probability; used to derive eta.

Returns:

Minimum amplitude SNR >= 0.

>>> 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

function det_snr_power¶

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

double det_snr_power (
    int dwell,
    double pd_min,
    double pfa
) 

Parameters:

• dwell Coherent integration depth; must be >= 1.
• pd_min Required detection probability.
• pfa False-alarm probability.

Returns:

Minimum power SNR >= 0.

>>> 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

function det_threshold¶

Threshold eta for a given false-alarm probability.

double det_threshold (
    double pfa
) 

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:

• pfa Desired false-alarm probability; must be in (0, 1).

Returns:

Threshold eta > 0.

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

function det_threshold_f¶

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.

double det_threshold_f (
    double pfa,
    int n
) 

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:

• pfa Tail probability budget, in (0, 1).
• n Degrees of freedom on each side (>= 1).

Returns:

The F(n, n) upper-pfa quantile; 0 on invalid input.

>>> 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

function det_threshold_gauss¶

Declare threshold for a Gaussian statistic sized by det_dwell_gauss.

double det_threshold_gauss (
    double mean,
    double pd,
    double pfa
) 

The crossover point that meets both budgets at once, in the statistic's own units:

thresh = Q_inv(pfa) * mean / (Q_inv(pfa) - Q_inv(pd))

Independent of the variance and of the look count those set how many looks are needed to reach this point, not where it is.

Parameters:

• mean H1 mean of one look, > 0.
• pd Required detection probability, in (0, 1).
• pfa Allowed false-alarm probability, in (0, 1) and below pd.

Returns:

Threshold in the statistic's units; 0.0 on invalid input.

>>> from doppler.detection import det_threshold_gauss
>>> round(det_threshold_gauss(mean=0.4, pd=0.99, pfa=1e-5), 4)
0.2588
>>> round(det_threshold_gauss(mean=0.8, pd=0.99, pfa=1e-5), 4)  # scales
0.5176

function det_threshold_noncoherent¶

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

double det_threshold_noncoherent (
    double pfa,
    int n_noncoh
) 

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:

• pfa Per-test false-alarm probability in (0, 1).
• n_noncoh Number of non-coherent looks; must be >= 1.

Returns:

Threshold eta_nc on the normalized statistic R.

>>> 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

function det_threshold_power¶

Power threshold p from Pfa for the power detector.

double det_threshold_power (
    double pfa
) 

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

p = -ln(Pfa)

Parameters:

• pfa Desired false-alarm probability; must be in (0, 1).

Returns:

Threshold p > 0.

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

function det_verify_count¶

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

int det_verify_count (
    double p_look,
    double p_target
) 

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).

That ~ is a BUDGET, and deliberately the conservative side of one: a consecutive-run detector's exact declare rate is p^n (1-p)/(1-p^n) (lockdet_core.h), which is lower, so sizing on p^n over-provisions n rather than under. The gap is ~p negligible where a detector is really sized, 10% at p = 0.1 so pick n here and predict what a caller will observe with det_verify_delay().

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:

• p_look Per-look probability (pfa or 1 - pd), in (0, 1).
• p_target Compound probability budget, in (0, 1).

Returns:

Smallest verify count n with p_look^n <= p_target.

>>> 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

function det_verify_delay¶

Expected looks until a run of n consecutive successes completes.

double det_verify_delay (
    double p_look,
    int n
) 

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:

• p_look Per-look success probability (e.g. pd), in [0, 1].
• n Run length (the verify count); clamped to >= 1.

Returns:

Expected number of looks to the first length-n run.

>>> 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

function marcum_q¶

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

double marcum_q (
    int m,
    double a,
    double b
) 

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):

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!.

The window is CENTRED on the Poisson mode k ~ u = a^2/2 and its half-width scales as 12*sqrt(u+1) + 60 terms, so the term count grows with a rather than being the fixed ~60 this comment used to claim: about 60 terms at a = 0, but ~187 at a = 15. That scaling is the whole point a Poisson(u) distribution's mass sits at k ~ u with spread ~sqrt(u), so a fixed window anchored at k = 0 misses it entirely once a is large, which is a real bug this code already carries a comment about (see marcum_q.c). Total cost: O(sqrt(u) + M).

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:

• m Integration order; must be >= 1.
• a Non-centrality parameter (signal strength). a = 0 for H0.
• b Threshold (same units as test_stat).

Returns:

Q_M(a, b) in [0, 1].

>>> 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

The documentation for this class was generated from the following file native/inc/detection/detection_core.h