355 lines
17 KiB
Python
355 lines
17 KiB
Python
"""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()
|