diff --git a/python_app/hardware_full/kamil_adc/processing.py b/python_app/hardware_full/kamil_adc/processing.py index 2180fe4..3cf0023 100644 --- a/python_app/hardware_full/kamil_adc/processing.py +++ b/python_app/hardware_full/kamil_adc/processing.py @@ -13,7 +13,10 @@ 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. + 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``). 2. **Amplitude normalization.** ``S = main / |reference|`` divides out the stimulus amplitude. Only the magnitude is removed; the reference phase is used @@ -47,6 +50,13 @@ _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: @@ -133,9 +143,45 @@ class KamilAdcSweepProcessor: Returns frequencies in *step order* (not sorted); see the module docstring for the calibration law. """ - phase = np.unwrap(np.angle(np.asarray(reference))) + phase = self._anchor_phase_to_calibration_branch( + 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/hardware_full/librevna_multi_device_driver/cycle_collection.py b/python_app/hardware_full/librevna_multi_device_driver/cycle_collection.py index 91f9db8..d887768 100644 --- a/python_app/hardware_full/librevna_multi_device_driver/cycle_collection.py +++ b/python_app/hardware_full/librevna_multi_device_driver/cycle_collection.py @@ -53,6 +53,32 @@ _MIN_CYCLE_ALIGNMENT_TOLERANCE_SECONDS = 0.5e-3 # S-parameter into noise. Mirrors the single-device C++ driver's reference guard. _MIN_REFERENCE_MAGNITUDE = 1e-12 +# A dropped USB frame leaves a single sweep point undelivered. Instead of +# abandoning the whole cycle and re-arming the hardware (a multi-second stall), +# we keep every delivered point and emit the few missing ones as zero. This caps +# that tolerance: once more than this fraction of a trace's points are +# unmeasurable, the loss is no longer a stray dropped frame but a genuine fault, +# so the cycle is rejected for proper recovery rather than returned as a +# mostly-zero trace. +_MAX_ZERO_FILLED_POINTS_FRACTION = 0.05 + + +class _DeviceCollectionState: + """Mutable per-device flag marking that this device's capture has finished. + + A device finishes when its final requested sweep cycle closes — either every + point of that cycle was delivered, or the free-running hardware sweep wrapped + back to point 0, proving the still-missing points are dropped frames that will + never arrive. The collector thread polls this flag so a gap ends the cycle at + the wrap instead of stalling until the no-progress timeout fires; the missing + points are zero-filled when the S-parameters are computed. + """ + + __slots__ = ("last_cycle_closed",) + + def __init__(self) -> None: + self.last_cycle_closed = False + class _CrossDeviceCycleSynchronizer: """Confirm every device thread anchors cycle 0 on the SAME physical sweep. @@ -158,7 +184,19 @@ def collect_complete_running_sweep_cycles( all_device_connections = [master_device_connection, *slave_device_connections] point_count = active_sweep_configuration.points - frequencies_hz = np.zeros(point_count, dtype=np.float64) + # Seed the frequency axis with the configured linear sweep grid so a point the + # master never delivers still carries its correct frequency instead of a 0 Hz + # hole. Every delivered point overwrites its slot with the device-reported + # frequency, which lands on this same grid. + if point_count > 1: + frequencies_hz = np.linspace( + float(active_sweep_configuration.start_hz), + float(active_sweep_configuration.stop_hz), + point_count, + dtype=np.float64, + ) + else: + frequencies_hz = np.array([float(active_sweep_configuration.start_hz)], dtype=np.float64) master_reference_measurements_by_port = { port: np.full((cycle_count, point_count), np.nan + 1j * np.nan, dtype=complex) for port in stimulus_ports @@ -254,14 +292,17 @@ def collect_complete_running_sweep_cycles( def collect_datapoints_from_device( device_connection: LibreVnaUsbBulkConnection, handle_datapoint: Callable[[ParsedVnaDatapoint], bool], + device_state: _DeviceCollectionState, ) -> None: - """Read datapoints from one device until enough are consumed or a timeout fires. + """Read datapoints from one device until its cycle closes or a timeout fires. Runs on a worker thread. Each datapoint is offered to ``handle_datapoint``, which returns whether it was consumed; only consumed datapoints count toward - progress and refresh the no-progress timeout. On any timeout or transport - error the error is recorded and ``stop_collection_requested`` is set so the - other collector threads also stop. + progress and refresh the no-progress timeout. The loop ends when every + requested point has arrived, when ``device_state.last_cycle_closed`` is set + (the sweep wrapped past the final cycle, so any remaining points are dropped + frames), or when a timeout/transport error is recorded — in which case + ``stop_collection_requested`` is set so the other collector threads stop too. """ datapoints_received = 0 expected_datapoint_count = cycle_count * point_count @@ -278,6 +319,13 @@ def collect_complete_running_sweep_cycles( if stop_collection_requested.is_set(): return + if device_state.last_cycle_closed: + # The sweep wrapped past the final requested cycle: its window is + # closed, so any points still missing are dropped frames, not points + # in flight. Finish now and let them be zero-filled downstream + # instead of waiting out the no-progress timeout. + return + now = time.monotonic() remaining_timeout_seconds = (last_consumed_timestamp + datapoint_timeout_seconds) - now if remaining_timeout_seconds <= 0: @@ -392,6 +440,7 @@ def collect_complete_running_sweep_cycles( def build_cycle_tracking_handler( cycle_aware_handler: Callable[[ParsedVnaDatapoint, int], None], + device_state: _DeviceCollectionState, ) -> Callable[[ParsedVnaDatapoint], bool]: """Wrap a cycle-aware handler with cross-device cycle tracking. @@ -446,6 +495,13 @@ def collect_complete_running_sweep_cycles( cycle_tracking_state["current_cycle_index"] += 1 current_cycle_index = cycle_tracking_state["current_cycle_index"] if current_cycle_index >= cycle_count: + # The sweep just wrapped past the final requested cycle, closing its + # collection window. Any of its points we never received are now + # confirmed dropped frames, not points still in flight, so flag the + # device done and let the collector thread stop immediately. The + # missing points are emitted as zero when the S-parameters are + # computed; delivered points keep their real values. + device_state.last_cycle_closed = True return False cycle_aware_handler(parsed_datapoint, current_cycle_index) @@ -498,7 +554,9 @@ def collect_complete_running_sweep_cycles( point_index, ] = port_receiver_value - def build_slave_datapoint_handler(slave_index: int) -> Callable[[ParsedVnaDatapoint], bool]: + def build_slave_datapoint_handler( + slave_index: int, device_state: _DeviceCollectionState + ) -> Callable[[ParsedVnaDatapoint], bool]: """Build a cycle-tracking datapoint handler for the given slave device. The slave's two receivers map to ports ``2*slave_index + 3`` and ``+ 4``, @@ -529,21 +587,31 @@ def collect_complete_running_sweep_cycles( point_index, ] = port_receiver_value - return build_cycle_tracking_handler(handle_slave_datapoint) + return build_cycle_tracking_handler(handle_slave_datapoint, device_state) + master_device_state = _DeviceCollectionState() collection_threads = [ threading.Thread( target=collect_datapoints_from_device, - args=(master_device_connection, build_cycle_tracking_handler(handle_master_datapoint)), + args=( + master_device_connection, + build_cycle_tracking_handler(handle_master_datapoint, master_device_state), + master_device_state, + ), daemon=True, name="collect-master", ) ] for slave_index, slave_device_connection in enumerate(slave_device_connections): + slave_device_state = _DeviceCollectionState() collection_threads.append( threading.Thread( target=collect_datapoints_from_device, - args=(slave_device_connection, build_slave_datapoint_handler(slave_index)), + args=( + slave_device_connection, + build_slave_datapoint_handler(slave_index, slave_device_state), + slave_device_state, + ), daemon=True, name=f"collect-slave{slave_index}", ) @@ -660,30 +728,61 @@ def calculate_last_cycle_s_parameters( raw_receiver_measurements_by_s_parameter: dict[str, np.ndarray], master_reference_measurements_by_port: dict[int, np.ndarray], ) -> dict[str, np.ndarray]: - """Convert raw receiver captures from the final sweep cycle into S-parameters.""" + """Convert raw receiver captures from the final sweep cycle into S-parameters. + + Points whose value never arrived (a dropped USB frame, still ``NaN``) or whose + incident reference is unusable (missing or a corrupt ~zero magnitude) cannot be + normalized, so they are emitted as exactly ``0``. Every delivered point keeps + its true normalized value, so no extra zeros are introduced. If more than + ``_MAX_ZERO_FILLED_POINTS_FRACTION`` of a trace is unmeasurable the loss is a + genuine fault rather than a stray dropped frame, and the cycle is rejected as + transient so the caller recovers instead of returning a mostly-zero trace. + """ s_parameters: dict[str, np.ndarray] = {} + zero_filled_points_by_trace: dict[str, int] = {} + point_count = 0 for s_parameter_name, raw_receiver_measurements in raw_receiver_measurements_by_s_parameter.items(): master_stimulus_port = int(s_parameter_name[-1]) master_reference_measurements = master_reference_measurements_by_port[master_stimulus_port] - if np.isnan(master_reference_measurements).any(): - missing_reference_count = int(np.isnan(master_reference_measurements).sum()) - raise TransientCollectionError( - f"Master port {master_stimulus_port} reference missing for {missing_reference_count} datapoints" - ) - if np.isnan(raw_receiver_measurements).any(): - missing_measurement_count = int(np.isnan(raw_receiver_measurements).sum()) - raise TransientCollectionError( - f"Measurement {s_parameter_name} missing for {missing_measurement_count} datapoints" - ) last_cycle_reference = master_reference_measurements[-1] - near_zero_reference_mask = np.abs(last_cycle_reference) <= _MIN_REFERENCE_MAGNITUDE - if near_zero_reference_mask.any(): - # A ~zero incident reference is a corrupt point: dividing by it would - # explode this S-parameter into noise. Reject the cycle for retry. + last_cycle_measurement = raw_receiver_measurements[-1] + point_count = last_cycle_measurement.size + + # Without a valid incident reference no trace at that point can be + # normalized; without this trace's receiver value there is nothing to + # normalize. Either way the point is unmeasurable and emitted as zero. + reference_unusable_mask = np.isnan(last_cycle_reference) | ( + np.abs(last_cycle_reference) <= _MIN_REFERENCE_MAGNITUDE + ) + measurement_missing_mask = np.isnan(last_cycle_measurement) + zero_fill_mask = reference_unusable_mask | measurement_missing_mask + + zero_filled_count = int(zero_fill_mask.sum()) + if point_count and zero_filled_count > _MAX_ZERO_FILLED_POINTS_FRACTION * point_count: raise TransientCollectionError( - f"Master port {master_stimulus_port} reference is ~zero for " - f"{int(near_zero_reference_mask.sum())} datapoints (corrupt incident signal)" + f"Trace {s_parameter_name} is unmeasurable at {zero_filled_count}/{point_count} " + f"points (missing reference/measurement); rejecting cycle for recovery" ) - s_parameters[s_parameter_name.lower()] = raw_receiver_measurements[-1] / last_cycle_reference + + # Substitute neutral operands so the division never produces NaN/inf, then + # force the unmeasurable points to zero. Delivered points are untouched. + safe_reference = np.where(reference_unusable_mask, 1.0 + 0j, last_cycle_reference) + safe_measurement = np.where(measurement_missing_mask, 0.0 + 0j, last_cycle_measurement) + s_parameter_values = safe_measurement / safe_reference + s_parameter_values[zero_fill_mask] = 0.0 + + if zero_filled_count: + zero_filled_points_by_trace[s_parameter_name] = zero_filled_count + s_parameters[s_parameter_name.lower()] = s_parameter_values + + if zero_filled_points_by_trace: + # One concise line per affected cycle (a rare dropped-frame event), so the + # gap-fill is visible at INFO without enabling DEBUG: this is the positive + # signal that the cycle was kept instead of triggering a re-arm stall. + logger.info( + "Zero-filled dropped sweep point(s) to keep the cycle (out of %d points/trace): %s", + point_count, + zero_filled_points_by_trace, + ) return s_parameters diff --git a/python_app/tests/test_kamil_adc_processing.py b/python_app/tests/test_kamil_adc_processing.py index 288885a..e010a50 100644 --- a/python_app/tests/test_kamil_adc_processing.py +++ b/python_app/tests/test_kamil_adc_processing.py @@ -135,6 +135,41 @@ 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. diff --git a/python_app/tests/test_librevna_multidevice.py b/python_app/tests/test_librevna_multidevice.py index cfc4055..6cd4645 100644 --- a/python_app/tests/test_librevna_multidevice.py +++ b/python_app/tests/test_librevna_multidevice.py @@ -225,7 +225,15 @@ class CycleAlignmentToleranceCalibrationTest(unittest.TestCase): class LastCycleSParameterTest(unittest.TestCase): - """Reference/measurement division + the corrupt-reference and NaN guards.""" + """Reference/measurement division, gap zero-filling, and the over-loss guard. + + A handful of dropped USB frames (points that never arrived, left as ``NaN``) + are emitted as zero so the cycle is kept instead of stalling on recovery, while + delivered points keep their true normalized value. Losing more than + ``_MAX_ZERO_FILLED_POINTS_FRACTION`` of a trace is treated as a genuine fault + and rejected as transient. The tiny-trace cases below cross that fraction with a + single missing point, so they still raise. + """ def test_divides_measurement_by_reference(self) -> None: reference = {1: np.array([[2 + 0j, 0 + 2j, 4 + 0j]], dtype=complex)} @@ -233,6 +241,49 @@ class LastCycleSParameterTest(unittest.TestCase): result = calculate_last_cycle_s_parameters(measurement, reference) np.testing.assert_allclose(result["s31"], np.array([2 + 0j, 2 + 0j, 0.5 + 0j])) + def test_single_dropped_measurement_is_zero_filled(self) -> None: + # One missing point in a realistically sized trace (< the loss fraction) is + # emitted as zero; every delivered point keeps its real value, no extras. + point_count = 200 + reference = {1: np.full((1, point_count), 2 + 0j, dtype=complex)} + measurement = {"S31": np.full((1, point_count), 4 + 0j, dtype=complex)} + measurement["S31"][0, 100] = np.nan + 1j * np.nan + result = calculate_last_cycle_s_parameters(measurement, reference) + self.assertEqual(result["s31"][100], 0) + np.testing.assert_allclose(np.delete(result["s31"], 100), 2 + 0j) + + def test_single_dropped_reference_is_zero_filled(self) -> None: + # A missing master reference cannot normalize that point for any trace, so + # the point is zeroed; the rest of the trace is untouched. + point_count = 200 + reference = {1: np.full((1, point_count), 2 + 0j, dtype=complex)} + reference[1][0, 50] = np.nan + 1j * np.nan + measurement = {"S31": np.full((1, point_count), 4 + 0j, dtype=complex)} + result = calculate_last_cycle_s_parameters(measurement, reference) + self.assertEqual(result["s31"][50], 0) + np.testing.assert_allclose(np.delete(result["s31"], 50), 2 + 0j) + + def test_single_near_zero_reference_is_zero_filled(self) -> None: + # A ~zero incident reference would explode the division into noise; one such + # point is zeroed rather than aborting the whole cycle, and stays finite. + point_count = 200 + reference = {1: np.full((1, point_count), 2 + 0j, dtype=complex)} + reference[1][0, 75] = 0 + 0j + measurement = {"S31": np.full((1, point_count), 4 + 0j, dtype=complex)} + result = calculate_last_cycle_s_parameters(measurement, reference) + self.assertEqual(result["s31"][75], 0) + self.assertTrue(np.isfinite(result["s31"]).all()) + + def test_excessive_loss_is_transient(self) -> None: + # Losing more than the tolerated fraction is a real fault, not a stray + # dropped frame: reject as transient so the caller recovers properly. + point_count = 200 + reference = {1: np.full((1, point_count), 2 + 0j, dtype=complex)} + measurement = {"S31": np.full((1, point_count), 4 + 0j, dtype=complex)} + measurement["S31"][0, :40] = np.nan + 1j * np.nan # 20% > 5% tolerance + with self.assertRaises(TransientCollectionError): + calculate_last_cycle_s_parameters(measurement, reference) + def test_nan_reference_is_transient(self) -> None: reference = {1: np.array([[2 + 0j, np.nan + 1j * np.nan, 4 + 0j]], dtype=complex)} measurement = {"S31": np.array([[4 + 0j, 4 + 0j, 2 + 0j]], dtype=complex)}