"""Signal processing for the Kamil ADC reference-channel acquisition. Each sweep arrives as two aligned complex arrays — the *main* signal and a *reference* signal — sampled at the same step indices (see :mod:`python_app.hardware_full.kamil_adc.protocol`). Turning that into a stable, comparable S21 trace is a fixed three-stage pipeline: 1. **Frequency axis from the reference phase.** The reference arm has a constant electrical delay, so its unwrapped phase is an affine function of frequency. Two fixed calibration anchors ``(phase0, freq0)`` and ``(phase1, freq1)`` — supplied by config, *never* derived from the live sweep — define that law:: 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``). 2. **Amplitude normalization.** ``S = main / |reference|`` divides out the stimulus amplitude. Only the magnitude is removed; the reference phase is used solely for the axis above. 3. **Crop + resample onto a fixed grid.** The floating per-sweep axis is resampled (linear, on real and imaginary parts) onto a single hardcoded grid ``linspace(band_start, band_stop, band_points)``. Every sweep then shares a byte-identical frequency axis, so traces can be averaged and subtracted. A sweep whose floated range does not fully span the band is *rejected* rather than edge-extrapolated. Everything here is pure and free of I/O so it can be unit-tested in isolation. """ from __future__ import annotations from dataclasses import dataclass import numpy as np from python_app.models.run_config_schema import KamilAdcModel # A reference sample whose magnitude falls to (effectively) zero carries no usable # amplitude or phase; such points are dropped before normalization rather than # producing a division blow-up. The threshold only guards genuine zeros — real # reference frames are emitted from settled measurements and sit far above it. _REFERENCE_AMPLITUDE_FLOOR = 1e-9 # Linear interpolation onto the band needs at least two distinct source samples; # 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: """Immutable phase→frequency calibration and output-grid definition. Built from :class:`~python_app.models.run_config_schema.KamilAdcModel` via :meth:`from_kamil_model`, but kept decoupled from the config types so the processing can be exercised with plain numbers in tests. """ phase0_rad: float freq0_hz: float phase1_rad: float freq1_hz: float band_start_hz: float band_stop_hz: float band_points: int reference_amplitude_floor: float = _REFERENCE_AMPLITUDE_FLOOR def __post_init__(self) -> None: if not np.isfinite([self.phase0_rad, self.phase1_rad, self.freq0_hz, self.freq1_hz]).all(): raise ValueError("Kamil ADC phase calibration anchors must be finite") if self.phase0_rad == self.phase1_rad: raise ValueError("Kamil ADC phase calibration anchors must use distinct phases") if self.freq0_hz == self.freq1_hz: raise ValueError("Kamil ADC phase calibration anchors must map to distinct frequencies") if not (np.isfinite(self.band_start_hz) and np.isfinite(self.band_stop_hz)): raise ValueError("Kamil ADC band edges must be finite") if self.band_stop_hz <= self.band_start_hz: raise ValueError("Kamil ADC band stop_hz must be greater than start_hz") if self.band_points < 2: raise ValueError("Kamil ADC band points must be >= 2") if self.reference_amplitude_floor <= 0.0: raise ValueError("Kamil ADC reference amplitude floor must be > 0") @classmethod def from_kamil_model(cls, kamil_adc: KamilAdcModel) -> KamilAdcProcessingParams: """Build parameters from the ``radar.kamil_adc`` config section.""" calibration = kamil_adc.phase_calibration band = kamil_adc.band return cls( phase0_rad=float(calibration.phase0_rad), freq0_hz=float(calibration.freq0_hz), phase1_rad=float(calibration.phase1_rad), freq1_hz=float(calibration.freq1_hz), band_start_hz=float(band.start_hz), band_stop_hz=float(band.stop_hz), band_points=int(band.points), ) @property def hz_per_rad(self) -> float: """Frequency change per radian of reference phase (the calibration slope).""" return (self.freq1_hz - self.freq0_hz) / (self.phase1_rad - self.phase0_rad) class KamilAdcSweepProcessor: """Turns aligned (main, reference) sweeps into S21 traces on a fixed grid. The output grid is computed once from the parameters and reused for every sweep, so all traces this processor emits share one identical frequency axis. """ __slots__ = ("_params", "_grid_hz") def __init__(self, params: KamilAdcProcessingParams) -> None: self._params = params self._grid_hz = np.linspace( params.band_start_hz, params.band_stop_hz, params.band_points, dtype=np.float64 ) @property def params(self) -> KamilAdcProcessingParams: return self._params @property def grid_hz(self) -> np.ndarray: """The fixed output frequency axis (float32), identical for every sweep.""" return self._grid_hz.astype(np.float32) def reference_frequency_axis(self, reference: np.ndarray) -> np.ndarray: """Map a reference signal's absolute unwrapped phase to frequency (Hz). 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))) ) 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. ``main`` and ``reference`` are equal-length complex arrays ordered by ascending step index. ``None`` is returned when the sweep is malformed (too few usable points) or does not fully cover the configured band. """ main = np.asarray(main, dtype=np.complex128) reference = np.asarray(reference, dtype=np.complex128) if main.size < _MIN_USABLE_POINTS or main.size != reference.size: return None # Frequency axis from the absolute unwrapped reference phase (step order). freqs = self.reference_frequency_axis(reference) # Amplitude-only normalization; drop points where the reference vanished. reference_amplitude = np.abs(reference) with np.errstate(divide="ignore", invalid="ignore"): s21 = main / reference_amplitude usable = ( (reference_amplitude > self._params.reference_amplitude_floor) & np.isfinite(freqs) & np.isfinite(s21.real) & np.isfinite(s21.imag) ) if int(np.count_nonzero(usable)) < _MIN_USABLE_POINTS: return None freqs = freqs[usable] s21 = s21[usable] # Sort onto a monotonically increasing axis (handles either sweep # direction) so interpolation and the coverage check are well defined. order = np.argsort(freqs, kind="stable") freqs = freqs[order] s21 = s21[order] # Reject sweeps that do not span the whole band: interpolating past the # measured edge would inject flat, non-physical points. if freqs[0] > self._params.band_start_hz or freqs[-1] < self._params.band_stop_hz: return None # Linear interpolation on real/imaginary parts. Because the band lies # within [freqs[0], freqs[-1]], np.interp never extrapolates here. real = np.interp(self._grid_hz, freqs, s21.real) imag = np.interp(self._grid_hz, freqs, s21.imag) return (real + 1j * imag).astype(np.complex64)