Lab055: audit received signal energy across symbol rates
This commit is contained in:
LittleSam129
384
experiments/lab055_energy_audit.py
Normal file
384
experiments/lab055_energy_audit.py
Normal file
@@ -0,0 +1,384 @@
|
||||
"""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()
|
||||
Reference in New Issue
Block a user