import argparse
import math
import os
import sys
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.patheffects as pe_fx
import matplotlib.ticker as ticker
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.gridspec import GridSpec
from matplotlib.lines import Line2D
import numpy as np
plt.rcParams.update({
"font.family": "serif",
"font.serif": ["Times New Roman", "DejaVu Serif"],
"font.size": 9,
"axes.titlesize": 9,
"axes.labelsize": 9,
"legend.fontsize": 8,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"axes.linewidth": 0.8,
"lines.linewidth": 1.2,
"grid.linewidth": 0.4,
"grid.alpha": 0.4,
"axes.grid": True,
"figure.dpi": 150,
"savefig.dpi": 150,
"savefig.bbox": "tight",
"savefig.pad_inches": 0.04,
"text.usetex": False,
})
C_ADMISSIBLE = "#2ca02c"
C_BOUNDARY = "#ff7f0e"
C_VIOLATION = "#d62728"
C_DSFB = "#1f77b4"
C_COMPARATOR = "#9467bd"
C_NEUTRAL = "#7f7f7f"
C_HIGHLIGHT = "#e377c2"
C_GOLD = "#d4af37"
def save(fig, idx: int, name: str, out_dir: Path, dpi: int = 150):
out_dir.mkdir(parents=True, exist_ok=True)
stub = f"fig_{idx:02d}_{name}"
for ext in ("pdf", "png"):
path = out_dir / f"{stub}.{ext}"
fig.savefig(path, dpi=dpi, format=ext)
print(f" Saved fig {idx:02d}: {stub}.*")
plt.close(fig)
def ordinal_pattern_3(a, b, c):
if a <= b <= c: return 0
if a <= c < b: return 1
if c < a <= b: return 2
if b < a <= c: return 3
if b <= c < a: return 4
return 5
def permutation_entropy(x, normalised=True):
counts = np.zeros(6, dtype=float)
for i in range(len(x) - 2):
p = ordinal_pattern_3(x[i], x[i+1], x[i+2])
counts[p] += 1
total = counts.sum()
if total == 0: return 0.0
probs = counts[counts > 0] / total
h = -np.sum(probs * np.log2(probs))
return h / math.log2(6) if normalised else h
def allan_deviation(x, tau):
n = len(x)
if n < 2*tau + 1: return float('nan')
diffs = [x[k + 2*tau] - 2*x[k + tau] + x[k] for k in range(n - 2*tau)]
avar = np.mean(np.array(diffs)**2) / (2 * tau**2)
return math.sqrt(max(avar, 0))
def arrhenius(T_c, ea_ev, alpha0=1.0):
kB = 8.617333e-5 return alpha0 * math.exp(-ea_ev / (kB * (T_c + 273.15)))
def cosine_sim(a, b):
a, b = np.array(a), np.array(b)
return float(np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b) + 1e-30))
def rng(seed=42):
return np.random.default_rng(seed)
def plot_fig21(out_dir, dpi):
r = rng()
t = np.arange(256)
noise = r.standard_normal(256) * 0.5 + 0.3
periodic = np.tile([0.1, 0.2, 0.3, 0.2], 64)
drifting = 0.05 + t * 0.001 + r.standard_normal(256) * 0.02
def rolling_pe(x, w=16):
out = []
for i in range(w, len(x)):
out.append(permutation_entropy(x[i-w:i]))
return np.array(out)
def rolling_shannon_approx(x, w=16, bins=8):
out = []
for i in range(w, len(x)):
seg = x[i-w:i]
counts, _ = np.histogram(seg, bins=bins)
probs = counts / w
probs = probs[probs > 0]
h = -np.sum(probs * np.log2(probs))
out.append(h / math.log2(bins))
return np.array(out)
fig, axes = plt.subplots(2, 1, figsize=(6.5, 4.5), sharex=True)
ax1, ax2 = axes
k = np.arange(len(rolling_pe(noise)))
for sig, label, col in [(noise, "Wide-sense stationary noise", C_DSFB),
(periodic, "Period-3 (hidden determinism)", C_VIOLATION),
(drifting, "Drifting baseline", C_BOUNDARY)]:
ax1.plot(k, rolling_pe(sig), label=label, color=col, lw=1.2)
ax2.plot(k, rolling_shannon_approx(sig), color=col, lw=1.2, linestyle="--")
ax1.axhline(0.70, color=C_VIOLATION, lw=0.8, ls=":", alpha=0.8)
ax1.axhline(0.92, color=C_ADMISSIBLE, lw=0.8, ls=":", alpha=0.8)
ax1.annotate("HiddenDeterminism | NPE < 0.70", xy=(5, 0.67), fontsize=7, color=C_VIOLATION)
ax1.annotate("StochasticNoise | NPE ≥ 0.92", xy=(5, 0.94), fontsize=7, color=C_ADMISSIBLE)
ax1.set_ylabel("Normalised PE")
ax1.set_title("Fig 21. Permutation Entropy (Bandt & Pompe 2002) vs Shannon Entropy — Rolling W=16")
ax1.legend(loc="upper right", fontsize=7)
ax1.set_ylim(-0.05, 1.05)
ax2.set_ylabel("Normalised Shannon\n(8-bin histogram)")
ax2.set_xlabel("Observation index")
ax2.set_ylim(-0.05, 1.05)
ax2.set_title("Shannon Entropy (reference, dashed)")
plt.tight_layout()
save(fig, 21, "permutation_entropy_vs_shannon", out_dir, dpi)
def plot_fig22(out_dir, dpi):
r = rng()
n = 30
def rat(x):
a = sum(1 for i in range(len(x)) for j in range(i+1, len(x)) if x[i] > x[j])
e = n*(n-1)/4
var = n*(2*n+5)*(n-1)/72
return (a - e) / math.sqrt(var)
windows = {
"Strictly increasing trend": np.linspace(0.05, 0.20, n),
"Strictly decreasing trend": np.linspace(0.20, 0.05, n),
"Sinusoidal (stationary)": 0.10 + 0.02*np.sin(np.linspace(0, 4*np.pi, n)),
"WGN (stationary)": 0.10 + r.standard_normal(n) * 0.02,
}
fig, axes = plt.subplots(2, 2, figsize=(7, 4.5))
axes = axes.ravel()
colors = [C_VIOLATION, C_BOUNDARY, C_ADMISSIBLE, C_DSFB]
for ax, (label, x), col in zip(axes, windows.items(), colors):
z = rat(x)
sig = abs(z) > 1.96
ax.plot(x, color=col, lw=1.3)
ax.set_title(f"{label}\nZ = {z:.2f} {'[TREND]' if sig else '[ok]'}", fontsize=8)
ax.axhline(np.mean(x), color=C_NEUTRAL, lw=0.8, ls="--", alpha=0.7)
ax.set_xlabel("Sample index", fontsize=8)
ax.set_ylabel("Residual norm", fontsize=8)
if sig:
for sp in ax.spines.values():
sp.set_edgecolor(C_VIOLATION); sp.set_linewidth(2)
fig.suptitle("Fig 22. Reverse Arrangements Test (Olmstead-Tukey 1947)\n"
"WSS pre-condition check for GUM calibration windows", fontsize=9)
plt.tight_layout()
save(fig, 22, "reverse_arrangements_test", out_dir, dpi)
def plot_fig23(out_dir, dpi):
snr_db = np.linspace(-10, 40, 200)
gammas = 10**(snr_db / 10)
rho_floor = 1.0 / np.sqrt(gammas)
n_obs = 100
crlb_phase = 1.0 / (n_obs * gammas)
crlb_freq = 6.0 / (n_obs**3 * gammas * (2*np.pi)**2)
rho_values = [0.05, 0.10, 0.20, 0.40]
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
ax1, ax2 = axes
ax1.semilogy(snr_db, rho_floor, color=C_NEUTRAL, lw=1.5, ls="--", label=r"$\rho_{\rm floor} = 1/\sqrt{\gamma}$")
for rho, col in zip(rho_values, [C_ADMISSIBLE, C_DSFB, C_BOUNDARY, C_VIOLATION]):
margin = np.array([rho / rf for rf in rho_floor])
ax1.axhline(rho, color=col, lw=0.9, alpha=0.7, ls=":")
above_3x = snr_db[margin >= 3.0]
if len(above_3x) > 0:
ax1.axvline(above_3x[0], color=col, lw=0.8, alpha=0.5)
ax1.fill_betweenx([1e-3, 2e0], -10, 0, alpha=0.07, color=C_VIOLATION, label="CRLB alert zone")
ax1.set_xlabel("SNR (dB)")
ax1.set_ylabel(r"Radius $\rho$")
ax1.set_title("CRLB Physics Floor vs ρ\n(Kay 1993; Rife-Boorstyn 1974)")
ax1.legend(fontsize=7)
ax1.set_ylim(1e-3, 2.0)
ax2.semilogy(snr_db, crlb_phase, color=C_DSFB, lw=1.5, label=r"CRLB$_\phi$ (N=100)")
ax2.semilogy(snr_db, crlb_freq, color=C_BOUNDARY, lw=1.5, label=r"CRLB$_f$ (N=100)")
ax2.set_xlabel("SNR (dB)")
ax2.set_ylabel("CRLB [rad² or Hz²]")
ax2.set_title(r"CRLB$_\phi$ = 1/(Nγ) CRLB$_f$ = 6/(N³γ(2π)²)")
ax2.legend(fontsize=7)
fig.suptitle("Fig 23. Cramér–Rao Lower Bound Floor and Admissibility Margin", y=1.01)
plt.tight_layout()
save(fig, 23, "crlb_floor_vs_rho", out_dir, dpi)
def plot_fig24(out_dir, dpi):
T = np.linspace(25, 200, 300)
drift_gaas = np.array([arrhenius(t, 1.6) for t in T])
drift_gan = np.array([arrhenius(t, 2.1) for t in T])
af_gaas = drift_gaas / drift_gaas[0]
af_gan = drift_gan / drift_gan[0]
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
ax1, ax2 = axes
ax1.semilogy(T, drift_gaas, color=C_VIOLATION, lw=1.5, label="GaAs pHEMT $E_a$ = 1.6 eV")
ax1.semilogy(T, drift_gan, color=C_DSFB, lw=1.5, label="GaN HEMT $E_a$ = 2.1 eV")
ax1.axvline(125, color=C_NEUTRAL, lw=0.8, ls="--", alpha=0.7)
ax1.annotate("125°C\noperating", xy=(125, drift_gaas[T>124][0]*1.3), fontsize=7, color=C_NEUTRAL)
ax1.set_xlabel("Junction temperature T (°C)")
ax1.set_ylabel("Normalised drift rate $k(T)$")
ax1.set_title("Arrhenius drift rate\n(Kayali 1999, JPL-96-25)")
ax1.legend(fontsize=7)
ax2.semilogy(T, af_gaas, color=C_VIOLATION, lw=1.5, label="GaAs pHEMT AF")
ax2.semilogy(T, af_gan, color=C_DSFB, lw=1.5, label="GaN HEMT AF")
ax2.axhline(10, color=C_NEUTRAL, lw=0.8, ls=":", alpha=0.7)
ax2.annotate("10× acceleration", xy=(60, 12), fontsize=7, color=C_NEUTRAL)
ax2.set_xlabel("Junction temperature T (°C)")
ax2.set_ylabel("Acceleration factor AF(T) = k(T)/k(25°C)")
ax2.set_title("Acceleration factor vs temperature")
ax2.legend(fontsize=7)
fig.suptitle("Fig 24. Arrhenius Thermal Model — PA Drift Rate (Physics-of-Failure)", y=1.01)
plt.tight_layout()
save(fig, 24, "arrhenius_pa_drift", out_dir, dpi)
def plot_fig25(out_dir, dpi):
r = rng()
n = 200
tau = 2
noise_x = r.standard_normal(n) * 0.15 + 0.3
t = np.arange(n)
periodic_x = 0.3 + 0.1*np.sin(2*np.pi*t/8) + r.standard_normal(n)*0.01
drift_x = 0.05 + t*0.001 + r.standard_normal(n)*0.02
def embed(x, tau):
return x[:-tau], x[tau:]
fig, axes = plt.subplots(1, 3, figsize=(7.5, 3), sharex=False, sharey=False)
configs = [
(noise_x, "Stochastic noise\n(D₂ ≈ 2.0)", C_DSFB),
(periodic_x, "Structured orbit\n(D₂ ≈ 0.9)", C_ADMISSIBLE),
(drift_x, "Drifting attractor\n(D₂ ≈ 0.3)", C_BOUNDARY),
]
for ax, (x, lab, col) in zip(axes, configs):
xd, xd2 = embed(x, tau)
sc = ax.scatter(xd, xd2, c=col, s=3, alpha=0.5)
ax.set_title(lab, fontsize=8)
ax.set_xlabel(r"$\|r(k)\|$", fontsize=8)
ax.set_ylabel(r"$\|r(k{-}\tau)\|$", fontsize=8)
ax.set_aspect("equal", adjustable="datalim")
fig.suptitle("Fig 25. Delay-Coordinate Phase Portrait (Takens 1981) — τ = 2 samples", y=1.02)
plt.tight_layout()
save(fig, 25, "delay_embedding_phase_portrait", out_dir, dpi)
def plot_fig26(out_dir, dpi):
r = rng()
n = 512; tau = 2
def gp_d2(x, radii):
xd, xd2 = x[:-tau], x[tau:]
pts = np.column_stack([xd, xd2])
n_pts = len(pts)
c = []
for rad in radii:
pairs = 0
for i in range(n_pts):
d = np.sqrt(np.sum((pts - pts[i])**2, axis=1))
pairs += np.sum(d < rad) - 1 c.append(pairs / (n_pts * (n_pts-1)))
return np.array(c)
radii = np.logspace(-2, 0, 40)
noise = r.standard_normal(n)*0.3 + 0.5
periodic = 0.5 + 0.2*np.sin(2*np.pi*np.arange(n)/7) + r.standard_normal(n)*0.005
c_noise = gp_d2(noise, radii)
c_per = gp_d2(periodic, radii)
def slope(r, c):
valid = (c > 0.01) & (c < 0.9)
if valid.sum() < 5: return float('nan')
return np.polyfit(np.log(r[valid]), np.log(c[valid]), 1)[0]
d2_noise = slope(radii, c_noise)
d2_per = slope(radii, c_per)
fig, ax = plt.subplots(figsize=(5.5, 3.5))
ax.loglog(radii, c_noise, color=C_DSFB, lw=1.3, label=f"WGN D₂ est. = {d2_noise:.2f}")
ax.loglog(radii, c_per, color=C_ADMISSIBLE, lw=1.3, label=f"Periodic D₂ est. = {d2_per:.2f}")
ax.set_xlabel("Radius $r$")
ax.set_ylabel("Correlation integral $C(r)$")
ax.set_title("Fig 26. Grassberger–Procaccia Correlation Dimension D₂\n(Physica D, 1983)")
ax.legend(fontsize=7)
plt.tight_layout()
save(fig, 26, "correlation_dimension_d2", out_dir, dpi)
def plot_fig27(out_dir, dpi):
r = rng()
def birth_death_1d(x, max_radius=0.5, n_radii=60):
n = len(x)
idx = np.argsort(x)
xs = x[idx]
parent = list(range(n))
birth = np.zeros(n)
death = np.full(n, fill_value=max_radius * 1.5)
def find(a):
while parent[a] != a: a = parent[a]
return a
dists_all = []
for i in range(n):
for j in range(i+1, n):
dists_all.append((abs(x[i]-x[j]), i, j))
dists_all.sort()
for d, i, j in dists_all:
ri, rj = find(i), find(j)
if ri != rj:
dying = max(ri, rj)
death[dying] = d
parent[ri] = rj
return birth, death
n = 50
noise_x = r.standard_normal(n)*0.1 + 0.3
cluster_x = np.concatenate([r.standard_normal(25)*0.01 + 0.1,
r.standard_normal(25)*0.01 + 0.4])
b_n, d_n = birth_death_1d(noise_x)
b_c, d_c = birth_death_1d(cluster_x)
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
ax1, ax2 = axes
for b, d, col, lab, ax in [
(b_n, d_n, C_DSFB, "WGN — high Betti₀", ax1),
(b_c, d_c, C_ADMISSIBLE, "2-cluster signal — low Betti₀", ax2),
]:
finite = d < 1.4
ax.scatter(b[finite], d[finite], c=col, s=15, zorder=3, label=f"{lab}")
lim = max(d[finite].max() if finite.sum() > 0 else 0.5, 0.5)
ax.plot([0, lim], [0, lim], 'k--', lw=0.8, alpha=0.4, label="birth=death")
ax.set_xlabel("Birth radius")
ax.set_ylabel("Death radius")
ax.set_title(lab, fontsize=8)
betti0 = np.sum(d == d.max())
ax.annotate(f"n_finite = {finite.sum()}", xy=(0.02, lim*0.82), fontsize=7)
ax.legend(fontsize=7)
fig.suptitle("Fig 27. TDA Persistence Diagram — Betti₀ Birth/Death (Edelsbrunner 2002)", y=1.02)
plt.tight_layout()
save(fig, 27, "tda_persistence_landscape", out_dir, dpi)
def plot_fig28(out_dir, dpi):
r = rng()
n = 40
radii = np.linspace(0, 0.6, 100)
def betti0(x, rad):
n = len(x)
parent = list(range(n))
def find(a):
while parent[a] != a: a = parent[a]
return a
def union(a, b):
ra, rb = find(a), find(b)
if ra != rb: parent[ra] = rb; return True
return False
dists = [(abs(x[i]-x[j]), i, j) for i in range(n) for j in range(i+1,n)]
dists.sort()
for d, i, j in dists:
if d <= rad: union(i, j)
return len({find(i) for i in range(n)})
env_wgn = r.standard_normal(n)*0.15 + 0.3
env_fhss = np.concatenate([np.full(10, 0.1), np.full(10, 0.4),
np.full(10, 0.1), np.full(10, 0.4)])
env_jammer = np.concatenate([np.full(20, 0.1), np.full(20, 0.45)])
configs = [
(env_wgn, "WGN (uniform spread)", C_DSFB),
(env_fhss, "FHSS (2 hop states)", C_BOUNDARY),
(env_jammer, "Jammer onset (2 levels)", C_VIOLATION),
]
fig, ax = plt.subplots(figsize=(6, 3.5))
for x, lab, col in configs:
b0s = [betti0(x, rad) for rad in radii]
ax.step(radii, b0s, color=col, lw=1.3, label=lab, where='post')
ax.set_xlabel("Filtration radius ε")
ax.set_ylabel("Betti₀ (connected components)")
ax.set_title("Fig 28. Betti₀ vs Filtration Radius — Topological Phase Transitions\n"
"Edelsbrunner et al. 2002; Bubenik 2015")
ax.legend(fontsize=7)
ax.set_xlim(0, 0.6)
ax.set_ylim(0, n+2)
plt.tight_layout()
save(fig, 28, "betti0_birth_death_events", out_dir, dpi)
def plot_fig29(out_dir, dpi):
r = rng()
n_total = 1000
admissible_time = 800
entropy = np.concatenate([
r.standard_normal(admissible_time)*0.01 + 0.05, np.linspace(0.05, 0.80, 100) + r.standard_normal(100)*0.02, r.standard_normal(100)*0.02 + 0.85, ])
threshold = 0.05
urgency = np.concatenate([np.zeros(900), np.ones(100)]) emit_flags = np.zeros(n_total, dtype=bool)
last = -999.0
for i in range(n_total):
delta = abs(entropy[i] - last)
if last < 0 or delta >= threshold or urgency[i] > 0.8:
emit_flags[i] = True
last = entropy[i]
emit_cumsum = np.cumsum(emit_flags)
naive_cumsum = np.arange(1, n_total+1)
efficiency = 1.0 - emit_cumsum / naive_cumsum
fig, axes = plt.subplots(3, 1, figsize=(6.5, 5), sharex=True)
ax1, ax2, ax3 = axes
t = np.arange(n_total)
ax1.plot(t, entropy, color=C_DSFB, lw=0.7)
ax1.axvline(800, color=C_VIOLATION, lw=0.8, ls="--")
ax1.axvline(900, color=C_BOUNDARY, lw=0.8, ls="--")
ax1.set_ylabel("Grammar\nentropy h")
ax1.set_title("Fig 29. Pragmatic Information Gating (Atlan & Cohen 1998)\n"
"SOSA Backplane Efficiency — >99% Admissible Suppression")
ax1.annotate("State\ntransition", xy=(800, 0.5), fontsize=7, color=C_VIOLATION)
emitted_t = t[emit_flags]
ax2.eventplot(emitted_t, lineoffsets=0.5, linelengths=0.8, color=C_HIGHLIGHT)
ax2.set_yticks([])
ax2.set_ylabel("Emitted\nevents")
ax3.plot(t, efficiency * 100, color=C_ADMISSIBLE, lw=1.2)
ax3.axhline(99, color=C_NEUTRAL, ls=":", lw=0.8)
ax3.set_ylabel("Backplane\neff. (%)")
ax3.set_xlabel("Observation index")
ax3.set_ylim(50, 101)
admissible_eff = efficiency[800] * 100
ax3.annotate(f"{admissible_eff:.1f}% @ time 800", xy=(750, admissible_eff - 4),
fontsize=7, color=C_ADMISSIBLE)
plt.tight_layout()
save(fig, 29, "pragmatic_gating_sosa_efficiency", out_dir, dpi)
def plot_fig30(out_dir, dpi):
r = rng()
taus = np.array([1, 2, 4, 8, 16, 32, 64, 128])
def avar_profile(h_white, h_flicker, h_rw):
avar = h_white/(2*taus) + h_flicker*2*math.log(2) + h_rw*(2*math.pi**2/3)*taus
return np.sqrt(avar)
ocxo = avar_profile(1e-9, 1e-11, 1e-17)
tcxo = avar_profile(1e-7, 1e-9, 1e-15)
mems = avar_profile(1e-5, 1e-7, 1e-13)
spoofed = avar_profile(1e-5, 1e-8, 1e-17)
fig, axes = plt.subplots(1, 2, figsize=(7.5, 3.5))
ax1, ax2 = axes
for sig, col, lab in [
(ocxo, C_ADMISSIBLE, "OCXO Class A"),
(tcxo, C_DSFB, "TCXO Grade B"),
(mems, C_BOUNDARY, "MEMS resonator"),
(spoofed, C_VIOLATION, "Spoofed clock"),
]:
ax1.loglog(taus, sig, marker='o', ms=4, color=col, lw=1.2, label=lab)
ax1.set_xlabel("Averaging time τ (samples)")
ax1.set_ylabel(r"Allan deviation $\sigma_y(\tau)$")
ax1.set_title("Oscillator class fingerprints\n(Allan 1966; IEEE Std 1139-2008)")
ax1.legend(fontsize=7)
sigs = np.array([ocxo, tcxo, mems, spoofed])
labels = ["OCXO", "TCXO", "MEMS", "Spoofed"]
n = len(sigs)
sim = np.zeros((n, n))
for i in range(n):
for j in range(n):
sim[i, j] = cosine_sim(sigs[i], sigs[j])
cmap = LinearSegmentedColormap.from_list("cosine", [C_VIOLATION, "#f7f7f7", C_ADMISSIBLE])
im = ax2.imshow(sim, vmin=0.5, vmax=1.0, cmap=cmap)
ax2.set_xticks(range(n)); ax2.set_xticklabels(labels, fontsize=8)
ax2.set_yticks(range(n)); ax2.set_yticklabels(labels, fontsize=8)
for i in range(n):
for j in range(n):
ax2.text(j, i, f"{sim[i,j]:.2f}", ha='center', va='center',
color='k' if sim[i,j] > 0.7 else 'w', fontsize=8)
plt.colorbar(im, ax=ax2, fraction=0.046)
ax2.set_title("Cosine similarity matrix\n(threshold = 0.95)")
ax2.axhline(2.5, color="k", lw=1.5)
ax2.axvline(2.5, color="k", lw=1.5)
fig.suptitle("Fig 30. Hardware DNA Fingerprinting via Allan Variance (Physical-Layer Auth.)", y=1.02)
plt.tight_layout()
save(fig, 30, "hardware_dna_allan_fingerprints", out_dir, dpi)
def plot_fig31(out_dir, dpi):
N_vals = np.logspace(1, 4, 200)
rho = 0.10
fig, ax = plt.subplots(figsize=(5.5, 3.5))
for snr_db, col, ls in [(-5, C_VIOLATION, "--"), (0, C_BOUNDARY, "-."),
(10, C_DSFB, "-"), (20, C_ADMISSIBLE, "-")]:
gamma = 10**(snr_db/10)
rho_floor = 1/math.sqrt(gamma)
margin = rho / rho_floor * np.ones_like(N_vals) rho_crlb = np.sqrt(1.0 / (N_vals * gamma))
m2 = rho / rho_crlb
ax.loglog(N_vals, m2, color=col, lw=1.2, ls=ls, label=f"SNR = {snr_db} dB")
ax.axhline(3.0, color=C_NEUTRAL, lw=1.0, ls=":", alpha=0.8)
ax.fill_between([10, 1e4], 0, 3, alpha=0.07, color=C_VIOLATION)
ax.annotate("CRLB alert (margin < 3×)", xy=(20, 2.4), fontsize=7, color=C_VIOLATION)
ax.set_xlabel("Number of calibration observations N")
ax.set_ylabel(r"Margin factor $\rho / \rho_{\rm CRLB}$")
ax.set_title("Fig 31. CRLB Admissibility Margin vs N\n(Increases as √N — more data tightens bound)")
ax.legend(fontsize=7)
ax.set_ylim(0.1, 500)
plt.tight_layout()
save(fig, 31, "crlb_margin_vs_n_obs", out_dir, dpi)
def plot_fig32(out_dir, dpi):
r = rng()
n = 256; tau = 2
t = np.arange(n)
noise_s = r.standard_normal(n)*0.1 + 0.3
period_s = 0.3 + 0.1*np.sin(2*np.pi*t/11) + 0.05*np.sin(2*np.pi*t/7)
transient = 0.3 + np.exp(-t/40)*0.2*np.sin(2*np.pi*t/9) + r.standard_normal(n)*0.01
def koopman_proxy(x, window=32):
out = []
for i in range(window, len(x)):
seg = x[i-window:i]
m = seg.mean()
v = seg.var()
out.append(v / m if m > 1e-6 else 0.0)
return np.array(out)
fig, axes = plt.subplots(2, 3, figsize=(8, 4.5))
configs = [
(noise_s, "WGN", C_DSFB),
(period_s, "Multi-tone", C_ADMISSIBLE),
(transient,"Transient", C_BOUNDARY),
]
for j, (x, lab, col) in enumerate(configs):
ax_ts = axes[0, j]
ax_kop = axes[1, j]
ax_ts.plot(t, x, color=col, lw=0.8, alpha=0.8)
ax_ts.set_title(lab, fontsize=8)
ax_ts.set_ylabel("‖r(k)‖", fontsize=8)
kp = koopman_proxy(x)
ax_kop.plot(np.arange(len(kp)), kp, color=col, lw=1.0)
ax_kop.axhline(1.0, color=C_NEUTRAL, lw=0.7, ls=":")
ax_kop.set_xlabel("Obs. index", fontsize=8)
ax_kop.set_ylabel("VM ratio", fontsize=8)
ax_kop.annotate(f"mean={kp.mean():.3f}", xy=(5, kp.max()*0.9), fontsize=7)
fig.suptitle("Fig 32. Koopman Proxy (Variance-to-Mean Ratio) — Stochastic vs Structured Modes\n"
"Mezić 2005; high VM → stochastic; low VM → structured Koopman modes")
plt.tight_layout()
save(fig, 32, "koopman_mode_proxy", out_dir, dpi)
def plot_fig33(out_dir, dpi):
r = rng()
norms = np.concatenate([
r.standard_normal(100)*0.01 + 0.05, np.linspace(0.05, 0.20, 100) + r.standard_normal(100)*0.005
])
norms = norms.clip(0, None)
scale = 65536.0
raw = np.round(norms * scale).astype(np.int32)
raw = raw.clip(-(2**31), 2**31 - 1)
deq = raw.astype(np.float64) / scale
error = deq - norms
bound = 2**-14
fig, axes = plt.subplots(2, 1, figsize=(6.5, 4), sharex=True)
ax1, ax2 = axes
t = np.arange(len(norms))
ax1.plot(t, norms, color=C_DSFB, lw=0.9, label="f32 original")
ax1.plot(t, deq, color=C_BOUNDARY, lw=0.7, ls="--", label="Q16.16 round-trip", alpha=0.8)
ax1.axvline(100, color=C_NEUTRAL, lw=0.8, ls=":")
ax1.annotate("Drift onset", xy=(101, 0.06), fontsize=7, color=C_NEUTRAL)
ax1.set_title("Fig 33. Bit-Exactness: f32 vs Q16.16 Round-Trip\n"
"(Residual norm quantisation error < 2⁻¹⁴)")
ax1.set_ylabel("‖r(k)‖")
ax1.legend(fontsize=7)
ax2.plot(t, np.abs(error), color=C_VIOLATION, lw=0.7, label="|error|")
ax2.axhline(bound, color=C_ADMISSIBLE, lw=0.9, ls="--", label=f"bound = 2⁻¹⁴ ≈ {bound:.2e}")
ax2.axhline(2**-16 * 0.5, color=C_NEUTRAL, lw=0.7, ls=":", label="Q16.16 resolution/2")
ax2.set_ylabel("|Q error|")
ax2.set_xlabel("Observation index")
ax2.legend(fontsize=7)
ax2.set_ylim(0, bound * 3)
pct_below = np.mean(np.abs(error) < bound) * 100
ax2.annotate(f"{pct_below:.1f}% samples below bound", xy=(5, bound*2.2), fontsize=7)
plt.tight_layout()
save(fig, 33, "bit_exactness_q16_error", out_dir, dpi)
def plot_fig34(out_dir, dpi):
taus = np.logspace(0, 3, 200)
def avar_full(h_white, h_flicker, h_rw, taus):
avar = h_white/(2*taus) + h_flicker*2*math.log(2) + h_rw*(2*math.pi**2/3)*taus
return np.sqrt(avar)
profiles = [
("OCXO Class A", 1e-9, 1e-11, 1e-17, C_ADMISSIBLE),
("TCXO Grade B", 1e-7, 1e-9, 1e-15, C_DSFB),
("MEMS resonator",1e-5, 1e-7, 1e-13, C_BOUNDARY),
]
fig, ax = plt.subplots(figsize=(6, 3.5))
for label, hw, hf, hrw, col in profiles:
sigma = avar_full(hw, hf, hrw, taus)
ax.loglog(taus, sigma, color=col, lw=1.4, label=label)
min_idx = np.argmin(sigma)
ax.scatter(taus[min_idx], sigma[min_idx], color=col, s=30, zorder=4)
slope_m1 = taus**-0.5 * 1e-4 slope_p1 = taus**0.5 * 3e-8 ax.loglog(taus, slope_m1, 'k:', lw=0.7, alpha=0.5, label=r"$\tau^{-1/2}$ (white FM)")
ax.loglog(taus, slope_p1, 'k--', lw=0.7, alpha=0.5, label=r"$\tau^{+1/2}$ (RW FM)")
ax.set_xlabel("Averaging time τ (samples)")
ax.set_ylabel(r"Allan deviation $\sigma_y(\tau)$")
ax.set_title("Fig 34. Allan Deviation — Flicker-Noise Regimes\n"
"(Allan 1966; IEEE Std 1193-2003; h-coefficient model)")
ax.legend(fontsize=7)
plt.tight_layout()
save(fig, 34, "allan_deviation_oscillator_classes", out_dir, dpi)
def plot_fig35(out_dir, dpi):
r = rng()
n = 500
t = np.arange(n)
jammer_on = 0.12 + 0.03*np.sin(2*np.pi*t/11)
noise_only = 0.12 + r.standard_normal(n)*0.015
modes = [
(noise_only, "No jammer (WGN)", C_DSFB),
(jammer_on + r.standard_normal(n)*0.01, "Cyclostationary jammer (11-sample period)", C_VIOLATION),
(jammer_on + r.standard_normal(n)*0.005, "Jammer + low noise (clear hidden det.)", C_BOUNDARY),
]
fig, axes = plt.subplots(2, 1, figsize=(6.5, 4.5), sharex=True)
ax1, ax2 = axes
w = 32
for x, lab, col in modes:
pe_roll = [permutation_entropy(x[i-w:i]) for i in range(w, n)]
ax2.plot(np.arange(len(pe_roll)), pe_roll, color=col, lw=0.9, label=lab)
ax1.plot(t, noise_only, color=C_DSFB, lw=0.7, alpha=0.7, label="WGN only")
ax1.plot(t, jammer_on + r.standard_normal(n)*0.01, color=C_VIOLATION, lw=0.7, alpha=0.7,
label="Jammer embedded")
ax1.set_title("Fig 35. Permutation Entropy — Cyclostationary Jammer Detection\n"
"(Bandt & Pompe 2002 — hidden determinism in RF streams)")
ax1.set_ylabel("‖r(k)‖")
ax1.legend(fontsize=7)
ax2.axhline(0.70, color=C_VIOLATION, ls=":", lw=0.8)
ax2.annotate("HiddenDeterminism threshold", xy=(5, 0.67), fontsize=7, color=C_VIOLATION)
ax2.set_ylabel("Normalised PE")
ax2.set_xlabel("Observation index")
ax2.set_ylim(0, 1.05)
ax2.legend(fontsize=7)
plt.tight_layout()
save(fig, 35, "pe_cyclostationary_jammer", out_dir, dpi)
def plot_fig36(out_dir, dpi):
n = 2000
t = np.arange(n)
events = [200, 450, 900, 1300, 1700]
naive_msgs = np.arange(1, n+1, dtype=float)
pragmatic_msgs = np.zeros(n)
count = 0
for i in range(n):
if any(abs(i - ev) < 2 for ev in events) or i % 50 == 0:
count += 1
pragmatic_msgs[i] = count
fig, axes = plt.subplots(2, 1, figsize=(6.5, 4.5), sharex=True)
ax1, ax2 = axes
ax1.plot(t, naive_msgs / n * 100, color=C_NEUTRAL, lw=1.0, label="Naive streaming")
ax1.plot(t, pragmatic_msgs / n * 100, color=C_DSFB, lw=1.3, label="DSFB pragmatic gate")
for ev in events:
ax1.axvline(ev, color=C_VIOLATION, lw=0.7, ls=":")
ax1.set_ylabel("Cumulative messages\n(% of samples)")
ax1.set_title("Fig 36. SOSA Backplane Utilisation: Event-Centric vs Naive Streaming\n"
"Pragmatic gate suppresses >97% of Admissible-state heartbeats")
ax1.legend(fontsize=7)
reduction = 1 - pragmatic_msgs / naive_msgs
ax2.plot(t, reduction * 100, color=C_ADMISSIBLE, lw=1.0)
ax2.set_ylabel("Backplane savings (%)")
ax2.set_xlabel("Observation index")
final_savings = reduction[-1] * 100
ax2.annotate(f"Final savings: {final_savings:.1f}%", xy=(n*0.7, final_savings - 5),
fontsize=8, color=C_ADMISSIBLE)
plt.tight_layout()
save(fig, 36, "sosa_backplane_event_vs_stream", out_dir, dpi)
def plot_fig37(out_dir, dpi):
r = rng()
taus = np.array([1, 2, 4, 8, 16, 32, 64, 128], dtype=float)
def avar_profile(h_white, h_flicker, h_rw):
avar = h_white/(2*taus) + h_flicker*2*math.log(2) + h_rw*(2*math.pi**2/3)*taus
return np.sqrt(avar)
registered = avar_profile(1e-7, 1e-9, 1e-15)
genuine_attempts = [registered * (1 + r.standard_normal(8)*0.015) for _ in range(20)]
mems = avar_profile(1e-5, 1e-7, 1e-13)
spoofed_attempts = [mems * (1 + r.standard_normal(8)*0.01) for _ in range(20)]
genuine_sims = [cosine_sim(x, registered) for x in genuine_attempts]
spoofed_sims = [cosine_sim(x, registered) for x in spoofed_attempts]
fig, axes = plt.subplots(1, 2, figsize=(7.5, 3.5))
ax1, ax2 = axes
ax1.loglog(taus, registered, color=C_DSFB, lw=2, label="Registered DNA", zorder=5)
for x in genuine_attempts[:5]:
ax1.loglog(taus, x, color=C_ADMISSIBLE, lw=0.7, alpha=0.5)
for x in spoofed_attempts[:5]:
ax1.loglog(taus, x, color=C_VIOLATION, lw=0.7, alpha=0.5)
ax1.loglog([], [], color=C_ADMISSIBLE, lw=1, label="Genuine attempts (×5)")
ax1.loglog([], [], color=C_VIOLATION, lw=1, label="Spoofed attempts (×5)")
ax1.set_xlabel("τ (samples)"); ax1.set_ylabel(r"$\sigma_y(\tau)$")
ax1.set_title("σ_y(τ) fingerprints", fontsize=8)
ax1.legend(fontsize=7)
ax2.scatter([1]*20, genuine_sims, color=C_ADMISSIBLE, s=30, label="Genuine", zorder=3)
ax2.scatter([2]*20, spoofed_sims, color=C_VIOLATION, s=30, label="Spoofed", zorder=3)
ax2.axhline(0.95, color=C_NEUTRAL, lw=0.9, ls="--", label="Auth. threshold = 0.95")
ax2.set_xticks([1, 2]); ax2.set_xticklabels(["Genuine", "Spoofed"])
ax2.set_ylabel("Cosine similarity vs registered DNA")
ax2.set_title("Authentication distribution", fontsize=8)
ax2.legend(fontsize=7)
ax2.set_ylim(0.5, 1.05)
fig.suptitle("Fig 37. Hardware DNA Authentication — Genuine vs Spoofed Clock Fingerprint", y=1.02)
plt.tight_layout()
save(fig, 37, "hardware_dna_auth_genuine_spoofed", out_dir, dpi)
def plot_fig38(out_dir, dpi):
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlim(0, 10); ax.set_ylim(0, 7)
ax.axis("off")
ax.set_title("Fig 38. DSFB Information-Physics Architecture\n"
"(Phase-4: TDA + Attractor + Pragmatic Gate + Hardware DNA)", fontsize=9)
def box(ax, x, y, w, h, label, sublabel="", color=C_DSFB, alpha=0.18):
rect = mpatches.FancyBboxPatch((x-w/2, y-h/2), w, h,
boxstyle="round,pad=0.05",
fc=color, ec=color, alpha=alpha, lw=1.5)
ax.add_patch(rect)
ax.text(x, y+0.05, label, ha='center', va='center', fontsize=7.5, fontweight='bold')
if sublabel:
ax.text(x, y-0.25, sublabel, ha='center', va='center', fontsize=6.5, color=C_NEUTRAL)
def arrow(ax, x1, y1, x2, y2):
ax.annotate("", xy=(x2, y2), xytext=(x1, y1),
arrowprops=dict(arrowstyle="->", color=C_NEUTRAL, lw=0.9))
box(ax, 2, 6.2, 3.0, 0.7, "IQ Residual Stream", "‖r(k)‖ from receiver", C_NEUTRAL)
box(ax, 5.5, 6.2, 2.0, 0.7, "Ingress", "f32 / Q16.16", C_DSFB)
arrow(ax, 3.5, 6.2, 4.5, 6.2)
box(ax, 1.5, 4.8, 2.8, 0.9, "GUM Uncertainty\nBudget", "CRLB floor check", "#d4af37")
box(ax, 5.0, 4.8, 2.8, 0.9, "WSS Verification\n+ RAT", "stationarity.rs", "#d4af37")
box(ax, 8.5, 4.8, 2.4, 0.9, "Permutation\nEntropy", "complexity.rs", "#d4af37")
for x2 in [1.5, 5.0, 8.5]: arrow(ax, 5.5, 5.85, x2, 5.25)
box(ax, 5.0, 3.3, 8.5, 1.0, "DSFB Grammar Engine", "envelope · grammar · heuristics · DSA · Lyapunov", C_DSFB)
for x2 in [1.5, 5.0, 8.5]: arrow(ax, x2, 4.35, 5.0, 3.8)
box(ax, 1.5, 1.8, 2.4, 0.9, "Attractor\nD₂ + Koopman", "attractor.rs", C_COMPARATOR)
box(ax, 4.2, 1.8, 2.4, 0.9, "TDA Betti₀\nPersistence", "tda.rs", C_COMPARATOR)
box(ax, 7.0, 1.8, 2.6, 0.9, "Physics Model\nArrhenius/Allan", "physics.rs", C_BOUNDARY)
for x2 in [1.5, 4.2, 7.0]: arrow(ax, 5.0, 2.8, x2, 2.25)
box(ax, 3.0, 0.5, 3.0, 0.7, "Pragmatic Gate", "SOSA backplane eff.", C_ADMISSIBLE)
box(ax, 7.5, 0.5, 2.8, 0.7, "DNA Auth.", "hardware fingerprint", C_ADMISSIBLE)
for x2 in [3.0, 7.5]: arrow(ax, 5.0, 1.35, x2, 0.85)
plt.tight_layout()
save(fig, 38, "information_physics_architecture", out_dir, dpi)
def plot_fig39(out_dir, dpi):
r = rng()
n = 300; tau = 3
t = np.arange(n)
scenarios = [
("WGN only", r.standard_normal(n)*0.1 + 0.3),
("FHSS hopper\n(2 freq)", np.where(t % 20 < 10, 0.15, 0.40) + r.standard_normal(n)*0.01),
("PA thermal drift", 0.05 + t*0.001 + r.standard_normal(n)*0.02),
("CW jammer onset", np.concatenate([r.standard_normal(150)*0.05+0.3,
r.standard_normal(150)*0.02+0.6])),
]
colors = [C_DSFB, C_ADMISSIBLE, C_BOUNDARY, C_VIOLATION]
fig, axes = plt.subplots(2, 4, figsize=(10, 5))
for j, ((lab, x), col) in enumerate(zip(scenarios, colors)):
ax_ts = axes[0, j]
ax_ph = axes[1, j]
ax_ts.plot(t, x, color=col, lw=0.7)
ax_ts.set_title(lab, fontsize=7.5)
ax_ts.set_xlabel("k", fontsize=7)
ax_ts.set_ylabel("‖r‖", fontsize=7)
xd, xd2 = x[:-tau], x[tau:]
ax_ph.scatter(xd, xd2, c=col, s=2, alpha=0.4)
ax_ph.set_xlabel(r"$x(k)$", fontsize=7)
ax_ph.set_ylabel(r"$x(k{-}\tau)$", fontsize=7)
ax_ph.set_aspect("equal", adjustable="datalim")
fig.suptitle("Fig 39. Delay-Coordinate Attractor Reconstruction (τ=3)\n"
"WGN → stochastic ball; FHSS → figure-8; drift → spiral; jammer → step",
fontsize=8.5)
plt.tight_layout()
save(fig, 39, "multimode_attractor_reconstruction", out_dir, dpi)
def plot_fig40(out_dir, dpi):
categories = [
"GUM\nUncertainty", "Physics\nModels", "Topological\nAnalysis",
"Perm. Entropy\n(TDS)", "Attractor\nReconstruction", "Hardware\nDNA",
"Pragmatic\nGating", "Bit-\nExactness", "Fixed-Point\nFPGA", "no_std\nno_alloc"
]
n = len(categories)
dsfb = [9.5, 9.0, 9.0, 9.5, 9.0, 9.5, 9.5, 9.5, 9.5, 10.0]
typical = [3.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0, 5.0, 3.0, 4.0]
mlclass = [4.0, 4.0, 6.0, 5.0, 5.0, 3.0, 2.0, 3.0, 1.0, 0.0]
angles = [i * 2 * np.pi / n for i in range(n)] + [0]
fig, ax = plt.subplots(figsize=(5.5, 5.5), subplot_kw=dict(polar=True))
for vals, col, lab in [(dsfb, C_DSFB, "dsfb-rf Phase 4"), (typical, C_NEUTRAL, "Typical SDR lib"),
(mlclass, C_COMPARATOR, "ML classifier")]:
vals_plot = vals + [vals[0]]
ax.plot(angles, vals_plot, color=col, lw=1.5, label=lab)
ax.fill(angles, vals_plot, color=col, alpha=0.08)
ax.set_xticks(angles[:-1])
ax.set_xticklabels(categories, size=7)
ax.set_ylim(0, 10)
ax.set_yticks([2, 4, 6, 8, 10]); ax.set_yticklabels(["2","4","6","8","10"], size=6)
ax.set_title("Fig 40. DSFB-RF Phase-4 Capability Radar\n"
"(Information-Physics Stack vs Alternatives)", pad=20)
ax.legend(loc="upper right", bbox_to_anchor=(1.35, 1.1), fontsize=7)
plt.tight_layout()
save(fig, 40, "capability_radar_phase4", out_dir, dpi)
FIGURES = {
21: plot_fig21, 22: plot_fig22, 23: plot_fig23, 24: plot_fig24,
25: plot_fig25, 26: plot_fig26, 27: plot_fig27, 28: plot_fig28,
29: plot_fig29, 30: plot_fig30, 31: plot_fig31, 32: plot_fig32,
33: plot_fig33, 34: plot_fig34, 35: plot_fig35, 36: plot_fig36,
37: plot_fig37, 38: plot_fig38, 39: plot_fig39, 40: plot_fig40,
}
def main():
parser = argparse.ArgumentParser(description="Generate DSFB-RF Phase-4 figures 21–40")
parser.add_argument("--fig", nargs="*", type=int,
help="Specific figure numbers (default: all 21-40)")
parser.add_argument("--dpi", type=int, default=150, help="Output DPI (default 150)")
parser.add_argument("--out", type=str, default="figs",
help="Output directory (default: ./figs)")
args = parser.parse_args()
script_dir = Path(__file__).parent
out_dir = script_dir / args.out
requested = set(args.fig) if args.fig else set(FIGURES.keys())
unknown = requested - set(FIGURES.keys())
if unknown:
print(f"Unknown figure numbers: {sorted(unknown)}", file=sys.stderr)
for idx in sorted(requested & set(FIGURES.keys())):
print(f"Generating fig {idx:02d}…")
FIGURES[idx](out_dir, args.dpi)
print(f"\nDone — {len(requested & set(FIGURES.keys()))} figures written to {out_dir}/")
if __name__ == "__main__":
main()