test(pps): analyze loop stability and telemetry
This commit is contained in:
@@ -70,6 +70,48 @@ python3 -m pps_end_to_end_validation.pps_end_to_end \
|
||||
--trusted-output-delay-ns 35
|
||||
```
|
||||
|
||||
## Closed-loop stability and dynamic tuning
|
||||
|
||||
The bit-level command above isolates the absolute epoch convention. The
|
||||
companion stability analysis connects that result to the fixed-point
|
||||
firmware/RTL PI model, replays the real GDOP and tracked-satellite history in
|
||||
`fusion/logs/gps_live_5min_20260713.csv`, and synthesizes only the telemetry
|
||||
that the repository did not capture:
|
||||
|
||||
```bash
|
||||
python3 -m pps_end_to_end_validation.stability_analysis
|
||||
python3 -m unittest \
|
||||
pps_end_to_end_validation/test_stability_analysis.py
|
||||
```
|
||||
|
||||
Generated results are in `output/stability/`:
|
||||
|
||||
- `stability_report.md`: real-data inventory, discrete poles, absolute phase
|
||||
wander, period jitter, dynamic parameter sweep, and recommendations;
|
||||
- `parameter_sweep.csv`: machine-readable results for every scenario/tuning;
|
||||
- `synthetic_telemetry.csv`: 5 Hz ABS_SAMPLE, REM_CODE_PHASE, corrected 125 MHz
|
||||
capture, HDOP/VDOP/TDOP/GDOP, C/N0, pseudorange noise, clock sigma/drift,
|
||||
oscillator and timing-error inputs.
|
||||
|
||||
The real capture contains 989 GDOP values and 990 tracked-satellite values
|
||||
(eight throughout), but no PPS phase/rate, ABS_SAMPLE, C/N0, TDOP/HDOP,
|
||||
pseudorange residual, oscillator temperature, or trusted-edge/TIC samples.
|
||||
The generated CSV labels its replacements as synthetic; it must not be read as
|
||||
board telemetry.
|
||||
|
||||
The current acquire, precision, and robust fixed-point poles are inside the
|
||||
unit circle at the 200 ms update rate, and the stress analysis observes no
|
||||
unlock. It also exposes an important metrics issue: cycle-to-cycle PPS period
|
||||
jitter can remain below 20 ns while absolute phase moves slowly through tens
|
||||
of nanoseconds. A physical 20 ns maximum therefore requires a common-clock
|
||||
time-interval-counter capture and cable/output-delay calibration; a p95 period
|
||||
jitter assertion or synthetic run cannot certify it.
|
||||
|
||||
The sweep includes the current coefficients, quieter/faster precision loops,
|
||||
a robust-noise guard, a larger Doppler-drift gain, and a quality-scheduled
|
||||
candidate. The generated report contains the current numerical comparison.
|
||||
No production constants are changed by the analysis.
|
||||
|
||||
## 20 ms diagnosis
|
||||
|
||||
The current tracking wrapper shifts the current bit into `shift_reg_n`, while
|
||||
|
||||
@@ -0,0 +1,25 @@
|
||||
scenario,tuning,runs,lock_p95_s,absolute_rms_median_ns,absolute_max_p95_ns,absolute_peak_to_peak_median_ns,period_rms_median_ns,period_p95_median_ns,period_max_worst_ns,wander_10s_peak_to_peak_median_ns,profile_transitions_total,unlocks_total
|
||||
real_geometry_board_noise,current,24,10.770000000000001,10.881176406988356,32.0,48.00000000000001,3.867084378182096,8.0,8.000000000000002,36.8,0,0
|
||||
real_geometry_board_noise,robust_noise_guard,24,10.770000000000001,10.881176406988356,32.0,48.00000000000001,3.867084378182096,8.0,8.000000000000002,36.8,0,0
|
||||
real_geometry_board_noise,quiet_precision,24,10.770000000000001,6.311779001549132,16.000000000000004,24.0,2.7197398634103065,7.999999999999999,8.000000000000005,20.0,0,0
|
||||
real_geometry_board_noise,fast_precision,24,10.770000000000001,11.482151620012546,40.0,56.0,4.4110037618184545,8.0,16.0,40.0,0,0
|
||||
real_geometry_board_noise,drift_gain_1_2,24,10.770000000000001,9.830551108590587,32.0,40.0,3.761748764716249,8.0,8.000000000000002,33.6,0,0
|
||||
real_geometry_board_noise,quality_scheduled_drift,24,10.770000000000001,9.830551108590587,32.0,40.0,3.761748764716249,8.0,8.000000000000002,33.6,0,0
|
||||
scope_like_multipath,current,24,10.770000000000001,22.943340809908015,62.79999999999998,88.0,4.727734084565672,8.000000000000002,16.0,80.80000000000001,0,0
|
||||
scope_like_multipath,robust_noise_guard,24,10.770000000000001,22.943340809908015,62.79999999999998,88.0,4.727734084565672,8.000000000000002,16.0,80.80000000000001,0,0
|
||||
scope_like_multipath,quiet_precision,24,10.770000000000001,17.147532219212493,40.0,63.99999999999999,3.7617487647162484,8.000000000000004,8.000000000000005,58.4,0,0
|
||||
scope_like_multipath,fast_precision,24,10.770000000000001,22.40714058089752,56.00000000000001,92.0,5.040539175116031,8.000000000000002,16.0,78.0,0,0
|
||||
scope_like_multipath,drift_gain_1_2,24,10.8,21.354092092763203,62.79999999999998,88.0,4.465388951979694,8.0,16.0,74.80000000000001,0,0
|
||||
scope_like_multipath,quality_scheduled_drift,24,10.8,21.354092092763203,62.79999999999998,88.0,4.465388951979694,8.0,16.0,74.80000000000001,0,0
|
||||
quality_transitions,current,24,10.770000000000001,18.351896116449183,96.0,80.0,4.4110037618184545,8.0,23.99999999999999,73.20000000000002,72,0
|
||||
quality_transitions,robust_noise_guard,24,10.770000000000001,17.481868755279223,80.0,80.0,4.374396022586886,8.0,16.00000000000001,67.6,72,0
|
||||
quality_transitions,quiet_precision,24,10.770000000000001,10.467771592924619,56.0,47.99999999999999,3.0800577533580347,7.999999999999999,16.0,42.0,72,0
|
||||
quality_transitions,fast_precision,24,10.770000000000001,19.145424296464835,88.00000000000001,88.0,4.861884148700504,8.0,24.000000000000004,73.6,72,0
|
||||
quality_transitions,drift_gain_1_2,24,10.770000000000001,16.291766736261422,83.39999999999999,76.0,4.089452052592408,8.0,16.000000000000004,63.599999999999994,72,0
|
||||
quality_transitions,quality_scheduled_drift,24,10.770000000000001,16.766578539927025,83.39999999999999,80.0,4.224793965697071,8.0,16.00000000000001,67.6,72,0
|
||||
thermal_0p05ppm,current,24,10.770000000000001,39.2873108676108,102.79999999999998,136.0,5.727479384554999,8.000000000000005,24.0,124.0,0,0
|
||||
thermal_0p05ppm,robust_noise_guard,24,10.770000000000001,39.2873108676108,102.79999999999998,136.0,5.727479384554999,8.000000000000005,24.0,124.0,0,0
|
||||
thermal_0p05ppm,quiet_precision,24,10.770000000000001,210.25812707207118,416.0,648.0,13.000191801484194,24.000000000000004,32.0,632.4,0,0
|
||||
thermal_0p05ppm,fast_precision,24,10.770000000000001,21.850277755888733,63.39999999999999,96.0,5.82476895848467,8.000000000000002,16.000000000000004,77.60000000000001,0,0
|
||||
thermal_0p05ppm,drift_gain_1_2,24,10.770000000000001,38.40828815100511,96.0,136.0,5.656836393570529,8.000000000000005,16.000000000000004,122.0,0,0
|
||||
thermal_0p05ppm,quality_scheduled_drift,24,10.770000000000001,38.40828815100511,96.0,136.0,5.656836393570529,8.000000000000005,16.000000000000004,122.0,0,0
|
||||
|
@@ -0,0 +1,70 @@
|
||||
# PPS stability analysis: current RTL + 180 ms firmware model
|
||||
|
||||
This report is generated by `stability_analysis.py`. The absolute LNAV epoch path and the stochastic control loop are evaluated separately so a 20 ms semantic error cannot be hidden by loop noise.
|
||||
|
||||
## Absolute epoch result
|
||||
|
||||
- Current RTL delivered the validated subframe after 9 following bits.
|
||||
- With 180 ms postponement, reference error was -0.000032 ms.
|
||||
- DUT minus trusted PPS was 32.01 ns. That remainder is the model's configured four-PCLK capture pipeline, not a measured board delay.
|
||||
- The production source still defines `POSTPONE_TOW` as 160 ms; this analysis uses the requested 180 ms candidate configuration.
|
||||
|
||||
## Repository hardware data
|
||||
|
||||
- Real GPS capture: 989 GDOP samples and 990 tracked-satellite samples over about five minutes.
|
||||
- GDOP min/median/p95/max: 1.86 / 1.93 / 3.15 / 3.15.
|
||||
- Tracked satellites: 8..8 (median 8).
|
||||
- Horizontal scatter RMS/p95/max: 8.58 / 15.15 / 23.77 m. No surveyed truth was recorded.
|
||||
- No real ABS_SAMPLE, REM_CODE_PHASE, TDOP/HDOP, C/N0, pseudorange residual, PPS phase/rate/integrator, temperature, or trusted-edge TIC samples exist in the repository. Those columns in the companion CSV are explicitly synthetic.
|
||||
|
||||
Synthetic timing uses the board-calibration values already encoded in the PPS tests: TDOP 1.50, GDOP 2.92, clock sigma 15.87 m (52.94 ns), C/N0 40.3/34.8 dB-Hz, drift noise 3.47 m/s, seven-satellite reference point, -291.75 m/s datum range rate, and 1.592 ppm local oscillator. Real GDOP and satellite count replace the geometry fields when available.
|
||||
|
||||
## Current-loop results
|
||||
|
||||
`period max` is the worst |edge[n]-edge[n-1]| across all seeds. `absolute p-p` and `10 s wander p-p` expose the slow left/right motion that period RMS hides.
|
||||
|
||||
| scenario | abs RMS median (ns) | abs p-p median (ns) | 10 s wander p-p (ns) | period RMS (ns) | period max (ns) | unlocks |
|
||||
|---|---:|---:|---:|---:|---:|---:|
|
||||
| real_geometry_board_noise | 10.88 | 48.00 | 36.80 | 3.87 | 8.00 | 0 |
|
||||
| scope_like_multipath | 22.94 | 88.00 | 80.80 | 4.73 | 16.00 | 0 |
|
||||
| quality_transitions | 18.35 | 80.00 | 73.20 | 4.41 | 24.00 | 0 |
|
||||
| thermal_0p05ppm | 39.29 | 136.00 | 124.00 | 5.73 | 24.00 | 0 |
|
||||
|
||||
The loop does not numerically oscillate or lose lock in these runs. The scope-like symptom is low-frequency phase wander: it can be tens of nanoseconds while one-second period jitter remains below 20 ns. Therefore a `p95 period jitter < 20 ns` assertion is not sufficient for an absolute ±20 ns requirement.
|
||||
|
||||
## Linear stability
|
||||
|
||||
The fixed-point coefficients give the following discrete poles at the 200 ms navigation update rate:
|
||||
|
||||
| profile | natural frequency | damping | poles | max magnitude |
|
||||
|---|---:|---:|---|---:|
|
||||
| acquire | 0.120 Hz | 0.90 | 0.852913+0.033242j, 0.852913-0.033242j | 0.853561 |
|
||||
| precision | 0.035 Hz | 0.90 | 0.959449+0.017031j, 0.959449-0.017031j | 0.959600 |
|
||||
| robust | 0.012 Hz | 1.00 | 0.986668+0.000000j, 0.982945+0.000000j | 0.986668 |
|
||||
|
||||
All magnitudes are below one. Profile changes preserve the PI integrator, and the RTL snapshots Kp/Ki/max-slew/lock limits together with the next commit, so the current acquire/precision/robust switching is structurally sound. Stable individual modes do not prove arbitrary-switching stability, but the asymmetric 8-update upgrade / 3-update downgrade hysteresis and the stress runs show no coefficient chatter or unlock.
|
||||
|
||||
## Dynamic parameter sweep
|
||||
|
||||
| tuning | worst period max (ns) | median abs RMS over scenarios (ns) | transitions | unlocks |
|
||||
|---|---:|---:|---:|---:|
|
||||
| current | 24.00 | 20.65 | 72 | 0 |
|
||||
| robust_noise_guard | 24.00 | 20.21 | 72 | 0 |
|
||||
| quiet_precision | 32.00 | 13.81 | 72 | 0 |
|
||||
| fast_precision | 24.00 | 20.50 | 72 | 0 |
|
||||
| drift_gain_1_2 | 16.00 | 18.82 | 72 | 0 |
|
||||
| quality_scheduled_drift | 16.00 | 19.06 | 72 | 0 |
|
||||
|
||||
The safest small candidate is `robust_noise_guard`: retain the current 0.035 Hz precision loop and 0.012 Hz robust loop, but reduce only the robust-profile WLS residual gain from 1/32 to 1/64. In the synthetic quality-transition case its worst period step changed from 24.00 to 16.00 ns. This is a candidate for hardware A/B testing, not enough evidence for a production constant change.
|
||||
|
||||
The best synthetic all-scenario period result is `quality_scheduled_drift`: use drift gain 1/2 only with a precision-quality Doppler solution, fall back to the current 1/4 gain in robust quality, and combine that with the 1/64 robust bias gain. Its worst period step was 16.00 ns with no unlocks. It is the most meaningful dynamic A/B candidate, but the repository has no real drift-covariance or TIC samples with which to validate the threshold.
|
||||
|
||||
Narrowing precision to 0.012 Hz improves stationary noise but lags even the explicit 0.05 ppm/120 s synthetic oscillator variation. Widening it improves oscillator tracking but passes more WLS/multipath motion. A geometry-only bandwidth scheduler cannot resolve that tradeoff.
|
||||
|
||||
## Recommended architecture update
|
||||
|
||||
1. Keep current acquire/precision/robust PI profiles until real TIC and oscillator-temperature data are available. Do not tune from scope appearance alone.
|
||||
2. Log the existing firmware status line at every solve during validation: ABS_SAMPLE, REM_CODE_PHASE, capture sequence/counter, TDOP/HDOP/GDOP, C/N0, pseudorange residual RMS, raw/filtered bias, innovation, bias drift, hardware phase, rate, integrator, profile, lock, FPGA/oscillator temperature, DUT edge and trusted edge.
|
||||
3. A/B test the simulated scheduled filter first: precision drift gain 1/2, robust drift gain 1/4, precision bias gain 1/32, and robust bias gain 1/64. Longer term, use estimated phase/drift covariance to interpolate those software gains within bounds. Keep RTL PI bandwidth changes hysteretic and preserve its integrator as today.
|
||||
4. If temperature/frequency acceleration is the main source, add oscillator temperature feed-forward or a phase/frequency/frequency-drift estimator. Switching PI bandwidth alone either passes WLS wander or lags the oscillator.
|
||||
5. Define acceptance separately: calibrated mean offset, absolute centered error/p-p wander, one-second period error, holdover, and re-lock. A physical max-20 ns claim requires a common-clock TIC run after cable/output-delay calibration; this synthetic test cannot certify it.
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,865 @@
|
||||
#!/usr/bin/env python3
|
||||
"""End-to-end PPS epoch, telemetry, and closed-loop stability analysis.
|
||||
|
||||
The bit-level model proves the current-RTL/180 ms epoch convention. This
|
||||
module then drives the fixed-point firmware/RTL loop model with the telemetry
|
||||
that is present in the repository and deterministic synthetic replacements for
|
||||
the timing fields that were not captured on the board.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import csv
|
||||
from dataclasses import dataclass, replace
|
||||
import math
|
||||
from pathlib import Path
|
||||
import random
|
||||
import re
|
||||
import statistics
|
||||
from typing import Iterable
|
||||
|
||||
from pps_python.adaptive_pps import (
|
||||
ACQUIRE_PROFILE,
|
||||
PROFILES,
|
||||
FixedPointAdaptivePps,
|
||||
LoopProfile,
|
||||
Observation,
|
||||
Profile,
|
||||
QualityMetrics,
|
||||
TRACKING_HZ,
|
||||
classify_quality,
|
||||
)
|
||||
|
||||
from .pps_end_to_end import DecoderMode, EndToEndConfig, run_end_to_end
|
||||
|
||||
|
||||
RF_HZ = 16_367_981.0
|
||||
CODE_FREQ_HZ = 1_023_000.0
|
||||
SPEED_OF_LIGHT_MPS = 299_792_458.0
|
||||
UPDATE_PERIOD_S = 0.2
|
||||
BOARD_CLOCK_SIGMA_M = 15.87
|
||||
BOARD_TDOP = 1.50
|
||||
BOARD_GDOP = 2.92
|
||||
BOARD_CN0_MEAN_DBHZ = 40.3
|
||||
BOARD_CN0_MIN_DBHZ = 34.8
|
||||
BOARD_DRIFT_SIGMA_MPS = 3.47
|
||||
BOARD_OSCILLATOR_PPM = 1.592
|
||||
BOARD_REFERENCE_RANGE_RATE_MPS = -291.75
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class RealTelemetry:
|
||||
gdop_time_s: tuple[float, ...]
|
||||
gdop: tuple[float, ...]
|
||||
tracked_time_s: tuple[float, ...]
|
||||
tracked: tuple[int, ...]
|
||||
horizontal_rms_m: float
|
||||
horizontal_p95_m: float
|
||||
horizontal_max_m: float
|
||||
|
||||
def gdop_at(self, time_s: float) -> float:
|
||||
return _linear_sample(self.gdop_time_s, self.gdop, time_s)
|
||||
|
||||
def tracked_at(self, time_s: float) -> int:
|
||||
if not self.tracked:
|
||||
return 8
|
||||
index = _previous_index(self.tracked_time_s, time_s)
|
||||
return self.tracked[index]
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Scenario:
|
||||
name: str
|
||||
multipath: bool = False
|
||||
geometry_transitions: bool = False
|
||||
thermal_amplitude_ppm: float = 0.0
|
||||
thermal_period_s: float = 120.0
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Tuning:
|
||||
name: str
|
||||
profiles: dict[Profile, LoopProfile]
|
||||
acquire: LoopProfile
|
||||
drift_gain: float
|
||||
robust_drift_gain: float | None = None
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class SyntheticEpoch:
|
||||
scenario: str
|
||||
time_s: float
|
||||
abs_sample: int
|
||||
rem_code_phase: int
|
||||
corrected_abs_sample: float
|
||||
capture_counter_125: int
|
||||
capture_correction_ns: float
|
||||
used_satellites: int
|
||||
gdop: float
|
||||
hdop: float
|
||||
vdop: float
|
||||
tdop: float
|
||||
cn0_mean_dbhz: float
|
||||
cn0_min_dbhz: float
|
||||
pseudorange_noise_rms_m: float
|
||||
clock_bias_sigma_m: float
|
||||
clock_drift_error_mps: float
|
||||
timing_error_ns: float
|
||||
oscillator_ppm: float
|
||||
quality_profile: str
|
||||
observation: Observation
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class RunMetrics:
|
||||
first_lock_s: float
|
||||
absolute_rms_ns: float
|
||||
absolute_p95_ns: float
|
||||
absolute_max_ns: float
|
||||
absolute_peak_to_peak_ns: float
|
||||
period_rms_ns: float
|
||||
period_p95_ns: float
|
||||
period_max_ns: float
|
||||
wander_10s_peak_to_peak_ns: float
|
||||
profile_transitions: int
|
||||
unlocks: int
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Aggregate:
|
||||
scenario: str
|
||||
tuning: str
|
||||
runs: int
|
||||
lock_p95_s: float
|
||||
absolute_rms_median_ns: float
|
||||
absolute_max_p95_ns: float
|
||||
absolute_peak_to_peak_median_ns: float
|
||||
period_rms_median_ns: float
|
||||
period_p95_median_ns: float
|
||||
period_max_worst_ns: float
|
||||
wander_10s_peak_to_peak_median_ns: float
|
||||
profile_transitions_total: int
|
||||
unlocks_total: int
|
||||
|
||||
|
||||
def _percentile(values: Iterable[float], fraction: float) -> float:
|
||||
ordered = sorted(values)
|
||||
if not ordered:
|
||||
return math.nan
|
||||
position = (len(ordered) - 1) * fraction
|
||||
lower = math.floor(position)
|
||||
upper = math.ceil(position)
|
||||
if lower == upper:
|
||||
return ordered[lower]
|
||||
weight = position - lower
|
||||
return ordered[lower] * (1.0 - weight) + ordered[upper] * weight
|
||||
|
||||
|
||||
def _previous_index(times: tuple[float, ...], time_s: float) -> int:
|
||||
if not times:
|
||||
return 0
|
||||
low = 0
|
||||
high = len(times)
|
||||
while low < high:
|
||||
middle = (low + high) // 2
|
||||
if times[middle] <= time_s:
|
||||
low = middle + 1
|
||||
else:
|
||||
high = middle
|
||||
return max(0, low - 1)
|
||||
|
||||
|
||||
def _linear_sample(
|
||||
times: tuple[float, ...], values: tuple[float, ...], time_s: float
|
||||
) -> float:
|
||||
if not values:
|
||||
return 1.93
|
||||
if time_s <= times[0]:
|
||||
return values[0]
|
||||
if time_s >= times[-1]:
|
||||
return values[-1]
|
||||
lower = _previous_index(times, time_s)
|
||||
upper = min(lower + 1, len(times) - 1)
|
||||
span = times[upper] - times[lower]
|
||||
if span <= 0.0:
|
||||
return values[lower]
|
||||
fraction = (time_s - times[lower]) / span
|
||||
return values[lower] + fraction * (values[upper] - values[lower])
|
||||
|
||||
|
||||
def load_real_telemetry(software_dir: Path) -> RealTelemetry:
|
||||
"""Load the only real GPS capture and its horizontal replay output."""
|
||||
log_path = software_dir / "fusion/logs/gps_live_5min_20260713.csv"
|
||||
comparison_path = (
|
||||
software_dir
|
||||
/ "fusion/output/model_compare_5min_20260713/model_comparison.csv"
|
||||
)
|
||||
gdop_time: list[float] = []
|
||||
gdop: list[float] = []
|
||||
tracked_time: list[float] = []
|
||||
tracked: list[int] = []
|
||||
if log_path.exists():
|
||||
with log_path.open(
|
||||
"r", encoding="utf-8", errors="replace", newline=""
|
||||
) as source:
|
||||
for row in csv.reader(source):
|
||||
if len(row) < 3:
|
||||
continue
|
||||
try:
|
||||
time_s = float(row[0])
|
||||
except ValueError:
|
||||
continue
|
||||
match = re.search(r"GDOP=([0-9.]+)", row[2])
|
||||
if match:
|
||||
gdop_time.append(time_s)
|
||||
gdop.append(float(match.group(1)))
|
||||
match = re.search(r"tracked satellites=(\d+)", row[2])
|
||||
if match:
|
||||
tracked_time.append(time_s)
|
||||
tracked.append(int(match.group(1)))
|
||||
|
||||
radii: list[float] = []
|
||||
if comparison_path.exists():
|
||||
with comparison_path.open("r", encoding="utf-8", newline="") as source:
|
||||
for row in csv.DictReader(source):
|
||||
try:
|
||||
radii.append(
|
||||
math.hypot(float(row["gps_e_m"]), float(row["gps_n_m"]))
|
||||
)
|
||||
except (KeyError, TypeError, ValueError):
|
||||
continue
|
||||
horizontal_rms = (
|
||||
math.sqrt(statistics.mean(value * value for value in radii))
|
||||
if radii
|
||||
else math.nan
|
||||
)
|
||||
return RealTelemetry(
|
||||
tuple(gdop_time),
|
||||
tuple(gdop),
|
||||
tuple(tracked_time),
|
||||
tuple(tracked),
|
||||
horizontal_rms,
|
||||
_percentile(radii, 0.95),
|
||||
max(radii, default=math.nan),
|
||||
)
|
||||
|
||||
|
||||
def scenario_matrix() -> tuple[Scenario, ...]:
|
||||
return (
|
||||
Scenario("real_geometry_board_noise"),
|
||||
Scenario("scope_like_multipath", multipath=True),
|
||||
Scenario("quality_transitions", geometry_transitions=True),
|
||||
Scenario("thermal_0p05ppm", thermal_amplitude_ppm=0.05),
|
||||
)
|
||||
|
||||
|
||||
def tuning_matrix() -> tuple[Tuning, ...]:
|
||||
current = dict(PROFILES)
|
||||
guarded = dict(PROFILES)
|
||||
guarded[Profile.ROBUST] = replace(
|
||||
guarded[Profile.ROBUST], software_bias_gain=1.0 / 64.0
|
||||
)
|
||||
quiet = dict(guarded)
|
||||
quiet[Profile.PRECISION] = replace(
|
||||
quiet[Profile.PRECISION],
|
||||
natural_frequency_hz=0.012,
|
||||
damping=1.0,
|
||||
software_bias_gain=1.0 / 16.0,
|
||||
)
|
||||
fast = dict(guarded)
|
||||
fast[Profile.PRECISION] = replace(
|
||||
fast[Profile.PRECISION], natural_frequency_hz=0.050, damping=1.0
|
||||
)
|
||||
return (
|
||||
Tuning("current", current, ACQUIRE_PROFILE, 1.0 / 4.0),
|
||||
Tuning("robust_noise_guard", guarded, ACQUIRE_PROFILE, 1.0 / 4.0),
|
||||
Tuning("quiet_precision", quiet, ACQUIRE_PROFILE, 1.0 / 4.0),
|
||||
Tuning("fast_precision", fast, ACQUIRE_PROFILE, 1.0 / 4.0),
|
||||
Tuning("drift_gain_1_2", guarded, ACQUIRE_PROFILE, 1.0 / 2.0),
|
||||
Tuning(
|
||||
"quality_scheduled_drift",
|
||||
guarded,
|
||||
ACQUIRE_PROFILE,
|
||||
1.0 / 2.0,
|
||||
1.0 / 4.0,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
def _weak_geometry(time_s: float) -> bool:
|
||||
return (80.0 <= time_s < 110.0) or (180.0 <= time_s < 205.0)
|
||||
|
||||
|
||||
def synthesize_epochs(
|
||||
telemetry: RealTelemetry,
|
||||
scenario: Scenario,
|
||||
*,
|
||||
duration_s: float,
|
||||
seed: int,
|
||||
) -> list[SyntheticEpoch]:
|
||||
"""Combine real GDOP/satellite count with documented board noise."""
|
||||
rng = random.Random(seed)
|
||||
result: list[SyntheticEpoch] = []
|
||||
count = round(duration_s / UPDATE_PERIOD_S)
|
||||
rf_initial = (1 << 32) - round(10.25 * RF_HZ) - 1000
|
||||
for update in range(count):
|
||||
time_s = update * UPDATE_PERIOD_S
|
||||
gdop = telemetry.gdop_at(time_s)
|
||||
satellites = max(4, min(8, telemetry.tracked_at(time_s)))
|
||||
tdop = gdop * (BOARD_TDOP / BOARD_GDOP)
|
||||
pdop = math.sqrt(max(0.0, gdop * gdop - tdop * tdop))
|
||||
vdop = 0.65 * pdop
|
||||
hdop = math.sqrt(max(0.0, pdop * pdop - vdop * vdop))
|
||||
cn0_mean = (
|
||||
BOARD_CN0_MEAN_DBHZ
|
||||
+ 1.2 * math.sin(2.0 * math.pi * time_s / 93.0)
|
||||
+ rng.gauss(0.0, 0.35)
|
||||
)
|
||||
cn0_min = cn0_mean - (
|
||||
BOARD_CN0_MEAN_DBHZ - BOARD_CN0_MIN_DBHZ
|
||||
)
|
||||
chi_square = 0.35
|
||||
velocity_residual = 0.23
|
||||
drift_sigma = BOARD_DRIFT_SIGMA_MPS
|
||||
if scenario.geometry_transitions and _weak_geometry(time_s):
|
||||
gdop = 8.0
|
||||
tdop = 5.0
|
||||
hdop = 4.3
|
||||
vdop = math.sqrt(max(0.0, gdop * gdop - tdop * tdop - hdop * hdop))
|
||||
satellites = 5
|
||||
cn0_mean = 27.0
|
||||
cn0_min = 18.0
|
||||
chi_square = 0.85
|
||||
velocity_residual = 1.8
|
||||
clock_sigma = 65.0
|
||||
drift_sigma = 5.0
|
||||
uere_sigma = clock_sigma / tdop
|
||||
else:
|
||||
uere_sigma = (
|
||||
(BOARD_CLOCK_SIGMA_M / BOARD_TDOP)
|
||||
* 10.0 ** ((BOARD_CN0_MEAN_DBHZ - cn0_mean) / 20.0)
|
||||
)
|
||||
clock_sigma = tdop * uere_sigma
|
||||
|
||||
pseudorange_errors = [
|
||||
rng.gauss(0.0, uere_sigma) for _ in range(satellites)
|
||||
]
|
||||
pseudorange_rms = math.sqrt(
|
||||
statistics.mean(value * value for value in pseudorange_errors)
|
||||
)
|
||||
long_omega = 2.0 * math.pi / 75.0
|
||||
short_omega = 2.0 * math.pi / 19.0
|
||||
true_datum_drift = (
|
||||
5.0 * math.sin(long_omega * time_s)
|
||||
+ 2.0 * math.sin(short_omega * time_s)
|
||||
)
|
||||
true_datum_offset = (
|
||||
5.0 * (1.0 - math.cos(long_omega * time_s)) / long_omega
|
||||
+ 2.0 * (1.0 - math.cos(short_omega * time_s)) / short_omega
|
||||
)
|
||||
timing_error_s = rng.gauss(0.0, clock_sigma / SPEED_OF_LIGHT_MPS)
|
||||
if scenario.multipath:
|
||||
timing_error_s += (
|
||||
25.0e-9 * math.sin(2.0 * math.pi * time_s / 67.0)
|
||||
+ 15.0e-9
|
||||
* math.sin(2.0 * math.pi * time_s / 29.0 + 0.8)
|
||||
)
|
||||
drift_error = true_datum_drift + rng.gauss(0.0, drift_sigma)
|
||||
oscillator_ppm = BOARD_OSCILLATOR_PPM + (
|
||||
scenario.thermal_amplitude_ppm
|
||||
* math.sin(2.0 * math.pi * time_s / scenario.thermal_period_s)
|
||||
)
|
||||
quality = QualityMetrics(
|
||||
tdop,
|
||||
gdop,
|
||||
clock_sigma,
|
||||
cn0_mean,
|
||||
cn0_min,
|
||||
chi_square,
|
||||
velocity_residual,
|
||||
satellites,
|
||||
)
|
||||
observation = Observation(
|
||||
time_s,
|
||||
timing_error_s,
|
||||
drift_error,
|
||||
quality,
|
||||
True,
|
||||
BOARD_REFERENCE_RANGE_RATE_MPS,
|
||||
1,
|
||||
true_datum_offset,
|
||||
)
|
||||
|
||||
continuous_abs = rf_initial + time_s * RF_HZ
|
||||
raw_unwrapped = math.ceil(continuous_abs - 1.0e-12)
|
||||
fractional_sample = raw_unwrapped - continuous_abs
|
||||
rem_code_phase = round(fractional_sample * CODE_FREQ_HZ)
|
||||
abs_sample = raw_unwrapped & 0xFFFF_FFFF
|
||||
corrected_abs = abs_sample - rem_code_phase / CODE_FREQ_HZ
|
||||
raw_capture_counter = math.floor(
|
||||
time_s * TRACKING_HZ * (1.0 + oscillator_ppm * 1.0e-6)
|
||||
)
|
||||
capture_correction_ticks = round(
|
||||
(rem_code_phase / CODE_FREQ_HZ) * (TRACKING_HZ / RF_HZ)
|
||||
)
|
||||
capture_counter = raw_capture_counter - capture_correction_ticks
|
||||
result.append(
|
||||
SyntheticEpoch(
|
||||
scenario.name,
|
||||
time_s,
|
||||
abs_sample,
|
||||
rem_code_phase,
|
||||
corrected_abs,
|
||||
capture_counter,
|
||||
capture_correction_ticks / TRACKING_HZ * 1.0e9,
|
||||
satellites,
|
||||
gdop,
|
||||
hdop,
|
||||
vdop,
|
||||
tdop,
|
||||
cn0_mean,
|
||||
cn0_min,
|
||||
pseudorange_rms,
|
||||
clock_sigma,
|
||||
drift_error,
|
||||
timing_error_s * 1.0e9,
|
||||
oscillator_ppm,
|
||||
classify_quality(quality).name.lower(),
|
||||
observation,
|
||||
)
|
||||
)
|
||||
return result
|
||||
|
||||
|
||||
def _run_loop(epochs: list[SyntheticEpoch], tuning: Tuning):
|
||||
model = FixedPointAdaptivePps(
|
||||
oscillator_ppm=BOARD_OSCILLATOR_PPM,
|
||||
profiles=tuning.profiles,
|
||||
acquire_profile=tuning.acquire,
|
||||
software_drift_gain=tuning.drift_gain,
|
||||
)
|
||||
for epoch in epochs:
|
||||
# Update at a solve boundary; the selected value applies to the
|
||||
# preceding synthetic interval, whose endpoint is this observation.
|
||||
model.oscillator_fraction = epoch.oscillator_ppm * 1.0e-6
|
||||
if tuning.robust_drift_gain is not None:
|
||||
target = classify_quality(epoch.observation.quality)
|
||||
model.software_drift_gain = (
|
||||
tuning.robust_drift_gain
|
||||
if target != Profile.PRECISION
|
||||
else tuning.drift_gain
|
||||
)
|
||||
model.update(epoch.observation)
|
||||
finish_s = epochs[-1].time_s + UPDATE_PERIOD_S if epochs else 0.0
|
||||
return model.finish(finish_s)
|
||||
|
||||
|
||||
def measure_trace(trace, *, after_s: float) -> RunMetrics:
|
||||
edge_errors = [
|
||||
error for _, error in trace.integer_second_edge_errors(after_s=after_s)
|
||||
]
|
||||
if len(edge_errors) < 12:
|
||||
raise RuntimeError("not enough locked integer-second PPS edges")
|
||||
center = statistics.median(edge_errors)
|
||||
centered = [value - center for value in edge_errors]
|
||||
periods = [
|
||||
edge_errors[index] - edge_errors[index - 1]
|
||||
for index in range(1, len(edge_errors))
|
||||
]
|
||||
moving_10s = [
|
||||
statistics.mean(centered[max(0, index - 9) : index + 1])
|
||||
for index in range(len(centered))
|
||||
]
|
||||
profiles = [
|
||||
profile
|
||||
for time_s, profile in zip(trace.time_s, trace.profile)
|
||||
if time_s >= after_s
|
||||
]
|
||||
profile_transitions = sum(
|
||||
profiles[index] != profiles[index - 1]
|
||||
for index in range(1, len(profiles))
|
||||
)
|
||||
return RunMetrics(
|
||||
trace.first_lock_time() if trace.first_lock_time() is not None else math.inf,
|
||||
math.sqrt(statistics.mean(value * value for value in centered)) * 1.0e9,
|
||||
_percentile((abs(value) for value in centered), 0.95) * 1.0e9,
|
||||
max(abs(value) for value in centered) * 1.0e9,
|
||||
(max(centered) - min(centered)) * 1.0e9,
|
||||
math.sqrt(statistics.mean(value * value for value in periods)) * 1.0e9,
|
||||
_percentile((abs(value) for value in periods), 0.95) * 1.0e9,
|
||||
max(abs(value) for value in periods) * 1.0e9,
|
||||
(max(moving_10s) - min(moving_10s)) * 1.0e9,
|
||||
profile_transitions,
|
||||
trace.unlock_events,
|
||||
)
|
||||
|
||||
|
||||
def aggregate(
|
||||
scenario: Scenario, tuning: Tuning, runs: list[RunMetrics]
|
||||
) -> Aggregate:
|
||||
return Aggregate(
|
||||
scenario.name,
|
||||
tuning.name,
|
||||
len(runs),
|
||||
_percentile((run.first_lock_s for run in runs), 0.95),
|
||||
statistics.median(run.absolute_rms_ns for run in runs),
|
||||
_percentile((run.absolute_max_ns for run in runs), 0.95),
|
||||
statistics.median(run.absolute_peak_to_peak_ns for run in runs),
|
||||
statistics.median(run.period_rms_ns for run in runs),
|
||||
statistics.median(run.period_p95_ns for run in runs),
|
||||
max(run.period_max_ns for run in runs),
|
||||
statistics.median(run.wander_10s_peak_to_peak_ns for run in runs),
|
||||
sum(run.profile_transitions for run in runs),
|
||||
sum(run.unlocks for run in runs),
|
||||
)
|
||||
|
||||
|
||||
def discrete_poles(profile: LoopProfile) -> tuple[complex, complex]:
|
||||
"""Poles of the exact sample/commit ordering, before slew saturation."""
|
||||
kp_q46, ki_q46, _ = profile.fixed_coefficients(UPDATE_PERIOD_S)
|
||||
kp = kp_q46 / (1 << 46)
|
||||
ki = ki_q46 / (1 << 46)
|
||||
ticks = UPDATE_PERIOD_S * TRACKING_HZ
|
||||
a = 1.0 - ticks * (kp + ki)
|
||||
b = ticks
|
||||
c = -ki
|
||||
d = 1.0
|
||||
trace = a + d
|
||||
determinant = a * d - b * c
|
||||
discriminant = complex(trace * trace - 4.0 * determinant, 0.0) ** 0.5
|
||||
return (trace + discriminant) / 2.0, (trace - discriminant) / 2.0
|
||||
|
||||
|
||||
def write_telemetry(path: Path, epochs: list[SyntheticEpoch]) -> None:
|
||||
fields = [
|
||||
"scenario", "time_s", "abs_sample", "rem_code_phase",
|
||||
"corrected_abs_sample", "capture_counter_125", "used_satellites",
|
||||
"capture_correction_ns",
|
||||
"gdop", "hdop", "vdop", "tdop", "cn0_mean_dbhz",
|
||||
"cn0_min_dbhz", "pseudorange_noise_rms_m", "clock_bias_sigma_m",
|
||||
"clock_drift_error_mps", "timing_error_ns", "oscillator_ppm",
|
||||
"quality_profile",
|
||||
]
|
||||
path.parent.mkdir(parents=True, exist_ok=True)
|
||||
with path.open("w", encoding="utf-8", newline="") as output:
|
||||
writer = csv.DictWriter(
|
||||
output, fieldnames=fields, lineterminator="\n"
|
||||
)
|
||||
writer.writeheader()
|
||||
for epoch in epochs:
|
||||
writer.writerow({field: getattr(epoch, field) for field in fields})
|
||||
|
||||
|
||||
def write_sweep(path: Path, results: list[Aggregate]) -> None:
|
||||
path.parent.mkdir(parents=True, exist_ok=True)
|
||||
with path.open("w", encoding="utf-8", newline="") as output:
|
||||
writer = csv.DictWriter(
|
||||
output,
|
||||
fieldnames=list(Aggregate.__annotations__),
|
||||
lineterminator="\n",
|
||||
)
|
||||
writer.writeheader()
|
||||
for item in results:
|
||||
writer.writerow(item.__dict__)
|
||||
|
||||
|
||||
def _fmt(value: float) -> str:
|
||||
return "n/a" if not math.isfinite(value) else f"{value:.2f}"
|
||||
|
||||
|
||||
def write_report(
|
||||
path: Path,
|
||||
telemetry: RealTelemetry,
|
||||
results: list[Aggregate],
|
||||
epoch_result,
|
||||
) -> None:
|
||||
gdop = telemetry.gdop
|
||||
tracked = telemetry.tracked
|
||||
current = [item for item in results if item.tuning == "current"]
|
||||
guarded = [
|
||||
item for item in results if item.tuning == "robust_noise_guard"
|
||||
]
|
||||
lines = [
|
||||
"# PPS stability analysis: current RTL + 180 ms firmware model",
|
||||
"",
|
||||
"This report is generated by `stability_analysis.py`. The absolute "
|
||||
"LNAV epoch path and the stochastic control loop are evaluated "
|
||||
"separately so a 20 ms semantic error cannot be hidden by loop noise.",
|
||||
"",
|
||||
"## Absolute epoch result",
|
||||
"",
|
||||
f"- Current RTL delivered the validated subframe after "
|
||||
f"{epoch_result.first_subframe_readout_bits} following bits.",
|
||||
f"- With 180 ms postponement, reference error was "
|
||||
f"{epoch_result.first_reference_error_ms:.6f} ms.",
|
||||
f"- DUT minus trusted PPS was "
|
||||
f"{epoch_result.pps.edges[0].dut_minus_trusted_s * 1.0e9:.2f} ns. "
|
||||
"That remainder is the model's configured four-PCLK capture pipeline, "
|
||||
"not a measured board delay.",
|
||||
"- The production source still defines `POSTPONE_TOW` as 160 ms; this "
|
||||
"analysis uses the requested 180 ms candidate configuration.",
|
||||
"",
|
||||
"## Repository hardware data",
|
||||
"",
|
||||
f"- Real GPS capture: {len(gdop)} GDOP samples and {len(tracked)} "
|
||||
"tracked-satellite samples over about five minutes.",
|
||||
f"- GDOP min/median/p95/max: {_fmt(min(gdop))} / "
|
||||
f"{_fmt(statistics.median(gdop))} / {_fmt(_percentile(gdop, 0.95))} / "
|
||||
f"{_fmt(max(gdop))}.",
|
||||
f"- Tracked satellites: {min(tracked)}..{max(tracked)} "
|
||||
f"(median {statistics.median(tracked):.0f}).",
|
||||
f"- Horizontal scatter RMS/p95/max: "
|
||||
f"{_fmt(telemetry.horizontal_rms_m)} / "
|
||||
f"{_fmt(telemetry.horizontal_p95_m)} / "
|
||||
f"{_fmt(telemetry.horizontal_max_m)} m. No surveyed truth was recorded.",
|
||||
"- No real ABS_SAMPLE, REM_CODE_PHASE, TDOP/HDOP, C/N0, pseudorange "
|
||||
"residual, PPS phase/rate/integrator, temperature, or trusted-edge TIC "
|
||||
"samples exist in the repository. Those columns in the companion CSV "
|
||||
"are explicitly synthetic.",
|
||||
"",
|
||||
"Synthetic timing uses the board-calibration values already encoded "
|
||||
"in the PPS tests: TDOP 1.50, GDOP 2.92, clock sigma 15.87 m "
|
||||
"(52.94 ns), C/N0 40.3/34.8 dB-Hz, drift noise 3.47 m/s, seven-satellite "
|
||||
"reference point, -291.75 m/s datum range rate, and 1.592 ppm local "
|
||||
"oscillator. Real GDOP and satellite count replace the geometry fields "
|
||||
"when available.",
|
||||
"",
|
||||
"## Current-loop results",
|
||||
"",
|
||||
"`period max` is the worst |edge[n]-edge[n-1]| across all seeds. "
|
||||
"`absolute p-p` and `10 s wander p-p` expose the slow left/right "
|
||||
"motion that period RMS hides.",
|
||||
"",
|
||||
"| scenario | abs RMS median (ns) | abs p-p median (ns) | "
|
||||
"10 s wander p-p (ns) | period RMS (ns) | period max (ns) | unlocks |",
|
||||
"|---|---:|---:|---:|---:|---:|---:|",
|
||||
]
|
||||
for item in current:
|
||||
lines.append(
|
||||
f"| {item.scenario} | {item.absolute_rms_median_ns:.2f} | "
|
||||
f"{item.absolute_peak_to_peak_median_ns:.2f} | "
|
||||
f"{item.wander_10s_peak_to_peak_median_ns:.2f} | "
|
||||
f"{item.period_rms_median_ns:.2f} | "
|
||||
f"{item.period_max_worst_ns:.2f} | {item.unlocks_total} |"
|
||||
)
|
||||
lines.extend(
|
||||
[
|
||||
"",
|
||||
"The loop does not numerically oscillate or lose lock in these "
|
||||
"runs. The scope-like symptom is low-frequency phase wander: it "
|
||||
"can be tens of nanoseconds while one-second period jitter remains "
|
||||
"below 20 ns. Therefore a `p95 period jitter < 20 ns` assertion is "
|
||||
"not sufficient for an absolute ±20 ns requirement.",
|
||||
"",
|
||||
"## Linear stability",
|
||||
"",
|
||||
"The fixed-point coefficients give the following discrete poles "
|
||||
"at the 200 ms navigation update rate:",
|
||||
"",
|
||||
"| profile | natural frequency | damping | poles | max magnitude |",
|
||||
"|---|---:|---:|---|---:|",
|
||||
]
|
||||
)
|
||||
profile_rows = [("acquire", ACQUIRE_PROFILE)] + [
|
||||
(name, PROFILES[profile])
|
||||
for name, profile in (
|
||||
("precision", Profile.PRECISION),
|
||||
("robust", Profile.ROBUST),
|
||||
)
|
||||
]
|
||||
for name, profile in profile_rows:
|
||||
poles = discrete_poles(profile)
|
||||
pole_text = ", ".join(
|
||||
f"{pole.real:.6f}{pole.imag:+.6f}j" for pole in poles
|
||||
)
|
||||
lines.append(
|
||||
f"| {name} | {profile.natural_frequency_hz:.3f} Hz | "
|
||||
f"{profile.damping:.2f} | {pole_text} | "
|
||||
f"{max(abs(pole) for pole in poles):.6f} |"
|
||||
)
|
||||
lines.extend(
|
||||
[
|
||||
"",
|
||||
"All magnitudes are below one. Profile changes preserve the PI "
|
||||
"integrator, and the RTL snapshots Kp/Ki/max-slew/lock limits "
|
||||
"together with the next commit, so the current acquire/precision/"
|
||||
"robust switching is structurally sound. Stable individual modes "
|
||||
"do not prove arbitrary-switching stability, but the asymmetric "
|
||||
"8-update upgrade / 3-update downgrade hysteresis and the stress "
|
||||
"runs show no coefficient chatter or unlock.",
|
||||
"",
|
||||
"## Dynamic parameter sweep",
|
||||
"",
|
||||
"| tuning | worst period max (ns) | median abs RMS over scenarios "
|
||||
"(ns) | transitions | unlocks |",
|
||||
"|---|---:|---:|---:|---:|",
|
||||
]
|
||||
)
|
||||
for tuning in tuning_matrix():
|
||||
subset = [item for item in results if item.tuning == tuning.name]
|
||||
median_scenario_rms = statistics.median(
|
||||
item.absolute_rms_median_ns for item in subset
|
||||
)
|
||||
lines.append(
|
||||
f"| {tuning.name} | "
|
||||
f"{max(item.period_max_worst_ns for item in subset):.2f} | "
|
||||
f"{median_scenario_rms:.2f} | "
|
||||
f"{sum(item.profile_transitions_total for item in subset)} | "
|
||||
f"{sum(item.unlocks_total for item in subset)} |"
|
||||
)
|
||||
current_transition = next(
|
||||
item for item in current if item.scenario == "quality_transitions"
|
||||
)
|
||||
guarded_transition = next(
|
||||
item for item in guarded if item.scenario == "quality_transitions"
|
||||
)
|
||||
scheduled = [
|
||||
item for item in results
|
||||
if item.tuning == "quality_scheduled_drift"
|
||||
]
|
||||
lines.extend(
|
||||
[
|
||||
"",
|
||||
"The safest small candidate is `robust_noise_guard`: retain the "
|
||||
"current 0.035 Hz precision loop and 0.012 Hz robust loop, but "
|
||||
"reduce only the robust-profile WLS residual gain from 1/32 to "
|
||||
"1/64. In the synthetic quality-transition case its worst period "
|
||||
f"step changed from {current_transition.period_max_worst_ns:.2f} "
|
||||
f"to {guarded_transition.period_max_worst_ns:.2f} ns. This is a "
|
||||
"candidate for hardware A/B testing, not enough evidence for a "
|
||||
"production constant change.",
|
||||
"",
|
||||
"The best synthetic all-scenario period result is "
|
||||
"`quality_scheduled_drift`: use drift gain 1/2 only with a "
|
||||
"precision-quality Doppler solution, fall back to the current 1/4 "
|
||||
"gain in robust quality, and combine that with the 1/64 robust "
|
||||
"bias gain. Its worst period step was "
|
||||
f"{max(item.period_max_worst_ns for item in scheduled):.2f} ns "
|
||||
"with no unlocks. It is the most meaningful dynamic A/B candidate, "
|
||||
"but the repository has no real drift-covariance or TIC samples "
|
||||
"with which to validate the threshold.",
|
||||
"",
|
||||
"Narrowing precision to 0.012 Hz improves stationary noise but "
|
||||
"lags even the explicit 0.05 ppm/120 s synthetic oscillator "
|
||||
"variation. Widening it improves oscillator tracking but passes "
|
||||
"more WLS/multipath motion. A geometry-only bandwidth scheduler "
|
||||
"cannot resolve that tradeoff.",
|
||||
"",
|
||||
"## Recommended architecture update",
|
||||
"",
|
||||
"1. Keep current acquire/precision/robust PI profiles until real "
|
||||
"TIC and oscillator-temperature data are available. Do not tune "
|
||||
"from scope appearance alone.",
|
||||
"2. Log the existing firmware status line at every solve during "
|
||||
"validation: ABS_SAMPLE, REM_CODE_PHASE, capture sequence/counter, "
|
||||
"TDOP/HDOP/GDOP, C/N0, pseudorange residual RMS, raw/filtered bias, "
|
||||
"innovation, bias drift, hardware phase, rate, integrator, profile, "
|
||||
"lock, FPGA/oscillator temperature, DUT edge and trusted edge.",
|
||||
"3. A/B test the simulated scheduled filter first: precision "
|
||||
"drift gain 1/2, robust drift gain 1/4, precision bias gain 1/32, "
|
||||
"and robust bias gain 1/64. Longer term, use estimated phase/drift "
|
||||
"covariance to interpolate those software gains within bounds. "
|
||||
"Keep RTL PI bandwidth changes hysteretic and preserve its "
|
||||
"integrator as today.",
|
||||
"4. If temperature/frequency acceleration is the main source, add "
|
||||
"oscillator temperature feed-forward or a phase/frequency/frequency-"
|
||||
"drift estimator. Switching PI bandwidth alone either passes WLS "
|
||||
"wander or lags the oscillator.",
|
||||
"5. Define acceptance separately: calibrated mean offset, absolute "
|
||||
"centered error/p-p wander, one-second period error, holdover, and "
|
||||
"re-lock. A physical max-20 ns claim requires a common-clock TIC "
|
||||
"run after cable/output-delay calibration; this synthetic test "
|
||||
"cannot certify it.",
|
||||
"",
|
||||
]
|
||||
)
|
||||
path.parent.mkdir(parents=True, exist_ok=True)
|
||||
path.write_text("\n".join(lines), encoding="utf-8")
|
||||
|
||||
|
||||
def parse_args() -> argparse.Namespace:
|
||||
parser = argparse.ArgumentParser(description=__doc__)
|
||||
parser.add_argument("--duration", type=float, default=300.0)
|
||||
parser.add_argument("--seeds", type=int, default=24)
|
||||
parser.add_argument(
|
||||
"--output-dir",
|
||||
type=Path,
|
||||
default=Path(__file__).resolve().parent / "output/stability",
|
||||
)
|
||||
parser.add_argument("--no-write", action="store_true")
|
||||
return parser.parse_args()
|
||||
|
||||
|
||||
def main() -> int:
|
||||
args = parse_args()
|
||||
if args.duration < 120.0 or args.seeds < 1:
|
||||
raise SystemExit("duration must be >= 120 s and seeds must be positive")
|
||||
software_dir = Path(__file__).resolve().parents[1]
|
||||
telemetry = load_real_telemetry(software_dir)
|
||||
if not telemetry.gdop or not telemetry.tracked:
|
||||
raise SystemExit("real GPS GDOP/satellite capture was not found")
|
||||
|
||||
epoch_result = run_end_to_end(
|
||||
EndToEndConfig(
|
||||
decoder_mode=DecoderMode.CURRENT_RTL,
|
||||
postpone_tow_s=0.180,
|
||||
)
|
||||
)
|
||||
analysis_start = max(60.0, args.duration / 3.0)
|
||||
results: list[Aggregate] = []
|
||||
telemetry_epochs: list[SyntheticEpoch] = []
|
||||
for scenario in scenario_matrix():
|
||||
seed_epochs = [
|
||||
synthesize_epochs(
|
||||
telemetry,
|
||||
scenario,
|
||||
duration_s=args.duration,
|
||||
seed=seed,
|
||||
)
|
||||
for seed in range(1, args.seeds + 1)
|
||||
]
|
||||
if scenario.name == "scope_like_multipath":
|
||||
telemetry_epochs = seed_epochs[0]
|
||||
for tuning in tuning_matrix():
|
||||
runs = [
|
||||
measure_trace(
|
||||
_run_loop(epochs, tuning), after_s=analysis_start
|
||||
)
|
||||
for epochs in seed_epochs
|
||||
]
|
||||
results.append(aggregate(scenario, tuning, runs))
|
||||
|
||||
print(
|
||||
"scenario tuning abs_rms "
|
||||
"abs_pp wander10 period_rms period_max unlock"
|
||||
)
|
||||
for item in results:
|
||||
print(
|
||||
f"{item.scenario:27s} {item.tuning:22s} "
|
||||
f"{item.absolute_rms_median_ns:7.2f} "
|
||||
f"{item.absolute_peak_to_peak_median_ns:7.2f} "
|
||||
f"{item.wander_10s_peak_to_peak_median_ns:8.2f} "
|
||||
f"{item.period_rms_median_ns:10.2f} "
|
||||
f"{item.period_max_worst_ns:10.2f} "
|
||||
f"{item.unlocks_total:6d}"
|
||||
)
|
||||
print(
|
||||
"current RTL + 180 ms: DUT-trusted = "
|
||||
f"{epoch_result.pps.edges[0].dut_minus_trusted_s * 1.0e9:.2f} ns"
|
||||
)
|
||||
if not args.no_write:
|
||||
write_telemetry(args.output_dir / "synthetic_telemetry.csv", telemetry_epochs)
|
||||
write_sweep(args.output_dir / "parameter_sweep.csv", results)
|
||||
write_report(
|
||||
args.output_dir / "stability_report.md",
|
||||
telemetry,
|
||||
results,
|
||||
epoch_result,
|
||||
)
|
||||
print(f"wrote {args.output_dir}")
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
@@ -0,0 +1,84 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from pathlib import Path
|
||||
import unittest
|
||||
|
||||
from pps_python.adaptive_pps import ACQUIRE_PROFILE, PROFILES, Profile
|
||||
|
||||
from pps_end_to_end_validation.stability_analysis import (
|
||||
discrete_poles,
|
||||
load_real_telemetry,
|
||||
measure_trace,
|
||||
scenario_matrix,
|
||||
synthesize_epochs,
|
||||
tuning_matrix,
|
||||
_run_loop,
|
||||
)
|
||||
|
||||
|
||||
SOFTWARE_DIR = Path(__file__).resolve().parents[1]
|
||||
|
||||
|
||||
class StabilityAnalysisTests(unittest.TestCase):
|
||||
@classmethod
|
||||
def setUpClass(cls):
|
||||
cls.telemetry = load_real_telemetry(SOFTWARE_DIR)
|
||||
|
||||
def test_repository_capture_is_loaded(self):
|
||||
self.assertEqual(len(self.telemetry.gdop), 989)
|
||||
self.assertEqual(len(self.telemetry.tracked), 990)
|
||||
self.assertEqual(set(self.telemetry.tracked), {8})
|
||||
self.assertAlmostEqual(min(self.telemetry.gdop), 1.8593929921)
|
||||
self.assertAlmostEqual(max(self.telemetry.gdop), 3.1537716348)
|
||||
self.assertGreater(self.telemetry.horizontal_rms_m, 8.0)
|
||||
self.assertLess(self.telemetry.horizontal_rms_m, 9.0)
|
||||
|
||||
def test_fixed_point_profiles_have_inside_unit_circle_poles(self):
|
||||
for profile in (
|
||||
ACQUIRE_PROFILE,
|
||||
PROFILES[Profile.PRECISION],
|
||||
PROFILES[Profile.ROBUST],
|
||||
):
|
||||
poles = discrete_poles(profile)
|
||||
self.assertTrue(all(abs(pole) < 1.0 for pole in poles))
|
||||
|
||||
def test_synthetic_telemetry_exercises_wrap_and_physical_fields(self):
|
||||
epochs = synthesize_epochs(
|
||||
self.telemetry,
|
||||
scenario_matrix()[0],
|
||||
duration_s=120.0,
|
||||
seed=1,
|
||||
)
|
||||
self.assertEqual(len(epochs), 600)
|
||||
self.assertTrue(any(epochs[i].abs_sample < epochs[i - 1].abs_sample
|
||||
for i in range(1, len(epochs))))
|
||||
self.assertTrue(all(0 <= epoch.rem_code_phase <= 1_023_000
|
||||
for epoch in epochs))
|
||||
self.assertTrue(all(epoch.used_satellites == 8 for epoch in epochs))
|
||||
self.assertTrue(all(math.isfinite(epoch.corrected_abs_sample)
|
||||
for epoch in epochs))
|
||||
self.assertTrue(all(epoch.capture_correction_ns >= 0.0
|
||||
for epoch in epochs))
|
||||
|
||||
def test_current_loop_locks_without_unlock_in_board_scenario(self):
|
||||
epochs = synthesize_epochs(
|
||||
self.telemetry,
|
||||
scenario_matrix()[0],
|
||||
duration_s=120.0,
|
||||
seed=4,
|
||||
)
|
||||
trace = _run_loop(epochs, tuning_matrix()[0])
|
||||
metrics = measure_trace(trace, after_s=60.0)
|
||||
self.assertLess(metrics.first_lock_s, 15.0)
|
||||
self.assertEqual(metrics.unlocks, 0)
|
||||
self.assertLessEqual(metrics.period_max_ns, 16.0 + 1.0e-9)
|
||||
# Slow absolute motion is deliberately measured separately.
|
||||
self.assertGreater(
|
||||
metrics.absolute_peak_to_peak_ns,
|
||||
metrics.period_max_ns,
|
||||
)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
@@ -11,7 +11,7 @@ from __future__ import annotations
|
||||
from dataclasses import dataclass
|
||||
from enum import IntEnum
|
||||
import math
|
||||
from typing import Iterable
|
||||
from typing import Iterable, Mapping
|
||||
|
||||
|
||||
TRACKING_HZ = 125_000_000
|
||||
@@ -242,17 +242,7 @@ class AdaptiveTrace:
|
||||
|
||||
def one_second_edge_jitter_rms(self, *, after_s: float = 0.0) -> float:
|
||||
"""Return 1 s PPS-edge jitter including the 125 MHz output grid."""
|
||||
selected = {
|
||||
round(time_s, 9): (
|
||||
round(-phase_s * TRACKING_HZ) / TRACKING_HZ
|
||||
)
|
||||
for time_s, phase_s, locked in zip(
|
||||
self.time_s, self.phase_error_s, self.locked
|
||||
)
|
||||
if locked
|
||||
and time_s >= after_s
|
||||
and abs(time_s - round(time_s)) < 1.0e-9
|
||||
}
|
||||
selected = dict(self.integer_second_edge_errors(after_s=after_s))
|
||||
differences = [
|
||||
edge_error - selected[round(time_s - 1.0, 9)]
|
||||
for time_s, edge_error in selected.items()
|
||||
@@ -264,6 +254,31 @@ class AdaptiveTrace:
|
||||
sum(value * value for value in differences) / len(differences)
|
||||
)
|
||||
|
||||
def integer_second_edge_errors(
|
||||
self,
|
||||
*,
|
||||
after_s: float = 0.0,
|
||||
locked_only: bool = True,
|
||||
) -> list[tuple[float, float]]:
|
||||
"""Return physical PPS edge error samples on the 125 MHz grid.
|
||||
|
||||
Period jitter alone hides low-frequency motion. Exposing the signed
|
||||
edge error lets stability analyses report absolute RMS, peak-to-peak
|
||||
wander and moving-average wander alongside cycle-to-cycle jitter.
|
||||
"""
|
||||
return [
|
||||
(
|
||||
round(time_s, 9),
|
||||
round(-phase_s * TRACKING_HZ) / TRACKING_HZ,
|
||||
)
|
||||
for time_s, phase_s, locked in zip(
|
||||
self.time_s, self.phase_error_s, self.locked
|
||||
)
|
||||
if (locked or not locked_only)
|
||||
and time_s >= after_s
|
||||
and abs(time_s - round(time_s)) < 1.0e-9
|
||||
]
|
||||
|
||||
def phase_rms(self, *, after_s: float = 0.0) -> float:
|
||||
values = [
|
||||
phase_s
|
||||
@@ -319,6 +334,10 @@ class FixedPointAdaptivePps:
|
||||
software_filter: bool = True,
|
||||
receiver_clock_bias_m: float = 1200.0,
|
||||
receiver_clock_drift_mps: float = 475.0,
|
||||
profiles: Mapping[Profile, LoopProfile] | None = None,
|
||||
acquire_profile: LoopProfile | None = None,
|
||||
software_drift_gain: float = 1.0 / 4.0,
|
||||
software_fallback_drift_gain: float = 1.0 / 16.0,
|
||||
) -> None:
|
||||
self.update_period_s = update_period_s
|
||||
self.oscillator_fraction = oscillator_ppm * 1.0e-6
|
||||
@@ -326,9 +345,24 @@ class FixedPointAdaptivePps:
|
||||
self.software_filter = software_filter
|
||||
self.receiver_clock_bias_m = receiver_clock_bias_m
|
||||
self.receiver_clock_drift_mps = receiver_clock_drift_mps
|
||||
self.profiles = dict(PROFILES if profiles is None else profiles)
|
||||
missing_profiles = set(Profile) - self.profiles.keys()
|
||||
if missing_profiles:
|
||||
raise ValueError(f"missing loop profiles: {sorted(missing_profiles)}")
|
||||
self.acquire_profile = (
|
||||
ACQUIRE_PROFILE if acquire_profile is None else acquire_profile
|
||||
)
|
||||
if not 0.0 < software_drift_gain <= 1.0:
|
||||
raise ValueError("software_drift_gain must be in (0, 1]")
|
||||
if not 0.0 < software_fallback_drift_gain <= 1.0:
|
||||
raise ValueError(
|
||||
"software_fallback_drift_gain must be in (0, 1]"
|
||||
)
|
||||
self.software_drift_gain = software_drift_gain
|
||||
self.software_fallback_drift_gain = software_fallback_drift_gain
|
||||
self.selector = QualitySelector()
|
||||
self.profile = Profile.HOLDOVER if adaptive else Profile.PRECISION
|
||||
self.parameters = PROFILES[self.profile]
|
||||
self.parameters = self.profiles[self.profile]
|
||||
self.kp_q46, self.ki_q46, self.max_slew_q32 = (
|
||||
self.parameters.fixed_coefficients(update_period_s)
|
||||
)
|
||||
@@ -376,10 +410,10 @@ class FixedPointAdaptivePps:
|
||||
|
||||
def _select_loop_parameters(self) -> None:
|
||||
parameters = (
|
||||
ACQUIRE_PROFILE
|
||||
self.acquire_profile
|
||||
if self.adaptive and
|
||||
self.profile != Profile.HOLDOVER and not self.locked
|
||||
else PROFILES[self.profile]
|
||||
else self.profiles[self.profile]
|
||||
)
|
||||
if parameters == self.parameters:
|
||||
return
|
||||
@@ -541,15 +575,16 @@ class FixedPointAdaptivePps:
|
||||
+ self.parameters.software_bias_gain * innovation_m
|
||||
)
|
||||
if observation.clock_drift_valid:
|
||||
self.software_clock_drift_mps += 0.25 * (
|
||||
self.software_clock_drift_mps += self.software_drift_gain * (
|
||||
measured_drift_mps - self.software_clock_drift_mps
|
||||
)
|
||||
else:
|
||||
if self.software_updates < self.parameters.software_settle_updates:
|
||||
self.software_clock_drift_mps = fallback_drift_mps
|
||||
else:
|
||||
self.software_clock_drift_mps += (1.0 / 16.0) * (
|
||||
fallback_drift_mps - self.software_clock_drift_mps
|
||||
self.software_clock_drift_mps += (
|
||||
self.software_fallback_drift_gain *
|
||||
(fallback_drift_mps - self.software_clock_drift_mps)
|
||||
)
|
||||
self.software_previous_measurement_m = measured_bias_m
|
||||
self.software_error_s = -(
|
||||
@@ -705,6 +740,10 @@ def run_observations(
|
||||
software_filter: bool = True,
|
||||
receiver_clock_bias_m: float = 1200.0,
|
||||
receiver_clock_drift_mps: float = 475.0,
|
||||
profiles: Mapping[Profile, LoopProfile] | None = None,
|
||||
acquire_profile: LoopProfile | None = None,
|
||||
software_drift_gain: float = 1.0 / 4.0,
|
||||
software_fallback_drift_gain: float = 1.0 / 16.0,
|
||||
) -> AdaptiveTrace:
|
||||
model = FixedPointAdaptivePps(
|
||||
oscillator_ppm=oscillator_ppm,
|
||||
@@ -712,6 +751,10 @@ def run_observations(
|
||||
software_filter=software_filter,
|
||||
receiver_clock_bias_m=receiver_clock_bias_m,
|
||||
receiver_clock_drift_mps=receiver_clock_drift_mps,
|
||||
profiles=profiles,
|
||||
acquire_profile=acquire_profile,
|
||||
software_drift_gain=software_drift_gain,
|
||||
software_fallback_drift_gain=software_fallback_drift_gain,
|
||||
)
|
||||
for observation in observations:
|
||||
model.update(observation)
|
||||
|
||||
Reference in New Issue
Block a user