The Polynomial-Phase Estimator — the reasoning¶
Scope: why PolynomialPhaseEstimator (ppe) has the shape it does — a
coherent two-dimensional search rather than a loop, one knob spanning
near-static Doppler through severe LEO chirp, and a transform four times
larger than its input. Not how to use it: that is the
API reference and the object's own certification,
src/doppler/dsss/tests/validation/ppe/results.md, which is where every
number quoted below was measured.
1. The problem: one shot, no second chance¶
A tracking loop earns its estimate over time. It can start wrong, because the error it makes on the first sample is fed back and walked out over the next thousand. That trade — accept a transient, buy accuracy — is the whole reason a loop exists, and it is unavailable to a burst receiver.
A burst is over before a loop has converged. BurstDemod gets one
frequency estimate per burst and there is no residual-walking stage behind it:
whatever the estimate is wrong by, the whole burst is demodulated wrong by. So
the estimator has to be right the first time, from the segment in hand, which
makes it feedforward and by-value — stateless, with reset() a documented
no-op, because an estimate depends only on the samples handed to that call.
That framing sets the bar. The question is not "does it converge" but "how close is a single shot, and at what SNR does it stop being close".
2. Why two dimensions¶
A signal with a constant frequency offset has linear phase. A signal seen across a satellite pass does not: the range rate itself changes, so the phase is quadratic and the frequency ramps. Fitting a constant frequency to a ramping one smears the energy across bins, and the smear is what destroys the estimate — not noise.
So the model is polynomial phase to second order: a frequency f
(cycles/sample) and a chirp rate r (cycles/sample²). The estimator searches
both:
For each chirp-rate hypothesis
r_i, dechirp the segment byexp(-j·π·r_i·m²)and take its FFT. The true(r, f)is the peak of the resulting (chirp-rate × frequency) surface.
Dechirping by the correct rate makes the phase linear again, which is exactly the condition under which an FFT concentrates the energy into one bin. Every wrong rate leaves residual curvature and a lower, wider peak. The search is therefore a matched filter over the two-parameter family, and it is fully coherent — which is what makes it optimal in the estimator sense, and what makes it hold at low SNR rather than falling off a cliff (§5).
max_rate = 0 is not a special case, it is the degenerate one. The rate
axis collapses to a single hypothesis, the search becomes a single FFT, and
the returned rate is forced to exactly 0.0 — a caller may test it for
equality. One knob therefore spans the whole range from near-static Doppler to
a severe LEO chirp, with the cost paid only where the dynamics demand it: the
work is linear in the number of rate hypotheses.
3. The grid is coarse; the refinement is what makes it accurate¶
Both axes are searched on a grid and then refined sub-grid by parabolic interpolation through the peak and its two neighbours. This is the decision that makes the object affordable, and it is worth being explicit about why, because the two axes get there differently.
The rate axis is genuinely coarse. At the certification's geometry, 53 hypotheses span ±5·10⁻⁵, a step of 1.92·10⁻⁶ — and the measured worst error across five true rates is 0.000 of a step. The refinement is doing the work; the grid only has to be fine enough to put the true rate inside the parabola's basin. The hypothesis count is forced odd, so a zero rate lands on a node rather than straddling two.
The frequency axis is refined the same way, but it is also zero-padded
4× before the transform — nfft = 4 · next_pow2(max_len). Padding buys no
information, so the reason is not resolution but conditioning: parabolic
interpolation assumes the peak is locally parabolic, and it is a better
assumption on a finely sampled main lobe than on a coarse one. The inputs here
are short — preamble partials, symbol streams — which is precisely where a
next-pow2 transform samples the main lobe with too few points for the parabola
to be a good fit.
Measured against the segment's own FFT bin (1/L — what a caller gets with no
refinement at all), the error across a full bin is 0.00e+00 bins,
noiseless, at every one of eight offsets.
The 4× is a memory decision too, and it was mis-documented.
nfftsizes three buffers, so the footprint is four times what "next pow2" suggests. The header said next-pow2 for as long as the object existed; a caller budgeting from it was out by 4×. Corrected during certification (finding F2) — and the reason it is called out here rather than quietly fixed is that the same 4× is what buys §3's accuracy. It is a trade, not an implementation detail.
4. The caller strips the modulation¶
The estimator searches for a tone under a chirp. A modulated stream has no tone, so something must remove the modulation first — and that something is the caller, deliberately, because only the caller knows which case it is in.
| case | what the caller does | cost |
|---|---|---|
| data-aided | multiply by the conjugate of the known symbols | none — full SNR retained |
| non-data-aided | raise an M-PSK stream to the M-th power (BPSK: square) | squaring loss, and the estimate comes back scaled |
The M-th-power trick is the one with a trap in it. Raising to the M-th power
multiplies the phase by M, which multiplies both f and r by M — so a
BPSK caller must halve both results. Putting that halving inside the estimator
would be wrong: the estimator cannot know whether it was handed a squared
stream or a clean tone, and a function that silently halves a correct answer is
worse than one that documents the convention.
5. The envelope¶
The header's claim was "matched-filter optimal, holds at low SNR". Nothing measured it until certification, and the envelope is what a caller sizing a burst preamble actually needs — how long a preamble must be for the estimate to be good enough is a question about this table, not about the algorithm.
| input SNR (dB) | median error (bins) | worst of 8 (bins) |
|---|---|---|
| +30 | 0.0005 | 0.0008 |
| +20 | 0.0004 | 0.0062 |
| +10 | 0.0041 | 0.0183 |
| 0 | 0.0192 | 0.0502 |
| −10 | 0.0486 | 0.1382 |
It degrades rather than breaks: at 0 dB input SNR the estimate is still inside 0.05 of a bin over eight noise draws, and at −10 dB it is inside 0.14. The estimate stays on the right peak and loses precision, which is the signature of a coherent search — the failure mode of an incoherent one is to pick a different peak entirely, and there is no graceful version of that.
snr_db is post-integration, and that matters more than it looks. The
field is a peak-to-mean taken after the coherent transform, so it carries the
processing gain: on identical input it grows with segment length (23.9 dB at
L=256, 29.9 dB at L=1024 — quadrupling adds 6.03 dB against the ideal 6.02).
A caller thresholding on it as though it were an input-referred SNR is
comparing an integrated number against an un-integrated one, and will read the
same signal as stronger simply for having handed over a longer segment.
6. Where it sits¶
ppe is the feedforward branch of the burst chain. Its counterpart is
BurstDespreader, which closes carrier and code
loops across the burst — better, when the burst is long enough to converge, and
unavailable when it is not. BurstDemod takes the feedforward branch for
exactly that reason, and drives ppe once per burst.
The composition it serves, and the search → refine → demod shape that
surrounds it, are on the
DsssBurstReceiver page. The naming — why
PolynomialPhaseEstimator rather than the shorter ChirpEstimator — is
§4.2 of the API Taxonomy.
7. The C surface¶
typedef struct
{
double freq_norm; /* frequency, cycles/sample, in [-0.5, 0.5) */
double rate_norm; /* chirp rate, cycles/sample^2 */
double snr_db; /* winning-row peak-to-mean (rough estimate) */
} ppe_result_t;
ppe_state_t *ppe_create (size_t max_len, double max_rate);
void ppe_destroy (ppe_state_t *state);
void ppe_reset (ppe_state_t *state); /* documented no-op */
ppe_result_t ppe_estimate (ppe_state_t *state, const float _Complex *x,
size_t n);
n must lie in [4, max_len]; out of range zeroes every field rather than
returning an error, so a caller that ignores the bound reads zeros rather than
garbage. The object is stateless and by-value — the measure-suite pattern —
and composes fft_core plus the spectral_core window and peak-finding free
functions rather than carrying its own.