"""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()