Add optional full-band Galileo E1 CBOC support

2026-07-31 02:56:04 +03:00
parent 9db955d02e
commit 3662ba7e53
28 changed files with 2065 additions and 101 deletions
@@ -0,0 +1,160 @@
+"""Full-band CBOC acquisition/tracking evidence at 32.736 Msps."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from galileo_rx.acquisition.acquisition_manager import AcquisitionManager
+from galileo_rx.config import AcquisitionConfig, TrackingConfig
+from galileo_rx.signal_profiles import galileo_e1_profile
+from galileo_rx.signals.code_provider import SyntheticCodeProvider
+from galileo_rx.signals.e1_modulation import E1Modulation
+from galileo_rx.signals.fullband_e1_generator import FullBandE1SignalGenerator
+from galileo_rx.signals.signal_generator import ChannelModelConfig, SatelliteSignalConfig
+from galileo_rx.tracking.pilot_aided_e1 import synchronize_e1c_for_tracking
+from galileo_rx.tracking.tracking_channel import TrackingChannel
+
+
+FS_HZ = 32.736e6
+IF_HZ = 2.25e6
+DOPPLER_HZ = 1250.0
+DELAY_CHIPS = 80.25
+
+
+def _run_case(transmit_modulation: E1Modulation, receiver_modulation: E1Modulation):
+    profile = galileo_e1_profile()
+    provider = SyntheticCodeProvider()
+    generator = FullBandE1SignalGenerator(
+        profile,
+        provider,
+        FS_HZ,
+        IF_HZ,
+        satellites=[
+            SatelliteSignalConfig(
+                1,
+                "E1C",
+                code_delay_chips=DELAY_CHIPS,
+                doppler_hz=DOPPLER_HZ,
+                carrier_phase_rad=0.3,
+                cn0_db_hz=55.0,
+                modulation=transmit_modulation,
+            )
+        ],
+        channel=ChannelModelConfig(seed=2026),
+    )
+    samples, truth = generator.generate(round(0.16 * FS_HZ))
+    acquisition_config = AcquisitionConfig(
+        modulation=receiver_modulation,
+        doppler_min_hz=1000.0,
+        doppler_max_hz=1500.0,
+        doppler_step_hz=250.0,
+        coherent_integration_ms=4.0,
+        noncoherent_accumulations=2,
+        detection_threshold=2.0,
+        zero_padding_mode="none",
+        parabolic_interpolation=False,
+    )
+    acquisition, _ = AcquisitionManager(
+        profile, provider, FS_HZ, acquisition_config
+    ).acquire_one(samples, "E1C", 1, IF_HZ)
+    synchronization = synchronize_e1c_for_tracking(
+        samples,
+        profile,
+        provider,
+        FS_HZ,
+        acquisition,
+        round(0.008 * FS_HZ),
+        search_primary_epochs=30,
+        residual_frequency_step_hz=1.0,
+        modulation=receiver_modulation,
+    )
+    tracking = TrackingChannel(
+        profile,
+        provider,
+        FS_HZ,
+        TrackingConfig(
+            modulation=receiver_modulation,
+            pull_in_duration_ms=20.0,
+            loss_of_lock_epochs=20,
+        ),
+        synchronization.acquisition,
+        synchronization.tracking_start_sample,
+        synchronization.sync.phase,
+    )
+    measurements = tracking.process(samples[synchronization.tracking_start_sample :])
+    return acquisition, synchronization.sync, measurements, truth
+
+
+@pytest.mark.parametrize(
+    ("transmit_modulation", "receiver_modulation"),
+    [
+        (E1Modulation.BOC11, E1Modulation.BOC11),
+        (E1Modulation.CBOC, E1Modulation.BOC11),
+        (E1Modulation.CBOC, E1Modulation.CBOC),
+    ],
+)
+def test_fullband_modulation_matrix_acquires_and_tracks(
+    transmit_modulation: E1Modulation,
+    receiver_modulation: E1Modulation,
+) -> None:
+    acquisition, secondary, measurements, truth = _run_case(
+        transmit_modulation, receiver_modulation
+    )
+    expected = truth.satellites[0]
+    assert acquisition.detected
+    assert acquisition.raw_doppler_hz == DOPPLER_HZ
+    assert abs(acquisition.code_phase_chips - expected.initial_code_delay_chips) <= 1.0 / 32.0
+    assert secondary.normalized_correlation > 0.98
+    assert secondary.sign_errors == 0
+    assert len(measurements) >= 35
+    assert measurements[-1].carrier_lock
+    assert measurements[-1].code_lock
+    assert not any(item.cycle_slip for item in measurements)
+    assert not any(item.sample_discontinuity for item in measurements)
+    tail = measurements[-10:]
+    assert np.sqrt(np.mean([item.pll_error_rad**2 for item in tail])) < 0.12
+    assert np.sqrt(np.mean([item.dll_error_chips**2 for item in tail])) < 0.12
+
+
+@pytest.mark.parametrize("intermediate_frequency_hz", [0.0, 2.25e6])
+def test_cboc_pcps_carrier_wipeoff_at_zero_and_nonzero_if(
+    intermediate_frequency_hz: float,
+) -> None:
+    profile = galileo_e1_profile()
+    provider = SyntheticCodeProvider()
+    samples, _ = FullBandE1SignalGenerator(
+        profile,
+        provider,
+        FS_HZ,
+        intermediate_frequency_hz,
+        satellites=[
+            SatelliteSignalConfig(
+                5,
+                "E1B",
+                code_delay_chips=19.5,
+                doppler_hz=-750.0,
+                modulation=E1Modulation.CBOC,
+            )
+        ],
+        channel=ChannelModelConfig(add_awgn=False),
+    ).generate(round(0.008 * FS_HZ))
+    result, _ = AcquisitionManager(
+        profile,
+        provider,
+        FS_HZ,
+        AcquisitionConfig(
+            modulation=E1Modulation.CBOC,
+            doppler_min_hz=-1000.0,
+            doppler_max_hz=-500.0,
+            doppler_step_hz=250.0,
+            coherent_integration_ms=4.0,
+            noncoherent_accumulations=2,
+            detection_threshold=2.0,
+            zero_padding_mode="none",
+            parabolic_interpolation=False,
+        ),
+    ).acquire_one(samples, "E1B", 5, intermediate_frequency_hz)
+    assert result.detected
+    assert result.raw_doppler_hz == -750.0
+    assert abs(result.code_phase_chips - 20.5) <= 1.0 / 32.0
@@ -0,0 +1,108 @@
+"""Unit evidence for phase-continuous Galileo E1 CBOC replicas."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from galileo_rx.signals.boc import SquareSubcarrierNco, boc_subcarrier
+from galileo_rx.signals.e1_modulation import (
+    CBOC_ALPHA,
+    CBOC_BETA,
+    E1Modulation,
+    e1_subcarrier_replica,
+    e1_transmitted_subcarrier_replica,
+)
+
+
+def _phases(sample_count: int, sampling_frequency_hz: float) -> np.ndarray:
+    return np.arange(sample_count, dtype=np.float64) * 1.023e6 / sampling_frequency_hz
+
+
+def test_boc11_subcarrier_sequence() -> None:
+    wave, ending = boc_subcarrier(8, 8.0, 1.0)
+    np.testing.assert_array_equal(wave, [1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0])
+    assert ending == 0.0
+
+
+def test_boc61_subcarrier_sequence() -> None:
+    wave, ending = boc_subcarrier(24, 24.0, 6.0)
+    np.testing.assert_array_equal(wave, np.tile([1.0, 1.0, -1.0, -1.0], 6))
+    assert ending == 0.0
+
+
+def test_subcarrier_phase_continuity_across_arbitrary_chunks() -> None:
+    fs = 31.25e6  # deliberately nonintegral samples per BOC(6,1) sign interval
+    total = 10_003
+    reference = SquareSubcarrierNco(fs, 6.138e6, 0.137).samples(total)
+    chunked_nco = SquareSubcarrierNco(fs, 6.138e6, 0.137)
+    chunks = [1, 7, 31, 2, 997, 4093, 3, 4869]
+    chunked = np.concatenate([chunked_nco.samples(count) for count in chunks])
+    np.testing.assert_array_equal(chunked, reference)
+
+
+def test_cboc_power_normalization_and_component_sign() -> None:
+    phase = _phases(4092 * 32, 32.736e6)
+    sc11 = e1_subcarrier_replica("E1B", E1Modulation.BOC11, phase)
+    sc61 = np.where((np.floor(2.0 * np.mod(6.0 * phase, 1.0) + 1e-12).astype(int) & 1) == 0, 1.0, -1.0)
+    e1b = e1_subcarrier_replica("E1B", E1Modulation.CBOC, phase, 6.0 * phase)
+    e1c = e1_subcarrier_replica("E1C", E1Modulation.CBOC, phase, 6.0 * phase)
+    assert abs(float(np.mean(e1b**2)) - 1.0) < 1e-12
+    assert abs(float(np.mean(e1c**2)) - 1.0) < 1e-12
+    np.testing.assert_allclose((e1b + e1c) / 2.0, CBOC_ALPHA * sc11, atol=1e-15)
+    np.testing.assert_allclose((e1b - e1c) / 2.0, CBOC_BETA * sc61, atol=1e-15)
+    transmitted_e1b = e1_transmitted_subcarrier_replica(
+        "E1B", E1Modulation.CBOC, phase, 6.0 * phase
+    )
+    transmitted_e1c = e1_transmitted_subcarrier_replica(
+        "E1C", E1Modulation.CBOC, phase, 6.0 * phase
+    )
+    # Equation (12): low-rate terms oppose, high-rate terms have equal sign.
+    np.testing.assert_allclose((transmitted_e1b - transmitted_e1c) / 2.0, CBOC_ALPHA * sc11, atol=1e-15)
+    np.testing.assert_allclose((transmitted_e1b + transmitted_e1c) / 2.0, CBOC_BETA * sc61, atol=1e-15)
+
+
+def test_cboc_spectrum_contains_the_boc61_lobe() -> None:
+    fs = 32.736e6
+    phase = _phases(4092 * 32, fs)
+    boc11 = e1_subcarrier_replica("E1B", E1Modulation.BOC11, phase)
+    cboc = e1_subcarrier_replica("E1B", E1Modulation.CBOC, phase, 6.0 * phase)
+    frequencies = np.fft.rfftfreq(phase.size, 1.0 / fs)
+    boc_psd = np.abs(np.fft.rfft(boc11)) ** 2
+    cboc_psd = np.abs(np.fft.rfft(cboc)) ** 2
+    high_lobe = (frequencies > 5.8e6) & (frequencies < 6.5e6)
+    fundamental = (frequencies > 0.8e6) & (frequencies < 1.25e6)
+    assert cboc_psd[high_lobe].sum() > 1_000.0 * boc_psd[high_lobe].sum()
+    assert cboc_psd[fundamental].sum() > cboc_psd[high_lobe].sum()
+
+
+def test_cboc_autocorrelation_side_structure_and_epl_symmetry() -> None:
+    fs = 32.736e6
+    samples_per_chip = 32
+    phase = _phases(4092 * samples_per_chip, fs)
+    cboc = e1_subcarrier_replica("E1B", E1Modulation.CBOC, phase, 6.0 * phase)
+
+    def circular_correlation(offset_samples: int) -> float:
+        return float(np.mean(cboc * np.roll(cboc, offset_samples)))
+
+    assert circular_correlation(0) > 0.999999
+    assert circular_correlation(samples_per_chip // 2) < -0.80
+    for offset in (2, 4, 6, 8, 10, 12):
+        assert abs(circular_correlation(offset) - circular_correlation(-offset)) < 2e-4
+    # E and L at +/-0.25 chip see equal magnitude at the correct prompt.
+    assert abs(abs(circular_correlation(8)) - abs(circular_correlation(-8))) < 2e-4
+
+
+def test_cboc_sample_rate_independence() -> None:
+    powers = []
+    quarter_chip_correlations = []
+    for fs in (32.736e6, 30.0e6, 27.5e6):
+        count = int(round(4092 * fs / 1.023e6))
+        phase = _phases(count, fs)
+        cboc = e1_subcarrier_replica("E1C", E1Modulation.CBOC, phase, 6.0 * phase)
+        powers.append(float(np.mean(cboc**2)))
+        quarter_late = e1_subcarrier_replica(
+            "E1C", E1Modulation.CBOC, phase + 0.25, 6.0 * (phase + 0.25)
+        )
+        quarter_chip_correlations.append(float(np.mean(cboc * quarter_late)))
+    assert max(abs(power - 1.0) for power in powers) < 2e-3
+    assert np.ptp(quarter_chip_correlations) < 0.04
@@ -0,0 +1,72 @@
+from __future__ import annotations
+
+import numpy as np
+
+from galileo_rx.signal_profiles import galileo_e1_profile
+from galileo_rx.signals.code_provider import SyntheticCodeProvider
+from galileo_rx.signals.e1_modulation import E1Modulation
+from galileo_rx.signals.fullband_e1_generator import FullBandE1SignalGenerator
+from galileo_rx.signals.signal_generator import ChannelModelConfig, SatelliteSignalConfig
+
+
+def _generator() -> FullBandE1SignalGenerator:
+    return FullBandE1SignalGenerator(
+        galileo_e1_profile(),
+        SyntheticCodeProvider(),
+        32.736e6,
+        2.25e6,
+        satellites=[
+            SatelliteSignalConfig(
+                3,
+                "E1C",
+                code_delay_chips=7.25,
+                doppler_hz=1375.0,
+                modulation=E1Modulation.CBOC,
+            )
+        ],
+        channel=ChannelModelConfig(add_awgn=False),
+    )
+
+
+def test_fullband_generator_is_streaming_and_double_precision() -> None:
+    one_shot, truth = _generator().generate(4099)
+    chunked_generator = _generator()
+    chunked = np.concatenate([chunked_generator.generate(count)[0] for count in (1, 17, 1000, 3, 3078)])
+    assert one_shot.dtype == np.complex128
+    np.testing.assert_allclose(chunked, one_shot, atol=2e-13, rtol=0.0)
+    assert truth.satellites[0].initial_code_delay_samples == 7.25 * 32.0 + 32.0
+
+
+def test_fullband_generator_supports_real_if_quantization() -> None:
+    generator = FullBandE1SignalGenerator(
+        galileo_e1_profile(),
+        SyntheticCodeProvider(),
+        32.736e6,
+        4.092e6,
+        True,
+        [SatelliteSignalConfig(1, "E1B", modulation=E1Modulation.CBOC)],
+        ChannelModelConfig(add_awgn=False, quantization_bits=8),
+    )
+    samples, _ = generator.generate(2048)
+    assert samples.dtype == np.float64
+    assert np.max(samples) <= 1.0 - 1.0 / 128.0
+    np.testing.assert_allclose(samples * 128.0, np.round(samples * 128.0), atol=1e-12)
+
+
+def test_fullband_generator_exposes_independent_code_doppler_truth() -> None:
+    generator = FullBandE1SignalGenerator(
+        galileo_e1_profile(),
+        SyntheticCodeProvider(),
+        satellites=[
+            SatelliteSignalConfig(
+                1,
+                "E1C",
+                doppler_hz=1000.0,
+                code_doppler_hz=0.75,
+                modulation=E1Modulation.CBOC,
+            )
+        ],
+        channel=ChannelModelConfig(add_awgn=False),
+    )
+    _, truth = generator.generate(128)
+    assert truth.satellites[0].initial_code_doppler_hz == 0.75