turov_av/SDR-Rover
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

CLOCK-001: validate independent-clock QPSK receiver

Browse Source
This commit is contained in:
LittleSam129
2026-09-17 15:05:44 +03:00
parent 6df5f8aa39
commit 3c31bac0da
12 changed files with 4458 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

354
tools/clock001.py Normal file
View File

@@ -0,0 +1,354 @@
"""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()