Initial Galileo receiver and SignalSim integration
This commit is contained in:
@@ -0,0 +1,8 @@
|
||||
"""Galileo receiver baseband reference model."""
|
||||
|
||||
from .config import ReceiverConfig, load_config
|
||||
from .signal_profiles import SignalProfile, galileo_e1_profile
|
||||
|
||||
__all__ = ["ReceiverConfig", "SignalProfile", "galileo_e1_profile", "load_config"]
|
||||
__version__ = "0.1.0"
|
||||
|
||||
@@ -0,0 +1,7 @@
|
||||
"""FFT parallel code-phase acquisition."""
|
||||
|
||||
from .acquisition_manager import AcquisitionManager
|
||||
from .fft_pcps import PcpsAcquisition, PcpsSearchGrid
|
||||
|
||||
__all__ = ["AcquisitionManager", "PcpsAcquisition", "PcpsSearchGrid"]
|
||||
|
||||
@@ -0,0 +1,67 @@
|
||||
"""PRN search orchestration and structured results."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.config import AcquisitionConfig
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.signals.code_provider import CodeProvider
|
||||
from galileo_rx.types import AcquisitionResult
|
||||
from .detector import boc_side_peak_suspected, parabolic_offset
|
||||
from .fft_pcps import PcpsAcquisition, PcpsSearchGrid
|
||||
from .metrics import calculate_metrics, estimate_acquisition_cn0_db_hz
|
||||
|
||||
|
||||
class AcquisitionManager:
|
||||
def __init__(self, profile: SignalProfile, provider: CodeProvider, sampling_frequency_hz: float, config: AcquisitionConfig) -> None:
|
||||
self.profile = profile
|
||||
self.sampling_frequency_hz = sampling_frequency_hz
|
||||
self.config = config
|
||||
self.engine = PcpsAcquisition(profile, provider, sampling_frequency_hz, config)
|
||||
|
||||
def acquire_one(self, samples: npt.NDArray[np.generic], component: str, prn: int, carrier_center_hz: float = 0.0) -> tuple[AcquisitionResult, PcpsSearchGrid]:
|
||||
grid = self.engine.search(samples, component, prn, carrier_center_hz)
|
||||
raw_doppler_index, raw_code_index = np.unravel_index(int(np.argmax(grid.power)), grid.power.shape)
|
||||
profile = grid.power[raw_doppler_index]
|
||||
samples_per_chip = self.sampling_frequency_hz / self.profile.code_rate_hz
|
||||
guard = max(1, int(np.ceil(self.config.peak_exclusion_chips * samples_per_chip)))
|
||||
metrics = calculate_metrics(profile, raw_code_index, guard)
|
||||
selected_metric = metrics[self.config.metric]
|
||||
code_offset = 0.0
|
||||
doppler_offset = 0.0
|
||||
if self.config.parabolic_interpolation:
|
||||
code_offset = parabolic_offset(profile[(raw_code_index - 1) % profile.size], profile[raw_code_index], profile[(raw_code_index + 1) % profile.size])
|
||||
if 0 < raw_doppler_index < grid.doppler_bins_hz.size - 1:
|
||||
doppler_offset = parabolic_offset(
|
||||
grid.power[raw_doppler_index - 1, raw_code_index],
|
||||
grid.power[raw_doppler_index, raw_code_index],
|
||||
grid.power[raw_doppler_index + 1, raw_code_index],
|
||||
)
|
||||
code_samples = float((raw_code_index + code_offset) % profile.size)
|
||||
doppler_hz = float(grid.doppler_bins_hz[raw_doppler_index] + doppler_offset * self.config.doppler_step_hz)
|
||||
cn0 = estimate_acquisition_cn0_db_hz(float(profile[raw_code_index]), metrics["noise_mean"], self.config.coherent_integration_ms * 1e-3)
|
||||
result = AcquisitionResult(
|
||||
detected=selected_metric >= self.config.detection_threshold,
|
||||
prn=prn,
|
||||
component=component,
|
||||
code_phase_samples=code_samples,
|
||||
code_phase_chips=code_samples / samples_per_chip,
|
||||
doppler_hz=doppler_hz,
|
||||
carrier_frequency_hz=carrier_center_hz + doppler_hz,
|
||||
peak_power=float(profile[raw_code_index]),
|
||||
second_peak_power=metrics["second_peak"],
|
||||
detection_metric=selected_metric,
|
||||
estimated_cn0_db_hz=cn0,
|
||||
coherent_integration_ms=self.config.coherent_integration_ms,
|
||||
noncoherent_accumulations=self.config.noncoherent_accumulations,
|
||||
raw_code_phase_samples=int(raw_code_index),
|
||||
raw_doppler_hz=float(grid.doppler_bins_hz[raw_doppler_index]),
|
||||
metrics=metrics,
|
||||
side_peak_suspected=boc_side_peak_suspected(profile, raw_code_index, samples_per_chip),
|
||||
)
|
||||
return result, grid
|
||||
|
||||
def acquire(self, samples: npt.NDArray[np.generic], component: str, prns: list[int], carrier_center_hz: float = 0.0) -> list[AcquisitionResult]:
|
||||
return [self.acquire_one(samples, component, prn, carrier_center_hz)[0] for prn in prns]
|
||||
@@ -0,0 +1,21 @@
|
||||
"""Detection decision and sub-bin interpolation."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def parabolic_offset(left: float, center: float, right: float) -> float:
|
||||
denominator = left - 2.0 * center + right
|
||||
if denominator == 0:
|
||||
return 0.0
|
||||
return float(np.clip(0.5 * (left - right) / denominator, -0.5, 0.5))
|
||||
|
||||
|
||||
def boc_side_peak_suspected(profile: np.ndarray, peak_index: int, samples_per_chip: float) -> bool:
|
||||
"""Flag a competitive BOC(1,1) half-chip neighbor; no false unambiguous claim."""
|
||||
offset = max(1, int(round(0.5 * samples_per_chip)))
|
||||
peak = profile[peak_index]
|
||||
neighbors = [profile[(peak_index - offset) % profile.size], profile[(peak_index + offset) % profile.size]]
|
||||
return bool(max(neighbors) > 0.75 * peak)
|
||||
|
||||
@@ -0,0 +1,138 @@
|
||||
"""Readable floating-point FFT parallel code-phase search."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.config import AcquisitionConfig
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.signals.code_provider import CodeProvider
|
||||
from galileo_rx.signals.galileo_e1 import e1_local_replica
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class PcpsSearchGrid:
|
||||
doppler_bins_hz: npt.NDArray[np.float64]
|
||||
power: npt.NDArray[np.float64]
|
||||
complex_sums: npt.NDArray[np.complex128] | None
|
||||
coherent_samples: int
|
||||
fft_length: int
|
||||
|
||||
|
||||
class PcpsAcquisition:
|
||||
def __init__(self, profile: SignalProfile, provider: CodeProvider, sampling_frequency_hz: float, config: AcquisitionConfig) -> None:
|
||||
self.profile = profile
|
||||
self.provider = provider
|
||||
self.sampling_frequency_hz = sampling_frequency_hz
|
||||
self.config = config
|
||||
config.validate()
|
||||
self._code_fft_cache: dict[tuple[str, int, int, int, tuple[int, ...], float], npt.NDArray[np.complex128]] = {}
|
||||
|
||||
@property
|
||||
def coherent_samples(self) -> int:
|
||||
exact = self.sampling_frequency_hz * self.config.coherent_integration_ms * 1e-3
|
||||
rounded = int(round(exact))
|
||||
if not np.isclose(exact, rounded, atol=1e-6):
|
||||
raise ValueError("coherent integration time is not an integer number of samples")
|
||||
return rounded
|
||||
|
||||
def _fft_length(self, sample_count: int) -> int:
|
||||
if self.config.fft_length is not None:
|
||||
if self.config.fft_length < sample_count:
|
||||
raise ValueError("fft_length cannot be shorter than coherent block")
|
||||
return self.config.fft_length
|
||||
if self.config.zero_padding_mode == "next_power_of_two":
|
||||
return 1 << (sample_count - 1).bit_length()
|
||||
return sample_count
|
||||
|
||||
def _replica_fft(self, component: str, prn: int, sample_count: int, fft_length: int, primary_signs: tuple[int, ...], doppler_hz: float) -> npt.NDArray[np.complex128]:
|
||||
effective_code_rate_hz = self.profile.code_rate_hz
|
||||
if self.config.code_doppler_aiding:
|
||||
effective_code_rate_hz *= 1.0 + doppler_hz / self.profile.carrier_frequency_hz
|
||||
key = (component, prn, sample_count, fft_length, primary_signs, effective_code_rate_hz)
|
||||
if key not in self._code_fft_cache:
|
||||
replica = e1_local_replica(
|
||||
self.profile, self.provider, component, prn, sample_count,
|
||||
self.sampling_frequency_hz, code_rate_hz=effective_code_rate_hz,
|
||||
)
|
||||
if primary_signs != (1,):
|
||||
primary_samples = int(round(self.profile.component(component).primary_code_period_s * self.sampling_frequency_hz))
|
||||
signs = np.repeat(np.asarray(primary_signs, dtype=np.float64), primary_samples)
|
||||
replica = replica * np.resize(signs, sample_count)
|
||||
self._code_fft_cache[key] = np.conjugate(np.fft.fft(replica, fft_length))
|
||||
return self._code_fft_cache[key]
|
||||
|
||||
def _sign_hypotheses(self, component: str, prn: int, primary_periods: int, first_primary_epoch: int) -> list[tuple[int, ...]]:
|
||||
spec = self.profile.component(component)
|
||||
if primary_periods <= 1:
|
||||
return [(1,)]
|
||||
if spec.channel_type == "pilot" and self.config.secondary_code_mode in {"wipeoff", "search"}:
|
||||
secondary = self.provider.secondary_code(self.profile, component, prn)
|
||||
if secondary is None:
|
||||
raise ValueError("secondary-code handling requested but provider returned no secondary code")
|
||||
phases = range(secondary.size) if self.config.secondary_code_mode == "search" else range(1)
|
||||
return [tuple(int(secondary[(phase + first_primary_epoch + epoch) % secondary.size]) for epoch in range(primary_periods)) for phase in phases]
|
||||
if spec.channel_type == "data" and self.config.bit_transition_mode == "hypothesis":
|
||||
if primary_periods > 8:
|
||||
raise ValueError("bit-transition hypothesis search is limited to 8 primary periods")
|
||||
return [tuple([1] + [1 if (mask >> bit) & 1 else -1 for bit in range(primary_periods - 1)]) for mask in range(1 << (primary_periods - 1))]
|
||||
return [(1,) * primary_periods]
|
||||
|
||||
def search(self, samples: npt.NDArray[np.generic], component: str, prn: int, carrier_center_hz: float = 0.0) -> PcpsSearchGrid:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
coherent_samples = self.coherent_samples
|
||||
required = coherent_samples * self.config.noncoherent_accumulations
|
||||
if array.size < required:
|
||||
raise ValueError(f"acquisition needs {required} samples, received {array.size}")
|
||||
fft_length = self._fft_length(coherent_samples)
|
||||
primary_samples = int(round(self.profile.component(component).primary_code_period_s * self.sampling_frequency_hz))
|
||||
primary_periods_exact = coherent_samples / primary_samples
|
||||
primary_periods = int(round(primary_periods_exact))
|
||||
if not np.isclose(primary_periods_exact, primary_periods, atol=1e-9):
|
||||
raise ValueError("coherent integration must contain a whole number of primary-code periods")
|
||||
bins = np.arange(
|
||||
self.config.doppler_min_hz,
|
||||
self.config.doppler_max_hz + 0.5 * self.config.doppler_step_hz,
|
||||
self.config.doppler_step_hz,
|
||||
dtype=np.float64,
|
||||
)
|
||||
# Code delay is periodic over one primary code even when coherent time
|
||||
# spans several signed primary epochs. Keeping one period prevents an
|
||||
# artificial duplicate peak from corrupting the second-peak metric.
|
||||
code_phase_bins = primary_samples
|
||||
power = np.zeros((bins.size, code_phase_bins), dtype=np.float64)
|
||||
complex_sums = np.zeros((bins.size, code_phase_bins), dtype=np.complex128) if self.config.preserve_complex_sums else None
|
||||
time_s = np.arange(coherent_samples, dtype=np.float64) / self.sampling_frequency_hz
|
||||
batch_size = 1 if self.config.parallelism == "sequential" else max(1, self.config.doppler_batch_size)
|
||||
for first_bin in range(0, bins.size, batch_size):
|
||||
batch_bins = bins[first_bin : first_bin + batch_size]
|
||||
carrier = np.exp(-2j * np.pi * (carrier_center_hz + batch_bins[:, None]) * time_s[None, :])
|
||||
previous_correlation: npt.NDArray[np.complex128] | None = None
|
||||
for accumulation in range(self.config.noncoherent_accumulations):
|
||||
start = accumulation * coherent_samples
|
||||
# Absolute sample time contributes only a common phase per coherent block,
|
||||
# which magnitude-squared noncoherent accumulation intentionally discards.
|
||||
mixed = array[start : start + coherent_samples][None, :] * carrier
|
||||
sample_fft = np.fft.fft(mixed, fft_length, axis=1)
|
||||
hypotheses = self._sign_hypotheses(component, prn, primary_periods, accumulation * primary_periods)
|
||||
code_ffts = np.asarray([[
|
||||
self._replica_fft(component, prn, coherent_samples, fft_length, signs, float(doppler_hz))
|
||||
for signs in hypotheses
|
||||
] for doppler_hz in batch_bins])
|
||||
stack = np.fft.ifft(sample_fft[:, None, :] * code_ffts, axis=-1)[:, :, :code_phase_bins]
|
||||
if stack.shape[1] == 1:
|
||||
correlation = stack[:, 0, :]
|
||||
else:
|
||||
winners = np.argmax(np.abs(stack) ** 2, axis=1)
|
||||
correlation = np.take_along_axis(stack, winners[:, None, :], axis=1)[:, 0, :]
|
||||
if self.config.bit_transition_mode == "differential" and previous_correlation is not None:
|
||||
power[first_bin : first_bin + batch_bins.size] += np.abs(correlation * np.conjugate(previous_correlation))
|
||||
else:
|
||||
power[first_bin : first_bin + batch_bins.size] += np.abs(correlation) ** 2
|
||||
if complex_sums is not None:
|
||||
complex_sums[first_bin : first_bin + batch_bins.size] += correlation
|
||||
previous_correlation = correlation
|
||||
return PcpsSearchGrid(bins, power, complex_sums, coherent_samples, fft_length)
|
||||
@@ -0,0 +1,41 @@
|
||||
"""Acquisition peak and noise metrics."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
def circular_guard_mask(length: int, peak_index: int, guard_samples: int) -> npt.NDArray[np.bool_]:
|
||||
indices = np.arange(length)
|
||||
circular_distance = np.minimum((indices - peak_index) % length, (peak_index - indices) % length)
|
||||
return circular_distance > guard_samples
|
||||
|
||||
|
||||
def calculate_metrics(power: npt.NDArray[np.float64], peak_index: int, guard_samples: int) -> dict[str, float]:
|
||||
valid = circular_guard_mask(power.size, peak_index, guard_samples)
|
||||
noise = power[valid]
|
||||
peak = float(power[peak_index])
|
||||
second = float(np.max(noise)) if noise.size else 0.0
|
||||
mean = float(np.mean(noise)) if noise.size else 0.0
|
||||
median = float(np.median(noise)) if noise.size else 0.0
|
||||
std = float(np.std(noise)) if noise.size else 0.0
|
||||
epsilon = np.finfo(np.float64).tiny
|
||||
return {
|
||||
"peak_to_second": peak / max(second, epsilon),
|
||||
"peak_to_mean": peak / max(mean, epsilon),
|
||||
"peak_to_sidelobe": (peak - mean) / max(second - mean, epsilon),
|
||||
"cfar": (peak - median) / max(std, epsilon),
|
||||
"noise_mean": mean,
|
||||
"noise_median": median,
|
||||
"noise_std": std,
|
||||
"second_peak": second,
|
||||
}
|
||||
|
||||
|
||||
def estimate_acquisition_cn0_db_hz(peak_power: float, noise_power: float, coherent_time_s: float) -> float | None:
|
||||
signal_power = peak_power - noise_power
|
||||
if signal_power <= 0 or noise_power <= 0 or coherent_time_s <= 0:
|
||||
return None
|
||||
return float(10.0 * np.log10(signal_power / noise_power / coherent_time_s))
|
||||
|
||||
@@ -0,0 +1,189 @@
|
||||
"""Command-line entry points for generation, acquisition, tracking, and comparison."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import logging
|
||||
from dataclasses import asdict
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
|
||||
from galileo_rx.acquisition.acquisition_manager import AcquisitionManager
|
||||
from galileo_rx.config import load_config
|
||||
from galileo_rx.diagnostics.reports import export_json, export_tracking_csv
|
||||
from galileo_rx.fixed_point.bittrue_pipeline import compare_acquisition
|
||||
from galileo_rx.fixed_point.fixed_acquisition import FixedPcpsAcquisition
|
||||
from galileo_rx.frontend.frontend_chain import FrontendChain
|
||||
from galileo_rx.io.sample_reader import read_samples
|
||||
from galileo_rx.io.sample_writer import write_samples
|
||||
from galileo_rx.signal_profiles import get_signal_profile
|
||||
from galileo_rx.signals.code_provider import SyntheticCodeProvider
|
||||
from galileo_rx.signals.galileo_e1 import e1_local_replica
|
||||
from galileo_rx.signals.signal_generator import ChannelModelConfig, SatelliteSignalConfig, SyntheticSignalGenerator
|
||||
from galileo_rx.tracking.tracking_channel import TrackingChannel
|
||||
from galileo_rx.types import AcquisitionResult
|
||||
from galileo_rx.verification.signalsim_validation import locate_signalsim, validate_signalsim
|
||||
|
||||
LOGGER = logging.getLogger("galileo_rx")
|
||||
|
||||
|
||||
def _load_and_frontend(config_path: str, input_path: str):
|
||||
config = load_config(config_path)
|
||||
raw = read_samples(input_path, config.rf.sample_format, config.rf.sample_endianness)
|
||||
frontend = FrontendChain(config.rf, config.frontend)
|
||||
return config, frontend.process(raw), frontend.output_sampling_frequency_hz
|
||||
|
||||
|
||||
def _result_from_json(path: str) -> AcquisitionResult:
|
||||
raw = json.loads(Path(path).read_text(encoding="utf-8"))
|
||||
if isinstance(raw, list):
|
||||
detected = [item for item in raw if item.get("detected")]
|
||||
if not detected:
|
||||
raise ValueError("acquisition result file contains no detected result")
|
||||
raw = detected[0]
|
||||
return AcquisitionResult(**raw)
|
||||
|
||||
|
||||
def command_generate(args: argparse.Namespace) -> None:
|
||||
config = load_config(args.config)
|
||||
profile = get_signal_profile(config.signal_profile)
|
||||
provider = SyntheticCodeProvider()
|
||||
satellite = SatelliteSignalConfig(
|
||||
prn=args.prn or config.prns[0], component=config.component, code_delay_chips=args.code_delay_chips,
|
||||
doppler_hz=args.doppler_hz, doppler_rate_hz_s=args.doppler_rate_hz_s,
|
||||
carrier_phase_rad=args.carrier_phase_rad, amplitude=args.amplitude, cn0_db_hz=args.cn0_db_hz,
|
||||
)
|
||||
channel = ChannelModelConfig(seed=args.seed, quantization_bits=args.quantization_bits, clip_amplitude=args.clip_amplitude)
|
||||
generator = SyntheticSignalGenerator(
|
||||
profile, provider, config.rf.sampling_frequency_hz, config.rf.intermediate_frequency_hz,
|
||||
config.rf.input_is_real_if, [satellite], channel,
|
||||
)
|
||||
sample_count = int(round(args.duration_ms * 1e-3 * config.rf.sampling_frequency_hz))
|
||||
samples, truth = generator.generate(sample_count)
|
||||
write_samples(args.output, samples, config.rf.sample_endianness)
|
||||
export_json(str(args.output) + ".truth.json", truth)
|
||||
LOGGER.info("wrote %d samples and truth metadata", sample_count)
|
||||
|
||||
|
||||
def command_acquire(args: argparse.Namespace) -> None:
|
||||
config, samples, fs = _load_and_frontend(args.config, args.input)
|
||||
profile = get_signal_profile(config.signal_profile)
|
||||
manager = AcquisitionManager(profile, SyntheticCodeProvider(), fs, config.acquisition)
|
||||
results = manager.acquire(samples, config.component, config.prns)
|
||||
payload = [asdict(result) for result in results]
|
||||
text = json.dumps(payload, indent=2)
|
||||
print(text)
|
||||
if args.output:
|
||||
Path(args.output).write_text(text, encoding="utf-8")
|
||||
|
||||
|
||||
def command_track(args: argparse.Namespace) -> None:
|
||||
config, samples, fs = _load_and_frontend(args.config, args.input)
|
||||
acquisition = _result_from_json(args.acquisition_result)
|
||||
channel = TrackingChannel(get_signal_profile(config.signal_profile), SyntheticCodeProvider(), fs, config.tracking, acquisition, args.start_sample)
|
||||
measurements = channel.process(samples[args.start_sample :])
|
||||
if args.output:
|
||||
export_tracking_csv(args.output, measurements)
|
||||
if measurements:
|
||||
print(json.dumps(asdict(measurements[-1]), indent=2, default=str))
|
||||
|
||||
|
||||
def command_run(args: argparse.Namespace) -> None:
|
||||
config, samples, fs = _load_and_frontend(args.config, args.input)
|
||||
profile = get_signal_profile(config.signal_profile)
|
||||
provider = SyntheticCodeProvider()
|
||||
acquisition_manager = AcquisitionManager(profile, provider, fs, config.acquisition)
|
||||
results = acquisition_manager.acquire(samples, config.component, config.prns)
|
||||
detected = [result for result in results if result.detected]
|
||||
if not detected:
|
||||
raise SystemExit("no configured PRN passed the acquisition threshold")
|
||||
consumed = acquisition_manager.engine.coherent_samples * config.acquisition.noncoherent_accumulations
|
||||
summary = []
|
||||
for result in detected:
|
||||
channel = TrackingChannel(profile, provider, fs, config.tracking, result, consumed)
|
||||
measurements = channel.process(samples[consumed:])
|
||||
summary.append({"acquisition": asdict(result), "last_tracking": asdict(measurements[-1]) if measurements else None})
|
||||
print(json.dumps(summary, indent=2, default=str))
|
||||
|
||||
|
||||
def command_compare_fixed(args: argparse.Namespace) -> None:
|
||||
config, samples, fs = _load_and_frontend(args.config, args.input)
|
||||
profile = get_signal_profile(config.signal_profile)
|
||||
provider = SyntheticCodeProvider()
|
||||
engine = AcquisitionManager(profile, provider, fs, config.acquisition)
|
||||
result, grid = engine.acquire_one(samples, config.component, config.prns[0])
|
||||
count = engine.engine.coherent_samples
|
||||
replica = e1_local_replica(profile, provider, config.component, config.prns[0], count, fs)
|
||||
fixed = FixedPcpsAcquisition().search(samples[:count], replica, fs, grid.doppler_bins_hz)
|
||||
reference = grid.power / config.acquisition.noncoherent_accumulations / count**2
|
||||
comparison = compare_acquisition(reference, fixed, config.acquisition.detection_threshold)
|
||||
print(json.dumps({"float_result": asdict(result), "comparison": comparison.to_dict()}, indent=2))
|
||||
|
||||
|
||||
def command_validate_signalsim(args: argparse.Namespace) -> None:
|
||||
source = Path(args.signalsim_path).resolve() if args.signalsim_path else locate_signalsim(Path.cwd())
|
||||
report, passed = validate_signalsim(source, args.config, args.output_dir, mode=args.mode)
|
||||
print(json.dumps({"overall_pass": passed, "passes": report["passes"], "report": str(Path(args.output_dir).resolve() / "signalsim_validation_report.json")}, indent=2))
|
||||
if not passed:
|
||||
raise SystemExit(1)
|
||||
|
||||
|
||||
def build_parser() -> argparse.ArgumentParser:
|
||||
parser = argparse.ArgumentParser(prog="galileo-rx")
|
||||
parser.add_argument("--log-level", default="INFO")
|
||||
commands = parser.add_subparsers(dest="command", required=True)
|
||||
generate = commands.add_parser("generate")
|
||||
generate.add_argument("--config", required=True)
|
||||
generate.add_argument("--output", required=True)
|
||||
generate.add_argument("--duration-ms", type=float, default=100.0)
|
||||
generate.add_argument("--prn", type=int)
|
||||
generate.add_argument("--code-delay-chips", type=float, default=100.0)
|
||||
generate.add_argument("--doppler-hz", type=float, default=1500.0)
|
||||
generate.add_argument("--doppler-rate-hz-s", type=float, default=0.0)
|
||||
generate.add_argument("--carrier-phase-rad", type=float, default=0.0)
|
||||
generate.add_argument("--amplitude", type=float, default=0.5)
|
||||
generate.add_argument("--cn0-db-hz", type=float, default=48.0)
|
||||
generate.add_argument("--seed", type=int, default=1)
|
||||
generate.add_argument("--quantization-bits", type=int)
|
||||
generate.add_argument("--clip-amplitude", type=float)
|
||||
generate.set_defaults(handler=command_generate)
|
||||
acquire = commands.add_parser("acquire")
|
||||
acquire.add_argument("--config", required=True)
|
||||
acquire.add_argument("--input", required=True)
|
||||
acquire.add_argument("--output")
|
||||
acquire.set_defaults(handler=command_acquire)
|
||||
track = commands.add_parser("track")
|
||||
track.add_argument("--config", required=True)
|
||||
track.add_argument("--input", required=True)
|
||||
track.add_argument("--acquisition-result", required=True)
|
||||
track.add_argument("--start-sample", type=int, default=0)
|
||||
track.add_argument("--output")
|
||||
track.set_defaults(handler=command_track)
|
||||
run = commands.add_parser("run")
|
||||
run.add_argument("--config", required=True)
|
||||
run.add_argument("--input", required=True)
|
||||
run.set_defaults(handler=command_run)
|
||||
compare = commands.add_parser("compare-fixed")
|
||||
compare.add_argument("--config", required=True)
|
||||
compare.add_argument("--input", required=True)
|
||||
compare.set_defaults(handler=command_compare_fixed)
|
||||
validate = commands.add_parser("validate-signalsim")
|
||||
validate.add_argument("--signalsim-path")
|
||||
validate.add_argument("--config", default="configs/signalsim_validation.yaml")
|
||||
validate.add_argument("--output-dir", default="reports/signalsim_validation")
|
||||
validate.add_argument("--mode", choices=["quick", "full"])
|
||||
validate.set_defaults(handler=command_validate_signalsim)
|
||||
return parser
|
||||
|
||||
|
||||
def main(argv: list[str] | None = None) -> None:
|
||||
parser = build_parser()
|
||||
args = parser.parse_args(argv)
|
||||
logging.basicConfig(level=getattr(logging, args.log_level.upper()), format="%(levelname)s %(name)s: %(message)s")
|
||||
args.handler(args)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,177 @@
|
||||
"""Typed configuration and strict YAML loading."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import asdict, dataclass, field, fields
|
||||
from pathlib import Path
|
||||
from typing import Any, Literal, TypeVar
|
||||
|
||||
import yaml
|
||||
|
||||
T = TypeVar("T")
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class RfConfig:
|
||||
sample_format: str = "complex128"
|
||||
input_is_real_if: bool = False
|
||||
sampling_frequency_hz: float = 16.368e6
|
||||
intermediate_frequency_hz: float = 4.092e6
|
||||
sample_bit_width: int = 8
|
||||
sample_endianness: Literal["little", "big"] = "little"
|
||||
|
||||
def validate(self) -> None:
|
||||
if self.sampling_frequency_hz <= 0:
|
||||
raise ValueError("sampling_frequency_hz must be positive")
|
||||
if self.input_is_real_if and not 0 < self.intermediate_frequency_hz < self.sampling_frequency_hz / 2:
|
||||
raise ValueError("real-IF intermediate_frequency_hz must lie between DC and Nyquist")
|
||||
if not 2 <= self.sample_bit_width <= 32:
|
||||
raise ValueError("sample_bit_width must be in [2, 32]")
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class AcquisitionConfig:
|
||||
doppler_min_hz: float = -10_000.0
|
||||
doppler_max_hz: float = 10_000.0
|
||||
doppler_step_hz: float = 250.0
|
||||
coherent_integration_ms: float = 4.0
|
||||
noncoherent_accumulations: int = 5
|
||||
detection_threshold: float = 2.5
|
||||
peak_exclusion_chips: float = 1.0
|
||||
fft_length: int | None = None
|
||||
zero_padding_mode: Literal["none", "next_power_of_two", "explicit"] = "next_power_of_two"
|
||||
metric: Literal["peak_to_second", "peak_to_mean", "cfar"] = "peak_to_second"
|
||||
parallelism: Literal["sequential", "doppler_vectorized", "prn_vectorized", "batched"] = "sequential"
|
||||
doppler_batch_size: int = 8
|
||||
preserve_complex_sums: bool = False
|
||||
bit_transition_mode: Literal["none", "hypothesis", "differential"] = "none"
|
||||
secondary_code_mode: Literal["none", "wipeoff", "search"] = "none"
|
||||
parabolic_interpolation: bool = True
|
||||
boc_side_peak_strategy: Literal["flag", "prefer_main", "none"] = "flag"
|
||||
code_doppler_aiding: bool = True
|
||||
|
||||
def validate(self) -> None:
|
||||
if self.doppler_max_hz < self.doppler_min_hz or self.doppler_step_hz <= 0:
|
||||
raise ValueError("invalid Doppler search range or step")
|
||||
if self.coherent_integration_ms <= 0 or self.noncoherent_accumulations < 1:
|
||||
raise ValueError("integration settings must be positive")
|
||||
if self.peak_exclusion_chips <= 0 or self.detection_threshold <= 0:
|
||||
raise ValueError("detector parameters must be positive")
|
||||
if self.zero_padding_mode == "explicit" and self.fft_length is None:
|
||||
raise ValueError("explicit zero padding requires fft_length")
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class TrackingConfig:
|
||||
integration_ms: float = 4.0
|
||||
early_late_spacing_chips: float = 0.25
|
||||
very_early_late_spacing_chips: float = 0.75
|
||||
dll_noise_bandwidth_hz: float = 1.0
|
||||
pll_noise_bandwidth_hz: float = 15.0
|
||||
fll_noise_bandwidth_hz: float = 20.0
|
||||
pll_order: int = 3
|
||||
dll_order: int = 2
|
||||
use_fll_assist: bool = True
|
||||
pull_in_duration_ms: float = 500.0
|
||||
damping_ratio: float = 0.707
|
||||
carrier_aiding: bool = True
|
||||
secondary_code_mode: Literal["none", "wipeoff", "search"] = "wipeoff"
|
||||
loss_of_lock_epochs: int = 10
|
||||
weak_signal_cn0_db_hz: float = 25.0
|
||||
|
||||
def validate(self) -> None:
|
||||
if self.integration_ms <= 0:
|
||||
raise ValueError("tracking integration_ms must be positive")
|
||||
if not 0 < self.early_late_spacing_chips < self.very_early_late_spacing_chips:
|
||||
raise ValueError("tap spacings must satisfy 0 < E/L < VE/VL")
|
||||
if self.pll_order not in (2, 3) or self.dll_order not in (1, 2):
|
||||
raise ValueError("supported PLL orders are 2/3 and DLL orders are 1/2")
|
||||
if min(self.dll_noise_bandwidth_hz, self.pll_noise_bandwidth_hz, self.fll_noise_bandwidth_hz) <= 0:
|
||||
raise ValueError("loop noise bandwidths must be positive")
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FrontendConfig:
|
||||
remove_dc: bool = False
|
||||
dc_alpha: float = 0.999
|
||||
agc_enabled: bool = False
|
||||
agc_target_rms: float = 0.5
|
||||
fir_taps: list[float] = field(default_factory=list)
|
||||
decimation: int = 1
|
||||
nco_phase_width: int = 32
|
||||
nco_lut_bits: int = 12
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedPointConfig:
|
||||
enabled: bool = False
|
||||
input_total_bits: int = 8
|
||||
input_fractional_bits: int = 7
|
||||
coefficient_total_bits: int = 16
|
||||
coefficient_fractional_bits: int = 15
|
||||
accumulator_total_bits: int = 32
|
||||
accumulator_fractional_bits: int = 20
|
||||
overflow: Literal["saturate", "wrap"] = "saturate"
|
||||
rounding: Literal["nearest", "truncate", "convergent"] = "nearest"
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class ReceiverConfig:
|
||||
signal_profile: str = "GALILEO_E1"
|
||||
component: Literal["E1B", "E1C"] = "E1C"
|
||||
prns: list[int] = field(default_factory=lambda: [1])
|
||||
rf: RfConfig = field(default_factory=RfConfig)
|
||||
acquisition: AcquisitionConfig = field(default_factory=AcquisitionConfig)
|
||||
tracking: TrackingConfig = field(default_factory=TrackingConfig)
|
||||
frontend: FrontendConfig = field(default_factory=FrontendConfig)
|
||||
fixed_point: FixedPointConfig = field(default_factory=FixedPointConfig)
|
||||
|
||||
def validate(self) -> None:
|
||||
self.rf.validate()
|
||||
self.acquisition.validate()
|
||||
self.tracking.validate()
|
||||
if not self.prns or any(prn <= 0 for prn in self.prns):
|
||||
raise ValueError("prns must contain positive identifiers")
|
||||
if self.acquisition.coherent_integration_ms > 4.0:
|
||||
if self.component == "E1B" and self.acquisition.bit_transition_mode == "none":
|
||||
raise ValueError("E1B coherent integration beyond 4 ms requires a bit-transition mode")
|
||||
if self.component == "E1C" and self.acquisition.secondary_code_mode == "none":
|
||||
raise ValueError("E1C coherent integration beyond 4 ms requires secondary-code handling")
|
||||
if self.tracking.integration_ms > 4.0:
|
||||
if self.component == "E1B":
|
||||
raise ValueError("E1B tracking integration beyond 4 ms is unsafe without symbol synchronization")
|
||||
if self.tracking.secondary_code_mode == "none":
|
||||
raise ValueError("E1C tracking beyond 4 ms requires secondary-code handling")
|
||||
|
||||
def to_dict(self) -> dict[str, Any]:
|
||||
return asdict(self)
|
||||
|
||||
|
||||
def _construct(cls: type[T], values: dict[str, Any] | None) -> T:
|
||||
values = values or {}
|
||||
known = {item.name for item in fields(cls)}
|
||||
unknown = set(values) - known
|
||||
if unknown:
|
||||
raise ValueError(f"unknown {cls.__name__} fields: {sorted(unknown)}")
|
||||
return cls(**values)
|
||||
|
||||
|
||||
def load_config(path: str | Path) -> ReceiverConfig:
|
||||
with Path(path).open("r", encoding="utf-8") as stream:
|
||||
raw = yaml.safe_load(stream) or {}
|
||||
known = {item.name for item in fields(ReceiverConfig)}
|
||||
unknown = set(raw) - known
|
||||
if unknown:
|
||||
raise ValueError(f"unknown ReceiverConfig fields: {sorted(unknown)}")
|
||||
config = ReceiverConfig(
|
||||
signal_profile=raw.get("signal_profile", "GALILEO_E1"),
|
||||
component=raw.get("component", "E1C"),
|
||||
prns=raw.get("prns", [1]),
|
||||
rf=_construct(RfConfig, raw.get("rf")),
|
||||
acquisition=_construct(AcquisitionConfig, raw.get("acquisition")),
|
||||
tracking=_construct(TrackingConfig, raw.get("tracking")),
|
||||
frontend=_construct(FrontendConfig, raw.get("frontend")),
|
||||
fixed_point=_construct(FixedPointConfig, raw.get("fixed_point")),
|
||||
)
|
||||
config.validate()
|
||||
return config
|
||||
@@ -0,0 +1,2 @@
|
||||
"""Optional plots and machine-readable reports."""
|
||||
|
||||
@@ -0,0 +1,96 @@
|
||||
"""Headless-safe diagnostics; matplotlib is imported only on demand."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def acquisition_map(power: np.ndarray, doppler_bins_hz: np.ndarray, title: str = "Acquisition map") -> Any:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
figure, axis = plt.subplots()
|
||||
image = axis.imshow(power, aspect="auto", origin="lower", extent=[0, power.shape[1], doppler_bins_hz[0], doppler_bins_hz[-1]])
|
||||
axis.set(xlabel="Code phase (samples)", ylabel="Doppler (Hz)", title=title)
|
||||
figure.colorbar(image, ax=axis, label="Correlation power")
|
||||
return figure
|
||||
|
||||
|
||||
def acquisition_profiles(power: np.ndarray, doppler_bins_hz: np.ndarray) -> Any:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
doppler_index, code_index = np.unravel_index(int(np.argmax(power)), power.shape)
|
||||
figure, axes = plt.subplots(1, 2)
|
||||
axes[0].plot(power[doppler_index])
|
||||
axes[0].set(xlabel="Code phase (samples)", ylabel="Power")
|
||||
axes[1].plot(doppler_bins_hz, power[:, code_index])
|
||||
axes[1].set(xlabel="Doppler (Hz)", ylabel="Power")
|
||||
return figure
|
||||
|
||||
|
||||
def acquisition_peak_neighborhood(power: np.ndarray, samples_per_chip: float, span_chips: float = 2.0) -> Any:
|
||||
"""Circular peak neighborhood including the BOC half-chip region."""
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
doppler_index, code_index = np.unravel_index(int(np.argmax(power)), power.shape)
|
||||
radius = max(1, int(np.ceil(span_chips * samples_per_chip)))
|
||||
offsets = np.arange(-radius, radius + 1)
|
||||
values = power[doppler_index, (code_index + offsets) % power.shape[1]]
|
||||
figure, axis = plt.subplots()
|
||||
axis.plot(offsets / samples_per_chip, values / np.max(values), marker=".")
|
||||
axis.axvline(-0.5, color="tab:red", linestyle="--", alpha=0.5)
|
||||
axis.axvline(0.5, color="tab:red", linestyle="--", alpha=0.5)
|
||||
axis.set(xlabel="Offset from peak (chips)", ylabel="Normalized power", title="Peak and BOC side-peak neighborhood")
|
||||
return figure
|
||||
|
||||
|
||||
def prompt_constellation(measurements: list[Any], ignored_epochs: int = 0) -> Any:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
selected = measurements[ignored_epochs:]
|
||||
figure, axis = plt.subplots()
|
||||
axis.scatter([item.prompt_i for item in selected], [item.prompt_q for item in selected], s=5, alpha=0.5)
|
||||
axis.set(xlabel="Prompt I", ylabel="Prompt Q", title="Prompt correlator constellation")
|
||||
axis.set_aspect("equal", adjustable="datalim")
|
||||
return figure
|
||||
|
||||
|
||||
def tracking_diagnostics(measurements: list[Any]) -> Any:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
time_s = [item.time_s for item in measurements]
|
||||
figure, axes = plt.subplots(4, 2, sharex=True)
|
||||
series = [
|
||||
("Prompt I/Q", [item.prompt_i for item in measurements], [item.prompt_q for item in measurements]),
|
||||
("Early/Prompt/Late magnitudes", [item.early_magnitude for item in measurements], [np.hypot(item.prompt_i, item.prompt_q) for item in measurements], [item.late_magnitude for item in measurements]),
|
||||
("PLL (rad)", [item.pll_error_rad for item in measurements], None),
|
||||
("FLL (Hz)", [item.fll_error_hz for item in measurements], None),
|
||||
("DLL (chips)", [item.dll_error_chips for item in measurements], None),
|
||||
("Carrier frequency (Hz)", [item.carrier_frequency_hz for item in measurements], None),
|
||||
("C/N0 (dB-Hz)", [np.nan if item.cn0_db_hz is None else item.cn0_db_hz for item in measurements], None),
|
||||
("Lock metrics", [item.carrier_lock_metric for item in measurements], [item.code_lock_metric for item in measurements]),
|
||||
]
|
||||
for axis, item in zip(axes.flat, series):
|
||||
title, first, *additional = item
|
||||
axis.plot(time_s, first)
|
||||
for second in additional:
|
||||
if second is not None:
|
||||
axis.plot(time_s, second)
|
||||
axis.set_title(title)
|
||||
axes[-1, 0].set_xlabel("Time (s)")
|
||||
axes[-1, 1].set_xlabel("Time (s)")
|
||||
figure.tight_layout()
|
||||
return figure
|
||||
|
||||
|
||||
def fixed_error(reference: np.ndarray, fixed: np.ndarray, overflow_indices: np.ndarray | None = None) -> Any:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
figure, axis = plt.subplots()
|
||||
axis.plot(np.asarray(fixed) - np.asarray(reference), label="fixed - float")
|
||||
if overflow_indices is not None and overflow_indices.size:
|
||||
axis.scatter(overflow_indices, np.zeros_like(overflow_indices), marker="x", label="overflow")
|
||||
axis.set(xlabel="Sample", ylabel="Error", title="Fixed-point error")
|
||||
axis.legend()
|
||||
return figure
|
||||
@@ -0,0 +1,32 @@
|
||||
"""Diagnostics export without a plotting dependency."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import json
|
||||
from dataclasses import asdict, is_dataclass
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def export_json(path: str | Path, value: Any) -> None:
|
||||
payload = asdict(value) if is_dataclass(value) else value
|
||||
Path(path).write_text(json.dumps(payload, indent=2, default=str), encoding="utf-8")
|
||||
|
||||
|
||||
def export_tracking_csv(path: str | Path, measurements: list[Any]) -> None:
|
||||
if not measurements:
|
||||
Path(path).write_text("", encoding="utf-8")
|
||||
return
|
||||
rows = [asdict(item) for item in measurements]
|
||||
with Path(path).open("w", newline="", encoding="utf-8") as stream:
|
||||
writer = csv.DictWriter(stream, fieldnames=list(rows[0]))
|
||||
writer.writeheader()
|
||||
writer.writerows(rows)
|
||||
|
||||
|
||||
def export_numpy(path: str | Path, **arrays: np.ndarray) -> None:
|
||||
np.savez(Path(path), **arrays)
|
||||
|
||||
@@ -0,0 +1,7 @@
|
||||
"""Bit-true signed fixed-point models."""
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .fixed_fft import FixedRadix2Fft
|
||||
from .quantizer import Quantizer
|
||||
|
||||
__all__ = ["FixedFormat", "OverflowMode", "RoundingMode", "Quantizer", "FixedRadix2Fft"]
|
||||
@@ -0,0 +1,65 @@
|
||||
"""Numerical comparison and word-length sweep utilities."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, asdict
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .fixed_acquisition import FixedPcpsAcquisition
|
||||
from .formats import FixedFormat
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class NumericComparison:
|
||||
absolute_error_mean: float
|
||||
relative_error_mean: float
|
||||
rms_error: float
|
||||
peak_error: float
|
||||
acquisition_decision_agreement: bool | None = None
|
||||
doppler_bin_agreement: bool | None = None
|
||||
code_phase_agreement: bool | None = None
|
||||
overflow_count: int = 0
|
||||
|
||||
def to_dict(self) -> dict[str, float | bool | None]:
|
||||
return asdict(self)
|
||||
|
||||
|
||||
def compare_arrays(reference: npt.ArrayLike, fixed: npt.ArrayLike, overflow_count: int = 0) -> NumericComparison:
|
||||
reference_array = np.asarray(reference, dtype=np.float64)
|
||||
fixed_array = np.asarray(fixed, dtype=np.float64)
|
||||
error = fixed_array - reference_array
|
||||
absolute = np.abs(error)
|
||||
relative = absolute / np.maximum(np.abs(reference_array), np.finfo(float).eps)
|
||||
return NumericComparison(float(np.mean(absolute)), float(np.mean(relative)), float(np.sqrt(np.mean(error**2))), float(np.max(absolute)), overflow_count=overflow_count)
|
||||
|
||||
|
||||
def compare_acquisition(reference_power: npt.ArrayLike, fixed_grid, threshold: float) -> NumericComparison:
|
||||
comparison = compare_arrays(reference_power, fixed_grid.power, fixed_grid.overflow_count)
|
||||
reference_peak = np.unravel_index(int(np.argmax(reference_power)), np.asarray(reference_power).shape)
|
||||
reference_metric = float(np.max(reference_power) / max(np.mean(reference_power), np.finfo(float).tiny))
|
||||
fixed_metric = float(np.max(fixed_grid.power) / max(np.mean(fixed_grid.power), np.finfo(float).tiny))
|
||||
comparison.acquisition_decision_agreement = bool((reference_metric >= threshold) == (fixed_metric >= threshold))
|
||||
comparison.doppler_bin_agreement = bool(reference_peak[0] == fixed_grid.peak_doppler_index)
|
||||
comparison.code_phase_agreement = bool(reference_peak[1] == fixed_grid.peak_code_index)
|
||||
return comparison
|
||||
|
||||
|
||||
def sweep_acquisition_word_lengths(samples: npt.ArrayLike, replica: npt.ArrayLike, sampling_frequency_hz: float, doppler_bins_hz: npt.ArrayLike, bit_widths: list[int]) -> list[dict[str, float | int | bool | None]]:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
code = np.asarray(replica)
|
||||
time_s = np.arange(array.size) / sampling_frequency_hz
|
||||
bins = np.asarray(doppler_bins_hz)
|
||||
reference = np.asarray([
|
||||
np.abs(np.fft.ifft(np.fft.fft(array * np.exp(-2j * np.pi * doppler * time_s)) * np.conjugate(np.fft.fft(code)))) ** 2
|
||||
for doppler in bins
|
||||
])
|
||||
reports = []
|
||||
for bits in bit_widths:
|
||||
model = FixedPcpsAcquisition(FixedFormat(bits, bits - 2), FixedFormat(bits + 8, bits), FixedFormat(48, 16, signed=False))
|
||||
grid = model.search(array, code, sampling_frequency_hz, bins)
|
||||
report = compare_acquisition(reference, grid, 2.5).to_dict()
|
||||
report["input_bits"] = bits
|
||||
reports.append(report)
|
||||
return reports
|
||||
@@ -0,0 +1,69 @@
|
||||
"""Quantized PCPS boundary model with explicit FFT scaling."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .quantizer import Quantizer, dequantize_complex, quantize_complex
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedAcquisitionGrid:
|
||||
power: npt.NDArray[np.float64]
|
||||
peak_doppler_index: int
|
||||
peak_code_index: int
|
||||
overflow_count: int
|
||||
|
||||
|
||||
class FixedPcpsAcquisition:
|
||||
"""Bit-true I/O/product/accumulator model around a normalized FFT primitive.
|
||||
|
||||
The transform primitive models an FFT IP configured with deterministic 1/N
|
||||
scaling. Vendor-specific butterfly truncation is intentionally delegated to
|
||||
the cycle/RTL FFT model selected during implementation.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
sample_format: FixedFormat = FixedFormat(10, 8),
|
||||
fft_format: FixedFormat = FixedFormat(20, 12),
|
||||
power_format: FixedFormat = FixedFormat(48, 16, signed=False),
|
||||
overflow: OverflowMode | str = OverflowMode.SATURATE,
|
||||
rounding: RoundingMode | str = RoundingMode.NEAREST,
|
||||
) -> None:
|
||||
self.input_q = Quantizer(sample_format, overflow, rounding)
|
||||
self.fft_q = Quantizer(fft_format, overflow, rounding)
|
||||
self.power_q = Quantizer(power_format, overflow, rounding)
|
||||
self.sample_format = sample_format
|
||||
self.fft_format = fft_format
|
||||
|
||||
def search(self, samples: npt.ArrayLike, local_replica: npt.ArrayLike, sampling_frequency_hz: float, doppler_bins_hz: npt.ArrayLike) -> FixedAcquisitionGrid:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
replica = np.asarray(local_replica, dtype=np.float64)
|
||||
if array.size != replica.size:
|
||||
raise ValueError("fixed acquisition samples and replica must have equal length")
|
||||
in_r, in_i = quantize_complex(array, self.input_q)
|
||||
quantized_samples = dequantize_complex(in_r, in_i, self.sample_format)
|
||||
code_i = self.input_q.quantize(replica)
|
||||
code = self.input_q.dequantize(code_i)
|
||||
code_spectrum = np.conjugate(np.fft.fft(code) / array.size)
|
||||
time_s = np.arange(array.size) / sampling_frequency_hz
|
||||
bins = np.asarray(doppler_bins_hz, dtype=np.float64)
|
||||
power = np.empty((bins.size, array.size), dtype=np.float64)
|
||||
for index, doppler_hz in enumerate(bins):
|
||||
mixed = quantized_samples * np.exp(-2j * np.pi * doppler_hz * time_s)
|
||||
spectrum = np.fft.fft(mixed) / array.size
|
||||
sr, si = quantize_complex(spectrum, self.fft_q)
|
||||
spectrum_q = dequantize_complex(sr, si, self.fft_format)
|
||||
product = spectrum_q * code_spectrum
|
||||
pr, pi = quantize_complex(product, self.fft_q)
|
||||
correlation = np.fft.ifft(dequantize_complex(pr, pi, self.fft_format)) * array.size
|
||||
power[index] = self.power_q.dequantize(self.power_q.quantize(np.abs(correlation) ** 2))
|
||||
peak = np.unravel_index(int(np.argmax(power)), power.shape)
|
||||
overflow_count = self.input_q.statistics.overflow_count + self.fft_q.statistics.overflow_count + self.power_q.statistics.overflow_count
|
||||
return FixedAcquisitionGrid(power, int(peak[0]), int(peak[1]), overflow_count)
|
||||
|
||||
@@ -0,0 +1,62 @@
|
||||
"""Integer mixer and correlator accumulation model."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .fixed_nco import FixedNco
|
||||
from .fixed_mixer import complex_mix_down
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .quantizer import Quantizer, quantize_complex
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedCorrelationResult:
|
||||
prompt_real: int
|
||||
prompt_imag: int
|
||||
fractional_bits: int
|
||||
overflow_count: int
|
||||
|
||||
@property
|
||||
def complex_value(self) -> complex:
|
||||
return complex(self.prompt_real, self.prompt_imag) / (1 << self.fractional_bits)
|
||||
|
||||
|
||||
class FixedCorrelator:
|
||||
def __init__(
|
||||
self,
|
||||
sampling_frequency_hz: float,
|
||||
carrier_frequency_hz: float,
|
||||
sample_format: FixedFormat = FixedFormat(8, 7),
|
||||
nco_format: FixedFormat = FixedFormat(16, 15),
|
||||
accumulator_format: FixedFormat = FixedFormat(40, 22),
|
||||
overflow: OverflowMode | str = OverflowMode.SATURATE,
|
||||
rounding: RoundingMode | str = RoundingMode.NEAREST,
|
||||
) -> None:
|
||||
self.sample_quantizer = Quantizer(sample_format, overflow, rounding)
|
||||
self.nco = FixedNco(sampling_frequency_hz, carrier_frequency_hz, output_format=nco_format)
|
||||
self.accumulator_quantizer = Quantizer(accumulator_format, overflow, rounding)
|
||||
self.sample_format = sample_format
|
||||
self.nco_format = nco_format
|
||||
self.accumulator_format = accumulator_format
|
||||
|
||||
def correlate_prompt(self, samples: npt.ArrayLike, replica: npt.ArrayLike) -> FixedCorrelationResult:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
local = np.asarray(replica, dtype=np.int64)
|
||||
if array.size != local.size or not set(np.unique(local).tolist()) <= {-1, 1}:
|
||||
raise ValueError("replica must be a same-length +/-1 sequence")
|
||||
sample_real, sample_imag = quantize_complex(array, self.sample_quantizer)
|
||||
cosine, sine = self.nco.generate(array.size)
|
||||
mixer_format = FixedFormat(
|
||||
min(62, self.sample_format.total_bits + self.nco_format.total_bits + 1),
|
||||
self.sample_format.fractional_bits + self.nco_format.fractional_bits,
|
||||
)
|
||||
mixed = complex_mix_down(sample_real, sample_imag, cosine, sine, self.sample_format, self.nco_format, mixer_format)
|
||||
real_sum = int(np.sum(mixed.real * local, dtype=object))
|
||||
imag_sum = int(np.sum(mixed.imag * local, dtype=object))
|
||||
real_q = int(self.accumulator_quantizer.requantize_integer([real_sum], mixer_format.fractional_bits)[0])
|
||||
imag_q = int(self.accumulator_quantizer.requantize_integer([imag_sum], mixer_format.fractional_bits)[0])
|
||||
return FixedCorrelationResult(real_q, imag_q, self.accumulator_format.fractional_bits, self.accumulator_quantizer.statistics.overflow_count)
|
||||
@@ -0,0 +1,100 @@
|
||||
"""Radix-2 stage-quantized FFT/IFFT oracle for RTL verification."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .quantizer import Quantizer, dequantize_complex, quantize_complex
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedFftResult:
|
||||
real: npt.NDArray[np.int64]
|
||||
imag: npt.NDArray[np.int64]
|
||||
fixed_format: FixedFormat
|
||||
overflow_count: int
|
||||
|
||||
def dequantize(self) -> npt.NDArray[np.complex128]:
|
||||
return dequantize_complex(self.real, self.imag, self.fixed_format)
|
||||
|
||||
|
||||
class FixedRadix2Fft:
|
||||
"""Decimation-in-time FFT with one-bit scaling at every stage.
|
||||
|
||||
Forward output equals ``numpy.fft.fft(x) / N`` apart from declared
|
||||
coefficient and butterfly quantization. The inverse uses the conjugation
|
||||
identity, so it has NumPy IFFT scaling.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
data_format: FixedFormat = FixedFormat(20, 15),
|
||||
twiddle_format: FixedFormat = FixedFormat(16, 15),
|
||||
overflow: OverflowMode | str = OverflowMode.SATURATE,
|
||||
rounding: RoundingMode | str = RoundingMode.CONVERGENT,
|
||||
) -> None:
|
||||
self.data_format = data_format
|
||||
self.twiddle_format = twiddle_format
|
||||
self.overflow = overflow
|
||||
self.rounding = rounding
|
||||
self.data_quantizer = Quantizer(data_format, overflow, rounding)
|
||||
self.twiddle_quantizer = Quantizer(twiddle_format, overflow, rounding)
|
||||
|
||||
@staticmethod
|
||||
def _bit_reversed_indices(length: int) -> npt.NDArray[np.int64]:
|
||||
bits = length.bit_length() - 1
|
||||
return np.asarray([int(f"{index:0{bits}b}"[::-1], 2) for index in range(length)], dtype=np.int64)
|
||||
|
||||
def _forward_integers(self, real: npt.NDArray[np.int64], imag: npt.NDArray[np.int64]) -> FixedFftResult:
|
||||
length = real.size
|
||||
if length == 0 or length & (length - 1):
|
||||
raise ValueError("fixed radix-2 FFT length must be a nonzero power of two")
|
||||
order = self._bit_reversed_indices(length)
|
||||
out_real = real[order].astype(np.int64, copy=True)
|
||||
out_imag = imag[order].astype(np.int64, copy=True)
|
||||
butterfly_q = Quantizer(self.data_format, self.overflow, self.rounding)
|
||||
size = 2
|
||||
while size <= length:
|
||||
half = size // 2
|
||||
angle = -2.0 * np.pi * np.arange(half) / size
|
||||
twiddle_real = self.twiddle_quantizer.quantize(np.cos(angle))
|
||||
twiddle_imag = self.twiddle_quantizer.quantize(np.sin(angle))
|
||||
product_fraction = self.data_format.fractional_bits + self.twiddle_format.fractional_bits
|
||||
for first in range(0, length, size):
|
||||
upper = slice(first, first + half)
|
||||
lower = slice(first + half, first + size)
|
||||
u_real = out_real[upper].copy()
|
||||
u_imag = out_imag[upper].copy()
|
||||
lower_real = out_real[lower].copy()
|
||||
lower_imag = out_imag[lower].copy()
|
||||
t_real = butterfly_q.requantize_integer(lower_real * twiddle_real - lower_imag * twiddle_imag, product_fraction)
|
||||
t_imag = butterfly_q.requantize_integer(lower_real * twiddle_imag + lower_imag * twiddle_real, product_fraction)
|
||||
# One-bit scale per stage bounds growth; source fractional+1
|
||||
# makes the integer requantizer perform a rounded divide by two.
|
||||
out_real[upper] = butterfly_q.requantize_integer(u_real + t_real, self.data_format.fractional_bits + 1)
|
||||
out_imag[upper] = butterfly_q.requantize_integer(u_imag + t_imag, self.data_format.fractional_bits + 1)
|
||||
out_real[lower] = butterfly_q.requantize_integer(u_real - t_real, self.data_format.fractional_bits + 1)
|
||||
out_imag[lower] = butterfly_q.requantize_integer(u_imag - t_imag, self.data_format.fractional_bits + 1)
|
||||
size *= 2
|
||||
# Q1.(W-1) maps the mathematical +1 twiddle to its maximum positive
|
||||
# code at table-build time; that coefficient approximation is not a
|
||||
# runtime datapath overflow.
|
||||
overflow_count = (
|
||||
self.data_quantizer.statistics.overflow_count
|
||||
+ butterfly_q.statistics.overflow_count
|
||||
)
|
||||
return FixedFftResult(out_real, out_imag, self.data_format, overflow_count)
|
||||
|
||||
def forward(self, values: npt.ArrayLike) -> FixedFftResult:
|
||||
real, imag = quantize_complex(values, self.data_quantizer)
|
||||
return self._forward_integers(real, imag)
|
||||
|
||||
def inverse(self, values: npt.ArrayLike) -> FixedFftResult:
|
||||
array = np.asarray(values, dtype=np.complex128)
|
||||
real, imag = quantize_complex(np.conjugate(array), self.data_quantizer)
|
||||
forward = self._forward_integers(real, imag)
|
||||
return FixedFftResult(forward.real, -forward.imag, forward.fixed_format, forward.overflow_count)
|
||||
@@ -0,0 +1,59 @@
|
||||
"""Bit-true tracking loop filter states and coefficients."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .quantizer import Quantizer
|
||||
from galileo_rx.tracking.loop_filters import second_order_natural_frequency_rad_s
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedSecondOrderLoopFilter:
|
||||
noise_bandwidth_hz: float
|
||||
integration_time_s: float
|
||||
damping_ratio: float = 0.707
|
||||
error_format: FixedFormat = FixedFormat(24, 18)
|
||||
coefficient_format: FixedFormat = FixedFormat(24, 16)
|
||||
state_format: FixedFormat = FixedFormat(32, 20)
|
||||
overflow: OverflowMode | str = OverflowMode.SATURATE
|
||||
rounding: RoundingMode | str = RoundingMode.CONVERGENT
|
||||
integrator_integer: int = 0
|
||||
overflow_count: int = 0
|
||||
_error_q: Quantizer = field(init=False, repr=False)
|
||||
_coefficient_q: Quantizer = field(init=False, repr=False)
|
||||
_state_q: Quantizer = field(init=False, repr=False)
|
||||
proportional_gain_integer: int = field(init=False)
|
||||
integral_step_gain_integer: int = field(init=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
if self.integration_time_s <= 0:
|
||||
raise ValueError("integration_time_s must be positive")
|
||||
self._error_q = Quantizer(self.error_format, self.overflow, self.rounding)
|
||||
self._coefficient_q = Quantizer(self.coefficient_format, self.overflow, self.rounding)
|
||||
self._state_q = Quantizer(self.state_format, self.overflow, self.rounding)
|
||||
wn = second_order_natural_frequency_rad_s(self.noise_bandwidth_hz, self.damping_ratio)
|
||||
self.proportional_gain_integer = int(self._coefficient_q.quantize([2 * self.damping_ratio * wn])[0])
|
||||
self.integral_step_gain_integer = int(self._coefficient_q.quantize([wn**2 * self.integration_time_s])[0])
|
||||
|
||||
def update(self, error: float) -> float:
|
||||
error_integer = int(self._error_q.quantize([error])[0])
|
||||
product_fraction = self.error_format.fractional_bits + self.coefficient_format.fractional_bits
|
||||
integral_delta = int(self._state_q.requantize_integer([error_integer * self.integral_step_gain_integer], product_fraction)[0])
|
||||
self.integrator_integer = int(self._state_q.requantize_integer(
|
||||
[self.integrator_integer + integral_delta], self.state_format.fractional_bits
|
||||
)[0])
|
||||
proportional = int(self._state_q.requantize_integer(
|
||||
[error_integer * self.proportional_gain_integer], product_fraction
|
||||
)[0])
|
||||
output_integer = int(self._state_q.requantize_integer(
|
||||
[proportional + self.integrator_integer], self.state_format.fractional_bits
|
||||
)[0])
|
||||
self.overflow_count = (
|
||||
self._error_q.statistics.overflow_count
|
||||
+ self._coefficient_q.statistics.overflow_count
|
||||
+ self._state_q.statistics.overflow_count
|
||||
)
|
||||
return output_integer / self.state_format.scale
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
"""Standalone bit-true complex mixer primitive."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
from .quantizer import Quantizer
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedComplex:
|
||||
real: npt.NDArray[np.int64]
|
||||
imag: npt.NDArray[np.int64]
|
||||
fixed_format: FixedFormat
|
||||
overflow_count: int
|
||||
|
||||
def dequantize(self) -> npt.NDArray[np.complex128]:
|
||||
return (self.real.astype(np.float64) + 1j * self.imag.astype(np.float64)) / self.fixed_format.scale
|
||||
|
||||
|
||||
def complex_mix_down(
|
||||
sample_real: npt.ArrayLike,
|
||||
sample_imag: npt.ArrayLike,
|
||||
carrier_cos: npt.ArrayLike,
|
||||
carrier_sin: npt.ArrayLike,
|
||||
sample_format: FixedFormat,
|
||||
carrier_format: FixedFormat,
|
||||
output_format: FixedFormat,
|
||||
overflow: OverflowMode | str = OverflowMode.SATURATE,
|
||||
rounding: RoundingMode | str = RoundingMode.NEAREST,
|
||||
) -> FixedComplex:
|
||||
"""Compute ``(I+jQ)*(cos-j*sin)`` from integer operands."""
|
||||
i = np.asarray(sample_real, dtype=np.int64)
|
||||
q = np.asarray(sample_imag, dtype=np.int64)
|
||||
c = np.asarray(carrier_cos, dtype=np.int64)
|
||||
s = np.asarray(carrier_sin, dtype=np.int64)
|
||||
if not (i.shape == q.shape == c.shape == s.shape):
|
||||
raise ValueError("fixed mixer operands must have equal shapes")
|
||||
product_fractional_bits = sample_format.fractional_bits + carrier_format.fractional_bits
|
||||
quantizer = Quantizer(output_format, overflow, rounding)
|
||||
real = quantizer.requantize_integer(i * c + q * s, product_fractional_bits)
|
||||
imag = quantizer.requantize_integer(q * c - i * s, product_fractional_bits)
|
||||
return FixedComplex(real, imag, output_format, quantizer.statistics.overflow_count)
|
||||
|
||||
@@ -0,0 +1,44 @@
|
||||
"""Bit-true phase word and quantized carrier LUT model."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .formats import FixedFormat
|
||||
from .quantizer import Quantizer
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FixedNco:
|
||||
sampling_frequency_hz: float
|
||||
frequency_hz: float
|
||||
phase_width: int = 32
|
||||
lut_bits: int = 10
|
||||
output_format: FixedFormat = FixedFormat(16, 15)
|
||||
phase_word: int = 0
|
||||
_cos_lut: npt.NDArray[np.int64] = field(init=False, repr=False)
|
||||
_sin_lut: npt.NDArray[np.int64] = field(init=False, repr=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
if self.sampling_frequency_hz <= 0 or self.lut_bits > self.phase_width:
|
||||
raise ValueError("invalid fixed NCO configuration")
|
||||
phase = 2 * np.pi * np.arange(1 << self.lut_bits) / (1 << self.lut_bits)
|
||||
quantizer = Quantizer(self.output_format)
|
||||
self._cos_lut = quantizer.quantize(np.cos(phase))
|
||||
self._sin_lut = quantizer.quantize(np.sin(phase))
|
||||
|
||||
@property
|
||||
def frequency_control_word(self) -> int:
|
||||
return int(round(self.frequency_hz / self.sampling_frequency_hz * (1 << self.phase_width)))
|
||||
|
||||
def generate(self, sample_count: int) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.int64]]:
|
||||
modulus = 1 << self.phase_width
|
||||
fcw = self.frequency_control_word
|
||||
words = np.asarray([(self.phase_word + index * fcw) % modulus for index in range(sample_count)], dtype=np.uint64)
|
||||
address = words >> (self.phase_width - self.lut_bits)
|
||||
self.phase_word = (self.phase_word + sample_count * fcw) % modulus
|
||||
return self._cos_lut[address], self._sin_lut[address]
|
||||
|
||||
@@ -0,0 +1,62 @@
|
||||
"""Fixed-point format declarations."""
|
||||
|
||||
from dataclasses import dataclass
|
||||
from enum import Enum
|
||||
|
||||
|
||||
class OverflowMode(str, Enum):
|
||||
SATURATE = "saturate"
|
||||
WRAP = "wrap"
|
||||
|
||||
|
||||
class RoundingMode(str, Enum):
|
||||
NEAREST = "nearest"
|
||||
TRUNCATE = "truncate"
|
||||
CONVERGENT = "convergent"
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class FixedFormat:
|
||||
total_bits: int
|
||||
fractional_bits: int
|
||||
signed: bool = True
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
if self.total_bits < 1 or self.fractional_bits < 0 or self.fractional_bits >= self.total_bits:
|
||||
raise ValueError("fixed format requires total_bits >= 1 and 0 <= fractional_bits < total_bits")
|
||||
|
||||
@property
|
||||
def minimum_integer(self) -> int:
|
||||
return -(1 << (self.total_bits - 1)) if self.signed else 0
|
||||
|
||||
@property
|
||||
def maximum_integer(self) -> int:
|
||||
return (1 << (self.total_bits - int(self.signed))) - 1
|
||||
|
||||
@property
|
||||
def scale(self) -> int:
|
||||
return 1 << self.fractional_bits
|
||||
|
||||
@property
|
||||
def lsb(self) -> float:
|
||||
return 1.0 / self.scale
|
||||
|
||||
|
||||
STARTING_WORD_LENGTHS: dict[str, FixedFormat] = {
|
||||
"input_samples": FixedFormat(8, 7),
|
||||
"mixer_sine_cosine": FixedFormat(16, 15),
|
||||
"mixer_output": FixedFormat(18, 15),
|
||||
"fft_input": FixedFormat(18, 15),
|
||||
"fft_twiddle": FixedFormat(16, 15),
|
||||
"fft_stage_output": FixedFormat(20, 15),
|
||||
"correlation_product": FixedFormat(36, 30),
|
||||
"ifft_output": FixedFormat(28, 15),
|
||||
"coherent_accumulator": FixedFormat(40, 15),
|
||||
"noncoherent_accumulator": FixedFormat(56, 20, signed=False),
|
||||
"carrier_phase_accumulator": FixedFormat(32, 0, signed=False),
|
||||
"code_phase_accumulator": FixedFormat(48, 32, signed=False),
|
||||
"loop_filter_state": FixedFormat(32, 20),
|
||||
"frequency_control_word": FixedFormat(32, 0),
|
||||
"discriminator_output": FixedFormat(24, 18),
|
||||
}
|
||||
|
||||
@@ -0,0 +1,71 @@
|
||||
"""Two's-complement quantization with observable overflow and error."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from .formats import FixedFormat, OverflowMode, RoundingMode
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class QuantizationStatistics:
|
||||
sample_count: int = 0
|
||||
overflow_count: int = 0
|
||||
squared_error_sum: float = 0.0
|
||||
peak_error: float = 0.0
|
||||
|
||||
@property
|
||||
def rms_error(self) -> float:
|
||||
return float(np.sqrt(self.squared_error_sum / self.sample_count)) if self.sample_count else 0.0
|
||||
|
||||
|
||||
class Quantizer:
|
||||
def __init__(self, fixed_format: FixedFormat, overflow: OverflowMode | str = OverflowMode.SATURATE, rounding: RoundingMode | str = RoundingMode.NEAREST) -> None:
|
||||
self.format = fixed_format
|
||||
self.overflow = OverflowMode(overflow)
|
||||
self.rounding = RoundingMode(rounding)
|
||||
self.statistics = QuantizationStatistics()
|
||||
|
||||
def _round(self, scaled: np.ndarray) -> np.ndarray:
|
||||
if self.rounding == RoundingMode.NEAREST:
|
||||
return np.floor(scaled + 0.5)
|
||||
if self.rounding == RoundingMode.TRUNCATE:
|
||||
return np.trunc(scaled)
|
||||
return np.rint(scaled) # IEEE ties-to-even is convergent rounding.
|
||||
|
||||
def quantize(self, values: npt.ArrayLike) -> npt.NDArray[np.int64]:
|
||||
array = np.asarray(values, dtype=np.float64)
|
||||
rounded = self._round(array * self.format.scale).astype(np.int64)
|
||||
overflow_mask = (rounded < self.format.minimum_integer) | (rounded > self.format.maximum_integer)
|
||||
self.statistics.overflow_count += int(np.count_nonzero(overflow_mask))
|
||||
if self.overflow == OverflowMode.SATURATE:
|
||||
quantized = np.clip(rounded, self.format.minimum_integer, self.format.maximum_integer)
|
||||
else:
|
||||
modulus = 1 << self.format.total_bits
|
||||
quantized = (rounded - self.format.minimum_integer) % modulus + self.format.minimum_integer
|
||||
restored = quantized.astype(np.float64) / self.format.scale
|
||||
error = restored - array
|
||||
self.statistics.sample_count += array.size
|
||||
self.statistics.squared_error_sum += float(np.sum(error**2))
|
||||
self.statistics.peak_error = max(self.statistics.peak_error, float(np.max(np.abs(error))) if error.size else 0.0)
|
||||
return quantized.astype(np.int64)
|
||||
|
||||
def dequantize(self, integers: npt.ArrayLike) -> npt.NDArray[np.float64]:
|
||||
return np.asarray(integers, dtype=np.float64) / self.format.scale
|
||||
|
||||
def requantize_integer(self, integers: npt.ArrayLike, source_fractional_bits: int) -> npt.NDArray[np.int64]:
|
||||
values = np.asarray(integers, dtype=np.float64) / float(1 << source_fractional_bits)
|
||||
return self.quantize(values)
|
||||
|
||||
|
||||
def quantize_complex(values: npt.ArrayLike, quantizer: Quantizer) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.int64]]:
|
||||
array = np.asarray(values, dtype=np.complex128)
|
||||
return quantizer.quantize(array.real), quantizer.quantize(array.imag)
|
||||
|
||||
|
||||
def dequantize_complex(real: npt.ArrayLike, imag: npt.ArrayLike, fixed_format: FixedFormat) -> npt.NDArray[np.complex128]:
|
||||
return (np.asarray(real, dtype=np.float64) + 1j * np.asarray(imag, dtype=np.float64)) / fixed_format.scale
|
||||
|
||||
@@ -0,0 +1,7 @@
|
||||
"""Streaming digital front-end blocks."""
|
||||
|
||||
from .frontend_chain import FrontendChain
|
||||
from .mixer import PhaseAccumulatorNco
|
||||
|
||||
__all__ = ["FrontendChain", "PhaseAccumulatorNco"]
|
||||
|
||||
@@ -0,0 +1,24 @@
|
||||
"""Block RMS normalization with bounded gain."""
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
class BlockAgc:
|
||||
def __init__(self, target_rms: float = 0.5, smoothing: float = 0.9, min_gain: float = 1e-3, max_gain: float = 1e3) -> None:
|
||||
if target_rms <= 0 or not 0 <= smoothing < 1:
|
||||
raise ValueError("invalid AGC parameters")
|
||||
self.target_rms = target_rms
|
||||
self.smoothing = smoothing
|
||||
self.min_gain = min_gain
|
||||
self.max_gain = max_gain
|
||||
self.gain = 1.0
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic]) -> npt.NDArray[np.complex128]:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
rms = float(np.sqrt(np.mean(np.abs(array) ** 2))) if array.size else 0.0
|
||||
desired = self.target_rms / max(rms, 1e-15)
|
||||
desired = float(np.clip(desired, self.min_gain, self.max_gain))
|
||||
self.gain = self.smoothing * self.gain + (1.0 - self.smoothing) * desired
|
||||
return array * self.gain
|
||||
|
||||
@@ -0,0 +1,25 @@
|
||||
"""First-order streaming DC blocker."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
class DcBlocker:
|
||||
def __init__(self, alpha: float = 0.999) -> None:
|
||||
if not 0 <= alpha < 1:
|
||||
raise ValueError("DC blocker alpha must be in [0, 1)")
|
||||
self.alpha = alpha
|
||||
self.previous_input = 0j
|
||||
self.previous_output = 0j
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic]) -> npt.NDArray[np.complex128]:
|
||||
output = np.empty(samples.size, dtype=np.complex128)
|
||||
for index, value in enumerate(np.asarray(samples, dtype=np.complex128)):
|
||||
current = value - self.previous_input + self.alpha * self.previous_output
|
||||
output[index] = current
|
||||
self.previous_input = value
|
||||
self.previous_output = current
|
||||
return output
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
"""Stateful FIR and integer decimation."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
from scipy.signal import lfilter
|
||||
|
||||
|
||||
class FirDecimator:
|
||||
def __init__(self, taps: list[float] | npt.NDArray[np.float64] | None = None, factor: int = 1) -> None:
|
||||
if factor < 1:
|
||||
raise ValueError("decimation factor must be >= 1")
|
||||
self.taps = np.asarray(taps if taps is not None and len(taps) else [1.0], dtype=np.float64)
|
||||
self.factor = factor
|
||||
self._state = np.zeros(max(0, self.taps.size - 1), dtype=np.complex128)
|
||||
self._decimation_phase = 0
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic]) -> npt.NDArray[np.complex128]:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
if self.taps.size == 1:
|
||||
filtered = array * self.taps[0]
|
||||
else:
|
||||
filtered, self._state = lfilter(self.taps, [1.0], array, zi=self._state)
|
||||
start = (-self._decimation_phase) % self.factor
|
||||
result = filtered[start::self.factor]
|
||||
self._decimation_phase = (self._decimation_phase + array.size) % self.factor
|
||||
return result
|
||||
|
||||
@@ -0,0 +1,42 @@
|
||||
"""Composable streaming RF-to-complex-baseband chain."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.config import FrontendConfig, RfConfig
|
||||
from .agc import BlockAgc
|
||||
from .dc_blocker import DcBlocker
|
||||
from .decimator import FirDecimator
|
||||
from .mixer import ComplexMixer, PhaseAccumulatorNco
|
||||
|
||||
|
||||
class FrontendChain:
|
||||
def __init__(self, rf: RfConfig, config: FrontendConfig) -> None:
|
||||
self.rf = rf
|
||||
self.config = config
|
||||
self.dc = DcBlocker(config.dc_alpha) if config.remove_dc else None
|
||||
self.agc = BlockAgc(config.agc_target_rms) if config.agc_enabled else None
|
||||
self.mixer = ComplexMixer(
|
||||
PhaseAccumulatorNco(rf.sampling_frequency_hz, rf.intermediate_frequency_hz, config.nco_phase_width, config.nco_lut_bits)
|
||||
)
|
||||
self.decimator = FirDecimator(config.fir_taps, config.decimation)
|
||||
|
||||
@property
|
||||
def output_sampling_frequency_hz(self) -> float:
|
||||
return self.rf.sampling_frequency_hz / self.config.decimation
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic]) -> npt.NDArray[np.complex128]:
|
||||
array = np.asarray(samples)
|
||||
if np.issubdtype(array.dtype, np.integer):
|
||||
array = array.astype(np.float64) / float(1 << (self.rf.sample_bit_width - 1))
|
||||
array = array.astype(np.complex128)
|
||||
if self.dc is not None:
|
||||
array = self.dc.process(array)
|
||||
if self.agc is not None:
|
||||
array = self.agc.process(array)
|
||||
if self.rf.input_is_real_if or self.rf.intermediate_frequency_hz != 0:
|
||||
array = self.mixer.process(array, downconvert=True)
|
||||
return self.decimator.process(array)
|
||||
|
||||
@@ -0,0 +1,66 @@
|
||||
"""Floating and LUT phase-accumulator carrier NCO."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class PhaseAccumulatorNco:
|
||||
sampling_frequency_hz: float
|
||||
frequency_hz: float
|
||||
phase_width: int = 32
|
||||
lut_bits: int = 12
|
||||
phase_word: int = 0
|
||||
_lut: npt.NDArray[np.complex128] = field(init=False, repr=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
if self.sampling_frequency_hz <= 0 or not 8 <= self.phase_width <= 64:
|
||||
raise ValueError("invalid NCO sampling frequency or phase width")
|
||||
if not 4 <= self.lut_bits <= min(self.phase_width, 20):
|
||||
raise ValueError("invalid LUT address width")
|
||||
size = 1 << self.lut_bits
|
||||
self._lut = np.exp(2j * np.pi * np.arange(size) / size)
|
||||
|
||||
@property
|
||||
def modulus(self) -> int:
|
||||
return 1 << self.phase_width
|
||||
|
||||
@property
|
||||
def frequency_control_word(self) -> int:
|
||||
return int(round(self.frequency_hz / self.sampling_frequency_hz * self.modulus))
|
||||
|
||||
@property
|
||||
def actual_frequency_hz(self) -> float:
|
||||
return self.frequency_control_word * self.sampling_frequency_hz / self.modulus
|
||||
|
||||
def set_frequency(self, frequency_hz: float) -> None:
|
||||
self.frequency_hz = frequency_hz
|
||||
|
||||
def generate(self, sample_count: int, conjugate: bool = False, use_lut: bool = False) -> npt.NDArray[np.complex128]:
|
||||
if sample_count < 0:
|
||||
raise ValueError("sample_count must be nonnegative")
|
||||
fcw = self.frequency_control_word
|
||||
words = (self.phase_word + np.arange(sample_count, dtype=object) * fcw)
|
||||
if use_lut:
|
||||
shift = self.phase_width - self.lut_bits
|
||||
addresses = np.asarray([(int(word) % self.modulus) >> shift for word in words], dtype=np.int64)
|
||||
values = self._lut[addresses]
|
||||
else:
|
||||
phases = np.asarray([int(word) % self.modulus for word in words], dtype=np.float64) / self.modulus
|
||||
values = np.exp(2j * np.pi * phases)
|
||||
self.phase_word = int((self.phase_word + sample_count * fcw) % self.modulus)
|
||||
return np.conjugate(values) if conjugate else values
|
||||
|
||||
|
||||
class ComplexMixer:
|
||||
def __init__(self, nco: PhaseAccumulatorNco, use_lut: bool = False) -> None:
|
||||
self.nco = nco
|
||||
self.use_lut = use_lut
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic], downconvert: bool = True) -> npt.NDArray[np.complex128]:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
return array * self.nco.generate(array.size, conjugate=downconvert, use_lut=self.use_lut)
|
||||
@@ -0,0 +1,14 @@
|
||||
"""Sample file I/O."""
|
||||
|
||||
from .sample_reader import SampleReader, read_samples
|
||||
from .sample_writer import write_samples
|
||||
from .signalsim_reader import SignalSimReader, SignalSimSampleFormat, SignalSimSampleSpec
|
||||
|
||||
__all__ = [
|
||||
"SampleReader",
|
||||
"SignalSimReader",
|
||||
"SignalSimSampleFormat",
|
||||
"SignalSimSampleSpec",
|
||||
"read_samples",
|
||||
"write_samples",
|
||||
]
|
||||
@@ -0,0 +1,35 @@
|
||||
"""Block-streaming raw sample input."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from pathlib import Path
|
||||
from typing import Iterator
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
def _dtype(sample_format: str, endianness: str) -> np.dtype:
|
||||
mapping = {"complex128": "c16", "float64": "f8", "complex64": "c8", "float32": "f4", "int8": "i1", "int16": "i2"}
|
||||
if sample_format not in mapping:
|
||||
raise ValueError(f"unsupported sample format {sample_format!r}")
|
||||
prefix = "<" if endianness == "little" else ">"
|
||||
return np.dtype(prefix + mapping[sample_format])
|
||||
|
||||
|
||||
def read_samples(path: str | Path, sample_format: str = "complex128", endianness: str = "little", count: int = -1, offset_samples: int = 0) -> npt.NDArray[np.generic]:
|
||||
dtype = _dtype(sample_format, endianness)
|
||||
return np.fromfile(Path(path), dtype=dtype, count=count, offset=offset_samples * dtype.itemsize)
|
||||
|
||||
|
||||
class SampleReader:
|
||||
def __init__(self, path: str | Path, sample_format: str = "complex128", endianness: str = "little") -> None:
|
||||
self.path = Path(path)
|
||||
self.dtype = _dtype(sample_format, endianness)
|
||||
|
||||
def blocks(self, block_samples: int) -> Iterator[npt.NDArray[np.generic]]:
|
||||
if block_samples <= 0:
|
||||
raise ValueError("block_samples must be positive")
|
||||
with self.path.open("rb") as stream:
|
||||
while block := np.fromfile(stream, dtype=self.dtype, count=block_samples):
|
||||
yield block
|
||||
@@ -0,0 +1,13 @@
|
||||
"""Raw sample output."""
|
||||
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
def write_samples(path: str | Path, samples: npt.NDArray[np.generic], endianness: str = "little") -> None:
|
||||
array = np.asarray(samples)
|
||||
byte_order = "<" if endianness == "little" else ">"
|
||||
array.astype(array.dtype.newbyteorder(byte_order), copy=False).tofile(Path(path))
|
||||
|
||||
@@ -0,0 +1,154 @@
|
||||
"""Strict adapters for globsky/SignalSim complex-IQ files.
|
||||
|
||||
All decoded samples are promoted to ``complex128``. ``scale_counts_per_unit``
|
||||
is explicit: use 1.0 for raw ADC counts, 25.0 for IQ8 physical amplitude at
|
||||
AGC gain 1, 3277.0 for IQ16, or the actual final AGC-adjusted scale.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from enum import Enum
|
||||
from pathlib import Path
|
||||
from typing import Iterator, Literal
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
class SignalSimSampleFormat(str, Enum):
|
||||
IQ8 = "IQ8"
|
||||
IQ16 = "IQ16"
|
||||
IQ4 = "IQ4"
|
||||
IQ2 = "IQ2"
|
||||
FLOAT32_IQ = "FLOAT32_IQ"
|
||||
FLOAT64_IQ = "FLOAT64_IQ"
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class SignalSimSampleSpec:
|
||||
sample_format: SignalSimSampleFormat | str = SignalSimSampleFormat.IQ8
|
||||
byte_order: Literal["little", "big"] = "little"
|
||||
iq_order: Literal["IQ", "QI"] = "IQ"
|
||||
scale_counts_per_unit: float = 1.0
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
object.__setattr__(self, "sample_format", SignalSimSampleFormat(self.sample_format))
|
||||
if self.scale_counts_per_unit <= 0.0 or not np.isfinite(self.scale_counts_per_unit):
|
||||
raise ValueError("scale_counts_per_unit must be finite and positive")
|
||||
if self.sample_format == SignalSimSampleFormat.IQ16 and self.byte_order != "little":
|
||||
# This is allowed for hand-crafted/converted data, but is not native SignalSim.
|
||||
pass
|
||||
|
||||
|
||||
def _decode_sign_magnitude(values: np.ndarray, sign_mask: int, magnitude_mask: int) -> np.ndarray:
|
||||
magnitude = (values & magnitude_mask).astype(np.float64)
|
||||
return np.where((values & sign_mask) != 0, -magnitude, magnitude)
|
||||
|
||||
|
||||
class SignalSimReader:
|
||||
"""Random-access and block-streaming SignalSim sample reader."""
|
||||
|
||||
def __init__(self, path: str | Path, spec: SignalSimSampleSpec | None = None, *, memory_map: bool = False) -> None:
|
||||
self.path = Path(path)
|
||||
self.spec = spec or SignalSimSampleSpec()
|
||||
self.memory_map = memory_map
|
||||
if not self.path.is_file():
|
||||
raise FileNotFoundError(self.path)
|
||||
|
||||
@property
|
||||
def sample_count(self) -> int:
|
||||
byte_count = self.path.stat().st_size
|
||||
fmt = self.spec.sample_format
|
||||
divisor = {
|
||||
SignalSimSampleFormat.IQ8: 2,
|
||||
SignalSimSampleFormat.IQ16: 4,
|
||||
SignalSimSampleFormat.IQ4: 1,
|
||||
SignalSimSampleFormat.FLOAT32_IQ: 8,
|
||||
SignalSimSampleFormat.FLOAT64_IQ: 16,
|
||||
}.get(fmt)
|
||||
if fmt == SignalSimSampleFormat.IQ2:
|
||||
return byte_count * 2
|
||||
assert divisor is not None
|
||||
if byte_count % divisor:
|
||||
raise ValueError(f"truncated {fmt.value} file: {byte_count} bytes is not a whole complex sample")
|
||||
return byte_count // divisor
|
||||
|
||||
def _raw_bytes(self) -> npt.NDArray[np.uint8]:
|
||||
if self.memory_map:
|
||||
return np.memmap(self.path, dtype=np.uint8, mode="r")
|
||||
return np.fromfile(self.path, dtype=np.uint8)
|
||||
|
||||
def read(self, count: int = -1, offset_samples: int = 0, *, allow_truncated: bool = False) -> npt.NDArray[np.complex128]:
|
||||
available = self.sample_count
|
||||
if offset_samples < 0 or offset_samples > available:
|
||||
raise ValueError(f"offset_samples must be within 0..{available}")
|
||||
requested = available - offset_samples if count < 0 else count
|
||||
if requested < 0:
|
||||
raise ValueError("count must be -1 or nonnegative")
|
||||
remaining = available - offset_samples
|
||||
if requested > remaining:
|
||||
if not allow_truncated:
|
||||
raise EOFError(f"requested {requested} samples at offset {offset_samples}, only {remaining} remain")
|
||||
requested = remaining
|
||||
if requested == 0:
|
||||
return np.empty(0, dtype=np.complex128)
|
||||
|
||||
fmt = self.spec.sample_format
|
||||
if fmt == SignalSimSampleFormat.IQ2:
|
||||
decoded = self._decode_iq2(offset_samples, requested)
|
||||
elif fmt == SignalSimSampleFormat.IQ4:
|
||||
raw = self._raw_bytes()[offset_samples : offset_samples + requested]
|
||||
i_values = _decode_sign_magnitude(raw >> 4, 0x8, 0x7)
|
||||
q_values = _decode_sign_magnitude(raw & 0xF, 0x8, 0x7)
|
||||
decoded = i_values + 1j * q_values
|
||||
else:
|
||||
scalar_dtype = {
|
||||
SignalSimSampleFormat.IQ8: np.dtype("i1"),
|
||||
SignalSimSampleFormat.IQ16: np.dtype(("<" if self.spec.byte_order == "little" else ">") + "i2"),
|
||||
SignalSimSampleFormat.FLOAT32_IQ: np.dtype(("<" if self.spec.byte_order == "little" else ">") + "f4"),
|
||||
SignalSimSampleFormat.FLOAT64_IQ: np.dtype(("<" if self.spec.byte_order == "little" else ">") + "f8"),
|
||||
}[fmt]
|
||||
scalar_offset = offset_samples * 2
|
||||
if self.memory_map:
|
||||
scalars = np.memmap(self.path, dtype=scalar_dtype, mode="r", offset=scalar_offset * scalar_dtype.itemsize, shape=(requested * 2,))
|
||||
else:
|
||||
scalars = np.fromfile(self.path, dtype=scalar_dtype, count=requested * 2, offset=scalar_offset * scalar_dtype.itemsize)
|
||||
first = scalars[0::2].astype(np.float64)
|
||||
second = scalars[1::2].astype(np.float64)
|
||||
decoded = first + 1j * second
|
||||
if self.spec.iq_order == "QI":
|
||||
decoded = decoded.imag + 1j * decoded.real
|
||||
return np.asarray(decoded / np.float64(self.spec.scale_counts_per_unit), dtype=np.complex128)
|
||||
|
||||
def _decode_iq2(self, offset_samples: int, count: int) -> npt.NDArray[np.complex128]:
|
||||
first_byte = offset_samples // 2
|
||||
last_sample = offset_samples + count
|
||||
last_byte = (last_sample + 1) // 2
|
||||
raw = self._raw_bytes()[first_byte:last_byte]
|
||||
output = np.empty(raw.size * 2, dtype=np.complex128)
|
||||
# Each 2-bit sign/magnitude value has levels 0 and +/-1. Samples 0/1
|
||||
# occupy the low/high nibbles respectively (I then Q within a nibble).
|
||||
for lane, shift in enumerate((0, 4)):
|
||||
nibble = raw >> shift
|
||||
i_values = _decode_sign_magnitude(nibble, 0x2, 0x1)
|
||||
q_values = _decode_sign_magnitude(nibble, 0x8, 0x4) / 4.0
|
||||
output[lane::2] = i_values + 1j * q_values
|
||||
begin = offset_samples & 1
|
||||
return output[begin : begin + count]
|
||||
|
||||
def blocks(self, block_samples: int, *, offset_samples: int = 0, count: int = -1) -> Iterator[npt.NDArray[np.complex128]]:
|
||||
if block_samples <= 0:
|
||||
raise ValueError("block_samples must be positive")
|
||||
remaining = self.sample_count - offset_samples if count < 0 else count
|
||||
cursor = offset_samples
|
||||
while remaining > 0:
|
||||
block_count = min(block_samples, remaining)
|
||||
block = self.read(block_count, cursor, allow_truncated=True)
|
||||
if block.size == 0:
|
||||
break
|
||||
yield block
|
||||
cursor += block.size
|
||||
remaining -= block.size
|
||||
if block.size < block_count:
|
||||
break
|
||||
@@ -0,0 +1,6 @@
|
||||
"""Cycle-oriented interfaces and hardware scheduling estimates."""
|
||||
|
||||
from .cycle_model import AcquisitionArchitectureEstimate, TrackingScheduler
|
||||
|
||||
__all__ = ["AcquisitionArchitectureEstimate", "TrackingScheduler"]
|
||||
|
||||
@@ -0,0 +1,71 @@
|
||||
# Receiver RTL module contract
|
||||
|
||||
Widths below are proposed starting points and must be confirmed by the sweep.
|
||||
Every stream transfers only on `valid && ready`; `last` marks an integration or
|
||||
FFT frame. Configuration registers are stable during a frame. Arithmetic data
|
||||
saturates and raises a sticky overflow flag unless natural NCO wrap is stated.
|
||||
|
||||
| Module | Input → output | Width / binary point | Latency & throughput | State / memory | Reset and verification |
|
||||
|---|---|---|---|---|---|
|
||||
| `sample_input` | ADC/file words → signed I/Q stream | 8–16 integer or Q1.(W-1) | 1, one sample/clock | optional CDC FIFO | clear valid/FIFO; replay exported vectors |
|
||||
| `dc_blocker` | I/Q → high-pass I/Q | 18/Q15 internal | 2–4, one/clock | previous x/y per component | zero state; impulse/DC tests |
|
||||
| `agc_or_scaler` | I/Q, gain → scaled I/Q | gain 18/Q16, out 18/Q15 | 2, one/clock | gain/RMS window | unity gain; float/fixed block comparison |
|
||||
| `ddc_nco` | signed FCW, phase load → cos/sin | phase 32 unsigned, LUT 16/Q15 | 1–3, one/clock | phase accumulator, LUT BRAM | programmed phase; FCW/continuity tests |
|
||||
| `complex_mixer` | sample, cos/sin → baseband | 8×16, out 18/Q15 | 3–5, one/clock | pipeline only | clear valid; tone/sign test |
|
||||
| `fir_decimator` | baseband → filtered/decimated | coeff 18/Q17, accum 36+ | tap dependent | delay line/coeff ROM/phase | clear delay/phase; SciPy impulse comparison |
|
||||
| `sample_buffer` | stream → framed random read | 18-bit I + 18-bit Q | one write/clock, read architecture dependent | `M*N_s` complex words, ping-pong | pointers zero; boundary/wrap tests |
|
||||
| `code_rom` | PRN/component/chip addr → +/-1 | 1 sign bit | 1 | validated code ROM only | no invented defaults; checksum against provider |
|
||||
| `code_nco` | FCW/phase → chip address/fraction | 48, 32 fractional | 1, one/sample | phase and epoch counter | phase load; noninteger-rate test |
|
||||
| `boc_subcarrier_nco` | code fraction → sign | phase fraction MSB | combinational/1 | none or phase state | defined sine-phase polarity test |
|
||||
| `fft_controller` | frame/config → FFT commands | address `ceil(log2 N)` | control dependent | frame/bin/PRN counters | idle reset; transaction assertions |
|
||||
| `streaming_fft` | complex frame → spectrum | 18/Q15 in, 20/Q15 stage, configurable shifts | IP-specific `C_FFT`, initiation one/clock | twiddle ROM, delay/transpose RAM | frame reset only; impulse/tone/random bit vectors |
|
||||
| `complex_multiplier` | spectra → product | 20/Q15 × 16/Q15, 36/Q30 full | 3–6, one/clock | pipeline | clear valid; exhaustive small-width test |
|
||||
| `code_fft_memory` | PRN/address → conjugated FFT | 20/Q15 complex | 1–2 | `P*N*40` bits or banked subset | contents retained; CRC and software comparison |
|
||||
| `correlation_ifft` | product frame → correlation | 20/Q15 → 28/Q15 | IP-specific | FFT RAM/twiddles | scaling-vector test |
|
||||
| `magnitude_squared` | complex → unsigned power | 28/Q15 → 56/Q30, then configured 48/Q16 | 4, one/clock | pipeline | saturation sticky; corner tests |
|
||||
| `noncoherent_accumulator` | power/address → sum | unsigned 56/Q20 | read-modify-write 2–3 | `D*N_s` words per active PRN/lane | clear at search; coherent/noncoherent vectors |
|
||||
| `peak_search` | power stream → top/second/noise | 56-bit power, wide mean accumulator | one scan, one/clock | top values, guard FIFO/circular rescan | clear per map; wraparound guard tests |
|
||||
| `acquisition_controller` | command → result stream | register/control | search latency | PRN/bin/epoch/hypothesis FSM | idle/abort; end-to-end generated vectors |
|
||||
| `tracking_correlator_bank` | wiped sample, five signs → products | 18/Q15 sample, sign code | 2–5, channel/sample throughput parameter | pipeline | valid reset; vector/cycle correlator equality |
|
||||
| `accumulate_and_dump` | tap products → five sums | 40/Q15 complex | one/sample, dump at epoch | ten accumulators/channel | zero on dump/reset; exact integer sums |
|
||||
| `tracking_discriminators` | tap sums → PLL/FLL/DLL | input 40/Q15, output 24/Q18 | LUT/CORDIC dependent | previous prompt per channel | zero history; signs/units/corner quadrants |
|
||||
| `loop_filter` | errors/config → corrections | 32/Q20 state, 24/Q18 coeff | 3–10 per update | integrators/channel | zero or software load; coefficient/step vectors |
|
||||
| `carrier_nco` | correction/FCW → carrier | phase 32 unsigned | one/sample | phase/FCW/channel | natural wrap; continuity and pull-in |
|
||||
| `tracking_channel_scheduler` | sample + state RAM → lane jobs | IDs/addresses | must meet scheduling inequality | active bitmap, round-robin pointer | idle; assertions for deadline/no collision |
|
||||
| `channel_state_ram` | scheduler read/write | implementation record ~512 bits/channel | dual-port 1–2 | all NCO/loop/lock/secondary state | configurable preserve/clear; context-switch test |
|
||||
| `measurement_output` | epoch state → host stream | packed 32-bit words | burst per epoch | small FIFO | clear FIFO; schema/golden CSV comparison |
|
||||
|
||||
`fir_decimator` is bypassable. `ddc_nco` can expose `(x,y,z,valid,ready)` instead
|
||||
of LUT output for a later CORDIC without changing mixer control. FFT RAM banking,
|
||||
scaling schedule, and vendor IP latency are deliberately parameters, not hidden
|
||||
assumptions.
|
||||
|
||||
For side-peak/unambiguous extensions, the acquisition controller accepts a
|
||||
replica-strategy identifier and the tracking bank already supplies VEVL. CBOC
|
||||
adds a modulation generator/combiner behind the same sign-replica interface.
|
||||
|
||||
## External-vector handoff contract
|
||||
|
||||
The SignalSim quick vector is the first external RTL oracle. The generated
|
||||
manifest contains raw IQ8 hex, sample timestamps, search limits, coherent and
|
||||
noncoherent counts, FFT length, carrier/code FCWs, acquired code phase and
|
||||
Doppler, E1B/E1C dump lengths, and E1C secondary phase. Symbolic keys become
|
||||
APB fields once the register address map is frozen.
|
||||
|
||||
The software handoff must program two independent quantities:
|
||||
|
||||
```text
|
||||
carrier_nco_frequency = configured_IF_center + acquired_Doppler
|
||||
code_aiding_Doppler = carrier_nco_frequency - configured_IF_center
|
||||
```
|
||||
|
||||
Using the full carrier NCO frequency for code aiding was observed to drive the
|
||||
code NCO about 1330 Hz fast at 2.046 MHz IF and lose lock after one epoch. RTL
|
||||
and firmware assertions should check that the aiding value lies within the
|
||||
configured Doppler search range, not the sampled IF range.
|
||||
|
||||
The E1B channel uses 16,368-sample (2 ms) pull-in dumps at 8.184 MHz so its
|
||||
BPSK-squared FLL has ±125 Hz unambiguous range. Two dumps are paired for a raw
|
||||
4 ms symbol. E1C uses 32,736-sample dumps plus a 25-state secondary counter.
|
||||
These measured modes must be representable by `MAX_INTEGRATION_SAMPLES`, the
|
||||
channel state RAM record, the scheduler, and measurement output.
|
||||
@@ -0,0 +1,74 @@
|
||||
"""Architecture sizing formulas and lossless scheduling checks."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from math import ceil, log2
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class AcquisitionArchitectureEstimate:
|
||||
sampling_frequency_hz: float
|
||||
coherent_integration_s: float
|
||||
noncoherent_accumulations: int
|
||||
number_of_prns: int
|
||||
number_of_doppler_bins: int
|
||||
fft_length: int
|
||||
clock_frequency_hz: float
|
||||
num_doppler_lanes: int = 1
|
||||
fft_cycles_per_transform: int | None = None
|
||||
|
||||
@property
|
||||
def samples_per_coherent_interval(self) -> int:
|
||||
return int(round(self.sampling_frequency_hz * self.coherent_integration_s))
|
||||
|
||||
@property
|
||||
def fft_operations_per_prn(self) -> int:
|
||||
# One local-code FFT may be precomputed; every bin/accumulation needs a
|
||||
# sample FFT and correlation IFFT.
|
||||
return 2 * self.number_of_doppler_bins * self.noncoherent_accumulations
|
||||
|
||||
@property
|
||||
def sample_buffer_complex_words(self) -> int:
|
||||
return self.samples_per_coherent_interval * self.noncoherent_accumulations
|
||||
|
||||
@property
|
||||
def code_fft_complex_words(self) -> int:
|
||||
return self.number_of_prns * self.fft_length
|
||||
|
||||
@property
|
||||
def complex_pointwise_multiplications(self) -> int:
|
||||
return self.number_of_prns * self.number_of_doppler_bins * self.noncoherent_accumulations * self.fft_length
|
||||
|
||||
@property
|
||||
def acquisition_time_s(self) -> float:
|
||||
fft_cycles = self.fft_cycles_per_transform or (self.fft_length + int(log2(self.fft_length)) + 16)
|
||||
bin_groups = ceil(self.number_of_doppler_bins / self.num_doppler_lanes)
|
||||
cycles = self.number_of_prns * bin_groups * self.noncoherent_accumulations * (2 * fft_cycles + self.fft_length)
|
||||
return cycles / self.clock_frequency_hz
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class TrackingScheduler:
|
||||
sampling_frequency_hz: float
|
||||
clock_frequency_hz: float
|
||||
num_tracking_channels: int
|
||||
num_correlator_lanes: int
|
||||
cycles_per_channel_sample: int
|
||||
|
||||
@property
|
||||
def available_cycles_per_input_sample(self) -> float:
|
||||
return self.clock_frequency_hz / self.sampling_frequency_hz
|
||||
|
||||
@property
|
||||
def required_cycles_per_input_sample(self) -> int:
|
||||
return ceil(self.num_tracking_channels / self.num_correlator_lanes) * self.cycles_per_channel_sample
|
||||
|
||||
@property
|
||||
def schedulable_without_loss(self) -> bool:
|
||||
return self.required_cycles_per_input_sample <= self.available_cycles_per_input_sample
|
||||
|
||||
@property
|
||||
def utilization(self) -> float:
|
||||
return self.required_cycles_per_input_sample / self.available_cycles_per_input_sample
|
||||
|
||||
@@ -0,0 +1,34 @@
|
||||
"""Vendor-neutral valid/ready transaction types."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class StreamTransaction:
|
||||
data: int | complex
|
||||
valid: bool = True
|
||||
ready: bool = True
|
||||
first: bool = False
|
||||
last: bool = False
|
||||
channel_id: int = 0
|
||||
user: int = 0
|
||||
|
||||
@property
|
||||
def transferred(self) -> bool:
|
||||
return self.valid and self.ready
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CorrelatorTransaction:
|
||||
sample_i: int
|
||||
sample_q: int
|
||||
carrier_cos: int
|
||||
carrier_sin: int
|
||||
code_taps: tuple[int, ...]
|
||||
channel_id: int
|
||||
epoch_last: bool
|
||||
valid: bool = True
|
||||
ready: bool = True
|
||||
|
||||
@@ -0,0 +1,76 @@
|
||||
"""Export integer transactions and manifests for HDL testbenches."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import json
|
||||
from pathlib import Path
|
||||
import csv
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.fixed_point.formats import FixedFormat
|
||||
from galileo_rx.fixed_point.quantizer import Quantizer, quantize_complex
|
||||
|
||||
|
||||
def export_complex_hex(
|
||||
directory: str | Path,
|
||||
stem: str,
|
||||
values: npt.ArrayLike,
|
||||
fixed_format: FixedFormat,
|
||||
) -> dict[str, object]:
|
||||
"""Write two's-complement I/Q hex files plus a machine-readable manifest."""
|
||||
output_directory = Path(directory)
|
||||
output_directory.mkdir(parents=True, exist_ok=True)
|
||||
quantizer = Quantizer(fixed_format)
|
||||
real, imag = quantize_complex(values, quantizer)
|
||||
mask = (1 << fixed_format.total_bits) - 1
|
||||
digits = (fixed_format.total_bits + 3) // 4
|
||||
real_path = output_directory / f"{stem}_i.hex"
|
||||
imag_path = output_directory / f"{stem}_q.hex"
|
||||
real_path.write_text("\n".join(f"{int(item) & mask:0{digits}x}" for item in real) + "\n", encoding="ascii")
|
||||
imag_path.write_text("\n".join(f"{int(item) & mask:0{digits}x}" for item in imag) + "\n", encoding="ascii")
|
||||
manifest = {
|
||||
"stem": stem,
|
||||
"sample_count": int(real.size),
|
||||
"total_bits": fixed_format.total_bits,
|
||||
"fractional_bits": fixed_format.fractional_bits,
|
||||
"signed": fixed_format.signed,
|
||||
"encoding": "two_complement_hex",
|
||||
"i_file": real_path.name,
|
||||
"q_file": imag_path.name,
|
||||
"overflow_count": quantizer.statistics.overflow_count,
|
||||
}
|
||||
(output_directory / f"{stem}.json").write_text(json.dumps(manifest, indent=2), encoding="utf-8")
|
||||
return manifest
|
||||
|
||||
|
||||
def export_timed_complex_hex(
|
||||
directory: str | Path,
|
||||
stem: str,
|
||||
values: npt.ArrayLike,
|
||||
fixed_format: FixedFormat,
|
||||
sampling_frequency_hz: float,
|
||||
register_values: dict[str, object],
|
||||
) -> dict[str, object]:
|
||||
"""Export RTL samples, explicit timestamps, and configuration metadata."""
|
||||
if sampling_frequency_hz <= 0.0:
|
||||
raise ValueError("sampling_frequency_hz must be positive")
|
||||
output_directory = Path(directory)
|
||||
manifest = export_complex_hex(output_directory, stem, values, fixed_format)
|
||||
timestamp_path = output_directory / f"{stem}_timestamps.csv"
|
||||
with timestamp_path.open("w", newline="", encoding="ascii") as stream:
|
||||
writer = csv.writer(stream)
|
||||
writer.writerow(["sample_index", "time_s"])
|
||||
for index in range(int(manifest["sample_count"])):
|
||||
writer.writerow([index, f"{index / sampling_frequency_hz:.17g}"])
|
||||
manifest.update(
|
||||
{
|
||||
"sampling_frequency_hz": sampling_frequency_hz,
|
||||
"sample_period_s": 1.0 / sampling_frequency_hz,
|
||||
"timestamp_file": timestamp_path.name,
|
||||
"register_values": register_values,
|
||||
}
|
||||
)
|
||||
(output_directory / f"{stem}.json").write_text(json.dumps(manifest, indent=2), encoding="utf-8")
|
||||
return manifest
|
||||
@@ -0,0 +1,58 @@
|
||||
"""Signal definitions consumed by generic acquisition and tracking engines."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import Literal
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class ComponentProfile:
|
||||
name: str
|
||||
channel_type: Literal["data", "pilot"]
|
||||
primary_code_length_chips: int
|
||||
primary_code_period_s: float
|
||||
secondary_code_length: int | None
|
||||
navigation_symbol_period_s: float | None
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class SignalProfile:
|
||||
name: str
|
||||
carrier_frequency_hz: float
|
||||
code_rate_hz: float
|
||||
modulation: str
|
||||
boc_subcarrier_rate_hz: float
|
||||
components: dict[str, ComponentProfile]
|
||||
|
||||
def component(self, name: str) -> ComponentProfile:
|
||||
try:
|
||||
return self.components[name]
|
||||
except KeyError as exc:
|
||||
raise ValueError(f"component {name!r} is not supported by profile {self.name}") from exc
|
||||
|
||||
|
||||
def galileo_e1_profile() -> SignalProfile:
|
||||
"""Galileo E1 parameters; code chips themselves come from a code provider."""
|
||||
return SignalProfile(
|
||||
name="GALILEO_E1",
|
||||
carrier_frequency_hz=1_575_420_000.0,
|
||||
code_rate_hz=1_023_000.0,
|
||||
modulation="BOC(1,1)",
|
||||
boc_subcarrier_rate_hz=1_023_000.0,
|
||||
components={
|
||||
"E1B": ComponentProfile("E1B", "data", 4092, 0.004, None, 0.004),
|
||||
"E1C": ComponentProfile("E1C", "pilot", 4092, 0.004, 25, None),
|
||||
},
|
||||
)
|
||||
|
||||
|
||||
PROFILES = {"GALILEO_E1": galileo_e1_profile}
|
||||
|
||||
|
||||
def get_signal_profile(name: str) -> SignalProfile:
|
||||
try:
|
||||
return PROFILES[name]()
|
||||
except KeyError as exc:
|
||||
raise ValueError(f"unknown signal profile {name!r}; available: {sorted(PROFILES)}") from exc
|
||||
|
||||
@@ -0,0 +1,7 @@
|
||||
"""Signal code, modulation, and generation."""
|
||||
|
||||
from .code_provider import CodeProvider, FileCodeProvider, SyntheticCodeProvider
|
||||
from .signal_generator import SatelliteSignalConfig, SyntheticSignalGenerator
|
||||
|
||||
__all__ = ["CodeProvider", "FileCodeProvider", "SyntheticCodeProvider", "SatelliteSignalConfig", "SyntheticSignalGenerator"]
|
||||
|
||||
@@ -0,0 +1,38 @@
|
||||
"""BOC modulation helpers with sample-continuous phase."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
def boc_subcarrier(
|
||||
sample_count: int,
|
||||
sampling_frequency_hz: float,
|
||||
subcarrier_frequency_hz: float,
|
||||
initial_phase_cycles: float = 0.0,
|
||||
) -> tuple[npt.NDArray[np.float64], float]:
|
||||
"""Generate sine-phased square-wave BOC subcarrier and ending phase in cycles."""
|
||||
if sample_count < 0 or sampling_frequency_hz <= 0 or subcarrier_frequency_hz <= 0:
|
||||
raise ValueError("invalid BOC generator parameters")
|
||||
phase = initial_phase_cycles + np.arange(sample_count, dtype=np.float64) * subcarrier_frequency_hz / sampling_frequency_hz
|
||||
wave = np.where(np.mod(phase, 1.0) < 0.5, 1.0, -1.0)
|
||||
ending_phase = float((initial_phase_cycles + sample_count * subcarrier_frequency_hz / sampling_frequency_hz) % 1.0)
|
||||
return wave, ending_phase
|
||||
|
||||
|
||||
def sample_primary_code(
|
||||
code: npt.NDArray[np.int8],
|
||||
sample_count: int,
|
||||
sampling_frequency_hz: float,
|
||||
code_rate_hz: float,
|
||||
initial_code_phase_chips: float = 0.0,
|
||||
) -> tuple[npt.NDArray[np.float64], float]:
|
||||
"""Nearest-chip code sampling supporting noninteger samples per chip."""
|
||||
if code.ndim != 1 or code.size == 0 or sampling_frequency_hz <= 0 or code_rate_hz <= 0:
|
||||
raise ValueError("invalid code sampling parameters")
|
||||
phases = initial_code_phase_chips + np.arange(sample_count, dtype=np.float64) * code_rate_hz / sampling_frequency_hz
|
||||
values = code[np.floor(phases).astype(np.int64) % code.size].astype(np.float64)
|
||||
ending_phase = float((initial_code_phase_chips + sample_count * code_rate_hz / sampling_frequency_hz) % code.size)
|
||||
return values, ending_phase
|
||||
|
||||
@@ -0,0 +1,82 @@
|
||||
"""Official-table integration point and deterministic synthetic test codes."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from abc import ABC, abstractmethod
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
|
||||
ChipArray = npt.NDArray[np.int8]
|
||||
|
||||
|
||||
class CodeProvider(ABC):
|
||||
@abstractmethod
|
||||
def primary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray:
|
||||
"""Return chips as +/-1 int8 values."""
|
||||
|
||||
@abstractmethod
|
||||
def secondary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray | None:
|
||||
"""Return one secondary-code period as +/-1, or None."""
|
||||
|
||||
|
||||
class SyntheticCodeProvider(CodeProvider):
|
||||
"""Deterministic random codes exclusively for simulation and verification.
|
||||
|
||||
These sequences are not Galileo ICD codes and cannot acquire live signals.
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def _chips(length: int, profile: str, component: str, prn: int, kind: str) -> ChipArray:
|
||||
seed = sum((index + 1) * ord(char) for index, char in enumerate(f"{profile}:{component}:{prn}:{kind}"))
|
||||
rng = np.random.default_rng(seed & 0xFFFFFFFF)
|
||||
return np.where(rng.integers(0, 2, length, dtype=np.int8) == 0, -1, 1).astype(np.int8)
|
||||
|
||||
def primary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray:
|
||||
spec = profile.component(component)
|
||||
return self._chips(spec.primary_code_length_chips, profile.name, component, prn, "primary")
|
||||
|
||||
def secondary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray | None:
|
||||
length = profile.component(component).secondary_code_length
|
||||
return None if length is None else self._chips(length, profile.name, component, prn, "secondary")
|
||||
|
||||
|
||||
class FileCodeProvider(CodeProvider):
|
||||
"""Load version-controlled code tables without embedding guessed sequences.
|
||||
|
||||
Files are ``<profile>_<component>_PRN<nnn>.txt`` and contain exactly one
|
||||
whitespace-separated row of either 0/1 or -1/+1 chips.
|
||||
"""
|
||||
|
||||
def __init__(self, primary_directory: str | Path, secondary_directory: str | Path) -> None:
|
||||
self.primary_directory = Path(primary_directory)
|
||||
self.secondary_directory = Path(secondary_directory)
|
||||
|
||||
@staticmethod
|
||||
def _load(path: Path, expected_length: int) -> ChipArray:
|
||||
if not path.exists():
|
||||
raise FileNotFoundError(f"official code table is not installed: {path}")
|
||||
raw = np.loadtxt(path, dtype=np.int8).reshape(-1)
|
||||
if raw.size != expected_length:
|
||||
raise ValueError(f"{path} has {raw.size} chips; expected {expected_length}")
|
||||
unique = set(raw.tolist())
|
||||
if unique <= {0, 1}:
|
||||
raw = 2 * raw - 1
|
||||
elif not unique <= {-1, 1}:
|
||||
raise ValueError(f"{path} must contain only 0/1 or -1/+1")
|
||||
return raw.astype(np.int8)
|
||||
|
||||
@staticmethod
|
||||
def _name(profile: SignalProfile, component: str, prn: int) -> str:
|
||||
return f"{profile.name.lower()}_{component.lower()}_PRN{prn:03d}.txt"
|
||||
|
||||
def primary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray:
|
||||
return self._load(self.primary_directory / self._name(profile, component, prn), profile.component(component).primary_code_length_chips)
|
||||
|
||||
def secondary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray | None:
|
||||
length = profile.component(component).secondary_code_length
|
||||
return None if length is None else self._load(self.secondary_directory / self._name(profile, component, prn), length)
|
||||
|
||||
@@ -0,0 +1,31 @@
|
||||
"""Profile-driven Galileo E1 local replica generation."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from .boc import boc_subcarrier, sample_primary_code
|
||||
from .code_provider import CodeProvider
|
||||
|
||||
|
||||
def e1_local_replica(
|
||||
profile: SignalProfile,
|
||||
provider: CodeProvider,
|
||||
component: str,
|
||||
prn: int,
|
||||
sample_count: int,
|
||||
sampling_frequency_hz: float,
|
||||
initial_code_phase_chips: float = 0.0,
|
||||
boc_enabled: bool = True,
|
||||
code_rate_hz: float | None = None,
|
||||
) -> npt.NDArray[np.float64]:
|
||||
primary = provider.primary_code(profile, component, prn)
|
||||
effective_code_rate_hz = profile.code_rate_hz if code_rate_hz is None else code_rate_hz
|
||||
code, _ = sample_primary_code(primary, sample_count, sampling_frequency_hz, effective_code_rate_hz, initial_code_phase_chips)
|
||||
if not boc_enabled:
|
||||
return code
|
||||
rate_scale = effective_code_rate_hz / profile.code_rate_hz
|
||||
subcarrier, _ = boc_subcarrier(sample_count, sampling_frequency_hz, profile.boc_subcarrier_rate_hz * rate_scale, initial_code_phase_chips % 1.0)
|
||||
return code * subcarrier
|
||||
@@ -0,0 +1,182 @@
|
||||
"""Streaming synthetic signal generator with exact state metadata."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Literal
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.types import GenerationTruth, SatelliteTruth
|
||||
from .code_provider import CodeProvider
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class SatelliteSignalConfig:
|
||||
prn: int
|
||||
component: Literal["E1B", "E1C"] = "E1C"
|
||||
code_delay_chips: float = 0.0
|
||||
doppler_hz: float = 0.0
|
||||
doppler_rate_hz_s: float = 0.0
|
||||
carrier_phase_rad: float = 0.0
|
||||
amplitude: float = 1.0
|
||||
cn0_db_hz: float | None = None
|
||||
navigation_bits: list[int] = field(default_factory=lambda: [1])
|
||||
navigation_bit_period_s: float = 0.004
|
||||
secondary_code_phase: int = 0
|
||||
|
||||
def validate(self, profile: SignalProfile) -> None:
|
||||
profile.component(self.component)
|
||||
if self.prn <= 0 or self.amplitude < 0:
|
||||
raise ValueError("PRN must be positive and amplitude nonnegative")
|
||||
if not self.navigation_bits or not set(self.navigation_bits) <= {-1, 1}:
|
||||
raise ValueError("navigation_bits must be a nonempty +/-1 list")
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class ChannelModelConfig:
|
||||
sampling_clock_offset_ppm: float = 0.0
|
||||
carrier_frequency_offset_hz: float = 0.0
|
||||
add_awgn: bool = True
|
||||
noise_std: float | None = None
|
||||
clip_amplitude: float | None = None
|
||||
quantization_bits: int | None = None
|
||||
seed: int = 1
|
||||
|
||||
|
||||
class SyntheticSignalGenerator:
|
||||
"""Generate real IF or analytic IQ while preserving time across blocks."""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
profile: SignalProfile,
|
||||
code_provider: CodeProvider,
|
||||
sampling_frequency_hz: float,
|
||||
intermediate_frequency_hz: float = 0.0,
|
||||
output_is_real_if: bool = False,
|
||||
satellites: list[SatelliteSignalConfig] | None = None,
|
||||
channel: ChannelModelConfig | None = None,
|
||||
) -> None:
|
||||
if sampling_frequency_hz <= 0:
|
||||
raise ValueError("sampling_frequency_hz must be positive")
|
||||
if output_is_real_if and not 0 < intermediate_frequency_hz < sampling_frequency_hz / 2:
|
||||
raise ValueError("real IF must lie strictly between DC and Nyquist")
|
||||
self.profile = profile
|
||||
self.code_provider = code_provider
|
||||
self.sampling_frequency_hz = sampling_frequency_hz
|
||||
self.intermediate_frequency_hz = intermediate_frequency_hz
|
||||
self.output_is_real_if = output_is_real_if
|
||||
self.satellites = satellites or []
|
||||
self.channel = channel or ChannelModelConfig()
|
||||
for satellite in self.satellites:
|
||||
satellite.validate(profile)
|
||||
self._sample_index = 0
|
||||
self._rng = np.random.default_rng(self.channel.seed)
|
||||
|
||||
@property
|
||||
def sample_index(self) -> int:
|
||||
return self._sample_index
|
||||
|
||||
def reset(self) -> None:
|
||||
self._sample_index = 0
|
||||
self._rng = np.random.default_rng(self.channel.seed)
|
||||
|
||||
@staticmethod
|
||||
def _sign_sequence(sequence: npt.NDArray[np.int8], indices: npt.NDArray[np.int64]) -> npt.NDArray[np.float64]:
|
||||
return sequence[np.mod(indices, sequence.size)].astype(np.float64)
|
||||
|
||||
def generate(self, sample_count: int) -> tuple[npt.NDArray[np.complex128] | npt.NDArray[np.float64], GenerationTruth]:
|
||||
if sample_count < 0:
|
||||
raise ValueError("sample_count must be nonnegative")
|
||||
start = self._sample_index
|
||||
nominal_t_s = (start + np.arange(sample_count, dtype=np.float64)) / self.sampling_frequency_hz
|
||||
signal_t_s = nominal_t_s * (1.0 + self.channel.sampling_clock_offset_ppm * 1e-6)
|
||||
output = np.zeros(sample_count, dtype=np.complex128)
|
||||
truths: list[SatelliteTruth] = []
|
||||
equivalent_noise_variance = 0.0
|
||||
|
||||
for satellite in self.satellites:
|
||||
primary = self.code_provider.primary_code(self.profile, satellite.component, satellite.prn)
|
||||
# Common line-of-sight dynamics scale carrier and spreading-code
|
||||
# frequencies by the same fractional Doppler. Keeping this coupling
|
||||
# makes carrier-aided code tracking physically self-consistent.
|
||||
code_phase_chips = (
|
||||
signal_t_s * self.profile.code_rate_hz
|
||||
+ signal_t_s * satellite.doppler_hz * self.profile.code_rate_hz / self.profile.carrier_frequency_hz
|
||||
+ 0.5 * signal_t_s**2 * satellite.doppler_rate_hz_s * self.profile.code_rate_hz / self.profile.carrier_frequency_hz
|
||||
- satellite.code_delay_chips
|
||||
)
|
||||
chip_indices = np.floor(code_phase_chips).astype(np.int64)
|
||||
spreading = self._sign_sequence(primary, chip_indices)
|
||||
subcarrier_phase = code_phase_chips
|
||||
subcarrier = np.where(np.mod(subcarrier_phase, 1.0) < 0.5, 1.0, -1.0)
|
||||
|
||||
component_spec = self.profile.component(satellite.component)
|
||||
if component_spec.channel_type == "data":
|
||||
bit_indices = np.floor(signal_t_s / satellite.navigation_bit_period_s).astype(np.int64)
|
||||
spreading *= self._sign_sequence(np.asarray(satellite.navigation_bits, dtype=np.int8), bit_indices)
|
||||
else:
|
||||
secondary = self.code_provider.secondary_code(self.profile, satellite.component, satellite.prn)
|
||||
if secondary is not None:
|
||||
epoch_indices = np.floor(signal_t_s / component_spec.primary_code_period_s).astype(np.int64)
|
||||
spreading *= self._sign_sequence(secondary, epoch_indices + satellite.secondary_code_phase)
|
||||
|
||||
phase_rad = (
|
||||
satellite.carrier_phase_rad
|
||||
+ 2.0 * np.pi * (self.intermediate_frequency_hz + satellite.doppler_hz + self.channel.carrier_frequency_offset_hz) * nominal_t_s
|
||||
+ np.pi * satellite.doppler_rate_hz_s * nominal_t_s**2
|
||||
)
|
||||
output += satellite.amplitude * spreading * subcarrier * np.exp(1j * phase_rad)
|
||||
if satellite.cn0_db_hz is not None:
|
||||
equivalent_noise_variance += satellite.amplitude**2 * self.sampling_frequency_hz / (10.0 ** (satellite.cn0_db_hz / 10.0))
|
||||
truths.append(
|
||||
SatelliteTruth(
|
||||
prn=satellite.prn,
|
||||
component=satellite.component,
|
||||
initial_code_delay_chips=satellite.code_delay_chips,
|
||||
initial_code_delay_samples=satellite.code_delay_chips * self.sampling_frequency_hz / self.profile.code_rate_hz,
|
||||
initial_doppler_hz=satellite.doppler_hz,
|
||||
doppler_rate_hz_s=satellite.doppler_rate_hz_s,
|
||||
initial_carrier_phase_rad=satellite.carrier_phase_rad,
|
||||
amplitude=satellite.amplitude,
|
||||
cn0_db_hz=satellite.cn0_db_hz,
|
||||
)
|
||||
)
|
||||
|
||||
if self.channel.add_awgn and sample_count:
|
||||
variance = self.channel.noise_std**2 if self.channel.noise_std is not None else equivalent_noise_variance
|
||||
if variance > 0:
|
||||
output += np.sqrt(variance / 2.0) * (
|
||||
self._rng.standard_normal(sample_count) + 1j * self._rng.standard_normal(sample_count)
|
||||
)
|
||||
result: npt.NDArray[np.complex128] | npt.NDArray[np.float64]
|
||||
if self.output_is_real_if:
|
||||
result = (np.sqrt(2.0) * output.real).astype(np.float64)
|
||||
else:
|
||||
result = output.astype(np.complex128)
|
||||
if self.channel.clip_amplitude is not None:
|
||||
limit = self.channel.clip_amplitude
|
||||
if np.iscomplexobj(result):
|
||||
result = (np.clip(result.real, -limit, limit) + 1j * np.clip(result.imag, -limit, limit)).astype(np.complex128)
|
||||
else:
|
||||
result = np.clip(result, -limit, limit).astype(np.float64)
|
||||
if self.channel.quantization_bits is not None:
|
||||
result = quantize_normalized(result, self.channel.quantization_bits)
|
||||
self._sample_index += sample_count
|
||||
truth = GenerationTruth(self.sampling_frequency_hz, self.intermediate_frequency_hz, self.output_is_real_if, start, sample_count, truths)
|
||||
return result, truth
|
||||
|
||||
|
||||
def quantize_normalized(samples: npt.NDArray[np.generic], bits: int) -> npt.NDArray[np.generic]:
|
||||
"""Quantize normalized signed samples to a symmetric fractional grid."""
|
||||
if not 2 <= bits <= 32:
|
||||
raise ValueError("quantization bits must be in [2, 32]")
|
||||
scale = float(1 << (bits - 1))
|
||||
maximum = 1.0 - 1.0 / scale
|
||||
if np.iscomplexobj(samples):
|
||||
real = np.round(np.clip(samples.real, -1.0, maximum) * scale) / scale
|
||||
imag = np.round(np.clip(samples.imag, -1.0, maximum) * scale) / scale
|
||||
return (real + 1j * imag).astype(samples.dtype)
|
||||
return (np.round(np.clip(samples, -1.0, maximum) * scale) / scale).astype(samples.dtype)
|
||||
@@ -0,0 +1,121 @@
|
||||
"""Code provider backed by an unmodified local SignalSim checkout.
|
||||
|
||||
The provider deliberately parses SignalSim's source at run time. It does not
|
||||
copy, guess, or silently replace Galileo code tables. This also makes the
|
||||
checkout commit part of the reproducibility boundary.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from functools import cached_property
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
import numpy as np
|
||||
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from .code_provider import ChipArray, CodeProvider
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class SignalSimConvention:
|
||||
"""Observed conventions for globsky/SignalSim commit 80b2621.
|
||||
|
||||
``code_sign`` describes table bit mapping (0 -> +1, 1 -> -1). SignalSim
|
||||
additionally multiplies E1B by -DataBit; that data-channel sign is kept out
|
||||
of the spreading-code provider and is reported separately.
|
||||
"""
|
||||
|
||||
prn_offset: int = 0
|
||||
code_phase_offset_chips: float = 0.0
|
||||
code_sign: int = 1
|
||||
boc_sign: int = 1
|
||||
doppler_sign: int = 1
|
||||
iq_order: str = "IQ"
|
||||
carrier_mixing_sign: int = -1
|
||||
|
||||
|
||||
class SignalSimCodeProvider(CodeProvider):
|
||||
"""Load Galileo E1B/E1C codes from a local SignalSim source tree."""
|
||||
|
||||
_E1_WORDS_PER_CODE = 128
|
||||
_E1_CODE_COUNT = 100
|
||||
_E1_CHIPS = 4092
|
||||
|
||||
def __init__(self, signalsim_path: str | Path) -> None:
|
||||
self.signalsim_path = Path(signalsim_path).expanduser().resolve()
|
||||
self.memory_code_path = self.signalsim_path / "src" / "MemoryCode.dat"
|
||||
self.pilot_bit_path = self.signalsim_path / "src" / "PilotBit.cpp"
|
||||
if not self.memory_code_path.is_file():
|
||||
raise FileNotFoundError(f"SignalSim E1 code table not found: {self.memory_code_path}")
|
||||
if not self.pilot_bit_path.is_file():
|
||||
raise FileNotFoundError(f"SignalSim pilot source not found: {self.pilot_bit_path}")
|
||||
|
||||
@cached_property
|
||||
def _e1_words(self) -> np.ndarray:
|
||||
text = self.memory_code_path.read_text(encoding="utf-8")
|
||||
match = re.search(
|
||||
r"E1MemoryCode\s*\[\s*100\s*\*\s*128\s*\]\s*=\s*\{(?P<body>.*?)\};",
|
||||
text,
|
||||
flags=re.DOTALL,
|
||||
)
|
||||
if match is None:
|
||||
raise ValueError(f"could not locate E1MemoryCode in {self.memory_code_path}")
|
||||
words = np.asarray(
|
||||
[int(token, 16) for token in re.findall(r"0x[0-9a-fA-F]+", match.group("body"))],
|
||||
dtype=np.uint32,
|
||||
)
|
||||
expected = self._E1_WORDS_PER_CODE * self._E1_CODE_COUNT
|
||||
if words.size != expected:
|
||||
raise ValueError(f"SignalSim E1 table has {words.size} words; expected {expected}")
|
||||
return words.reshape(self._E1_CODE_COUNT, self._E1_WORDS_PER_CODE)
|
||||
|
||||
@staticmethod
|
||||
def _unpack_msb_first(words: np.ndarray) -> ChipArray:
|
||||
# GetMemorySequence() divides a code into four 32-word sectors. It
|
||||
# consumes 1023 MSB-first bits from each 1024-bit sector and deliberately
|
||||
# skips the last padding bit before advancing to the next sector.
|
||||
sectors = words.reshape(4, 32)
|
||||
sector_indices = np.arange(1023, dtype=np.int64)
|
||||
bits = np.concatenate(
|
||||
[(sector[sector_indices >> 5] >> (31 - (sector_indices & 31))) & 1 for sector in sectors]
|
||||
)
|
||||
# SignalSim's modulator maps a stored 0 to +1 and stored 1 to -1.
|
||||
return np.where(bits == 0, 1, -1).astype(np.int8)
|
||||
|
||||
def primary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray:
|
||||
if profile.name != "GALILEO_E1":
|
||||
raise ValueError(f"SignalSim E1 provider cannot serve profile {profile.name!r}")
|
||||
if not 1 <= prn <= 50:
|
||||
raise ValueError("SignalSim Galileo E1 PRN must be in 1..50")
|
||||
normalized = component.upper()
|
||||
if normalized == "E1B":
|
||||
table_index = prn - 1
|
||||
elif normalized == "E1C":
|
||||
table_index = prn + 49
|
||||
else:
|
||||
raise ValueError(f"SignalSim E1 provider does not support component {component!r}")
|
||||
expected = profile.component(normalized).primary_code_length_chips
|
||||
if expected != self._E1_CHIPS:
|
||||
raise ValueError(f"profile requests {expected} E1 chips; SignalSim table contains {self._E1_CHIPS}")
|
||||
return self._unpack_msb_first(self._e1_words[table_index])
|
||||
|
||||
def secondary_code(self, profile: SignalProfile, component: str, prn: int) -> ChipArray | None:
|
||||
del prn # SignalSim uses one common 25-symbol E1C sequence.
|
||||
normalized = component.upper()
|
||||
if normalized == "E1B":
|
||||
return None
|
||||
if normalized != "E1C":
|
||||
raise ValueError(f"SignalSim E1 provider does not support component {component!r}")
|
||||
length = profile.component(normalized).secondary_code_length
|
||||
text = self.pilot_bit_path.read_text(encoding="utf-8")
|
||||
match = re.search(r"SecondaryCodeE1\s*\[\s*1\s*\]\s*=\s*\{\s*(0x[0-9a-fA-F]+)\s*\}", text)
|
||||
if match is None:
|
||||
raise ValueError(f"could not locate SecondaryCodeE1 in {self.pilot_bit_path}")
|
||||
if length != 25:
|
||||
raise ValueError(f"profile requests secondary length {length}; SignalSim implements 25")
|
||||
word = int(match.group(1), 16)
|
||||
# PilotBit.cpp addresses the word LSB first by secondary-code position.
|
||||
bits = np.asarray([(word >> index) & 1 for index in range(length)], dtype=np.int8)
|
||||
return np.where(bits == 0, 1, -1).astype(np.int8)
|
||||
@@ -0,0 +1,7 @@
|
||||
"""Stateful single- and multi-channel tracking."""
|
||||
|
||||
from .tracking_channel import TrackingChannel
|
||||
from .tracking_manager import TrackingManager
|
||||
|
||||
__all__ = ["TrackingChannel", "TrackingManager"]
|
||||
|
||||
@@ -0,0 +1,31 @@
|
||||
"""Moment-based prompt C/N0 estimator."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from collections import deque
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class Cn0Estimator:
|
||||
def __init__(self, integration_time_s: float, window_epochs: int = 50, *, data_sign_insensitive: bool = False) -> None:
|
||||
self.integration_time_s = integration_time_s
|
||||
self.data_sign_insensitive = data_sign_insensitive
|
||||
self.values: deque[complex] = deque(maxlen=window_epochs)
|
||||
|
||||
def update(self, prompt: complex) -> float | None:
|
||||
self.values.append(prompt)
|
||||
if len(self.values) < max(5, self.values.maxlen // 4):
|
||||
return None
|
||||
values = np.asarray(self.values)
|
||||
if self.data_sign_insensitive:
|
||||
# A locked BPSK data channel lies on the I axis. Decision-directed
|
||||
# polarity wipeoff retains the carrier phasor while removing data.
|
||||
decisions = np.where(values.real >= 0.0, 1.0, -1.0)
|
||||
values = values * decisions
|
||||
mean_power = float(np.mean(np.abs(values) ** 2))
|
||||
coherent_power = abs(np.mean(values)) ** 2
|
||||
noise_power = mean_power - coherent_power
|
||||
if coherent_power <= 0 or noise_power <= 0:
|
||||
return None
|
||||
return float(10 * np.log10(coherent_power / noise_power / self.integration_time_s))
|
||||
@@ -0,0 +1,75 @@
|
||||
"""Vector and sample-oriented Galileo correlator banks."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.signals.code_provider import CodeProvider
|
||||
from galileo_rx.types import CorrelatorOutputs
|
||||
from .nco import CarrierNco, CodeNco
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CorrelatorBank:
|
||||
profile: SignalProfile
|
||||
provider: CodeProvider
|
||||
component: str
|
||||
prn: int
|
||||
carrier_nco: CarrierNco
|
||||
code_nco: CodeNco
|
||||
early_late_spacing_chips: float
|
||||
very_early_late_spacing_chips: float
|
||||
_code: npt.NDArray[np.int8] = field(init=False, repr=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
self._code = self.provider.primary_code(self.profile, self.component, self.prn)
|
||||
|
||||
def _replica(self, phases_chips: np.ndarray, offset_chips: float) -> np.ndarray:
|
||||
tap_phase = phases_chips + offset_chips
|
||||
code = self._code[np.floor(tap_phase).astype(np.int64) % self._code.size]
|
||||
boc = np.where(np.mod(tap_phase, 1.0) < 0.5, 1.0, -1.0)
|
||||
return code * boc
|
||||
|
||||
def correlate(self, samples: npt.NDArray[np.generic], secondary_sign: int = 1) -> CorrelatorOutputs:
|
||||
array = np.asarray(samples, dtype=np.complex128)
|
||||
carrier_phase = self.carrier_nco.phases(array.size)
|
||||
code_phase = self.code_nco.phases(array.size)
|
||||
wiped = array * np.exp(-1j * carrier_phase) * secondary_sign
|
||||
taps = {
|
||||
"very_early": -self.very_early_late_spacing_chips,
|
||||
"early": -self.early_late_spacing_chips,
|
||||
"prompt": 0.0,
|
||||
"late": self.early_late_spacing_chips,
|
||||
"very_late": self.very_early_late_spacing_chips,
|
||||
}
|
||||
values = {name: complex(np.sum(wiped * self._replica(code_phase, offset))) for name, offset in taps.items()}
|
||||
return CorrelatorOutputs(**values)
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CycleCorrelatorState:
|
||||
"""One-sample transaction model used as an RTL comparison oracle."""
|
||||
|
||||
accumulators: dict[str, complex]
|
||||
sample_count: int = 0
|
||||
|
||||
|
||||
class SampleStreamCorrelator:
|
||||
def __init__(self, bank: CorrelatorBank) -> None:
|
||||
self.bank = bank
|
||||
self.state = CycleCorrelatorState({name: 0j for name in ("very_early", "early", "prompt", "late", "very_late")})
|
||||
|
||||
def accept(self, sample: complex, secondary_sign: int = 1) -> None:
|
||||
outputs = self.bank.correlate(np.asarray([sample]), secondary_sign)
|
||||
for name in self.state.accumulators:
|
||||
self.state.accumulators[name] += getattr(outputs, name)
|
||||
self.state.sample_count += 1
|
||||
|
||||
def dump(self) -> CorrelatorOutputs:
|
||||
result = CorrelatorOutputs(**self.state.accumulators)
|
||||
self.state = CycleCorrelatorState({name: 0j for name in self.state.accumulators})
|
||||
return result
|
||||
@@ -0,0 +1,51 @@
|
||||
"""Tracking discriminators with explicit output units."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def pll_two_quadrant_rad(prompt: complex) -> float:
|
||||
"""Decision-directed carrier phase error in radians, range +/-pi/2."""
|
||||
return float(np.arctan(prompt.imag / prompt.real)) if prompt.real != 0 else float(np.sign(prompt.imag) * np.pi / 2)
|
||||
|
||||
|
||||
def pll_four_quadrant_rad(prompt: complex) -> float:
|
||||
"""Pilot carrier phase error in radians, range +/-pi."""
|
||||
return float(np.arctan2(prompt.imag, prompt.real))
|
||||
|
||||
|
||||
def fll_phase_difference_hz(previous_prompt: complex | None, prompt: complex, integration_time_s: float) -> float:
|
||||
"""Cross-product/phase-difference FLL error in hertz."""
|
||||
if previous_prompt is None or integration_time_s <= 0 or previous_prompt == 0 or prompt == 0:
|
||||
return 0.0
|
||||
return float(np.angle(prompt * np.conjugate(previous_prompt)) / (2.0 * np.pi * integration_time_s))
|
||||
|
||||
|
||||
def fll_bpsk_phase_difference_hz(previous_prompt: complex | None, prompt: complex, integration_time_s: float) -> float:
|
||||
"""BPSK-data-insensitive phase-difference FLL error in hertz.
|
||||
|
||||
Squaring the inter-epoch phasor removes a possible pi-radian data-symbol
|
||||
transition. The unambiguous range is consequently +/-1/(4*T).
|
||||
"""
|
||||
if previous_prompt is None or integration_time_s <= 0 or previous_prompt == 0 or prompt == 0:
|
||||
return 0.0
|
||||
difference = prompt * np.conjugate(previous_prompt)
|
||||
return float(np.angle(difference * difference) / (4.0 * np.pi * integration_time_s))
|
||||
|
||||
|
||||
def dll_normalized_envelope_chips(early: complex, late: complex, spacing_chips: float) -> float:
|
||||
"""Normalized E-minus-L envelope, calibrated to chip units near lock."""
|
||||
early_mag, late_mag = abs(early), abs(late)
|
||||
return float(spacing_chips * (early_mag - late_mag) / max(early_mag + late_mag, np.finfo(float).tiny))
|
||||
|
||||
|
||||
def dll_noncoherent_power_chips(early: complex, late: complex, spacing_chips: float) -> float:
|
||||
early_power, late_power = abs(early) ** 2, abs(late) ** 2
|
||||
return float(spacing_chips * (early_power - late_power) / max(early_power + late_power, np.finfo(float).tiny))
|
||||
|
||||
|
||||
def dll_dot_product_chips(early: complex, late: complex, prompt: complex, spacing_chips: float) -> float:
|
||||
numerator = np.real((early - late) * np.conjugate(prompt))
|
||||
denominator = max(abs(prompt) ** 2, np.finfo(float).tiny)
|
||||
return float(0.5 * spacing_chips * numerator / denominator)
|
||||
@@ -0,0 +1,39 @@
|
||||
"""Lock metrics and hysteretic failure counting."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
|
||||
from galileo_rx.types import CorrelatorOutputs
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class LockStatus:
|
||||
carrier_metric: float
|
||||
code_metric: float
|
||||
locked: bool
|
||||
cycle_slip: bool
|
||||
side_peak_suspected: bool
|
||||
|
||||
|
||||
class LockDetector:
|
||||
def __init__(self, failure_limit: int = 10) -> None:
|
||||
self.failure_limit = failure_limit
|
||||
self.failures = 0
|
||||
self.previous_phase: float | None = None
|
||||
|
||||
def update(self, correlators: CorrelatorOutputs) -> LockStatus:
|
||||
prompt_power = abs(correlators.prompt) ** 2
|
||||
carrier_metric = float(abs(correlators.prompt.real) / max(abs(correlators.prompt), np.finfo(float).tiny))
|
||||
code_metric = float(prompt_power / max(0.5 * (abs(correlators.early) ** 2 + abs(correlators.late) ** 2), np.finfo(float).tiny))
|
||||
phase = float(np.angle(correlators.prompt))
|
||||
cycle_slip = self.previous_phase is not None and abs(float(np.angle(np.exp(1j * (phase - self.previous_phase))))) > 0.9 * np.pi
|
||||
self.previous_phase = phase
|
||||
locked = carrier_metric > 0.55 and code_metric > 0.5
|
||||
self.failures = 0 if locked else self.failures + 1
|
||||
vevl = max(abs(correlators.very_early), abs(correlators.very_late))
|
||||
side_peak = vevl > 1.25 * max(abs(correlators.prompt), np.finfo(float).tiny)
|
||||
return LockStatus(carrier_metric, code_metric, locked, cycle_slip, side_peak)
|
||||
|
||||
@@ -0,0 +1,70 @@
|
||||
"""Loop filters derived from noise bandwidth and damping parameters."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
|
||||
def second_order_natural_frequency_rad_s(noise_bandwidth_hz: float, damping_ratio: float) -> float:
|
||||
"""Convert one-sided equivalent noise bandwidth to natural radian frequency."""
|
||||
if noise_bandwidth_hz <= 0 or damping_ratio <= 0:
|
||||
raise ValueError("bandwidth and damping ratio must be positive")
|
||||
return 8.0 * damping_ratio * noise_bandwidth_hz / (4.0 * damping_ratio**2 + 1.0)
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class SecondOrderLoopFilter:
|
||||
noise_bandwidth_hz: float
|
||||
damping_ratio: float = 0.707
|
||||
integrator: float = 0.0
|
||||
proportional_gain: float = field(init=False)
|
||||
integral_gain: float = field(init=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
wn = second_order_natural_frequency_rad_s(self.noise_bandwidth_hz, self.damping_ratio)
|
||||
self.proportional_gain = 2.0 * self.damping_ratio * wn
|
||||
self.integral_gain = wn**2
|
||||
|
||||
def update(self, error: float, integration_time_s: float) -> float:
|
||||
if integration_time_s <= 0:
|
||||
raise ValueError("integration_time_s must be positive")
|
||||
self.integrator += self.integral_gain * error * integration_time_s
|
||||
return self.proportional_gain * error + self.integrator
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class ThirdOrderLoopFilter:
|
||||
"""Type-III filter using a repeated stable pole derived from Bn.
|
||||
|
||||
The continuous correction is ``3*w*e + 3*w^2*integral(e) +
|
||||
w^3*double_integral(e)``. This explicit pole placement provides a
|
||||
deterministic starting point for later loop optimization.
|
||||
"""
|
||||
|
||||
noise_bandwidth_hz: float
|
||||
damping_ratio: float = 0.707
|
||||
first_integrator: float = 0.0
|
||||
second_integrator: float = 0.0
|
||||
omega_rad_s: float = field(init=False)
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
self.omega_rad_s = second_order_natural_frequency_rad_s(self.noise_bandwidth_hz, self.damping_ratio)
|
||||
|
||||
def update(self, error: float, integration_time_s: float) -> float:
|
||||
if integration_time_s <= 0:
|
||||
raise ValueError("integration_time_s must be positive")
|
||||
self.first_integrator += error * integration_time_s
|
||||
self.second_integrator += self.first_integrator * integration_time_s
|
||||
w = self.omega_rad_s
|
||||
return 3.0 * w * error + 3.0 * w**2 * self.first_integrator + w**3 * self.second_integrator
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class FirstOrderFllFilter:
|
||||
noise_bandwidth_hz: float
|
||||
output_hz: float = 0.0
|
||||
|
||||
def update(self, error_hz: float, integration_time_s: float) -> float:
|
||||
alpha = min(1.0, 2.0 * self.noise_bandwidth_hz * integration_time_s)
|
||||
self.output_hz += alpha * (error_hz - self.output_hz)
|
||||
return self.output_hz
|
||||
@@ -0,0 +1,33 @@
|
||||
"""Carrier and fractional-chip code NCO states."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CarrierNco:
|
||||
sampling_frequency_hz: float
|
||||
frequency_hz: float
|
||||
phase_rad: float = 0.0
|
||||
|
||||
def phases(self, sample_count: int) -> np.ndarray:
|
||||
phase = self.phase_rad + 2.0 * np.pi * self.frequency_hz / self.sampling_frequency_hz * np.arange(sample_count)
|
||||
self.phase_rad = float((self.phase_rad + 2.0 * np.pi * self.frequency_hz / self.sampling_frequency_hz * sample_count) % (2.0 * np.pi))
|
||||
return phase
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CodeNco:
|
||||
sampling_frequency_hz: float
|
||||
code_length_chips: int
|
||||
frequency_hz: float
|
||||
phase_chips: float = 0.0
|
||||
|
||||
def phases(self, sample_count: int) -> np.ndarray:
|
||||
phase = self.phase_chips + self.frequency_hz / self.sampling_frequency_hz * np.arange(sample_count)
|
||||
self.phase_chips = float((self.phase_chips + self.frequency_hz / self.sampling_frequency_hz * sample_count) % self.code_length_chips)
|
||||
return phase
|
||||
|
||||
@@ -0,0 +1,133 @@
|
||||
"""Navigation-symbol recovery and pilot secondary-code synchronization."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class SymbolComparison:
|
||||
recovered_symbols: int
|
||||
compared_symbols: int
|
||||
alignment_offset_symbols: int
|
||||
sign_inverted: bool
|
||||
symbol_errors: int
|
||||
symbol_error_rate: float
|
||||
|
||||
|
||||
@dataclass(frozen=True, slots=True)
|
||||
class SecondaryCodeSyncResult:
|
||||
phase: int
|
||||
residual_frequency_hz: float
|
||||
normalized_correlation: float
|
||||
compared_periods: int
|
||||
sign_errors: int
|
||||
|
||||
|
||||
def recover_bpsk_symbols(
|
||||
prompt_values: npt.ArrayLike,
|
||||
epochs_per_symbol: int,
|
||||
*,
|
||||
ignored_epochs: int = 0,
|
||||
) -> tuple[npt.NDArray[np.int8], npt.NDArray[np.float64]]:
|
||||
"""Accumulate decision-directed prompt I and return +/-1 raw symbols."""
|
||||
prompts = np.asarray(prompt_values, dtype=np.complex128).reshape(-1)
|
||||
if epochs_per_symbol <= 0 or ignored_epochs < 0:
|
||||
raise ValueError("epochs_per_symbol must be positive and ignored_epochs nonnegative")
|
||||
usable = prompts[ignored_epochs:]
|
||||
symbol_count = usable.size // epochs_per_symbol
|
||||
coherent_i = usable[: symbol_count * epochs_per_symbol].real.reshape(symbol_count, epochs_per_symbol).sum(axis=1)
|
||||
symbols = np.where(coherent_i >= 0.0, 1, -1).astype(np.int8)
|
||||
return symbols, np.asarray(coherent_i, dtype=np.float64)
|
||||
|
||||
|
||||
def compare_symbol_streams(
|
||||
recovered: npt.ArrayLike,
|
||||
expected: npt.ArrayLike,
|
||||
*,
|
||||
max_alignment_offset_symbols: int,
|
||||
minimum_compared_symbols: int = 1,
|
||||
) -> SymbolComparison:
|
||||
"""Compare +/-1 streams with bounded offset and one global sign ambiguity."""
|
||||
rx = np.asarray(recovered, dtype=np.int8).reshape(-1)
|
||||
truth = np.asarray(expected, dtype=np.int8).reshape(-1)
|
||||
if max_alignment_offset_symbols < 0 or minimum_compared_symbols <= 0:
|
||||
raise ValueError("alignment bound must be nonnegative and comparison length positive")
|
||||
if not set(rx.tolist()) <= {-1, 1} or not set(truth.tolist()) <= {-1, 1}:
|
||||
raise ValueError("symbol streams must contain only +/-1")
|
||||
candidates: list[tuple[int, int, int, bool, int]] = []
|
||||
for offset in range(-max_alignment_offset_symbols, max_alignment_offset_symbols + 1):
|
||||
rx_start = max(0, -offset)
|
||||
truth_start = max(0, offset)
|
||||
count = min(rx.size - rx_start, truth.size - truth_start)
|
||||
if count < minimum_compared_symbols:
|
||||
continue
|
||||
lhs = rx[rx_start : rx_start + count]
|
||||
rhs = truth[truth_start : truth_start + count]
|
||||
for inverted in (False, True):
|
||||
errors = int(np.count_nonzero(lhs != (-rhs if inverted else rhs)))
|
||||
candidates.append((errors, -count, offset, inverted, count))
|
||||
if not candidates:
|
||||
raise ValueError("streams do not overlap for the requested bounded alignment")
|
||||
errors, negative_count, offset, inverted, count = min(candidates)
|
||||
del negative_count
|
||||
return SymbolComparison(rx.size, count, offset, inverted, errors, errors / count)
|
||||
|
||||
|
||||
def synchronize_secondary_code(
|
||||
prompt_values: npt.ArrayLike,
|
||||
secondary_code: npt.ArrayLike,
|
||||
integration_time_s: float,
|
||||
*,
|
||||
residual_frequency_min_hz: float = -125.0,
|
||||
residual_frequency_max_hz: float = 125.0,
|
||||
residual_frequency_step_hz: float = 1.0,
|
||||
) -> SecondaryCodeSyncResult:
|
||||
"""Jointly search secondary phase and residual carrier frequency.
|
||||
|
||||
The correlation is complex, so an unknown constant carrier phase does not
|
||||
affect the metric. Frequency search prevents residual phase rotation from
|
||||
destroying the 25-period coherent sum.
|
||||
"""
|
||||
prompts = np.asarray(prompt_values, dtype=np.complex128).reshape(-1)
|
||||
secondary = np.asarray(secondary_code, dtype=np.int8).reshape(-1)
|
||||
if prompts.size == 0 or secondary.size == 0 or integration_time_s <= 0:
|
||||
raise ValueError("nonempty prompts/secondary and positive integration time are required")
|
||||
if not set(secondary.tolist()) <= {-1, 1}:
|
||||
raise ValueError("secondary code must contain only +/-1")
|
||||
frequencies = np.arange(
|
||||
residual_frequency_min_hz,
|
||||
residual_frequency_max_hz + 0.5 * residual_frequency_step_hz,
|
||||
residual_frequency_step_hz,
|
||||
dtype=np.float64,
|
||||
)
|
||||
if residual_frequency_step_hz <= 0 or frequencies.size == 0:
|
||||
raise ValueError("invalid residual-frequency search")
|
||||
epochs = np.arange(prompts.size, dtype=np.float64)
|
||||
carrier = np.exp(-2j * np.pi * frequencies[:, None] * integration_time_s * epochs[None, :])
|
||||
derotated = carrier * prompts[None, :]
|
||||
denominator = max(float(np.sum(np.abs(prompts))), np.finfo(np.float64).tiny)
|
||||
best_metric = -1.0
|
||||
best_phase = 0
|
||||
best_frequency_index = 0
|
||||
best_sum = 0j
|
||||
for phase in range(secondary.size):
|
||||
signs = secondary[(phase + np.arange(prompts.size)) % secondary.size]
|
||||
sums = derotated @ signs.astype(np.float64)
|
||||
index = int(np.argmax(np.abs(sums)))
|
||||
metric = float(abs(sums[index]) / denominator)
|
||||
if metric > best_metric:
|
||||
best_metric = metric
|
||||
best_phase = phase
|
||||
best_frequency_index = index
|
||||
best_sum = complex(sums[index])
|
||||
best_frequency = float(frequencies[best_frequency_index])
|
||||
signs = secondary[(best_phase + np.arange(prompts.size)) % secondary.size]
|
||||
corrected = prompts * np.exp(-2j * np.pi * best_frequency * integration_time_s * epochs)
|
||||
corrected *= np.exp(-1j * np.angle(best_sum))
|
||||
decisions = np.where(corrected.real >= 0.0, 1, -1)
|
||||
errors = int(np.count_nonzero(decisions != signs))
|
||||
return SecondaryCodeSyncResult(best_phase, best_frequency, best_metric, prompts.size, errors)
|
||||
@@ -0,0 +1,174 @@
|
||||
"""Independent Galileo component tracking channel state machine."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.config import TrackingConfig
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.signals.code_provider import CodeProvider
|
||||
from galileo_rx.types import AcquisitionResult, CorrelatorOutputs, TrackingMeasurement, TrackingState
|
||||
from .cn0_estimator import Cn0Estimator
|
||||
from .correlator import CorrelatorBank
|
||||
from .discriminators import (
|
||||
dll_normalized_envelope_chips,
|
||||
fll_bpsk_phase_difference_hz,
|
||||
fll_phase_difference_hz,
|
||||
pll_four_quadrant_rad,
|
||||
pll_two_quadrant_rad,
|
||||
)
|
||||
from .lock_detector import LockDetector
|
||||
from .loop_filters import FirstOrderFllFilter, SecondOrderLoopFilter, ThirdOrderLoopFilter
|
||||
from .nco import CarrierNco, CodeNco
|
||||
|
||||
|
||||
class TrackingChannel:
|
||||
def __init__(
|
||||
self,
|
||||
profile: SignalProfile,
|
||||
provider: CodeProvider,
|
||||
sampling_frequency_hz: float,
|
||||
config: TrackingConfig,
|
||||
acquisition: AcquisitionResult,
|
||||
start_sample_index: int = 0,
|
||||
secondary_code_phase_at_start: int | None = None,
|
||||
) -> None:
|
||||
config.validate()
|
||||
self.profile = profile
|
||||
self.provider = provider
|
||||
self.sampling_frequency_hz = sampling_frequency_hz
|
||||
self.config = config
|
||||
self.acquisition = acquisition
|
||||
self.component_spec = profile.component(acquisition.component)
|
||||
self.integration_samples = int(round(config.integration_ms * 1e-3 * sampling_frequency_hz))
|
||||
if not np.isclose(self.integration_samples, config.integration_ms * 1e-3 * sampling_frequency_hz, atol=1e-6):
|
||||
raise ValueError("tracking integration is not an integer number of samples")
|
||||
initial_code_phase = (
|
||||
-acquisition.code_phase_chips + start_sample_index * profile.code_rate_hz / sampling_frequency_hz
|
||||
) % self.component_spec.primary_code_length_chips
|
||||
initial_carrier_phase = 2.0 * np.pi * acquisition.carrier_frequency_hz * start_sample_index / sampling_frequency_hz
|
||||
self.carrier_nco = CarrierNco(sampling_frequency_hz, acquisition.carrier_frequency_hz, initial_carrier_phase % (2 * np.pi))
|
||||
initial_code_frequency_hz = profile.code_rate_hz * (
|
||||
1.0 + acquisition.doppler_hz / profile.carrier_frequency_hz
|
||||
)
|
||||
self.code_nco = CodeNco(
|
||||
sampling_frequency_hz,
|
||||
self.component_spec.primary_code_length_chips,
|
||||
initial_code_frequency_hz,
|
||||
initial_code_phase,
|
||||
)
|
||||
self.correlator = CorrelatorBank(
|
||||
profile, provider, acquisition.component, acquisition.prn, self.carrier_nco, self.code_nco,
|
||||
config.early_late_spacing_chips, config.very_early_late_spacing_chips,
|
||||
)
|
||||
self.pll_filter = ThirdOrderLoopFilter(config.pll_noise_bandwidth_hz, config.damping_ratio) if config.pll_order == 3 else SecondOrderLoopFilter(config.pll_noise_bandwidth_hz, config.damping_ratio)
|
||||
self.dll_filter = SecondOrderLoopFilter(config.dll_noise_bandwidth_hz, config.damping_ratio)
|
||||
self.fll_filter = FirstOrderFllFilter(config.fll_noise_bandwidth_hz)
|
||||
self.lock_detector = LockDetector(config.loss_of_lock_epochs)
|
||||
self.cn0_estimator = Cn0Estimator(
|
||||
config.integration_ms * 1e-3,
|
||||
data_sign_insensitive=self.component_spec.channel_type == "data",
|
||||
)
|
||||
self.state = TrackingState.ACQUIRED
|
||||
self.epoch = 0
|
||||
self.previous_prompt: complex | None = None
|
||||
self.measurements: list[TrackingMeasurement] = []
|
||||
self.last_correlators = CorrelatorOutputs()
|
||||
self.nominal_carrier_hz = acquisition.carrier_frequency_hz
|
||||
# Acquisition reports the NCO frequency as IF center plus Doppler. Code
|
||||
# aiding must use range-rate Doppler only; feeding the IF itself into
|
||||
# f_code/f_carrier would be a MHz-scale unit-boundary error.
|
||||
self.carrier_center_hz = acquisition.carrier_frequency_hz - acquisition.doppler_hz
|
||||
self._pll_reference_frequency_hz = acquisition.carrier_frequency_hz
|
||||
self._pll_active = False
|
||||
self.secondary_epoch_offset = (
|
||||
int(secondary_code_phase_at_start)
|
||||
if secondary_code_phase_at_start is not None
|
||||
else int(np.floor(start_sample_index / (self.component_spec.primary_code_period_s * sampling_frequency_hz)))
|
||||
)
|
||||
|
||||
def _secondary_sign(self) -> int:
|
||||
if self.component_spec.channel_type != "pilot" or self.config.secondary_code_mode == "none":
|
||||
return 1
|
||||
secondary = self.provider.secondary_code(self.profile, self.acquisition.component, self.acquisition.prn)
|
||||
if secondary is None:
|
||||
return 1
|
||||
return int(secondary[(self.secondary_epoch_offset + self.epoch) % secondary.size])
|
||||
|
||||
def process_epoch(self, samples: npt.NDArray[np.generic]) -> TrackingMeasurement:
|
||||
if len(samples) != self.integration_samples:
|
||||
raise ValueError(f"tracking epoch requires {self.integration_samples} samples")
|
||||
dt_s = self.config.integration_ms * 1e-3
|
||||
correlators = self.correlator.correlate(samples, self._secondary_sign())
|
||||
self.last_correlators = correlators
|
||||
pll_error = pll_two_quadrant_rad(correlators.prompt) if self.component_spec.channel_type == "data" else pll_four_quadrant_rad(correlators.prompt)
|
||||
fll_error = (
|
||||
fll_bpsk_phase_difference_hz(self.previous_prompt, correlators.prompt, dt_s)
|
||||
if self.component_spec.channel_type == "data"
|
||||
else fll_phase_difference_hz(self.previous_prompt, correlators.prompt, dt_s)
|
||||
)
|
||||
dll_error = dll_normalized_envelope_chips(correlators.early, correlators.late, self.config.early_late_spacing_chips)
|
||||
elapsed_ms = (self.epoch + 1) * self.config.integration_ms
|
||||
if self.state == TrackingState.ACQUIRED:
|
||||
self.state = TrackingState.FLL_PULL_IN if self.config.use_fll_assist else TrackingState.PLL_TRACK
|
||||
use_fll = self.config.use_fll_assist and elapsed_ms <= self.config.pull_in_duration_ms
|
||||
if use_fll:
|
||||
correction_hz = self.fll_filter.update(fll_error, dt_s)
|
||||
self.carrier_nco.frequency_hz += correction_hz
|
||||
self.state = TrackingState.FLL_PULL_IN
|
||||
else:
|
||||
if not self._pll_active:
|
||||
self._pll_reference_frequency_hz = self.carrier_nco.frequency_hz
|
||||
self._pll_active = True
|
||||
correction_hz = self.pll_filter.update(pll_error, dt_s) / (2.0 * np.pi)
|
||||
self.carrier_nco.frequency_hz = self._pll_reference_frequency_hz + correction_hz
|
||||
self.state = TrackingState.PLL_TRACK
|
||||
code_correction_hz = self.dll_filter.update(dll_error, dt_s)
|
||||
carrier_doppler_hz = self.carrier_nco.frequency_hz - self.carrier_center_hz
|
||||
aiding_hz = (
|
||||
carrier_doppler_hz / self.profile.carrier_frequency_hz * self.profile.code_rate_hz
|
||||
if self.config.carrier_aiding
|
||||
else 0.0
|
||||
)
|
||||
# Positive E-L means the local phase must be retarded for the phase
|
||||
# convention used by CodeNco, hence the negative actuator sign.
|
||||
self.code_nco.frequency_hz = self.profile.code_rate_hz + aiding_hz - code_correction_hz
|
||||
cn0 = self.cn0_estimator.update(correlators.prompt)
|
||||
lock = self.lock_detector.update(correlators)
|
||||
reacquire = self.lock_detector.failures >= self.config.loss_of_lock_epochs
|
||||
if reacquire:
|
||||
self.state = TrackingState.REACQUISITION
|
||||
elif cn0 is not None and cn0 < self.config.weak_signal_cn0_db_hz:
|
||||
self.state = TrackingState.WEAK_SIGNAL
|
||||
measurement = TrackingMeasurement(
|
||||
epoch=self.epoch,
|
||||
time_s=(self.epoch + 1) * dt_s,
|
||||
state=self.state,
|
||||
prompt_i=float(correlators.prompt.real),
|
||||
prompt_q=float(correlators.prompt.imag),
|
||||
early_magnitude=float(abs(correlators.early)),
|
||||
late_magnitude=float(abs(correlators.late)),
|
||||
pll_error_rad=pll_error,
|
||||
fll_error_hz=fll_error,
|
||||
dll_error_chips=dll_error,
|
||||
carrier_frequency_hz=self.carrier_nco.frequency_hz,
|
||||
code_frequency_hz=self.code_nco.frequency_hz,
|
||||
carrier_phase_rad=self.carrier_nco.phase_rad,
|
||||
code_phase_chips=self.code_nco.phase_chips,
|
||||
cn0_db_hz=cn0,
|
||||
carrier_lock_metric=lock.carrier_metric,
|
||||
code_lock_metric=lock.code_metric,
|
||||
cycle_slip=lock.cycle_slip,
|
||||
side_peak_suspected=lock.side_peak_suspected,
|
||||
reacquisition_requested=reacquire,
|
||||
)
|
||||
self.previous_prompt = correlators.prompt
|
||||
self.measurements.append(measurement)
|
||||
self.epoch += 1
|
||||
return measurement
|
||||
|
||||
def process(self, samples: npt.NDArray[np.generic]) -> list[TrackingMeasurement]:
|
||||
array = np.asarray(samples)
|
||||
epochs = array.size // self.integration_samples
|
||||
return [self.process_epoch(array[index * self.integration_samples : (index + 1) * self.integration_samples]) for index in range(epochs)]
|
||||
@@ -0,0 +1,30 @@
|
||||
"""Multi-channel scheduler with separate E1B/E1C logical channels."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy.typing as npt
|
||||
|
||||
from galileo_rx.config import TrackingConfig
|
||||
from galileo_rx.signal_profiles import SignalProfile
|
||||
from galileo_rx.signals.code_provider import CodeProvider
|
||||
from galileo_rx.types import AcquisitionResult, TrackingMeasurement
|
||||
from .tracking_channel import TrackingChannel
|
||||
|
||||
|
||||
class TrackingManager:
|
||||
def __init__(self, profile: SignalProfile, provider: CodeProvider, sampling_frequency_hz: float, config: TrackingConfig) -> None:
|
||||
self.profile = profile
|
||||
self.provider = provider
|
||||
self.sampling_frequency_hz = sampling_frequency_hz
|
||||
self.config = config
|
||||
self.channels: dict[tuple[int, str], TrackingChannel] = {}
|
||||
|
||||
def add_channel(self, acquisition: AcquisitionResult, start_sample_index: int = 0) -> TrackingChannel:
|
||||
key = (acquisition.prn, acquisition.component)
|
||||
channel = TrackingChannel(self.profile, self.provider, self.sampling_frequency_hz, self.config, acquisition, start_sample_index)
|
||||
self.channels[key] = channel
|
||||
return channel
|
||||
|
||||
def process_epoch(self, samples: npt.NDArray, keys: list[tuple[int, str]] | None = None) -> dict[tuple[int, str], TrackingMeasurement]:
|
||||
selected = keys if keys is not None else list(self.channels)
|
||||
return {key: self.channels[key].process_epoch(samples) for key in selected}
|
||||
@@ -0,0 +1,102 @@
|
||||
"""Shared, unit-explicit data transfer objects."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from enum import Enum
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
ComplexArray = npt.NDArray[np.complexfloating[Any, Any]]
|
||||
FloatArray = npt.NDArray[np.floating[Any]]
|
||||
|
||||
|
||||
class TrackingState(str, Enum):
|
||||
IDLE = "IDLE"
|
||||
ACQUIRED = "ACQUIRED"
|
||||
FLL_PULL_IN = "FLL_PULL_IN"
|
||||
FLL_PLL_TRANSITION = "FLL_PLL_TRANSITION"
|
||||
PLL_TRACK = "PLL_TRACK"
|
||||
WEAK_SIGNAL = "WEAK_SIGNAL"
|
||||
LOSS_OF_LOCK = "LOSS_OF_LOCK"
|
||||
REACQUISITION = "REACQUISITION"
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class AcquisitionResult:
|
||||
detected: bool
|
||||
prn: int
|
||||
component: str
|
||||
code_phase_samples: float
|
||||
code_phase_chips: float
|
||||
doppler_hz: float
|
||||
carrier_frequency_hz: float
|
||||
peak_power: float
|
||||
second_peak_power: float
|
||||
detection_metric: float
|
||||
estimated_cn0_db_hz: float | None
|
||||
coherent_integration_ms: float
|
||||
noncoherent_accumulations: int
|
||||
raw_code_phase_samples: int = 0
|
||||
raw_doppler_hz: float = 0.0
|
||||
metrics: dict[str, float] = field(default_factory=dict)
|
||||
side_peak_suspected: bool = False
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class CorrelatorOutputs:
|
||||
very_early: complex = 0j
|
||||
early: complex = 0j
|
||||
prompt: complex = 0j
|
||||
late: complex = 0j
|
||||
very_late: complex = 0j
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class TrackingMeasurement:
|
||||
epoch: int
|
||||
time_s: float
|
||||
state: TrackingState
|
||||
prompt_i: float
|
||||
prompt_q: float
|
||||
early_magnitude: float
|
||||
late_magnitude: float
|
||||
pll_error_rad: float
|
||||
fll_error_hz: float
|
||||
dll_error_chips: float
|
||||
carrier_frequency_hz: float
|
||||
code_frequency_hz: float
|
||||
carrier_phase_rad: float
|
||||
code_phase_chips: float
|
||||
cn0_db_hz: float | None
|
||||
carrier_lock_metric: float
|
||||
code_lock_metric: float
|
||||
cycle_slip: bool = False
|
||||
side_peak_suspected: bool = False
|
||||
reacquisition_requested: bool = False
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class SatelliteTruth:
|
||||
prn: int
|
||||
component: str
|
||||
initial_code_delay_chips: float
|
||||
initial_code_delay_samples: float
|
||||
initial_doppler_hz: float
|
||||
doppler_rate_hz_s: float
|
||||
initial_carrier_phase_rad: float
|
||||
amplitude: float
|
||||
cn0_db_hz: float | None
|
||||
|
||||
|
||||
@dataclass(slots=True)
|
||||
class GenerationTruth:
|
||||
sampling_frequency_hz: float
|
||||
intermediate_frequency_hz: float
|
||||
input_is_real_if: bool
|
||||
start_sample: int
|
||||
sample_count: int
|
||||
satellites: list[SatelliteTruth]
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
"""External-vector verification flows."""
|
||||
@@ -0,0 +1,591 @@
|
||||
"""Executable end-to-end validation against a local globsky/SignalSim tree."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import asdict, replace
|
||||
from hashlib import sha256
|
||||
import json
|
||||
import math
|
||||
from pathlib import Path
|
||||
import subprocess
|
||||
import time
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import yaml
|
||||
|
||||
from galileo_rx.acquisition.acquisition_manager import AcquisitionManager
|
||||
from galileo_rx.config import AcquisitionConfig, TrackingConfig
|
||||
from galileo_rx.diagnostics.plots import (
|
||||
acquisition_map,
|
||||
acquisition_peak_neighborhood,
|
||||
acquisition_profiles,
|
||||
prompt_constellation,
|
||||
tracking_diagnostics,
|
||||
)
|
||||
from galileo_rx.diagnostics.reports import export_tracking_csv
|
||||
from galileo_rx.fixed_point.fixed_acquisition import FixedPcpsAcquisition
|
||||
from galileo_rx.fixed_point.formats import FixedFormat
|
||||
from galileo_rx.io.signalsim_reader import SignalSimReader, SignalSimSampleSpec
|
||||
from galileo_rx.rtl_model.vector_export import export_timed_complex_hex
|
||||
from galileo_rx.signal_profiles import galileo_e1_profile
|
||||
from galileo_rx.signals.galileo_e1 import e1_local_replica
|
||||
from galileo_rx.signals.signalsim_codes import SignalSimCodeProvider
|
||||
from galileo_rx.tracking.correlator import CorrelatorBank
|
||||
from galileo_rx.tracking.nco import CarrierNco, CodeNco
|
||||
from galileo_rx.tracking.symbol_sync import recover_bpsk_symbols, synchronize_secondary_code
|
||||
from galileo_rx.tracking.tracking_channel import TrackingChannel
|
||||
|
||||
|
||||
def locate_signalsim(start: str | Path) -> Path:
|
||||
start_path = Path(start).resolve()
|
||||
roots = [start_path, *start_path.parents[:3], Path.home() / "Desktop", Path.home() / "projects", Path.home() / "src"]
|
||||
names = ("SignalSim", "signalsim", "signal-sim")
|
||||
candidates: list[Path] = []
|
||||
for root in roots:
|
||||
for name in names:
|
||||
candidates.extend((root / name, root / "signalsim" / name))
|
||||
for candidate in candidates:
|
||||
if (candidate / "IFdataGen" / "CMakeLists.txt").is_file() and (candidate / "src" / "PrnGenerate.cpp").is_file():
|
||||
return candidate.resolve()
|
||||
raise FileNotFoundError("could not locate a SignalSim source tree from the workspace or common project directories")
|
||||
|
||||
|
||||
def _git(source: Path, *args: str) -> str:
|
||||
return subprocess.check_output(["git", "-C", str(source), *args], text=True).strip()
|
||||
|
||||
|
||||
def _hash_file(path: Path) -> str:
|
||||
digest = sha256()
|
||||
with path.open("rb") as stream:
|
||||
for block in iter(lambda: stream.read(1024 * 1024), b""):
|
||||
digest.update(block)
|
||||
return digest.hexdigest()
|
||||
|
||||
|
||||
def _save_figure(figure: Any, path: Path) -> None:
|
||||
figure.savefig(path, dpi=140, bbox_inches="tight")
|
||||
import matplotlib.pyplot as plt
|
||||
plt.close(figure)
|
||||
|
||||
|
||||
def _tracking_summary(measurements: list[Any], if_hz: float, expected_doppler_hz: float) -> dict[str, Any]:
|
||||
tail_count = min(200, len(measurements))
|
||||
tail = measurements[-tail_count:]
|
||||
doppler = np.asarray([item.carrier_frequency_hz - if_hz for item in tail], dtype=np.float64)
|
||||
dll = np.asarray([item.dll_error_chips for item in tail], dtype=np.float64)
|
||||
prompt = np.asarray([complex(item.prompt_i, item.prompt_q) for item in tail], dtype=np.complex128)
|
||||
return {
|
||||
"epoch_count": len(measurements),
|
||||
"tail_epoch_count": tail_count,
|
||||
"final_state": measurements[-1].state.value,
|
||||
"reacquisition_epoch_count": sum(item.reacquisition_requested for item in measurements),
|
||||
"loss_of_lock_flag_count": sum(item.state.value in {"LOSS_OF_LOCK", "REACQUISITION"} for item in measurements),
|
||||
"estimated_doppler_hz": float(np.mean(doppler)),
|
||||
"doppler_error_hz": float(np.mean(doppler) - expected_doppler_hz),
|
||||
"doppler_jitter_hz_rms": float(np.std(doppler)),
|
||||
"dll_error_chips_rms": float(np.sqrt(np.mean(dll * dll))),
|
||||
"mean_prompt_magnitude": float(np.mean(np.abs(prompt))),
|
||||
"mean_carrier_lock_metric": float(np.mean([item.carrier_lock_metric for item in tail])),
|
||||
"mean_code_lock_metric": float(np.mean([item.code_lock_metric for item in tail])),
|
||||
"median_cn0_db_hz": float(np.nanmedian([np.nan if item.cn0_db_hz is None else item.cn0_db_hz for item in tail])),
|
||||
}
|
||||
|
||||
|
||||
def validate_signalsim(source_path: str | Path, config_path: str | Path, output_dir: str | Path, *, mode: str | None = None) -> tuple[dict[str, Any], bool]:
|
||||
source = Path(source_path).resolve()
|
||||
config_file = Path(config_path).resolve()
|
||||
project_root = config_file.parent.parent
|
||||
cfg = yaml.safe_load(config_file.read_text(encoding="utf-8"))
|
||||
selected_mode = mode or cfg.get("mode", "quick")
|
||||
if selected_mode not in {"quick", "full"}:
|
||||
raise ValueError("SignalSim validation mode must be quick or full")
|
||||
out = Path(output_dir).resolve()
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
vector_cfg = cfg["vector"]
|
||||
samples_path = project_root / vector_cfg["samples"]
|
||||
if not samples_path.is_file():
|
||||
generator = project_root / "tools" / "generate_signalsim_vectors.py"
|
||||
subprocess.run([str(generator), "--signalsim-path", str(source), "--project-root", str(project_root)], check=True)
|
||||
fs = float(vector_cfg["sampling_frequency_hz"])
|
||||
if_hz = float(vector_cfg["intermediate_frequency_hz"])
|
||||
expected_doppler_hz = float(vector_cfg["expected_doppler_hz"])
|
||||
prn = int(vector_cfg["prn"])
|
||||
reader = SignalSimReader(samples_path, SignalSimSampleSpec(vector_cfg["sample_format"]), memory_map=True)
|
||||
samples = reader.read()
|
||||
profile = galileo_e1_profile()
|
||||
provider = SignalSimCodeProvider(source)
|
||||
acq_cfg = AcquisitionConfig(**cfg["acquisition"], parallelism="doppler_vectorized")
|
||||
manager = AcquisitionManager(profile, provider, fs, acq_cfg)
|
||||
acquisition_samples = manager.engine.coherent_samples * acq_cfg.noncoherent_accumulations
|
||||
acquisition: dict[str, Any] = {}
|
||||
grids: dict[str, Any] = {}
|
||||
for component in ("E1B", "E1C"):
|
||||
started = time.perf_counter()
|
||||
result, grid = manager.acquire_one(samples[:acquisition_samples], component, prn, if_hz)
|
||||
elapsed = time.perf_counter() - started
|
||||
item = asdict(result)
|
||||
item.update({
|
||||
"expected_prn": prn,
|
||||
"expected_doppler_hz": expected_doppler_hz,
|
||||
"raw_doppler_error_hz": result.raw_doppler_hz - expected_doppler_hz,
|
||||
"interpolated_doppler_error_hz": result.doppler_hz - expected_doppler_hz,
|
||||
"expected_code_phase_samples": vector_cfg.get("expected_code_phase_samples"),
|
||||
"code_phase_error_samples": None,
|
||||
"processing_time_s": elapsed,
|
||||
"overflow_count": 0,
|
||||
})
|
||||
acquisition[component] = item
|
||||
grids[component] = grid
|
||||
_save_figure(acquisition_map(grid.power, grid.doppler_bins_hz, f"SignalSim {component} PRN {prn}"), out / f"acquisition_map_{component.lower()}.png")
|
||||
_save_figure(acquisition_profiles(grid.power, grid.doppler_bins_hz), out / f"acquisition_profiles_{component.lower()}.png")
|
||||
_save_figure(acquisition_peak_neighborhood(grid.power, fs / profile.code_rate_hz), out / f"acquisition_peak_neighborhood_{component.lower()}.png")
|
||||
|
||||
e1b_result = manager.acquire_one(samples[:acquisition_samples], "E1B", prn, if_hz)[0]
|
||||
primary_samples = int(round(fs * profile.component("E1B").primary_code_period_s))
|
||||
boundary_start = int(round(e1b_result.code_phase_samples + math.ceil((acquisition_samples - e1b_result.code_phase_samples) / primary_samples) * primary_samples))
|
||||
tracking_cfg = cfg["tracking"]
|
||||
e1b_channel = TrackingChannel(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
TrackingConfig(
|
||||
integration_ms=float(tracking_cfg["e1b_integration_ms"]),
|
||||
pull_in_duration_ms=float(tracking_cfg["e1b_pull_in_ms"]),
|
||||
pll_noise_bandwidth_hz=float(tracking_cfg["pll_noise_bandwidth_hz"]),
|
||||
fll_noise_bandwidth_hz=float(tracking_cfg["fll_noise_bandwidth_hz"]),
|
||||
dll_noise_bandwidth_hz=float(tracking_cfg["dll_noise_bandwidth_hz"]),
|
||||
secondary_code_mode="none",
|
||||
loss_of_lock_epochs=40,
|
||||
),
|
||||
e1b_result,
|
||||
boundary_start,
|
||||
)
|
||||
e1b_measurements = e1b_channel.process(samples[boundary_start:])
|
||||
e1b_summary = _tracking_summary(e1b_measurements, if_hz, expected_doppler_hz)
|
||||
export_tracking_csv(out / "tracking_e1b.csv", e1b_measurements)
|
||||
_save_figure(tracking_diagnostics(e1b_measurements), out / "tracking_e1b.png")
|
||||
_save_figure(prompt_constellation(e1b_measurements, ignored_epochs=150), out / "prompt_constellation_e1b.png")
|
||||
|
||||
ignored_epochs = int(round(0.300 / (float(tracking_cfg["e1b_integration_ms"]) * 1e-3)))
|
||||
prompt_e1b = np.asarray([complex(item.prompt_i, item.prompt_q) for item in e1b_measurements])
|
||||
symbols, soft_symbols = recover_bpsk_symbols(prompt_e1b, 2, ignored_epochs=ignored_epochs)
|
||||
sync_bits = np.asarray([0, 1, 0, 1, 1, 0, 0, 0, 0, 0], dtype=np.int8)
|
||||
sync_signs = (2 * sync_bits - 1).astype(np.int8)
|
||||
sync_matches: list[dict[str, Any]] = []
|
||||
for inversion in (1, -1):
|
||||
for index in range(max(0, symbols.size - sync_signs.size + 1)):
|
||||
errors = int(np.count_nonzero(inversion * symbols[index : index + 10] != sync_signs))
|
||||
if errors == 0:
|
||||
sync_matches.append({"symbol_index": index, "sign_inverted": inversion == -1, "errors": errors})
|
||||
symbol_report = {
|
||||
"status": "raw_symbols_recovered_page_decoding_not_implemented",
|
||||
"comparison_scope": "known_10_symbol_INAV_sync_words_only",
|
||||
"recovered_symbols": int(symbols.size),
|
||||
"compared_symbols": len(sync_matches) * 10,
|
||||
"initial_ignored_tracking_epochs": ignored_epochs,
|
||||
"initial_ignored_time_s": ignored_epochs * float(tracking_cfg["e1b_integration_ms"]) * 1e-3,
|
||||
"bounded_alignment": "code boundary from PCPS; sync search within remaining 1.876 s",
|
||||
"sync_matches": sync_matches,
|
||||
"symbol_errors": 0 if sync_matches else None,
|
||||
"symbol_error_rate": 0.0 if sync_matches else None,
|
||||
"median_soft_symbol_magnitude": float(np.median(np.abs(soft_symbols))),
|
||||
"frame_or_page_sync": len(sync_matches) >= 2,
|
||||
"fec_decoding": "not_implemented",
|
||||
"crc_validation": "not_implemented",
|
||||
}
|
||||
|
||||
e1c_result = manager.acquire_one(samples[:acquisition_samples], "E1C", prn, if_hz)[0]
|
||||
e1c_boundary = int(round(e1c_result.code_phase_samples + math.ceil((acquisition_samples - e1c_result.code_phase_samples) / primary_samples) * primary_samples))
|
||||
open_epoch_count = min(100, (samples.size - e1c_boundary) // primary_samples)
|
||||
carrier = CarrierNco(fs, e1c_result.carrier_frequency_hz, 2 * np.pi * e1c_result.carrier_frequency_hz * e1c_boundary / fs % (2 * np.pi))
|
||||
code = CodeNco(
|
||||
fs,
|
||||
profile.component("E1C").primary_code_length_chips,
|
||||
profile.code_rate_hz * (1 + e1c_result.doppler_hz / profile.carrier_frequency_hz),
|
||||
(-e1c_result.code_phase_chips + e1c_boundary * profile.code_rate_hz / fs) % 4092,
|
||||
)
|
||||
bank = CorrelatorBank(profile, provider, "E1C", prn, carrier, code, 0.25, 0.75)
|
||||
open_prompts = np.asarray([
|
||||
bank.correlate(samples[e1c_boundary + epoch * primary_samples : e1c_boundary + (epoch + 1) * primary_samples]).prompt
|
||||
for epoch in range(open_epoch_count)
|
||||
], dtype=np.complex128)
|
||||
secondary = provider.secondary_code(profile, "E1C", prn)
|
||||
assert secondary is not None
|
||||
sync = synchronize_secondary_code(open_prompts, secondary, 0.004, residual_frequency_step_hz=0.5)
|
||||
epochs = np.arange(min(25, open_prompts.size), dtype=np.float64)
|
||||
corrected = open_prompts[: epochs.size] * np.exp(-2j * np.pi * sync.residual_frequency_hz * 0.004 * epochs)
|
||||
wipe = secondary[(sync.phase + np.arange(epochs.size)) % secondary.size]
|
||||
wiped_gain_db = float(10 * np.log10(abs(np.sum(corrected * wipe)) ** 2 / (epochs.size * np.mean(np.abs(corrected) ** 2))))
|
||||
unwiped_gain_db = float(10 * np.log10(max(abs(np.sum(corrected)) ** 2, np.finfo(np.float64).tiny) / (epochs.size * np.mean(np.abs(corrected) ** 2))))
|
||||
adjusted = replace(
|
||||
e1c_result,
|
||||
doppler_hz=e1c_result.doppler_hz + sync.residual_frequency_hz,
|
||||
carrier_frequency_hz=e1c_result.carrier_frequency_hz + sync.residual_frequency_hz,
|
||||
)
|
||||
e1c_channel = TrackingChannel(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
TrackingConfig(
|
||||
integration_ms=float(tracking_cfg["e1c_integration_ms"]),
|
||||
pull_in_duration_ms=float(tracking_cfg["e1c_pull_in_ms"]),
|
||||
pll_noise_bandwidth_hz=float(tracking_cfg["pll_noise_bandwidth_hz"]),
|
||||
fll_noise_bandwidth_hz=float(tracking_cfg["fll_noise_bandwidth_hz"]),
|
||||
dll_noise_bandwidth_hz=float(tracking_cfg["dll_noise_bandwidth_hz"]),
|
||||
secondary_code_mode="wipeoff",
|
||||
loss_of_lock_epochs=20,
|
||||
),
|
||||
adjusted,
|
||||
e1c_boundary,
|
||||
secondary_code_phase_at_start=sync.phase,
|
||||
)
|
||||
e1c_measurements = e1c_channel.process(samples[e1c_boundary:])
|
||||
e1c_summary = _tracking_summary(e1c_measurements, if_hz, expected_doppler_hz)
|
||||
export_tracking_csv(out / "tracking_e1c.csv", e1c_measurements)
|
||||
_save_figure(tracking_diagnostics(e1c_measurements), out / "tracking_e1c.png")
|
||||
_save_figure(prompt_constellation(e1c_measurements, ignored_epochs=75), out / "prompt_constellation_e1c.png")
|
||||
secondary_report = asdict(sync) | {
|
||||
"sequence_source": str(provider.pilot_bit_path),
|
||||
"sequence_classification": "SignalSim_implemented_sequence_not_independently_certified_against_ICD",
|
||||
"wiped_25_period_coherent_gain_db": wiped_gain_db,
|
||||
"unwiped_25_period_coherent_gain_db": unwiped_gain_db,
|
||||
"wipeoff_gain_improvement_db": wiped_gain_db - unwiped_gain_db,
|
||||
}
|
||||
|
||||
# Bit-true boundary comparison on the same external IQ8 coherent block.
|
||||
# The floating reference remains float64/complex128; only the explicit
|
||||
# input/FFT/power boundaries are quantized.
|
||||
fixed_samples = samples[:primary_samples] / np.float64(128.0)
|
||||
fixed_bins_hz = np.arange(1500.0, 2125.0 + 0.5 * 125.0, 125.0, dtype=np.float64)
|
||||
fixed_replica = e1_local_replica(
|
||||
profile,
|
||||
provider,
|
||||
"E1C",
|
||||
prn,
|
||||
primary_samples,
|
||||
fs,
|
||||
code_rate_hz=profile.code_rate_hz * (1.0 + 1750.0 / profile.carrier_frequency_hz),
|
||||
)
|
||||
fixed_grid = FixedPcpsAcquisition(
|
||||
FixedFormat(10, 8),
|
||||
FixedFormat(32, 20),
|
||||
FixedFormat(64, 32, signed=False),
|
||||
).search(fixed_samples, fixed_replica, fs, if_hz + fixed_bins_hz)
|
||||
time_axis_s = np.arange(primary_samples, dtype=np.float64) / fs
|
||||
reference_power = np.asarray([
|
||||
np.abs(
|
||||
np.fft.ifft(
|
||||
np.fft.fft(fixed_samples * np.exp(-2j * np.pi * (if_hz + doppler_hz) * time_axis_s))
|
||||
* np.conjugate(np.fft.fft(fixed_replica))
|
||||
)
|
||||
) ** 2
|
||||
for doppler_hz in fixed_bins_hz
|
||||
], dtype=np.float64)
|
||||
reference_peak = np.unravel_index(int(np.argmax(reference_power)), reference_power.shape)
|
||||
reference_normalized = reference_power / np.max(reference_power)
|
||||
fixed_normalized = fixed_grid.power / np.max(fixed_grid.power)
|
||||
fixed_error = fixed_normalized - reference_normalized
|
||||
fixed_report = {
|
||||
"reference_dtype": "complex128/float64",
|
||||
"input_format": {"total_bits": 10, "fractional_bits": 8},
|
||||
"fft_boundary_format": {"total_bits": 32, "fractional_bits": 20},
|
||||
"power_format": {"total_bits": 64, "fractional_bits": 32, "signed": False},
|
||||
"reference_doppler_bin_hz": float(fixed_bins_hz[reference_peak[0]]),
|
||||
"fixed_doppler_bin_hz": float(fixed_bins_hz[fixed_grid.peak_doppler_index]),
|
||||
"reference_code_phase_samples": int(reference_peak[1]),
|
||||
"fixed_code_phase_samples": fixed_grid.peak_code_index,
|
||||
"doppler_bin_agreement": bool(fixed_grid.peak_doppler_index == reference_peak[0]),
|
||||
"code_phase_agreement": bool(fixed_grid.peak_code_index == reference_peak[1]),
|
||||
"acquisition_decision_agreement": bool((np.max(reference_power) / np.mean(reference_power) >= 2.5) == (np.max(fixed_grid.power) / np.mean(fixed_grid.power) >= 2.5)),
|
||||
"normalized_absolute_error_mean": float(np.mean(np.abs(fixed_error))),
|
||||
"normalized_rms_error": float(np.sqrt(np.mean(fixed_error * fixed_error))),
|
||||
"normalized_peak_error": float(np.max(np.abs(fixed_error))),
|
||||
"overflow_count": fixed_grid.overflow_count,
|
||||
}
|
||||
|
||||
carrier_fcw = int(round((if_hz + expected_doppler_hz) / fs * 2**32)) & 0xFFFFFFFF
|
||||
code_fcw = int(round(profile.code_rate_hz * (1 + expected_doppler_hz / profile.carrier_frequency_hz) / fs * 2**48))
|
||||
rtl_registers = {
|
||||
"SAMPLE_WIDTH": 8,
|
||||
"SAMPLE_FRACTIONAL_BITS": 0,
|
||||
"INPUT_REPRESENTATION": "COMPLEX_IQ",
|
||||
"IQ_ORDER": "IQ",
|
||||
"SAMPLING_FREQUENCY_HZ": fs,
|
||||
"INTERMEDIATE_FREQUENCY_HZ": if_hz,
|
||||
"PRN": prn,
|
||||
"COMPONENT": "E1C",
|
||||
"ACQ_DOPPLER_MIN_HZ": acq_cfg.doppler_min_hz,
|
||||
"ACQ_DOPPLER_MAX_HZ": acq_cfg.doppler_max_hz,
|
||||
"ACQ_DOPPLER_STEP_HZ": acq_cfg.doppler_step_hz,
|
||||
"ACQ_COHERENT_SAMPLES": manager.engine.coherent_samples,
|
||||
"ACQ_NONCOHERENT_COUNT": acq_cfg.noncoherent_accumulations,
|
||||
"ACQ_FFT_LENGTH": manager.engine.coherent_samples,
|
||||
"TRACK_E1B_INTEGRATION_SAMPLES": int(round(fs * float(tracking_cfg["e1b_integration_ms"]) * 1e-3)),
|
||||
"TRACK_E1C_INTEGRATION_SAMPLES": primary_samples,
|
||||
"CARRIER_NCO_PHASE_WIDTH": 32,
|
||||
"CARRIER_NCO_FCW": carrier_fcw,
|
||||
"CODE_NCO_PHASE_WIDTH": 48,
|
||||
"CODE_NCO_FCW": code_fcw,
|
||||
"E1C_SECONDARY_PHASE": sync.phase,
|
||||
"ACQ_TO_TRACK_CODE_PHASE_SAMPLES": acquisition["E1C"]["code_phase_samples"],
|
||||
"ACQ_TO_TRACK_DOPPLER_HZ": expected_doppler_hz,
|
||||
"NOTE": "Symbolic names; APB addresses remain to be assigned by the RTL IP register map",
|
||||
}
|
||||
rtl_manifest = export_timed_complex_hex(
|
||||
out / "rtl_vectors",
|
||||
"signalsim_acquisition_input_iq8",
|
||||
samples[:acquisition_samples],
|
||||
FixedFormat(8, 0),
|
||||
fs,
|
||||
rtl_registers,
|
||||
)
|
||||
|
||||
criteria = cfg["pass_criteria"]
|
||||
acquisition_pass = all(
|
||||
acquisition[component]["detected"]
|
||||
and acquisition[component]["prn"] == prn
|
||||
and abs(acquisition[component]["raw_doppler_error_hz"]) <= criteria["maximum_raw_doppler_error_hz"]
|
||||
for component in ("E1B", "E1C")
|
||||
)
|
||||
def tracking_pass(summary: dict[str, Any]) -> bool:
|
||||
return (
|
||||
summary["final_state"] == "PLL_TRACK"
|
||||
and summary["reacquisition_epoch_count"] == 0
|
||||
and abs(summary["doppler_error_hz"]) <= criteria["maximum_tracking_doppler_error_hz"]
|
||||
and summary["doppler_jitter_hz_rms"] <= criteria["maximum_tracking_doppler_jitter_hz"]
|
||||
and summary["dll_error_chips_rms"] <= criteria["maximum_dll_rms_chips"]
|
||||
)
|
||||
full_vector_report: dict[str, Any] = {}
|
||||
full_vector_pass = True
|
||||
if selected_mode == "full":
|
||||
definitions = {
|
||||
"nominal_cn0_single_prn": {1: 1828.0},
|
||||
"negative_doppler": {27: -2637.0},
|
||||
"multiple_prns": {1: 1828.0, 4: 2171.0, 19: -367.0, 21: -712.0, 27: -2637.0},
|
||||
}
|
||||
generator = project_root / "tools" / "generate_signalsim_vectors.py"
|
||||
for vector_name, expected_by_prn in definitions.items():
|
||||
metadata_path = project_root / "verification" / "vectors" / "signalsim" / vector_name / "metadata.yaml"
|
||||
if not metadata_path.is_file():
|
||||
subprocess.run([str(generator), "--signalsim-path", str(source), "--project-root", str(project_root), "--vector", vector_name], check=True)
|
||||
metadata = yaml.safe_load(metadata_path.read_text(encoding="utf-8"))
|
||||
path = project_root / metadata["output"]
|
||||
vector_samples = SignalSimReader(path, SignalSimSampleSpec("IQ8"), memory_map=True).read(acquisition_samples)
|
||||
broad_cfg = AcquisitionConfig(
|
||||
doppler_min_hz=-4000.0,
|
||||
doppler_max_hz=3000.0,
|
||||
doppler_step_hz=250.0,
|
||||
coherent_integration_ms=4.0,
|
||||
noncoherent_accumulations=5,
|
||||
detection_threshold=2.5,
|
||||
zero_padding_mode="none",
|
||||
parallelism="doppler_vectorized",
|
||||
)
|
||||
broad_manager = AcquisitionManager(profile, provider, fs, broad_cfg)
|
||||
vector_results: dict[int, Any] = {}
|
||||
for vector_prn, vector_truth_hz in expected_by_prn.items():
|
||||
vector_result, _ = broad_manager.acquire_one(vector_samples, "E1C", vector_prn, if_hz)
|
||||
item = asdict(vector_result) | {
|
||||
"expected_doppler_hz": vector_truth_hz,
|
||||
"raw_doppler_error_hz": vector_result.raw_doppler_hz - vector_truth_hz,
|
||||
}
|
||||
item["pass"] = bool(vector_result.detected and abs(item["raw_doppler_error_hz"]) <= 250.0)
|
||||
full_vector_pass &= item["pass"]
|
||||
vector_results[vector_prn] = item
|
||||
full_vector_report[vector_name] = {
|
||||
"sha256": _hash_file(path),
|
||||
"size_bytes": path.stat().st_size,
|
||||
"acquisition": vector_results,
|
||||
}
|
||||
|
||||
# Negative-Doppler tracking is checked over the long vector. Transient
|
||||
# lock-detector assertions during pull-in are reported, while the pass
|
||||
# condition applies after convergence.
|
||||
negative_path = project_root / yaml.safe_load(
|
||||
(project_root / "verification" / "vectors" / "signalsim" / "negative_doppler" / "metadata.yaml").read_text(encoding="utf-8")
|
||||
)["output"]
|
||||
negative_samples = SignalSimReader(negative_path, SignalSimSampleSpec("IQ8"), memory_map=True).read()
|
||||
negative_acq_manager = AcquisitionManager(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
AcquisitionConfig(doppler_min_hz=-4000.0, doppler_max_hz=-1500.0, doppler_step_hz=250.0, noncoherent_accumulations=5, zero_padding_mode="none"),
|
||||
)
|
||||
negative_acq, _ = negative_acq_manager.acquire_one(negative_samples[:acquisition_samples], "E1B", 27, if_hz)
|
||||
negative_start = int(round(negative_acq.code_phase_samples + math.ceil((acquisition_samples - negative_acq.code_phase_samples) / primary_samples) * primary_samples))
|
||||
negative_channel = TrackingChannel(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
TrackingConfig(integration_ms=2.0, pull_in_duration_ms=300.0, pll_noise_bandwidth_hz=10.0, fll_noise_bandwidth_hz=20.0, dll_noise_bandwidth_hz=1.0, secondary_code_mode="none", loss_of_lock_epochs=40),
|
||||
negative_acq,
|
||||
negative_start,
|
||||
)
|
||||
negative_measurements = negative_channel.process(negative_samples[negative_start:])
|
||||
negative_summary = _tracking_summary(negative_measurements, if_hz, -2637.0)
|
||||
post_pull_in = negative_measurements[500:]
|
||||
negative_summary["post_pull_in_reacquisition_epoch_count"] = sum(item.reacquisition_requested for item in post_pull_in)
|
||||
negative_tracking_pass = bool(
|
||||
negative_summary["final_state"] == "PLL_TRACK"
|
||||
and negative_summary["post_pull_in_reacquisition_epoch_count"] == 0
|
||||
and abs(negative_summary["doppler_error_hz"]) <= criteria["maximum_tracking_doppler_error_hz"]
|
||||
and negative_summary["doppler_jitter_hz_rms"] <= criteria["maximum_tracking_doppler_jitter_hz"]
|
||||
and negative_summary["dll_error_chips_rms"] <= criteria["maximum_dll_rms_chips"]
|
||||
)
|
||||
negative_summary["pass"] = negative_tracking_pass
|
||||
full_vector_report["negative_doppler"]["tracking_e1b"] = negative_summary
|
||||
full_vector_pass &= negative_tracking_pass
|
||||
|
||||
# Nominal-C/N0 uses the pilot because secondary wipeoff removes the
|
||||
# 100 ms overlay transition from the carrier loop.
|
||||
nominal_path = project_root / yaml.safe_load(
|
||||
(project_root / "verification" / "vectors" / "signalsim" / "nominal_cn0_single_prn" / "metadata.yaml").read_text(encoding="utf-8")
|
||||
)["output"]
|
||||
nominal_samples = SignalSimReader(nominal_path, SignalSimSampleSpec("IQ8"), memory_map=True).read()
|
||||
nominal_acq_manager = AcquisitionManager(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
AcquisitionConfig(doppler_min_hz=1000.0, doppler_max_hz=2500.0, doppler_step_hz=250.0, noncoherent_accumulations=5, zero_padding_mode="none"),
|
||||
)
|
||||
nominal_acq, _ = nominal_acq_manager.acquire_one(nominal_samples[:acquisition_samples], "E1C", 1, if_hz)
|
||||
nominal_start = int(round(nominal_acq.code_phase_samples + math.ceil((acquisition_samples - nominal_acq.code_phase_samples) / primary_samples) * primary_samples))
|
||||
nominal_open_count = min(100, (nominal_samples.size - nominal_start) // primary_samples)
|
||||
nominal_carrier = CarrierNco(fs, nominal_acq.carrier_frequency_hz, 2 * np.pi * nominal_acq.carrier_frequency_hz * nominal_start / fs % (2 * np.pi))
|
||||
nominal_code = CodeNco(
|
||||
fs,
|
||||
4092,
|
||||
profile.code_rate_hz * (1 + nominal_acq.doppler_hz / profile.carrier_frequency_hz),
|
||||
(-nominal_acq.code_phase_chips + nominal_start * profile.code_rate_hz / fs) % 4092,
|
||||
)
|
||||
nominal_bank = CorrelatorBank(profile, provider, "E1C", 1, nominal_carrier, nominal_code, 0.25, 0.75)
|
||||
nominal_prompts = np.asarray([
|
||||
nominal_bank.correlate(nominal_samples[nominal_start + epoch * primary_samples : nominal_start + (epoch + 1) * primary_samples]).prompt
|
||||
for epoch in range(nominal_open_count)
|
||||
], dtype=np.complex128)
|
||||
nominal_sync = synchronize_secondary_code(nominal_prompts, secondary, 0.004, residual_frequency_step_hz=0.5)
|
||||
nominal_adjusted = replace(
|
||||
nominal_acq,
|
||||
doppler_hz=nominal_acq.doppler_hz + nominal_sync.residual_frequency_hz,
|
||||
carrier_frequency_hz=nominal_acq.carrier_frequency_hz + nominal_sync.residual_frequency_hz,
|
||||
)
|
||||
nominal_channel = TrackingChannel(
|
||||
profile,
|
||||
provider,
|
||||
fs,
|
||||
TrackingConfig(integration_ms=4.0, pull_in_duration_ms=200.0, pll_noise_bandwidth_hz=10.0, fll_noise_bandwidth_hz=10.0, dll_noise_bandwidth_hz=1.0, secondary_code_mode="wipeoff", loss_of_lock_epochs=20),
|
||||
nominal_adjusted,
|
||||
nominal_start,
|
||||
secondary_code_phase_at_start=nominal_sync.phase,
|
||||
)
|
||||
nominal_measurements = nominal_channel.process(nominal_samples[nominal_start:])
|
||||
nominal_summary = _tracking_summary(nominal_measurements, if_hz, 1828.0)
|
||||
nominal_summary["secondary_sync"] = asdict(nominal_sync)
|
||||
nominal_summary["post_pull_in_reacquisition_epoch_count"] = sum(item.reacquisition_requested for item in nominal_measurements[75:])
|
||||
nominal_tracking_pass = bool(
|
||||
nominal_summary["final_state"] == "PLL_TRACK"
|
||||
and nominal_summary["post_pull_in_reacquisition_epoch_count"] == 0
|
||||
and abs(nominal_summary["doppler_error_hz"]) <= criteria["maximum_tracking_doppler_error_hz"]
|
||||
and nominal_summary["doppler_jitter_hz_rms"] <= criteria["maximum_tracking_doppler_jitter_hz"]
|
||||
and nominal_summary["dll_error_chips_rms"] <= criteria["maximum_dll_rms_chips"]
|
||||
and nominal_sync.normalized_correlation >= criteria["minimum_secondary_correlation"]
|
||||
and nominal_sync.sign_errors == 0
|
||||
)
|
||||
nominal_summary["pass"] = nominal_tracking_pass
|
||||
full_vector_report["nominal_cn0_single_prn"]["tracking_e1c"] = nominal_summary
|
||||
full_vector_pass &= nominal_tracking_pass
|
||||
|
||||
passes = {
|
||||
"acquisition": acquisition_pass,
|
||||
"e1b_tracking": tracking_pass(e1b_summary),
|
||||
"e1b_raw_symbol_sync": len(sync_matches) >= int(criteria["minimum_inav_sync_words"]),
|
||||
"e1c_secondary_sync": sync.normalized_correlation >= criteria["minimum_secondary_correlation"] and sync.sign_errors <= criteria["maximum_secondary_sign_errors"],
|
||||
"e1c_tracking": tracking_pass(e1c_summary),
|
||||
"fixed_point_acquisition": bool(
|
||||
fixed_report["doppler_bin_agreement"]
|
||||
and fixed_report["code_phase_agreement"]
|
||||
and fixed_report["acquisition_decision_agreement"]
|
||||
and fixed_report["overflow_count"] == 0
|
||||
),
|
||||
}
|
||||
if selected_mode == "full":
|
||||
passes["full_vector_set"] = full_vector_pass
|
||||
report = {
|
||||
"simulator": {
|
||||
"name": "SignalSim",
|
||||
"path": str(source),
|
||||
"remote": _git(source, "remote", "get-url", "origin"),
|
||||
"commit": _git(source, "rev-parse", "HEAD"),
|
||||
"git_status": _git(source, "status", "--short", "--branch"),
|
||||
},
|
||||
"vector": {
|
||||
**vector_cfg,
|
||||
"path": str(samples_path),
|
||||
"size_bytes": samples_path.stat().st_size,
|
||||
"complex_sample_count": reader.sample_count,
|
||||
"sha256": _hash_file(samples_path),
|
||||
"internal_dtype": "complex128",
|
||||
},
|
||||
"mode": selected_mode,
|
||||
"acquisition": acquisition,
|
||||
"tracking": {"E1B": e1b_summary, "E1C": e1c_summary},
|
||||
"e1b_symbols": symbol_report,
|
||||
"e1c_secondary": secondary_report,
|
||||
"fixed_point_acquisition": fixed_report,
|
||||
"rtl_vector_manifest": rtl_manifest,
|
||||
"full_vectors": full_vector_report,
|
||||
"pass_criteria": criteria,
|
||||
"passes": passes,
|
||||
"overall_pass": all(passes.values()),
|
||||
"known_failures": [
|
||||
"SignalSim E1 output cannot select E1B-only or E1C-only",
|
||||
"SignalSim E1 uses BOC(1,1), not full CBOC",
|
||||
"SignalSim has no native real-IF output",
|
||||
"SignalSim IQ2 writer omits the output-byte assignment at this commit",
|
||||
"Full I/NAV receiver-side FEC/page/CRC decoding is not implemented",
|
||||
"SignalSim does not emit exact sample-side code-delay or navigation-symbol truth",
|
||||
],
|
||||
}
|
||||
(out / "signalsim_validation_report.json").write_text(json.dumps(report, indent=2, default=str) + "\n", encoding="utf-8")
|
||||
(out / "signalsim_e1b_symbol_comparison.json").write_text(json.dumps(symbol_report, indent=2) + "\n", encoding="utf-8")
|
||||
reports_root = project_root / "reports"
|
||||
reports_root.mkdir(parents=True, exist_ok=True)
|
||||
acquisition_report = {
|
||||
"simulator": report["simulator"],
|
||||
"vector": report["vector"],
|
||||
"acquisition": acquisition,
|
||||
"pass": acquisition_pass,
|
||||
"truth_limitations": "SignalSim CLI logs rounded Doppler but does not export exact code delay",
|
||||
}
|
||||
(reports_root / "signalsim_acquisition_report.json").write_text(json.dumps(acquisition_report, indent=2, default=str) + "\n", encoding="utf-8")
|
||||
(reports_root / "signalsim_e1b_symbol_comparison.json").write_text(json.dumps(symbol_report, indent=2) + "\n", encoding="utf-8")
|
||||
markdown = [
|
||||
"# SignalSim acquisition and tracking validation",
|
||||
"",
|
||||
f"Overall {selected_mode}-mode result: **{'PASS' if report['overall_pass'] else 'FAIL'}**",
|
||||
"",
|
||||
"| Check | Result |",
|
||||
"|---|---:|",
|
||||
*[f"| {name} | {'PASS' if value else 'FAIL'} |" for name, value in passes.items()],
|
||||
"",
|
||||
f"E1B tracking Doppler: {e1b_summary['estimated_doppler_hz']:.3f} Hz (error {e1b_summary['doppler_error_hz']:+.3f} Hz).",
|
||||
f"E1C tracking Doppler: {e1c_summary['estimated_doppler_hz']:.3f} Hz (error {e1c_summary['doppler_error_hz']:+.3f} Hz).",
|
||||
f"E1C secondary phase: {sync.phase}; normalized correlation {sync.normalized_correlation:.4f}; sign errors {sync.sign_errors}/{sync.compared_periods}.",
|
||||
f"E1B exact I/NAV sync words: {len(sync_matches)}. Full page decoding is not implemented.",
|
||||
]
|
||||
(out / "signalsim_validation_report.md").write_text("\n".join(markdown) + "\n", encoding="utf-8")
|
||||
acquisition_markdown = [
|
||||
"# SignalSim acquisition report",
|
||||
"",
|
||||
"| Component | Detected | Raw Doppler (Hz) | Truth (Hz) | Code phase (samples) | Peak/second |",
|
||||
"|---|---:|---:|---:|---:|---:|",
|
||||
*[
|
||||
f"| {component} | {item['detected']} | {item['raw_doppler_hz']:.1f} | {expected_doppler_hz:.1f} | {item['code_phase_samples']:.3f} | {item['detection_metric']:.3f} |"
|
||||
for component, item in acquisition.items()
|
||||
],
|
||||
"",
|
||||
"SignalSim does not export exact code-delay truth. Agreement between independently acquired E1B and E1C code phases is used as a cross-component consistency check; no code-phase error is fabricated.",
|
||||
]
|
||||
(reports_root / "signalsim_acquisition_report.md").write_text("\n".join(acquisition_markdown) + "\n", encoding="utf-8")
|
||||
return report, bool(report["overall_pass"])
|
||||
Reference in New Issue
Block a user