110 lines
5.1 KiB
Python
110 lines
5.1 KiB
Python
"""Только программные нарушения IQ; эталоны сюда не передаются приёмнику."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import json
|
||
import math
|
||
from pathlib import Path
|
||
import statistics
|
||
from typing import cast
|
||
|
||
import numpy as np
|
||
from scipy.optimize import brentq
|
||
from scipy.signal import lfilter
|
||
|
||
from protocol.video003_modem import SAMPLE_RATE_HZ, samples_per_symbol
|
||
from protocol.video005_iq_waveform import NOMINAL_CARRIER_HZ, SYMBOL_RATE
|
||
|
||
|
||
ROOT = Path(__file__).resolve().parents[1]
|
||
LAB066_RESULT = ROOT / "data/processed/lab066/rf_result.json"
|
||
|
||
|
||
def impair(
|
||
iq: np.ndarray,
|
||
*,
|
||
timing_offset_symbols: float = 0.0,
|
||
cfo_hz: float = 0.0,
|
||
phase_degrees: float = 0.0,
|
||
ebn0_db: float | None = None,
|
||
seed: int = 20260917,
|
||
symbol_rate: int = SYMBOL_RATE,
|
||
) -> np.ndarray:
|
||
"""Сдвиг времени, CFO, фаза и AWGN с фиксированной затравкой."""
|
||
if not all(math.isfinite(value) for value in (timing_offset_symbols, cfo_hz, phase_degrees)):
|
||
raise ValueError("параметры канала должны быть конечными")
|
||
if ebn0_db is not None and not math.isfinite(ebn0_db):
|
||
raise ValueError("Eb/N0 должно быть конечным")
|
||
transmitted = np.asarray(iq, dtype=np.complex64)
|
||
received = np.empty_like(transmitted)
|
||
delay = timing_offset_symbols * samples_per_symbol(symbol_rate)
|
||
# Ограничиваем рабочие массивы: полная видеозапись занимает миллионы IQ-отсчётов.
|
||
random_generator = np.random.default_rng(seed)
|
||
sigma = (math.sqrt(1 / (2 * (2 / 3) * 10 ** (ebn0_db / 10)) / 2)
|
||
if ebn0_db is not None else 0.0)
|
||
for start in range(0, transmitted.size, 1_000_000):
|
||
stop = min(start + 1_000_000, transmitted.size)
|
||
indexes = np.arange(start, stop, dtype=np.float64)
|
||
if delay:
|
||
positions = indexes - delay
|
||
lower = np.floor(positions).astype(np.int64)
|
||
valid = (lower >= 0) & (lower + 1 < transmitted.size)
|
||
part = np.zeros(stop - start, dtype=np.complex64)
|
||
fraction = positions[valid] - lower[valid]
|
||
part[valid] = (transmitted[lower[valid]] * (1 - fraction)
|
||
+ transmitted[lower[valid] + 1] * fraction)
|
||
else:
|
||
part = transmitted[start:stop].copy()
|
||
if cfo_hz or phase_degrees:
|
||
phase = math.radians(phase_degrees) + 2 * math.pi * cfo_hz * indexes / SAMPLE_RATE_HZ
|
||
part *= np.exp(1j * phase).astype(np.complex64)
|
||
if ebn0_db is not None:
|
||
part += (sigma * (random_generator.standard_normal(stop - start)
|
||
+ 1j * random_generator.standard_normal(stop - start))
|
||
).astype(np.complex64)
|
||
received[start:stop] = part
|
||
return received
|
||
|
||
|
||
def lab066_like(iq: np.ndarray, *, seed: int = 20260917) -> tuple[np.ndarray, dict]:
|
||
"""Диагностический суррогат по двум измеренным сводным характеристикам.
|
||
|
||
Частотная рябь и добавочный остаток не воспроизводят точную OTA-передачу.
|
||
"""
|
||
report = json.loads(LAB066_RESULT.read_text(encoding="utf-8"))
|
||
rows = report["physical_diagnostics"]["rows"]
|
||
ripple_db = statistics.median(row["band_gain_nonuniformity_p95_p05_db"] for row in rows)
|
||
useful_to_residual_db = statistics.mean(row["matched_residual_db"] for row in rows)
|
||
frequencies = np.linspace(NOMINAL_CARRIER_HZ - 337_500,
|
||
NOMINAL_CARRIER_HZ + 337_500, 2048)
|
||
delay_samples = 8
|
||
|
||
def modeled_ripple(alpha: float) -> float:
|
||
response = np.abs(1 + alpha * np.exp(-2j * np.pi * frequencies
|
||
* delay_samples / SAMPLE_RATE_HZ))
|
||
return 20 * math.log10(float(np.percentile(response, 95)
|
||
/ np.percentile(response, 5)))
|
||
|
||
alpha = cast(float, brentq(lambda value: modeled_ripple(value) - ripple_db,
|
||
0.0, 0.8, full_output=False))
|
||
coefficients = np.zeros(delay_samples + 1)
|
||
coefficients[0], coefficients[-1] = 1.0, alpha
|
||
filtered = cast(np.ndarray, lfilter(coefficients, [1.0],
|
||
np.asarray(iq, dtype=np.complex64)))
|
||
residual_ratio = 10 ** (useful_to_residual_db / 10)
|
||
sigma = math.sqrt(1 / (2 * residual_ratio))
|
||
random_generator = np.random.default_rng(seed)
|
||
result = filtered + sigma * (
|
||
random_generator.standard_normal(filtered.size)
|
||
+ 1j * random_generator.standard_normal(filtered.size)
|
||
)
|
||
return result.astype(np.complex64), {
|
||
"label": "LAB066-LIKE DIAGNOSTIC CHANNEL MODEL",
|
||
"measured_band_p95_p05_db": ripple_db,
|
||
"modeled_band_p95_p05_db": modeled_ripple(alpha),
|
||
"measured_matched_residual_db": useful_to_residual_db,
|
||
"filter_delay_samples": delay_samples,
|
||
"filter_second_tap": alpha,
|
||
"limitation": "FIR ripple and white residual surrogate do not reproduce OTA propagation",
|
||
}
|