"""Integer CUSUM drift monitor: a bit-exact model of the embedded C. The [embedded companion](embedded.qmd) puts a CUSUM on the block-energy stream a CFAR burst detector is already computing, so that a slow change in the noise floor -- the very thing the CFAR is designed to be immune to -- stops being invisible. On the target hardware that detector is integer arithmetic end to end. This module is the reference implementation of exactly that integer pipeline, in Python, using only operations the C can perform: uint64 block energy -> ln() from the bit length plus a Q12 Horner polynomial -> Q8 standardization against a stored baseline -> Q8 two-sided CUSUM with a saturating accumulator -> re-baselining from a ring buffer of recent log-energies, stored Q12 >> 4 because a raw Q12 log-energy is about 2.3x too large for the int16 ring the firmware can afford It exists so the page's claims about the firmware are *measured* rather than asserted, and so the tests can pin them. Every shift here is the shift the C performs, including the ones that truncate toward negative infinity; ``>>`` on a negative int is an arithmetic shift in Python and on every compiler this workshop targets, so the two agree bit for bit. Two of those shifts are easy to leave out and were, until the review battery of 2026-08-01 caught them. The ring stores ``ln_q12 >> 4`` and reconstructs with ``<< 4``, which costs a floor-truncation per entry and so biases the re-estimated baseline *downward* by roughly 1-2% of a log-energy sigma; that bias exists in the firmware and is worth knowing about, because it tilts the detector toward upward alarms after every re-baseline. And ``filled`` is reset to zero after a re-baseline, so the detector refuses to re-baseline again until the ring holds only post-change blocks. Modelling the first without the second, or either without the other, reproduces the C's numbers by accident rather than by construction: see ``_vv/change-detection.md``. Why integers at all ------------------- The CUSUM itself would be perfectly happy in float. The reason to go integer is the logarithm: turning a multiplicative change in noise *power* into an additive change the CUSUM can accumulate needs a ln() per block, and ``logf`` from libm is the most expensive thing in the loop by an order of magnitude. Replacing it with an exponent lookup and a cubic costs accuracy that is measurably irrelevant here (see ``LN_MAX_ABS_ERR_NATS`` and the tests). The one hazard worth naming --------------------------- ``ln(0)`` is not representable, and a block energy of exactly zero is not hypothetical: a muted microphone, a disconnected input, or a DMA buffer read before the first conversion all produce it. ``int_ln_q12`` therefore floors its argument at 1 rather than returning a nonsense value, and the accumulator saturates rather than wrapping, because those two guards catch different halves of the same failure. This is the same pair of measures the voice-pitch capstone needed for the log of a degenerate spectral bin. """ import numpy as np # ln(m) for m in [0.5, 1), degree-3 fit, coefficients in Q12. # Fitted at Chebyshev nodes; see the tests for the accuracy pin. LN_C3, LN_C2, LN_C1, LN_C0 = 3530, -11802, 17202, -8929 LN2_Q12 = 2839 # round(ln(2) * 4096) #: Worst-case error of ``int_ln_q12`` over the range a 12-bit ADC can #: produce, in nats. Measured, not estimated; the tests re-measure it. LN_MAX_ABS_ERR_NATS = 2.4e-3 Q8_MAX = 32767 # int16 ceiling for the saturating adds def int_ln_q12(energy: int) -> int: """Natural log of a positive integer, in Q12, using integers only. Writes ``e = m * 2^p`` with ``m`` in [0.5, 1) by taking ``p`` from the bit length, then evaluates a cubic in ``m``. The C form is p = 64 - __builtin_clzll(e); m = (p >= 15) ? (e >> (p - 15)) : (e << (15 - p)); /* Q15 */ Parameters ---------- energy : int Block energy, >= 0. Zero is floored to one (see module docstring): ln(0) has no representation and a zero-energy block is a real hardware state, not a hypothetical. Returns ------- int ln(energy) in Q12 (i.e. multiplied by 4096), truncated. """ e = int(energy) if e < 0: raise ValueError("energy must be non-negative") if e == 0: e = 1 # the floor; ln -> 0, not nonsense p = e.bit_length() m_q15 = (e >> (p - 15)) if p >= 15 else (e << (15 - p)) t = LN_C3 t = ((t * m_q15) >> 15) + LN_C2 t = ((t * m_q15) >> 15) + LN_C1 t = ((t * m_q15) >> 15) + LN_C0 return p * LN2_Q12 + t class IntDriftMonitor: """Two-sided integer CUSUM on a stream of integer block energies. Parameters ---------- k_q8, h_q8 : int Reference value and threshold, in Q8 units of the in-control log-energy standard deviation. inv_sd_q8 : int ``round(256 / sd)`` with ``sd`` the in-control standard deviation of the log-energy in nats; converts a Q12 log difference into a Q8 standardized deviate. mu_q12 : int Baseline log-energy in Q12 nats, from a training run. rebase_n : int, default 128 Ring-buffer length used to re-estimate the baseline after an alarm. Must be a power of two so the division is a shift. Longer averages the estimate down; too long and the window reaches back across the change that just fired, which costs far more than it saves (measured on the embedded page). rebaseline : bool, default True Re-estimate the baseline after each alarm. Turning it off leaves a detector that alarms forever once the floor has moved, which is only useful for isolating the arithmetic from the adaptation when comparing against a float reference. Notes ----- ``push`` returns +1 for an upward alarm, -1 for downward, 0 for no alarm, and resets both arms on an alarm. """ def __init__(self, k_q8: int, h_q8: int, inv_sd_q8: int, mu_q12: int, rebase_n: int = 128, rebaseline: bool = True): if rebase_n & (rebase_n - 1): raise ValueError("rebase_n must be a power of two") if h_q8 <= 0 or k_q8 < 0: raise ValueError("h_q8 must be positive and k_q8 non-negative") self.k_q8, self.h_q8 = int(k_q8), int(h_q8) self.inv_sd_q8 = int(inv_sd_q8) self.mu_q12 = int(mu_q12) self.rebase_n = int(rebase_n) self.rebaseline = bool(rebaseline) self.rebase_shift = int(rebase_n).bit_length() - 1 self.hist = [0] * rebase_n self.head = 0 self.filled = 0 self.hi = 0 self.lo = 0 self.peak = 0 # high-water mark, for the range claim def push(self, energy: int) -> int: """Feed one block energy; return +1, -1 or 0.""" ln_q12 = int_ln_q12(energy) self.hist[self.head] = ln_q12 >> 4 # Q8 nats: fits int16 self.head = (self.head + 1) & (self.rebase_n - 1) if self.filled < self.rebase_n: self.filled += 1 z_q8 = ((ln_q12 - self.mu_q12) * self.inv_sd_q8) >> 12 self.hi = min(max(0, self.hi + z_q8 - self.k_q8), Q8_MAX) self.lo = min(max(0, self.lo - z_q8 - self.k_q8), Q8_MAX) self.peak = max(self.peak, self.hi, self.lo) if self.hi >= self.h_q8 or self.lo >= self.h_q8: direction = 1 if self.hi >= self.h_q8 else -1 self.hi = self.lo = 0 if self.rebaseline and self.filled >= self.rebase_n: s = sum(self.hist) >> self.rebase_shift self.mu_q12 = s << 4 # Q8 back to Q12 self.filled = 0 # refill before next return direction return 0 def run(self, energies) -> list: """Push a sequence; return ``(index, direction)`` for each alarm.""" out = [] for i, e in enumerate(energies): d = self.push(e) if d: out.append((i, d)) return out def adc_block_energies(n_blocks: int, sigma_lsb: float, rng, block_n: int = 64, dc_lsb: int = 2048) -> np.ndarray: """Block energies as the CFAR front end actually delivers them. This returns the CFAR's *integer energy numerator* N = M * sum(x^2) - (sum x)^2 = M * sum (x - xbar)^2, which is the quantity ``energy_push`` computes in [the CFAR firmware](../detection-theory/embedded.qmd) one line before it divides by ``M`` and casts to float. Taking it at that point is what lets the drift monitor stay in integers and keeps the promise that it adds no sampling and no buffers of its own. Two properties make this the right interface, and both are pinned by tests rather than asserted here: * **It is exactly invariant to the DC offset.** A real 12-bit ADC sits at mid-rail, and a *raw* sum of squares is then about 99.8% pedestal: a 1 dB change in the noise moves its log by 0.5% of the CUSUM's per-block slack, so a monitor fed raw sums of squares is deaf to the very drift it exists to catch. Removing the block mean cancels the offset identically in integer arithmetic, which is why ``dc_lsb`` may be set to anything without changing a single returned value. * **It costs one degree of freedom.** ``N`` is chi-square with ``M - 1``, not ``M``, so the in-control spread is ``polygamma(1, (M - 1) / 2)``. The sibling CFAR page makes the same point about its own thresholds. The constant factor ``M`` is deliberate: it keeps the value an exact integer, and in log space it is an additive constant that the stored baseline absorbs, so it never has to be divided out. Parameters ---------- n_blocks : int Number of blocks. sigma_lsb : float Noise standard deviation in ADC counts. rng : numpy.random.Generator block_n : int, default 64 Samples per block. dc_lsb : int, default 2048 Mid-rail pedestal in counts. Present so the invariance above is testable; it cannot affect the result. Returns ------- ndarray of int64 ``M * sum(x^2) - (sum x)^2`` per block. """ x = (np.rint(sigma_lsb * rng.standard_normal((n_blocks, block_n))) + int(dc_lsb)).astype(np.int64) s = x.sum(axis=1) return block_n * np.sum(x ** 2, axis=1) - s ** 2