#!/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}")