"""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()