Python Utilities API¶

Square-clip a complex sample: clip the real and imaginary parts independently to [-lin, lin] (a square region in the IQ plane).

Parameters:

Returns:

Examples:

>>> from doppler.util import square_clip
>>> square_clip(0.5+0.25j, 1.0)   # within bounds, passed through
(0.5+0.25j)
>>> square_clip(2.0+0.5j, 1.0)    # real clipped, imag unchanged
(1+0.5j)
>>> square_clip(3.0-4.0j, 1.0)    # both components clipped
(1-1j)
>>> square_clip(0.5+0.5j, 0.25)   # smaller threshold clips both
(0.25+0.25j)
>>> square_clip(-2.0+0.0j, 1.0)   # negative real clipped
(-1+0j)

Saturate a value into [lo, hi], total over every double including NaN and both infinities. The NaN destination is a parameter because which end is safe is domain knowledge: a gain control guarding a measured power wants the ceiling, a lock statistic wants the floor. Use it at the boundary where an untrusted value first becomes persistent state -- the input of an EMA, accumulator or integrator.

fmin/fmax are not enough for this job. A plain fmin(fmax(v, lo), hi) propagates NaN on some platforms and silently returns a bound on others, and a hand-written v > hi ? hi : v leaves NaN untouched, because every comparison against NaN is false. This function has no fall-through: a value that is neither inside the interval, nor below it, nor above it can only be NaN.

Parameters:

Returns:

Examples:

>>> from doppler.util import saturate
>>> saturate(0.5, 0.0, 1.0, 1.0)     # inside the interval
0.5
>>> saturate(2.0, 0.0, 1.0, 1.0)     # above the ceiling
1.0
>>> saturate(-3.0, 0.0, 1.0, 1.0)    # below the floor
0.0
>>> saturate(float("inf"), 0.0, 1.0, 1.0)   # infinity is just above
1.0
>>> saturate(float("nan"), 0.0, 1.0, 1.0)   # NaN takes the caller's end
1.0
>>> saturate(float("nan"), 0.0, 1.0, 0.0)   # ... which may be the other
0.0

One step of a first-order exponential moving average, state + alpha*(x - state). The canonical EMA for the library: it was written out four times in two different algebraic forms before this existed, and duplicated implementations drift. The incremental form is the more accurate of the two everywhere the library operates, by a margin that grows as the average lengthens; alpha == 1 (pass-through) and alpha == 0 (frozen) are both exact. NOT total in x -- a non-finite observation poisons the state permanently, because an EMA remembers, so saturate() belongs on this function's input.

The canonical EMA for the whole library. It was written out four times before this existed — agc (power detector), async_dsss_receiver (the lock_num/lock_den pair), acc_trace (ACC_TRACE_EXP) and the recursion det_ema_alpha sizes — in two different algebraic forms, which are identical on paper and not in floating point. Duplicated implementations drift; this is the one.

Why this form, and not alpha*x + (1-alpha)*state¶

Both were measured against a 60-digit reference over 5000 steps. The incremental form written here is the more accurate one everywhere the library actually operates, by a margin that grows as the average gets longer — which is the direction a narrow-band estimator moves:

The other form wins exactly one case, and it is a boundary rather than a regime: at alpha == 1 it returns x bit-exactly while the incremental form does not (measured inexact for 9.6% of random (state, x) pairs, because state + 1*(x - state) rounds twice). That case is real — det_ema_alpha returns exactly 1.0 for "no gain requested, so no averaging" — so it is handled explicitly below rather than paid for at every alpha.

Parameters:

Returns:

Examples:

>>> from doppler.util import ema_step
>>> ema_step(0.0, 1.0, 0.5)          # halfway to the observation
0.5
>>> ema_step(2.0, 2.0, 0.25)         # at its fixed point, no motion
2.0
>>> ema_step(1.0, 7.0, 1.0)          # alpha 1 is exact pass-through
7.0
>>> ema_step(1.0, 7.0, 0.0)          # alpha 0 freezes the state
1.0

The EMA coefficient that advances d samples in one step, 1 - (1 - alpha)^d. A decimated loop updates once per chunk of d samples and must not thereby change its own time constant. Computed through expm1/log1p because the direct expression cancels catastrophically for small alpha -- 26865 ulps off at alpha 1e-5, d 1 -- and being exact at d == 1 is what lets the decimated and per-sample paths be compared bit-for-bit.

A decimated loop updates its average once per chunk of d samples and must not thereby change its own time constant. Compounding the pole exactly is what makes decim a performance knob instead of a retune.

Why expm1/log1p rather than the direct expression¶

1.0 - pow(1.0 - alpha, d) cancels catastrophically for small alpha, and the damage is worst exactly where a narrow-band estimator lives. Measured at d == 1, where the answer must be alpha itself:

dp_agc_steps used the repeated-multiply form and had this defect; it now forms BOTH its per-chunk coefficients with this function. Being exact at d == 1 is the property that lets a caller set decim = 1 and get bit-for-bit the undecimated recursion, so the decimated and per-sample paths can be compared at all.

Parameters:

Returns:

Examples:

>>> from doppler.util import ema_alpha_decim
>>> ema_alpha_decim(0.05, 1)         # d == 1 returns alpha exactly
0.05
>>> round(ema_alpha_decim(0.05, 8), 12)
0.336579568711
>>> ema_alpha_decim(1.0, 4)          # pass-through stays pass-through
1.0
>>> ema_alpha_decim(0.0, 8)          # frozen stays frozen
0.0

1 - (1 - p)^x for p in [0, 1] and real x >= 0, computed through expm1/log1p so it stays accurate when p is small. The one kernel behind two library quantities: the EMA coefficient that advances d samples (ema_alpha_decim, x = d) and the per-cell false-alarm probability that splits a search's Pfa over n independent cells (det_pfa_cell, x = 1/n).

Written directly, 1 - pow(1 - p, x) loses everything 1 - p rounded away: at p = 1e-5 it is 26865 ulps off. -expm1(x * log1p(-p)) is the same quantity with nothing cancelled.

Two library quantities are this one expression, and both call it:

• the EMA coefficient that advances d samples in one step, ema_alpha_decim(alpha, d) (x = d);
• the per-cell false-alarm probability that splits a search's pfa over n independent cells, det_pfa_cell(pfa, n) (x = 1/n, Šidák).

Parameters:

Returns:

Examples:

>>> from doppler.util import complement_power
>>> complement_power(0.05, 1.0)          # one trial is p exactly
0.05
>>> round(complement_power(0.5, 2.0), 12)  # 1 - 0.25
0.75
>>> round(complement_power(1e-3, 1 / 1000) * 1e6, 6)  # Sidak split
1.0005
>>> complement_power(0.3, 0.0)
0.0

Normalized sinc, sin(pi u)/(pi u), with sinc(0) = 1: the amplitude response of a rectangular window, so the straddle loss of a signal u bins off a DFT bin's centre.

The amplitude response of a rectangular window, which makes it the straddle loss of every correlator and DFT: a signal u bins off a bin's centre keeps sinc(u) of its amplitude in that bin.

Parameters:

Returns:

Examples:

>>> from doppler.util import sinc
>>> sinc(0.0)
1.0
>>> round(sinc(0.5), 12)                 # half a bin: 2/pi
0.636619772368
>>> abs(sinc(1.0)) < 1e-15               # the first null
True

The mean of sinc(u) over u in [0, umax]: the average amplitude loss of a signal whose offset from the nearest bin centre is uniform over umax bins. 1 for umax <= 0. 64-interval Simpson (simpson_weights) over segments of at most half a bin: within 3e-10 at any umax.

The average amplitude loss of a signal whose offset from the nearest bin centre is uniform over umax bins: the scalloping a Pd model averages over, where sinc(umax) would be only the worst case. 64-interval Simpson (dp_simpson_weights()) over segments of at most half a bin: within 3e-10 at any umax, far below any model this feeds.

Parameters:

Returns:

Examples:

>>> from doppler.util import mean_sinc
>>> mean_sinc(0.0)
1.0
>>> round(mean_sinc(0.5), 9)             # uniform over half a bin
0.8726543

Fill w with composite Simpson weights for the MEAN of a function over an interval: sum(w[i] * f(a + i*(b - a)/(n - 1))) approximates the mean of f over [a, b], for n = len(w) odd and at least 3. The weights sum to 1. Raises for any other length.

With n = w_len points, sum(w[i] * f(a + i*(b - a)/(n - 1))) is the mean of f over [a, b] (multiply by b - a for the integral). The weights are 1, 4, 2, 4, ..., 2, 4, 1 over 3 (n - 1) and sum to 1. Exact for any cubic; the error falls as (n - 1)^-4 for a smooth f.

Parameters:

Examples:

>>> import numpy as np
>>> from doppler.util import simpson_weights
>>> w = np.empty(5)
>>> simpson_weights(w)
>>> w * 12                               # 1, 4, 2, 4, 1 over 12
array([1., 4., 2., 4., 1.])
>>> u = np.linspace(0.0, 1.0, 5)
>>> round(float(w @ u**3), 12)           # mean of u^3 over [0, 1]
0.25

Fill u with the midpoint-rule nodes on [0, 1], u[k] = (k + 1/2)/n for n = len(u): the points a uniform average over n equal cells is evaluated at, each weighted 1/n.

The points a uniform average over n equal cells is evaluated at, each weighted 1/n. Scale to [a, b] as a + (b - a) * u[k].

Parameters:

Examples:

>>> import numpy as np
>>> from doppler.util import midpoint_nodes
>>> u = np.empty(4)
>>> midpoint_nodes(u)
>>> u
array([0.125, 0.375, 0.625, 0.875])

Fill z and p with the n-point Gauss-Hermite rule for a STANDARD NORMAL: sum(p[i] * f(z[i])) approximates E[f(Z)], Z ~ N(0, 1), exactly for any polynomial f of degree up to 2n - 1. For X ~ N(mu, sigma^2) evaluate f(mu + sigma * z[i]). Nodes ascend and the weights sum to 1. z and p must be the same length n >= 1; raises otherwise.

sum(p[i] * f(z[i])) approximates E[f(Z)], Z ~ N(0, 1), and is exact for any polynomial f of degree up to 2n - 1. For X ~ N(mu, sigma^2), evaluate f(mu + sigma * z[i]). The nodes ascend and are symmetric about 0; the weights sum to 1.

The nodes are the roots of the probabilists' Hermite polynomial He_n, found by Newton's method on its orthonormal recurrence h[k+1] = (z h[k] - sqrt(k) h[k-1]) / sqrt(k+1), which cannot overflow the way He_n and n! do. Each starts from the classical asymptotic guesses (Numerical Recipes' gauher). The weight of a root is 1 / (n h[n-1](z)^2).

Parameters:

Examples:

>>> import numpy as np
>>> from doppler.util import gauss_hermite
>>> z, p = np.empty(2), np.empty(2)
>>> gauss_hermite(z, p)
>>> z, p                                 # +-1, each half
(array([-1.,  1.]), array([0.5, 0.5]))
>>> z, p = np.empty(5), np.empty(5)
>>> gauss_hermite(z, p)
>>> round(float(p @ z**4), 12)           # E[Z^4] = 3
3.0

The path from here to to: where a script writes an artifact.

The start, here, defaults to the directory of the file that called fromto, so a bare file name lands beside the calling script whatever directory the run started in. The destination is joined onto the start with pathlib's /, which returns its right operand unchanged when that is absolute. So one call covers both cases a script has: no destination given (its default name, beside itself) and a destination given (used as given).

Parameters:

Returns:

Raises:

Examples:

A relative destination lands under the start, an absolute one where it says, and no destination is the start itself:

>>> from doppler.util import fromto
>>> fromto("dll_lock.png", here="/runs").as_posix()
'/runs/dll_lock.png'
>>> fromto("/tmp/dll_lock.png", here="/runs").as_posix()
'/tmp/dll_lock.png'
>>> fromto(here="/runs").as_posix()
'/runs'

With no here, the start is the directory of the file making the call:

>>> from pathlib import Path
>>> fromto("x.png") == Path(__file__).resolve().parent / "x.png"
True