#!/usr/bin/env python3 # -*- coding: utf-8 -*- """N 点最小二乘标定求解器:**全仓唯一实现**(v2.2 起固件侧拟合已删除)。 为什么是唯一实现:标定是台架一次性离线动作——注入已知量 → 读屏上 `meas:` → 在 PC 上解出 r/b → 经 104 写入设备。设备本体只需要接受校正后的参数,不需要 自己拟合曲线。固件里曾有过一份 C 实现,生产代码零调用(被 linker 移除), 留着只会变成第二个真相源,v2.2 已整段删除。所以**本脚本的整数算术就是口径**, `--selftest` 是它的回归锚点(原 C 侧金标向量已搬进来)。 把台架注入的一串 (屏上读数, 注入已知量) 点拟合成 y = a*x + c(普通最小二乘), 再折回标定层能落盘的 (r_q16, b): r_new_q16 = round(a * r_old_q16) # 一次项 = 端子↔引脚变比 c_pin_mV = c * 65536 / r_new_q16 # 截距从二次值域折回引脚域 b_new = b_old + round(c_pin_mV * 65536 * 1000 / k_q16) 推导(disp 单位 mV):折算模型 disp = (code + b) × (k_q16/65536) × 1e-3 × (r_q16/65536),把 disp_new = a×disp_old + c 的 code 一次项与常数项分别对齐, 一次项给出 r_new = a·r_old,常数项给出上式 b_new(b_old = 写之前设备里的零点, 必须作为入参传进来,因为读数是在那个零点下读的)。 为什么用两点之外的算法:两点没有冗余,单点读数一抖解出的 r 就整条平移;N 点最小二乘 把残差摊到所有点上,并且能顺手判「非线性/削顶/接线反」——这正是两点的盲区。 全程 int64 定点(不用 float 做累加):结果确定、可跨平台复现,同一组输入永远 解出同一个 r_q16。整数路径的圆整次序是本脚本语义的一部分,改一处就要同步 --selftest 的期望值。 用法(仓库根目录,解释器用 .venv): 自检(回归锚点,改动算法后必跑): .venv/Scripts/python.exe test/host/cal_fit.py --selftest 单路(x = `meas:` 行原样读数 mV/mA;y = 继保仪注入值 V/A,内部 ×1000): .venv/Scripts/python.exe test/host/cal_fit.py --sem ua \\ --x 5021.5,10041.77,25000,50024.09,99933.02 --y 5,10,25,50,100 批量(读执行卡 之间的长表,每行一个点): .venv/Scripts/python.exe test/host/cal_fit.py --from-doc \\ docs/测试说明/专项验证/ADC多点标定执行流程.md 公共参数: --r-old 本次读数采集时设备里的变比(浮点):默认 ua/un=1000.0、ia=29.4118; 也可 ua=1000,un=1000,ia=29.4118 形式逐路给 --b-old 本次读数采集时设备里的零点(整数码值,默认 0):裸值或 ua=0,ia=-2 形式 ⚠ **同一路的 x 必须在同一个 r 下采集**:拿旧 r 时代的读数和新 r 时代的读数混着 做直线拟合,斜率会被系统性掰弯(表现是解出的 r 差几十个百分点,而合法带仍然放行)。 |r_new/r_old - 1| > 10% 时打 WARN 就是给这件事提个醒。 判据(任一不满足即 REJECT,打印原因,不算可写): n 不在 [3, 8] 冗余不足或超上限 den <= 0 示值全同 / 点位没拉开 x 与 y 相邻增量方向相反 接线反、通道串、削顶的典型征兆 a <= 0 零斜率/负斜率 r_new_q16 不在 [1.0, 32767.99] 与 run_par_cal_r_set() 的粗守同带 |b_new| > 8192 零点越带(先查前端,不用系数凑) 最大残差 > 满量程跨度 2% 非线性/削顶 输出纯 ASCII(中文只在本注释里);exit 0 = 全部解出且合法,2 = 有拒绝/未填行; --selftest 的 exit 0 = 全部断言通过,1 = 有失败。 创建日期:2026-11-13(两点 → N 点最小二乘;同日 v2.2 固件侧拟合删除,本脚本转为唯一实现) """ import argparse import re import sys # ---- 判据参数(本脚本是唯一实现,这些值即为口径;改动会使 GOLDEN 期望值失效)---- # ⚠ 前 5 个是**独立于固件**的算法参数(原先与已删除的 ADC_CAL_FIT_* 同值, # 现在以本文件为准);后 4 个是设备侧合法带,必须与 App/algo/adc_cal.h / # run_par.c 保持一致——那两处仍是设备侧安全闸(写参数时的拦截),不随拟合删除。 K_Q16 = 10000000 # 该路前端斜率(UA/UN/IA 三路同源,同值) B_MAX = 8192 # 设备侧 ADC_CAL_B_MAX(零点幅值上限,码值) R_MIN = 65536 # 设备侧 ADC_CAL_R_MIN(比值 1.0) R_MAX = 2147483647 # 设备侧 ADC_CAL_R_MAX(int32 上限) MILLI_MAX = 200000 # 单点值域上界(mV / mA),±200 V / ±200 A C_MAX = 100000 # 截距 c 上界(mV),±100 V N_MIN = 3 # 点数下限(少于 3 点无冗余) N_MAX = 8 # 点数上限 RES_PCT = 2 # 最大残差 ÷ 满量程跨度上限(百分比) # 比值上界的两位小数**向下取整**得 32767.99,不能用 R_MAX/65536.0 直接 %.2f: # 2147483647/65536 = 32767.99998474,四舍五入会印成 32768.00,与 adc_cal.h / # run_par.c 的"1.0 ~ 32767.99"口径差一个数(提示文案里的带边界必须和真带一致)。 R_RATIO_MAX = (R_MAX * 100 // 65536) / 100.0 # 32767.99(仅用于提示文案) # 语义槽位(App/algo/adc_input.h 枚举序:UA0/UB1/UC2/UN3/IA4/IB5/IC6/IN7) SEM = { "ua": dict(slot=0, r_old=1000.0), # 电压族 TPAFE "un": dict(slot=3, r_old=1000.0), # 电压族 TPAFE "ia": dict(slot=4, r_old=29.4118), # 电流族 TPAFE } # 拒绝原因码(负数):本脚本的对外契约,--selftest 与 REASON 文案都按这张表 ERR_ARG = -1 ERR_POINTS = -2 ERR_DEN = -3 ERR_MONO = -4 ERR_SLOPE = -5 ERR_R_BAND = -6 ERR_B_BAND = -7 ERR_RESID = -8 REASON = { ERR_ARG: "bad argument (|x|,|y| < " + str(MILLI_MAX) + ", k_q16 in (0,1e8], r_old in band)", ERR_POINTS: "n outside [%d,%d]: got %%d point(s)" % (N_MIN, N_MAX), ERR_DEN: "x span collapsed (den<=0): readings not spread", ERR_MONO: "non-monotonic: adjacent x/y increments oppose", ERR_SLOPE: "slope <= 0: readings flat or reversed vs injected value", ERR_R_BAND: "r_new out of band [1.0, %.2f]" % R_RATIO_MAX, ERR_B_BAND: "|b_new| beyond %d codes" % B_MAX, ERR_RESID: "residual > %d%% of y span (non-linearity / clipping)" % RES_PCT, } # ------------------------------------------------------------ 整数算术(唯一实现) # v2.2 起这里不再是"镜像":App/algo/adc_cal.c 的拟合算术已整段删除,本文件就是口径。 def div_trunc(num, den): """整数除法:向零截断(Python // 向负无穷取整,不能直接用)。""" q = abs(num) // abs(den) return q if (num < 0) == (den < 0) else -q def div_round(num, den): """四舍五入(half away from zero);den 恒 > 0。 这是本仓库拟合算术**唯一的圆整口径**,不要换成 `(num + den/2)//den` 那种 写法:Python 的 // 对负数是向负无穷取整,负数会在 -0.5 处折向 0, 与期望的 half-away-from-zero 不同。先判符号再取绝对值,不依赖语言行为。 """ if num >= 0: return (num + den // 2) // den return -((-num + den // 2) // den) def round_int(v): """浮点 → 整数:二分之一远离零圆整(入参量化口径,与 to_milli 配对)。""" return int(v + 0.5) if v >= 0 else -int(-v + 0.5) def to_milli(v): """读数 → 整数毫单位(mV / mA 域,即 fit() 的直接入参)。""" return round_int(v) def fit(xs, ys, k_q16=K_Q16, r_old_q16=None, b_old=0): """N 点最小二乘标定求解(唯一实现)。返回 (rc, info);rc == 0 时 info 是解。 xs / ys 已是整数 mV / mA;r_old_q16 = 本轮读数采集时的变比 Q16 原值。 整数路径决定 rc 与 r_q16/b;info 里的 a/res%/R2 是报告用浮点量,不参与判定。 """ n = len(xs) if r_old_q16 is None: return ERR_ARG, None if (k_q16 <= 0) or (k_q16 > 100000000): return ERR_ARG, None if (r_old_q16 < R_MIN) or (r_old_q16 > R_MAX): return ERR_ARG, None if (n < N_MIN) or (n > N_MAX): return ERR_POINTS, None for v in list(xs) + list(ys): if abs(v) >= MILLI_MAX: return ERR_ARG, None # 一趟:累加 + 单调性(相邻点方向必须一致) sx = sy = sxy = sxx = 0 for i in range(n): x = xs[i] y = ys[i] sx += x sy += y sxy += x * y sxx += x * x if i > 0: dx = x - xs[i - 1] dy = y - ys[i - 1] if (dx > 0 and dy < 0) or (dx < 0 and dy > 0): return ERR_MONO, None num = n * sxy - sx * sy den = n * sxx - sx * sx if den <= 0: return ERR_DEN, None a_q16 = div_round(num * 65536, den) if a_q16 <= 0: return ERR_SLOPE, None # 先按除法把 a_q16 夹到"r_new 不会越上界"的范围,再做 a_q16 × r_old—— # 顺序反了会在极端点位下把 int64 顶穿(a_q16 可到 1e17 量级) if a_q16 > (R_MAX * 65536) // r_old_q16: return ERR_R_BAND, None r_new_q16 = (a_q16 * r_old_q16 + 32768) >> 16 if (r_new_q16 < R_MIN) or (r_new_q16 > R_MAX): return ERR_R_BAND, None # 截距 c:y = a·x + c ⇒ c = (Σy − a·Σx) / n。两遍式(用已算好的 a_q16)。 # 若改成从 num/den 代数展开的 (Σy·den − num·Σx)/(n·den),分子峰值到 1.2e19, # 超出 int64 上界——Python 的大整数不会溢出,但这条路径要能直接照抄成 # 定点实现,所以这里**刻意保持 int64 可容纳的两遍式**,别改成展开式 c_milli = div_trunc(sy * 65536 - a_q16 * sx, n * 65536) if c_milli > C_MAX or c_milli < -C_MAX: return ERR_B_BAND, None # 截距域离谱 = b 必然越带(且护住下一步的 ×2^32) # 截距折回引脚域,保留 Q16 亚毫伏精度(截到整 mV 会带来约 6.5 码 / mV 的 b # 量化,把真正的细分零点修正抹平);两步走,别一次乘完(会顶穿 int64) c_pin_q16 = div_round(c_milli * 65536 * 65536, r_new_q16) db = div_round(c_pin_q16 * 1000, k_q16) b_new = b_old + db if (b_new > B_MAX) or (b_new < -B_MAX): return ERR_B_BAND, None # 残差:|y_i − (a·x_i + c)| 与满量程跨度比 y_min = min(ys) y_max = max(ys) y_span = y_max - y_min res_max = 0 for i in range(n): pred = div_round(a_q16 * xs[i], 65536) + c_milli res = ys[i] - pred if res < 0: res = -res if res > res_max: res_max = res if (res_max * 100) > (y_span * RES_PCT): return ERR_RESID, None # 报告用浮点量(不参与判定) a_f = a_q16 / 65536.0 y_mean = float(sy) / n ss_res = 0.0 ss_tot = 0.0 for i in range(n): d = ys[i] - (a_f * xs[i] + c_milli) t = ys[i] - y_mean ss_res += d * d ss_tot += t * t r2 = (1.0 - ss_res / ss_tot) if ss_tot > 0.0 else 0.0 return 0, dict(n=n, a_q16=a_q16, a=a_f, c_milli=c_milli, r_new_q16=r_new_q16, r_new=r_new_q16 / 65536.0, c_pin_mV=c_pin_q16 / 65536.0, db=db, b=b_new, res_max=res_max, res_pct=(100.0 * res_max / y_span) if y_span else 0.0, r2=r2, y_span=y_span, dev_pct=100.0 * (r_new_q16 / float(r_old_q16) - 1.0)) # ---------------------------------------------------------------- 输入解析 def parse_floats(text): return [float(s) for s in text.split(",") if s.strip() != ""] def parse_map(text, defaults, conv): """--r-old / --b-old:裸值作用于所有路,或 "ua=…,ia=…" 逐路给。""" if text is None: return dict(defaults) if "=" not in text: v = conv(text) return {k: v for k in defaults} out = dict(defaults) for part in text.split(","): part = part.strip() if not part: continue k, v = part.split("=", 1) k = k.strip().lower() if k not in defaults: raise ValueError("unknown sem %r (use ua/un/ia)" % k) out[k] = conv(v) return out def parse_cal_table(path): """提取 ~ 之间的长表。 列序:路 | 序号 | 屏上读数(mV/mA) | 注入已知量(V/A),每行一个点。 返回 {sem: [(idx, x_s, y_s), ...]};标记不存在返回 None。 两格全空 = 该点本轮没测(调用方跳过);只填一格 = 半填错误(报 INCOMPLETE)。 """ with open(path, encoding="utf-8") as f: text = f.read() m = re.search(r"(.*?)", text, re.S) if m is None: return None groups = {} for line in m.group(1).splitlines(): line = line.strip() if not line.startswith("|"): continue cells = [c.strip() for c in line.strip("|").split("|")] sem = cells[0].lower() if sem not in SEM: continue # 表头行 / 分隔行自然跳过 if len(cells) < 4: raise ValueError("row for %r has %d columns; need 4 " "[sem|idx|reading|injected]" % (sem, len(cells))) if len(cells) > 4: # 典型来源:没清空的**旧两点表**(6 列)。静默按前四列解会解出垃圾, # 所以这里明说——比让用户去猜"怎么 REJECT 了"省一轮台架 raise ValueError("row for %r has %d columns (>4): that is the old " "two-point table; the N-point table needs exactly " "4 columns [sem|idx|reading|injected]" % (sem, len(cells))) groups.setdefault(sem, []).append((cells[1], cells[2], cells[3])) return groups def parse_cell(s): if s in ("", "-", "—"): return None return float(s) # ---------------------------------------------------------------- 求解与输出 def solve_and_print(sem, xs, ys, r_old_q16, b_old): """解一路并打印判读;成功返回 (sem, slot, r_new, b),拒绝返回 None。""" ch = SEM[sem] rc, info = fit(xs, ys, K_Q16, r_old_q16, b_old) if rc != 0: reason = REASON[rc] if rc == ERR_POINTS: reason = reason % len(xs) print("sem=%s slot=%d n=%d REJECT %s" % (sem, ch["slot"], len(xs), reason)) return None print("sem=%s slot=%d n=%d OK r=%.6f r_q16=%d b=%d a=%.9f res%%=%.3f " "R2=%.8f" % (sem, ch["slot"], info["n"], info["r_new"], info["r_new_q16"], info["b"], info["a"], info["res_pct"], info["r2"])) if abs(info["dev_pct"]) > 10.0: print("sem=%s WARN r_new dev %+.2f%% vs r_old=%.6f - all points must be " "read at the same r, and units are mV vs V" % (sem, info["dev_pct"], r_old_q16 / 65536.0)) return (sem, ch["slot"], info["r_new"], info["b"]) def write_all(done): """末尾合并写入命令:先全部 r,再全部 b(cal_write.py 的两段顺序契约)。""" parts = ["--r 0x%04X=%.6f" % (0x8050 + s, r) for _, s, r, _ in done] parts += ["--b 0x%04X=%d" % (0x8048 + s, b) for _, s, _, b in done] print("WRITE-ALL: test/host/cal_write.py " + " ".join(parts)) def run_single(a): missing = [n for n, v in (("--sem", a.sem), ("--x", a.x), ("--y", a.y)) if v is None] if missing: print("ERROR: single-route mode needs %s (or use --from-doc)" % ",".join(missing)) return 2 try: xv = parse_floats(a.x) yv = parse_floats(a.y) r_map = parse_map(a.r_old, {s: SEM[s]["r_old"] for s in SEM}, float) b_map = parse_map(a.b_old, {s: 0 for s in SEM}, int) except ValueError as exc: print("ERROR: %s" % exc) return 2 if len(xv) != len(yv): print("ERROR: --x has %d points, --y has %d" % (len(xv), len(yv))) return 2 xs = [to_milli(v) for v in xv] # 屏上读数本来就是 mV/mA ys = [to_milli(v * 1000.0) for v in yv] # 注入已知量 V/A → mV/mA r_old_q16 = round_int(r_map[a.sem] * 65536.0) res = solve_and_print(a.sem, xs, ys, r_old_q16, b_map[a.sem]) if res is None: return 2 write_all([res]) return 0 def run_from_doc(a): try: groups = parse_cal_table(a.from_doc) except OSError as exc: print("BAD-DOC: %s" % exc) return 2 except ValueError as exc: print("BAD-TABLE: %s" % exc) return 2 if groups is None: print("NO-MARKER: no table in %s" % a.from_doc) return 2 if not groups: print("EMPTY: CALDATA table has no ua/un/ia rows") return 2 try: r_map = parse_map(a.r_old, {s: SEM[s]["r_old"] for s in SEM}, float) b_map = parse_map(a.b_old, {s: 0 for s in SEM}, int) except ValueError as exc: print("ERROR: %s" % exc) return 2 done = [] failed = 0 for sem in sorted(groups): xs, ys = [], [] for idx, x_s, y_s in groups[sem]: xv = parse_cell(x_s) yv = parse_cell(y_s) if xv is None and yv is None: continue # 整行空 = 本轮没测这个点 if xv is None or yv is None: print("sem=%s idx=%s INCOMPLETE: need both reading and injected " "value (or leave both empty); skipped" % (sem, idx)) failed += 1 continue xs.append(to_milli(xv)) ys.append(to_milli(yv * 1000.0)) if not xs: print("sem=%s SKIP: no filled points" % sem) continue r_old_q16 = round_int(r_map[sem] * 65536.0) res = solve_and_print(sem, xs, ys, r_old_q16, b_map[sem]) if res is None: failed += 1 else: done.append(res) if done: write_all(done) return 0 if (not failed and done) else 2 # ---------------------------------------------------------------- selftest 回归 # 金标向量:原固件单测(v2.1 的 tc_00b_fit)里的那一组,v2.2 固件侧拟合删除后 # 搬到这里。输入与期望值**全部写死**,不在这里现算——现算等于自证,锚点就没意义了。 # # 数据怎么来的(推导留档,便于将来重构时重新生成): # code = 500/700/900/1100/1300,真值 r_true=600.0、b_true=+12, # k_q16=1e7(TPAFE 152.587890625 µV/码),r_old=1000.0,b_old=0 # 引脚电压 = code × 0.152587890625 mV → x(旧 r 下的屏上读数)= 该值 × 1000/1 # y(注入已知量)= (code + 12) × 0.152587890625 × 600.0 GOLDEN = { "x": [76294, 106812, 137329, 167847, 198364], "y": [46875, 65186, 83496, 101807, 120117], "k_q16": K_Q16, "r_old_q16": 1000 * 65536, "b_old": 0, "expect_rc": 0, "expect_r_q16": 39322000, "expect_b": 12, } # 拒绝路径断言:(名称, xs, ys, r_old_q16, b_old, 期望 rc, 期望原因前缀) NG_CASES = [ ("n<3", [1000, 2000], [1000, 2000], 1000 * 65536, 0, ERR_POINTS, "[3,8]"), ("x all same", [5000, 5000, 5000], [1000, 2000, 3000], 1000 * 65536, 0, ERR_DEN, "x span collapsed"), ("mono reversed", [1000, 2000, 3000, 4000, 5000], [1000, 2000, 3000, 2000, 1000], 1000 * 65536, 0, ERR_MONO, "non-monotonic"), ("r out of band", [1000, 101000, 199999], [1000, 1001, 1002], 1000 * 65536, 0, ERR_R_BAND, "out of band"), ] def run_selftest(): """跑金标向量 + 拒绝路径断言。exit 0 = 全过,1 = 有失败。""" fails = 0 total = 0 def chk(cond, label): nonlocal fails, total total += 1 if cond: print("[PASS] %s" % label) else: print("[FAIL] %s" % label) fails += 1 g = GOLDEN rc, info = fit(g["x"], g["y"], g["k_q16"], g["r_old_q16"], g["b_old"]) chk(rc == g["expect_rc"], "golden: rc == %d (got %d)" % (g["expect_rc"], rc)) if info is None: chk(False, "golden: solver returned no solution") else: chk(info["r_new_q16"] == g["expect_r_q16"], "golden: r_q16 == %d (got %d)" % (g["expect_r_q16"], info["r_new_q16"])) chk(info["b"] == g["expect_b"], "golden: b == %d (got %d)" % (g["expect_b"], info["b"])) for name, xs, ys, ro, bo, want_rc, want_txt in NG_CASES: rc, _ = fit(xs, ys, K_Q16, ro, bo) chk(rc == want_rc, "reject %-14s rc == %d (got %d)" % (name, want_rc, rc)) chk(want_txt in REASON.get(want_rc, ""), "reject %-14s message mentions %r" % (name, want_txt)) print("=== cal_fit selftest: %d/%d PASSED (failed=%d) ===" % (total - fails, total, fails)) return 0 if fails == 0 else 1 def main(): ap = argparse.ArgumentParser( description="N-point least-squares calibration solver " "(sole implementation; firmware no longer fits)") ap.add_argument("--selftest", action="store_true", help="run the golden-vector regression (exit 0 = all passed)") ap.add_argument("--sem", choices=sorted(SEM), default=None) ap.add_argument("--x", default=None, help="comma-separated on-screen readings (mV or mA)") ap.add_argument("--y", default=None, help="comma-separated injected values (V or A; x1000 inside)") ap.add_argument("--r-old", default=None, help="r in the device while readings were taken " "(default ua/un=1000.0, ia=29.4118; or ua=..,ia=..)") ap.add_argument("--b-old", default=None, help="b in the device while readings were taken " "(codes, default 0; or ua=..,ia=..)") ap.add_argument("--from-doc", default=None, help="read the CALDATA-marked long table from a doc") a = ap.parse_args() if a.selftest: return run_selftest() if a.from_doc is not None: return run_from_doc(a) return run_single(a) if __name__ == "__main__": sys.exit(main())