Delay-Doppler Search¶

doppler.spectral.DelayDoppler answers one question about a stretch of samples: where is this code, and how fast is it rotating? You give it one period of a code (a PN code, a Zadoff-Chu sequence, any sequence with a sharp autocorrelation) and a frame of the received stream. It hands back a delay-Doppler map: a table with a row for every Doppler offset and a column for every delay, whose largest cell is the best guess for both.

Use it when you need the map itself: to see two emitters at once, to read a Doppler to a fraction of a bin, to rank many candidate frames, or to search a carrier offset that is larger than a code bin. If you only want bursts found in a stream and handed back, start with Acquiring a Repeated Preamble instead; it runs this same search for you.

What a map is¶

A frame is ny code periods of nx samples each, back to back, so ny * nx samples in all. Each period is correlated against the code (that is the delay axis, nx lags). The slow-time FFT then runs down every delay column across the ny periods (that is the Doppler axis).

doppler_step is 1 / (ny_out * nx) cycles per sample, so multiply by the sample rate for Hz. The map comes back flat (ny_out * nx values, row by row); reshape it to (ny_out, nx).

One map¶

A code delayed by 5 samples and rotating at 3 Doppler bins, in a frame of 16 periods:

import numpy as np
from doppler.spectral import DelayDoppler

rng = np.random.default_rng(0)
NX, NY = 64, 16                       # code length, periods per frame
DELAY, DOPPLER_BINS = 5, 3
code = rng.choice([-1.0, 1.0], NX).astype(np.complex64)

n = np.arange(NY * NX)
f = DOPPLER_BINS / (NY * NX)          # cycles per sample
x = (np.tile(np.roll(code, DELAY), NY)
     * np.exp(2j * np.pi * f * n)).astype(np.complex64)

dd = DelayDoppler(code, ny=NY)
m = dd.execute(x).reshape(NY, NX)     # a dump on every frame (dwell 1)

row, col = np.unravel_index(np.argmax(np.abs(m)), m.shape)
print(int(row) - NY // 2, int(col))   # (Doppler bins, delay)
assert (int(row) - NY // 2, int(col)) == (3, 5)

The peak sits at (3, 5): three Doppler bins from zero, five samples late. Its height is about ny * nx for a unit-amplitude signal, a little less when the carrier turns during a code period (that rotation is not corrected inside the period; the bank below is what searches it).

Three knobs shape the map:

• ny_out zero-pads the slow-time samples to ny_out > ny rows, a finer Doppler grid over the same response. A Doppler between two native bins lands on its fractional row, and the peak at a bin centre is unchanged.
• dwell coherently sums that many frames before dumping. execute returns None until the last frame of a dwell, and the map after it.
• nthreads splits the work; the map is bit-identical at any count.

# A Doppler of 3.5 bins is between two native rows. A 4x finer grid finds it.
f = 3.5 / (NY * NX)
x35 = (np.tile(np.roll(code, DELAY), NY)
       * np.exp(2j * np.pi * f * n)).astype(np.complex64)

fine = DelayDoppler(code, ny=NY, ny_out=4 * NY)
mf = fine.execute(x35).reshape(fine.ny_out, NX)
r, c = np.unravel_index(np.argmax(np.abs(mf)), mf.shape)
print((int(r) - fine.ny_out // 2) * fine.doppler_step * (NY * NX), int(c))
assert int(r) - fine.ny_out // 2 == 14 and int(c) == DELAY  # 3.5 bins

If you hold the output yourself, execute(x, out=buf) writes the map into buf (it must hold execute_max_out() values) and returns a view of it; the next dump into the same buf overwrites that view. Without out= every dump is a new array.

Batches: many frames, one call¶

A receiver rarely has one frame. It tries a frame under several hypotheses, or ranks the candidate frames from several antennas, and keeps the best. Looping over execute works but leaves the loop in Python and, with a dwell above 1, mixes the frames into one dwell. A batch is H frames back to back (H * ny * nx samples) handed over in one call:

• the object's dwell is never touched, so a batch can run in the middle of a dwell;
• map h is exactly what a fresh dwell=1 object returns for frame h, bit for bit;
• a batch that is empty, or not a whole number of frames, is a ValueError with nothing written.

What comes back depends on how you want the H maps combined:

execute_coherent and execute_power_sum differ exactly where it matters: frames x and -x cancel in the first and add in the second. execute_peaks and execute_power_max differ too, and the difference is the point: the best frame by its own peak-to-mean is not always the frame that owns the strongest cell of the pointwise maximum map, so ask for the statistic you mean.

Here one frame is tried under eight carrier hypotheses (each removes a candidate offset from a copy of the frame), and the batch picks the hypothesis that left the signal most compact:

import numpy as np

from doppler.spectral import DelayDoppler

NX, NY, DELAY = 128, 16, 37  # code length, periods per frame, code delay
OFFSET_BINS = 2.4  # true carrier offset, in code bins (1 / NX cycles/sample)
STEP_BINS, H = 0.5, 8  # hypothesis spacing and count

rng = np.random.default_rng(2683)
code = rng.choice(np.array([-1.0, 1.0]), NX).astype(np.complex64)

m = np.arange(NY * NX)
rx = np.tile(np.roll(code, DELAY), NY) * np.exp(
    2j * np.pi * OFFSET_BINS / NX * m
)
removed = STEP_BINS * np.arange(H)
batch = np.concatenate(
    [rx * np.exp(-2j * np.pi * r / NX * m) for r in removed]
).astype(np.complex64)

dd = DelayDoppler(code, ny=NY, nthreads=0)  # 0: every core
peak, row, col, mean = dd.execute_peaks(batch)
ratio = peak / mean
best = int(np.argmax(ratio))

print("hypothesis  removed (bins)  peak/mean  delay")
for h in range(H):
    print(
        f"  {h:2d}          {removed[h]:4.1f}        {ratio[h]:8.1f}   "
        f"{col[h]:3d}"
    )
print(
    f"best: hypothesis {best} ({removed[best]:.1f} bins removed; "
    f"true offset {OFFSET_BINS})"
)

power, which = dd.execute_power_max(batch)
cell = int(np.argmax(power))
print(
    f"strongest cell over all maps: delay {cell % NX}, "
    f"from hypothesis {which[cell]}"
)

execute_all takes out= as a 2-D complex64 array of shape (rows, n_out) with at least H rows; the others return new arrays.

Unknown carrier: the hypothesis bank¶

A map corrects no carrier within a code period, so a carrier that is d code bins off costs the Dirichlet factor of d: about 0.9 dB at a quarter bin, 3.9 dB at half a bin, 19 dB at 0.9 bin, and at a whole bin the correlation falls into its null (the delay is not even found). When the offset is unknown and can be larger than a bin, you search it, and a batch of corrected frame copies is the slow way to do that. The bank does the same search without the copies: every carrier shares the frame's row transforms.

The bank speaks in a carrier grid. With L fractions per bin, carrier g is an offset of g / (L * nx) cycles per sample. bank_power_max makes one map per carrier and keeps the strongest power per cell, and which carrier it was. bank_detect is the whole detector in one call:

1. Detect on a coarse grid (L_d, say half bins) and take every cell that clears t_candidate and is a local maximum.
2. Refine each candidate at its own cell on a fine grid (L_r, say quarter bins), which recovers the loss the coarse grid's straddle cost.
3. Confirm: keep the candidates whose refined power clears t_confirm.

The thresholds are powers in the map's own units, so they come from noise-only frames, as any detector's do. The example below measures the strongest coarse cell in 50 noise frames and sets the candidate and confirm thresholds a little above it:

import numpy as np

from doppler.spectral import DelayDoppler

NX, NY, DELAY = 128, 16, 37  # code length, periods per frame, code delay
CARRIER_BINS = 1.3  # the carrier, in code bins (1 / NX cycles/sample)
SNR_DB = -15.0  # per sample; the map's gain is NY * NX = 33 dB
NO_RATES = np.empty(0, np.complex64)

rng = np.random.default_rng(2683)
code = rng.choice(np.array([-1.0, 1.0]), NX).astype(np.complex64)


def noise() -> np.ndarray:
    n = rng.standard_normal(NY * NX) + 1j * rng.standard_normal(NY * NX)
    return (n / np.sqrt(2)).astype(np.complex64)


m = np.arange(NY * NX)
signal = np.tile(np.roll(code, DELAY), NY) * np.exp(
    2j * np.pi * CARRIER_BINS / NX * m
)
x = (10 ** (SNR_DB / 20) * signal + noise()).astype(np.complex64)

dd = DelayDoppler(code, ny=NY)
COARSE = (-4, 5)  # carriers g / (2 NX): -2 to +2.5 bins, on half bins

# Thresholds from noise alone: the strongest coarse cell in 50 noise frames,
# with margin, for the candidate; the confirm a little above it.
floor = max(
    float(dd.bank_power_max(noise(), 2, *COARSE, NO_RATES)[0].max())
    for _ in range(50)
)
t_candidate, t_confirm = 1.2 * floor, 1.5 * floor

power, g, rate, row, col = dd.bank_detect(
    x, 2, 4, *COARSE, t_candidate, t_confirm, NO_RATES
)
coarse_peak = float(dd.bank_power_max(x, 2, *COARSE, NO_RATES)[0].max())

print("record  delay  Doppler row  carrier (bins)  refined / coarse peak")
for k in range(len(power)):
    print(
        f"  {k}      {col[k]:3d}      {row[k]:3d}          {g[k] / 4:+5.2f}"
        f"           {power[k] / coarse_peak:.3f}"
    )

One record comes back: the code at delay 37, on the quarter-bin carrier nearest 1.3 bins, stronger than the coarse grid saw it. Records arrive in cell order, one entry per record in each of the five arrays (power, g, rate, row, col). rates (empty here) adds per-period weights, a second hypothesis axis for a rate or a window.

See also¶

• Python: Spectral API: DelayDoppler and Corr2D, every method and property
• Gallery: Delay-Doppler Search: the map, the coarse straddle loss and the refine that recovers it, drawn
• The runnable examples delay_doppler_batch_demo.py and delay_doppler_bank_demo.py, each with a C twin in native/examples/
• Acquiring a Repeated Preamble: the burst-level search that uses this map
• Design: batches and hypotheses: why the batch and the bank are shaped this way (for those who want the reasoning)