"""Tests for the detection module. The important ones are the Monte Carlo pins: every closed-form Pfa/Pd in the module is checked against a simulation of the detector it describes, per the arc's standing rule that a detection-statistic distribution quoted from memory ships only after a numeric check. """ import numpy as np import pytest from scipy.stats import chi2, norm from detection import (block_energies, ca_cfar, cfar_factor, cfar_pd, energy_detector_pd, energy_detector_threshold, incoherent_pd, matched_filter_pd, roc_empirical) class TestMatchedFilterPd: def test_monte_carlo_pins_the_closed_form(self): # Correlator statistic sum(y s) against its exact threshold. rng = np.random.default_rng(3) n, sigma = 64, 1.3 s = 0.6 * np.sin(2 * np.pi * 0.11 * np.arange(n)) es = np.sum(s**2) enr = es / sigma**2 pfa = 0.05 gamma = sigma * np.sqrt(es) * norm.isf(pfa) trials = 40000 w = sigma * rng.standard_normal((trials, n)) t0 = w @ s t1 = (w + s) @ s assert abs(np.mean(t0 > gamma) - pfa) < 0.005 assert abs(np.mean(t1 > gamma) - matched_filter_pd(pfa, enr)) < 0.01 def test_pd_at_zero_enr_is_pfa(self): assert abs(matched_filter_pd(0.1, 0.0) - 0.1) < 1e-12 def test_pd_is_monotone_in_enr_and_pfa(self): assert matched_filter_pd(0.01, 16.0) > matched_filter_pd(0.01, 4.0) assert matched_filter_pd(0.1, 4.0) > matched_filter_pd(0.01, 4.0) def test_shape_does_not_matter_only_energy(self): # Two very different waveforms with equal energy, each run # through the actual correlator: the measured Pd of both must # sit on the one closed form that only sees the energy. rng = np.random.default_rng(16) n, pfa, es = 64, 0.05, 12.0 chirp = np.sin(2 * np.pi * (0.05 + 0.2 * np.arange(n) / n) * np.arange(n)) front = np.zeros(n) front[:4] = 1.0 # all energy in 4 samples closed = matched_filter_pd(pfa, es) gamma = np.sqrt(es) * norm.isf(pfa) for shape in (chirp, front): s = shape * np.sqrt(es / np.sum(shape**2)) t1 = (rng.standard_normal((30000, n)) + s) @ s assert abs(np.mean(t1 > gamma) - closed) < 0.012 def test_invalid_inputs_raise(self): with pytest.raises(ValueError): matched_filter_pd(0.0, 1.0) with pytest.raises(ValueError): matched_filter_pd(1.0, 1.0) with pytest.raises(ValueError): matched_filter_pd(0.1, -1.0) class TestEnergyDetector: def test_threshold_hits_the_target_pfa(self): rng = np.random.default_rng(4) n, sigma2, pfa = 32, 1.7, 0.02 gamma = energy_detector_threshold(pfa, n, sigma2) e = np.sum(sigma2 * rng.standard_normal((60000, n))**2, axis=1) assert abs(np.mean(e > gamma) - pfa) < 0.003 def test_monte_carlo_pins_pd(self): rng = np.random.default_rng(5) n, sigma = 32, 1.0 s = 0.45 * np.cos(2 * np.pi * 0.17 * np.arange(n) + 0.3) enr = np.sum(s**2) / sigma**2 pfa = 0.05 gamma = energy_detector_threshold(pfa, n, sigma**2) e1 = np.sum((s + sigma * rng.standard_normal((40000, n)))**2, axis=1) assert abs(np.mean(e1 > gamma) - energy_detector_pd(pfa, enr, n)) < 0.01 def test_always_below_the_matched_filter(self): # Same pfa, same energy: not knowing the waveform never helps. for n in (8, 64, 512): for enr in (1.0, 9.0, 25.0): assert (energy_detector_pd(0.01, enr, n) < matched_filter_pd(0.01, enr)) def test_gap_grows_with_window_length(self): # At fixed enr the energy detector decays toward pfa as n # grows (more noise chi-square mass, same signal energy). pds = [energy_detector_pd(0.01, 9.0, n) for n in (8, 64, 512)] assert pds[0] > pds[1] > pds[2] def test_zero_enr_is_pfa(self): assert abs(energy_detector_pd(0.03, 0.0, 16) - 0.03) < 1e-12 def test_invalid_inputs_raise(self): with pytest.raises(ValueError): energy_detector_threshold(0.5, 0) with pytest.raises(ValueError): energy_detector_threshold(0.5, 8, 0.0) with pytest.raises(ValueError): energy_detector_pd(0.5, 1.0, 0) class TestIncoherentPd: def test_rayleigh_pfa_closed_form(self): # Threshold from chi2_isf(pfa, 2) must equal the exp form # gamma^2 = -2 sigma_z^2 ln(pfa). pfa = 0.01 assert abs(chi2.isf(pfa, 2) - (-2 * np.log(pfa))) < 1e-9 def test_monte_carlo_pins_pd(self): rng = np.random.default_rng(6) sz, a, pfa = 0.8, 2.4, 0.02 lam = a**2 / sz**2 gamma2 = sz**2 * chi2.isf(pfa, 2) trials = 60000 z0 = sz * (rng.standard_normal(trials) + 1j * rng.standard_normal(trials)) z1 = z0 + a assert abs(np.mean(np.abs(z0)**2 > gamma2) - pfa) < 0.003 assert abs(np.mean(np.abs(z1)**2 > gamma2) - incoherent_pd(pfa, lam)) < 0.01 def test_below_coherent_detection_at_same_lambda(self): # Knowing the phase (coherent, 1-D Gaussian test with # deflection sqrt(lam)) always beats thresholding |z|. for lam in (4.0, 9.0, 25.0): assert incoherent_pd(0.01, lam) < matched_filter_pd(0.01, lam) def test_invalid_inputs_raise(self): with pytest.raises(ValueError): incoherent_pd(0.0, 1.0) with pytest.raises(ValueError): incoherent_pd(0.1, -0.1) class TestCfarFactor: def test_pfa_is_constant_across_noise_levels(self): # The CFAR property itself: same alpha, same measured Pfa at # wildly different sigma. This is the module's central claim. rng = np.random.default_rng(7) m, n_ref, pfa = 32, 16, 0.02 alpha = cfar_factor(pfa, m, n_ref) trials = 40000 for sigma in (0.1, 1.0, 10.0): cut = np.sum(sigma * rng.standard_normal((trials, m))**2 * sigma, axis=1) ref = np.sum(sigma**2 * rng.standard_normal( (trials, n_ref * m))**2, axis=1) / n_ref assert abs(np.mean(cut > alpha * ref) - pfa) < 0.004 def test_reduces_to_classic_radar_form_for_two_dof(self): # For n_block = 2 (square-law / exponential cells) the F # inverse survival function must reproduce the textbook # CA-CFAR result Pfa = (1 + alpha/N)^(-N). for pfa in (1e-1, 1e-2, 1e-4): for n_ref in (4, 16, 64): alpha = cfar_factor(pfa, 2, n_ref) assert abs((1 + alpha / n_ref)**(-n_ref) - pfa) < 1e-10 def test_approaches_known_sigma_threshold_as_nref_grows(self): # Infinite reference cells = perfectly known noise level: the # implied per-block threshold alpha * E[E_ref] must approach # the chi-square threshold of the clairvoyant detector. m, pfa = 64, 1e-3 known = chi2.isf(pfa, m) implied = [cfar_factor(pfa, m, n_ref) * m for n_ref in (8, 64, 512)] gaps = [abs(x - known) for x in implied] assert gaps[0] > gaps[1] > gaps[2] assert gaps[2] / known < 0.02 def test_invalid_inputs_raise(self): with pytest.raises(ValueError): cfar_factor(0.0, 8, 8) with pytest.raises(ValueError): cfar_factor(0.1, 0, 8) with pytest.raises(ValueError): cfar_factor(0.1, 8, 0) class TestCfarPd: def test_monte_carlo_pins_pd(self): rng = np.random.default_rng(8) m, n_ref, pfa = 32, 16, 0.05 sigma = 1.4 s = 0.55 * np.cos(2 * np.pi * 0.21 * np.arange(m)) enr = np.sum(s**2) / sigma**2 alpha = cfar_factor(pfa, m, n_ref) trials = 40000 cut = np.sum((s + sigma * rng.standard_normal((trials, m)))**2, axis=1) ref = np.sum(sigma**2 * rng.standard_normal( (trials, n_ref * m))**2, axis=1) / n_ref assert abs(np.mean(cut > alpha * ref) - cfar_pd(pfa, enr, m, n_ref)) < 0.01 def test_cfar_loss_is_positive_and_shrinks_with_nref(self): # At the same pfa and enr, CFAR detects less than the # known-sigma energy detector, and the gap closes as the # reference window grows. m, pfa, enr = 64, 1e-3, 40.0 ideal = energy_detector_pd(pfa, enr, m) pds = [cfar_pd(pfa, enr, m, n_ref) for n_ref in (2, 8, 64)] assert all(pd < ideal for pd in pds) assert pds[0] < pds[1] < pds[2] assert ideal - pds[2] < 0.02 def test_zero_enr_is_pfa(self): assert abs(cfar_pd(0.03, 0.0, 16, 8) - 0.03) < 1e-12 class TestRocEmpirical: def test_identical_distributions_lie_on_the_chance_line(self): rng = np.random.default_rng(9) pfa, pd = roc_empirical(rng.standard_normal(20000), rng.standard_normal(20000)) assert np.max(np.abs(pd - pfa)) < 0.03 def test_separated_distributions_reach_the_corner(self): rng = np.random.default_rng(10) pfa, pd = roc_empirical(rng.standard_normal(5000), rng.standard_normal(5000) + 10.0) assert np.any((pd > 0.999) & (pfa < 0.001)) def test_matches_gaussian_closed_form(self): # The empirical curve of the mean-shift problem must sit on # Pd = Q(Q^-1(Pfa) - d). rng = np.random.default_rng(11) d = 2.0 pfa, pd = roc_empirical(rng.standard_normal(60000), rng.standard_normal(60000) + d) mid = (pfa > 0.01) & (pfa < 0.5) closed = norm.sf(norm.isf(pfa[mid]) - d) assert np.max(np.abs(pd[mid] - closed)) < 0.02 def test_empty_raises(self): with pytest.raises(ValueError): roc_empirical(np.array([]), np.ones(4)) class TestBlockEnergies: def test_partitions_and_sums_squares(self): x = np.arange(7, dtype=float) e = block_energies(x, 3) assert e.shape == (2,) assert abs(e[0] - (0 + 1 + 4)) < 1e-12 assert abs(e[1] - (9 + 16 + 25)) < 1e-12 def test_mean_removal_loses_one_dof(self): # Mean-removed block energies are sigma^2 chi2_{M-1}: their # mean is sigma^2 (M-1), not sigma^2 M, and the M-1 threshold # hits the target Pfa where the M one misses low. rng = np.random.default_rng(12) m, sigma2, pfa = 16, 2.3, 0.05 e = block_energies(sigma2**0.5 * rng.standard_normal(m * 40000), m, remove_mean=True) assert abs(e.mean() / (sigma2 * (m - 1)) - 1) < 0.01 gamma = energy_detector_threshold(pfa, m - 1, sigma2) assert abs(np.mean(e > gamma) - pfa) < 0.004 gamma_wrong = energy_detector_threshold(pfa, m, sigma2) assert np.mean(e > gamma_wrong) < 0.8 * pfa def test_dc_pedestal_is_removed(self): e_raw = block_energies(np.full(32, 5.0), 8) e_ac = block_energies(np.full(32, 5.0), 8, remove_mean=True) assert e_raw[0] > 100.0 and np.all(e_ac < 1e-20) def test_invalid_inputs_raise(self): with pytest.raises(ValueError): block_energies(np.ones(4), 0) with pytest.raises(ValueError): block_energies(np.ones(4), 8) with pytest.raises(ValueError): block_energies(np.ones(4), 1, remove_mean=True) class TestCaCfar: def test_warmup_produces_no_verdict(self): det, thr = ca_cfar(np.ones(20), alpha=2.0, n_ref=8, n_guard=2) assert not det[:10].any() assert np.isnan(thr[:10]).all() assert np.isfinite(thr[10:]).all() def test_detects_a_burst_and_only_the_burst(self): rng = np.random.default_rng(13) m, n_ref, pfa = 32, 16, 1e-3 x = rng.standard_normal(m * 200) burst = slice(m * 100, m * 101) x[burst] += 1.2 * np.sin(2 * np.pi * 0.13 * np.arange(m)) e = block_energies(x, m) det, _ = ca_cfar(e, cfar_factor(pfa, m, n_ref), n_ref, n_guard=1) assert det[100] assert det.sum() <= 2 # the burst, at most one false alarm def test_pfa_holds_when_noise_level_steps(self): # The reason CFAR exists: after the noise floor jumps 10 dB # and the reference window refills, the false-alarm rate is # the same as before the jump. A fixed threshold set for the # quiet half would fire constantly in the loud half. rng = np.random.default_rng(14) m, n_ref, pfa = 32, 16, 0.01 alpha = cfar_factor(pfa, m, n_ref) n_blocks = 30000 sigma = np.ones(n_blocks) sigma[n_blocks // 2:] = 10**0.5 # +10 dB in power x = np.repeat(sigma, m) * rng.standard_normal(n_blocks * m) e = block_energies(x, m) det, _ = ca_cfar(e, alpha, n_ref, n_guard=1) quiet = det[n_ref + 1:n_blocks // 2] loud = det[n_blocks // 2 + n_ref + 1:] assert abs(quiet.mean() - pfa) < 0.004 assert abs(loud.mean() - pfa) < 0.004 gamma_fixed = energy_detector_threshold(pfa, m, 1.0) assert np.mean(e[n_blocks // 2:] > gamma_fixed) > 0.5 def test_guard_blocks_shield_the_reference(self): # A signal spanning two blocks: with a guard the second block # is judged against clean noise; without, the first block # contaminates its own reference estimate is avoided only by # the guard. Assert the guard version's threshold for the # second burst block is not inflated by the first. rng = np.random.default_rng(15) m, n_ref = 32, 8 x = rng.standard_normal(m * 40) x[m * 20:m * 22] += 3.0 # two consecutive loud blocks e = block_energies(x, m) _, thr_guard = ca_cfar(e, 2.0, n_ref, n_guard=1) _, thr_no = ca_cfar(e, 2.0, n_ref, n_guard=0) # For the second burst block (index 21), the no-guard # reference includes loud block 20; the guarded one does not. assert thr_no[21] > 1.5 * thr_guard[21] def test_invalid_inputs_raise(self): with pytest.raises(ValueError): ca_cfar(np.ones(10), 0.0, 4) with pytest.raises(ValueError): ca_cfar(np.ones(10), 1.0, 0) with pytest.raises(ValueError): ca_cfar(np.ones(10), 1.0, 4, n_guard=-1)