"""Lab059: пределы preamble-only восстановления несущей по сохранённым IQ.""" from __future__ import annotations import csv from dataclasses import asdict, dataclass, replace import hashlib import json import math from pathlib import Path import sys from typing import Iterable if __package__ in (None, ""): sys.path.insert(0, str(Path(__file__).resolve().parents[1])) import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np from experiments import lab054_hardware_control as lab054 from experiments import lab057_symbol_diagnosis as lab057 from experiments import lab058_abab_symbol_quality as lab058 ROOT = Path(__file__).resolve().parents[1] OUTPUT = ROOT / "data/processed/lab059" METHODS = ("P0", "P1", "P2", "P3") SAMPLE_RATE = 8_000_000 PREAMBLE_SYMBOLS = 1280 PROBE_SYMBOLS = 1424 FULL_BLOCK_SYMBOLS = 20 P2_BLOCK_SYMBOLS = 16 P2_TRAIN = np.arange(32, 608) P2_VALIDATION = np.arange(672, 1248) BURST_DURATION_SECONDS = 0.2048 MONTE_CARLO_TRIALS = 300 SYNTHETIC_FREQUENCIES_HZ = np.array( [-100, -20, -10, -5, -2, -1, -0.5, -0.25, 0, 0.25, 0.5, 1, 2, 5, 10, 20, 100], dtype=np.float64, ) QUALITY_LEVELS = ( ("LAB058_A_EVM", 0.317390), ("INTERMEDIATE_EVM", 0.55), ("LAB058_B_EVM", 0.845350), ) CATASTROPHIC_ERROR_HZ = 1.0 / (2.0 * BURST_DURATION_SECONDS) @dataclass(frozen=True) class CarrierEstimate: """Оценка несущей с фазой, отнесённой к t=0.""" method: str frequency_hz: float gain: complex estimation_symbols: int estimation_duration_seconds: float phase_fit_rmse_rad: float mean_block_coherence: float unambiguous_abs_hz: float validation_phase_rmse_rad: float | None = None validation_errors: int | None = None validation_symbols: int | None = None def serializable(value): """Преобразовать numpy и complex в строгий JSON.""" if isinstance(value, dict): return {str(key): serializable(item) for key, item in value.items()} if isinstance(value, (list, tuple)): return [serializable(item) for item in value] if isinstance(value, np.ndarray): return serializable(value.tolist()) if isinstance(value, complex): return {"real": float(value.real), "imag": float(value.imag)} if isinstance(value, np.generic): return serializable(value.item()) return value def write_json(path: Path, value) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text( json.dumps(serializable(value), ensure_ascii=False, indent=2, allow_nan=False) + "\n", encoding="utf-8", ) 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() def phase_slope_to_hz(slope_rad_per_second: float, *, squared_bpsk: bool) -> float: """Перевести наклон фазы в Гц с учётом удвоения при BPSK².""" divisor = 4.0 * np.pi if squared_bpsk else 2.0 * np.pi return float(slope_rad_per_second / divisor) def phase_drift_degrees(frequency_error_hz: float, duration_seconds: float) -> float: """Накопленный фазовый уход от ошибки частоты.""" return float(360.0 * frequency_error_hz * duration_seconds) def maximum_frequency_error_hz(phase_limit_degrees: float, duration_seconds: float) -> float: """Допустимая ошибка частоты для заданного фазового ухода.""" return float(phase_limit_degrees / (360.0 * duration_seconds)) def awgn_cfo_crlb_standard_deviation_hz( evm_rms: float, symbol_count: int, symbol_rate: float, ) -> float: """CRLB СКО CFO для известных символов, неизвестной фазы и complex AWGN. Здесь EVM² = E|n|²/|h|². Это модельная проверка масштаба, а не описание частотно-селективного OTA-остатка. """ variance = ( 6.0 * evm_rms**2 * symbol_rate**2 / ((2.0 * np.pi) ** 2 * symbol_count * (symbol_count**2 - 1)) ) return float(np.sqrt(variance)) def _huber_phase_line( block_values: np.ndarray, block_times_seconds: np.ndarray, ) -> tuple[float, float, float, float]: """Устойчиво оценить phase(t)=phi0+slope*t заранее фиксированным IRLS.""" values = np.asarray(block_values, dtype=np.complex128) times = np.asarray(block_times_seconds, dtype=np.float64) if values.ndim != 1 or times.ndim != 1 or len(values) != len(times) or len(values) < 3: raise ValueError("Недостаточно фазовых блоков") phases = np.unwrap(np.angle(values)) base_weights = np.maximum(np.abs(values), 1e-12) center_time = float(np.sum(base_weights * times) / np.sum(base_weights)) centered = times - center_time design = np.column_stack((np.ones(len(times)), centered)) weights = base_weights.copy() coefficients = np.zeros(2, dtype=np.float64) for _iteration in range(12): root = np.sqrt(weights / max(float(np.max(weights)), 1e-30)) updated, *_unused = np.linalg.lstsq( design * root[:, None], phases * root, rcond=None ) residual = phases - design @ updated median = float(np.median(residual)) scale = 1.4826 * float(np.median(np.abs(residual - median))) + 1e-9 limit = 1.345 * scale robust = np.ones_like(residual) outside = np.abs(residual) > limit robust[outside] = limit / np.abs(residual[outside]) weights = base_weights * robust if np.allclose(updated, coefficients, rtol=0.0, atol=1e-12): coefficients = updated break coefficients = updated residual = phases - design @ coefficients rmse = float(np.sqrt(np.sum(base_weights * residual**2) / np.sum(base_weights))) phase_at_zero = float(coefficients[0] - coefficients[1] * center_time) coherence = float(abs(np.sum(values)) / max(float(np.sum(abs(values))), 1e-30)) return phase_at_zero, float(coefficients[1]), rmse, coherence def _block_values( values: np.ndarray, indexes: np.ndarray, symbol_rate: int, block_symbol_count: int, ) -> tuple[np.ndarray, np.ndarray]: """Собрать неперекрывающиеся когерентные блоки по исходной шкале времени.""" values = np.asarray(values, dtype=np.complex128) indexes = np.asarray(indexes, dtype=np.int64) block_count = len(values) // block_symbol_count if block_count < 3: raise ValueError("Слишком мало когерентных блоков") usable = block_count * block_symbol_count reshaped_values = values[:usable].reshape(block_count, block_symbol_count) reshaped_indexes = indexes[:usable].reshape(block_count, block_symbol_count) return ( np.mean(reshaped_values, axis=1), np.mean((reshaped_indexes + 0.5) / symbol_rate, axis=1), ) def _gain_at_zero( symbols: np.ndarray, known_symbols: np.ndarray, indexes: np.ndarray, symbol_rate: int, frequency_hz: float, ) -> complex: times = (indexes + 0.5) / symbol_rate return complex( np.mean( symbols * np.conj(known_symbols) * np.exp(-2j * np.pi * frequency_hz * times) ) ) def _validation_metrics( symbols: np.ndarray, known_symbols: np.ndarray, indexes: np.ndarray, symbol_rate: int, frequency_hz: float, gain: complex, ) -> tuple[float, int]: times = (indexes + 0.5) / symbol_rate normalized = symbols * np.exp(-2j * np.pi * frequency_hz * times) / gain phase_error = np.angle(normalized * np.conj(known_symbols)) errors = int(np.count_nonzero((normalized.real >= 0) != (known_symbols.real >= 0))) return float(np.sqrt(np.mean(phase_error**2))), errors def estimate_carrier( preamble_symbols: np.ndarray, symbol_rate: int, method: str, ) -> CarrierEstimate: """Оценить частоту и фазу только по 1280 известным символам преамбулы.""" received = np.asarray(preamble_symbols, dtype=np.complex128) known = lab057.EXPECTED.astype(np.complex128) if received.shape != (PREAMBLE_SYMBOLS,): raise ValueError("Оценивателю требуется ровно 1280 символов преамбулы") all_indexes = np.arange(PREAMBLE_SYMBOLS, dtype=np.int64) if method == "P0": statistics = lab057.squared_cfo_statistics(received, symbol_rate) frequency = float(statistics["cfo_hz"]) gain = _gain_at_zero(received, known, all_indexes, symbol_rate, frequency) phase_rmse, validation_errors = _validation_metrics( received[P2_VALIDATION], known[P2_VALIDATION], P2_VALIDATION, symbol_rate, frequency, gain, ) blocks = (received[:1024] ** 2).reshape(2, 512).mean(axis=1) return CarrierEstimate( method, frequency, gain, 1024, 1024 / symbol_rate, 0.0, float(abs(np.sum(blocks)) / max(float(np.sum(abs(blocks))), 1e-30)), symbol_rate / 2048.0, phase_rmse, validation_errors, len(P2_VALIDATION), ) if method == "P2": indexes = P2_TRAIN despread = received[indexes] * np.conj(known[indexes]) blocks, times = _block_values(despread, indexes, symbol_rate, P2_BLOCK_SYMBOLS) _phase0, slope, rmse, coherence = _huber_phase_line(blocks, times) frequency = phase_slope_to_hz(slope, squared_bpsk=False) gain = _gain_at_zero(received[indexes], known[indexes], indexes, symbol_rate, frequency) phase_rmse, validation_errors = _validation_metrics( received[P2_VALIDATION], known[P2_VALIDATION], P2_VALIDATION, symbol_rate, frequency, gain, ) return CarrierEstimate( method, frequency, gain, len(indexes), (indexes[-1] - indexes[0] + 1) / symbol_rate, rmse, coherence, symbol_rate / (2.0 * P2_BLOCK_SYMBOLS), phase_rmse, validation_errors, len(P2_VALIDATION), ) if method == "P1": transformed = received * np.conj(known) squared_bpsk = False elif method == "P3": transformed = received**2 squared_bpsk = True else: raise ValueError(f"Неизвестный метод {method}") blocks, times = _block_values( transformed, all_indexes, symbol_rate, FULL_BLOCK_SYMBOLS ) _phase0, slope, rmse, coherence = _huber_phase_line(blocks, times) frequency = phase_slope_to_hz(slope, squared_bpsk=squared_bpsk) gain = _gain_at_zero(received, known, all_indexes, symbol_rate, frequency) phase_rmse, validation_errors = _validation_metrics( received[P2_VALIDATION], known[P2_VALIDATION], P2_VALIDATION, symbol_rate, frequency, gain, ) ambiguity_divisor = 4.0 if squared_bpsk else 2.0 return CarrierEstimate( method, frequency, gain, PREAMBLE_SYMBOLS, PREAMBLE_SYMBOLS / symbol_rate, rmse, coherence, symbol_rate / (ambiguity_divisor * FULL_BLOCK_SYMBOLS), phase_rmse, validation_errors, len(P2_VALIDATION), ) def estimate_from_burst( burst_symbols: np.ndarray, symbol_rate: int, method: str, ) -> CarrierEstimate: """Явная граница: оцениватель видит только prefix-преамбулу.""" return estimate_carrier(np.asarray(burst_symbols)[:PREAMBLE_SYMBOLS], symbol_rate, method) def compensate( symbols: np.ndarray, symbol_rate: int, estimate: CarrierEstimate, ) -> np.ndarray: """Один раз компенсировать CFO и фазу по единой шкале времени.""" times = (np.arange(len(symbols), dtype=np.float64) + 0.5) / symbol_rate return symbols * np.exp(-2j * np.pi * estimate.frequency_hz * times) / estimate.gain def _extract_symbols( matched: np.ndarray, nominal: int, timing_samples: int, count: int, samples_per_symbol: int, ) -> np.ndarray: values = matched[ nominal + timing_samples: nominal + timing_samples + count * samples_per_symbol: samples_per_symbol ].astype(np.complex128) if len(values) != count: raise ValueError("Недостаточно отфильтрованных символов") return values def local_preamble_symbols( raw_values: np.ndarray, frontend: lab057.Frontend, timing_samples: int, selected_indexes: np.ndarray, ) -> np.ndarray: """RRC-точки, чьи raw DC и FIR-опора ограничены выбранной преамбулой.""" indexes = np.asarray(selected_indexes, dtype=np.int64) if indexes.ndim != 1 or not len(indexes) or indexes.min() < 0 or indexes.max() >= PREAMBLE_SYMBOLS: raise ValueError("Некорректные индексы локальной преамбулы") positions = ( frontend.nominal + int(timing_samples) + indexes * frontend.samples_per_symbol ) low = int(positions.min() - len(frontend.taps) + 1) high = int(positions.max() + 1) if low < 0 or high > len(raw_values): raise ValueError("Сырая опора преамбулы выходит за окно") local = np.asarray(raw_values[low:high], dtype=np.complex64) rotation = lab057.packet.lab043._carrier_downconversion_vector(low, len(local)) filtered = lab057.oaconvolve( (local - np.mean(local)) * rotation, frontend.taps, mode="full", ) output = np.zeros(PREAMBLE_SYMBOLS, dtype=np.complex128) output[indexes] = filtered[positions - low] return output def real_carrier_estimates( raw_values: np.ndarray, frontend: lab057.Frontend, timing_samples: int, stock_symbols: np.ndarray, symbol_rate: int, methods: Iterable[str], ) -> dict[str, CarrierEstimate]: """P0 остаётся stock; P1–P3 не читают raw после своей преамбульной опоры.""" requested = tuple(methods) estimates: dict[str, CarrierEstimate] = {} if "P0" in requested: estimates["P0"] = estimate_from_burst(stock_symbols, symbol_rate, "P0") full_local = None if any(method in requested for method in ("P1", "P2", "P3")): full_local = local_preamble_symbols( raw_values, frontend, timing_samples, np.arange(PREAMBLE_SYMBOLS) ) for method in ("P1", "P3"): if method in requested: estimates[method] = estimate_carrier(full_local, symbol_rate, method) if "P2" in requested: train_local = local_preamble_symbols( raw_values, frontend, timing_samples, P2_TRAIN ) estimate = estimate_carrier(train_local, symbol_rate, "P2") validation_rmse, validation_errors = _validation_metrics( full_local[P2_VALIDATION], lab057.EXPECTED[P2_VALIDATION], P2_VALIDATION, symbol_rate, estimate.frequency_hz, estimate.gain, ) estimates["P2"] = replace( estimate, validation_phase_rmse_rad=validation_rmse, validation_errors=validation_errors, ) return {method: estimates[method] for method in requested} def _estimate_dict(estimate: CarrierEstimate) -> dict: return asdict(estimate) def _lab056_inputs() -> tuple[list[dict], dict[str, float]]: metadata = json.loads((ROOT / "data/raw/lab056/tx_metadata.json").read_text(encoding="utf-8")) analysis = json.loads((ROOT / "data/processed/lab056/rf_analysis.json").read_text(encoding="utf-8")) lab058_summary = json.loads((ROOT / "data/processed/lab058/summary.json").read_text(encoding="utf-8")) measured = {row["label"]: row for row in analysis["segments"]} oracle = {row["label"]: row["oracle"]["cfo_hz"] for row in lab058_summary["segments"]} inputs = [] for segment in metadata["segments"]: inputs.append({ "label": segment["label"], "symbol_rate": segment["symbol_rate"], "samples": segment["samples"], "rx_start": measured[segment["label"]]["measurement"]["start"], }) return inputs, oracle def analyze_lab056(methods: Iterable[str]) -> list[dict]: """Применить оценки к полным A1/B1/A2/B2, truth только для метрик после fit.""" reader = lab058.Cs8File(ROOT / "data/raw/lab056/rx_stream.cs8") inputs, oracle = _lab056_inputs() rows: list[dict] = [] for item in inputs: rate = item["symbol_rate"] expected = lab058.expected_symbols(rate) values = reader.read(item["rx_start"], item["samples"]) frontend = lab057.Frontend(rate) stock = frontend.extract(values, len(expected)) if not stock.passed: raise RuntimeError(f"Lab056 {item['label']}: штатное извлечение {stock.reason}") matched = frontend.filter(values) symbols = _extract_symbols( matched, frontend.nominal, stock.timing_offset_samples, len(expected), frontend.samples_per_symbol, ) estimates = real_carrier_estimates( values, frontend, stock.timing_offset_samples, symbols, rate, methods ) if "P0" in estimates: if not np.isclose(estimates["P0"].frequency_hz, stock.carrier["cfo_hz"], atol=1e-10): raise AssertionError("P0 Lab056 не воспроизвёл штатную оценку CFO") selected = np.arange(len(expected)) >= PREAMBLE_SYMBOLS for method, estimate in estimates.items(): result = lab058.evaluate( {"symbols": symbols, "h": estimate.gain, "cfo_hz": estimate.frequency_hz}, expected, rate, selected, ) reference = float(oracle[item["label"]]) rows.append({ "dataset": "Lab056", "label": item["label"], "frame": None, "symbol_rate": rate, "method": method, "timing_samples": int(stock.timing_offset_samples), "frequency_hz": estimate.frequency_hz, "diagnostic_reference_hz": reference, "frequency_error_hz": estimate.frequency_hz - reference, "estimation_symbols": estimate.estimation_symbols, "estimation_duration_seconds": estimate.estimation_duration_seconds, "phase_drift_to_burst_end_degrees": phase_drift_degrees( estimate.frequency_hz - reference, BURST_DURATION_SECONDS ), "preamble_errors": int(np.count_nonzero( (compensate(symbols[:PREAMBLE_SYMBOLS], rate, estimate).real >= 0) != lab057.PREAMBLE )), "validation_errors": estimate.validation_errors, "validation_symbols": estimate.validation_symbols, "validation_phase_rmse_rad": estimate.validation_phase_rmse_rad, "header_raw_errors": None, "header_voted_errors": None, "HEADER_ACCEPTED": None, "FEC_STARTED": None, "CRC_PASS": None, "PACKET_OUTPUT": None, "symbol_errors": result["errors"], "symbols": result["symbols"], "symbol_error_rate": result["ber"], "evm_rms": result["evm_rms"], "matched_residual_db": result["matched_residual_db"], "estimate": _estimate_dict(estimate), }) return rows def _packet_sources() -> tuple[list[dict], list[dict], lab058.Cs8File, lab058.Cs8File]: mapping = json.loads((ROOT / "data/processed/lab057/candidate_mapping.json").read_text(encoding="utf-8")) lab054_items = [row for row in mapping if row["true_start_frame"]] audit = json.loads((ROOT / "data/processed/lab055/energy_audit_v2.json").read_text(encoding="utf-8")) old_directory = Path(audit["experiments"]["Lab047"]["directory"]) old = json.loads((old_directory / "independent_reprocess_result_corrected.json").read_text(encoding="utf-8")) lab047_items = [row for row in old["strict"]["outcomes"] if row["accepted"]] return ( lab054_items, lab047_items, lab058.Cs8File(ROOT / "data/raw/lab054/rx_stream.cs8"), lab058.Cs8File(old_directory / "rx_stream.cs8"), ) def _packet_rows( reader: lab058.Cs8File, start: int, symbol_rate: int, dataset: str, frame_number: int, expected_header: np.ndarray, methods: Iterable[str], ) -> list[dict]: frontend = lab057.Frontend(symbol_rate) probe_count = PROBE_SYMBOLS probe_values = reader.read(start, (probe_count + 42) * frontend.samples_per_symbol) stock = frontend.extract(probe_values, probe_count) if not stock.passed: raise RuntimeError(f"{dataset} frame {frame_number}: {stock.reason}") matched = frontend.filter(probe_values) symbols = _extract_symbols( matched, frontend.nominal, stock.timing_offset_samples, probe_count, frontend.samples_per_symbol, ) estimates = real_carrier_estimates( probe_values, frontend, stock.timing_offset_samples, symbols, symbol_rate, methods ) if "P0" in estimates and not np.isclose( estimates["P0"].frequency_hz, stock.carrier["cfo_hz"], atol=1e-10 ): raise AssertionError(f"P0 {dataset} не воспроизвёл штатную оценку") rows = [] for method, estimate in estimates.items(): normalized_probe = compensate(symbols, symbol_rate, estimate) header = lab057.physical.decode_phy_header_soft( normalized_probe.real[PREAMBLE_SYMBOLS:PROBE_SYMBOLS] ) full_soft = normalized_probe.real[PREAMBLE_SYMBOLS:] if header["passed"]: packet_bytes = int(header["packet_bytes"]) count = ( PREAMBLE_SYMBOLS + 144 + lab057.physical.conv.encoded_bit_count(packet_bytes * 8) ) full_values = reader.read(start, (count + 42) * frontend.samples_per_symbol) if method == "P0": full_stock = frontend.extract(full_values, count) if not full_stock.passed: raise RuntimeError(f"{dataset} P0 full extract: {full_stock.reason}") full_soft = full_stock.normalized_real[PREAMBLE_SYMBOLS:] else: full_matched = frontend.filter(full_values) full_symbols = _extract_symbols( full_matched, frontend.nominal, stock.timing_offset_samples, count, frontend.samples_per_symbol, ) full_soft = compensate(full_symbols, symbol_rate, estimate).real[PREAMBLE_SYMBOLS:] decoded, _packet = lab057.decode_body(full_soft) decisions = normalized_probe.real >= 0 expected_header_bool = expected_header.astype(bool) voted = normalized_probe.real[PREAMBLE_SYMBOLS:PROBE_SYMBOLS].reshape(48, 3).sum(axis=1) >= 0 rows.append({ "dataset": dataset, "label": f"{dataset}_{frame_number:02d}", "frame": frame_number, "symbol_rate": symbol_rate, "method": method, "timing_samples": int(stock.timing_offset_samples), "frequency_hz": estimate.frequency_hz, "diagnostic_reference_hz": None, "frequency_error_hz": None, "estimation_symbols": estimate.estimation_symbols, "estimation_duration_seconds": estimate.estimation_duration_seconds, "phase_drift_to_burst_end_degrees": phase_drift_degrees( estimate.frequency_hz, BURST_DURATION_SECONDS ), "preamble_errors": int(np.count_nonzero( decisions[:PREAMBLE_SYMBOLS] != lab057.PREAMBLE )), "validation_errors": estimate.validation_errors, "validation_symbols": estimate.validation_symbols, "validation_phase_rmse_rad": estimate.validation_phase_rmse_rad, "header_raw_errors": int(np.count_nonzero( decisions[PREAMBLE_SYMBOLS:PROBE_SYMBOLS] != expected_header_bool )), "header_voted_errors": int(np.count_nonzero( voted != expected_header_bool[::3] )), "HEADER_ACCEPTED": bool(decoded["HEADER_ACCEPTED"]), "FEC_STARTED": bool(decoded["FEC_STARTED"]), "CRC_PASS": decoded["CRC_RESULT"] == lab057.physical.CRC_PASS, "PACKET_OUTPUT": bool(decoded["PACKET_OUTPUT"]), "symbol_errors": None, "symbols": None, "symbol_error_rate": None, "evm_rms": None, "matched_residual_db": None, "estimate": _estimate_dict(estimate), }) return rows def analyze_packets(methods: Iterable[str]) -> list[dict]: lab054_items, lab047_items, lab054_reader, lab047_reader = _packet_sources() rows: list[dict] = [] plans = lab054.plans() for item in lab054_items: frame = int(item["true_start_frame"]) rows.extend(_packet_rows( lab054_reader, int(item["refined_start"]), 500_000, "Lab054", frame, plans[frame - 1].plan.phy_header_bits, methods, )) expected_header_047 = lab057.physical.encode_phy_header(268) for frame, item in enumerate(lab047_items, 1): rows.extend(_packet_rows( lab047_reader, int(item["refined_start_sample"]), 20_000, "Lab047", frame, expected_header_047, methods, )) return rows def _clean_synthetic_preamble(symbol_rate: int) -> tuple[np.ndarray, dict]: """Сформировать чистые точки через фактические Fs и TX/RX RRC.""" values, truth = lab057.synthetic_waveform( symbol_rate, delay_symbols=0.0, cfo_hz=0.0, phase_rad=0.0, es_n0_db=None ) frontend = lab057.Frontend(symbol_rate) probe = values[: (PROBE_SYMBOLS + 42) * frontend.samples_per_symbol] stock = frontend.extract(probe, PROBE_SYMBOLS) matched = frontend.filter(probe) symbols = _extract_symbols( matched, frontend.nominal, stock.timing_offset_samples, PREAMBLE_SYMBOLS, frontend.samples_per_symbol, ) indexes = np.arange(PREAMBLE_SYMBOLS) gain = _gain_at_zero(symbols, lab057.EXPECTED, indexes, symbol_rate, 0.0) normalized = symbols / gain return normalized, { "sample_rate_hz": SAMPLE_RATE, "symbol_rate": symbol_rate, "samples_per_symbol": frontend.samples_per_symbol, "rrc_taps": len(frontend.taps), "tx_truth": truth, "clean_evm_rms": float(np.sqrt(np.mean(abs(normalized - lab057.EXPECTED) ** 2))), } def synthetic_study() -> tuple[list[dict], list[dict], dict, dict]: """Воспроизводимый Monte Carlo с фиксированными частотами, EVM и seed.""" aggregate_rows: list[dict] = [] by_cfo_rows: list[dict] = [] distributions: dict[tuple[int, str, str], np.ndarray] = {} clean_metadata = {} for symbol_rate in (20_000, 500_000): clean, clean_info = _clean_synthetic_preamble(symbol_rate) clean_metadata[str(symbol_rate)] = clean_info times = (np.arange(PREAMBLE_SYMBOLS) + 0.5) / symbol_rate for quality_index, (quality_name, evm_rms) in enumerate(QUALITY_LEVELS): method_errors = {method: [] for method in METHODS} method_validation = {method: [] for method in METHODS} for frequency_index, true_frequency in enumerate(SYNTHETIC_FREQUENCIES_HZ): seed = 59_000_000 + symbol_rate + quality_index * 100_000 + frequency_index * 1_000 random_generator = np.random.default_rng(seed) errors_by_method = {method: [] for method in METHODS} validation_by_method = {method: [] for method in METHODS} for _trial in range(MONTE_CARLO_TRIALS): phase = random_generator.uniform(-np.pi, np.pi) noise = evm_rms / np.sqrt(2.0) * ( random_generator.normal(size=PREAMBLE_SYMBOLS) + 1j * random_generator.normal(size=PREAMBLE_SYMBOLS) ) received = ( clean * np.exp(1j * (phase + 2 * np.pi * true_frequency * times)) + noise ) for method in METHODS: estimate = estimate_carrier(received, symbol_rate, method) error = estimate.frequency_hz - true_frequency errors_by_method[method].append(error) validation_by_method[method].append(estimate.validation_phase_rmse_rad) for method in METHODS: errors = np.asarray(errors_by_method[method], dtype=np.float64) validation = np.asarray(validation_by_method[method], dtype=np.float64) method_errors[method].extend(errors.tolist()) method_validation[method].extend(validation.tolist()) by_cfo_rows.append({ "symbol_rate": symbol_rate, "quality": quality_name, "target_evm_rms": evm_rms, "true_frequency_hz": float(true_frequency), "method": method, "trials": MONTE_CARLO_TRIALS, **_error_statistics(errors), "validation_phase_rmse_rad": float(np.mean(validation)), }) for method in METHODS: errors = np.asarray(method_errors[method], dtype=np.float64) validation = np.asarray(method_validation[method], dtype=np.float64) distributions[(symbol_rate, quality_name, method)] = errors aggregate_rows.append({ "symbol_rate": symbol_rate, "quality": quality_name, "target_evm_rms": evm_rms, "method": method, "trials": len(errors), "cfo_values": len(SYNTHETIC_FREQUENCIES_HZ), **_error_statistics(errors), "validation_phase_rmse_rad": float(np.mean(validation)), }) return aggregate_rows, by_cfo_rows, distributions, clean_metadata def _error_statistics(errors: np.ndarray) -> dict: absolute = np.abs(errors) return { "bias_hz": float(np.mean(errors)), "rmse_hz": float(np.sqrt(np.mean(errors**2))), "median_abs_error_hz": float(np.median(absolute)), "p95_abs_error_hz": float(np.percentile(absolute, 95)), "p99_abs_error_hz": float(np.percentile(absolute, 99)), "catastrophic_error_rate": float(np.mean(absolute > CATASTROPHIC_ERROR_HZ)), "probability_abs_error_le_0p25_hz": float(np.mean(absolute <= 0.25)), "probability_abs_error_le_0p5_hz": float(np.mean(absolute <= 0.5)), "probability_abs_error_le_1p0_hz": float(np.mean(absolute <= 1.0)), } def _csv_value(value): if isinstance(value, complex): return json.dumps(serializable(value), ensure_ascii=False) if isinstance(value, (dict, list, tuple)): return json.dumps(serializable(value), ensure_ascii=False, separators=(",", ":")) return value def write_csv(path: Path, rows: list[dict], exclude: tuple[str, ...] = ()) -> None: fields: list[str] = [] for row in rows: fields.extend(key for key in row if key not in fields and key not in exclude) with path.open("w", newline="", encoding="utf-8-sig") as stream: writer = csv.DictWriter(stream, fieldnames=fields) writer.writeheader() for row in rows: writer.writerow({key: _csv_value(row.get(key)) for key in fields}) def summarize_real(rows: list[dict]) -> dict: summary = {} for dataset in ("Lab056", "Lab054", "Lab047"): summary[dataset] = {} for method in METHODS: selected = [row for row in rows if row["dataset"] == dataset and row["method"] == method] if not selected: continue if dataset == "Lab056": b_rows = [row for row in selected if row["label"].startswith("B")] a_rows = [row for row in selected if row["label"].startswith("A")] def mean_or_none(items: list[dict], key: str) -> float | None: return float(np.mean([row[key] for row in items])) if items else None summary[dataset][method] = { "segments": len(selected), "mean_abs_frequency_error_hz": float(np.mean([abs(row["frequency_error_hz"]) for row in selected])), "A_frequency_errors_hz": [row["frequency_error_hz"] for row in a_rows], "B_frequency_errors_hz": [row["frequency_error_hz"] for row in b_rows], "B_mean_abs_frequency_error_hz": ( float(np.mean([abs(row["frequency_error_hz"]) for row in b_rows])) if b_rows else None ), "B_mean_symbol_error_rate": mean_or_none(b_rows, "symbol_error_rate"), "A_mean_symbol_error_rate": mean_or_none(a_rows, "symbol_error_rate"), "B_mean_evm_rms": mean_or_none(b_rows, "evm_rms"), "B_mean_matched_residual_db": mean_or_none(b_rows, "matched_residual_db"), } else: summary[dataset][method] = { "frames": len(selected), "header_accepted": sum(bool(row["HEADER_ACCEPTED"]) for row in selected), "fec_started": sum(bool(row["FEC_STARTED"]) for row in selected), "crc_pass": sum(bool(row["CRC_PASS"]) for row in selected), "packet_output": sum(bool(row["PACKET_OUTPUT"]) for row in selected), "preamble_errors": sum(int(row["preamble_errors"]) for row in selected), "header_raw_errors": sum(int(row["header_raw_errors"]) for row in selected), "header_voted_errors": sum(int(row["header_voted_errors"]) for row in selected), "mean_frequency_hz": float(np.mean([row["frequency_hz"] for row in selected])), } return summary def phase_limit_table() -> list[dict]: return [ { "phase_limit_degrees": limit, "max_frequency_error_hz": maximum_frequency_error_hz(limit, BURST_DURATION_SECONDS), } for limit in (15, 30, 45, 60, 90) ] def mathematics_audit() -> dict: return { "active_input": "matched-filter symbol centers after nominal 50 kHz downconversion", "p0_squares_bpsk": True, "p0_symbol_ranges_half_open": [[0, 512], [512, 1024]], "p0_discarded_preamble_symbols_half_open": [1024, 1280], "block_center_symbol_times": [256.0, 768.0], "center_distance_symbols": 512, "phase_difference_formula": "angle(block1)-angle(block0), unwrapped", "frequency_formula": "delta_phase/(4*pi*(512/Rs))", "factor_two_after_squaring": True, "compensation_sign": "negative complex exponential", "reference_time": "gain phase at t=0; symbol k center at (k+0.5)/Rs", "p0_unambiguous_abs_hz": {"20000": 20_000 / 2048, "500000": 500_000 / 2048}, "double_compensation": False, } def choose_verdict(real_summary: dict, synthetic_rows: list[dict]) -> tuple[str, str, str]: b_quality = { row["method"]: row for row in synthetic_rows if row["symbol_rate"] == 500_000 and row["quality"] == "LAB058_B_EVM" } stable = { method: b_quality[method]["probability_abs_error_le_0p5_hz"] >= 0.95 for method in METHODS } best = min(METHODS, key=lambda method: real_summary["Lab056"][method]["B_mean_symbol_error_rate"]) diagnostic = 0.03495896 approaches_diagnostic = real_summary["Lab056"][best]["B_mean_symbol_error_rate"] <= diagnostic + 0.02 lab047_preserved = real_summary["Lab047"][best]["packet_output"] == 36 a_preserved = real_summary["Lab056"][best]["A_mean_symbol_error_rate"] <= 0.001 if stable[best] and approaches_diagnostic and lab047_preserved and a_preserved: return ( "PREAMBLE_ONLY_SUFFICIENT", f"{best} стабилен в Monte Carlo и сохраняет A/Lab047", f"подготовить внедрение {best} в штатный receiver", ) if max(row["probability_abs_error_le_0p5_hz"] for row in b_quality.values()) >= 0.5: return ( "PREAMBLE_ONLY_MARGINAL", "улучшение есть, но субгерцовая точность на качестве B достигается нестабильно", "проверить decision-directed carrier tracking through burst", ) return ( "PREAMBLE_ONLY_INSUFFICIENT", "ни один preamble-only метод не обеспечивает устойчивую ошибку менее 0,5 Гц на качестве B", "проверить decision-directed carrier tracking through burst", ) def make_plots( real_rows: list[dict], synthetic_rows: list[dict], distributions: dict, ) -> None: colors = {"P0": "#6c757d", "P1": "#2878b5", "P2": "#e17c05", "P3": "#54a24b"} b_rows = [row for row in real_rows if row["dataset"] == "Lab056" and row["label"].startswith("B")] figure, axis = plt.subplots(figsize=(9, 5), constrained_layout=True) x = np.arange(2) width = 0.18 for index, method in enumerate(METHODS): selected = sorted((row for row in b_rows if row["method"] == method), key=lambda row: row["label"]) axis.bar(x + (index - 1.5) * width, [row["frequency_error_hz"] for row in selected], width, label=method, color=colors[method]) axis.axhline(0, color="black", linewidth=0.8) axis.set(xticks=x, xticklabels=["B1", "B2"], ylabel="Ошибка оценки, Гц", title="Lab056 B: preamble-only CFO относительно whole-burst reference") axis.legend(ncol=4); axis.grid(axis="y", alpha=0.25) figure.savefig(OUTPUT / "frequency_error_by_method.png", dpi=160) plt.close(figure) figure, axis = plt.subplots(figsize=(9, 5), constrained_layout=True) bins = np.linspace(-35, 35, 141) for method in METHODS: values = distributions[(500_000, "LAB058_B_EVM", method)] axis.hist(values, bins=bins, density=True, histtype="step", linewidth=1.5, label=method, color=colors[method]) axis.axvspan(-0.5, 0.5, color="green", alpha=0.12, label="±0,5 Гц") axis.set(xlabel="Ошибка оценки частоты, Гц", ylabel="Плотность", title="Synthetic 500 ксимв/с, качество Lab058 B") axis.legend(); axis.grid(alpha=0.2) figure.savefig(OUTPUT / "synthetic_500k_error_distribution.png", dpi=160) plt.close(figure) figure, axis = plt.subplots(figsize=(9, 5), constrained_layout=True) time = np.linspace(0, BURST_DURATION_SECONDS, 300) for method in METHODS: errors = [row["frequency_error_hz"] for row in b_rows if row["method"] == method] axis.plot(time * 1000, 360 * float(np.mean(errors)) * time, label=method, color=colors[method]) for limit in (-30, 30): axis.axhline(limit, color="black", linestyle="--", linewidth=0.8) axis.set(xlabel="Время от начала burst, мс", ylabel="Фазовый уход, градусы", title="Средняя ошибка B1/B2, экстраполированная на burst") axis.legend(); axis.grid(alpha=0.25) figure.savefig(OUTPUT / "burst_phase_drift.png", dpi=160) plt.close(figure) figure, axis = plt.subplots(figsize=(9, 5), constrained_layout=True) for index, method in enumerate(METHODS): selected = sorted((row for row in b_rows if row["method"] == method), key=lambda row: row["label"]) axis.bar(x + (index - 1.5) * width, [100 * row["symbol_error_rate"] for row in selected], width, label=method, color=colors[method]) axis.axhline(3.495896, color="black", linestyle="--", label="whole-burst reference") axis.set(xticks=x, xticklabels=["B1", "B2"], ylabel="Ошибочные символы, %", title="Lab056 B: полный burst после preamble-only оценки") axis.legend(ncol=3); axis.grid(axis="y", alpha=0.25) figure.savefig(OUTPUT / "lab056_b_symbol_errors.png", dpi=160) plt.close(figure) lab054_rows = [row for row in real_rows if row["dataset"] == "Lab054"] figure, axis = plt.subplots(figsize=(9, 5), constrained_layout=True) accepted = [sum(row["HEADER_ACCEPTED"] for row in lab054_rows if row["method"] == method) for method in METHODS] crc = [sum(row["CRC_PASS"] for row in lab054_rows if row["method"] == method) for method in METHODS] positions = np.arange(len(METHODS)) axis.bar(positions - 0.18, accepted, 0.36, label="HEADER_ACCEPTED") axis.bar(positions + 0.18, crc, 0.36, label="CRC_PASS") axis.set(xticks=positions, xticklabels=METHODS, ylabel="Кадры из 16", title="Lab054: строгий результат по preamble-only методам", ylim=(0, 16.5)) if max(accepted + crc) == 0: axis.text(1.5, 8.0, "Все методы: 0/16 HEADER и CRC", ha="center", va="center") axis.set_xlim(-0.6, 3.6) axis.legend(); axis.grid(axis="y", alpha=0.25) figure.savefig(OUTPUT / "lab054_result_by_method.png", dpi=160) plt.close(figure) def run_functional_tests() -> None: """Быстрые проверки формул без файлов IQ.""" assert maximum_frequency_error_hz(30, BURST_DURATION_SECONDS) == 0.40690104166666663 assert np.isclose(phase_slope_to_hz(2 * np.pi * 7, squared_bpsk=False), 7) assert np.isclose(phase_slope_to_hz(4 * np.pi * -7, squared_bpsk=True), -7) assert set(P2_TRAIN).isdisjoint(P2_VALIDATION) assert len(P2_TRAIN) == len(P2_VALIDATION) == 576 def main() -> None: OUTPUT.mkdir(parents=True, exist_ok=True) run_functional_tests() source_paths = [ ROOT / "data/raw/lab056/rx_stream.cs8", ROOT / "data/raw/lab056/tx_abab.cs8", ROOT / "data/raw/lab054/rx_stream.cs8", ROOT / "data/raw/lab054/tx_stream.cs8", ] audit = json.loads((ROOT / "data/processed/lab055/energy_audit_v2.json").read_text(encoding="utf-8")) lab047_directory = Path(audit["experiments"]["Lab047"]["directory"]) source_paths.append(lab047_directory / "rx_stream.cs8") integrity_before = {str(path): {"bytes": path.stat().st_size, "sha256": sha256(path)} for path in source_paths} # P0 действительно завершается и сохраняется до вычисления P1/P2/P3. baseline_rows = analyze_lab056(("P0",)) + analyze_packets(("P0",)) baseline = { "algorithm": "unchanged active CURRENT estimator", "rows": baseline_rows, "summary": summarize_real(baseline_rows), "mathematics": mathematics_audit(), } write_json(OUTPUT / "baseline_p0.json", baseline) real_rows = analyze_lab056(METHODS) + analyze_packets(METHODS) real_summary = summarize_real(real_rows) synthetic_rows, by_cfo_rows, distributions, clean_metadata = synthetic_study() verdict, reason, next_step = choose_verdict(real_summary, synthetic_rows) integrity_after = {str(path): {"bytes": path.stat().st_size, "sha256": sha256(path)} for path in source_paths} if integrity_before != integrity_after: raise AssertionError("Сохранённые исходные IQ изменились") phase_limits = phase_limit_table() summary = { "status": "COMPLETED", "lab058_commit": "82f17499ace5e655991240348401b9700c7c3349", "rf_tx": 0, "hardware_access": 0, "preamble_symbols": PREAMBLE_SYMBOLS, "preamble_duration_seconds": {"20000": PREAMBLE_SYMBOLS / 20_000, "500000": PREAMBLE_SYMBOLS / 500_000}, "burst_duration_seconds": BURST_DURATION_SECONDS, "phase_limits": phase_limits, "catastrophic_error_threshold_hz": CATASTROPHIC_ERROR_HZ, "mathematics": mathematics_audit(), "p2_split": {"train_half_open": [32, 608], "validation_half_open": [672, 1248]}, "real": real_summary, "synthetic": { "trials_per_frequency_quality_rate": MONTE_CARLO_TRIALS, "frequency_values_hz": SYNTHETIC_FREQUENCIES_HZ, "quality_levels": [{"name": name, "evm_rms": value} for name, value in QUALITY_LEVELS], "clean_rrc_controls": clean_metadata, "awgn_crlb_standard_deviation_hz": { "formula": "sqrt(6*EVM^2*Rs^2/((2*pi)^2*N*(N^2-1)))", "assumptions": "known symbols, unknown constant phase, complex AWGN", "500000_lab058_b_evm": awgn_cfo_crlb_standard_deviation_hz( 0.845350, PREAMBLE_SYMBOLS, 500_000 ), }, "aggregate": synthetic_rows, }, "verdict": verdict, "main_reason": reason, "recommended_next_single_step": next_step, "future_options": { "decision_directed": "первый кандидат: не меняет формат и не знает payload truth", "costas_loop": "альтернатива для слежения по всему burst; требует настройки полосы петли", "periodic_pilots": "наиболее явная опора, но требует изменения формата и overhead", }, "physical_limit": { "lab058_b_diagnostic_symbol_error_rate": 0.03495896, "lab058_b_evm_rms": 0.845350, "lab058_b_matched_residual_db": 1.464846, "lab058_b_channel_ripple_db": 5.471382, }, "lab054_historical_result": "NOT_A, 0/16", "lab059_commit": "NOT_PERFORMED", "push": "NOT_PERFORMED", "integrity": {"before": integrity_before, "after": integrity_after, "passed": True}, } write_csv(OUTPUT / "real_results.csv", real_rows, exclude=("estimate",)) write_csv(OUTPUT / "synthetic_results.csv", synthetic_rows) write_csv(OUTPUT / "synthetic_by_cfo.csv", by_cfo_rows) write_json(OUTPUT / "summary.json", summary) make_plots(real_rows, synthetic_rows, distributions) print(json.dumps(serializable({ "status": summary["status"], "verdict": verdict, "main_reason": reason, "real": real_summary, "synthetic_500k_b": [ row for row in synthetic_rows if row["symbol_rate"] == 500_000 and row["quality"] == "LAB058_B_EVM" ], "rf_tx": 0, }), ensure_ascii=False, indent=2)) if __name__ == "__main__": main()