313 lines
18 KiB
Python
313 lines
18 KiB
Python
"""Lab056: программная подготовка ABAB и диагностика сохранённого RX."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import json
|
||
from pathlib import Path
|
||
import sys
|
||
|
||
if __package__ in (None, ""):
|
||
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
||
|
||
import numpy as np
|
||
from scipy import signal
|
||
|
||
from experiments import lab043_pluto_to_rtlsdr as shaping
|
||
from experiments import lab055_energy_audit as energy
|
||
|
||
SAMPLE_RATE = 8_000_000
|
||
CARRIER = 50_000
|
||
GAP_SAMPLES = 800_000
|
||
RX_SAMPLES = 24_000_000
|
||
RX_LEAD_SECONDS = .5
|
||
RAW_DIRECTORY = energy.PROJECT_ROOT / "data/raw/lab056"
|
||
OUTPUT = energy.PROJECT_ROOT / "data/processed/lab056"
|
||
|
||
|
||
def save_new(path: Path, value: dict) -> None:
|
||
text = json.dumps(value, ensure_ascii=False, indent=2, allow_nan=False) + "\n"
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
with path.open("x", encoding="utf-8") as stream:
|
||
stream.write(text)
|
||
|
||
|
||
def build_burst(symbol_rate: int) -> np.ndarray:
|
||
"""Прежние PRBS7/BPSK/RRC, 204,8 мс символов, RMS всего burst=44."""
|
||
if symbol_rate not in (20000, 500000):
|
||
raise ValueError("Разрешены только два исследуемых режима")
|
||
symbols = 2. * shaping.pilot_prbs7(4096 * (symbol_rate // 20000)) - 1.
|
||
samples_per_symbol = SAMPLE_RATE // symbol_rate
|
||
waveform = shaping.shape_frame_like_lab042(symbols, samples_per_symbol)
|
||
waveform *= np.exp(2j * np.pi * CARRIER * np.arange(len(waveform)) / SAMPLE_RATE)
|
||
waveform *= 44 / np.sqrt(np.mean(abs(waveform)**2))
|
||
values = np.empty(2 * len(waveform), dtype=np.int8)
|
||
values[::2] = np.clip(np.rint(waveform.real), -128, 127).astype(np.int8)
|
||
values[1::2] = np.clip(np.rint(waveform.imag), -128, 127).astype(np.int8)
|
||
return values
|
||
|
||
|
||
def prepare(directory: Path = RAW_DIRECTORY) -> dict:
|
||
"""Создать единственный конечный TX-файл; аппаратная частота не назначена."""
|
||
directory.mkdir(parents=True, exist_ok=True)
|
||
tx_path, metadata_path = directory / "tx_abab.cs8", directory / "tx_metadata.json"
|
||
if tx_path.exists() or metadata_path.exists():
|
||
raise FileExistsError("Артефакты Lab056 уже существуют")
|
||
bursts = {rate: build_burst(rate) for rate in (20000, 500000)}
|
||
segments = []
|
||
cursor = GAP_SAMPLES
|
||
with tx_path.open("xb") as stream:
|
||
stream.write(bytes(2 * GAP_SAMPLES))
|
||
for label, rate in (("A1", 20000), ("B1", 500000), ("A2", 20000), ("B2", 500000)):
|
||
values = bursts[rate]
|
||
count = len(values) // 2
|
||
complex_values = values[::2].astype(np.float32) + 1j * values[1::2].astype(np.float32)
|
||
whole, active = energy.Moments(), energy.Moments()
|
||
whole.add(complex_values)
|
||
trim = 21 * (SAMPLE_RATE // rate)
|
||
active.add(complex_values[trim:-trim])
|
||
segments.append({"label": label, "symbol_rate": rate, "start": cursor,
|
||
"samples": count, "active_start": cursor + trim, "active_samples": count - 2 * trim,
|
||
"whole": whole.result(rate), "active": active.result(rate)})
|
||
stream.write(values.tobytes())
|
||
stream.write(bytes(2 * GAP_SAMPLES))
|
||
cursor += count + GAP_SAMPLES
|
||
result = {"stage": "PREPARATION_ONLY", "rf_tx": 0, "rf_frequency_hz": None,
|
||
"physical_path": "OTA_REQUESTED_FREQUENCY_AND_GEOMETRY_PENDING", "sample_rate_hz": SAMPLE_RATE,
|
||
"carrier_hz": CARRIER, "tx_samples": cursor, "tx_bytes": 2 * cursor,
|
||
"tx_seconds": cursor / SAMPLE_RATE, "rx_samples": RX_SAMPLES,
|
||
"rx_seconds": RX_SAMPLES / SAMPLE_RATE, "rx_lead_seconds": RX_LEAD_SECONDS,
|
||
"gap_seconds": GAP_SAMPLES / SAMPLE_RATE, "gap_is_zero_digital_iq_not_rf_off": True,
|
||
"tx_sha256": energy.sha256(tx_path), "segments": segments,
|
||
"whole_power_difference_db": float(10 * np.log10(segments[1]["whole"]["power_code2"] / segments[0]["whole"]["power_code2"])),
|
||
"active_power_difference_db": float(10 * np.log10(segments[1]["active"]["power_code2"] / segments[0]["active"]["power_code2"]))}
|
||
if abs(result["whole_power_difference_db"]) > .01 or any(s["whole"]["clipped_samples_fraction"] for s in segments):
|
||
raise RuntimeError("Цифровая подготовка не прошла проверку мощности/клиппинга")
|
||
save_new(metadata_path, result)
|
||
return result
|
||
|
||
|
||
def locate_sequence(tx: energy.Cs8File, rx: energy.Cs8File) -> dict:
|
||
"""Диагностическое согласование полной ABAB; TX reference разрешён.
|
||
|
||
Грубая сетка 32 отсчёта использует всю последовательность, затем начало
|
||
каждого burst уточняется прежним методом Lab055 на обучающей плитке.
|
||
Прореживание без антиалиасинга здесь только для корреляционного поиска,
|
||
а не для измерений мощности или спектра.
|
||
"""
|
||
stride = 32
|
||
def decimated(source):
|
||
chunks = [source.read(start, min(262144, source.size - start))[::stride]
|
||
for start in range(0, source.size, 262144)]
|
||
return np.concatenate(chunks).astype(np.complex128)
|
||
reference, received = decimated(tx), decimated(rx)
|
||
if len(received) < len(reference):
|
||
raise ValueError("RX короче TX")
|
||
reference -= reference.mean()
|
||
received -= received.mean()
|
||
cumulative = np.concatenate(([0.], np.cumsum(abs(received)**2)))
|
||
window_energy = cumulative[len(reference):] - cumulative[:-len(reference)]
|
||
reference_energy = np.sum(abs(reference)**2)
|
||
time = np.arange(len(reference)) * stride / SAMPLE_RATE
|
||
best = (-1., 0, 0.)
|
||
for frequency in np.arange(-2., 2.01, .25):
|
||
rotated = reference * np.exp(2j * np.pi * frequency * time)
|
||
correlation = signal.fftconvolve(received, np.conj(rotated[::-1]), mode="valid")
|
||
scores = abs(correlation) / np.sqrt(np.maximum(reference_energy * window_energy, 1e-30))
|
||
index = int(np.argmax(scores))
|
||
if scores[index] > best[0]:
|
||
best = (float(scores[index]), index * stride, float(frequency))
|
||
# Это качество диагностического совпадения, не критерий PHY/CRC.
|
||
if best[0] < .05:
|
||
raise ValueError("ABAB не локализована достаточно уверенно; измерение UNKNOWN")
|
||
return {"start": best[1], "correlation": best[0], "coarse_cfo_hz": best[2], "stride": stride}
|
||
|
||
|
||
class OffsetSource:
|
||
def __init__(self, source: energy.Cs8File, offset: int):
|
||
self.source, self.offset = source, offset
|
||
|
||
def read(self, start: int, count: int) -> np.ndarray:
|
||
return self.source.read(self.offset + start, count)
|
||
|
||
|
||
def window_power(values: np.ndarray, symbol_rate: int) -> dict:
|
||
moments = energy.Moments()
|
||
moments.add(values)
|
||
result = moments.result(symbol_rate)
|
||
centered = values - values.mean()
|
||
frequencies, density = energy.spectrum(centered)
|
||
for width in (27000, 675000):
|
||
result[f"band_{width}_power"] = energy.band_power(frequencies, density, CARRIER - width / 2, CARRIER + width / 2)
|
||
size = 32000 # Неперекрывающиеся блоки по 4 мс для оценки разброса мощности.
|
||
powers = [float(np.mean(abs(centered[start:start + size])**2)) for start in range(0, len(centered) - size + 1, size)]
|
||
result["power_standard_error"] = float(np.std(powers, ddof=1) / np.sqrt(len(powers))) if len(powers) > 1 else None
|
||
return result
|
||
|
||
|
||
def difference_db(new: float | None, old: float | None) -> float | None:
|
||
if new is None or old is None or new <= 0 or old <= 0:
|
||
return None
|
||
return float(10 * np.log10(new / old))
|
||
|
||
|
||
def complex_gain(reference: np.ndarray, received: np.ndarray, cfo_hz: float) -> dict:
|
||
"""Вывести комплексный h прежней модели; те же чётные обучающие плитки."""
|
||
boundaries = np.linspace(0, len(reference), 9, dtype=int)
|
||
train = np.zeros(len(reference), dtype=bool)
|
||
for tile in range(0, 8, 2):
|
||
train[boundaries[tile]:boundaries[tile + 1]] = True
|
||
rotated = reference * np.exp(2j * np.pi * cfo_hz * np.arange(len(reference)) / SAMPLE_RATE)
|
||
x, y = rotated[train], received[train].astype(np.complex128)
|
||
x, y = x - x.mean(), y - y.mean()
|
||
gain = np.vdot(x, y) / max(float(np.sum(abs(x)**2)), 1e-30)
|
||
return {"real": float(gain.real), "imag": float(gain.imag),
|
||
"phase_rad_at_active_start": float(np.angle(gain))}
|
||
|
||
|
||
def comparison_summary(rows: list[dict]) -> dict:
|
||
"""Сравнить повторы отдельно; X/Z используют средние линейные мощности."""
|
||
by_label = {row["label"]: row for row in rows}
|
||
repeatability = []
|
||
for mode in ("A", "B"):
|
||
first, second = by_label[mode + "1"], by_label[mode + "2"]
|
||
repeatability.append({"comparison": mode + "2/" + mode + "1",
|
||
"raw_rx_power_difference_db": difference_db(second["raw_rx"]["power_code2"], first["raw_rx"]["power_code2"]),
|
||
"background_difference_db": difference_db(second["background_power"], first["background_power"]),
|
||
"received_power_difference_db": difference_db(second["background_subtracted_signal_power"], first["background_subtracted_signal_power"]),
|
||
"channel_gain_difference_db": difference_db(second["channel_gain_power"], first["channel_gain_power"]),
|
||
"matched_component_difference_db": difference_db(second["measurement"]["heldout_signal_power"], first["measurement"]["heldout_signal_power"])})
|
||
def mean(mode, key):
|
||
values = [by_label[mode + str(number)][key] for number in (1, 2)]
|
||
return None if any(value is None or value <= 0 for value in values) else float(np.mean(values))
|
||
a_power, b_power = mean("A", "background_subtracted_signal_power"), mean("B", "background_subtracted_signal_power")
|
||
return {"repeatability": repeatability,
|
||
"averaging": "RATIO_OF_MEAN_LINEAR_POWERS_B1_B2_OVER_A1_A2",
|
||
"received_power_difference_db": difference_db(b_power, a_power),
|
||
"channel_gain_difference_db": difference_db(mean("B", "channel_gain_power"), mean("A", "channel_gain_power")),
|
||
"energy_per_symbol_difference_db": difference_db(None if b_power is None else b_power / 500000,
|
||
None if a_power is None else a_power / 20000),
|
||
"energy_basis": "BACKGROUND_SUBTRACTED_RECEIVED_POWER_DIVIDED_BY_SYMBOL_RATE"}
|
||
|
||
|
||
def analyze(directory: Path, rx_path: Path) -> dict:
|
||
metadata = energy.read_json(directory / "tx_metadata.json")
|
||
tx_path = directory / "tx_abab.cs8"
|
||
if energy.sha256(tx_path) != metadata["tx_sha256"]:
|
||
raise ValueError("TX SHA не совпал")
|
||
tx, rx = energy.Cs8File(tx_path), energy.Cs8File(rx_path)
|
||
alignment = locate_sequence(tx, rx)
|
||
rows = []
|
||
for segment in metadata["segments"]:
|
||
rate, count = segment["symbol_rate"], segment["samples"]
|
||
start = alignment["start"] + segment["start"]
|
||
measurement = energy.analyze_rx(OffsetSource(tx, segment["start"]), rx, [start], count,
|
||
SAMPLE_RATE // rate)["frames"][0]
|
||
start = measurement["start"]
|
||
trim = 21 * (SAMPLE_RATE // rate)
|
||
received_active = rx.read(start + trim, count - 2 * trim)
|
||
raw = window_power(received_active, rate)
|
||
gain = complex_gain(tx.read(segment["start"] + trim, count - 2 * trim), received_active, measurement["cfo_hz"])
|
||
if not np.isclose(gain["real"]**2 + gain["imag"]**2, measurement["gain_abs"]**2, rtol=1e-10, atol=1e-15):
|
||
raise ValueError("Вывод комплексного h расходится с прежней моделью Lab055")
|
||
windows = [window_power(rx.read(begin, 480000), rate) for begin in
|
||
(start - 640000, start + count + 160000)]
|
||
background = float(np.mean([w["ac_power_code2"] for w in windows]))
|
||
excess = raw["ac_power_code2"] - background
|
||
uncertainty = float(np.sqrt(raw["power_standard_error"]**2 +
|
||
sum(w["power_standard_error"]**2 for w in windows) / 4 +
|
||
((windows[0]["ac_power_code2"] - windows[1]["ac_power_code2"]) / 2)**2))
|
||
stable = excess > 3 * uncertainty
|
||
rows.append({"label": segment["label"], "symbol_rate": rate, "raw_rx": raw,
|
||
"local_zero_iq_background_windows": windows, "background_power": background,
|
||
"background_subtracted_raw": excess, "background_subtracted_uncertainty": uncertainty,
|
||
"background_subtracted_stable": stable, "background_subtracted_signal_power": excess if stable else None,
|
||
"channel_gain_complex": gain, "channel_gain_power": measurement["gain_abs"]**2, "measurement": measurement})
|
||
by_label = {r["label"]: r for r in rows}
|
||
pairs = []
|
||
for number in (1, 2):
|
||
a, b = by_label[f"A{number}"], by_label[f"B{number}"]
|
||
pairs.append({"pair": number,
|
||
"received_power_difference_db": difference_db(b["background_subtracted_signal_power"], a["background_subtracted_signal_power"]),
|
||
"channel_gain_difference_db": difference_db(b["channel_gain_power"], a["channel_gain_power"]),
|
||
"matched_component_difference_db": difference_db(b["measurement"]["heldout_signal_power"], a["measurement"]["heldout_signal_power"]),
|
||
"matched_component_symbol_energy_difference_db": difference_db(b["measurement"]["heldout_signal_power"] / b["symbol_rate"], a["measurement"]["heldout_signal_power"] / a["symbol_rate"])})
|
||
# До начала и после конца всего цифрового файла: отдельные внешние окна.
|
||
outside = []
|
||
for begin in (alignment["start"] - 800000, alignment["start"] + tx.size + 400000):
|
||
if 0 <= begin and begin + 400000 <= rx.size:
|
||
outside.append({"start": begin, **window_power(rx.read(begin, 400000), 20000)})
|
||
return {"scope": "TX_REFERENCE_DIAGNOSTICS_NOT_INDEPENDENT_PACKET_RX", "rf_tx": 0,
|
||
"tx_sha256": metadata["tx_sha256"], "rx_sha256": energy.sha256(rx_path),
|
||
"alignment": alignment, "segments": rows, "pairs": pairs, "comparison_summary": comparison_summary(rows),
|
||
"outside_waveform_background": outside,
|
||
"background_note": "Внутренние gaps: нулевой цифровой IQ при работающем TX; это не RF OFF. Внешние окна также требуют сверки времён запуска/остановки по журналам.",
|
||
"power_units": "RX code squared, not calibrated watts/dBm",
|
||
"residual_note": "Шум и неописанные моделью искажения; RX не нормируется для выравнивания A/B"}
|
||
|
||
|
||
def software_capture(directory: Path, b_gain: float, name: str) -> dict:
|
||
"""Синтетический CS8 с известным усилением, DC и CFO; не аппаратная запись."""
|
||
metadata = energy.read_json(directory / "tx_metadata.json")
|
||
tx = energy.Cs8File(directory / "tx_abab.cs8")
|
||
path = directory / f"synthetic_{name}.cs8"
|
||
random_generator = np.random.default_rng(56056)
|
||
lead = 4_160_013
|
||
with path.open("xb") as stream:
|
||
for low in range(0, RX_SAMPLES, 262144):
|
||
count = min(262144, RX_SAMPLES - low)
|
||
values = (1 + 2j + .6 * (random_generator.normal(size=count) + 1j * random_generator.normal(size=count)))
|
||
first, last = max(low, lead), min(low + count, lead + tx.size)
|
||
if first < last:
|
||
source = tx.read(first - lead, last - first)
|
||
gains = np.full(len(source), .03)
|
||
indexes = np.arange(first - lead, last - lead)
|
||
for segment in metadata["segments"]:
|
||
if segment["label"].startswith("B"):
|
||
gains[(indexes >= segment["start"]) & (indexes < segment["start"] + segment["samples"])] = b_gain
|
||
values[first - low:last - low] += source * gains * np.exp(1j * (.7 + 2 * np.pi * .2 * np.arange(first, last) / SAMPLE_RATE))
|
||
packed = np.empty(2 * count, dtype=np.int8)
|
||
packed[::2] = np.clip(np.rint(values.real), -128, 127)
|
||
packed[1::2] = np.clip(np.rint(values.imag), -128, 127)
|
||
stream.write(packed.tobytes())
|
||
result = analyze(directory, path)
|
||
expected = 20 * np.log10(b_gain / .03)
|
||
errors = [abs(p["channel_gain_difference_db"] - expected) for p in result["pairs"]]
|
||
result.update(synthetic=True, expected_gain_difference_db=float(expected),
|
||
known_sequence_start=lead, passed=bool(max(errors) < .3))
|
||
if not result["passed"]:
|
||
raise RuntimeError(f"Синтетическая проверка не прошла: {errors}")
|
||
return result
|
||
|
||
|
||
def run_functional_tests() -> None:
|
||
assert difference_db(.25, 1) == 10 * np.log10(.25)
|
||
assert difference_db(None, 1) is None
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description=__doc__)
|
||
parser.add_argument("stage", choices=("prepare", "software", "analyze"))
|
||
parser.add_argument("--rx", type=Path)
|
||
parser.add_argument("--output", type=Path)
|
||
args = parser.parse_args()
|
||
if args.stage == "prepare":
|
||
result = prepare()
|
||
save_new(OUTPUT / "preparation.json", result)
|
||
elif args.stage == "software":
|
||
for name, gain in (("equal", .03), ("minus6db", .015)):
|
||
result = software_capture(RAW_DIRECTORY, gain, name)
|
||
save_new(OUTPUT / f"software_{name}.json", result)
|
||
print(name, result["pairs"], flush=True)
|
||
return
|
||
else:
|
||
if args.rx is None or args.output is None:
|
||
parser.error("Для сохранённого RX нужны --rx и новый --output")
|
||
result = analyze(RAW_DIRECTORY, args.rx)
|
||
save_new(args.output, result)
|
||
print(json.dumps(result, ensure_ascii=False, indent=2))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|