"""Офлайн-диагностика асимметрии ошибок BPSK Lab043. Скрипт читает только сохранённые NPY/JSON, проверяет SHA-256 и не импортирует аппаратные библиотеки. Результаты предназначены для диагностики, а не для автоматического изменения алгоритма приёмника. """ from __future__ import annotations import argparse import csv import hashlib import json import math from dataclasses import dataclass from pathlib import Path import numpy as np from matplotlib import pyplot as plt from scipy.signal import csd, fftconvolve, welch from experiments import lab043_pluto_to_rtlsdr as lab043 from protocol import bpsk_radio as radio VALID_CAPTURE_NAMES = ( "lab043_prbs11_s_20260819_134831.json", "lab043_prbs11_s_20260819_150001.json", ) TRIPLETS = tuple(f"{value:03b}" for value in range(8)) @dataclass(frozen=True) class CaptureAnalysis: name: str metadata: dict raw_bpsk: np.ndarray coarse_corrected: np.ndarray matched: np.ndarray payload_bits: np.ndarray payload_symbols: np.ndarray payload_normalized: np.ndarray sample_phase: int start_symbol: int marker_score: float residual_cfo_hz: float sha256: str def sha256_file(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 _known_sequence_phase_fit( received: np.ndarray, expected: np.ndarray, symbol_rate: float, block_symbols: int = 64, ) -> tuple[np.ndarray, float, float]: """Убрать только общую линейную фазу по известной PRBS для диагностики.""" usable = len(received) // block_symbols * block_symbols despread = (received[:usable] * np.conj(expected[:usable])).reshape( -1, block_symbols ) block_sums = np.sum(despread, axis=1) centers = ( np.arange(len(block_sums), dtype=np.float64) * block_symbols + (block_symbols - 1) / 2.0 ) phases = np.unwrap(np.angle(block_sums)) weights = np.maximum(np.abs(block_sums), 1e-12) design = np.column_stack((np.ones(len(centers)), centers)) coefficients = np.zeros(2, dtype=np.float64) for _ in range(12): root = np.sqrt(weights) coefficients, *_ = np.linalg.lstsq( design * root[:, None], phases * root, rcond=None ) residual = phases - design @ coefficients center = float(np.median(residual)) scale = 1.4826 * float(np.median(np.abs(residual - center))) + 1e-12 limit = 1.345 * scale robust = np.ones_like(residual) outlier = np.abs(residual) > limit robust[outlier] = limit / np.abs(residual[outlier]) weights = np.maximum(np.abs(block_sums), 1e-12) * robust indexes = np.arange(len(received), dtype=np.float64) corrected = received * np.exp(-1j * (coefficients[0] + coefficients[1] * indexes)) residual_cfo_hz = float(coefficients[1] * symbol_rate / (2.0 * np.pi)) fit_rmse = float(np.sqrt(np.mean((phases - design @ coefficients) ** 2))) return corrected, residual_cfo_hz, fit_rmse def _normalize_symbols(received: np.ndarray, expected: np.ndarray) -> tuple[np.ndarray, complex]: gain = np.vdot(expected, received) / (np.vdot(expected, expected) + 1e-15) return received / gain, complex(gain) def load_capture( raw_dir: Path, metadata_name: str, *, require_processing_success: bool = True, ) -> CaptureAnalysis: metadata_path = raw_dir / metadata_name metadata = json.loads(metadata_path.read_text(encoding="utf-8")) iq_path = raw_dir / metadata["iq_file"] actual_sha = sha256_file(iq_path) if actual_sha != metadata["iq_sha256"]: raise RuntimeError(f"SHA-256 не совпал: {iq_path}") if metadata.get("capture_status") != "captured": raise RuntimeError(f"Захват не имеет статуса captured: {metadata_path}") analysis = metadata.get("analysis", {}) if metadata.get("processing_status") != "success" or ( require_processing_success and not analysis.get("processing_success", False) ): raise RuntimeError(f"Обработка захвата неуспешна: {metadata_path}") continuity = analysis.get("continuity", {}) if continuity.get("confirmed_discontinuity_count") != 0: raise RuntimeError(f"В захвате есть подтверждённые разрывы: {metadata_path}") raw = np.load(iq_path, mmap_mode="r") sections = analysis["sections"] raw_bpsk = np.asarray( raw[sections["bpsk_slice_start_sample"] : sections["bpsk_slice_end_sample"]], dtype=np.complex128, ) fs = float(metadata["actual_sample_rate_hz"]) coarse = radio.apply_coarse_frequency_correction( raw_bpsk, analysis["calibration"]["carrier_offset_hz"], fs ) taps = radio.root_raised_cosine_taps( lab043.LAB042_RRC_ROLLOFF, lab043.SAMPLES_PER_SYMBOL, lab043.LAB042_RRC_SPAN_SYMBOLS, ) matched = fftconvolve(coarse, taps, mode="full") marker_bits, marker_symbols = radio.build_frame_marker() found = radio.find_radio_frame(matched, marker_symbols, lab043.SAMPLES_PER_SYMBOL) pilot_count = int(metadata["waveform"].get("pilot_symbol_count", 0)) payload_count = int(metadata["waveform"]["transmitted_prbs_bit_count"]) required = len(marker_bits) + pilot_count + payload_count start = int(found["start_symbol_index"]) symbols = found["symbol_samples"][start : start + required] if len(symbols) != required: raise RuntimeError(f"Кадр обрезан в {metadata_path}") marker_carrier = radio.estimate_carrier_parameters( symbols[: len(marker_bits)], marker_symbols ) symbols = radio.correct_phase_and_frequency( symbols, float(marker_carrier["constant_phase"]), 0.0 ) payload_bits = lab043.prbs11(payload_count, lab043.PRBS11_INITIAL_STATE) expected = radio.bpsk_modulate(payload_bits) payload = symbols[len(marker_bits) + pilot_count :] phase_corrected, residual_cfo_hz, _ = _known_sequence_phase_fit( payload, expected, lab043.SYMBOL_RATE ) normalized, _ = _normalize_symbols(phase_corrected, expected) return CaptureAnalysis( name=metadata_path.stem, metadata=metadata, raw_bpsk=raw_bpsk, coarse_corrected=coarse, matched=matched, payload_bits=payload_bits, payload_symbols=payload, payload_normalized=normalized, sample_phase=int(found["sample_phase"]), start_symbol=start, marker_score=float(found["score"]), residual_cfo_hz=residual_cfo_hz, sha256=actual_sha, ) def triplet_rows( source: str, bits: np.ndarray, values: np.ndarray, *, raw_values: np.ndarray | None = None, fitted_gain: complex = 1.0 + 0.0j, ) -> list[dict]: decisions = (values.real >= 0.0).astype(np.uint8) indexes = np.arange(1, len(bits) - 1) rows: list[dict] = [] for triplet in TRIPLETS: pattern = np.fromiter((int(bit) for bit in triplet), dtype=np.uint8) selected = indexes[ (bits[indexes - 1] == pattern[0]) & (bits[indexes] == pattern[1]) & (bits[indexes + 1] == pattern[2]) ] group = values[selected] raw_group = group if raw_values is None else raw_values[selected] ideal = 1.0 if pattern[1] else -1.0 errors = int(np.count_nonzero(decisions[selected] != bits[selected])) if len(group): phase_errors = np.angle(group * ideal) circular_phase = float(np.angle(np.mean(np.exp(1j * phase_errors)))) metrics = { "ber": errors / len(selected), "mean_i": float(np.mean(group.real)), "mean_q": float(np.mean(group.imag)), "median_i": float(np.median(group.real)), "median_q": float(np.median(group.imag)), "std_i": float(np.std(group.real)), "std_q": float(np.std(group.imag)), "mean_amplitude": float(np.mean(np.abs(group))), "mean_phase_error_deg": math.degrees(circular_phase), "mean_distance_to_ideal": float(np.mean(np.abs(group - ideal))), "decision_cross_fraction": errors / len(selected), "raw_mean_i": float(np.mean(raw_group.real)), "raw_mean_q": float(np.mean(raw_group.imag)), } else: metrics = { key: float("nan") for key in ( "ber", "mean_i", "mean_q", "median_i", "median_q", "std_i", "std_q", "mean_amplitude", "mean_phase_error_deg", "mean_distance_to_ideal", "decision_cross_fraction", "raw_mean_i", "raw_mean_q", ) } rows.append( { "source": source, "triplet": triplet, "count": int(len(selected)), "bit_errors": errors, **metrics, "fitted_gain_i": float(fitted_gain.real), "fitted_gain_q": float(fitted_gain.imag), } ) return rows def tx_stage_rows(bits: np.ndarray) -> list[dict]: marker_bits, marker_symbols = radio.build_frame_marker() pilot_bits = lab043.pilot_prbs7() all_symbols = np.concatenate( (marker_symbols, radio.bpsk_modulate(pilot_bits), radio.bpsk_modulate(bits)) ) expected = radio.bpsk_modulate(bits) taps = radio.root_raised_cosine_taps( lab043.LAB042_RRC_ROLLOFF, lab043.SAMPLES_PER_SYMBOL, lab043.LAB042_RRC_SPAN_SYMBOLS, ) upsampled = np.zeros(len(all_symbols) * lab043.SAMPLES_PER_SYMBOL, dtype=np.complex128) upsampled[:: lab043.SAMPLES_PER_SYMBOL] = all_symbols guard = np.zeros( lab043.LAB042_GUARD_SYMBOL_COUNT * lab043.SAMPLES_PER_SYMBOL, dtype=np.complex128, ) rrc = fftconvolve(np.concatenate((guard, upsampled, guard)), taps, mode="full") payload_symbol_offset = len(marker_bits) + len(pilot_bits) first = ( len(guard) + (len(taps) - 1) // 2 + payload_symbol_offset * lab043.SAMPLES_PER_SYMBOL ) positions = first + np.arange(len(bits)) * lab043.SAMPLES_PER_SYMBOL rrc_symbols = rrc[positions] scaled = rrc_symbols * lab043.LAB042_TX_AMPLITUDE pluto_complex64 = (scaled * (2**14)).astype(np.complex64) stages = ( ("идеальные символы", expected), ("после TX RRC", rrc_symbols), ("после масштаба 0.70", scaled), ("после complex64 и масштаба 2^14", pluto_complex64), ) rows: list[dict] = [] for name, raw_values in stages: normalized, gain = _normalize_symbols(raw_values, expected) rows.extend( triplet_rows(name, bits, normalized, raw_values=raw_values, fitted_gain=gain) ) return rows def _fractional_delay(values: np.ndarray, delay_samples: float) -> np.ndarray: indexes = np.arange(len(values), dtype=np.float64) source = indexes - delay_samples real = np.interp(source, indexes, values.real, left=0.0, right=0.0) imag = np.interp(source, indexes, values.imag, left=0.0, right=0.0) return real + 1j * imag def ideal_channel_rows(bits: np.ndarray) -> list[dict]: marker_bits, marker_symbols = radio.build_frame_marker() pilot_bits = lab043.pilot_prbs7() pilot_symbols = radio.bpsk_modulate(pilot_bits) expected = radio.bpsk_modulate(bits) all_symbols = np.concatenate((marker_symbols, pilot_symbols, expected)) waveform = lab043.shape_frame_like_lab042(all_symbols) taps = radio.root_raised_cosine_taps( lab043.LAB042_RRC_ROLLOFF, lab043.SAMPLES_PER_SYMBOL, lab043.LAB042_RRC_SPAN_SYMBOLS, ) scale = 0.9999965 variants: list[tuple[str, np.ndarray]] = [("без искажений", waveform)] variants.append(("только SRO -3.5 ppm", lab043.radio.resample_for_clock_scale(waveform, scale))) indexes = np.arange(len(waveform), dtype=np.float64) cfo = waveform * np.exp(1j * 2.0 * np.pi * 0.5 * indexes / lab043.SAMPLE_RATE_HZ) variants.append(("только остаточная CFO +0.5 Гц", cfo)) variants.append(("только дробный тайминг +0.45 отсчёта", _fractional_delay(waveform, 0.45))) variants.append(("только gain 0.8 exp(j0.4)", waveform * 0.8 * np.exp(0.4j))) combined = lab043.radio.resample_for_clock_scale(waveform, scale) combined = _fractional_delay(combined, 0.45) combined_indexes = np.arange(len(combined), dtype=np.float64) combined *= 0.8 * np.exp( 0.4j + 1j * 2.0 * np.pi * 0.5 * combined_indexes / lab043.SAMPLE_RATE_HZ ) variants.append(("SRO+CFO+тайминг+gain", combined)) rows: list[dict] = [] for name, channel in variants: matched = fftconvolve(channel, taps, mode="full") found = radio.find_radio_frame(matched, marker_symbols, lab043.SAMPLES_PER_SYMBOL) start = int(found["start_symbol_index"]) required = len(all_symbols) symbols = found["symbol_samples"][start : start + required] marker_carrier = radio.estimate_carrier_parameters( symbols[: len(marker_bits)], marker_symbols ) symbols = radio.correct_phase_and_frequency( symbols, float(marker_carrier["constant_phase"]), 0.0 ) payload = symbols[len(marker_bits) + len(pilot_bits) :] normalized, gain = _normalize_symbols(payload, expected) variant_rows = triplet_rows( name, bits, normalized, raw_values=payload, fitted_gain=gain ) for row in variant_rows: row["sample_phase"] = int(found["sample_phase"]) row["marker_score"] = float(found["score"]) rows.extend(variant_rows) return rows def rrc_rows() -> tuple[list[dict], dict]: taps = radio.root_raised_cosine_taps( lab043.LAB042_RRC_ROLLOFF, lab043.SAMPLES_PER_SYMBOL, lab043.LAB042_RRC_SPAN_SYMBOLS, ) full = np.convolve(taps, taps) center = len(full) // 2 normalized = full / full[center] rows = [] for offset in range(-5, 6): rows.append( { "symbol_offset": offset, "sample_offset": offset * lab043.SAMPLES_PER_SYMBOL, "h_normalized": float(normalized[center + offset * lab043.SAMPLES_PER_SYMBOL]), } ) audit = { "tx_tap_count": len(taps), "rx_tap_count": len(taps), "rolloff": lab043.LAB042_RRC_ROLLOFF, "span_symbols": lab043.LAB042_RRC_SPAN_SYMBOLS, "samples_per_symbol": lab043.SAMPLES_PER_SYMBOL, "symmetry_max_abs_error": float(np.max(np.abs(taps - taps[::-1]))), "tap_energy": float(np.sum(taps**2)), "tx_group_delay_samples": (len(taps) - 1) // 2, "rx_group_delay_samples": (len(taps) - 1) // 2, "total_group_delay_samples": len(taps) - 1, "total_group_delay_symbols": (len(taps) - 1) / lab043.SAMPLES_PER_SYMBOL, "sampling_index_rule": "find_radio_frame: phase 0..119, затем одинаковый индекс phase+n*120 для обоих знаков", "one_sample_sign_dependent_shift": False, } return rows, audit def phase_sweep_rows(capture: CaptureAnalysis) -> list[dict]: marker_bits, marker_symbols = radio.build_frame_marker() pilot_count = int(capture.metadata["waveform"].get("pilot_symbol_count", 0)) expected = radio.bpsk_modulate(capture.payload_bits) rows: list[dict] = [] for phase in range(lab043.SAMPLES_PER_SYMBOL): symbol_samples = capture.matched[phase:: lab043.SAMPLES_PER_SYMBOL] correlation = np.correlate(symbol_samples, marker_symbols, mode="valid") energy = np.convolve( np.abs(symbol_samples) ** 2, np.ones(len(marker_symbols)), mode="valid" ) score = np.abs(correlation) / np.sqrt(energy * len(marker_symbols) + 1e-15) start = int(np.argmax(score)) required = len(marker_bits) + pilot_count + len(expected) symbols = symbol_samples[start : start + required] if len(symbols) != required: continue marker_carrier = radio.estimate_carrier_parameters( symbols[: len(marker_bits)], marker_symbols ) symbols = radio.correct_phase_and_frequency( symbols, float(marker_carrier["constant_phase"]), 0.0 ) payload = symbols[len(marker_bits) + pilot_count :] normalized, _ = _normalize_symbols(payload, expected) decisions = (normalized.real >= 0.0).astype(np.uint8) errors = decisions != capture.payload_bits indexes = np.arange(1, len(expected) - 1) triplet_labels = np.asarray( [ f"{capture.payload_bits[i-1]}{capture.payload_bits[i]}{capture.payload_bits[i+1]}" for i in indexes ] ) error_indexes = indexes[errors[indexes]] error_labels = triplet_labels[errors[indexes]] negative = normalized[capture.payload_bits == 0].real positive = normalized[capture.payload_bits == 1].real rows.append( { "capture": capture.name, "phase": phase, "marker_start_symbol": start, "marker_score": float(score[start]), "bit_errors": int(np.count_nonzero(errors)), "ber": float(np.mean(errors)), "errors_000": int(np.count_nonzero(error_labels == "000")), "errors_001": int(np.count_nonzero(error_labels == "001")), "errors_other": int( len(error_indexes) - np.count_nonzero(error_labels == "000") - np.count_nonzero(error_labels == "001") ), "negative_mean_i": float(np.mean(negative)), "positive_mean_i": float(np.mean(positive)), "negative_p95_i": float(np.percentile(negative, 95)), "positive_p05_i": float(np.percentile(positive, 5)), "eye_opening_i_p05_p95": float( np.percentile(positive, 5) - np.percentile(negative, 95) ), "negative_std_i": float(np.std(negative)), "positive_std_i": float(np.std(positive)), } ) return rows def dc_iq_rows(capture: CaptureAnalysis) -> list[dict]: expected = radio.bpsk_modulate(capture.payload_bits) values = capture.payload_normalized plus = values[capture.payload_bits == 1] minus = values[capture.payload_bits == 0] plus_center = complex(np.mean(plus)) minus_center = complex(np.mean(minus)) dc = (plus_center + minus_center) / 2.0 residual = np.concatenate((plus - plus_center, minus - minus_center)) covariance = np.cov(np.column_stack((residual.real, residual.imag)), rowvar=False) eigenvalues = np.linalg.eigvalsh(covariance) ellipse_ratio = math.sqrt(float(eigenvalues[-1] / max(eigenvalues[0], 1e-15))) phase_separation = abs( math.degrees(np.angle(np.exp(1j * (np.angle(plus_center) - np.angle(minus_center))))) ) return [ { "capture": capture.name, "raw_bpsk_mean_i": float(np.mean(capture.raw_bpsk).real), "raw_bpsk_mean_q": float(np.mean(capture.raw_bpsk).imag), "coarse_before_rrc_mean_i": float(np.mean(capture.coarse_corrected).real), "coarse_before_rrc_mean_q": float(np.mean(capture.coarse_corrected).imag), "matched_mean_i": float(np.mean(capture.matched).real), "matched_mean_q": float(np.mean(capture.matched).imag), "positive_center_i": plus_center.real, "positive_center_q": plus_center.imag, "negative_center_i": minus_center.real, "negative_center_q": minus_center.imag, "symbol_dc_i": dc.real, "symbol_dc_q": dc.imag, "positive_distance_from_zero": abs(plus_center), "negative_distance_from_zero": abs(minus_center), "positive_decision_margin_i": plus_center.real, "negative_decision_margin_i": -minus_center.real, "residual_std_i": float(np.std(residual.real)), "residual_std_q": float(np.std(residual.imag)), "iq_std_ratio": float(np.std(residual.real) / max(np.std(residual.imag), 1e-15)), "iq_correlation": float(np.corrcoef(residual.real, residual.imag)[0, 1]), "ellipse_axis_ratio": ellipse_ratio, "cluster_phase_separation_deg": phase_separation, "negative_is_closer_to_boundary": bool(-minus_center.real < plus_center.real), } ] def run_length_rows(capture: CaptureAnalysis) -> list[dict]: bits = capture.payload_bits values = capture.payload_normalized decisions = (values.real >= 0.0).astype(np.uint8) rows: list[dict] = [] start = 0 runs: list[tuple[int, int, int]] = [] for index in range(1, len(bits) + 1): if index == len(bits) or bits[index] != bits[start]: runs.append((start, index, int(bits[start]))) start = index grouped: dict[tuple[int, str, str], list[int]] = {} for first, end, bit in runs: length = end - first category = str(length) if length <= 6 else "7+" for index in range(first, end): if length == 1: position = "single" elif index == first: position = "first" elif index == end - 1: position = "last" else: position = "middle" grouped.setdefault((bit, category, position), []).append(index) for (bit, category, position), selected_list in sorted(grouped.items()): selected = np.asarray(selected_list, dtype=np.int64) group = values[selected] ideal = 1.0 if bit else -1.0 phase = np.angle(group * ideal) rows.append( { "capture": capture.name, "bit": bit, "run_length": category, "position": position, "symbol_count": len(selected), "mean_i": float(np.mean(group.real)), "mean_q": float(np.mean(group.imag)), "mean_amplitude": float(np.mean(np.abs(group))), "mean_signed_amplitude": float(np.mean(group.real * ideal)), "mean_phase_error_deg": math.degrees( float(np.angle(np.mean(np.exp(1j * phase)))) ), "ber": float(np.mean(decisions[selected] != bits[selected])), } ) return rows def spectrum_rows(captures: list[CaptureAnalysis], bits: np.ndarray) -> list[dict]: marker_bits, marker_symbols = radio.build_frame_marker() pilot_bits = lab043.pilot_prbs7() tx_symbols = np.concatenate( (marker_symbols, radio.bpsk_modulate(pilot_bits), radio.bpsk_modulate(bits)) ) tx = lab043.shape_frame_like_lab042(tx_symbols) sources: list[tuple[str, np.ndarray, float]] = [("чистый TX waveform", tx, lab043.SAMPLE_RATE_HZ)] sources.extend( (f"{capture.name}: RX после грубой CFO", capture.coarse_corrected, float(capture.metadata["actual_sample_rate_hz"])) for capture in captures ) rows: list[dict] = [] for name, values, fs in sources: frequencies, psd = welch( values, fs=fs, window="hann", nperseg=min(65536, len(values)), return_onesided=False, scaling="density", ) frequencies = np.fft.fftshift(frequencies) psd = np.fft.fftshift(psd) total = float(np.sum(psd)) + 1e-30 near_reference = (np.abs(frequencies) >= 2_000) & (np.abs(frequencies) <= 10_000) dc_index = int(np.argmin(np.abs(frequencies))) row = { "source": name, "mean_i": float(np.mean(values).real), "mean_q": float(np.mean(values).imag), "mean_amplitude": float(np.mean(np.abs(values))), "dc_bin_to_2_10khz_median_db": float( 10.0 * np.log10((psd[dc_index] + 1e-30) / (np.median(psd[near_reference]) + 1e-30)) ), } for band in (100, 500, 1_000, 5_000, 20_000): row[f"power_fraction_abs_f_below_{band}_hz"] = float( np.sum(psd[np.abs(frequencies) < band]) / total ) rows.append(row) return rows def symbol_frequency_rows(captures: list[CaptureAnalysis]) -> list[dict]: """Оценить передаточную характеристику по известным символам без коррекции АЧХ.""" rows: list[dict] = [] bands = ((10, 100), (100, 500), (500, 1_000), (1_000, 5_000), (5_000, 9_000)) for capture in captures: expected = radio.bpsk_modulate(capture.payload_bits).astype(np.complex128) received = capture.payload_normalized frequencies, input_psd = welch( expected, fs=lab043.SYMBOL_RATE, window="hann", nperseg=1024, noverlap=512, return_onesided=False, scaling="density", ) _, cross_psd = csd( expected, received, fs=lab043.SYMBOL_RATE, window="hann", nperseg=1024, noverlap=512, return_onesided=False, scaling="density", ) response = np.abs(cross_psd / (input_psd + 1e-30)) for low, high in bands: selected = (np.abs(frequencies) >= low) & (np.abs(frequencies) < high) rows.append( { "capture": capture.name, "band_low_hz": low, "band_high_hz": high, "mean_transfer_magnitude": float(np.mean(response[selected])), "median_transfer_magnitude": float(np.median(response[selected])), "mean_transfer_db": float( 20.0 * np.log10(np.mean(response[selected]) + 1e-30) ), } ) return rows def fitted_symbol_channel_rows( captures: list[CaptureAnalysis], half_span: int = 8 ) -> tuple[list[dict], list[dict], list[dict]]: """Оценить на первой половине и проверить на второй gain/DC/FIR-модели.""" summaries: list[dict] = [] coefficient_rows: list[dict] = [] triplets: list[dict] = [] lags = np.arange(-half_span, half_span + 1) for capture in captures: expected = radio.bpsk_modulate(capture.payload_bits).astype(np.complex128) centers = np.arange(half_span, len(expected) - half_span) fir_matrix = np.column_stack([expected[centers - lag] for lag in lags]) target = capture.payload_normalized[centers] split = len(centers) // 2 train = slice(0, split) validation = slice(split, None) models = { "gain": fir_matrix[:, [half_span]], "gain+DC": np.column_stack((fir_matrix[:, half_span], np.ones(len(centers)))), "FIR": fir_matrix, "FIR+DC": np.column_stack((fir_matrix, np.ones(len(centers)))), } validation_bits = capture.payload_bits[centers[validation]] for model_name, matrix in models.items(): coefficients, *_ = np.linalg.lstsq(matrix[train], target[train], rcond=None) predicted = matrix[validation] @ coefficients actual = target[validation] decisions = (predicted.real >= 0.0).astype(np.uint8) errors = int(np.count_nonzero(decisions != validation_bits)) actual_rmse = float(np.sqrt(np.mean(np.abs(actual - predicted) ** 2))) summaries.append( { "capture": capture.name, "model": model_name, "training_symbols": split, "validation_symbols": len(predicted), "validation_rmse": actual_rmse, "predicted_bit_errors": errors, "predicted_ber": errors / len(predicted), "predicted_mean_i": float(np.mean(predicted).real), "predicted_mean_q": float(np.mean(predicted).imag), } ) source = f"{capture.name}: {model_name}, проверочная половина" triplets.extend(triplet_rows(source, validation_bits, predicted)) if model_name == "gain": coefficient_rows.append( { "capture": capture.name, "model": model_name, "term": "gain", "lag_symbols": 0, "coefficient_i": float(coefficients[0].real), "coefficient_q": float(coefficients[0].imag), "magnitude": float(abs(coefficients[0])), } ) elif model_name == "gain+DC": for term, coefficient in zip(("gain", "DC"), coefficients): coefficient_rows.append( { "capture": capture.name, "model": model_name, "term": term, "lag_symbols": 0 if term == "gain" else "", "coefficient_i": float(coefficient.real), "coefficient_q": float(coefficient.imag), "magnitude": float(abs(coefficient)), } ) else: fir_coefficients = coefficients[: len(lags)] for lag, coefficient in zip(lags, fir_coefficients): coefficient_rows.append( { "capture": capture.name, "model": model_name, "term": "tap", "lag_symbols": int(lag), "coefficient_i": float(coefficient.real), "coefficient_q": float(coefficient.imag), "magnitude": float(abs(coefficient)), } ) if model_name == "FIR+DC": coefficient = coefficients[-1] coefficient_rows.append( { "capture": capture.name, "model": model_name, "term": "DC", "lag_symbols": "", "coefficient_i": float(coefficient.real), "coefficient_q": float(coefficient.imag), "magnitude": float(abs(coefficient)), } ) return summaries, coefficient_rows, triplets def write_csv(path: Path, rows: list[dict]) -> None: if not rows: raise ValueError(f"Нет строк для {path}") with path.open("w", encoding="utf-8-sig", newline="") as stream: writer = csv.DictWriter(stream, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) def save_plots(output_dir: Path, phase_rows: list[dict], run_rows: list[dict]) -> None: figure, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True) for capture in sorted({row["capture"] for row in phase_rows}): selected = [row for row in phase_rows if row["capture"] == capture] axes[0].plot([row["phase"] for row in selected], [row["ber"] for row in selected], label=capture) axes[1].plot( [row["phase"] for row in selected], [row["eye_opening_i_p05_p95"] for row in selected], label=capture, ) axes[0].set_ylabel("Доля ошибок") axes[0].set_title("Lab043: зависимость от фазы символьной выборки") axes[0].grid(True, alpha=0.3) axes[0].legend(fontsize=8) axes[1].set_xlabel("Фаза, отсчётов из 120") axes[1].set_ylabel("Численное раскрытие глаза по I") axes[1].grid(True, alpha=0.3) figure.tight_layout() figure.savefig(output_dir / "phase_sweep.png", dpi=150) plt.close(figure) figure, axis = plt.subplots(figsize=(10, 5)) categories = ["1", "2", "3", "4", "5", "6", "7+"] for capture in sorted({row["capture"] for row in run_rows}): for bit in (0, 1): values = [] for category in categories: chosen = [ row for row in run_rows if row["capture"] == capture and row["bit"] == bit and row["run_length"] == category ] count = sum(row["symbol_count"] for row in chosen) values.append( sum(row["mean_signed_amplitude"] * row["symbol_count"] for row in chosen) / count if count else np.nan ) axis.plot(categories, values, marker="o", label=f"{capture}, бит {bit}") axis.set_title("Lab043: уровень символов и длина серии") axis.set_xlabel("Длина серии одинаковых битов") axis.set_ylabel("Средняя знаковая амплитуда") axis.grid(True, alpha=0.3) axis.legend(fontsize=8) figure.tight_layout() figure.savefig(output_dir / "run_length_dependence.png", dpi=150) plt.close(figure) def _summary_report( captures: list[CaptureAnalysis], capture_triplets: list[dict], tx_triplets: list[dict], ideal_triplets: list[dict], phase_rows: list[dict], dc_rows: list[dict], run_rows: list[dict], rrc_audit: dict, rrc_impulse: list[dict], spectra: list[dict], symbol_frequency: list[dict], channel_models: list[dict], channel_coefficients: list[dict], ) -> str: lines = ["Lab043: офлайн-диагностика асимметрии BPSK 000/001", ""] lines.append("Исходные данные") for capture in captures: lines.append( f"- {capture.name}: SHA-256 {capture.sha256}; корреляция маркера " f"{capture.marker_score:.6f}; диагностическая остаточная CFO " f"{capture.residual_cfo_hz:+.6f} Гц." ) lines.append("") lines.append("Ошибки по тройкам после общей линейной фазовой коррекции по известной PRBS") for capture in captures: rows = [row for row in capture_triplets if row["source"] == capture.name] total = sum(row["bit_errors"] for row in rows) target = sum(row["bit_errors"] for row in rows if row["triplet"] in {"000", "001"}) lines.append(f"- {capture.name}: всего {total}, в 000/001 {target} ({target / max(total, 1):.1%}).") for row in rows: lines.append( f" {row['triplet']}: n={row['count']}, ошибок={row['bit_errors']}, " f"BER={row['ber']:.6f}, mean={row['mean_i']:+.5f}{row['mean_q']:+.5f}j." ) lines.append("") lines.append("Цифровой TX и идеальный канал") for source in dict.fromkeys(row["source"] for row in tx_triplets): rows = [row for row in tx_triplets if row["source"] == source] spread = max(row["mean_amplitude"] for row in rows) - min(row["mean_amplitude"] for row in rows) lines.append( f"- {source}: ошибок {sum(row['bit_errors'] for row in rows)}, " f"размах средних амплитуд по тройкам {spread:.9f}." ) for source in dict.fromkeys(row["source"] for row in ideal_triplets): rows = [row for row in ideal_triplets if row["source"] == source] target = sum(row["bit_errors"] for row in rows if row["triplet"] in {"000", "001"}) total = sum(row["bit_errors"] for row in rows) lines.append(f"- Идеальный канал, {source}: ошибок {total}, из них 000/001 {target}.") lines.append("") lines.append("RRC") lines.append(json.dumps(rrc_audit, ensure_ascii=False, indent=2)) lines.append("h[kT]: " + ", ".join(f"h[{row['symbol_offset']:+d}T]={row['h_normalized']:+.9e}" for row in rrc_impulse)) lines.append("") lines.append("Фаза выборки") for capture in captures: rows = [row for row in phase_rows if row["capture"] == capture.name] selected = next(row for row in rows if row["phase"] == capture.sample_phase) best = min(rows, key=lambda row: row["bit_errors"]) lines.append( f"- {capture.name}: выбранная фаза {capture.sample_phase}, ошибок {selected['bit_errors']} " f"(000={selected['errors_000']}, 001={selected['errors_001']}, прочие={selected['errors_other']}), " f"глаз={selected['eye_opening_i_p05_p95']:+.5f}; минимум обзора на фазе {best['phase']}: " f"{best['bit_errors']} ошибок, глаз={best['eye_opening_i_p05_p95']:+.5f}. Фаза не изменялась." ) lines.append("") lines.append("DC, I/Q и расстояние до границы") for row in dc_rows: lines.append( f"- {row['capture']}: symbol DC={row['symbol_dc_i']:+.5f}{row['symbol_dc_q']:+.5f}j; " f"margin -1={row['negative_decision_margin_i']:.5f}, +1={row['positive_decision_margin_i']:.5f}; " f"ellipse={row['ellipse_axis_ratio']:.3f}, corr(I,Q)={row['iq_correlation']:+.3f}." ) lines.append("") lines.append("Серии одинаковых символов") for capture in captures: for bit in (0, 1): compact = [] for category in ("1", "2", "3", "4", "5", "6", "7+"): selected = [row for row in run_rows if row["capture"] == capture.name and row["bit"] == bit and row["run_length"] == category] count = sum(row["symbol_count"] for row in selected) if count: amplitude = sum(row["mean_signed_amplitude"] * row["symbol_count"] for row in selected) / count errors = sum(row["ber"] * row["symbol_count"] for row in selected) compact.append(f"{category}: A={amplitude:.4f}, BER={errors/count:.4f}") lines.append(f"- {capture.name}, бит {bit}: " + "; ".join(compact)) lines.append("") lines.append("Низкие частоты") for row in spectra: lines.append( f"- {row['source']}: mean={row['mean_i']:+.3e}{row['mean_q']:+.3e}j, " f"P(|f|<500 Гц)={row['power_fraction_abs_f_below_500_hz']:.6e}, " f"DC/median(2..10 кГц)={row['dc_bin_to_2_10khz_median_db']:+.2f} дБ." ) lines.append("Передаточная характеристика на символьной скорости по известной PRBS") for row in symbol_frequency: lines.append( f"- {row['capture']}, {row['band_low_hz']}..{row['band_high_hz']} Гц: " f"|H|={row['mean_transfer_magnitude']:.5f} ({row['mean_transfer_db']:+.2f} дБ)." ) lines.append("Диагностические модели символьного канала, обучение/проверка 50/50") for row in channel_models: lines.append( f"- {row['capture']}, {row['model']}: RMSE={row['validation_rmse']:.5f}, " f"предсказано ошибок {row['predicted_bit_errors']}/{row['validation_symbols']} " f"({row['predicted_ber']:.6f})." ) for capture in captures: dc = next( row for row in channel_coefficients if row["capture"] == capture.name and row["model"] == "FIR+DC" and row["term"] == "DC" ) taps = [ row for row in channel_coefficients if row["capture"] == capture.name and row["model"] == "FIR+DC" and row["term"] == "tap" ] center = next(row for row in taps if row["lag_symbols"] == 0) off_center = sum(row["magnitude"] for row in taps if row["lag_symbols"] != 0) lines.append( f"- {capture.name}: FIR+DC: DC={dc['coefficient_i']:+.5f}{dc['coefficient_q']:+.5f}j, " f"|h0|={center['magnitude']:.5f}, сумма |h[k!=0]|={off_center:.5f}." ) lines.append("") lines.append( "Интерпретация: таблицы являются измерениями; никакая коррекция DC, новый порог, " "эквалайзер или новая фаза автоматически не применялись." ) lines.append( "Классификация: B. Ошибки 000/001 отсутствуют в полном цифровом TX-RRC -> RX-RRC " "тракте и появляются только в реальном IQ. Промежуточная зависимость амплитуды после " "одного TX-RRC полностью компенсируется согласованным RX-RRC в идеальном канале." ) lines.append( "Наиболее вероятный механизм: совместное действие комплексного смещения середины " "кластеров и длинной памяти/ISI аппаратного тракта. FIR+DC-модель, обученная на первой " "половине и проверенная на второй, воспроизводит часть ошибок 000 без шума, но не весь " "измеренный BER; место внутри Pluto TX, кабеля/аттенюатора или RTL RX по этим данным не разделяется." ) lines.append( "AC coupling отдельно не доказан: падение уровня длинных серий устойчиво, но оценённая " "низкочастотная АЧХ различается между двумя S и не образует однозначной монотонной " "высокочастотной характеристики." ) lines.append( "Следующий диагностический шаг, пока не выполненный: сравнительный приём того же S " "другим приёмником либо loopback-разделение TX/RX с заранее фиксированными параметрами; " "после локализации отдельно проверить DC/АЧХ." ) return "\n".join(lines) + "\n" def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--raw-dir", type=Path, default=Path("data/raw/lab043")) parser.add_argument( "--output-dir", type=Path, default=Path("data/processed/lab043_diagnostics") ) args = parser.parse_args() args.output_dir.mkdir(parents=True, exist_ok=True) captures = [load_capture(args.raw_dir, name) for name in VALID_CAPTURE_NAMES] bits = captures[0].payload_bits capture_triplets: list[dict] = [] for capture in captures: capture_triplets.extend( triplet_rows(capture.name, capture.payload_bits, capture.payload_normalized) ) tx_triplets = tx_stage_rows(bits) ideal_triplets = ideal_channel_rows(bits) rrc_impulse, rrc_audit = rrc_rows() phase_rows = [row for capture in captures for row in phase_sweep_rows(capture)] dc_rows = [row for capture in captures for row in dc_iq_rows(capture)] run_rows = [row for capture in captures for row in run_length_rows(capture)] spectra = spectrum_rows(captures, bits) symbol_frequency = symbol_frequency_rows(captures) channel_models, channel_coefficients, channel_triplets = fitted_symbol_channel_rows(captures) write_csv(args.output_dir / "capture_triplets.csv", capture_triplets) write_csv(args.output_dir / "tx_stages_triplets.csv", tx_triplets) write_csv(args.output_dir / "ideal_channel_triplets.csv", ideal_triplets) write_csv(args.output_dir / "rrc_impulse_response.csv", rrc_impulse) write_csv(args.output_dir / "phase_sweep.csv", phase_rows) write_csv(args.output_dir / "dc_iq_metrics.csv", dc_rows) write_csv(args.output_dir / "run_length_metrics.csv", run_rows) write_csv(args.output_dir / "spectrum_metrics.csv", spectra) write_csv(args.output_dir / "symbol_frequency_response.csv", symbol_frequency) write_csv(args.output_dir / "symbol_channel_models.csv", channel_models) write_csv(args.output_dir / "symbol_channel_coefficients.csv", channel_coefficients) write_csv(args.output_dir / "symbol_channel_triplets.csv", channel_triplets) (args.output_dir / "rrc_audit.json").write_text( json.dumps(rrc_audit, ensure_ascii=False, indent=2) + "\n", encoding="utf-8" ) (args.output_dir / "input_manifest.json").write_text( json.dumps( [ { "metadata_file": f"{capture.name}.json", "iq_file": capture.metadata["iq_file"], "verified_sha256": capture.sha256, } for capture in captures ], ensure_ascii=False, indent=2, ) + "\n", encoding="utf-8", ) save_plots(args.output_dir, phase_rows, run_rows) report = _summary_report( captures, capture_triplets, tx_triplets, ideal_triplets, phase_rows, dc_rows, run_rows, rrc_audit, rrc_impulse, spectra, symbol_frequency, channel_models, channel_coefficients, ) (args.output_dir / "report.txt").write_text(report, encoding="utf-8") print(report) return 0 if __name__ == "__main__": raise SystemExit(main())