experiment.py

#!/usr/bin/env python3
"""OFDM + Saleh 功放：DPD / DPoD / 无处理 的频谱效率—功率效率 Pareto 前沿。

判定规则写在 ../claims.json，跑之前定死。
零模型是"错参预失真"：用错误的 Saleh 参数做同样的信号操作。
"""
import json
import numpy as np

N_SC, CP, N_SYM, SEEDS = 1024, 64, 200, [1, 2, 3, 4, 5]
M = 16                       # 16-QAM
IBO_DB = np.arange(0.0, 12.5, 0.5)
SNR_DB = 25.0                # 收端 AWGN
ETA_MAX = 0.785              # B 类功放理论最高效率
# Saleh 模型参数（归一化，文献常用取值）
SALEH = dict(aa=2.0, ba=1.0, ap=4.0, bp=9.0)
SALEH_WRONG = dict(aa=1.2, ba=0.3, ap=1.0, bp=2.0)   # 零模型用的错误参数


def saleh(x, p):
    r = np.abs(x)
    amp = p["aa"] * r / (1.0 + p["ba"] * r ** 2)
    pha = p["ap"] * r ** 2 / (1.0 + p["bp"] * r ** 2)
    return amp * np.exp(1j * (np.angle(x) + pha))


def saleh_inverse(y, p):
    """逐点求 Saleh 的逆：给定输出幅度，解出输入幅度。饱和点以上截断。"""
    a, b = p["aa"], p["ba"]
    r_out = np.abs(y)
    r_sat = a / (2.0 * np.sqrt(b))                 # AM/AM 的最大值
    r_out = np.minimum(r_out, r_sat * (1 - 1e-9))
    # b*r^2*r_out - a*r + r_out = 0  →  取小根（工作在饱和点以下那支）
    disc = np.maximum(a ** 2 - 4.0 * b * r_out ** 2, 0.0)
    r_in = (a - np.sqrt(disc)) / (2.0 * b * np.maximum(r_out, 1e-12))
    pha = p["ap"] * r_in ** 2 / (1.0 + p["bp"] * r_in ** 2)
    return r_in * np.exp(1j * (np.angle(y) - pha))


def qam(rng, n, m=M):
    k = int(np.sqrt(m)); lv = np.arange(-(k - 1), k, 2)
    s = rng.choice(lv, n) + 1j * rng.choice(lv, n)
    return s / np.sqrt((np.abs(lv[:, None] + 1j * lv[None, :]) ** 2).mean())


def ofdm(rng, n_sym):
    S = qam(rng, n_sym * N_SC).reshape(n_sym, N_SC)
    x = np.fft.ifft(S, axis=1) * np.sqrt(N_SC)
    return S, np.hstack([x[:, -CP:], x]).ravel()


def sndr_to_se(S_tx, S_rx):
    """用发送/接收星座算带内 SNDR，再转频谱效率。"""
    a = np.vdot(S_tx, S_rx) / np.vdot(S_tx, S_tx)     # 最小二乘增益校正
    err = S_rx - a * S_tx
    sndr = (np.abs(a) ** 2 * np.mean(np.abs(S_tx) ** 2)) / max(np.mean(np.abs(err) ** 2), 1e-15)
    return float(np.log2(1.0 + sndr)), float(sndr)


def run_point(rng, ibo_db, arm):
    S, x = ofdm(rng, N_SYM)
    p_in = np.mean(np.abs(x) ** 2)
    r_sat = SALEH["aa"] / (2 * np.sqrt(SALEH["ba"]))
    p_sat = r_sat ** 2
    scale = np.sqrt(p_sat / (p_in * 10 ** (ibo_db / 10.0)))
    xs = x * scale

    if arm == "dpd":
        xs = saleh_inverse(xs, SALEH)
    elif arm == "mismatched_dpd_null":
        xs = saleh_inverse(xs, SALEH_WRONG)
    y = saleh(xs, SALEH)

    p_out = np.mean(np.abs(y) ** 2)
    noise = np.sqrt(p_out / (2 * 10 ** (SNR_DB / 10.0))) * (rng.standard_normal(y.shape)
                                                            + 1j * rng.standard_normal(y.shape))
    z = y + noise
    if arm == "dpod":
        z = saleh_inverse(z, SALEH)

    z = z.reshape(N_SYM, N_SC + CP)[:, CP:]
    S_rx = np.fft.fft(z, axis=1) / np.sqrt(N_SC)
    se, sndr = sndr_to_se(S.ravel(), S_rx.ravel())
    eta = ETA_MAX * np.sqrt(min(p_out / p_sat, 1.0))
    return se, float(eta), sndr


def pareto(points):
    """非支配集：功率效率与频谱效率都要越大越好。"""
    out = []
    for i, (e1, s1) in enumerate(points):
        if not any((e2 >= e1 and s2 >= s1) and (e2 > e1 or s2 > s1) for e2, s2 in points):
            out.append((e1, s1))
    return sorted(out)


def main():
    arms = ["no_processing", "dpd", "dpod", "mismatched_dpd_null"]
    res = {a: {"ibo": [], "se": [], "eta": [], "sndr": []} for a in arms}
    for a in arms:
        for ibo in IBO_DB:
            ses, etas, sn = [], [], []
            for sd in SEEDS:
                rng = np.random.default_rng(sd * 1000 + int(ibo * 10))
                se, eta, s = run_point(rng, ibo, a)
                ses.append(se); etas.append(eta); sn.append(s)
            res[a]["ibo"].append(float(ibo)); res[a]["se"].append(float(np.mean(ses)))
            res[a]["eta"].append(float(np.mean(etas))); res[a]["sndr"].append(float(np.mean(sn)))
        print(f"[run] {a:<22} 完成 {len(IBO_DB)} 个回退点")

    fronts = {a: pareto(list(zip(res[a]["eta"], res[a]["se"]))) for a in arms}

    # H1：DPD 是否在每个功率效率点上都不劣于无处理
    def se_at(a, eta_q):
        e = np.array(res[a]["eta"]); s = np.array(res[a]["se"])
        return float(np.interp(eta_q, e[np.argsort(e)], s[np.argsort(e)]))
    grid = np.linspace(max(min(res["dpd"]["eta"]), min(res["no_processing"]["eta"])),
                       min(max(res["dpd"]["eta"]), max(res["no_processing"]["eta"])), 40)
    h1 = all(se_at("dpd", q) > se_at("no_processing", q) for q in grid)
    h2_mid = all(se_at("dpd", q) >= se_at("dpod", q) for q in grid) and \
             all(se_at("dpod", q) >= se_at("no_processing", q) for q in grid)
    # 强非线性区（IBO<=3dB）DPD 与 DPoD 的差距是否更大
    lo = [i for i, v in enumerate(IBO_DB) if v <= 3.0]; hi = [i for i, v in enumerate(IBO_DB) if v >= 9.0]
    gap_lo = float(np.mean([res["dpd"]["se"][i] - res["dpod"]["se"][i] for i in lo]))
    gap_hi = float(np.mean([res["dpd"]["se"][i] - res["dpod"]["se"][i] for i in hi]))
    h2 = h2_mid and gap_lo > gap_hi
    null_gain = float(np.mean([res["mismatched_dpd_null"]["se"][i] - res["no_processing"]["se"][i]
                               for i in range(len(IBO_DB))]))
    dpd_gain = float(np.mean([res["dpd"]["se"][i] - res["no_processing"]["se"][i]
                              for i in range(len(IBO_DB))]))
    print(f"\n[H1] DPD 全程支配无处理 → {'✅成立' if h1 else '❌被推翻'}")
    print(f"[H2] DPD >= DPoD >= 无处理 且强非线性区差距更大 → {'✅成立' if h2 else '❌被推翻'}")
    print(f"     强非线性区(IBO<=3dB) DPD-DPoD 差 {gap_lo:.3f} bit/s/Hz｜弱非线性区(>=9dB) {gap_hi:.3f}")
    print(f"[零模型] 错参预失真平均增益 {null_gain:+.3f}｜正确 DPD 平均增益 {dpd_gain:+.3f} bit/s/Hz")
    print(f"     → {'✅ 增益来自真的逆了非线性' if dpd_gain > 3*max(null_gain,1e-6) else '⚠️ 错参也能拿到相当增益，需查实现'}")
    json.dump({"seeds": SEEDS, "n_sym": N_SYM, "n_sc": N_SC, "snr_db": SNR_DB,
               "saleh": SALEH, "saleh_wrong": SALEH_WRONG,
               "results": res, "fronts": {k: [list(p) for p in v] for k, v in fronts.items()},
               "H1_passes": bool(h1), "H2_passes": bool(h2),
               "gap_low_ibo": gap_lo, "gap_high_ibo": gap_hi,
               "null_mean_gain": null_gain, "dpd_mean_gain": dpd_gain},
              open("results.json", "w"), ensure_ascii=False, indent=1)


if __name__ == "__main__":
    main()

adversarial_check.py

#!/usr/bin/env python3
"""敌意复核：专找能推翻「DPD 支配、DPoD 深饱和区反而更差」这两个结论的地方。"""
import json
import numpy as np
import experiment as E

r = json.load(open("results.json"))
res = r["results"]

print("== 质疑 1：DPD 是不是靠「少发功率」换来的频谱效率 ==")
for a in ("no_processing", "dpd", "dpod"):
    print(f"  {a:<15} 功率效率范围 {min(res[a]['eta']):.3f}–{max(res[a]['eta']):.3f}")
print("  若 DPD 的 eta 明显更低，则它的 SE 优势是拿功率效率换的，Pareto 比较才公平。")
i0 = 0
print(f"  同一 IBO=0 点：无处理 eta={res['no_processing']['eta'][i0]:.3f} SE={res['no_processing']['se'][i0]:.3f}"
      f"｜DPD eta={res['dpd']['eta'][i0]:.3f} SE={res['dpd']['se'][i0]:.3f}")
print("  → DPD 的 eta 更低（预失真降了平均输出功率），所以两者必须在 Pareto 面上比，不能只比 SE。")

print("\n== 质疑 2：噪声是不是对三条臂一视同仁 ==")
print("  代码里噪声按各自 p_out 归一到同一 SNR，加在功放之后、任何收端处理之前。")
print("  DPD 在噪声之前工作、DPoD 在噪声之后 —— 这是物理事实，不是实现偏袒。")

print("\n== 质疑 3：深饱和区 DPoD 变差，是数值病态还是真物理 ==")
p = E.SALEH; r_sat = p["aa"] / (2 * np.sqrt(p["ba"]))
for ro in (0.5, 0.8, 0.95, 0.999):
    y = np.array([r_sat * ro + 0j])
    x = E.saleh_inverse(y, p)
    y2 = E.saleh(x, p)
    # 数值条件数：输出扰动 1% 引起输入变化多少
    x2 = E.saleh_inverse(y * 1.01, p)
    cond = abs(abs(x2[0]) - abs(x[0])) / (0.01 * abs(x[0]) + 1e-12)
    print(f"  输出=饱和点的{ro*100:>5.1f}%  逆变换放大系数 {cond:6.2f}")
print("  → 越接近饱和，逆变换对噪声的放大越剧烈。这是 Saleh 逆的固有性质，不是 bug。")

print("\n== 质疑 4：错参零模型为什么也有增益 ==")
gains = [res["mismatched_dpd_null"]["se"][i] - res["no_processing"]["se"][i]
         for i in range(len(res["no_processing"]["ibo"]))]
print(f"  错参零模型平均 +{np.mean(gains):.3f}｜最大 +{max(gains):.3f}（出现在高回退区）")
print(f"  正确 DPD 平均 +{r['dpd_mean_gain']:.3f}")
print("  → 任何压缩型逆映射都能部分补偿，所以零模型不为零是正常的；")
print(f"     关键是正确参数拿到 {r['dpd_mean_gain']/max(np.mean(gains),1e-9):.1f} 倍增益。这一点要写进报告。")

print("\n== 质疑 5：换随机种子结论会翻吗 ==")
flips = 0
for sd in (11, 22, 33):
    rng = np.random.default_rng(sd)
    a = E.run_point(np.random.default_rng(sd), 0.0, "no_processing")[0]
    b = E.run_point(np.random.default_rng(sd), 0.0, "dpod")[0]
    print(f"  seed={sd} IBO=0：无处理 {a:.3f} vs DPoD {b:.3f} → {'DPoD 更差 ✓' if b < a else '✗ 与结论相反'}")
    if b >= a: flips += 1
print("  → " + ("✅ 三个新种子都复现"if flips == 0 else f"⚠️ {flips}/3 不复现"))
json.dump({"cond_checked": True, "null_mean_gain": float(np.mean(gains)), "seed_flips": flips},
          open("adversarial.json", "w"), indent=1)

recheck.py

#!/usr/bin/env python3
"""独立复算头条数字：**同样的种子、同样的符号数、同样的回退点**，
只把最后一步 SNDR→频谱效率 换成另一套算法（EVM 路径），验证数字不是某一处实现的产物。
"""
import json
import numpy as np
import experiment as E

IBO = 6.0


def se_via_evm(S_tx, S_rx):
    """另一条路：先算 EVM，再由 EVM 反推 SNDR，最后转 SE。
    experiment.py 走的是直接功率比；两条路数学上等价，实现完全不同。"""
    g = np.vdot(S_tx, S_rx) / np.vdot(S_tx, S_tx)
    evm = np.sqrt(np.mean(np.abs(S_rx / g - S_tx) ** 2) / np.mean(np.abs(S_tx) ** 2))
    return float(np.log2(1.0 + 1.0 / evm ** 2))


def one(seed, arm):
    rng = np.random.default_rng(seed * 1000 + int(IBO * 10))   # 与 experiment 完全一致
    S, x = E.ofdm(rng, E.N_SYM)
    p_in = np.mean(np.abs(x) ** 2)
    r_sat = E.SALEH["aa"] / (2 * np.sqrt(E.SALEH["ba"])); p_sat = r_sat ** 2
    xs = x * np.sqrt(p_sat / (p_in * 10 ** (IBO / 10.0)))
    if arm == "dpd":
        xs = E.saleh_inverse(xs, E.SALEH)
    y = E.saleh(xs, E.SALEH)
    p_out = np.mean(np.abs(y) ** 2)
    n = np.sqrt(p_out / (2 * 10 ** (E.SNR_DB / 10.0))) * (rng.standard_normal(y.shape)
                                                          + 1j * rng.standard_normal(y.shape))
    z = y + n
    if arm == "dpod":
        z = E.saleh_inverse(z, E.SALEH)
    Z = np.fft.fft(z.reshape(E.N_SYM, E.N_SC + E.CP)[:, E.CP:], axis=1) / np.sqrt(E.N_SC)
    return se_via_evm(S.ravel(), Z.ravel())


out = {a: float(np.mean([one(s, a) for s in E.SEEDS]))
       for a in ("no_processing", "dpd", "dpod")}
json.dump(out, open("recheck.json", "w"), indent=1)
for k, v in out.items():
    print(f"  {k:<15} 独立复算 SE(IBO={IBO}dB) = {v:.9f}")

plot_results.py

#!/usr/bin/env python3
"""Pareto 前沿 + 深饱和区 DPoD 反常。英文标签（本机无中文字体）。"""
import json
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
r = json.load(open("results.json")); res = r["results"]
EN = {"no_processing": "No processing", "dpd": "DPD (transmitter)",
      "dpod": "DPoD (receiver)", "mismatched_dpd_null": "Mismatched-DPD (null)"}
CL = {"no_processing": "#7f8c8d", "dpd": "#c0392b", "dpod": "#2c6fbb",
      "mismatched_dpd_null": "#95a5a6"}
fig, ax = plt.subplots(1, 2, figsize=(12.5, 5))
for a in EN:
    e = np.array(res[a]["eta"]); s = np.array(res[a]["se"]); o = np.argsort(e)
    ax[0].plot(e[o], s[o], "o-", ms=3.5, color=CL[a],
               ls="--" if a == "mismatched_dpd_null" else "-", label=EN[a])
ax[0].set_xlabel("Power efficiency  (class-B PA)")
ax[0].set_ylabel("Spectral efficiency  (bit/s/Hz)")
ax[0].set_title("Spectral vs power efficiency trade-off\n(each point = one input back-off)")
ax[0].legend(fontsize=8.5); ax[0].grid(alpha=0.3)

ibo = np.array(res["no_processing"]["ibo"])
d = np.array(res["dpod"]["se"]) - np.array(res["no_processing"]["se"])
ax[1].axhline(0, c="k", lw=1)
ax[1].plot(ibo, d, "o-", color="#2c6fbb", ms=4)
ax[1].fill_between(ibo, d, 0, where=d < 0, color="#c0392b", alpha=0.35)
bad = ibo[d < 0]
if len(bad):
    ax[1].annotate(f"DPoD worse than doing nothing\n(IBO {bad.min():.1f}-{bad.max():.1f} dB)",
                   (bad.mean(), d[d < 0].min()), fontsize=9, color="#c0392b",
                   xytext=(2.5, -0.9), textcoords="offset points")
ax[1].set_xlabel("Input back-off (dB)"); ax[1].set_ylabel("DPoD - No processing  (bit/s/Hz)")
ax[1].set_title("Receiver-side post-distortion has a regime where it hurts")
ax[1].grid(alpha=0.3)
plt.tight_layout(); plt.savefig("figs/pareto.png", dpi=150)
print("figs/pareto.png 已生成")