From da7c23f00fbc3d9d95fea93e87262a18f2106da5 Mon Sep 17 00:00:00 2001 From: ameoli Date: Mon, 26 Oct 2026 22:00:14 +0100 Subject: [PATCH] Validate fractional-lag recovery with independent acquisition noise --- docs/theory.md | 46 +++++++ tests/test_demodulation.py | 12 +- tests/test_fractional_lag_validation.py | 171 ++++++++++++++++++++++++ 3 files changed, 226 insertions(+), 3 deletions(-) create mode 100644 tests/test_fractional_lag_validation.py diff --git a/docs/theory.md b/docs/theory.md index a6f1c42..38c313d 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -448,6 +448,52 @@ or when optimization fails. This is a minimal validity check, not a confidence score. Tests use known phase ramps, both shift signs, periodic boundaries, and the fractional shift corresponding to a 5 K temperature change. +#### Fractional-lag validation (#105) + +The implementation uses fixed 8x coarse localization followed by continuous +refinement; it does not expose the configurable padding factor proposed in +the original issue. Validation therefore checks the shipped estimator against +known shifts, rather than treating a padding-factor sweep as an implemented +feature. + +`tests/test_fractional_lag_validation.py` checks 39 lags from -0.95 to +0.95 +bins and four cases straddling the periodic boundaries, for three Rayleigh +seeds and window sizes 64 and 488. The circular lag error must be below +`1e-7` bins. The legacy parabolic method is retained as a negative comparison: +its maximum error over the 39 interior lags exceeds 0.05 bins for each profile. +These are noiseless numerical checks, not physical resolution claims. + +The same file tests full-pipeline recovery of +/-17 microstrain and +/-5 K +on a 0.2--0.8 m segment of a 1 m fiber, for two seeds and both detector modes. +Reference sweep 0 and measurement sweep 1 share the fiber seed. Undoing the +known perturbation must recover the same complex fiber fingerprint. Subtracting +each matching noiseless receiver output gives nonzero receiver-noise residuals +whose absolute Pearson correlation must be below 0.03. Changing the global seed +would change the fiber too, so it is not a substitute for independent sweeps. + +This test uses a linear 40 nm sweep at 1550 nm, 131072 samples at 20 MS/s, +10 mW laser power, an 8 MHz fourth-order receiver filter, shot noise, +`1e-12 W/sqrt(Hz)` thermal NEP and `1e-9 A` dark current. The ADC has 18 bits, +a 0.2 V span and 1000 ohm transimpedance; clipping is checked explicitly. +For non-overlapping 1 cm gauges wholly inside the segment, the absolute median +bias must be below 5% of the imposed magnitude and the window RMSE below 10%. +The median offset before the segment must also be below 5% of that magnitude. +The existing strain robustness tests additionally use independent sweeps at +2 mm, 1 cm and 5 cm gauges, with shot noise on and off. + +The new tests record numerical errors and recovery/noise metrics as pytest +JUnit properties. For example: + +```bash +pytest tests/test_fractional_lag_validation.py -o junit_family=xunit1 --junitxml=validation.xml +``` + +These finite, reproducible cases cover the regression's receiver-noise regime. +They do not establish accuracy under fiber decorrelation, phase noise, RIN, +clipping, arbitrary SNR, or experimental instrument conditions. The original +issue's configurable-factor proposal is superseded by validation of the fixed +coarse search plus continuous refinement described above. + ### Quality metrics Spatial (`analysis/spatial_metrics.py`): peak SNR, integrated reflected power per region, baseline noise floor, dead-zone. Strain (`analysis/strain_quality.py`): RMS error vs ground truth, white and 1/f noise slopes, non-overlapping Allan variance. diff --git a/tests/test_demodulation.py b/tests/test_demodulation.py index 59a3f64..1b41ec8 100644 --- a/tests/test_demodulation.py +++ b/tests/test_demodulation.py @@ -604,7 +604,8 @@ def test_pipeline_strain_recovery(self): class TestStrainRecoveryRobustness: """Strain recovery via windowed sub-spectrum xcorr must hold across realistic detector / gauge combinations. Sweeps shot noise on/off - against three gauges (2 mm, 1 cm, 5 cm). Part of #5.""" + against three gauges (2 mm, 1 cm, 5 cm), with independent acquisition + noise on the same fiber. Part of #5 and #105.""" @staticmethod def _run_pair(shot_on): @@ -617,10 +618,15 @@ def _run_pair(shot_on): "thermal_nep": 0.0, "dark_current": 0.0}, "source": {**CFG["source"], "linewidth": 0.0}} - cfg_str = {**cfg_ref, "strain": {"segments": + cfg_str = {**cfg_ref, + "simulation": {**cfg_ref["simulation"], "n_sweeps": 2}, + "strain": {"segments": [{"start": 0.3, "end": 0.7, "epsilon": eps_true}]}} acq_ref = run_campaign(cfg_ref)[0] - acq_str = run_campaign(cfg_str)[0] + # Sweep index selects fresh detector/ADC/laser noise; the common base + # seed preserves the Rayleigh realization. Do not compare sweep 0 + # from both campaigns, which reuses the same noise streams. + acq_str = run_campaign(cfg_str)[1] return acq_ref, acq_str, eps_true @pytest.mark.parametrize("shot_on", [False, True]) diff --git a/tests/test_fractional_lag_validation.py b/tests/test_fractional_lag_validation.py new file mode 100644 index 0000000..414d99e --- /dev/null +++ b/tests/test_fractional_lag_validation.py @@ -0,0 +1,171 @@ +"""Validation of the shipped continuous estimator, including issue #105. + +Known phase ramps test numerical accuracy without a second peak estimator +as the oracle. Campaign pairs keep the fiber seed fixed and select different +sweeps so receiver noise cannot cancel between reference and measurement. +""" + +from copy import deepcopy + +import numpy as np +import pytest + +from helpers import CFG +from pyofdr.analysis.demodulation import ( + reflectogram_from_acq, + windowed_xcorr_strain, +) +from pyofdr.core.acquisition import Acquisition +from pyofdr.core.campaign import run_campaign +from pyofdr.detection.detector import Detector +from pyofdr.detection.filter import AntiAliasFilter +from pyofdr.utils.constants import C + + +@pytest.mark.parametrize("seed", [7, 42, 105]) +@pytest.mark.parametrize("window", [64, 488]) +def test_continuous_accuracy_across_fractional_lags(seed, window, record_property): + """Fixed 8x localization + continuous refinement recovers known lags. + + This validates the selected algorithm instead of exposing a padding-factor + option just to implement the alternative originally proposed in #105. + """ + rng = np.random.default_rng(seed) + reference = rng.normal(size=window) + 1j * rng.normal(size=window) + lags = np.r_[np.linspace(-0.95, 0.95, 39), + -window / 2 - 0.03, -window / 2 + 0.03, + window / 2 - 0.03, window / 2 + 0.03] + errors = {"continuous": [], "parabolic": []} + for lag in lags: + measured = reference * np.exp( + -2j * np.pi * lag * np.arange(window) / window) + for method in errors: + _, strain = windowed_xcorr_strain( + measured, reference, dz=1.0, gauge_length=window, + stride=window, sweep_range_hz=window, center_freq=1.0, + p_e=0.0, peak_method=method, + ) + recovered_lag = -strain[0] + assert -window / 2 <= recovered_lag < window / 2 + error = (recovered_lag - lag + window / 2) % window - window / 2 + errors[method].append(abs(error)) + + record_property("max_continuous_error_bins", max(errors["continuous"])) + record_property("max_parabolic_interior_error_bins", max(errors["parabolic"][:39])) + assert max(errors["continuous"]) < 1e-7 + # The regression must exercise the original bias, not only easy grid peaks. + assert max(errors["parabolic"][:39]) > 0.05 + + +def _campaign_config(seed, balanced): + cfg = deepcopy(CFG) + cfg["simulation"] = {"seed": seed, "n_sweeps": 2} + # 131072 samples: enough to resolve the 1 m fiber without Nyquist aliasing, + # but much cheaper than the default two million samples per sweep. + cfg["source"].update(sweep_duration=131072 / 20e6, + linewidth=0.0, rin_dB_per_Hz=None) + cfg["detection"].update( + balanced=balanced, bandwidth=8e6, filter_order=4, + shot_noise=True, thermal_nep=1e-12, dark_current=1e-9, + ) + cfg["adc"].update(sample_rate=20e6, bits=18, voltage_range=0.2, + input_impedance=1000.0) + cfg["strain"] = {"photoelastic_coefficient": 0.22, "segments": []} + cfg["temperature"] = {"thermal_expansion": 5.5e-7, + "thermo_optic": 6.5e-6, "segments": []} + return cfg + + +def _receiver_noise(acq, cfg): + """Subtract the matching noiseless receiver output, before ADC rounding.""" + quiet = deepcopy(cfg) + quiet["detection"].update(shot_noise=False, thermal_nep=0.0, + dark_current=0.0) + baseline = Acquisition(dt=acq.dt, photocurrent_main=acq.photocurrent_main) + baseline = Detector(quiet).process(baseline) + baseline = AntiAliasFilter(quiet).process(baseline) + residual = (acq.analog_main - baseline.analog_main)[0] + # Exclude filter settling and edges when checking noise correlation. + return residual[len(residual) // 10:-len(residual) // 10] + + +@pytest.mark.parametrize("seed", [7, 105]) +@pytest.mark.parametrize("balanced", [False, True]) +@pytest.mark.parametrize("kind,magnitude", [("strain", 17e-6), + ("temperature", 5.0)]) +@pytest.mark.parametrize("sign", [-1, 1]) +def test_recovery_with_independent_acquisition_noise(seed, balanced, + kind, magnitude, sign, + record_property): + cfg = _campaign_config(seed, balanced) + perturbed = deepcopy(cfg) + imposed = sign * magnitude + field = "epsilon" if kind == "strain" else "delta_T" + perturbed[kind]["segments"] = [{"start": 0.2, "end": 0.8, field: imposed}] + + # Same fiber seed, different sweep indices. Changing the global seed would + # also change the Rayleigh fingerprint and would test a different problem. + reference = run_campaign(cfg)[0] + measured = run_campaign(perturbed)[1] + assert (reference.sweep_index, measured.sweep_index) == (0, 1) + + p_e = cfg["strain"]["photoelastic_coefficient"] + thermal = (cfg["temperature"]["thermal_expansion"] + + cfg["temperature"]["thermo_optic"]) + coefficient = 1 - p_e if kind == "strain" else thermal + applied_field = (measured.strain_field if kind == "strain" + else measured.temperature_field) + wavelength = cfg["source"]["center_wavelength"] + phase = (4 * np.pi * cfg["fiber"]["n_core"] / wavelength + * coefficient * np.cumsum(applied_field) * measured.dz) + # Undo only the known physical perturbation: the complex fingerprint, + # not just its magnitude, must match the reference. + np.testing.assert_allclose( + measured.fiber_profile * np.exp(-1j * phase), + reference.fiber_profile, rtol=1e-12, atol=0, + ) + + noise_ref = _receiver_noise(reference, cfg) + noise_meas = _receiver_noise(measured, perturbed) + assert np.std(noise_ref) > 0 + assert np.std(noise_meas) > 0 + # The original same-sweep setup produces almost identical receiver noise; + # this assertion ensures it cannot silently return to these regressions. + noise_correlation = np.corrcoef(noise_ref, noise_meas)[0, 1] + record_property("receiver_noise_correlation", noise_correlation) + assert abs(noise_correlation) < 0.03 + + for acq in (reference, measured): + # Reject clipping as an explanation for estimator errors. + assert np.max(np.abs(acq.analog_main)) < cfg["adc"]["voltage_range"] / 2 + H_ref, _ = reflectogram_from_acq(reference) + H_meas, _ = reflectogram_from_acq(measured) + gauge = 0.01 + centers, apparent_strain = windowed_xcorr_strain( + H_meas, H_ref, reference.fft_dz, gauge_length=gauge, stride=gauge, + sweep_range_hz=reference.fft_sweep_range_hz, + center_freq=C / wavelength, p_e=p_e, + ) + recovered = (apparent_strain if kind == "strain" + else apparent_strain * (1 - p_e) / thermal) + # Non-overlapping windows wholly within the perturbed segment; exclude + # its edges where one gauge would average perturbed and quiet material. + inside = (centers > 0.2 + gauge) & (centers < 0.8 - gauge) + outside = (centers > 0.03) & (centers < 0.2 - gauge) + assert inside.sum() >= 40 + assert outside.sum() >= 10 + assert np.isfinite(recovered[inside | outside]).all() + errors = recovered[inside] - imposed + relative_bias = abs(np.median(errors)) / magnitude + relative_rmse = np.sqrt(np.mean(errors**2)) / magnitude + quiet_offset = abs(np.median(recovered[outside])) / magnitude + record_property("imposed_value", imposed) + record_property("recovered_median", np.median(recovered[inside])) + record_property("relative_median_bias", relative_bias) + record_property("relative_window_rmse", relative_rmse) + record_property("relative_quiet_offset", quiet_offset) + # Acceptance limits for this stated receiver/SNR regime, not universal + # uncertainty guarantees: median bias <5%, window RMSE <10% of input. + assert relative_bias < 0.05 + assert relative_rmse < 0.10 + assert quiet_offset < 0.05