385 lines
21 KiB
Python
385 lines
21 KiB
Python
"""Lab055: энергия сохранённых IQ; аппаратных зависимостей и команд нет."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
from dataclasses import dataclass
|
||
import hashlib
|
||
import json
|
||
from pathlib import Path
|
||
import sys
|
||
|
||
if __package__ in (None, ""):
|
||
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
|
||
|
||
import numpy as np
|
||
from scipy import signal
|
||
|
||
from protocol.bpsk_radio import root_raised_cosine_taps
|
||
|
||
SAMPLE_RATE = 8_000_000
|
||
CARRIER = 50_000
|
||
SPECTRUM_SIZE = 65536
|
||
PROJECT_ROOT = Path(__file__).resolve().parents[1]
|
||
OLD_DIRECTORY = Path(r"C:\Users\user\Documents\SDR_Rover_Lab047_raw\20260903_115747")
|
||
NEW_DIRECTORY = PROJECT_ROOT / "data/raw/lab054"
|
||
|
||
|
||
def read_json(path: Path) -> dict:
|
||
return json.loads(path.read_text(encoding="utf-8"))
|
||
|
||
|
||
def sha256(path: Path) -> str:
|
||
"""Хеш без загрузки файла целиком."""
|
||
digest = hashlib.sha256()
|
||
with path.open("rb") as stream:
|
||
while block := stream.read(1024 * 1024):
|
||
digest.update(block)
|
||
return digest.hexdigest()
|
||
|
||
|
||
class Cs8File:
|
||
"""Только чтение interleaved signed IQ; в память извлекается одно окно."""
|
||
|
||
def __init__(self, path: Path):
|
||
if path.stat().st_size == 0 or path.stat().st_size % 2:
|
||
raise ValueError("Нужен непустой CS8-файл с полными IQ-парами")
|
||
self.values = np.memmap(path, mode="r", dtype=np.int8)
|
||
self.size = self.values.size // 2
|
||
|
||
def read(self, start: int, count: int) -> np.ndarray:
|
||
if start < 0 or count <= 0 or start + count > self.size:
|
||
raise ValueError("Окно вне сохранённой записи")
|
||
values = self.values[2 * start:2 * (start + count)]
|
||
return values[::2].astype(np.float32) + 1j * values[1::2].astype(np.float32)
|
||
|
||
|
||
@dataclass
|
||
class Moments:
|
||
"""Точные накопители моментов цифровых отсчётов, без массивов всей записи."""
|
||
|
||
count: int = 0
|
||
total: complex = 0j
|
||
energy: float = 0.0
|
||
peak_power: float = 0.0
|
||
max_i: float = 0.0
|
||
max_q: float = 0.0
|
||
clipped_components: int = 0
|
||
clipped_samples: int = 0
|
||
|
||
def add(self, values: np.ndarray) -> None:
|
||
power = values.real.astype(float)**2 + values.imag.astype(float)**2
|
||
clip_i = (values.real == -128) | (values.real == 127)
|
||
clip_q = (values.imag == -128) | (values.imag == 127)
|
||
self.count += values.size
|
||
self.total += complex(np.sum(values, dtype=np.complex128))
|
||
self.energy += float(np.sum(power))
|
||
self.peak_power = max(self.peak_power, float(np.max(power)))
|
||
self.max_i = max(self.max_i, float(np.max(abs(values.real))))
|
||
self.max_q = max(self.max_q, float(np.max(abs(values.imag))))
|
||
self.clipped_components += int(np.sum(clip_i) + np.sum(clip_q))
|
||
self.clipped_samples += int(np.sum(clip_i | clip_q))
|
||
|
||
def result(self, symbol_rate: int) -> dict:
|
||
mean = self.total / self.count
|
||
power = self.energy / self.count
|
||
return {"samples": self.count, "mean_i": mean.real, "mean_q": mean.imag,
|
||
"power_code2": power, "ac_power_code2": power - abs(mean)**2,
|
||
"rms_code": float(np.sqrt(power)), "peak_magnitude_code": float(np.sqrt(self.peak_power)),
|
||
"max_abs_i": self.max_i, "max_abs_q": self.max_q,
|
||
"papr_db": float(10 * np.log10(self.peak_power / power)),
|
||
"clipped_components_fraction": self.clipped_components / (2 * self.count),
|
||
"clipped_samples_fraction": self.clipped_samples / self.count,
|
||
"symbol_energy_code2_seconds": power / symbol_rate}
|
||
|
||
|
||
def spectrum(values: np.ndarray, sample_rate: int = SAMPLE_RATE) -> tuple[np.ndarray, np.ndarray]:
|
||
"""Двусторонняя PSD комплексного сигнала; интеграл имеет единицы код²."""
|
||
size = min(SPECTRUM_SIZE, values.size)
|
||
frequencies, density = signal.welch(values, fs=sample_rate, window="hann",
|
||
nperseg=size, noverlap=size // 2, detrend=False, return_onesided=False, scaling="density")
|
||
return frequencies, density.astype(float)
|
||
|
||
|
||
def band_power(frequencies: np.ndarray, density: np.ndarray, low: float, high: float) -> float:
|
||
width = abs(float(frequencies[1] - frequencies[0]))
|
||
return float(np.sum(density[(frequencies >= low) & (frequencies < high)]) * width)
|
||
|
||
|
||
def spread(values: list[float]) -> dict:
|
||
return {"min": float(np.min(values)), "median": float(np.median(values)),
|
||
"max": float(np.max(values)), "mean": float(np.mean(values))}
|
||
|
||
|
||
def filter_metrics(samples_per_symbol: int) -> tuple[np.ndarray, dict]:
|
||
taps = root_raised_cosine_taps(.35, samples_per_symbol, 10)
|
||
energy = float(np.sum(abs(taps)**2))
|
||
dc_gain = float(abs(np.sum(taps))**2)
|
||
return taps, {"rolloff": .35, "span_symbols": 10, "taps": len(taps),
|
||
"group_delay_samples": (len(taps) - 1) / 2,
|
||
"energy": energy, "dc_power_gain": dc_gain,
|
||
"equivalent_noise_bandwidth_hz": SAMPLE_RATE * energy / dc_gain}
|
||
|
||
|
||
def matched_noise(frequencies: np.ndarray, density: np.ndarray, taps: np.ndarray,
|
||
carrier: float = CARRIER) -> float:
|
||
"""Интеграл PSD*|H|² для фактического единично-энергетического RRC."""
|
||
shifted = taps * np.exp(2j * np.pi * carrier * np.arange(len(taps)) / SAMPLE_RATE)
|
||
response = np.fft.fft(shifted, len(frequencies))
|
||
return float(np.sum(density * abs(response)**2) * SAMPLE_RATE / len(frequencies))
|
||
|
||
|
||
def align_training(reference: np.ndarray, received: np.ndarray, margin: int) -> int:
|
||
"""Задержка по обучающему окну; проверочные участки не участвуют."""
|
||
reference = reference - reference.mean()
|
||
received = received - received.mean()
|
||
correlation = signal.fftconvolve(received, np.conj(reference[::-1]), mode="valid")
|
||
energy = abs(received)**2
|
||
cumulative = np.concatenate(([0.], np.cumsum(energy, dtype=float)))
|
||
windows = cumulative[len(reference):] - cumulative[:-len(reference)]
|
||
scores = abs(correlation) / np.sqrt(np.maximum(windows, 1e-30))
|
||
return int(np.argmax(scores)) - margin
|
||
|
||
|
||
def fit_channel(reference: np.ndarray, received: np.ndarray, sample_rate: int,
|
||
train_parity: int = 0, spectral_taps: np.ndarray | None = None) -> dict:
|
||
"""Модель y=d+g*x*exp(j2πft); 8 плиток, чередующиеся обучение/проверка.
|
||
|
||
Частота сначала оценивается по фазе четырёх обучающих проекций.
|
||
g и d оцениваются только по обучающим отсчётам. Проверочный остаток
|
||
содержит шум, частотную/временную ошибку и прочие искажения модели.
|
||
"""
|
||
count = len(reference)
|
||
if count != len(received) or count < 128:
|
||
raise ValueError("Несопоставимые или слишком короткие окна")
|
||
boundaries = np.linspace(0, count, 9, dtype=int)
|
||
train = np.zeros(count, dtype=bool)
|
||
times, gains, weights = [], [], []
|
||
for tile in range(train_parity, 8, 2):
|
||
low, high = boundaries[tile:tile + 2]
|
||
train[low:high] = True
|
||
x = reference[low:high].astype(np.complex128)
|
||
y = received[low:high].astype(np.complex128)
|
||
x -= x.mean()
|
||
y -= y.mean()
|
||
energy = float(np.vdot(x, x).real)
|
||
gain = np.vdot(x, y) / max(energy, 1e-30)
|
||
gains.append(gain)
|
||
times.append((low + high - 1) / (2 * sample_rate))
|
||
weights.append(abs(gain)**2 * energy)
|
||
times = np.asarray(times)
|
||
matrix = np.column_stack((np.ones(4), times - times.mean()))
|
||
roots = np.sqrt(np.maximum(weights, 1e-30))
|
||
phase = np.unwrap(np.angle(gains))
|
||
fit = np.linalg.lstsq(matrix * roots[:, None], phase * roots, rcond=None)[0]
|
||
frequency = float(fit[1] / (2 * np.pi))
|
||
model_reference = reference * np.exp(2j * np.pi * frequency * np.arange(count) / sample_rate)
|
||
x = model_reference[train]
|
||
y = received[train].astype(np.complex128)
|
||
x_mean, y_mean = x.mean(), y.mean()
|
||
gain = np.vdot(x - x_mean, y - y_mean) / max(float(np.sum(abs(x - x_mean)**2)), 1e-30)
|
||
dc = y_mean - gain * x_mean
|
||
result = {"cfo_hz": frequency, "gain_abs": float(abs(gain)),
|
||
"dc_i": float(dc.real), "dc_q": float(dc.imag)}
|
||
for label, mask in (("train", train), ("heldout", ~train)):
|
||
predicted = gain * model_reference[mask]
|
||
actual = received[mask] - dc
|
||
residual = actual - predicted
|
||
result[label + "_signal_power"] = float(np.mean(abs(predicted)**2))
|
||
result[label + "_residual_power"] = float(np.mean(abs(residual)**2))
|
||
result[label + "_actual_ac_power"] = float(np.mean(abs(actual)**2))
|
||
result[label + "_cross_power"] = float(2 * np.mean((predicted.conj() * residual).real))
|
||
if spectral_taps is not None:
|
||
# Отдельные плитки: не склеиваем разорванные окна перед оценкой PSD.
|
||
spectral_rows = []
|
||
for tile in range(1 - train_parity, 8, 2):
|
||
low, high = boundaries[tile:tile + 2]
|
||
residual = received[low:high] - dc - gain * model_reference[low:high]
|
||
frequencies, density = spectrum(residual, sample_rate)
|
||
spectral_rows.append({
|
||
**{f"heldout_residual_band_{width}_power": band_power(
|
||
frequencies, density, CARRIER - width / 2, CARRIER + width / 2)
|
||
for width in (27000, 675000, 1000000)},
|
||
"heldout_residual_matched_power": matched_noise(frequencies, density, spectral_taps)})
|
||
result.update({key: float(np.mean([row[key] for row in spectral_rows])) for key in spectral_rows[0]})
|
||
return result
|
||
|
||
|
||
def analyze_tx(path: Path, frame_samples: int, frames: int, symbol_rate: int) -> dict:
|
||
source = Cs8File(path)
|
||
samples_per_symbol = SAMPLE_RATE // symbol_rate
|
||
if source.size != frame_samples * frames:
|
||
raise ValueError("Размер TX не соответствует метаданным")
|
||
moments = {"whole": Moments(), "active": Moments()}
|
||
rows = []
|
||
for frame in range(frames):
|
||
values = source.read(frame * frame_samples, frame_samples)
|
||
active = values[21 * samples_per_symbol:-21 * samples_per_symbol]
|
||
row = {"frame": frame + 1}
|
||
for label, array in (("whole", values), ("active", active)):
|
||
moments[label].add(array)
|
||
local = Moments()
|
||
local.add(array)
|
||
row[label] = local.result(symbol_rate)
|
||
# Прямоугольная периодограмма полного окна сохраняет Parseval,
|
||
# включая защитные интервалы. Welch недовзвешивает края файла.
|
||
frequencies, density = signal.periodogram(array, fs=SAMPLE_RATE,
|
||
window="boxcar", detrend=False, return_onesided=False, scaling="density")
|
||
row[label]["occupied_band_power_code2"] = band_power(
|
||
frequencies, density, CARRIER - .675 * symbol_rate, CARRIER + .675 * symbol_rate)
|
||
rows.append(row)
|
||
return {"whole": moments["whole"].result(symbol_rate),
|
||
"active": moments["active"].result(symbol_rate), "frames": rows,
|
||
"occupied_band_power_whole": spread([r["whole"]["occupied_band_power_code2"] for r in rows]),
|
||
"occupied_band_power_active": spread([r["active"]["occupied_band_power_code2"] for r in rows])}
|
||
|
||
|
||
def analyze_background(source: Cs8File, seconds: list[float], samples_per_symbol: int) -> dict:
|
||
taps, metrics = filter_metrics(samples_per_symbol)
|
||
rows = []
|
||
for start_seconds in seconds:
|
||
values = source.read(round(start_seconds * SAMPLE_RATE), SAMPLE_RATE // 10)
|
||
moments = Moments()
|
||
moments.add(values)
|
||
row = {"start_seconds": start_seconds, **moments.result(SAMPLE_RATE // samples_per_symbol)}
|
||
values -= values.mean()
|
||
frequencies, density = spectrum(values)
|
||
for width in (27000, 675000, 1000000):
|
||
power = band_power(frequencies, density, CARRIER - width / 2, CARRIER + width / 2)
|
||
row[f"band_{width}_power"] = power
|
||
row[f"band_{width}_density"] = power / width
|
||
row["matched_noise_power"] = matched_noise(frequencies, density, taps)
|
||
row["equivalent_n0_code2_per_hz"] = row["matched_noise_power"] / (SAMPLE_RATE * metrics["energy"])
|
||
row["carrier_noise_model"] = {str(carrier): matched_noise(frequencies, density, taps, carrier)
|
||
for carrier in (50000, 250000, 500000, 750000, 1000000)}
|
||
rows.append(row)
|
||
return {"window_seconds": .1, "rrc": metrics, "windows": rows,
|
||
"summary": {key: spread([r[key] for r in rows]) for key in (
|
||
"mean_i", "mean_q", "ac_power_code2", "band_27000_power", "band_675000_power",
|
||
"band_1000000_power", "band_27000_density", "band_675000_density",
|
||
"matched_noise_power", "equivalent_n0_code2_per_hz")}}
|
||
|
||
|
||
def analyze_rx(tx: Cs8File, rx: Cs8File, starts: list[int], frame_samples: int,
|
||
samples_per_symbol: int) -> dict:
|
||
rows = []
|
||
taps, _ = filter_metrics(samples_per_symbol)
|
||
for frame, start in enumerate(starts):
|
||
reference = tx.read(frame * frame_samples, frame_samples)
|
||
trim = 21 * samples_per_symbol
|
||
training_count = (frame_samples - 2 * trim) // 8
|
||
margin = samples_per_symbol
|
||
# Поиск остаётся внутри итоговой обучающей плитки при любом lag±margin.
|
||
if training_count <= 4 * margin:
|
||
raise ValueError("Обучающая плитка слишком коротка для защищённого поиска")
|
||
train_reference = reference[trim + 2 * margin:trim + training_count - 2 * margin]
|
||
train_received = rx.read(start + trim + margin, training_count - 2 * margin)
|
||
offset = align_training(train_reference, train_received, margin)
|
||
start += offset
|
||
received = rx.read(start, frame_samples)
|
||
active_reference = reference[trim:-trim]
|
||
active_received = received[trim:-trim]
|
||
# Задержка обучена на плитке 0: только чётные плитки обучают основную модель.
|
||
fitted = fit_channel(active_reference, active_received, SAMPLE_RATE, spectral_taps=taps)
|
||
moments = Moments()
|
||
moments.add(active_received)
|
||
row = {"frame": frame + 1, "start": start, "alignment_adjustment": offset,
|
||
"rx_active": moments.result(SAMPLE_RATE // samples_per_symbol), **fitted}
|
||
row["effective_signal_to_residual_db"] = float(10 * np.log10(
|
||
fitted["heldout_signal_power"] / fitted["heldout_residual_power"]))
|
||
rows.append(row)
|
||
print(f"RX frame {frame + 1}/{len(starts)}", flush=True)
|
||
return {"frames": rows, "summary": {key: spread([r[key] for r in rows]) for key in (
|
||
"alignment_adjustment", "cfo_hz", "gain_abs", "heldout_signal_power",
|
||
"train_residual_power", "heldout_residual_power", "heldout_actual_ac_power",
|
||
"heldout_cross_power", "effective_signal_to_residual_db", "heldout_residual_band_27000_power",
|
||
"heldout_residual_band_675000_power", "heldout_residual_band_1000000_power",
|
||
"heldout_residual_matched_power")}}
|
||
|
||
|
||
def audit(old_directory: Path = OLD_DIRECTORY, new_directory: Path = NEW_DIRECTORY) -> dict:
|
||
"""Все входы читаются; результат возвращается отдельно от исторических файлов."""
|
||
old_meta = read_json(old_directory / "tx_metadata.json")
|
||
new_meta = read_json(new_directory / "tx_metadata.json")
|
||
old_capture = read_json(old_directory / "capture_result.json")
|
||
official_path = PROJECT_ROOT / "data/processed/lab054/rf_summary.json"
|
||
official = read_json(official_path)
|
||
protected = sorted(set([p for directory in (old_directory, new_directory,
|
||
official_path.parent) for p in directory.iterdir() if p.is_file()]))
|
||
hashes = {str(path): sha256(path) for path in protected}
|
||
expected = {old_directory / "tx_stream.cs8": old_meta["tx_stream"]["sha256"],
|
||
old_directory / "rx_stream.cs8": old_capture["rx_sha256"],
|
||
new_directory / "tx_stream.cs8": new_meta["tx_sha256"],
|
||
new_directory / "rx_stream.cs8": official["artifact_sha256"]["rx_stream.cs8"]}
|
||
for path, expected_hash in expected.items():
|
||
if hashes[str(path)] != expected_hash:
|
||
raise ValueError(f"SHA не совпал: {path}")
|
||
if official["result"] != "NOT_A" or official["exact_packets"] != 0:
|
||
raise ValueError("Историческое состояние Lab054 изменилось; нужен новый аудит")
|
||
accepted = [r for r in read_json(old_directory / "independent_reprocess_result_corrected.json")["strict"]["outcomes"] if r["accepted"]]
|
||
old_starts = [r["refined_start_sample"] for r in accepted]
|
||
# Только диагностические исходные положения из завершённого анализа Lab054.
|
||
new_starts = [4161351, 4253606, 4345863, 4438119, 4530375, 4622631, 4714887, 4807143,
|
||
4899399, 4991655, 5083911, 5176167, 5268423, 5360679, 5452935, 5545191]
|
||
specifications = [("Lab047", old_directory, 20000, old_starts,
|
||
[.05, .20, .35, .50, 12.05, 12.4, 12.75, 13.1, 13.45, 13.8]),
|
||
("Lab054", new_directory, 500000, new_starts,
|
||
[.02, .12, .22, .32, 1., 1.25, 1.5, 1.75, 2., 2.5, 2.8])]
|
||
result = {"rf_tx": 0, "official_lab054": "NOT_A", "independent_receive": False,
|
||
"sample_rate_hz": SAMPLE_RATE, "carrier_hz": CARRIER,
|
||
"input_sha256": hashes, "experiments": {}}
|
||
for name, directory, symbol_rate, starts, background_seconds in specifications:
|
||
frame_samples = 5766 * (SAMPLE_RATE // symbol_rate)
|
||
tx, rx = Cs8File(directory / "tx_stream.cs8"), Cs8File(directory / "rx_stream.cs8")
|
||
print(name, "TX", flush=True)
|
||
tx_result = analyze_tx(directory / "tx_stream.cs8", frame_samples, len(starts), symbol_rate)
|
||
background = analyze_background(rx, background_seconds, SAMPLE_RATE // symbol_rate)
|
||
rx_result = analyze_rx(tx, rx, starts, frame_samples, SAMPLE_RATE // symbol_rate)
|
||
power = rx_result["summary"]["heldout_signal_power"]["median"]
|
||
n0 = background["summary"]["equivalent_n0_code2_per_hz"]["median"]
|
||
result["experiments"][name] = {"directory": str(directory), "symbol_rate": symbol_rate,
|
||
"frame_samples": frame_samples, "tx_bytes": tx.size * 2, "rx_bytes": rx.size * 2,
|
||
"tx": tx_result, "background": background, "rx": rx_result,
|
||
"conditional_es_n0_equivalent_db": float(10 * np.log10(power / symbol_rate / n0))}
|
||
old, new = result["experiments"].values()
|
||
ratio_db = lambda new_value, old_value: float(10 * np.log10(new_value / old_value))
|
||
result["changes_db"] = {
|
||
"tx_active_power": ratio_db(new["tx"]["active"]["power_code2"], old["tx"]["active"]["power_code2"]),
|
||
"rx_fitted_signal_power": ratio_db(new["rx"]["summary"]["heldout_signal_power"]["median"], old["rx"]["summary"]["heldout_signal_power"]["median"]),
|
||
"rate_penalty": ratio_db(new["symbol_rate"], old["symbol_rate"]),
|
||
"equivalent_n0": ratio_db(new["background"]["summary"]["equivalent_n0_code2_per_hz"]["median"], old["background"]["summary"]["equivalent_n0_code2_per_hz"]["median"])}
|
||
result["changes_db"]["conditional_es_n0_equivalent"] = new["conditional_es_n0_equivalent_db"] - old["conditional_es_n0_equivalent_db"]
|
||
changed = [str(path) for path in protected if sha256(path) != hashes[str(path)]]
|
||
if changed:
|
||
raise RuntimeError(f"Исходные артефакты изменились: {changed}")
|
||
result["source_integrity_after_analysis"] = True
|
||
return result
|
||
|
||
|
||
def run_functional_tests() -> None:
|
||
"""Минимальная проверка единиц и моментов без экспериментальных данных."""
|
||
moments = Moments()
|
||
moments.add(np.full(32, 3 + 4j))
|
||
assert moments.result(5)["power_code2"] == 25
|
||
assert moments.result(5)["symbol_energy_code2_seconds"] == 5
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description=__doc__)
|
||
parser.add_argument("--output", type=Path, default=PROJECT_ROOT / "data/processed/lab055/energy_audit.json")
|
||
arguments = parser.parse_args()
|
||
target = arguments.output.resolve()
|
||
allowed = (PROJECT_ROOT / "data/processed/lab055").resolve()
|
||
if not target.is_relative_to(allowed) or target.exists():
|
||
raise ValueError("Нужен новый файл только внутри data/processed/lab055")
|
||
result = audit()
|
||
target.parent.mkdir(parents=True, exist_ok=True)
|
||
with target.open("x", encoding="utf-8") as stream:
|
||
json.dump(result, stream, ensure_ascii=False, indent=2, allow_nan=False)
|
||
stream.write("\n")
|
||
print(json.dumps(result["changes_db"], indent=2))
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|