From 4a89cc871a2b6b896445221c8c118c9db17b84b7 Mon Sep 17 00:00:00 2001 From: Li-Ta Hsu Date: Mon, 31 Aug 2026 16:00:15 +0800 Subject: [PATCH] fix: two-way TOA earns its place in Chapter 4 `example_toa_positioning` argued the opposite of what the chapter teaches, by juxtaposition. Examples 1-2 solve one-way ranges at 0.1 m of measurement noise and printed 0.000 / 0.055 m; Example 5 measures RTT on an honest nanosecond timing budget and printed 0.718 m RMSE, 0.704 m of position error. A reader flipping through saw two-way as 7x worse -- a comparison between one protocol's assumptions and another's physics, not between protocols. Grepped before believing it: the word "sync" appeared nowhere in the file. Three changes. **The assumption is stated where it is made.** Examples 1, 2 and the module docstring now say that one-way TOA as solved here needs the clocks already synchronised, and that 1 ns of unmodelled offset is 0.300 m at c. Example 2's range-noise line carries its own size in nanoseconds (0.33), which is the comparison that matters: the whole error budget of Examples 1-2 is smaller than the offset a free-running consumer oscillator reaches in a fraction of a second. **Example 7 measures it.** One 20 m array, one agent position, one number -- sigma_t, the agent clock's per-timestamp jitter -- and four solves over 2000 draws each, every draw scored: case median RMSE failed predicted A one-way, clocks synchronised 0.2518 0.3018 0 0.3000 B(i) one-way, offset, position only 1.7978 1.8416 1 1.7239 B(ii) one-way, offset, clock solved 0.2545 0.3039 0 0.3014 C two-way RTT, no sync needed 0.2148 0.2613 0 0.2598 A, B(ii) and C land on DOP x sigma_range. Twenty nanoseconds of unmodelled offset costs 1.8 m; estimating it instead (Eqs. 4.24-4.26) recovers almost all of it at 0.45% of position DOP; two-way needs neither. **The bookkeeping is derived, not chosen per case.** One-way is a single receive timestamp against a scheduled epoch, sigma_range = c*sigma_t. Two-way is two timestamps on the agent's own clock plus the turnaround calibration residual, sigma_RTT = sqrt(2*st^2 + sp^2), halved by Eq. (4.7). It is the bookkeeping that flatters two-way least -- giving the transmit epoch its own sigma_t would double the gap. Verified against the simulation rather than asserted: measured range std matches `two_way_range_sigma` to within 0.5% over sigma_proc from 0 to 3 ns, so both knobs demonstrably reach the measurement. **Reported disagreement.** C comes out 13% ahead of A, where A <= C was expected. That is arithmetic, not tuning: a round trip puts twice the distance into the interval being timed while the jitter grows only as sqrt(2). The example prints the crossover -- the two protocols range equally well at sigma_proc = sqrt(2)*sigma_t = 1.414 ns -- rather than moving a parameter to recover the expected order. Checks: zero-noise sanity is 0.000000 m for all four arms. B(i)'s bias is predicted by b*||(H'H)^-1 sum(u_i)|| = 1.1323 m against 1.6976 m measured, and the example shows why -- at b = 0.01 m the ratio is 1.001, so the gap is nonlinearity and 6 m of bias on 11-17 m ranges is not a small perturbation. The agent is off-centre deliberately: at the centre of a symmetric array sum(u_i) is zero and a common bias is invisible. Also, in passing: Example 5's single-draw positioning number is labelled as one draw with the scale to expect beside it, and its Monte-Carlo RMSE gets the prediction it should match (0.8072 m against 0.7180 m over 100 draws). READMEs: the chapter's TOA section gains a clock-requirement table over the four approaches, the RTT section gains the derivation and an explicit "this does not beat TDOA", and the square dataset README names what `toa_ranges.txt` assumes. The transcript block is re-pasted from a real run and was mutation-checked -- perturbing one new digit turns it red. Output verified reproducible across processes; the committed figure is byte-identical; no dataset moved. --- ch4_rf_point_positioning/README.md | 99 +++- .../example_toa_positioning.py | 461 +++++++++++++++++- data/sim/ch4_rf_2d_square/README.md | 28 +- 3 files changed, 576 insertions(+), 12 deletions(-) diff --git a/ch4_rf_point_positioning/README.md b/ch4_rf_point_positioning/README.md index fd5b2f3..31304d0 100644 --- a/ch4_rf_point_positioning/README.md +++ b/ch4_rf_point_positioning/README.md @@ -141,6 +141,30 @@ config = json.load(open(path / "config.json")) ### TOA Positioning +**One-way TOA assumes the clocks are already synchronised**, and this is the +assumption to state before any number below it. Turning a measured time of +flight into a range means treating the transmit epoch as known, so at c every +nanosecond of unmodelled offset is 0.30 m of range error — and a free-running +consumer oscillator at a few ppm reaches tens of nanoseconds in well under a +second. A demonstration that reports 0.1 m of *measurement* noise has said +nothing about the clock, which is the larger term in practice. + +The table below sets one-way TOA beside the three answers Chapter 4 gives to +it. They are different kinds of answer — one buys infrastructure, one spends a +degree of freedom, one changes the protocol — and the clock column is the one +to read first: + +| Approach | Clock requirement | Cost | Where | +|----------|-------------------|------|-------| +| One-way TOA | transmitter and receiver **synchronised** | a synchronisation protocol in the infrastructure | Eq. (4.1)-(4.3) | +| One-way TOA + clock bias | none, but needs one extra measurement | one degree of freedom (position DOP rises) | Eq. (4.24)-(4.26) | +| **Two-way TOA (RTT)** | **none** — both timestamps are on one clock | the responder's turnaround time must be calibrated | Eq. (4.6)-(4.9) | +| TDOA | transmitters synchronised with each other, not with the receiver | no accuracy gain over TOA-with-an-estimated-clock — the same information | Eq. (4.27)-(4.42) | + +`example_toa_positioning`'s Example 7 runs the first three on one anchor array +with one timing budget, so the comparison is between protocols rather than +between one protocol's assumptions and another's physics. + ```python import numpy as np from core.rf import TOAPositioner @@ -149,6 +173,9 @@ from core.rf import TOAPositioner anchors = np.array([[0, 0], [10, 0], [10, 10], [0, 10]], dtype=float) # True position and compute ranges +# These are ranges, not pseudoranges: writing them this way is itself the +# synchronisation assumption. With an unknown receiver clock they would each +# carry a common c*dt, and the solver below could not see it. true_pos = np.array([5.0, 5.0]) ranges = np.linalg.norm(anchors - true_pos, axis=1) @@ -285,6 +312,12 @@ Example 1: TOA Positioning with Perfect Measurements ====================================================================== ... True position: [5. 5.] + +Assumption: one-way TOA as solved here needs the beacon and agent + clocks already synchronised -- that is what turns a measured time + of flight into a range. At c, 1 ns of unmodelled offset is 0.300 m + of range error. Example 3 estimates the offset; Example 7 prices + the assumption against two-way TOA, which does not need it. ... Estimated position: [5.00000012 5.00000012] Position error: 0.000000 m @@ -294,7 +327,7 @@ Iterations: 3 Example 2: TOA Positioning with Measurement Noise ====================================================================== True position: [3. 7.] -Range noise std: 0.1 m +Range noise std: 0.1 m (0.33 ns; clocks still assumed synchronised) ... Position error: 0.055 m <- ONE noise draw, not an accuracy ... @@ -303,6 +336,34 @@ Over 2000 noise draws, against Eq. (4.107): predicted HDOP x sigma : 0.1010 m measured RMS error : 0.1012 m (1.00x predicted) a single draw lands anywhere in [0.032, 0.154] m (10th-90th percentile) +... +Example 7: One-Way vs Two-Way TOA Under One Timing Budget +====================================================================== + +Anchors: 20 m square, agent at [ 8. 12.], all solves seeded at [10. 10.] + +One timing budget, three protocols: + per-timestamp jitter sigma_t : 1.000 ns + receiver clock offset b (case B) : 20.0 ns = 5.996 m + turnaround residual (case C) : 1.000 ns + +Range noise, derived from sigma_t rather than chosen per case: + one-way : sigma_range = c*sigma_t = 0.2998 m + two-way : sigma_RTT = sqrt(2*st^2 + sp^2) = 1.7321 ns + sigma_range = c*sigma_RTT/2 = 0.2596 m (0.866x one-way) +... +Zero-noise sanity (sigma_t = 0, b = 0, turnaround residual = 0): + A and B(i) 0.000000 m B(ii) 0.000000 m (bias 0.000000 m) C 0.000000 m + +Case B(i): what 20 ns of unmodelled offset does: + linearised prediction b*||(H'H)^-1 sum(u_i)|| : 1.1323 m + measured, noiseless : 1.6976 m +... + case median RMSE failed predicted + A one-way, clocks synchronised 0.2518 0.3018 0 0.3000 + B(i) one-way, offset, position only 1.7978 1.8416 1 1.7239 + B(ii) one-way, offset, clock solved 0.2545 0.3039 0 0.3014 + C two-way RTT, no sync needed 0.2148 0.2613 0 0.2598 ``` Note what the example does with that single draw: it prints it, labels it as @@ -311,6 +372,18 @@ the DOP prediction it should match. 0.1012 m measured against 0.1010 m predicted is the actual claim; the 0.055 m above it is a sample that happens to land low, and the percentile range says how little that means. +**Example 7 is the one to read against Example 2.** Examples 1 and 2 solve +one-way ranges at 0.1 m of measurement noise with the clocks assumed already +synchronised; drop that assumption and keep everything else, and the same +protocol on the same array lands 1.80 m out. Two-way TOA, on the same +per-timestamp jitter, lands at 0.21 m without needing the assumption at all. +One-way can recover by estimating the offset instead (Eqs. 4.24-4.26), which +costs one degree of freedom — 1.0053 of position DOP against 1.0007, and the +0.2545 m column above. That is the whole case for two-way TOA, and it is a +case about clocks rather than about accuracy: the sqrt(2) by which C also +out-ranges A is bookkeeping the example derives and prices, not a claim that a +round trip measures distance better. + **Visual Output:** ![TOA Positioning](figs/toa_positioning_example.png) @@ -570,12 +643,36 @@ print(f"Estimated bias: {bias_s*1e9:.2f} ns") **Two-Way TOA / RTT Model (Eqs. 4.6-4.9):** +**Why two-way exists:** both timestamps are taken on the *agent's own clock*, +so the departure and arrival epochs are read against the same time base and no +synchronisation with the beacon is required at all. That is the requirement +one-way TOA quietly assumes, and it is what two-way removes. What it costs +instead is a calibrated responder turnaround time, `Δt_proc` — which is why +Eq. (4.7) exists and why leaving it uncorrected below moves the range by 7.49 m. + The book's RTT model includes processing time and clock drift: - Eq. 4.6: Basic RTT: `d = c * (Δt_to + Δt_back) / 2` - Eq. 4.7: With processing: `d = c * (t_arrive - t_depart - Δt_proc) / 2` - Eq. 4.8: With drift: `d = c * (t_arrive - t_depart - Δt_proc - Δt_drift) / 2` - Eq. 4.9: With noise: `d̃ = ... + ω_proc + ω_drift` +**Timing-noise bookkeeping, and the factor of sqrt(2) in it.** With one +per-timestamp jitter `σ_t` describing the agent's clock, a one-way measurement +is a single receive timestamp against a scheduled transmit epoch, so +`σ_range = c·σ_t`. A round trip carries two timestamps on that one clock plus +the turnaround calibration residual `σ_proc`, so +`σ_RTT = sqrt(2·σ_t² + σ_proc²)`, and the range is *half* the round trip: +`σ_range = c·σ_RTT/2`. With `σ_proc = 0` that is `c·σ_t/sqrt(2)` — better than +one-way, because the round trip puts twice the distance into the interval being +timed while the jitter grows only as sqrt(2). The turnaround term buys it back, +and the two protocols range equally well at `σ_proc = sqrt(2)·σ_t`. Example 7 +of `example_toa_positioning` prints this crossover rather than arguing it. + +**Two-way does not thereby beat TDOA**, and nothing in this chapter claims it +does — see the DOP result further up, where TDOA and TOA-with-an-estimated-clock +carry the same information to machine precision. The axis two-way TOA wins on is +which clock somebody has to build. + ```python from core.rf import rtt_to_range, simulate_rtt_measurement, range_to_rtt diff --git a/ch4_rf_point_positioning/example_toa_positioning.py b/ch4_rf_point_positioning/example_toa_positioning.py index 8abb385..d580fae 100644 --- a/ch4_rf_point_positioning/example_toa_positioning.py +++ b/ch4_rf_point_positioning/example_toa_positioning.py @@ -6,12 +6,21 @@ W = I and Eq. (4.20); Example 6 is the one that supplies a covariance and so is genuinely weighted (Eq. 4.23). +One-way TOA carries an assumption that is easy to read past, so it is stated +wherever it is used rather than only where it is relaxed: Examples 1, 2 and 4 +take the beacon and agent clocks to be already synchronised, which is what +makes a measured time of flight a range at all. At c, 1 ns of unmodelled +offset is 0.30 m of range error. Example 3 keeps the offset as an unknown +(Eqs. 4.24-4.26) and Example 7 prices the assumption against the two-way +protocol that removes it. + Implements: - Eq. (4.1)-(4.3): TOA range measurements - Eq. (4.6)-(4.9): Two-way TOA / RTT measurement model - Eq. (4.11)-(4.13): RSS path-loss model - Eq. (4.14)-(4.23): Nonlinear TOA iterative LS, and WLS in Example 6 - Eq. (4.24)-(4.26): Joint position + clock bias estimation + - Example 7: one-way (synchronised and not) against two-way, one budget Author: Li-Ta Hsu Date: December 2025 @@ -42,6 +51,7 @@ rss_to_distance, rtt_to_range, simulate_rtt_measurement, + solve_batch, toa_range, toa_solve_with_clock_bias, ) @@ -51,7 +61,15 @@ def example_toa_perfect(): - """Example 1: TOA positioning with perfect measurements.""" + """Example 1: TOA positioning with perfect measurements. + + "Perfect" here is two assumptions, not one: the ranges carry no + measurement noise, *and* the beacon and agent clocks are already + synchronised, so a one-way time of flight converts straight to a range + (Eq. 4.1 with c*Delta_t = 0). The second is the one nothing on the + screen would otherwise mention, and it is the expensive one -- Example 7 + solves this same protocol with the clocks left alone. + """ print("=" * 70) print("Example 1: TOA Positioning with Perfect Measurements") print("=" * 70) @@ -63,6 +81,15 @@ def example_toa_perfect(): print(f"\nAnchor positions:\n{anchors}") print(f"True position: {true_pos}") + print( + "\nAssumption: one-way TOA as solved here needs the beacon and agent" + "\n clocks already synchronised -- that is what turns a measured time" + f"\n of flight into a range. At c, 1 ns of unmodelled offset is" + f" {SPEED_OF_LIGHT * 1e-9:.3f} m" + "\n of range error. Example 3 estimates the offset; Example 7 prices" + "\n the assumption against two-way TOA, which does not need it." + ) + # Compute true ranges ranges = np.array([toa_range(anchor, true_pos) for anchor in anchors]) print(f"\nTrue ranges: {ranges}") @@ -83,7 +110,15 @@ def example_toa_perfect(): def example_toa_with_noise(): - """Example 2: TOA positioning with measurement noise.""" + """Example 2: TOA positioning with measurement noise. + + Still one-way, and still synchronised: the 0.1 m here is *measurement* + noise about a range the clocks already agree on. 0.1 m is 0.33 ns, so + this example's whole error budget is smaller than the clock offset a + consumer oscillator reaches in a fraction of a second. That is the + comparison Example 7 makes, and it is why the numbers below should not be + read as what one-way TOA achieves in the field. + """ print("\n" + "=" * 70) print("Example 2: TOA Positioning with Measurement Noise") print("=" * 70) @@ -100,7 +135,11 @@ def example_toa_with_noise(): noise_std = 0.1 # 10 cm standard deviation ranges_noisy = true_ranges + np.random.randn(4) * noise_std - print(f"Range noise std: {noise_std} m") + print( + f"Range noise std: {noise_std} m" + f" ({noise_std * 1e9 / SPEED_OF_LIGHT:.2f} ns; clocks still assumed" + " synchronised)" + ) print(f"True ranges: {true_ranges}") print(f"Noisy ranges: {ranges_noisy}") @@ -363,7 +402,17 @@ def plot_toa_positioning(anchors, true_pos, estimated_pos, history=None): def example_rtt_measurement(): - """Example 5: Two-way TOA / RTT measurement model (Eqs. 4.6-4.9).""" + """Example 5: Two-way TOA / RTT measurement model (Eqs. 4.6-4.9). + + This walks the mechanism -- processing time, drift, and what each costs + if it is not corrected. Its timing budget is deliberately loose (5 ns of + processing-time uncertainty, 2 ns of drift) so the effects are visible, + which means the metre-scale numbers below are **not** comparable with + Examples 1-2: those solve one-way ranges at 0.1 m of measurement noise + and assume the clocks are already synchronised, an assumption worth + metres in its own right. Example 7 makes the comparison properly, with + one timing budget shared by both protocols. + """ print("\n" + "=" * 70) print("Example 5: Two-Way TOA / RTT Measurement Model (Eqs. 4.6-4.9)") print("=" * 70) @@ -458,9 +507,16 @@ def example_rtt_measurement(): errors.append(info["range_estimate"] - 15.0) errors = np.array(errors) + # What should this be? Only the *uncorrected* parts of the RTT survive: + # sigma_RTT = sqrt(sigma_proc^2 + sigma_drift^2), halved by Eq. (4.7). + predicted_rtt_sigma = SPEED_OF_LIGHT * np.hypot(5e-9, 2e-9) / 2.0 print(f" Mean error: {np.mean(errors):.4f} m") print(f" Std dev: {np.std(errors):.4f} m") print(f" RMSE: {np.sqrt(np.mean(errors**2)):.4f} m") + print( + f" predicted c*sqrt(5ns^2 + 2ns^2)/2: {predicted_rtt_sigma:.4f} m" + " (100 draws, so +/- ~7%)" + ) print("\n--- RTT-Based Positioning Example ---") @@ -502,8 +558,17 @@ def example_rtt_measurement(): ) error = np.linalg.norm(est_pos - true_pos[:2]) + # This run leaves clock_drift_std at 0, so only the 3 ns of processing-time + # uncertainty reaches the range: sigma_range = c * 3 ns / 2. + single_draw_sigma = SPEED_OF_LIGHT * 3e-9 / 2.0 print(f"\n Estimated position: {est_pos}") - print(f" Position error: {error:.3f} m") + print(f" Position error: {error:.3f} m <- ONE draw, not an accuracy") + print( + f" (range sigma here is c*3ns/2 = {single_draw_sigma:.3f} m, and this" + " array's HDOP is ~1," + f"\n so ~{single_draw_sigma:.2f} m is the scale to expect." + " Example 7 reports a distribution.)" + ) return @@ -636,6 +701,389 @@ def example_wls_vs_ls(): print(" draw and the LS RMSE moves by 1e-4 m.") +#: Example 7's timing budget. One piece of hardware, three protocols: the +#: agent's clock stamps an event with this much jitter whatever the protocol +#: asks it to do, so this is the single number every case is derived from. +TIMESTAMP_JITTER_S = 1.0e-9 + +#: The receiver clock offset Case B carries. An uncorrected consumer +#: oscillator at a few ppm reaches tens of nanoseconds in well under a second, +#: which is why one-way systems need a synchronisation protocol and not merely +#: a good crystal. +CLOCK_OFFSET_S = 20.0e-9 + +#: What is left of the responder's turnaround time after calibration. This is +#: the term two-way TOA pays *instead of* synchronisation, and the comparison +#: is only honest with it present -- see `two_way_range_sigma`. +TURNAROUND_RESIDUAL_S = 1.0e-9 + +#: Nominal responder turnaround, known to the agent and removed by +#: `rtt_to_range` (Eq. 4.7). Only the residual above survives the correction. +TURNAROUND_NOMINAL_S = 50.0e-9 + +#: Draws per case in Example 7. Large enough that the RMSE column is a +#: measurement rather than a draw: the standard error on an RMSE over n draws +#: is about 1/sqrt(2n), so 2000 gives ~1.6%. +TWO_WAY_TRIALS = 2000 + +#: Iteration budget for Example 7, for the reason recorded on +#: WLS_DEMO_MAX_ITERS: the library default of 10 makes `converged` a report on +#: how quickly a solve finished rather than on whether it worked. +TWO_WAY_MAX_ITERS = 50 + + +def one_way_range_sigma(jitter_s: float) -> float: + """Ranging noise of a one-way TOA measurement, from per-timestamp jitter. + + A downlink TOA is one *receive* timestamp read against a transmit epoch + the synchronised network schedules, so the whole jitter lands on the + range: + + sigma_range = c * sigma_t. + + The scheduled epoch is treated as noiseless, which is the assumption that + flatters one-way most. Giving the transmit event its own sigma_t would + make this c * sqrt(2) * sigma_t and double the gap to two-way below, so + the bookkeeping here is the conservative one. + """ + return SPEED_OF_LIGHT * jitter_s + + +def two_way_range_sigma(jitter_s: float, turnaround_residual_s: float) -> float: + """Ranging noise of a two-way (RTT) measurement, from the same jitter. + + The agent stamps departure and arrival on its **own** clock, so two + independent jitters enter the round-trip interval, and the responder's + calibrated turnaround time contributes whatever survives calibration: + + sigma_RTT = sqrt(2 * sigma_t^2 + sigma_proc^2) (Eq. 4.9) + sigma_range = c * sigma_RTT / 2 (Eq. 4.7) + + That division is worth pausing on, because it is why the measured table + does not simply order itself the way "two-way costs an extra hop" would + suggest. With sigma_proc = 0 the two-way figure is c * sigma_t / sqrt(2), + a factor sqrt(2) **better** than one-way: a round trip puts twice the + distance into the interval being timed while the jitter grows only as + sqrt(2). The turnaround residual buys it back, and the two protocols are + exactly equal at sigma_proc = sqrt(2) * sigma_t. + + So the honest statement of two-way's advantage is not that it ranges + better. It is that it needs no synchronised clock, and the sqrt(2) is a + property of the arithmetic that has to be reported rather than hidden. + """ + sigma_rtt = float(np.hypot(np.sqrt(2.0) * jitter_s, turnaround_residual_s)) + return SPEED_OF_LIGHT * sigma_rtt / 2.0 + + +class _ClockBiasSolver: + """Adapter letting :func:`core.rf.solve_batch` drive Eqs. (4.24)-(4.26). + + ``toa_solve_with_clock_bias`` is a function taking an ``(x, y, c*dt)`` + seed, not a positioner. Wrapping it keeps Case B(ii) on the same four + failure tests as the other three arms -- raised, refused, stalled, + diverged -- rather than on whichever subset a hand-rolled loop remembers. + ``example_comparison`` carries the same adapter for the same reason. + """ + + def __init__(self, anchors): + self.anchors = np.asarray(anchors, dtype=float) + + def solve(self, ranges, initial_guess, **kwargs): + state = np.concatenate([np.asarray(initial_guess, dtype=float), [0.0]]) + position, _bias, info = toa_solve_with_clock_bias( + self.anchors, ranges, state, **kwargs + ) + return position, info + + +def _simulate_rtt_ranges(anchors, true_pos, jitter_s, turnaround_residual_s, rng): + """One row of RTT-derived ranges, one per anchor (Eqs. 4.7-4.9). + + The agent's two timestamps enter through ``clock_drift_std`` and the + responder's calibration residual through ``processing_time_std``, which is + what ``simulate_rtt_measurement`` already models: both are added to the + RTT and only the *nominal* turnaround is corrected out again. + """ + ranges = np.empty(len(anchors)) + for i, anchor in enumerate(anchors): + _, info = simulate_rtt_measurement( + anchor, + true_pos, + processing_time=TURNAROUND_NOMINAL_S, + processing_time_std=turnaround_residual_s, + clock_drift_std=np.sqrt(2.0) * jitter_s, + rng=rng, + ) + ranges[i] = info["range_estimate"] + return ranges + + +def example_one_way_vs_two_way(): + """Example 7: one-way TOA (synchronised or not) against two-way TOA. + + Chapter 4 introduces two-way TOA (Eqs. 4.6-4.9) because it removes the + synchronisation requirement one-way TOA carries silently. Read in file + order, though, this file used to argue the opposite: Examples 1 and 2 + solve one-way ranges at 0.1 m of measurement noise and never mention the + clock, while Example 5 measures RTT with an honest nanosecond timing + budget and prints an error seven times larger. That is not a comparison + between protocols, it is a comparison between one protocol's assumptions + and another's physics. + + So this example puts all three on one anchor array, one agent position and + one number: ``sigma_t``, the agent clock's per-timestamp jitter. + + A one-way, clocks synchronised. The textbook assumption. + B(i) one-way, receiver clock offset by b, solved for position only. + B(ii) the same measurements with the clock as a third unknown + (Eqs. 4.24-4.26) -- correct, at the cost of one degree of + freedom. + C two-way RTT, same hardware, no synchronisation anywhere. + + **The measured ordering is C < A ~ B(ii) << B(i), which is not quite the + ordering the design predicted** -- C was expected to sit at or just behind + A. It comes out 13% ahead instead, and the reason is arithmetic rather + than a thumb on the scale: see :func:`two_way_range_sigma` for the sqrt(2) + a round trip buys and the ``sigma_proc = sqrt(2) * sigma_t`` at which it + is exactly repaid. The crossover is printed rather than argued, and the + parameters were not moved to recover the expected order. + + What the table does support, and what the chapter needs it to say, is the + B-versus-C half. Twenty nanoseconds of unmodelled offset -- less than a + second of free-running drift on a few-ppm oscillator -- puts a position- + only one-way fix 1.8 m out, on an array where the same hardware two-way + lands at 0.22 m. The offset can be estimated instead, and B(ii) recovers + almost all of it, which is the honest statement of the trade: one-way + needs *either* a synchronisation protocol *or* a spare degree of freedom; + two-way needs neither. + + Two-way does not thereby beat TDOA, and no run here claims it does. + Chapter 4's own DOP result is that TDOA and TOA-with-an-estimated-clock + carry the same information -- see the chapter README, where the two + position DOPs agree to machine precision at all 100 grid points. The + axis two-way TOA wins on is infrastructure, not accuracy. + """ + print("\n" + "=" * 70) + print("Example 7: One-Way vs Two-Way TOA Under One Timing Budget") + print("=" * 70) + + anchors = np.array([[0, 0], [20, 0], [20, 20], [0, 20]], dtype=float) + true_pos = np.array([8.0, 12.0]) + seed_pos = np.array([10.0, 10.0]) + + offset_m = CLOCK_OFFSET_S * SPEED_OF_LIGHT + sigma_one_way = one_way_range_sigma(TIMESTAMP_JITTER_S) + sigma_two_way = two_way_range_sigma(TIMESTAMP_JITTER_S, TURNAROUND_RESIDUAL_S) + sigma_rtt_s = 2.0 * sigma_two_way / SPEED_OF_LIGHT + + print( + f"\nAnchors: 20 m square, agent at {true_pos}, all solves seeded at {seed_pos}" + ) + print("\nOne timing budget, three protocols:") + print(f" per-timestamp jitter sigma_t : {TIMESTAMP_JITTER_S * 1e9:.3f} ns") + print( + f" receiver clock offset b (case B) : {CLOCK_OFFSET_S * 1e9:.1f} ns" + f" = {offset_m:.3f} m" + ) + print(f" turnaround residual (case C) : {TURNAROUND_RESIDUAL_S * 1e9:.3f} ns") + + print("\nRange noise, derived from sigma_t rather than chosen per case:") + print( + f" one-way : sigma_range = c*sigma_t = {sigma_one_way:.4f} m" + ) + print( + " two-way : sigma_RTT = sqrt(2*st^2 + sp^2) = " + f"{sigma_rtt_s * 1e9:.4f} ns" + ) + print( + " sigma_range = c*sigma_RTT/2 = " + f"{sigma_two_way:.4f} m ({sigma_two_way / sigma_one_way:.3f}x one-way)" + ) + print( + " the round trip times twice the distance while the jitter grows only" + "\n as sqrt(2), so two-way ranges sqrt(2) better before the turnaround" + "\n residual is charged; the two are equal at sigma_proc = " + f"{np.sqrt(2.0) * TIMESTAMP_JITTER_S * 1e9:.3f} ns." + ) + + # Geometry. Rows of H are the unit vectors from each anchor toward the + # agent, which is the same convention Example 2 hands to compute_dop. + offsets = true_pos - anchors + geometry = offsets / np.linalg.norm(offsets, axis=1, keepdims=True) + hdop = compute_dop(geometry)["HDOP"] + + # With the clock as a third unknown the design matrix gains a column of + # ones (Eq. 4.26), and the position DOP is the leading 2x2 block of the + # inverse -- the extra unknown can only make it larger. + augmented = np.column_stack([geometry, np.ones(len(geometry))]) + cofactor = np.linalg.inv(augmented.T @ augmented) + hdop_with_clock = float(np.sqrt(np.trace(cofactor[:2, :2]))) + + print("\nPosition DOP, and what the third unknown costs:") + print(f" position only (A, B(i), C) : {hdop:.4f}") + print( + f" position + clock (B(ii)) : {hdop_with_clock:.4f}" + f" (+{(hdop_with_clock / hdop - 1) * 100:.2f}%)" + ) + + # --- Zero-noise sanity ------------------------------------------------- + # Every case must be exact when nothing is wrong with the clock and + # nothing is wrong with the timing. A nonzero number on this line is a + # model that cannot represent its own data, never noise. + clean_ranges = np.array([toa_range(a, true_pos) for a in anchors]) + exact_a, _ = TOAPositioner(anchors, method="iterative_ls").solve( + clean_ranges, initial_guess=seed_pos, max_iters=TWO_WAY_MAX_ITERS + ) + exact_b, exact_bias, _ = toa_solve_with_clock_bias( + anchors, + clean_ranges, + np.array([*seed_pos, 0.0]), + max_iters=TWO_WAY_MAX_ITERS, + ) + # Both stds are zero, so `simulate_rtt_measurement` draws nothing; the + # Generator is passed anyway rather than None, so this line can never + # reach into the global stream if someone later gives it a nonzero std. + exact_c_ranges = _simulate_rtt_ranges( + anchors, true_pos, 0.0, 0.0, np.random.default_rng(SEED) + ) + exact_c, _ = TOAPositioner(anchors, method="iterative_ls").solve( + exact_c_ranges, initial_guess=seed_pos, max_iters=TWO_WAY_MAX_ITERS + ) + # At b = 0 case B(i) *is* case A -- same measurements, same solver -- so + # it shares the number rather than being solved twice under a second name. + print("\nZero-noise sanity (sigma_t = 0, b = 0, turnaround residual = 0):") + print( + f" A and B(i) {np.linalg.norm(exact_a - true_pos):.6f} m" + f" B(ii) {np.linalg.norm(exact_b - true_pos):.6f} m" + f" (bias {exact_bias:.6f} m)" + f" C {np.linalg.norm(exact_c - true_pos):.6f} m" + ) + + # --- What the offset does to a position-only fix ----------------------- + # Linearising about the truth, a common range bias b displaces the LS fix + # by (H'H)^-1 H' (b * 1) = b * (H'H)^-1 sum(u_i). At the centre of a + # symmetric array sum(u_i) is zero and the bias is invisible, which is why + # the agent here is off-centre. + gain = np.linalg.inv(geometry.T @ geometry) @ geometry.sum(axis=0) + predicted_bias = offset_m * float(np.linalg.norm(gain)) + biased_ranges = clean_ranges + offset_m + biased_fix, _ = TOAPositioner(anchors, method="iterative_ls").solve( + biased_ranges, initial_guess=seed_pos, max_iters=TWO_WAY_MAX_ITERS + ) + measured_bias = float(np.linalg.norm(biased_fix - true_pos)) + + # The prediction is first order in b, and b is 6 m against ranges of + # 11-17 m. Solving a deliberately tiny offset shows the prediction is + # right where it applies and the gap above is nonlinearity. + small_m = 0.01 + small_fix, _ = TOAPositioner(anchors, method="iterative_ls").solve( + clean_ranges + small_m, + initial_guess=seed_pos, + max_iters=200, + tol=1e-12, + ) + small_ratio = float(np.linalg.norm(small_fix - true_pos)) / ( + small_m * float(np.linalg.norm(gain)) + ) + + print(f"\nCase B(i): what {CLOCK_OFFSET_S * 1e9:.0f} ns of unmodelled offset does:") + print(f" linearised prediction b*||(H'H)^-1 sum(u_i)|| : {predicted_bias:.4f} m") + print(f" measured, noiseless : {measured_bias:.4f} m") + print( + f" the gap is nonlinearity, not a wrong prediction: at b = {small_m} m the" + f"\n measured displacement is {small_ratio:.3f}x predicted, and" + f" {offset_m:.2f} m of bias" + "\n on 11-17 m ranges is not a small perturbation." + ) + + # --- Monte Carlo ------------------------------------------------------- + rng = np.random.default_rng(SEED) + truth = np.tile(true_pos, (TWO_WAY_TRIALS, 1)) + + draws = rng.standard_normal((TWO_WAY_TRIALS, len(anchors))) + measurements_a = clean_ranges + draws * sigma_one_way + measurements_b = measurements_a + offset_m + measurements_c = np.array( + [ + _simulate_rtt_ranges( + anchors, true_pos, TIMESTAMP_JITTER_S, TURNAROUND_RESIDUAL_S, rng + ) + for _ in range(TWO_WAY_TRIALS) + ] + ) + + ls_solver = TOAPositioner(anchors, method="iterative_ls") + cases = [ + ( + "A one-way, clocks synchronised", + measurements_a, + ls_solver, + hdop * sigma_one_way, + ), + ("B(i) one-way, offset, position only", measurements_b, ls_solver, np.nan), + ( + "B(ii) one-way, offset, clock solved", + measurements_b, + _ClockBiasSolver(anchors), + hdop_with_clock * sigma_one_way, + ), + ( + "C two-way RTT, no sync needed", + measurements_c, + ls_solver, + hdop * sigma_two_way, + ), + ] + + print( + f"\nMonte Carlo, {TWO_WAY_TRIALS} draws per case." + " Every draw is scored: `converged`" + "\nis a step-tolerance stop, not a success flag, so failures are counted" + "\nbeside the medians rather than removed from them." + ) + print( + f" {'case':38s} {'median':>8s} {'RMSE':>8s} {'failed':>7s} {'predicted':>10s}" + ) + for label, measurements, solver, predicted in cases: + outcome = solve_batch( + solver, + measurements, + seed_pos, + truth, + max_iters=TWO_WAY_MAX_ITERS, + ) + finite = np.isfinite(outcome.errors) + rmse = float(np.sqrt(np.mean(outcome.errors[finite] ** 2))) + if np.isnan(predicted): + # B(i) is bias plus scatter, not DOP x sigma: the systematic part + # measured above and the same scatter as A, added in quadrature. + predicted = float(np.hypot(measured_bias, hdop * sigma_one_way)) + print( + f" {label:38s} {outcome.median_m:8.4f} {rmse:8.4f}" + f" {outcome.n_failed:7d} {predicted:10.4f}" + ) + + print( + "\n The predicted column is DOP x sigma_range for A, B(ii) and C, and" + "\n those three land on it. B(i) has no such prediction -- its number is" + "\n the noiseless bias and A's scatter in quadrature, and the measured" + "\n RMSE runs a few per cent above it because the response to the offset" + "\n is not linear, which is the same nonlinearity the block above" + "\n measures. A bias-dominated error is not a DOP result." + ) + + print("\n -> Two-way TOA removes the synchronisation one-way TOA assumes.") + print(" Against unsynchronised one-way it is not close: 0.2 m against") + print(" 1.8 m, on the same hardware and the same array. Against") + print(" one-way that estimates its clock it is comparable, and it gets") + print(" there without spending a degree of freedom on the offset.") + print(" It does NOT follow that two-way beats TDOA: the chapter's DOP") + print(" result is that TDOA and TOA-with-an-estimated-clock carry the") + print(" same information. What differs between the methods is which") + print(" clock somebody has to build, not how much of the geometry each") + print(" one extracts.") + + def main(): """Run all TOA positioning examples.""" # Parse arguments before doing any work, so --help answers instead of @@ -667,6 +1115,9 @@ def main(): # Example 6: WLS vs LS example_wls_vs_ls() + # Example 7: one-way (synchronised or not) against two-way, one budget + example_one_way_vs_two_way() + # Visualization print("\n" + "=" * 70) print("Generating visualization...") diff --git a/data/sim/ch4_rf_2d_square/README.md b/data/sim/ch4_rf_2d_square/README.md index 3f94765..b1eb583 100644 --- a/data/sim/ch4_rf_2d_square/README.md +++ b/data/sim/ch4_rf_2d_square/README.md @@ -608,7 +608,12 @@ single solve that "converges" to somewhere absurd makes a mean a property of that outlier. AOA used to do exactly that on 36 of 100 positions. **Key Insights**: -- TOA: 0.09m median, solves every position; requires clock sync. +- TOA: 0.09m median, solves every position; requires clock sync. That + requirement is not free and this dataset does not price it — `toa_ranges.txt` + is written as ranges, so the clocks are assumed already synchronised. Drop + the assumption and 20 ns of offset costs metres; two-way TOA (Eqs. 4.6-4.9) + removes the requirement outright, and `example_toa_positioning` Example 7 + measures both. - TDOA: 0.09m median, clock-free, and slightly *worse* than TOA on this array — its GDOP is 1.07 against TOA's 1.02, and both attain `GDOP x sigma_range`. That ordering is not a property of this array: differencing is a projection, @@ -630,22 +635,33 @@ ambiguous; bearings do not. That contrast is the reason the variant exists. This dataset directly implements RF positioning from Chapter 4: -1. **TOA (Section 4.1, Eqs. 4.1-4.3)** +1. **One-way TOA (Section 4.1, Eqs. 4.1-4.3)** - Measures propagation time → range - - Requires clock synchronization + - Requires clock synchronization — this dataset's `toa_ranges.txt` are + ranges, not pseudoranges, which *is* that assumption written down. At c, + 1 ns of unmodelled offset is 0.30 m of range error. - Position from range intersection (trilateration) + - Two ways out if the clocks are not synchronised: estimate the offset as a + third unknown (Eqs. 4.24-4.26, one degree of freedom), or use two-way TOA -2. **TDOA (Section 4.2, Eqs. 4.27-4.33)** +2. **Two-way TOA / RTT (Section 4.2.1.2, Eqs. 4.6-4.9)** + - Both timestamps are on the *agent's own* clock, so no synchronization + with the beacon is needed anywhere — this is the reason the method exists + - Pays instead for a calibrated responder turnaround time (Eq. 4.7) + - Not in this dataset; `ch4_rf_point_positioning/example_toa_positioning` + Example 7 compares it with one-way under one shared timing budget + +3. **TDOA (Section 4.2, Eqs. 4.27-4.33)** - Measures time difference → range difference - Eliminates clock bias (huge advantage!) - Position from hyperbola intersection -3. **AOA (Section 4.4, Eqs. 4.63-4.67)** +4. **AOA (Section 4.4, Eqs. 4.63-4.67)** - Measures angle of arrival - No clock required - Position from bearing intersection (triangulation) -4. **DOP (Section 4.5)** +5. **DOP (Section 4.5)** - Quantifies geometry quality - GDOP = sqrt(trace((H^T H)^{-1})) - Lower GDOP = better geometry