"""CLOCK-001: программная проверка QPSK при независимых тактах TX и RX.""" from __future__ import annotations import argparse from dataclasses import dataclass import hashlib import json import math from pathlib import Path import statistics import sys import time import numpy as np if __package__ in (None, ""): sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from protocol.packet import parse_packet from protocol.video003_modem import SAMPLE_RATE_HZ, build_frame from protocol.video005_iq_waveform import ( GUARD_SAMPLES, NOMINAL_CARRIER_HZ, SYMBOL_RATE, rrc_taps, shape_packet, ) from protocol.video005_receiver import IqStreamResult, receive_iq_stream from tools.qpsk_rf_common import GAP_SAMPLES, SEED ROOT = Path(__file__).resolve().parents[1] OUTPUT = ROOT / "data/processed/clock001/results.json" PAYLOAD_SIZES = (256, 512, 1012) OFFSETS_PPM = (0, -1, 1, -5, 5, -10, 10, -20, 20, -50, 50, -100, 100) PRIMARY_LIMIT_PPM = 50 PACKETS_PER_CASE = 3 @dataclass(frozen=True) class ClockFixture: iq: np.ndarray expected_packets: dict[int, bytes] starts: dict[int, int] symbol_counts: dict[int, int] payload_bytes: int def payload_for_case(payload_bytes: int, sequence_number: int) -> bytes: """Расширяет фиксированную RF-001 SHA-последовательность до длины пакета.""" if not 0 <= sequence_number <= 0xFFFF or not 1 <= payload_bytes <= 1012: raise ValueError("неверный номер или размер пакета") prefix = SEED.to_bytes(8, "big") + sequence_number.to_bytes(2, "big") blocks = math.ceil(payload_bytes / hashlib.sha256().digest_size) return b"".join( hashlib.sha256(prefix + block.to_bytes(2, "big")).digest() for block in range(blocks) )[:payload_bytes] def build_fixture(payload_bytes: int, packet_count: int = PACKETS_PER_CASE) -> ClockFixture: """Собирает поток с RF-001 паузой; границы сохраняет только для оценки.""" if packet_count < 1 or packet_count > 100: raise ValueError("число пакетов должно быть от 1 до 100") separator = GAP_SAMPLES - 2 * GUARD_SAMPLES if separator < 0: raise AssertionError("пауза короче защитных интервалов") waveforms = [] expected_packets = {} starts = {} symbol_counts = {} position = 0 for sequence_number in range(packet_count): frame = build_frame(payload_for_case(payload_bytes, sequence_number), "QPSK", "2/3", sequence_number=sequence_number) waveform = shape_packet(frame) starts[sequence_number] = position expected_packets[sequence_number] = frame.packet symbol_counts[sequence_number] = len(frame.symbols) waveforms.append(waveform) position += len(waveform) + (separator if sequence_number < packet_count - 1 else 0) iq = np.zeros(position, dtype=np.complex64) for sequence_number, waveform in enumerate(waveforms): start = starts[sequence_number] iq[start:start + len(waveform)] = waveform return ClockFixture(iq, expected_packets, starts, symbol_counts, payload_bytes) def sample_independent_clock( tx_iq: np.ndarray, sample_clock_offset_ppm: float, *, cfo_hz: float = 0.0, initial_phase_degrees: float = 0.0, ebn0_db: float | None = None, seed: int = SEED, ) -> np.ndarray: """Выбирает TX-волну в моменты АЦП Fs_RX=Fs_TX·(1+ppm·1e-6). Дробная интерполяция меняет временную шкалу по всему потоку. CFO задаётся отдельно в герцах физического времени RX; затем добавляется AWGN. """ if (tx_iq.ndim != 1 or not np.iscomplexobj(tx_iq) or tx_iq.size < 2 or not np.all(np.isfinite(tx_iq))): raise ValueError("нужен конечный одномерный комплексный TX IQ") if not all(math.isfinite(value) for value in ( sample_clock_offset_ppm, cfo_hz, initial_phase_degrees)): raise ValueError("параметры такта и несущей должны быть конечными") if abs(sample_clock_offset_ppm) > 1000: raise ValueError("модель ограничена диапазоном ±1000 ppm") if ebn0_db is not None and not math.isfinite(ebn0_db): raise ValueError("Eb/N0 должно быть конечным") scale = 1 + sample_clock_offset_ppm * 1e-6 output_count = math.floor((tx_iq.size - 1) * scale) + 1 received = np.empty(output_count, dtype=np.complex64) # Обрабатываем длинные потоки частями, сохраняя единую временную шкалу. for start in range(0, output_count, 1_000_000): stop = min(start + 1_000_000, output_count) positions = np.arange(start, stop, dtype=np.float64) / scale lower = np.floor(positions).astype(np.int64) fraction = positions - lower upper = np.minimum(lower + 1, tx_iq.size - 1) part = (tx_iq[lower] * (1 - fraction) + tx_iq[upper] * fraction).astype(np.complex64) if cfo_hz or initial_phase_degrees: physical_time = np.arange(start, stop, dtype=np.float64) / (SAMPLE_RATE_HZ * scale) phase = math.radians(initial_phase_degrees) + 2 * math.pi * cfo_hz * physical_time part *= np.exp(1j * phase).astype(np.complex64) received[start:stop] = part if ebn0_db is not None: sigma = math.sqrt(1 / (2 * (2 / 3) * 10 ** (ebn0_db / 10)) / 2) random_generator = np.random.default_rng(seed) for start in range(0, output_count, 1_000_000): stop = min(start + 1_000_000, output_count) received[start:stop] += (sigma * ( random_generator.standard_normal(stop - start) + 1j * random_generator.standard_normal(stop - start) )).astype(np.complex64) return received def _median(values: list[float]) -> float | None: return float(statistics.median(values)) if values else None def summarize_case( result: IqStreamResult, fixture: ClockFixture, *, sample_clock_offset_ppm: float, cfo_hz: float, initial_phase_degrees: float, ebn0_db: float | None, receiver_time_s: float, received_samples: int, ) -> dict: """Сравнивает строгий выход приёмника с эталоном после декодирования.""" scale = 1 + sample_clock_offset_ppm * 1e-6 tx_filter_delay = (len(rrc_taps()) - 1) / 2 rx_filter_delay = tx_filter_delay samples_per_symbol = SAMPLE_RATE_HZ // SYMBOL_RATE effective_cfo_hz = (NOMINAL_CARRIER_HZ + cfo_hz) / scale - NOMINAL_CARRIER_HZ exact_sequences: set[int] = set() seen: set[int] = set() duplicates: list[int] = [] unexpected: list[int] = [] initial_errors = [] final_errors = [] drift_values = [] cfo_errors = [] for packet in result.packets: sequence_number = parse_packet(packet.packet).sequence_number if sequence_number in seen: duplicates.append(sequence_number) seen.add(sequence_number) expected = fixture.expected_packets.get(sequence_number) if expected is None or packet.packet != expected: unexpected.append(sequence_number) continue exact_sequences.add(sequence_number) tx_preamble = fixture.starts[sequence_number] + GUARD_SAMPLES + tx_filter_delay ideal_rx_start = tx_preamble * scale + rx_filter_delay initial_error = (packet.start_sample - ideal_rx_start) / samples_per_symbol symbol_intervals = fixture.symbol_counts[sequence_number] - 1 final_error = initial_error - symbol_intervals * sample_clock_offset_ppm * 1e-6 initial_errors.append(initial_error) final_errors.append(final_error) drift_values.append(final_error - initial_error) cfo_errors.append(packet.carrier_offset_hz - effective_cfo_hz) packet_count = len(fixture.expected_packets) modeled_drift = -(next(iter(fixture.symbol_counts.values())) - 1) * sample_clock_offset_ppm * 1e-6 # Приёмник выдаёт только CRC-valid пакеты. Отдельной стадии "пакет до CRC" # в его API нет; число корреляционных кандидатов указано отдельно. return { "sample_clock_offset_ppm": sample_clock_offset_ppm, "cfo_hz": cfo_hz, "effective_cfo_hz_at_nominal_rx_rate": effective_cfo_hz, "initial_phase_degrees": initial_phase_degrees, "ebn0_db": ebn0_db, "payload_bytes": fixture.payload_bytes, "tx_packets": packet_count, "correlation_candidates": result.candidate_count, "detected_packets": len(result.packets), "crc_valid_packets": len(result.packets), "exact_packets": len(exact_sequences), "packet_success_rate": len(exact_sequences) / packet_count, "exact_sequences": sorted(exact_sequences), "missing_sequences": sorted(fixture.expected_packets.keys() - exact_sequences), "duplicate_sequences": duplicates, "unexpected_or_wrong_sequences": unexpected, "initial_timing_error_symbols_median": _median(initial_errors), "final_timing_error_symbols_median": _median(final_errors), "timing_drift_across_packet_symbols_median": _median(drift_values), "injected_timing_drift_across_packet_symbols": modeled_drift, "maximum_absolute_final_timing_error_symbols": ( max(map(abs, final_errors)) if final_errors else None ), "cfo_estimate_error_hz_median": _median(cfo_errors), "timing_error_method": "accepted-packet start versus injected time base; final error extrapolated with fixed 16-sample step", "reacquisition_per_packet": True, "continuous_timing_tracking_used": False, "receiver_time_s": receiver_time_s, "received_duration_s": received_samples / (SAMPLE_RATE_HZ * scale), "receiver_realtime_factor": receiver_time_s / (received_samples / (SAMPLE_RATE_HZ * scale)), } def evaluate_case(fixture: ClockFixture, *, ppm: float, cfo_hz: float = 0.0, phase_degrees: float = 0.0, ebn0_db: float | None = None, seed: int = SEED) -> dict: rx_iq = sample_independent_clock( fixture.iq, ppm, cfo_hz=cfo_hz, initial_phase_degrees=phase_degrees, ebn0_db=ebn0_db, seed=seed, ) started = time.perf_counter() result = receive_iq_stream(rx_iq) elapsed = time.perf_counter() - started return summarize_case(result, fixture, sample_clock_offset_ppm=ppm, cfo_hz=cfo_hz, initial_phase_degrees=phase_degrees, ebn0_db=ebn0_db, receiver_time_s=elapsed, received_samples=len(rx_iq)) def run() -> dict: """Снимает отдельную ppm-сетку, короткие комбинации и прогретую скорость.""" fixtures = {size: build_fixture(size) for size in PAYLOAD_SIZES} grid = [evaluate_case(fixtures[size], ppm=ppm) for size in PAYLOAD_SIZES for ppm in OFFSETS_PPM] combinations = [ (5, 20, 15, 9), (5, 20, 15, 12), (10, 100, 90, 9), (10, 100, 90, 12), (20, 100, -45, 9), (20, 100, -45, 12), ] combined = [evaluate_case(fixtures[256], ppm=ppm, cfo_hz=cfo, phase_degrees=phase, ebn0_db=ebn0, seed=SEED + index) for index, (ppm, cfo, phase, ebn0) in enumerate(combinations)] single_long_fixture = build_fixture(1012, packet_count=1) single_long_stress = [evaluate_case(single_long_fixture, ppm=ppm) for ppm in (-100, 100)] from tools.video005_benchmark import build_radio_fixture benchmark_packets, clean_benchmark_iq, _ = build_radio_fixture() benchmark_iq = sample_independent_clock(clean_benchmark_iq, 20, cfo_hz=20, initial_phase_degrees=15, ebn0_db=9) benchmark_times = [] benchmark_exact = [] for _ in range(3): started = time.perf_counter() received = receive_iq_stream(benchmark_iq) benchmark_times.append(time.perf_counter() - started) observed = {parse_packet(packet.packet).sequence_number: packet.packet for packet in received.packets} benchmark_exact.append(sum(observed.get(index) == packet for index, packet in enumerate(benchmark_packets))) duration = len(benchmark_iq) / (SAMPLE_RATE_HZ * (1 + 20e-6)) benchmark_all_exact = all(count == len(benchmark_packets) for count in benchmark_exact) primary = [case for case in grid if abs(case["sample_clock_offset_ppm"]) <= PRIMARY_LIMIT_PPM] primary_all_exact = all(case["exact_packets"] == case["tx_packets"] for case in primary) stress = [case for case in grid if abs(case["sample_clock_offset_ppm"]) > PRIMARY_LIMIT_PPM] stress_all_exact = all(case["exact_packets"] == case["tx_packets"] for case in stress) all_combined_exact = all(case["exact_packets"] == case["tx_packets"] for case in combined) symmetric_limits = {} for size in PAYLOAD_SIZES: maximum = 0 for offset in (1, 5, 10, 20, 50, 100): pair = [case for case in grid if case["payload_bytes"] == size and abs(case["sample_clock_offset_ppm"]) == offset] if len(pair) != 2 or any(case["exact_packets"] != case["tx_packets"] for case in pair): break maximum = offset symmetric_limits[str(size)] = maximum return { "status": ("PASS" if primary_all_exact and stress_all_exact and all_combined_exact and benchmark_all_exact else "PARTIAL"), "experiment": "CLOCK-001 INDEPENDENT TX/RX CLOCK READINESS", "seed": SEED, "sample_rate_nominal_hz": SAMPLE_RATE_HZ, "symbol_rate": SYMBOL_RATE, "fec_rate": "2/3", "packets_per_case": PACKETS_PER_CASE, "payload_sizes": PAYLOAD_SIZES, "offsets_ppm": OFFSETS_PPM, "primary_limit_ppm": PRIMARY_LIMIT_PPM, "grid": grid, "combined": combined, "single_long_stress": single_long_stress, "primary_grid_all_exact": primary_all_exact, "stress_grid_all_exact": stress_all_exact, "combined_all_exact": all_combined_exact, "benchmark_all_exact": benchmark_all_exact, "maximum_symmetric_ppm_with_all_exact_by_payload": symmetric_limits, "maximum_symmetric_ppm_with_all_payloads_exact": min(symmetric_limits.values()), "per_packet_reacquisition_sufficient_for_primary_grid": primary_all_exact, "continuous_tracking_required_for_primary_grid": not primary_all_exact, "intra_packet_timing_correction_needed_to_cover_long_packet_stress": not stress_all_exact, "benchmark": { "description": "VIDEO-005 32 long packets; independent clocks +20 ppm, CFO 20 Hz, phase 15 deg, Eb/N0 9 dB", "iq_duration_s": duration, "warm_receiver_times_s": benchmark_times, "warm_realtime_factor": statistics.median(benchmark_times) / duration, "exact_packets_each_run": benchmark_exact, "tx_packets": len(benchmark_packets), }, "rf_tx": 0, "sdr_hardware": "NOT_ACCESSED", "hardware_status": "NOT_RUN", } def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--numba-packages", type=Path) parser.add_argument("--report", type=Path, default=OUTPUT) args = parser.parse_args() if args.numba_packages is not None: if not args.numba_packages.is_dir(): parser.error("каталог Numba недоступен") sys.path.insert(0, str(args.numba_packages.resolve())) if args.report.exists(): parser.error("отчёт уже существует; выберите новый путь") report = run() args.report.parent.mkdir(parents=True, exist_ok=True) args.report.write_text(json.dumps(report, ensure_ascii=False, indent=2) + "\n", encoding="utf-8") print(f"CLOCK-001: {report['status']}; primary_all_exact={report['primary_grid_all_exact']}; " f"combined_all_exact={report['combined_all_exact']}; " f"realtime_factor={report['benchmark']['warm_realtime_factor']:.3f}; " "RF TX=0") if __name__ == "__main__": main()