Files
SDR-Rover/experiments/lab043_bpsk_asymmetry_diagnostics.py

1035 lines
46 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
"""Офлайн-диагностика асимметрии ошибок 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())