245 lines
8.7 KiB
Python
245 lines
8.7 KiB
Python
from __future__ import annotations
|
|
|
|
import math
|
|
import re
|
|
from dataclasses import dataclass
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
|
|
|
|
HEADROOM_OK_DB = 6.0
|
|
HEADROOM_LIMIT_DB = 0.0
|
|
DEFAULT_HEADROOM_EPS_W = 1e-12
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class ExcursionSanityResult:
|
|
ok: bool
|
|
median_ratio: float
|
|
expected_ratio: float
|
|
relative_error: float
|
|
message: str
|
|
|
|
|
|
def headroom_status(headroom_db: float) -> str:
|
|
if not math.isfinite(headroom_db):
|
|
return "DISABLED"
|
|
if headroom_db <= HEADROOM_LIMIT_DB:
|
|
return "LIMIT EXCEEDED"
|
|
if headroom_db <= HEADROOM_OK_DB:
|
|
return "WARNING"
|
|
return "OK"
|
|
|
|
|
|
class WinISDLimitCurve:
|
|
"""Frequency-dependent WinISD maximum-power protection curve."""
|
|
|
|
def __init__(
|
|
self,
|
|
frequency_hz: np.ndarray,
|
|
pmax_w: np.ndarray,
|
|
source: Path | None = None,
|
|
eps_w: float = DEFAULT_HEADROOM_EPS_W,
|
|
):
|
|
if frequency_hz.ndim != 1 or pmax_w.ndim != 1:
|
|
raise ValueError("WinISD limit arrays must be one-dimensional")
|
|
if len(frequency_hz) != len(pmax_w):
|
|
raise ValueError("WinISD limit arrays must have the same length")
|
|
if len(frequency_hz) < 2:
|
|
raise ValueError("WinISD limit curve needs at least two points")
|
|
if eps_w <= 0.0:
|
|
raise ValueError("eps_w must be positive")
|
|
|
|
valid = np.isfinite(frequency_hz) & np.isfinite(pmax_w)
|
|
valid &= frequency_hz > 0.0
|
|
valid &= pmax_w > eps_w
|
|
if np.count_nonzero(valid) < 2:
|
|
raise ValueError("WinISD limit curve has fewer than two valid points")
|
|
|
|
frequency_hz = frequency_hz[valid].astype(np.float64)
|
|
pmax_w = pmax_w[valid].astype(np.float64)
|
|
order = np.argsort(frequency_hz)
|
|
frequency_hz = frequency_hz[order]
|
|
pmax_w = pmax_w[order]
|
|
|
|
unique_frequency, unique_indices = np.unique(frequency_hz, return_index=True)
|
|
self.frequency_hz = unique_frequency
|
|
self.pmax_w = pmax_w[unique_indices]
|
|
self.source = source
|
|
self.eps_w = float(eps_w)
|
|
|
|
@classmethod
|
|
def from_file(cls, path: str | Path, eps_w: float = DEFAULT_HEADROOM_EPS_W) -> "WinISDLimitCurve":
|
|
path = Path(path)
|
|
frequency_hz, values = load_two_column_curve(path)
|
|
return cls(frequency_hz, values, source=path, eps_w=eps_w)
|
|
|
|
def interpolate(self, frequencies_hz: np.ndarray) -> np.ndarray:
|
|
frequencies_hz = np.asarray(frequencies_hz, dtype=np.float64)
|
|
pmax = np.full_like(frequencies_hz, np.nan, dtype=np.float64)
|
|
valid = np.isfinite(frequencies_hz) & (frequencies_hz > 0.0)
|
|
if not np.any(valid):
|
|
return pmax
|
|
|
|
interp_f = np.clip(frequencies_hz[valid], self.frequency_hz[0], self.frequency_hz[-1])
|
|
pmax[valid] = np.interp(
|
|
np.log(interp_f),
|
|
np.log(self.frequency_hz),
|
|
self.pmax_w,
|
|
left=self.pmax_w[0],
|
|
right=self.pmax_w[-1],
|
|
)
|
|
return pmax
|
|
|
|
def compute_headroom(self, p_bin: np.ndarray, f_bins: np.ndarray) -> dict:
|
|
p_bin = np.asarray(p_bin, dtype=np.float64)
|
|
f_bins = np.asarray(f_bins, dtype=np.float64)
|
|
if p_bin.shape != f_bins.shape:
|
|
raise ValueError("p_bin and f_bins must have the same shape")
|
|
|
|
pmax = self.interpolate(f_bins)
|
|
protected_power = np.maximum(p_bin, 0.0)
|
|
valid = np.isfinite(pmax) & (pmax > self.eps_w) & np.isfinite(protected_power)
|
|
valid &= np.isfinite(f_bins) & (f_bins > 0.0)
|
|
|
|
utilization = np.zeros_like(protected_power, dtype=np.float64)
|
|
headroom = np.full_like(protected_power, np.inf, dtype=np.float64)
|
|
utilization[valid] = protected_power[valid] / pmax[valid]
|
|
headroom[valid] = 10.0 * np.log10(pmax[valid] / np.maximum(protected_power[valid], self.eps_w))
|
|
|
|
if not np.any(valid):
|
|
worst = float("nan")
|
|
critical_index = 0
|
|
critical_frequency = float("nan")
|
|
critical_power = float("nan")
|
|
critical_pmax = float("nan")
|
|
else:
|
|
valid_indices = np.flatnonzero(valid)
|
|
critical_index = int(valid_indices[np.argmin(headroom[valid])])
|
|
worst = float(headroom[critical_index])
|
|
critical_frequency = float(f_bins[critical_index])
|
|
critical_power = float(protected_power[critical_index])
|
|
critical_pmax = float(pmax[critical_index])
|
|
|
|
return {
|
|
"utilization_bin": utilization,
|
|
"headroom_db_bin": headroom,
|
|
"worst_headroom_db": worst,
|
|
"critical_frequency_hz": critical_frequency,
|
|
"critical_power_w": critical_power,
|
|
"critical_pmax_w": critical_pmax,
|
|
"over_limit_bool": bool(math.isfinite(worst) and worst <= HEADROOM_LIMIT_DB),
|
|
"status": headroom_status(worst),
|
|
}
|
|
|
|
|
|
def load_two_column_curve(path: str | Path) -> tuple[np.ndarray, np.ndarray]:
|
|
path = Path(path)
|
|
try:
|
|
text = path.read_text(encoding="utf-8-sig", errors="replace")
|
|
except OSError as exc:
|
|
raise ValueError(f"Could not read WinISD curve '{path}': {exc}") from exc
|
|
|
|
rows: list[tuple[float, float]] = []
|
|
for line in text.splitlines():
|
|
stripped = line.strip()
|
|
if not stripped or stripped.startswith(("#", "*")):
|
|
continue
|
|
values = _parse_numeric_line(stripped)
|
|
if len(values) >= 2:
|
|
rows.append((values[0], values[1]))
|
|
|
|
if not rows:
|
|
raise ValueError(f"No two-column numeric data found in WinISD curve '{path}'")
|
|
|
|
data = np.asarray(rows, dtype=np.float64)
|
|
return data[:, 0], data[:, 1]
|
|
|
|
|
|
def check_excursion_scaling(
|
|
excursion_10w_path: str | Path,
|
|
excursion_100w_path: str | Path,
|
|
tolerance: float = 0.10,
|
|
) -> ExcursionSanityResult:
|
|
f10, x10 = load_two_column_curve(excursion_10w_path)
|
|
f100, x100 = load_two_column_curve(excursion_100w_path)
|
|
valid_10 = np.isfinite(f10) & np.isfinite(x10) & (f10 > 0.0) & (x10 > 0.0)
|
|
valid_100 = np.isfinite(f100) & np.isfinite(x100) & (f100 > 0.0) & (x100 > 0.0)
|
|
f10 = f10[valid_10]
|
|
x10 = x10[valid_10]
|
|
f100 = f100[valid_100]
|
|
x100 = x100[valid_100]
|
|
if len(f10) < 2 or len(f100) < 2:
|
|
raise ValueError("Excursion sanity check needs at least two valid rows per file")
|
|
|
|
order = np.argsort(f10)
|
|
f10 = f10[order]
|
|
x10 = x10[order]
|
|
in_range = (f100 >= f10[0]) & (f100 <= f10[-1])
|
|
if np.count_nonzero(in_range) < 2:
|
|
raise ValueError("Excursion files do not overlap enough in frequency")
|
|
|
|
x10_interp = np.interp(np.log(f100[in_range]), np.log(f10), x10)
|
|
ratios = x100[in_range] / x10_interp
|
|
ratios = ratios[np.isfinite(ratios) & (ratios > 0.0)]
|
|
if len(ratios) == 0:
|
|
raise ValueError("Excursion sanity check produced no valid ratios")
|
|
|
|
expected = math.sqrt(10.0)
|
|
median_ratio = float(np.median(ratios))
|
|
relative_error = abs(median_ratio - expected) / expected
|
|
ok = relative_error <= tolerance
|
|
if ok:
|
|
message = (
|
|
f"Excursion scaling OK: median ratio {median_ratio:.3f}, "
|
|
f"expected {expected:.3f}."
|
|
)
|
|
else:
|
|
message = (
|
|
f"Excursion scaling warning: median ratio {median_ratio:.3f}, "
|
|
f"expected {expected:.3f}, error {relative_error * 100.0:.1f}%."
|
|
)
|
|
return ExcursionSanityResult(ok, median_ratio, expected, relative_error, message)
|
|
|
|
|
|
def run_self_test() -> None:
|
|
frequencies = np.array([20.0, 100.0, 1000.0], dtype=np.float64)
|
|
limit = WinISDLimitCurve(frequencies, np.full_like(frequencies, 100.0))
|
|
bins = np.array([20.0, 50.0, 100.0], dtype=np.float64)
|
|
|
|
cases = (
|
|
(50.0, 3.010299956639812, False),
|
|
(100.0, 0.0, True),
|
|
(200.0, -3.010299956639812, True),
|
|
)
|
|
for power, expected_headroom, expected_over_limit in cases:
|
|
p_bin = np.array([0.0, power, 0.0], dtype=np.float64)
|
|
result = limit.compute_headroom(p_bin, bins)
|
|
actual = result["worst_headroom_db"]
|
|
if not math.isclose(actual, expected_headroom, abs_tol=1e-6):
|
|
raise AssertionError(f"headroom for {power:g} W: expected {expected_headroom}, got {actual}")
|
|
if result["over_limit_bool"] is not expected_over_limit:
|
|
raise AssertionError(
|
|
f"over-limit for {power:g} W: expected {expected_over_limit}, got {result['over_limit_bool']}"
|
|
)
|
|
|
|
|
|
def _parse_numeric_line(line: str) -> list[float]:
|
|
if ";" in line:
|
|
parts = line.split(";")
|
|
else:
|
|
parts = re.split(r"[\t, ]+", line)
|
|
|
|
values: list[float] = []
|
|
for part in parts:
|
|
token = part.strip()
|
|
if not token:
|
|
continue
|
|
token = token.replace(",", ".")
|
|
try:
|
|
values.append(float(token))
|
|
except ValueError:
|
|
return []
|
|
return values
|