550 lines
30 KiB
Python
550 lines
30 KiB
Python
"""Lab058: офлайн-анализ качества BPSK-символов записи Lab056."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import csv
|
||
import hashlib
|
||
import json
|
||
from pathlib import Path
|
||
import sys
|
||
|
||
if __package__ in (None, ""):
|
||
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
||
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
from scipy import signal
|
||
|
||
from experiments import lab043_pluto_to_rtlsdr as shaping
|
||
from experiments import lab057_symbol_diagnosis as lab057
|
||
|
||
|
||
PROJECT_ROOT = Path(__file__).resolve().parents[1]
|
||
RAW = PROJECT_ROOT / "data/raw/lab056"
|
||
LAB056 = PROJECT_ROOT / "data/processed/lab056"
|
||
OUTPUT = PROJECT_ROOT / "data/processed/lab058"
|
||
TX_PATH = RAW / "tx_abab.cs8"
|
||
RX_PATH = RAW / "rx_stream.cs8"
|
||
TX_SHA256 = "5b34bbe1a8e6d037d2d6d367a74e9eb84e8940db7da61600ffaa7f7497d4f600"
|
||
RX_SHA256 = "e336dd653fa554058052fcf851e9f2462f66b8cf7c710815ca2b9a09ffc23ad4"
|
||
SAMPLE_RATE = 8_000_000
|
||
CARRIER_HZ = 50_000
|
||
TILE_COUNT = 64
|
||
TIMING_GRID = np.linspace(-0.5, 0.5, 33)
|
||
ISI_OFFSETS = np.arange(-5, 6)
|
||
AWGN_SEED = 5801
|
||
AWGN_COMPLEX_VARIANCE = 400.0
|
||
EXPECTED_ENERGY_DIFFERENCE_DB = -13.870505
|
||
|
||
|
||
def sha256(path: Path) -> str:
|
||
digest = hashlib.sha256()
|
||
with path.open("rb") as stream:
|
||
for block in iter(lambda: stream.read(1024 * 1024), b""):
|
||
digest.update(block)
|
||
return digest.hexdigest()
|
||
|
||
|
||
class Cs8File:
|
||
"""Read-only interleaved signed-I/Q file."""
|
||
|
||
def __init__(self, path: Path):
|
||
self.path = path
|
||
self.raw = np.memmap(path, np.int8, mode="r")
|
||
if self.raw.size % 2:
|
||
raise ValueError(f"Нечётный размер cs8: {path}")
|
||
self.size = self.raw.size // 2
|
||
|
||
def read(self, start: int, count: int) -> np.ndarray:
|
||
if start < 0 or count < 0 or start + count > self.size:
|
||
raise ValueError(f"Окно вне {self.path}: {start}+{count}>{self.size}")
|
||
values = self.raw[2 * start:2 * (start + count)]
|
||
return values[::2].astype(np.float32) + 1j * values[1::2].astype(np.float32)
|
||
|
||
|
||
def expected_symbols(symbol_rate: int) -> np.ndarray:
|
||
count = 4096 * (symbol_rate // 20_000)
|
||
bits = shaping.pilot_prbs7(count)
|
||
if not np.array_equal(bits[:1280], lab057.PREAMBLE):
|
||
raise AssertionError("PRBS Lab056 не совпала со штатной преамбулой")
|
||
return 2.0 * bits.astype(np.float64) - 1.0
|
||
|
||
|
||
def tile_ids(count: int) -> np.ndarray:
|
||
return np.minimum(np.arange(count, dtype=np.int64) * TILE_COUNT // count, TILE_COUNT - 1)
|
||
|
||
|
||
def training_mask(count: int) -> np.ndarray:
|
||
return tile_ids(count) % 2 == 0
|
||
|
||
|
||
def guarded_tile_masks(count: int, margin_symbols: int = 8) -> tuple[np.ndarray, np.ndarray]:
|
||
"""Return disjoint even/odd tile interiors outside FIR/interpolation support."""
|
||
train = np.zeros(count, dtype=bool)
|
||
validation = np.zeros(count, dtype=bool)
|
||
for identifier in range(TILE_COUNT):
|
||
low = (identifier * count + TILE_COUNT - 1) // TILE_COUNT + margin_symbols
|
||
high = ((identifier + 1) * count + TILE_COUNT - 1) // TILE_COUNT - margin_symbols
|
||
if low >= high:
|
||
raise ValueError("Плитка слишком коротка для защитного края")
|
||
(train if identifier % 2 == 0 else validation)[low:high] = True
|
||
return train, validation
|
||
|
||
|
||
def oracle_matched_filter(values: np.ndarray, taps: np.ndarray,
|
||
symbol_rate: int, symbol_count: int) -> np.ndarray:
|
||
"""Matched filter whose DC nuisance estimate reads only guarded TRAIN intervals."""
|
||
samples_per_symbol = SAMPLE_RATE // symbol_rate
|
||
train, _ = guarded_tile_masks(symbol_count)
|
||
centers = (16 + 5) * samples_per_symbol + np.arange(symbol_count) * samples_per_symbol
|
||
raw_train = np.zeros(len(values), dtype=bool)
|
||
for index in np.flatnonzero(train):
|
||
raw_train[centers[index]:centers[index] + samples_per_symbol] = True
|
||
dc = np.mean(values[raw_train])
|
||
indexes = np.arange(len(values))
|
||
baseband = (values - dc) * np.exp(-2j * np.pi * CARRIER_HZ * indexes / SAMPLE_RATE)
|
||
return signal.oaconvolve(baseband, taps, mode="full")
|
||
|
||
|
||
def fit_cfo_gain(samples: np.ndarray, expected: np.ndarray, indexes: np.ndarray,
|
||
symbol_rate: int, total_symbol_count: int | None = None) -> tuple[float, complex, float]:
|
||
"""Fit one CFO and h using only the supplied indexes."""
|
||
corrected_sign = samples * expected
|
||
if total_symbol_count is None:
|
||
total_symbol_count = len(expected)
|
||
identifiers = np.minimum(indexes * TILE_COUNT // total_symbol_count, TILE_COUNT - 1)
|
||
block_values, block_times = [], []
|
||
for identifier in np.unique(identifiers):
|
||
select = identifiers == identifier
|
||
value = np.mean(corrected_sign[select])
|
||
if abs(value) > 0:
|
||
block_values.append(value)
|
||
block_times.append(float(np.mean((indexes[select] + 0.5) / symbol_rate)))
|
||
block_values = np.asarray(block_values)
|
||
block_times = np.asarray(block_times)
|
||
if len(block_times) == 0:
|
||
return 0.0, 0j, 0.0
|
||
phase = np.unwrap(np.angle(block_values))
|
||
weights = np.maximum(abs(block_values), 1e-12)
|
||
matrix = np.column_stack((np.ones(len(block_times)), block_times - block_times[0]))
|
||
coefficients = np.linalg.lstsq(matrix * weights[:, None], phase * weights, rcond=None)[0]
|
||
cfo_hz = float(coefficients[1] / (2 * np.pi))
|
||
times = (indexes + 0.5) / symbol_rate
|
||
h = complex(np.mean(corrected_sign * np.exp(-2j * np.pi * cfo_hz * times)))
|
||
coherence = float(abs(h) ** 2 / max(float(np.mean(abs(samples) ** 2)), 1e-30))
|
||
return cfo_hz, h, coherence
|
||
|
||
|
||
def oracle_estimate(matched: np.ndarray, nominal: int, expected: np.ndarray,
|
||
symbol_rate: int) -> dict:
|
||
samples_per_symbol = SAMPLE_RATE // symbol_rate
|
||
train = np.flatnonzero(guarded_tile_masks(len(expected))[0])
|
||
stride = max(1, len(train) // 8192)
|
||
probe = train[::stride]
|
||
candidates = []
|
||
for fraction in TIMING_GRID:
|
||
positions = nominal + fraction * samples_per_symbol + probe * samples_per_symbol
|
||
samples = lab057.interpolate(matched, positions)
|
||
cfo_hz, h, coherence = fit_cfo_gain(samples, expected[probe], probe, symbol_rate,
|
||
len(expected))
|
||
candidates.append((coherence, float(fraction), cfo_hz, h))
|
||
_, fraction, _, _ = max(candidates, key=lambda row: row[0])
|
||
positions = nominal + fraction * samples_per_symbol + np.arange(len(expected)) * samples_per_symbol
|
||
symbols = lab057.interpolate(matched, positions)
|
||
train_samples = symbols[train]
|
||
cfo_hz, h, coherence = fit_cfo_gain(train_samples, expected[train], train, symbol_rate,
|
||
len(expected))
|
||
return {"timing_fraction": fraction, "timing_samples": fraction * samples_per_symbol,
|
||
"cfo_hz": cfo_hz, "h": h, "train_coherence": coherence,
|
||
"symbols": symbols}
|
||
|
||
|
||
def receiver_like(frontend: lab057.Frontend, values: np.ndarray,
|
||
expected: np.ndarray) -> dict:
|
||
window = frontend.extract(values, len(expected))
|
||
if not window.passed:
|
||
raise RuntimeError(f"Штатное извлечение не выполнено: {window.reason}")
|
||
matched = frontend.filter(values)
|
||
start = frontend.nominal + window.timing_offset_samples
|
||
symbols = matched[start:start + len(expected) * frontend.samples_per_symbol:
|
||
frontend.samples_per_symbol].astype(np.complex128)
|
||
carrier = window.carrier
|
||
h = carrier["abs_h"] * np.exp(1j * carrier["phase_rad"])
|
||
return {"timing_fraction": window.timing_offset_samples / frontend.samples_per_symbol,
|
||
"timing_samples": float(window.timing_offset_samples),
|
||
"cfo_hz": float(carrier["cfo_hz"]), "h": complex(h),
|
||
"train_coherence": None, "symbols": symbols, "matched": matched}
|
||
|
||
|
||
def evaluate(estimate: dict, expected: np.ndarray, symbol_rate: int,
|
||
selected: np.ndarray) -> dict:
|
||
indexes = np.flatnonzero(selected)
|
||
symbols = estimate["symbols"][selected]
|
||
truth = expected[selected]
|
||
times = (indexes + 0.5) / symbol_rate
|
||
derotated = symbols * np.exp(-2j * np.pi * estimate["cfo_hz"] * times)
|
||
h = estimate["h"]
|
||
normalized = derotated / h
|
||
predicted = h * truth
|
||
residual = derotated - predicted
|
||
decisions = normalized.real >= 0
|
||
bits = truth > 0
|
||
errors = int(np.count_nonzero(decisions != bits))
|
||
matched_power = float(np.mean(abs(predicted) ** 2))
|
||
residual_power = float(np.mean(abs(residual) ** 2))
|
||
plus, minus = normalized.real[truth > 0], normalized.real[truth < 0]
|
||
pooled = float(np.sqrt((np.var(plus) + np.var(minus)) / 2))
|
||
separation = (float("inf") if pooled == 0 else
|
||
float((np.mean(plus) - np.mean(minus)) / pooled))
|
||
margin = truth * normalized.real
|
||
return {"errors": errors, "symbols": int(len(indexes)), "ber": errors / len(indexes),
|
||
"plus_mean": float(np.mean(plus)), "plus_std": float(np.std(plus)),
|
||
"minus_mean": float(np.mean(minus)), "minus_std": float(np.std(minus)),
|
||
"cloud_separation_sigma": separation,
|
||
"evm_rms": float(np.sqrt(residual_power / matched_power)),
|
||
"matched_component_power": matched_power, "residual_power": residual_power,
|
||
"matched_residual_db": (float("inf") if residual_power == 0 else
|
||
float(10 * np.log10(matched_power / residual_power))),
|
||
"margin_median": float(np.median(margin)),
|
||
"margin_p10": float(np.percentile(margin, 10)),
|
||
"margin_p01": float(np.percentile(margin, 1)),
|
||
"normalized": normalized, "truth": truth, "indexes": indexes}
|
||
|
||
|
||
def eye_metrics(matched: np.ndarray, nominal: int, estimate: dict, expected: np.ndarray,
|
||
symbol_rate: int) -> tuple[dict, np.ndarray, np.ndarray]:
|
||
samples_per_symbol = SAMPLE_RATE // symbol_rate
|
||
valid = np.flatnonzero(guarded_tile_masks(len(expected))[1])
|
||
valid = valid[(valid >= 2) & (valid < len(expected) - 2)]
|
||
if len(valid) > 4096:
|
||
valid = valid[np.linspace(0, len(valid) - 1, 4096, dtype=int)]
|
||
phases = np.linspace(-0.5, 0.5, 41)
|
||
openings, traces = [], []
|
||
for phase in phases:
|
||
positions = (nominal + (estimate["timing_fraction"] + phase) * samples_per_symbol
|
||
+ valid * samples_per_symbol)
|
||
values = lab057.interpolate(matched, positions)
|
||
times = (valid + 0.5 + phase) / symbol_rate
|
||
normalized = values * np.exp(-2j * np.pi * estimate["cfo_hz"] * times) / estimate["h"]
|
||
truth = expected[valid]
|
||
openings.append(float(np.mean(normalized.real[truth > 0]) -
|
||
np.mean(normalized.real[truth < 0])))
|
||
traces.append(normalized.real)
|
||
openings = np.asarray(openings)
|
||
center = int(np.argmin(abs(phases)))
|
||
left = int(np.argmin(abs(phases + 0.1)))
|
||
right = int(np.argmin(abs(phases - 0.1)))
|
||
metrics = {"opening": float(openings[center]),
|
||
"opening_at_minus_0p1": float(openings[left]),
|
||
"opening_at_plus_0p1": float(openings[right]),
|
||
"drop_minus_0p1_percent": float(100 * (1 - openings[left] / openings[center])),
|
||
"drop_plus_0p1_percent": float(100 * (1 - openings[right] / openings[center]))}
|
||
return metrics, phases, np.asarray(traces).T
|
||
|
||
|
||
def isi_metrics(estimate: dict, expected: np.ndarray, symbol_rate: int) -> tuple[dict, np.ndarray]:
|
||
indexes = np.arange(5, len(expected) - 5)
|
||
train = guarded_tile_masks(len(expected))[0][indexes]
|
||
indexes = indexes[train]
|
||
if len(indexes) > 50000:
|
||
indexes = indexes[np.linspace(0, len(indexes) - 1, 50000, dtype=int)]
|
||
times = (indexes + 0.5) / symbol_rate
|
||
y = estimate["symbols"][indexes] * np.exp(-2j * np.pi * estimate["cfo_hz"] * times)
|
||
matrix = np.column_stack([expected[indexes - offset] for offset in ISI_OFFSETS])
|
||
taps = np.linalg.lstsq(matrix, y, rcond=None)[0]
|
||
main = float(abs(taps[5]) ** 2)
|
||
adjacent = float(np.sum(abs(np.delete(taps, 5)) ** 2))
|
||
return {"main_tap_power": main, "adjacent_taps_power": adjacent,
|
||
"isi_ratio_db": float(10 * np.log10(adjacent / main)),
|
||
"offsets": ISI_OFFSETS.tolist(),
|
||
"taps": [{"real": float(x.real), "imag": float(x.imag)} for x in taps]}, taps
|
||
|
||
|
||
def channel_metrics(tx: np.ndarray, rx: np.ndarray, symbol_rate: int,
|
||
cfo_hz: float) -> tuple[dict, dict]:
|
||
indexes = np.arange(len(tx))
|
||
carrier = np.exp(-2j * np.pi * CARRIER_HZ * indexes / SAMPLE_RATE)
|
||
x = (tx - np.mean(tx)) * carrier
|
||
y = (rx - np.mean(rx)) * carrier * np.exp(-2j * np.pi * cfo_hz * indexes / SAMPLE_RATE)
|
||
frequency, pxx = signal.welch(x, fs=SAMPLE_RATE, window="hann", nperseg=8192,
|
||
noverlap=4096, return_onesided=False)
|
||
_, pyy = signal.welch(y, fs=SAMPLE_RATE, window="hann", nperseg=8192,
|
||
noverlap=4096, return_onesided=False)
|
||
_, pxy = signal.csd(x, y, fs=SAMPLE_RATE, window="hann", nperseg=8192,
|
||
noverlap=4096, return_onesided=False)
|
||
order = np.argsort(frequency)
|
||
frequency, pxx, pyy, pxy = frequency[order], pxx[order], pyy[order], pxy[order]
|
||
transfer = pxy / np.maximum(pxx, 1e-30)
|
||
coherence = abs(pxy) ** 2 / np.maximum(pxx * pyy, 1e-30)
|
||
bandwidth = symbol_rate * (1 + 0.35) / 2
|
||
mask = ((abs(frequency) <= bandwidth) &
|
||
(pxx >= np.max(pxx[abs(frequency) <= bandwidth]) * 0.01))
|
||
f, h, coh = frequency[mask], transfer[mask], coherence[mask]
|
||
gain_db = 20 * np.log10(np.maximum(abs(h), 1e-30))
|
||
phase = np.unwrap(np.angle(h))
|
||
phase_fit = np.polyfit(f, phase, 1)
|
||
phase_residual = phase - np.polyval(phase_fit, f)
|
||
group_delay = -np.gradient(phase, f) / (2 * np.pi)
|
||
stable = bool(len(f) >= 20 and np.median(coh) >= 0.2)
|
||
group_delay_stable = bool(len(f) >= 64 and np.median(coh) >= 0.9)
|
||
metrics = {"stable": stable, "bins": int(len(f)),
|
||
"median_coherence": float(np.median(coh)),
|
||
"amplitude_ripple_p95_p05_db": float(np.percentile(gain_db, 95) - np.percentile(gain_db, 5)),
|
||
"max_min_gain_db": float(np.max(gain_db) - np.min(gain_db)),
|
||
"gain_slope_db_per_mhz": float(np.polyfit(f / 1e6, gain_db, 1)[0]),
|
||
"common_group_delay_seconds": float(-phase_fit[0] / (2 * np.pi)),
|
||
"residual_phase_ripple_p95_p05_rad": float(np.percentile(phase_residual, 95) - np.percentile(phase_residual, 5)),
|
||
"group_delay_stable": group_delay_stable,
|
||
"group_delay_variation_p95_p05_seconds": float(np.percentile(group_delay, 95) - np.percentile(group_delay, 5)) if group_delay_stable else None,
|
||
"band_edge_hz": bandwidth, "tx_psd_threshold_db": -20.0}
|
||
arrays = {"frequency": f, "gain_db": gain_db, "phase_residual": phase_residual,
|
||
"coherence": coh}
|
||
return metrics, arrays
|
||
|
||
|
||
def serializable(value):
|
||
if isinstance(value, dict):
|
||
return {key: serializable(item) for key, item in value.items()
|
||
if key not in ("matched", "normalized", "truth", "indexes")}
|
||
if isinstance(value, list):
|
||
return [serializable(item) for item in value]
|
||
if isinstance(value, complex):
|
||
return {"real": float(value.real), "imag": float(value.imag)}
|
||
if isinstance(value, np.generic):
|
||
return value.item()
|
||
return value
|
||
|
||
|
||
def analyze_values(values: np.ndarray, expected: np.ndarray, symbol_rate: int) -> tuple[dict, dict]:
|
||
frontend = lab057.Frontend(symbol_rate)
|
||
receiver = receiver_like(frontend, values, expected)
|
||
receiver_select = np.arange(len(expected)) >= 1280
|
||
receiver_result = evaluate(receiver, expected, symbol_rate, receiver_select)
|
||
oracle_matched = oracle_matched_filter(values, frontend.taps, symbol_rate, len(expected))
|
||
oracle = oracle_estimate(oracle_matched, frontend.nominal, expected, symbol_rate)
|
||
validation = guarded_tile_masks(len(expected))[1]
|
||
oracle_result = evaluate(oracle, expected, symbol_rate, validation)
|
||
eye, phases, traces = eye_metrics(oracle_matched, frontend.nominal, oracle,
|
||
expected, symbol_rate)
|
||
isi, taps = isi_metrics(oracle, expected, symbol_rate)
|
||
receiver_public = {key: value for key, value in receiver.items() if key not in ("symbols", "matched")}
|
||
oracle_public = {key: value for key, value in oracle.items() if key != "symbols"}
|
||
result = {"receiver_like": {**serializable(receiver_public), "metrics": serializable(receiver_result)},
|
||
"oracle": {**serializable(oracle_public), "metrics": serializable(oracle_result),
|
||
"eye": eye, "isi": isi}}
|
||
working = {"frontend": frontend, "receiver": receiver, "oracle": oracle,
|
||
"oracle_matched": oracle_matched,
|
||
"receiver_result": receiver_result, "oracle_result": oracle_result,
|
||
"eye_phases": phases, "eye_traces": traces, "isi_taps": taps}
|
||
return result, working
|
||
|
||
|
||
def geometric_db(new: float, old: float) -> float:
|
||
return float(20 * np.log10(new / old))
|
||
|
||
|
||
def run_functional_tests() -> None:
|
||
"""Быстрые проверки формул без raw IQ и оборудования."""
|
||
expected = expected_symbols(20_000)
|
||
indexes = np.arange(len(expected))
|
||
h = 0.4 + 0.3j
|
||
estimate = {"symbols": h * expected, "h": h, "cfo_hz": 0.0}
|
||
metrics = evaluate(estimate, expected, 20_000, np.ones(len(expected), bool))
|
||
if metrics["errors"] != 0 or metrics["evm_rms"] > 1e-12:
|
||
raise AssertionError("Точный BPSK control не прошёл")
|
||
isi, _ = isi_metrics(estimate, expected, 20_000)
|
||
if isi["adjacent_taps_power"] > isi["main_tap_power"] * 1e-20:
|
||
raise AssertionError("Точный одно-tap control ошибочно содержит ISI")
|
||
|
||
|
||
def make_plots(rows: list[dict], work: dict[str, dict]) -> None:
|
||
OUTPUT.mkdir(parents=True, exist_ok=True)
|
||
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
|
||
for axis, mode in zip(axes, ("A", "B")):
|
||
for label in (mode + "1", mode + "2"):
|
||
metric = work[label]["oracle_result"]
|
||
values, truth = metric["normalized"], metric["truth"]
|
||
select = np.linspace(0, len(values) - 1, min(5000, len(values)), dtype=int)
|
||
axis.scatter(values.real[select], values.imag[select], s=2, alpha=.2, label=label)
|
||
axis.set(title=f"{mode}: oracle, held-out symbols", xlabel="I / h", ylabel="Q / h", xlim=(-4, 4), ylim=(-4, 4))
|
||
axis.grid(True, alpha=.25); axis.legend()
|
||
fig.tight_layout(); fig.savefig(OUTPUT / "constellation_a_vs_b.png", dpi=160); plt.close(fig)
|
||
|
||
for mode in ("A", "B"):
|
||
fig, axis = plt.subplots(figsize=(8, 4.5))
|
||
for label in (mode + "1", mode + "2"):
|
||
phases, traces = work[label]["eye_phases"], work[label]["eye_traces"]
|
||
select = np.linspace(0, len(traces) - 1, min(180, len(traces)), dtype=int)
|
||
axis.plot(phases, traces[select].T, alpha=.05, color="C0" if label.endswith("1") else "C1")
|
||
axis.axvline(0, color="black", lw=.8); axis.axhline(0, color="black", lw=.8)
|
||
axis.set(title=f"Eye diagram {mode}, oracle alignment", xlabel="Смещение, доли символа", ylabel="Re(r/h)")
|
||
axis.grid(True, alpha=.2); fig.tight_layout(); fig.savefig(OUTPUT / f"eye_{mode.lower()}.png", dpi=160); plt.close(fig)
|
||
|
||
labels = [row["label"] for row in rows]
|
||
figure, axes = plt.subplots(1, 3, figsize=(12, 4))
|
||
axes[0].bar(labels, [100 * row["oracle"]["metrics"]["ber"] for row in rows]); axes[0].set_ylabel("BER, %")
|
||
axes[1].bar(labels, [row["oracle"]["metrics"]["evm_rms"] for row in rows]); axes[1].set_ylabel("RMS EVM")
|
||
axes[2].bar(labels, [row["oracle"]["metrics"]["matched_residual_db"] for row in rows]); axes[2].set_ylabel("Matched/residual, dB")
|
||
for axis in axes: axis.grid(True, axis="y", alpha=.25)
|
||
figure.tight_layout(); figure.savefig(OUTPUT / "symbol_quality_comparison.png", dpi=160); plt.close(figure)
|
||
|
||
fig, axes = plt.subplots(2, 2, figsize=(12, 7), sharex="col")
|
||
for column, mode in enumerate(("A", "B")):
|
||
for label in (mode + "1", mode + "2"):
|
||
arrays = work[label]["channel_arrays"]
|
||
axes[0, column].plot(arrays["frequency"] / 1000, arrays["gain_db"], alpha=.8, label=label)
|
||
axes[1, column].plot(arrays["frequency"] / 1000, arrays["phase_residual"], alpha=.8, label=label)
|
||
axes[0, column].set(title=f"{mode}: усреднённая H(f)", ylabel="|H(f)|, dB")
|
||
axes[1, column].set(xlabel="Частота baseband, кГц", ylabel="Остаточная фаза, rad")
|
||
for row in range(2): axes[row, column].grid(True, alpha=.25); axes[row, column].legend()
|
||
fig.tight_layout(); fig.savefig(OUTPUT / "channel_transfer.png", dpi=160); plt.close(fig)
|
||
|
||
fig, axes = plt.subplots(2, 2, figsize=(10, 7), sharex=True)
|
||
for axis, label in zip(axes.flat, labels):
|
||
taps = work[label]["isi_taps"]
|
||
axis.stem(ISI_OFFSETS, 20 * np.log10(np.maximum(abs(taps / taps[5]), 1e-8)), basefmt=" ")
|
||
axis.set(title=label, ylabel="Tap/main, dB"); axis.grid(True, alpha=.25)
|
||
axes[1, 0].set_xlabel("Смещение, символы"); axes[1, 1].set_xlabel("Смещение, символы")
|
||
fig.tight_layout(); fig.savefig(OUTPUT / "diagnostic_isi.png", dpi=160); plt.close(fig)
|
||
|
||
|
||
def main() -> None:
|
||
OUTPUT.mkdir(parents=True, exist_ok=True)
|
||
metadata = json.loads((RAW / "tx_metadata.json").read_text(encoding="utf-8"))
|
||
analysis = json.loads((LAB056 / "rf_analysis.json").read_text(encoding="utf-8"))
|
||
integrity = {"tx": {"path": str(TX_PATH.relative_to(PROJECT_ROOT)), "bytes": TX_PATH.stat().st_size,
|
||
"sha256": sha256(TX_PATH), "expected_sha256": TX_SHA256},
|
||
"rx": {"path": str(RX_PATH.relative_to(PROJECT_ROOT)), "bytes": RX_PATH.stat().st_size,
|
||
"sha256": sha256(RX_PATH), "expected_sha256": RX_SHA256}}
|
||
if integrity["tx"]["sha256"] != TX_SHA256 or integrity["rx"]["sha256"] != RX_SHA256:
|
||
raise RuntimeError("SHA-256 исходных данных Lab056 не совпал")
|
||
tx_file, rx_file = Cs8File(TX_PATH), Cs8File(RX_PATH)
|
||
measured = {row["label"]: row for row in analysis["segments"]}
|
||
rows, work = [], {}
|
||
for segment in metadata["segments"]:
|
||
label, rate = segment["label"], segment["symbol_rate"]
|
||
expected = expected_symbols(rate)
|
||
tx = tx_file.read(segment["start"], segment["samples"])
|
||
rx_start = measured[label]["measurement"]["start"]
|
||
rx = rx_file.read(rx_start, segment["samples"])
|
||
result, working = analyze_values(rx, expected, rate)
|
||
trim = 21 * (SAMPLE_RATE // rate)
|
||
channel, arrays = channel_metrics(tx[trim:-trim], rx[trim:-trim], rate,
|
||
result["oracle"]["cfo_hz"])
|
||
result["oracle"]["channel"] = channel
|
||
row = {"label": label, "symbol_rate": rate, "samples_per_symbol": SAMPLE_RATE // rate,
|
||
"rx_start": rx_start, **result}
|
||
rows.append(row); working["channel_arrays"] = arrays; work[label] = working
|
||
|
||
controls = {"tx": {}, "awgn": {}}
|
||
for rate, label in ((20_000, "A"), (500_000, "B")):
|
||
segment = next(item for item in metadata["segments"] if item["label"] == label + "1")
|
||
tx = tx_file.read(segment["start"], segment["samples"])
|
||
expected = expected_symbols(rate)
|
||
tx_result, _ = analyze_values(tx, expected, rate)
|
||
controls["tx"][label] = tx_result
|
||
random_generator = np.random.default_rng(AWGN_SEED)
|
||
noise = np.sqrt(AWGN_COMPLEX_VARIANCE / 2) * (
|
||
random_generator.normal(size=len(tx)) + 1j * random_generator.normal(size=len(tx)))
|
||
awgn_result, _ = analyze_values((tx + noise).astype(np.complex64), expected, rate)
|
||
controls["awgn"][label] = awgn_result
|
||
controls["awgn"].update({"seed": AWGN_SEED, "complex_noise_variance_per_sample": AWGN_COMPLEX_VARIANCE,
|
||
"invariant": "same Fs, TX RMS and sampled complex-noise variance"})
|
||
|
||
background_start = analysis["alignment"]["start"] - metadata["segments"][1]["samples"] - 100_000
|
||
background_values = rx_file.read(background_start, metadata["segments"][1]["samples"])
|
||
background_result, _ = analyze_values(background_values, expected_symbols(500_000), 500_000)
|
||
controls["background"] = {"rx_start": background_start, **background_result}
|
||
if any(controls["tx"][mode]["oracle"]["metrics"]["ber"] != 0 for mode in ("A", "B")):
|
||
raise RuntimeError("TX waveform control BER не равен нулю")
|
||
if not (controls["awgn"]["B"]["oracle"]["metrics"]["evm_rms"] >
|
||
controls["awgn"]["A"]["oracle"]["metrics"]["evm_rms"]):
|
||
raise RuntimeError("AWGN control не показал ожидаемого ухудшения EVM")
|
||
background_oracle = controls["background"]["oracle"]
|
||
if not (0.4 < background_oracle["metrics"]["ber"] < 0.6 and
|
||
background_oracle["train_coherence"] < 0.01):
|
||
raise RuntimeError("Background ошибочно выглядит как качественный BPSK")
|
||
|
||
def mean_metric(mode: str, path: tuple[str, ...]) -> float:
|
||
values = []
|
||
for row in rows:
|
||
if row["label"].startswith(mode):
|
||
value = row
|
||
for key in path: value = value[key]
|
||
values.append(value)
|
||
return float(np.mean(values))
|
||
|
||
primary_paths = {
|
||
"symbol_ber": ("oracle", "metrics", "ber"), "evm_rms": ("oracle", "metrics", "evm_rms"),
|
||
"matched_residual_db": ("oracle", "metrics", "matched_residual_db"),
|
||
"decision_margin_median": ("oracle", "metrics", "margin_median"),
|
||
"eye_opening": ("oracle", "eye", "opening"), "isi_ratio_db": ("oracle", "isi", "isi_ratio_db"),
|
||
"channel_ripple_db": ("oracle", "channel", "amplitude_ripple_p95_p05_db")}
|
||
comparison = {}
|
||
for name, path in primary_paths.items():
|
||
a, b = mean_metric("A", path), mean_metric("B", path)
|
||
comparison[name] = {"a_mean": a, "b_mean": b, "b_minus_a": b - a}
|
||
repeatability = {}
|
||
for mode in ("A", "B"):
|
||
first = next(row for row in rows if row["label"] == mode + "1")
|
||
second = next(row for row in rows if row["label"] == mode + "2")
|
||
repeatability[mode] = {}
|
||
for name, path in primary_paths.items():
|
||
one, two = first, second
|
||
for key in path: one = one[key]; two = two[key]
|
||
entry = {"absolute_second_minus_first": float(two - one),
|
||
"relative_second_minus_first": None if one == 0 else float((two - one) / abs(one))}
|
||
if name in ("evm_rms", "eye_opening") and one > 0 and two > 0:
|
||
entry["difference_db"] = geometric_db(two, one)
|
||
elif name == "symbol_ber" and one > 0 and two > 0:
|
||
entry["difference_db"] = float(10 * np.log10(two / one))
|
||
elif name.endswith("_db"):
|
||
entry["difference_db"] = float(two - one)
|
||
else: entry["difference_db"] = None
|
||
repeatability[mode][name] = entry
|
||
|
||
summary = {"lab": "Lab058", "status": "COMPLETED", "scope": "OFFLINE_SAVED_LAB056_ONLY",
|
||
"rf_tx": 0, "hardware_access": 0, "integrity": integrity,
|
||
"method": {"tiles": TILE_COUNT, "training_tiles": "even", "validation_tiles": "odd",
|
||
"guard_symbols_each_tile_edge": 8,
|
||
"oracle_dc_source": "guarded_even_training_raw_intervals_only",
|
||
"timing_grid_symbol": TIMING_GRID.tolist(), "isi_offsets": ISI_OFFSETS.tolist()},
|
||
"segments": serializable(rows), "comparison": comparison,
|
||
"repeatability": repeatability, "controls": serializable(controls),
|
||
"known_energy_per_symbol_difference_db": EXPECTED_ENERGY_DIFFERENCE_DB,
|
||
"lab054_historical_result": "NOT_A, 0/16", "lab058_commit": "NOT_PERFORMED",
|
||
"push": "NOT_PERFORMED"}
|
||
(OUTPUT / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2, allow_nan=False) + "\n", encoding="utf-8")
|
||
|
||
flat_rows = []
|
||
for row in rows:
|
||
for mode in ("receiver_like", "oracle"):
|
||
metric = row[mode]["metrics"]
|
||
flat_rows.append({"label": row["label"], "symbol_rate": row["symbol_rate"], "analysis_mode": mode,
|
||
**{key: metric[key] for key in ("errors", "symbols", "ber", "plus_mean", "plus_std",
|
||
"minus_mean", "minus_std", "cloud_separation_sigma", "evm_rms",
|
||
"matched_component_power", "residual_power", "matched_residual_db",
|
||
"margin_median", "margin_p10", "margin_p01")},
|
||
"timing_fraction": row[mode]["timing_fraction"], "cfo_hz": row[mode]["cfo_hz"],
|
||
"eye_opening": row.get("oracle", {}).get("eye", {}).get("opening") if mode == "oracle" else None,
|
||
"isi_ratio_db": row.get("oracle", {}).get("isi", {}).get("isi_ratio_db") if mode == "oracle" else None,
|
||
"channel_ripple_db": row.get("oracle", {}).get("channel", {}).get("amplitude_ripple_p95_p05_db") if mode == "oracle" else None})
|
||
with (OUTPUT / "segment_metrics.csv").open("w", newline="", encoding="utf-8") as stream:
|
||
writer = csv.DictWriter(stream, fieldnames=list(flat_rows[0])); writer.writeheader(); writer.writerows(flat_rows)
|
||
with (OUTPUT / "channel_response.csv").open("w", newline="", encoding="utf-8") as stream:
|
||
writer = csv.writer(stream); writer.writerow(["label", "frequency_hz", "gain_db", "phase_residual_rad", "coherence"])
|
||
for label, item in work.items():
|
||
arrays = item["channel_arrays"]
|
||
writer.writerows(zip([label] * len(arrays["frequency"]), arrays["frequency"], arrays["gain_db"], arrays["phase_residual"], arrays["coherence"]))
|
||
make_plots(rows, work)
|
||
print(json.dumps({"status": summary["status"], "comparison": comparison,
|
||
"receiver_ber": {r["label"]: r["receiver_like"]["metrics"]["ber"] for r in rows},
|
||
"oracle_ber": {r["label"]: r["oracle"]["metrics"]["ber"] for r in rows},
|
||
"controls": {"tx_ber": {k: v["oracle"]["metrics"]["ber"] for k, v in controls["tx"].items()},
|
||
"awgn_ber": {k: v["oracle"]["metrics"]["ber"] for k, v in controls["awgn"].items() if isinstance(v, dict)},
|
||
"background_ber": controls["background"]["oracle"]["metrics"]["ber"]}}, indent=2))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|