774 lines
44 KiB
Python
774 lines
44 KiB
Python
"""Lab057: изолированная диагностика сохранённых символов, без аппаратных команд."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import csv
|
||
from collections import Counter
|
||
from dataclasses import asdict, dataclass
|
||
from functools import lru_cache
|
||
import hashlib
|
||
import json
|
||
from pathlib import Path
|
||
import sys
|
||
from types import SimpleNamespace
|
||
|
||
import numpy as np
|
||
from scipy.signal import oaconvolve
|
||
|
||
if __package__ in (None, ""):
|
||
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
||
|
||
from experiments import lab053_packet_processing as packet
|
||
from experiments import lab053_streaming_receiver as streaming
|
||
from experiments import lab054_hardware_control as historical
|
||
from experiments import lab055_energy_audit as energy
|
||
|
||
physical = packet.physical
|
||
radio = packet.previous.iq.radio
|
||
ROOT = Path(__file__).resolve().parents[1]
|
||
OUTPUT = ROOT / "data/processed/lab057"
|
||
TRAIN = np.arange(32, 608)
|
||
VALIDATION = np.arange(672, 1248)
|
||
VARIANTS = ("S0", "D0", "D1", "D2", "D3")
|
||
FS = 8_000_000
|
||
PROBE_SYMBOLS = 1424
|
||
PREAMBLE = physical.lab044.preamble_bits()
|
||
EXPECTED = radio.bpsk_modulate(PREAMBLE)
|
||
|
||
|
||
def plain(value):
|
||
"""JSON без потери комплексных статистик."""
|
||
if isinstance(value, dict):
|
||
return {str(key): plain(item) for key, item in value.items()}
|
||
if isinstance(value, (list, tuple)):
|
||
return [plain(item) for item in value]
|
||
if isinstance(value, np.ndarray):
|
||
return plain(value.tolist())
|
||
if isinstance(value, complex):
|
||
return {"real": value.real, "imag": value.imag}
|
||
if isinstance(value, np.generic):
|
||
return plain(value.item())
|
||
return value
|
||
|
||
|
||
def save(name: str, value: dict | list) -> None:
|
||
"""Писать исключительно новый результат Lab057."""
|
||
OUTPUT.mkdir(parents=True, exist_ok=True)
|
||
(OUTPUT / name).write_text(json.dumps(plain(value), ensure_ascii=False, indent=2), encoding="utf-8")
|
||
|
||
|
||
class ObservedProcessor(packet.PacketProcessor):
|
||
"""Штатный обработчик с наблюдением первого probe; вычисления не подменяются."""
|
||
|
||
def __init__(self):
|
||
super().__init__(viterbi=physical.conv.viterbi_decode_soft)
|
||
self.probes = []
|
||
|
||
def extract(self, samples, frame_start, symbol_count):
|
||
result = super().extract(samples, frame_start, symbol_count)
|
||
self.probes.append(result)
|
||
return result
|
||
|
||
def decode(self, samples, candidate, number=1):
|
||
self.probes = []
|
||
outcome, row = super().decode(samples, candidate, number)
|
||
probe = self.probes[0] if self.probes else None
|
||
row["probe"] = None if probe is None else {
|
||
"timing": probe.timing_offset_samples, "carrier": probe.carrier,
|
||
"soft": probe.normalized_real.tolist(), "passed": probe.passed,
|
||
}
|
||
row["fec_started"] = outcome.fec_attempted
|
||
row["crc_result"] = outcome.crc
|
||
row["header_accepted"] = outcome.phy_header_pass
|
||
row["extract_calls"] = len(self.probes)
|
||
return outcome, row
|
||
|
||
|
||
def reproduce_baseline() -> dict:
|
||
"""Полный штатный поток получает только RX, без начала, числа кадров и TX."""
|
||
rows = []
|
||
|
||
def observe(row, outcome):
|
||
rows.append(dict(row, candidate_id=len(rows) + 1,
|
||
reason=outcome.reason, accepted=outcome.accepted))
|
||
|
||
receiver = streaming.StreamingReceiver(processor=ObservedProcessor(), packet_observer=observe)
|
||
received = []
|
||
with (ROOT / "data/raw/lab054/rx_stream.cs8").open("rb") as source:
|
||
while block := source.read(2 * 65536):
|
||
values = np.frombuffer(block, np.int8)
|
||
samples = (values[::2].astype(np.float32) + 1j * values[1::2]).astype(np.complex64)
|
||
received.extend(receiver.push_iq_block(samples))
|
||
received.extend(receiver.finish())
|
||
expected = energy.read_json(ROOT / "data/processed/lab054/independent_result.json")
|
||
comparisons = [(row["absolute_refined"], row["reason"]) for row in rows]
|
||
old = [(row["start"], row["reason"]) for row in expected["rows"]]
|
||
result = {"independent_stream": True, "tx_truth_used": False, "rf_tx": 0,
|
||
"packets": len(received), "candidates": len(rows),
|
||
"reasons": dict(Counter(row["reason"] for row in rows)),
|
||
"fec_started": sum(row["fec_started"] for row in rows),
|
||
"crc_started": sum(row["crc_result"] != "NOT_REACHED" for row in rows),
|
||
"historical_order_and_starts_equal": comparisons == old, "rows": rows}
|
||
result["passed"] = (comparisons == old and len(received) == 0 and len(rows) == 100
|
||
and result["fec_started"] == result["crc_started"] == 0
|
||
and set(result["reasons"]) == {"PHY_HEADER_MAGIC_MISMATCH"})
|
||
save("baseline.json", result)
|
||
if not result["passed"]:
|
||
raise RuntimeError("Baseline отличается от исторического: альтернативы запрещены")
|
||
return result
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class Estimate:
|
||
timing_samples: float
|
||
cfo_hz: float
|
||
gain: complex
|
||
train_coherence: float
|
||
frequency_boundary: bool = False
|
||
timing_boundary: bool = False
|
||
|
||
|
||
def interpolate(values: np.ndarray, positions: np.ndarray) -> np.ndarray:
|
||
"""Оконная sinc, ±8 исходных отсчётов; целые позиции сохраняются точно."""
|
||
positions = np.asarray(positions, np.float64)
|
||
centers = np.floor(positions).astype(np.int64)
|
||
offsets = np.arange(-7, 9)
|
||
indexes = centers[..., None] + offsets
|
||
if indexes.min() < 0 or indexes.max() >= len(values):
|
||
raise ValueError("Опора интерполяции вне окна")
|
||
distance = positions[..., None] - indexes
|
||
weights = np.sinc(distance) * np.where(np.abs(distance) < 8, .5 + .5 * np.cos(np.pi * distance / 8), 0)
|
||
weights /= weights.sum(axis=-1, keepdims=True)
|
||
return np.sum(values[indexes] * weights, axis=-1)
|
||
|
||
|
||
@lru_cache(maxsize=2)
|
||
def frequency_grid(symbol_rate: int) -> tuple[np.ndarray, np.ndarray]:
|
||
"""Один заранее зафиксированный диапазон в циклах/символ."""
|
||
frequencies = np.linspace(-250, 250, 501) * (symbol_rate / 500_000)
|
||
times = (TRAIN + .5) / symbol_rate
|
||
rotations = np.exp(-2j * np.pi * frequencies[:, None] * times)
|
||
return frequencies, rotations
|
||
|
||
|
||
def fit_gain(symbols: np.ndarray, symbol_rate: int, cfo_hz: float) -> complex:
|
||
"""Единственная оценка комплексного коэффициента по TRAIN, фаза при t=0."""
|
||
times = (TRAIN + .5) / symbol_rate
|
||
return complex(np.mean(symbols * EXPECTED[TRAIN] * np.exp(-2j * np.pi * cfo_hz * times)))
|
||
|
||
|
||
def coherence(symbols: np.ndarray, gain: complex) -> float:
|
||
"""Нормированный TRAIN-критерий, без проверочных и информационных битов."""
|
||
return float(abs(gain) ** 2 / max(float(np.mean(abs(symbols) ** 2)), 1e-30))
|
||
|
||
|
||
def frequency_search(symbols: np.ndarray, symbol_rate: int) -> tuple[float, complex, float, bool]:
|
||
"""Сетка и единственное локальное уточнение; диапазон не расширяется."""
|
||
frequencies, rotations = frequency_grid(symbol_rate)
|
||
power = abs(rotations @ (symbols * EXPECTED[TRAIN])) ** 2
|
||
index = int(np.argmax(power))
|
||
boundary = index in (0, len(frequencies) - 1)
|
||
offset = 0.0
|
||
if not boundary:
|
||
left, middle, right = power[index - 1:index + 2]
|
||
denominator = left - 2 * middle + right
|
||
if denominator < 0:
|
||
offset = float(np.clip(.5 * (left - right) / denominator, -.5, .5))
|
||
frequency = float(frequencies[index] + offset * (frequencies[1] - frequencies[0]))
|
||
gain = fit_gain(symbols, symbol_rate, frequency)
|
||
return frequency, gain, coherence(symbols, gain), boundary
|
||
|
||
|
||
def estimate_variant(matched: np.ndarray, nominal: int, symbol_rate: int,
|
||
variant: str, stock_timing: float, stock_cfo: float) -> Estimate:
|
||
"""Оцениватель не получает ни header/payload, ни VALIDATION/CRC."""
|
||
samples_per_symbol = FS // symbol_rate
|
||
if variant == "D3":
|
||
offsets = np.linspace(-.5, .5, 33) * samples_per_symbol
|
||
else:
|
||
offsets = np.array([stock_timing])
|
||
best = None
|
||
for offset in offsets:
|
||
positions = nominal + offset + TRAIN * samples_per_symbol
|
||
symbols = interpolate(matched, positions)
|
||
if variant in ("D2", "D3"):
|
||
frequency, gain, score, boundary = frequency_search(symbols, symbol_rate)
|
||
else:
|
||
frequency = stock_cfo if variant == "D0" else 0.0
|
||
gain = fit_gain(symbols, symbol_rate, frequency)
|
||
score, boundary = coherence(symbols, gain), False
|
||
candidate = Estimate(float(offset), frequency, gain, score, boundary,
|
||
variant == "D3" and abs(offset) == samples_per_symbol / 2)
|
||
if best is None or candidate.train_coherence > best.train_coherence:
|
||
best = candidate
|
||
return best
|
||
|
||
|
||
def compensate(symbols: np.ndarray, symbol_rate: int, estimate: Estimate) -> np.ndarray:
|
||
"""Один знак компенсации, единая шкала времени для всех вариантов."""
|
||
if abs(estimate.gain) == 0:
|
||
raise ValueError("Нулевой коэффициент")
|
||
times = (np.arange(len(symbols)) + .5) / symbol_rate
|
||
return symbols * np.exp(-2j * np.pi * estimate.cfo_hz * times) / estimate.gain
|
||
|
||
|
||
def squared_cfo_statistics(symbols: np.ndarray, symbol_rate: int) -> dict:
|
||
"""Разложить ровно активную штатную двухблочную оценку, без иной модели."""
|
||
blocks = (symbols[:1024] ** 2).reshape(2, 512).mean(axis=1)
|
||
phase = np.unwrap(np.angle(blocks))
|
||
delta = float(phase[1] - phase[0])
|
||
interval = 512 / symbol_rate
|
||
return {"symbol_ranges": [[0, 512], [512, 1024]], "block_values": blocks,
|
||
"block_abs": abs(blocks), "phase_difference_rad": delta,
|
||
"block_center_times_seconds": [256 / symbol_rate, 768 / symbol_rate],
|
||
"delta_time_seconds": interval, "cfo_hz": delta / (4 * np.pi * interval),
|
||
"unambiguous_abs_hz": symbol_rate / 2048,
|
||
"formula": "delta_phase / (4*pi*delta_time)",
|
||
"input_rate_hz": symbol_rate, "discarded_symbols": [1024, 1280]}
|
||
|
||
|
||
class Frontend:
|
||
"""Доступ к фактическому RX-фильтру и штатному извлечению нужной скорости."""
|
||
|
||
def __init__(self, symbol_rate: int):
|
||
self.symbol_rate = symbol_rate
|
||
self.samples_per_symbol = FS // symbol_rate
|
||
self.processor = packet.PacketProcessor(viterbi=physical.conv.viterbi_decode_soft)
|
||
self.taps = packet.bind(packet.lab043._matched_filter_taps,
|
||
SAMPLES_PER_SYMBOL=self.samples_per_symbol)()
|
||
self.nominal = 16 * self.samples_per_symbol + len(self.taps) - 1
|
||
|
||
def filter(self, values: np.ndarray) -> np.ndarray:
|
||
if self.symbol_rate == 500_000:
|
||
return self.processor.filter(self.processor.shift(values))
|
||
return packet.lab043._matched_waveform(values, 0, len(values))[0]
|
||
|
||
def extract(self, values: np.ndarray, count: int):
|
||
if self.symbol_rate == 500_000:
|
||
return self.processor.extract(values, 0, count)
|
||
return packet.previous.iq._extract_soft_window(values, 0, count)
|
||
|
||
def train_support(self) -> tuple[int, int]:
|
||
"""Сырая опора всех TRAIN-проб: FIR, sinc и полсимвола по краям."""
|
||
margin = int(np.ceil(self.samples_per_symbol / 2)) + 8
|
||
low = self.nominal + int(TRAIN[0]) * self.samples_per_symbol - len(self.taps) + 1 - margin
|
||
high = self.nominal + int(TRAIN[-1]) * self.samples_per_symbol + margin + 1
|
||
return low, high
|
||
|
||
def estimation_filter(self, values: np.ndarray) -> np.ndarray:
|
||
"""Те же RRC/carrier; оценивание не читает IQ вне сырой опоры TRAIN.
|
||
|
||
Локальный DC здесь — nuisance preprocessing оценивателя. Решения символов
|
||
используют исходный штатный filter(), поэтому D0 не меняет вход решения.
|
||
Короткая свёртка исключает даже округления FFT от проверочного участка.
|
||
"""
|
||
low, high = self.train_support()
|
||
local = values[low:high]
|
||
rotation = packet.lab043._carrier_downconversion_vector(low, len(local))
|
||
filtered = oaconvolve((local - np.mean(local)) * rotation, self.taps, mode="full")
|
||
matched = np.zeros(len(values) + len(self.taps) - 1, dtype=filtered.dtype)
|
||
matched[low:low + len(filtered)] = filtered
|
||
return matched
|
||
|
||
|
||
def decode_body(soft: np.ndarray) -> tuple[dict, bytes | None]:
|
||
"""Существующая проверка header/FEC/CRC; наблюдение реальных вызовов."""
|
||
calls = Counter()
|
||
|
||
def viterbi(values, count):
|
||
calls["fec"] += 1
|
||
return physical.conv.viterbi_decode_soft(values, count)
|
||
|
||
def parse(values):
|
||
calls["crc"] += 1
|
||
return physical.parse_packet(values)
|
||
|
||
decoder = packet.bind(physical.analyze_post_preamble_soft,
|
||
conv=SimpleNamespace(encoded_bit_count=physical.conv.encoded_bit_count,
|
||
viterbi_decode_soft=viterbi), parse_packet=parse)
|
||
reception = decoder(soft)
|
||
result = reception.analysis
|
||
return {"HEADER_ACCEPTED": bool(result["phy_header_pass"]),
|
||
"FEC_STARTED": bool(calls["fec"]),
|
||
"FEC_RESULT": ("PASS" if result["fec_decode_pass"] else "FAIL") if calls["fec"] else "NOT_REACHED",
|
||
"CRC_STARTED": bool(calls["crc"]), "CRC_RESULT": result["crc"],
|
||
"PACKET_OUTPUT": result["crc"] == physical.CRC_PASS,
|
||
"protocol_header_pass": bool(result["protocol_header_pass"]),
|
||
"decoder_analysis": result}, reception.packet_bytes
|
||
|
||
|
||
def error_metrics(normalized: np.ndarray, symbols: np.ndarray, estimate: Estimate,
|
||
symbol_rate: int, expected_header: np.ndarray | None) -> dict:
|
||
"""Оценка ошибки только после заморозки параметров; эталон не идёт оценивателю."""
|
||
decisions = normalized.real >= 0
|
||
errors = decisions[:1280] != PREAMBLE
|
||
header = normalized.real[1280:1424]
|
||
voted = header.reshape(48, 3).sum(axis=1) >= 0
|
||
times = (VALIDATION + .5) / symbol_rate
|
||
predicted = estimate.gain * EXPECTED[VALIDATION] * np.exp(2j * np.pi * estimate.cfo_hz * times)
|
||
residual = symbols[VALIDATION] - predicted
|
||
result = {"preamble_errors": int(errors.sum()), "train_errors": int(errors[TRAIN].sum()),
|
||
"validation_errors": int(errors[VALIDATION].sum()),
|
||
"train_symbols": len(TRAIN), "validation_symbols": len(VALIDATION),
|
||
"validation_residual_power": float(np.mean(abs(residual) ** 2)),
|
||
"validation_model_power": float(np.mean(abs(predicted) ** 2)),
|
||
"validation_actual_power": float(np.mean(abs(symbols[VALIDATION]) ** 2)),
|
||
"decoded_header_hex": np.packbits(voted).tobytes().hex(),
|
||
"header_raw_errors": None, "header_voted_errors": None, "magic_errors": None}
|
||
if expected_header is not None:
|
||
result.update(header_raw_errors=int(np.count_nonzero(decisions[1280:1424] != expected_header)),
|
||
header_voted_errors=int(np.count_nonzero(voted != expected_header[::3])),
|
||
magic_errors=int(np.count_nonzero(voted[:16] != expected_header[::3][:16])))
|
||
return result
|
||
|
||
|
||
def analyze_frame(reader, start: int, symbol_rate: int, group: str, frame_number: int,
|
||
expected_header: np.ndarray | None, stock_probe: dict | None = None,
|
||
*, write_packets: bool = True) -> tuple[list, dict, dict]:
|
||
"""Диагностика при известном начале; длина тела берётся только из header."""
|
||
frontend = Frontend(symbol_rate)
|
||
samples_per_symbol = frontend.samples_per_symbol
|
||
values = reader.read(start, (PROBE_SYMBOLS + 42) * samples_per_symbol)
|
||
window = frontend.extract(values, PROBE_SYMBOLS)
|
||
matched = frontend.filter(values)
|
||
estimation_matched = frontend.estimation_filter(values)
|
||
nominal = frontend.nominal
|
||
if stock_probe is not None:
|
||
if not np.array_equal(window.normalized_real, stock_probe["soft"]):
|
||
raise AssertionError("Диагностический probe не совпал со штатным baseline")
|
||
stock_symbols = matched[nominal + window.timing_offset_samples:
|
||
nominal + window.timing_offset_samples + PROBE_SYMBOLS * samples_per_symbol:
|
||
samples_per_symbol].astype(np.complex128)
|
||
statistics = squared_cfo_statistics(stock_symbols[:1280], symbol_rate)
|
||
if not np.isclose(statistics["cfo_hz"], window.carrier["cfo_hz"], atol=1e-8):
|
||
raise AssertionError("Разложение CFO не воспроизвело штатную ветвь")
|
||
results, observations = [], {}
|
||
for variant in VARIANTS:
|
||
if not window.passed and variant != "D3":
|
||
results.append({"group": group, "frame": frame_number, "variant": variant,
|
||
"start_sample": start, "known_start_diagnostic": True, "symbol_rate": symbol_rate,
|
||
"timing_samples": window.timing_offset_samples,
|
||
"timing_symbols": window.timing_offset_samples / samples_per_symbol,
|
||
"cfo_hz": window.carrier.get("cfo_hz"), "phase_rad": window.carrier.get("phase_rad"),
|
||
"phase_epoch": "t=0, symbol k at (k+0.5)/Rs", "gain_abs": window.carrier.get("abs_h"),
|
||
"train_coherence": None, "frequency_boundary": False, "timing_boundary": False,
|
||
"header_reason": window.reason, "accepted_packet_bytes": None,
|
||
"preamble_errors": None, "train_errors": None, "validation_errors": None,
|
||
"train_symbols": len(TRAIN), "validation_symbols": len(VALIDATION),
|
||
"header_raw_errors": None, "header_voted_errors": None, "magic_errors": None,
|
||
"HEADER_ACCEPTED": False, "FEC_STARTED": False, "FEC_RESULT": "NOT_REACHED",
|
||
"CRC_STARTED": False, "CRC_RESULT": "NOT_REACHED", "PACKET_OUTPUT": False,
|
||
"diagnostic_available": False})
|
||
continue
|
||
if variant == "S0":
|
||
gain = window.carrier["abs_h"] * np.exp(1j * window.carrier["phase_rad"])
|
||
estimate = Estimate(window.timing_offset_samples, window.carrier["cfo_hz"], gain,
|
||
coherence(stock_symbols[TRAIN], gain))
|
||
symbols = stock_symbols
|
||
else:
|
||
estimate = estimate_variant(estimation_matched, nominal, symbol_rate, variant,
|
||
window.timing_offset_samples, window.carrier["cfo_hz"])
|
||
symbols = interpolate(matched, nominal + estimate.timing_samples +
|
||
np.arange(PROBE_SYMBOLS) * samples_per_symbol)
|
||
normalized = compensate(symbols, symbol_rate, estimate)
|
||
if variant == "S0" and not np.allclose(normalized.real, window.normalized_real, atol=1e-12, rtol=1e-12):
|
||
raise AssertionError("S0 отличается от штатной нормировки")
|
||
header = physical.decode_phy_header_soft(normalized.real[1280:])
|
||
result = {"group": group, "frame": frame_number, "variant": variant, "start_sample": start,
|
||
"diagnostic_available": True,
|
||
"known_start_diagnostic": True, "symbol_rate": symbol_rate,
|
||
"timing_samples": estimate.timing_samples,
|
||
"timing_symbols": estimate.timing_samples / samples_per_symbol,
|
||
"cfo_hz": estimate.cfo_hz, "phase_rad": float(np.angle(estimate.gain)),
|
||
"phase_epoch": "t=0, symbol k at (k+0.5)/Rs",
|
||
"estimator_preprocessing": "STOCK" if variant == "S0" else "TRAIN_SUPPORT_DC",
|
||
"decision_preprocessing": "STOCK_GLOBAL_DC",
|
||
"gain_abs": abs(estimate.gain), "train_coherence": estimate.train_coherence,
|
||
"frequency_boundary": estimate.frequency_boundary, "timing_boundary": estimate.timing_boundary,
|
||
"header_reason": header["reason"], "accepted_packet_bytes": header["packet_bytes"],
|
||
**error_metrics(normalized, symbols, estimate, symbol_rate, expected_header)}
|
||
body_soft = normalized.real[1280:]
|
||
if header["passed"]:
|
||
packet_bytes = int(header["packet_bytes"])
|
||
count = 1280 + 144 + physical.conv.encoded_bit_count(packet_bytes * 8)
|
||
length = (count + 42) * samples_per_symbol
|
||
if start + length > reader.size:
|
||
result.update(HEADER_ACCEPTED=True, FEC_STARTED=False, FEC_RESULT="NOT_REACHED",
|
||
CRC_STARTED=False, CRC_RESULT="NOT_REACHED", PACKET_OUTPUT=False,
|
||
body_reason="TRUNCATED_FULL_IQ_WINDOW")
|
||
results.append(result)
|
||
continue
|
||
full_values = reader.read(start, length)
|
||
if variant == "S0":
|
||
full_window = frontend.extract(full_values, count)
|
||
body_soft = full_window.normalized_real[1280:]
|
||
else:
|
||
full_matched = frontend.filter(full_values)
|
||
full_symbols = interpolate(full_matched, nominal + estimate.timing_samples +
|
||
np.arange(count) * samples_per_symbol)
|
||
body_soft = compensate(full_symbols, symbol_rate, estimate).real[1280:]
|
||
decoded, packet_bytes = decode_body(body_soft)
|
||
result.update(decoded)
|
||
if result["PACKET_OUTPUT"] and write_packets:
|
||
name = f"packet_{group}_{frame_number:02d}_{variant}.bin"
|
||
(OUTPUT / name).write_bytes(packet_bytes)
|
||
result["packet_file"] = name
|
||
result["packet_sha256"] = hashlib.sha256(packet_bytes).hexdigest()
|
||
observations[variant] = {"normalized": normalized, "symbols": symbols,
|
||
"validation_residual": (symbols[VALIDATION] - estimate.gain * EXPECTED[VALIDATION] *
|
||
np.exp(2j * np.pi * estimate.cfo_hz * (VALIDATION + .5) / symbol_rate))}
|
||
results.append(result)
|
||
return results, observations, statistics
|
||
|
||
|
||
def candidate_mapping(baseline: dict, frames: list) -> list:
|
||
"""Связать независимые кандидаты с TX только после завершения baseline."""
|
||
mapping = []
|
||
chosen = {min(range(len(baseline["rows"])), key=lambda i: abs(baseline["rows"][i]["absolute_refined"] - frame["start"])):
|
||
frame["frame"] for frame in frames}
|
||
for index, row in enumerate(baseline["rows"]):
|
||
nearest = min(frames, key=lambda frame: abs(frame["start"] - row["absolute_refined"]))
|
||
delta = row["absolute_refined"] - nearest["start"]
|
||
mapping.append({"candidate_id": row["candidate_id"], "coarse_start": row["base_sample"] + row["coarse"],
|
||
"refined_start": row["absolute_refined"], "nearest_frame": nearest["frame"],
|
||
"difference_from_diagnostic_start": delta,
|
||
"true_start_frame": chosen.get(index), "inside_known_tx": abs(delta) < 92256,
|
||
"reason": row["reason"]})
|
||
if len(chosen) != 16 or any(abs(row["difference_from_diagnostic_start"]) > 1 for row in mapping if row["true_start_frame"]):
|
||
raise AssertionError("Не воспроизведены 16 настоящих начал")
|
||
return mapping
|
||
|
||
|
||
def summarize(rows: list) -> dict:
|
||
"""Все варианты раздельно; лучшего по CRC не выбираем."""
|
||
output = {}
|
||
for group in sorted({row["group"] for row in rows}):
|
||
output[group] = {}
|
||
for variant in VARIANTS:
|
||
selected = [row for row in rows if row["group"] == group and row["variant"] == variant]
|
||
output[group][variant] = {"frames": len(selected),
|
||
"header_accepted": sum(row["HEADER_ACCEPTED"] for row in selected),
|
||
"fec_started": sum(row["FEC_STARTED"] for row in selected),
|
||
"crc_pass": sum(row["CRC_RESULT"] == physical.CRC_PASS for row in selected),
|
||
"diagnostic_available": sum(row.get("diagnostic_available", True) for row in selected),
|
||
"train_errors": sum(row["train_errors"] or 0 for row in selected),
|
||
"validation_errors": sum(row["validation_errors"] or 0 for row in selected),
|
||
"header_voted_errors": (sum(row["header_voted_errors"] or 0 for row in selected)
|
||
if any(row["header_voted_errors"] is not None for row in selected) else None),
|
||
"frequency_boundaries": sum(row["frequency_boundary"] for row in selected),
|
||
"timing_boundaries": sum(row["timing_boundary"] for row in selected)}
|
||
return output
|
||
|
||
|
||
def write_tables(rows: list) -> None:
|
||
"""Машинная полная таблица и плоский CSV с состояниями строгого декодера."""
|
||
save("frame_results.json", rows)
|
||
fields = [key for key in rows[0] if key != "decoder_analysis"]
|
||
fields += ["packet_file", "packet_sha256", "candidate_id"]
|
||
fields = list(dict.fromkeys(fields))
|
||
with (OUTPUT / "frame_results.csv").open("w", newline="", encoding="utf-8-sig") as stream:
|
||
writer = csv.DictWriter(stream, fields, extrasaction="ignore")
|
||
writer.writeheader()
|
||
writer.writerows(rows)
|
||
|
||
|
||
def run_diagnosis(baseline: dict) -> dict:
|
||
"""Фиксированный набор сохранённых кадров и контролей; без аппаратуры."""
|
||
audit = energy.read_json(ROOT / "data/processed/lab055/energy_audit_v2.json")
|
||
mapping = candidate_mapping(baseline, audit["experiments"]["Lab054"]["rx"]["frames"])
|
||
save("candidate_mapping.json", mapping)
|
||
save("other_candidates.json", [row for row, mapped in zip(baseline["rows"], mapping) if not mapped["true_start_frame"]])
|
||
rx = energy.Cs8File(ROOT / "data/raw/lab054/rx_stream.cs8")
|
||
plans = historical.plans()
|
||
rows, diagnostics, observations = [], [], {}
|
||
for mapped in mapping:
|
||
if not mapped["true_start_frame"]:
|
||
continue
|
||
number = mapped["true_start_frame"]
|
||
selected, observed, statistics = analyze_frame(rx, mapped["refined_start"], 500_000, "Lab054", number,
|
||
plans[number - 1].plan.phy_header_bits, baseline["rows"][mapped["candidate_id"] - 1]["probe"])
|
||
for row in selected:
|
||
row["candidate_id"] = mapped["candidate_id"]
|
||
rows.extend(selected)
|
||
diagnostics.append({"group": "Lab054", "frame": number, **statistics})
|
||
for variant, values in observed.items():
|
||
for kind, array in values.items():
|
||
observations[f"frame{number:02d}_{variant}_{kind}"] = array
|
||
print(f"Lab054 frame {number}/16", flush=True)
|
||
tx = energy.Cs8File(ROOT / "data/raw/lab054/tx_stream.cs8")
|
||
for number, plan in enumerate(plans, 1):
|
||
selected, _, _ = analyze_frame(tx, (number - 1) * 92256, 500_000, "TX054", number, plan.plan.phy_header_bits)
|
||
rows.extend(selected)
|
||
old_directory = Path(audit["experiments"]["Lab047"]["directory"])
|
||
old = energy.read_json(old_directory / "independent_reprocess_result_corrected.json")
|
||
old_rx = energy.Cs8File(old_directory / "rx_stream.cs8")
|
||
accepted = [outcome for outcome in old["strict"]["outcomes"] if outcome["accepted"]]
|
||
for number, outcome in enumerate(accepted, 1):
|
||
start = outcome["refined_start_sample"]
|
||
# Все исторические кадры имеют 268 байт; эталон используется только в error_metrics.
|
||
expected_header = physical.encode_phy_header(268)
|
||
selected, _, statistics = analyze_frame(old_rx, start, 20_000, "RX047", number, expected_header)
|
||
rows.extend(selected)
|
||
diagnostics.append({"group": "RX047", "frame": number, **statistics})
|
||
print(f"Lab047 frame {number}/36", flush=True)
|
||
for number in range(16):
|
||
selected, _, _ = analyze_frame(rx, 8_000_000 + number * 100_000, 500_000,
|
||
"BACKGROUND054", number + 1, None)
|
||
rows.extend(selected)
|
||
write_tables(rows)
|
||
save("squared_cfo_statistics.json", diagnostics)
|
||
np.savez_compressed(OUTPUT / "symbol_observations.npz", **observations)
|
||
summary = {"rf_tx": 0, "historical_lab054": "NOT_A, 0/16", "baseline_reproduced": baseline["passed"],
|
||
"train": [32, 608], "validation": [672, 1248], "known_start_diagnostics": True,
|
||
"groups": summarize(rows)}
|
||
save("summary.json", summary)
|
||
return summary
|
||
|
||
|
||
class ArrayFile:
|
||
"""Синтетический контроль с тем же интерфейсом, без файла и оборудования."""
|
||
|
||
def __init__(self, values: np.ndarray):
|
||
self.values = values
|
||
self.size = len(values)
|
||
|
||
def read(self, start: int, count: int) -> np.ndarray:
|
||
if start < 0 or start + count > self.size:
|
||
raise ValueError("Синтетическое окно вне массива")
|
||
return self.values[start:start + count]
|
||
|
||
|
||
def synthetic_waveform(symbol_rate: int, *, delay_symbols: float = .171875,
|
||
cfo_hz: float = 0, phase_rad: float = .7,
|
||
es_n0_db: float | None = None, seed: int = 5701) -> tuple[np.ndarray, dict]:
|
||
"""Физические параметры задаёт генератор; оцениватель их не получает.
|
||
|
||
При сравнении скоростей равен Es/N0 (по средней цифровой энергии),
|
||
дисперсия комплексного белого шума масштабируется с SPS.
|
||
"""
|
||
mode = historical.phy.PhyMode(symbol_rate)
|
||
waveform = historical.phy.build_parameterized_tx_complex(historical.plans()[0].plan, mode)
|
||
indexes = np.arange(len(waveform))
|
||
baseband = waveform * np.exp(-2j * np.pi * 50_000 * indexes / FS)
|
||
delay_samples = delay_symbols * mode.samples_per_symbol
|
||
padding = int(np.ceil(abs(delay_samples))) + 20
|
||
delayed = interpolate(np.pad(baseband, padding), padding + indexes - delay_samples)
|
||
values = delayed * np.exp(1j * (phase_rad + 2 * np.pi * (50_000 + cfo_hz) * indexes / FS))
|
||
signal_power = float(np.mean(abs(values) ** 2))
|
||
variance = 0.0
|
||
if es_n0_db is not None:
|
||
variance = signal_power * mode.samples_per_symbol / 10 ** (es_n0_db / 10)
|
||
random_generator = np.random.default_rng(seed)
|
||
noise = (random_generator.normal(size=len(values)) + 1j * random_generator.normal(size=len(values)))
|
||
values = values + np.sqrt(variance / 2) * noise
|
||
return values.astype(np.complex64), {"symbol_rate": symbol_rate, "delay_samples": delay_samples,
|
||
"delay_symbols": delay_symbols, "cfo_hz": cfo_hz, "phase_rad": phase_rad,
|
||
"phase_epoch": "raw n=0, before RX filtering", "es_n0_db": es_n0_db,
|
||
"complex_noise_variance": variance, "signal_power": signal_power, "seed": seed,
|
||
"comparison_condition": "equal digital Es/N0; noise variance scales with SPS"}
|
||
|
||
|
||
def synthetic_controls() -> dict:
|
||
"""Небольшие фиксированные контроли фазы, CFO, дробной задержки и шума."""
|
||
symbol_checks, waveform_checks = [], []
|
||
for symbol_rate in (20_000, 500_000):
|
||
times = (np.arange(1280) + .5) / symbol_rate
|
||
true_cfo = 37.5 * symbol_rate / 500_000
|
||
true_phase = .7
|
||
symbols = EXPECTED * np.exp(1j * (true_phase + 2 * np.pi * true_cfo * times))
|
||
frequency, gain, score, boundary = frequency_search(symbols[TRAIN], symbol_rate)
|
||
stock = squared_cfo_statistics(symbols, symbol_rate)
|
||
symbol_checks.append({"symbol_rate": symbol_rate, "true_cfo": true_cfo, "estimated_cfo": frequency,
|
||
"true_phase": true_phase, "estimated_phase": float(np.angle(gain)),
|
||
"stock_cfo": stock["cfo_hz"], "unambiguous_abs_hz": stock["unambiguous_abs_hz"],
|
||
"passed": abs(frequency - true_cfo) < .01 and abs(np.angle(gain) - true_phase) < .001})
|
||
outside = symbol_rate / 2048 * 1.2
|
||
aliased = squared_cfo_statistics(EXPECTED * np.exp(2j * np.pi * outside * times), symbol_rate)
|
||
symbol_checks.append({"symbol_rate": symbol_rate, "true_cfo": outside,
|
||
"estimated_cfo": aliased["cfo_hz"], "outside_stock_unambiguous_range": True,
|
||
"expected_alias_hz": outside - symbol_rate / 1024,
|
||
"passed": np.isclose(aliased["cfo_hz"], outside - symbol_rate / 1024)})
|
||
for case, noise_db, residual_cfo in (("fractional_clean", None, 0.0), ("fractional_noise", 12.0, true_cfo)):
|
||
values, truth = synthetic_waveform(symbol_rate, es_n0_db=noise_db, cfo_hz=residual_cfo, seed=5701)
|
||
frontend = Frontend(symbol_rate)
|
||
probe = values[:(PROBE_SYMBOLS + 42) * frontend.samples_per_symbol]
|
||
stock_window = frontend.extract(probe, PROBE_SYMBOLS)
|
||
estimation_matched = frontend.estimation_filter(probe)
|
||
output_matched = frontend.filter(probe)
|
||
for variant in VARIANTS[1:]:
|
||
estimate = estimate_variant(estimation_matched, frontend.nominal, symbol_rate, variant,
|
||
stock_window.timing_offset_samples, stock_window.carrier["cfo_hz"])
|
||
symbols = interpolate(output_matched, frontend.nominal + estimate.timing_samples +
|
||
np.arange(PROBE_SYMBOLS) * frontend.samples_per_symbol)
|
||
normalized = compensate(symbols, symbol_rate, estimate)
|
||
waveform_checks.append({"case": case, "variant": variant, "truth": truth,
|
||
"estimate": asdict(estimate), "timing_error_symbols":
|
||
estimate.timing_samples / frontend.samples_per_symbol - truth["delay_symbols"],
|
||
"cfo_error_hz": estimate.cfo_hz - residual_cfo,
|
||
"train_errors": int(np.count_nonzero((normalized.real[TRAIN] >= 0) != PREAMBLE[TRAIN])),
|
||
"validation_errors": int(np.count_nonzero((normalized.real[VALIDATION] >= 0) != PREAMBLE[VALIDATION]))})
|
||
# Полный контроль утечки файла-эталона, без подмены исходных файлов проекта.
|
||
values, _ = synthetic_waveform(500_000, delay_symbols=0, cfo_hz=0, es_n0_db=None)
|
||
reader = ArrayFile(values)
|
||
expected = historical.plans()[0].plan.phy_header_bits
|
||
first, _, _ = analyze_frame(reader, 0, 500_000, "LEAKAGE_CONTROL", 1, expected, write_packets=False)
|
||
changed, _, _ = analyze_frame(reader, 0, 500_000, "LEAKAGE_CONTROL", 1, 1 - expected, write_packets=False)
|
||
parameter_keys = ("timing_samples", "cfo_hz", "phase_rad", "gain_abs", "HEADER_ACCEPTED", "CRC_RESULT")
|
||
reference_invariance = all(all(a[key] == b[key] for key in parameter_keys) for a, b in zip(first, changed))
|
||
raw = energy.Cs8File(ROOT / "data/raw/lab054/rx_stream.cs8")
|
||
values = raw.read(4161352, (PROBE_SYMBOLS + 42) * 16)
|
||
altered = values.copy()
|
||
frontend = Frontend(500_000)
|
||
low, high = frontend.train_support()
|
||
altered[:low] += 10 + 7j
|
||
altered[high:] += 10 + 7j
|
||
original = estimate_variant(frontend.estimation_filter(values), frontend.nominal, 500_000, "D3", 0, 0)
|
||
perturbation = estimate_variant(frontend.estimation_filter(altered), frontend.nominal, 500_000, "D3", 99, 99)
|
||
clean_delay_passed = all(abs(row["timing_error_symbols"]) <= 1/32 and row["validation_errors"] == 0
|
||
for row in waveform_checks if row["variant"] == "D3" and row["case"] == "fractional_clean")
|
||
result = {"symbol_checks": symbol_checks, "waveform_checks": waveform_checks,
|
||
"clean_fractional_delay_passed": clean_delay_passed,
|
||
"expected_bits_invariance": reference_invariance,
|
||
"raw_outside_train_invariance": original == perturbation,
|
||
"raw_train_support_500k": [low, high],
|
||
"passed": all(row["passed"] for row in symbol_checks) and clean_delay_passed and reference_invariance and original == perturbation}
|
||
save("synthetic_controls.json", result)
|
||
if not result["passed"]:
|
||
raise AssertionError("Контроли синхронизации/утечки не прошли")
|
||
return result
|
||
|
||
|
||
def make_plots() -> None:
|
||
"""Созвездия, фазы, расположение ошибок и остаток всех 16 кадров."""
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
|
||
rows = energy.read_json(OUTPUT / "frame_results.json")
|
||
rows = [row for row in rows if row["group"] == "Lab054"]
|
||
with np.load(OUTPUT / "symbol_observations.npz") as arrays:
|
||
colors = dict(zip(VARIANTS, ("#666666", "#a76b25", "#14835b", "#2465bc", "#9135a3")))
|
||
for variant in VARIANTS:
|
||
figure, axes = plt.subplots(4, 4, figsize=(14, 12), constrained_layout=True)
|
||
for number, axis in enumerate(axes.flat, 1):
|
||
symbols = arrays[f"frame{number:02d}_{variant}_normalized"][VALIDATION]
|
||
axis.scatter(symbols.real, symbols.imag, s=3, alpha=.4, color=colors[variant])
|
||
axis.scatter([-1, 1], [0, 0], marker="x", color="black")
|
||
axis.set(title=f"Кадр {number}", xlabel="I / |h|", ylabel="Q / |h|")
|
||
axis.grid(alpha=.2)
|
||
figure.suptitle(f"Lab054, {variant}: созвездие VALIDATION, решения после штатного RRC")
|
||
figure.savefig(OUTPUT / f"constellation_{variant}.png", dpi=130)
|
||
plt.close(figure)
|
||
figure, axes = plt.subplots(4, 4, figsize=(15, 12), constrained_layout=True)
|
||
for number, axis in enumerate(axes.flat, 1):
|
||
for variant in VARIANTS:
|
||
corrected = arrays[f"frame{number:02d}_{variant}_normalized"][:1280] * EXPECTED
|
||
# Круговое среднее фазы, без переоценки компенсации.
|
||
unit = corrected / np.maximum(abs(corrected), 1e-30)
|
||
phase = np.angle(unit.reshape(40, 32).mean(axis=1))
|
||
axis.plot((np.arange(40) * 32 + 16) / 500, phase, color=colors[variant], label=variant, linewidth=.8)
|
||
axis.axvspan(32/500, 608/500, color="green", alpha=.06)
|
||
axis.axvspan(672/500, 1248/500, color="blue", alpha=.06)
|
||
axis.set(title=f"Кадр {number}", xlabel="Время преамбулы, мс", ylabel="Ошибка фазы, рад", ylim=(-np.pi, np.pi))
|
||
axes.flat[0].legend(fontsize=8, ncol=3)
|
||
figure.suptitle("Фаза относительно PRBS7; круговое среднее каждых 32 символов")
|
||
figure.savefig(OUTPUT / "phase_errors.png", dpi=130)
|
||
plt.close(figure)
|
||
errors, residuals, labels = [], [], []
|
||
for number in range(1, 17):
|
||
for variant in VARIANTS:
|
||
symbols = arrays[f"frame{number:02d}_{variant}_normalized"][:1280]
|
||
errors.append((symbols.real >= 0) != PREAMBLE)
|
||
residuals.append(abs(arrays[f"frame{number:02d}_{variant}_validation_residual"]) ** 2)
|
||
labels.append(f"{number:02d}/{variant}")
|
||
figure, axes = plt.subplots(1, 2, figsize=(16, 15), constrained_layout=True)
|
||
axes[0].imshow(errors, aspect="auto", interpolation="nearest", cmap="Greys", vmin=0, vmax=1)
|
||
plot = axes[1].imshow(10*np.log10(np.maximum(residuals, 1e-12)), aspect="auto", interpolation="nearest", cmap="viridis")
|
||
for axis in axes:
|
||
axis.set_yticks(np.arange(80), labels, fontsize=6)
|
||
axis.set_ylabel("Кадр / вариант")
|
||
axes[0].set(title="Ошибочные решения преамбулы", xlabel="Номер символа преамбулы")
|
||
axes[1].set(title="Мощность остатка VALIDATION, дБ относительно 1 код²", xlabel="Индекс внутри VALIDATION")
|
||
figure.colorbar(plot, ax=axes[1], shrink=.5)
|
||
figure.savefig(OUTPUT / "errors_and_residual.png", dpi=130)
|
||
plt.close(figure)
|
||
figure, axes = plt.subplots(2, 2, figsize=(13, 8), constrained_layout=True)
|
||
for variant in VARIANTS:
|
||
selected = sorted((row for row in rows if row["variant"] == variant), key=lambda row: row["frame"])
|
||
for axis, key in zip(axes.flat, ("validation_errors", "header_voted_errors", "cfo_hz", "timing_symbols")):
|
||
axis.plot(range(1, 17), [row[key] for row in selected], marker=".", label=variant, color=colors[variant])
|
||
axis.set_xlabel("Номер кадра")
|
||
axis.grid(alpha=.2)
|
||
for axis, title in zip(axes.flat, ("Ошибки VALIDATION из 576", "Ошибки заголовка из 48", "Оценка CFO, Гц", "Timing, доля символа")):
|
||
axis.set_title(title)
|
||
axes.flat[0].legend(ncol=5)
|
||
figure.savefig(OUTPUT / "comparison.png", dpi=140)
|
||
plt.close(figure)
|
||
|
||
|
||
def verify_integrity() -> dict:
|
||
"""Проверить исходники и старые артефакты; журнал разрешает только дописывание."""
|
||
initial = energy.read_json(OUTPUT / "initial_state.json")
|
||
changed_sources, changed_data = [], []
|
||
for name, expected in initial["source_files"].items():
|
||
if energy.sha256(ROOT / name) != expected["sha256"]:
|
||
changed_sources.append(name)
|
||
for name, expected in initial["historical_files"].items():
|
||
path = Path(name)
|
||
if path.name in ("PROJECT_LOG.md", "LABS_INDEX.md"):
|
||
actual = hashlib.sha256(path.read_bytes()[:expected["bytes"]]).hexdigest()
|
||
else:
|
||
actual = energy.sha256(path)
|
||
if actual != expected["sha256"]:
|
||
changed_data.append(name)
|
||
loaded = {}
|
||
for module in list(sys.modules.values()):
|
||
filename = getattr(module, "__file__", None)
|
||
if filename:
|
||
path = Path(filename).resolve()
|
||
if path.is_relative_to(ROOT) and path.suffix == ".py":
|
||
loaded[path.relative_to(ROOT).as_posix()] = energy.sha256(path)
|
||
result = {"source_changes": changed_sources, "historical_changes": changed_data,
|
||
"loaded_processing_sources_sha256": loaded,
|
||
"loaded_versions": {name: {"version": getattr(sys.modules.get(name), "__version__", None),
|
||
"file": getattr(sys.modules.get(name), "__file__", None)}
|
||
for name in ("numpy", "scipy", "matplotlib", "numba")},
|
||
"passed": not changed_sources and not changed_data}
|
||
save("integrity.json", result)
|
||
if not result["passed"]:
|
||
raise AssertionError("Изменились исходные файлы вне Lab057")
|
||
return result
|
||
|
||
|
||
def run_functional_tests() -> None:
|
||
"""Быстрые свойства без IQ-файлов и оборудования."""
|
||
assert len(TRAIN) == len(VALIDATION) == 576
|
||
assert VALIDATION[0] - TRAIN[-1] > 20
|
||
assert set(TRAIN).isdisjoint(VALIDATION)
|
||
times = (np.arange(1280) + .5) / 500_000
|
||
symbols = EXPECTED * np.exp(1j * (.7 + 2 * np.pi * 31.25 * times))
|
||
frequency, gain, score, boundary = frequency_search(symbols[TRAIN], 500_000)
|
||
assert abs(frequency - 31.25) < .01 and score > .99999 and not boundary
|
||
assert abs(np.angle(gain) - .7) < .001
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description=__doc__)
|
||
parser.add_argument("--baseline-only", action="store_true")
|
||
arguments = parser.parse_args()
|
||
run_functional_tests()
|
||
result = reproduce_baseline()
|
||
print("BASELINE", result["passed"], result["reasons"], "RF TX=0", flush=True)
|
||
if not arguments.baseline_only:
|
||
summary = run_diagnosis(result)
|
||
summary["controls_passed"] = synthetic_controls()["passed"]
|
||
make_plots()
|
||
summary["integrity_passed"] = verify_integrity()["passed"]
|
||
save("summary.json", summary)
|
||
print(json.dumps(plain(summary), ensure_ascii=False, indent=2))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|