awe/radar_system
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added diagnostics

Browse Source
This commit is contained in:
Ayzen
2026-06-22 16:30:55 +03:00
parent 7c381facaf
commit d7ddca6e76
3 changed files with 320 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files
python_app/scripts/kamil_adc_calibrate.py
+318
View File
@@ -37,6 +37,10 @@ import time
import numpy as np
from python_app.hardware_full.kamil_adc import KamilAdcService, apply_kamil_adc_laser_control
from python_app.hardware_full.kamil_adc.processing import (
KamilAdcProcessingParams,
KamilAdcSweepProcessor,
)
from python_app.models.run_config_model import RunConfigModel
logger = logging.getLogger("kamil_adc_calibrate")
@@ -200,6 +204,293 @@ def _report(result: dict) -> None:
print("=" * 72 + "\n")
_TWO_PI = 2.0 * np.pi
def _capture_raw_sweeps(
service: KamilAdcService, *, warmup: int, sweeps: int
) -> list[tuple[np.ndarray, np.ndarray, np.ndarray]]:
"""Capture full ``(steps, main, reference)`` arrays after discarding `warmup`.
Unlike :func:`_capture_reference_phases`, this keeps the complete complex main
and reference samples (and their step indices) so the diagnosis can study the
raw phase, amplitude, and the pass-through behaviour — not just the endpoints.
"""
for index in range(warmup):
service.read_raw_sweep()
if index == 0:
logger.info("Warming up (%d sweeps) while the sweep settles...", warmup)
captured: list[tuple[np.ndarray, np.ndarray, np.ndarray]] = []
for index in range(sweeps):
raw = service.read_raw_sweep()
if raw.reference.size >= 2 and raw.main.size == raw.reference.size:
captured.append(
(
np.asarray(raw.steps),
raw.main.astype(np.complex128),
raw.reference.astype(np.complex128),
)
)
if (index + 1) % 50 == 0:
logger.info("Captured %d/%d sweeps", index + 1, sweeps)
return captured
def _diagnose(
sweeps: list[tuple[np.ndarray, np.ndarray, np.ndarray]], config: RunConfigModel
) -> dict:
"""Characterise how the reference phase / frequency axis moves across sweeps.
The reference arm has a fixed electrical delay, so each sweep's unwrapped phase
is a straight ramp; the sweep start frequency wandering shifts that ramp
vertically (the *physical float*). When ``np.angle(ref[0])`` crosses the +/-pi
cut, ``np.unwrap`` re-anchors the whole ramp one 2*pi turn away (the *spurious
branch wrap*) — a ~``2*pi*slope`` jump of the entire frequency axis that wrecks
the band mapping for that sweep.
This separates the two: it reconstructs the anchor phase per sweep, unwraps it
*across* sweeps (the candidate fix), and reports how that removes discrete 2*pi
steps while leaving the slow physical float intact. It also runs the live
pass-through (``KamilAdcSweepProcessor``) to flag which sweeps actually come out
distorted, and correlates those with the detected branch wraps.
"""
kamil = config.radar.kamil_adc
cal = kamil.phase_calibration
phase0, phase1 = float(cal.phase0_rad), float(cal.phase1_rad)
freq0, freq1 = float(cal.freq0_hz), float(cal.freq1_hz)
slope = (freq1 - freq0) / (phase1 - phase0) # Hz per rad
band_start, band_stop = float(kamil.band.start_hz), float(kamil.band.stop_hz)
n = len(sweeps)
anchor_phase = np.empty(n) # unwrap(angle(ref))[0] == angle(ref[0]) in (-pi, pi]
stop_phase = np.empty(n) # unwrap(angle(ref))[-1]
span_rad = np.empty(n) # stop_phase - anchor_phase (fixed-delay invariant)
n_points = np.empty(n, dtype=np.int64)
f_lo = np.empty(n) # min reconstructed frequency (config calibration)
f_hi = np.empty(n) # max reconstructed frequency
ref_amp_med = np.empty(n)
for i, (_steps, _main, reference) in enumerate(sweeps):
phase = np.unwrap(np.angle(reference))
anchor_phase[i] = phase[0]
stop_phase[i] = phase[-1]
span_rad[i] = phase[-1] - phase[0]
n_points[i] = reference.size
freqs = freq0 + (phase - phase0) * slope
f_lo[i] = float(np.min(freqs))
f_hi[i] = float(np.max(freqs))
ref_amp_med[i] = float(np.median(np.abs(reference)))
# --- Candidate fix: unwrap the anchor phase ACROSS the sweep-time axis ------
# The physical float moves the anchor smoothly; a branch wrap injects a +/-2*pi
# step. Unwrapping along sweep index removes the discrete steps and keeps the
# slow float. branch_turns is the integer turns each sweep was wrapped by.
anchor_unwrapped = np.unwrap(anchor_phase)
branch_turns = np.round((anchor_unwrapped - anchor_phase) / _TWO_PI).astype(np.int64)
axis_shift = branch_turns * _TWO_PI * slope # Hz the whole axis was displaced
f_lo_corr = f_lo + axis_shift
f_hi_corr = f_hi + axis_shift
# Branch wraps: consecutive anchor steps above pi are spurious 2*pi jumps.
anchor_step = np.diff(anchor_phase)
wrap_indices = np.nonzero(np.abs(anchor_step) > np.pi)[0] + 1 # sweep index after wrap
covers_raw = (f_lo <= band_start) & (f_hi >= band_stop)
covers_corr = (f_lo_corr <= band_start) & (f_hi_corr >= band_stop)
passthrough = _passthrough_flatness(sweeps, kamil)
return {
"sweeps": n,
"slope_hz_per_rad": slope,
"axis_2pi_shift_hz": _TWO_PI * slope,
"band_start_hz": band_start,
"band_stop_hz": band_stop,
"anchor_phase": anchor_phase,
"anchor_unwrapped": anchor_unwrapped,
"branch_turns": branch_turns,
"n_branch_wraps": int(wrap_indices.size),
"wrap_indices": wrap_indices,
"span_rad": span_rad,
"n_points": n_points,
"f_lo": f_lo,
"f_hi": f_hi,
"f_lo_corr": f_lo_corr,
"f_hi_corr": f_hi_corr,
"axis_shift_hz": axis_shift,
"ref_amp_med": ref_amp_med,
"coverage_raw": float(np.mean(covers_raw)),
"coverage_corr": float(np.mean(covers_corr)),
# Largest sweep-to-sweep step of the *unwrapped* anchor: the fix is only
# valid if this stays below pi (otherwise the cross-sweep unwrap is itself
# ambiguous). Reported so we can confirm the float really is slow.
"max_unwrapped_anchor_step_rad": float(np.max(np.abs(np.diff(anchor_unwrapped))))
if n > 1
else 0.0,
**passthrough,
}
def _passthrough_flatness(
sweeps: list[tuple[np.ndarray, np.ndarray, np.ndarray]], kamil
) -> dict:
"""Run the live processor and measure how flat each pass-through trace is.
Builds the same :class:`KamilAdcSweepProcessor` the runtime uses, processes
every sweep onto the fixed grid, then measures each trace's deviation from a
robust (median) reference trace. A branch-wrapped sweep lands on a frequency
axis offset by ~one 2*pi turn, so after resampling it beats against the
reference — showing up as a large phase-residual std. This is the observable
symptom ("constant offset + oscillations") tied back to the raw-phase analysis.
"""
processor = KamilAdcSweepProcessor(KamilAdcProcessingParams.from_kamil_model(kamil))
traces = [processor.process(main, reference) for _steps, main, reference in sweeps]
passed = np.array([trace is not None for trace in traces])
n = len(traces)
phase_residual_std = np.full(n, np.nan)
amp_residual_std = np.full(n, np.nan)
stack = np.array([trace for trace in traces if trace is not None], dtype=np.complex128)
if stack.shape[0] >= 3:
# Robust per-grid-point reference: median of real/imag over passing sweeps.
reference_trace = np.median(stack.real, axis=0) + 1j * np.median(stack.imag, axis=0)
safe_reference = np.where(np.abs(reference_trace) < 1e-12, 1.0 + 0j, reference_trace)
passed_index = 0
for i, trace in enumerate(traces):
if trace is None:
continue
ratio = trace.astype(np.complex128) / safe_reference
phase_residual_std[i] = float(np.std(np.angle(ratio)))
amp_residual_std[i] = float(np.std(np.abs(ratio)))
passed_index += 1
return {
"passthrough_passed_fraction": float(np.mean(passed)),
"passthrough_phase_residual_std": phase_residual_std,
"passthrough_amp_residual_std": amp_residual_std,
}
def _report_diagnosis(result: dict) -> None:
"""Print the detailed reference-phase / frequency-axis diagnosis."""
n = result["sweeps"]
print("\n" + "=" * 72)
print(f"Kamil ADC reference-phase DIAGNOSIS over {n} sweeps")
print("=" * 72)
print(
f"calibration slope = {result['slope_hz_per_rad'] / 1e6:.4f} MHz/rad\n"
f"one 2*pi branch wrap = {result['axis_2pi_shift_hz'] / 1e6:.2f} MHz axis shift"
)
print("-" * 72)
print(
f"anchor phase angle(ref[0]): {_summarize(result['anchor_phase'])} rad\n"
f"phase span (stop - start) : {_summarize(result['span_rad'])} rad "
f"(std {np.std(result['span_rad']):.4f} -> fixed-delay {'OK' if np.std(result['span_rad']) < 1.0 else 'SUSPECT'})\n"
f"points per sweep : {_summarize(result['n_points'].astype(float))}"
)
print("-" * 72)
print(
f"branch wraps detected : {result['n_branch_wraps']} "
f"(at sweep indices {result['wrap_indices'].tolist()})"
)
print(f"branch turns range : {_summarize(result['branch_turns'].astype(float))} turns")
print(
f"max sweep-to-sweep step of UNWRAPPED anchor = "
f"{result['max_unwrapped_anchor_step_rad']:.4f} rad "
f"({'< pi: cross-sweep unwrap is unambiguous' if result['max_unwrapped_anchor_step_rad'] < np.pi else '>= pi: AMBIGUOUS, fix may misstep'})"
)
print("-" * 72)
print("reconstructed band edges (GHz), config calibration:")
print(f" raw start {_summarize(result['f_lo'] / 1e9)} stop {_summarize(result['f_hi'] / 1e9)}")
print(f" corrected start {_summarize(result['f_lo_corr'] / 1e9)} stop {_summarize(result['f_hi_corr'] / 1e9)}")
print(
f"band [{result['band_start_hz'] / 1e9:.3f}, {result['band_stop_hz'] / 1e9:.3f}] GHz covered: "
f"raw {result['coverage_raw'] * 100:.1f}% -> corrected {result['coverage_corr'] * 100:.1f}%"
)
print("-" * 72)
phase_res = result["passthrough_phase_residual_std"]
finite = phase_res[np.isfinite(phase_res)]
print(f"pass-through: {result['passthrough_passed_fraction'] * 100:.1f}% of sweeps covered the band")
if finite.size:
print(f" phase-residual std vs median trace: {_summarize(finite)} rad")
# Sweeps whose pass-through trace is far from flat (the visible glitches).
bad = np.nonzero(phase_res > (np.median(finite) + 5.0 * (np.std(finite) + 1e-9)))[0]
print(f" distorted sweeps (>5 sigma residual): {bad.tolist()}")
print(f" do they coincide with branch wraps? wraps at {result['wrap_indices'].tolist()}")
print("-" * 72)
print(
"HOW TO READ THIS:\n"
" * phase span std small -> fixed reference delay confirmed (model holds).\n"
" * 'branch wraps' are the suspected glitches: angle(ref[0]) crossing +/-pi.\n"
" * If the UNWRAPPED-anchor max step stays < pi, the sweep-to-sweep float is\n"
" slow enough that a CROSS-SWEEP unwrap can separate the real drift from the\n"
" spurious 2*pi wrap. That is the candidate fix.\n"
" * Expect: 'corrected' band edges move CONTINUOUSLY (no ~157 MHz steps) while\n"
" 'raw' jumps by one wrap at each detected wrap index, and the distorted\n"
" pass-through sweeps line up with those wrap indices."
)
print("=" * 72 + "\n")
def _dump_diagnosis(
path: Path,
result: dict,
sweeps: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
) -> None:
"""Save per-sweep metrics and the full raw sweeps to an ``.npz`` for offline study."""
arrays = {key: value for key, value in result.items() if isinstance(value, np.ndarray)}
arrays["raw_steps"] = np.array([steps for steps, _m, _r in sweeps], dtype=object)
arrays["raw_main"] = np.array([main for _s, main, _r in sweeps], dtype=object)
arrays["raw_reference"] = np.array([reference for _s, _m, reference in sweeps], dtype=object)
np.savez(path, **arrays)
print(f"Saved diagnosis data ({len(sweeps)} sweeps) to {path}")
def _plot_diagnosis(path: Path, result: dict) -> None:
"""Render the key diagnostic plots to a PNG (no-op if matplotlib is missing)."""
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except Exception as exc: # noqa: BLE001 — plotting is optional
logger.warning("Plotting skipped (matplotlib unavailable: %s)", exc)
return
sweep_index = np.arange(result["sweeps"])
fig, axes = plt.subplots(3, 1, figsize=(11, 12), sharex=True)
axes[0].plot(sweep_index, result["anchor_phase"], ".", label="angle(ref[0]) raw (wrapped)")
axes[0].plot(sweep_index, result["anchor_unwrapped"], "-", label="unwrapped across sweeps (fix)")
for wrap in result["wrap_indices"]:
axes[0].axvline(wrap, color="r", alpha=0.3)
axes[0].set_ylabel("anchor phase [rad]")
axes[0].legend(loc="best")
axes[0].set_title("Reference anchor phase: raw wraps vs cross-sweep unwrap")
axes[1].plot(sweep_index, result["f_lo"] / 1e9, ".", label="start raw")
axes[1].plot(sweep_index, result["f_lo_corr"] / 1e9, "-", label="start corrected")
axes[1].plot(sweep_index, result["f_hi"] / 1e9, ".", label="stop raw")
axes[1].plot(sweep_index, result["f_hi_corr"] / 1e9, "-", label="stop corrected")
axes[1].axhline(result["band_start_hz"] / 1e9, color="k", ls="--", alpha=0.5)
axes[1].axhline(result["band_stop_hz"] / 1e9, color="k", ls="--", alpha=0.5)
axes[1].set_ylabel("reconstructed band edge [GHz]")
axes[1].legend(loc="best")
axes[2].plot(sweep_index, result["passthrough_phase_residual_std"], ".", label="phase residual std")
for wrap in result["wrap_indices"]:
axes[2].axvline(wrap, color="r", alpha=0.3, label="_branch wrap")
axes[2].set_ylabel("pass-through residual [rad]")
axes[2].set_xlabel("sweep index")
axes[2].legend(loc="best")
fig.tight_layout()
fig.savefig(path, dpi=110)
plt.close(fig)
print(f"Saved diagnosis plot to {path}")
def _apply(config_path: Path, result: dict) -> None:
"""Write the migrated collector args + calibrated anchors + points to the config."""
payload = json.loads(config_path.read_text(encoding="utf-8"))
@@ -228,6 +519,17 @@ def main() -> int:
parser.add_argument("--warmup", type=int, default=10, help="Sweeps to discard first (default 10)")
parser.add_argument("--apply", action="store_true", help="Write the calibration back to --config")
parser.add_argument("--no-laser", action="store_true", help="Skip laser setup (already running)")
parser.add_argument(
"--diagnose",
action="store_true",
help="Investigate reference-phase branch wraps instead of calibrating",
)
parser.add_argument(
"--dump", type=Path, default=None, help="(--diagnose) save raw sweeps + metrics to this .npz"
)
parser.add_argument(
"--plot", type=Path, default=None, help="(--diagnose) render diagnostic plots to this .png"
)
args = parser.parse_args()
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(name)s: %(message)s")
@@ -246,6 +548,22 @@ def main() -> int:
service = KamilAdcService(config)
logger.info("Opening collector: %s", " ".join(service.command))
_open_with_retry(service)
if args.diagnose:
try:
sweeps = _capture_raw_sweeps(service, warmup=args.warmup, sweeps=args.sweeps)
finally:
service.close()
if len(sweeps) < max(2, args.sweeps // 2):
raise SystemExit(f"Only {len(sweeps)} usable sweeps captured; check the reference signal")
result = _diagnose(sweeps, config)
_report_diagnosis(result)
if args.dump is not None:
_dump_diagnosis(args.dump, result, sweeps)
if args.plot is not None:
_plot_diagnosis(args.plot, result)
return 0
try:
phases = _capture_reference_phases(service, warmup=args.warmup, sweeps=args.sweeps)
finally: