From d7ddca6e767117e1b8340a0f95d8e8a6a61b0b95 Mon Sep 17 00:00:00 2001 From: Ayzen Date: Mon, 22 Jun 2026 16:30:55 +0300 Subject: [PATCH] added diagnostics --- .../hardware_full/kamil_adc/processing.py | 50 +-- python_app/scripts/kamil_adc_calibrate.py | 318 ++++++++++++++++++ python_app/tests/test_kamil_adc_processing.py | 35 -- 3 files changed, 320 insertions(+), 83 deletions(-) diff --git a/python_app/hardware_full/kamil_adc/processing.py b/python_app/hardware_full/kamil_adc/processing.py index 3cf0023..2180fe4 100644 --- a/python_app/hardware_full/kamil_adc/processing.py +++ b/python_app/hardware_full/kamil_adc/processing.py @@ -13,10 +13,7 @@ comparable S21 trace is a fixed three-stage pipeline: f(phase) = freq0 + (phase - phase0) * (freq1 - freq0) / (phase1 - phase0) Trigger jitter shifts every sample's absolute phase together, so the measured - band floats from sweep to sweep around the fixed calibration. When that float - carries the unwrap anchor (sample 0) across the +/-pi branch cut, a stray sweep - is offset by a whole 2*pi turn; it is snapped back onto the branch nearest the - calibration before mapping (see ``_anchor_phase_to_calibration_branch``). + band floats from sweep to sweep around the fixed calibration. 2. **Amplitude normalization.** ``S = main / |reference|`` divides out the stimulus amplitude. Only the magnitude is removed; the reference phase is used @@ -50,13 +47,6 @@ _REFERENCE_AMPLITUDE_FLOOR = 1e-9 # a sweep yielding fewer usable points is malformed and rejected. _MIN_USABLE_POINTS = 2 -# One full turn of phase. ``np.unwrap`` reconstructs each sweep's phase ramp but -# anchors it to the raw ``np.angle`` of the first sample, which lives on the -# (-pi, pi] branch. Trigger jitter occasionally lands that anchor on the far side -# of the +/-pi branch cut for a stray sweep or two, rigidly offsetting the whole -# ramp by exactly this much before it settles back onto the physical branch. -_PHASE_BRANCH_PERIOD_RAD = 2.0 * np.pi - @dataclass(frozen=True, slots=True) class KamilAdcProcessingParams: @@ -143,45 +133,9 @@ class KamilAdcSweepProcessor: Returns frequencies in *step order* (not sorted); see the module docstring for the calibration law. """ - phase = self._anchor_phase_to_calibration_branch( - np.unwrap(np.angle(np.asarray(reference))) - ) + phase = np.unwrap(np.angle(np.asarray(reference))) return self._params.freq0_hz + (phase - self._params.phase0_rad) * self._params.hz_per_rad - def _anchor_phase_to_calibration_branch(self, phase: np.ndarray) -> np.ndarray: - """Collapse a stray 2*pi branch excursion back onto the physical branch. - - ``np.unwrap`` reconstructs a continuous phase ramp but pins its absolute - level to the raw angle of the first sample, which lives on the (-pi, pi] - branch. Trigger jitter occasionally lands that anchor on the wrong side of - the +/-pi cut, rigidly shifting the whole sweep by one - :data:`_PHASE_BRANCH_PERIOD_RAD` (~157 MHz on the rig) until it settles back - a sweep or two later. Such an excursion would otherwise wreck the frequency - axis, the band-coverage check, and the normalization. - - The calibration's ``phase0_rad`` is the expected first-sample phase (its - median across many sweeps), so the physical branch is the one nearest it. - We round the first sample onto that branch and shift the whole ramp by the - same whole number of turns. This is: - - * **Stateless** — each sweep is judged only against the fixed calibration, - so a glitch can never propagate into, or latch, later sweeps. - * **Self-correcting** — a glitched sweep is pulled back onto the band and - yields usable data instead of being rejected. - * **Span-invariant** — it keys on the first sample (a fixed sweep start), - not on how much band the sweep happens to span. - - Genuine sweep-to-sweep float (well under pi against a calibration centered - on its median) rounds to zero turns and is left untouched. A float that - ever drifts past pi is a recalibration concern, not a per-sweep glitch. - """ - if phase.size == 0: - return phase - branch_turns = np.round((phase[0] - self._params.phase0_rad) / _PHASE_BRANCH_PERIOD_RAD) - if branch_turns: - phase = phase - branch_turns * _PHASE_BRANCH_PERIOD_RAD - return phase - def process(self, main: np.ndarray, reference: np.ndarray) -> np.ndarray | None: """Return the S21 trace resampled onto the fixed grid, or ``None`` to reject. diff --git a/python_app/scripts/kamil_adc_calibrate.py b/python_app/scripts/kamil_adc_calibrate.py index 0e799ea..5c3bcc4 100644 --- a/python_app/scripts/kamil_adc_calibrate.py +++ b/python_app/scripts/kamil_adc_calibrate.py @@ -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: diff --git a/python_app/tests/test_kamil_adc_processing.py b/python_app/tests/test_kamil_adc_processing.py index e010a50..288885a 100644 --- a/python_app/tests/test_kamil_adc_processing.py +++ b/python_app/tests/test_kamil_adc_processing.py @@ -135,41 +135,6 @@ class SweepProcessorTest(unittest.TestCase): np.testing.assert_allclose(result.real, expected, atol=1e-3) np.testing.assert_allclose(result.imag, 0.0, atol=1e-3) - def test_recovers_sweep_offset_by_a_full_branch(self) -> None: - """A stray +/-2*pi branch jump on the unwrap anchor must be snapped back. - - The same physical sweep, offset by one full turn (as happens when trigger - jitter carries sample 0 across the +/-pi cut), must map to the SAME frequency - axis as the unshifted sweep instead of sliding ~one branch off the band. - """ - processor = self._processor() - phase = np.linspace(0.0, 100.0, 401) # freq [2.0, 4.0] GHz, spans the band - baseline = processor.reference_frequency_axis(_reference(phase)) - for turns in (+1, -1, +2): - shifted = processor.reference_frequency_axis(_reference(phase + turns * 2.0 * np.pi)) - np.testing.assert_allclose(shifted, baseline, atol=1e-3) - - def test_branch_recovery_keeps_a_glitched_sweep_usable(self) -> None: - """A branch-jumped sweep is pulled back onto the band, not rejected.""" - processor = self._processor() - phase = np.linspace(0.0, 100.0, 401) + 2.0 * np.pi # one full turn off - ref = _reference(phase) - result = processor.process(np.abs(ref).astype(np.complex128), ref) - self.assertIsNotNone(result) - np.testing.assert_allclose(np.abs(result), 1.0, atol=1e-3) - - def test_leaves_genuine_sub_branch_float_untouched(self) -> None: - """A real None: # phase 0 -> 2 GHz, phase -100 -> 4 GHz (negative slope). Phase ramp # 0 -> -100 therefore sweeps frequency UP across the band.