Validate fractional-lag recovery with independent acquisition noise

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.

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])

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