实现编码与突发错误纠错)
简介RS码Reed-Solomon码作为应用广泛的信道编码方案在通信与数据存储领域承担着纠错重任。资源内提供了一套完整的RS码MATLAB实现面向通信工程、电子信息类学生及科研人员可辅助理解伽罗华域运算、RS编码、译码及AWGN信道下的仿真流程。压缩包共18个文件以15个.m函数脚本为主另有3个.asv自动备份文件整体仅7KB代码轻量清晰。文件覆盖GF(2^m)域生成、多项式乘加、BPSK调制映射、AWGN信道模拟、RS编码、Chien搜索与Forney译码等核心模块便于逐段研读与复用。目前已有275人学习下载适合需要结合理论进行仿真实验或完成课程设计的读者。通过运行demo可直观观察RS码在加性高斯白噪声环境下的纠错表现为后续深入通信系统研究打下基础。1. 突发误码与爆浆错误RS 码仍是信道编码的主心骨在真实信道里错误很少是均匀撒开的一两颗“椒盐”更多是连续一段符号整体被打翻。光盘上的划痕、闪存的坏块、无线信号被汽车点火干扰后的持续衰落这类突发错误让只能纠单个比特的码型立刻失去作用。Reed-Solomon 码RS 码恰恰是站在“符号”而不是“比特”的粒度上做编码译码的一个 8 比特符号整体被改掉在 RS 码眼里只是一个错误符号这与信道编码中常见的单比特纠错码是完全不同的纠错模型。也正因为如此RS 码从光盘、QR 码一路活到了 5G 的 eCPRI 接口和深空通信 CCSDS 标准里始终是信道编码里兜底的那一层。这篇文章不讨论某份现成源码而是从 GF(2^m) 的运算法则出发完整写好一个 RS(255,223) 的编码器与译码器再把它们放进突发错误模型里看参数边界在哪里。2. RS 码的数学底盘GF(2^m) 运算与 RS(n,k) 参数解读2.1 伽罗瓦域四个运算加法用 XOR乘法查表RS 码的“符号”不是普通整数而是有限域 GF(2^m) 里的元素。最常见的 m8 意味着一个符号占 8 bitGF(2^8) 里一共有 256 个合法元素每个元素都对应一个 8 bit 多项式。域上的加法就是按位异或这一点比整数加法简单得多乘法则要靠一个本原多项式做模约减。我一般会提前把“指数-元素”和“元素-指数”两张表算好之后所有乘法都变成查表跑 RS 译码时才不会在每次乘除上浪费周期。# 构建 GF(2^8) 的两张表exp 输出 α^ilog 给出元素的对数 def gf_init(prim0x11d, field256): exp [0] * (2 * field) log [0] * field x 1 for i in range(field - 1): exp[i] x log[x] i x 1 # 左移一位相当于乘以生成元 α if x field: # 溢出到第 9 位按本原多项式约减 x ^ prim # 把表扩到两倍后续 exp[log[a]log[b]] 就不必每次都取模 for i in range(field - 1, 2 * field - 1): exp[i] exp[i - (field - 1)] return exp, log def gf_mul(x, y, exp, log): if x 0 or y 0: return 0 return exp[log[x] log[y]] def gf_inv(x, exp, log): if x 0: raise ValueError(0 没有逆元) return exp[255 - log[x]]这里0x11d是 RS 常用的本原多项式 x^8 x^4 x^3 x^2 1。代码里x field判断左移后是否出现了 x^8 这一项有就异或本原多项式相当于做一次模 2 多项式除法。表建好之后gf_mul的代价只有一次数组寻址gf_inv用255 - log[x]是因为非零元素构成 255 阶循环群求逆即找指数互为相反数的元素。这套查表结构是后面生成多项式、校验子计算、Chien 搜索能跑得动的前提。2.2 最小距离与生成多项式的关系RS(n,k) 码的码字被设计成生成多项式 g(x) 的倍式。给定 2t 个校验符号生成多项式写作g(x) (x α^0)(x α^1)...(x α^(2t-1))只要发送端给出的码字是 g(x) 的倍式接收端用同样一组根去验证就能知道哪里不对。这也直接决定了码的最小距离 d_min n - k 1即码字之间最少可以相差多少个符号。RS 是一个最大距离可分码MDS 码这意味着同样冗余下它能达到的 d_min 已经顶到了理论天花板。信道编码选 RS 时看重的正是“已知 t 个错误符号必定能纠”的这个硬保证而不是那种靠概率才纠得动的软保证。2.3 RS(n,k) 参数表t 与码率的折中参数含义RS(255,223) 的取值m每个符号的比特数8GF(2^8)n码字包含的符号数255 2^8 - 1k一个码字里的信息符号数223n-k校验符号数32t可纠错误符号上限(n-k)/2 16d_min码最小距离n-k1 33码率 R有效信息占比k/n 0.8745这里最容易看走眼的是 t (n-k)/2RS 是按符号纠错的所以“t16”表示一帧里最多能容忍 16 个错误符号而不是 16 个错误比特。一个符号在信道里可能有好几个比特同时翻转但只要这 8 个比特同属一个错误符号计数仍然只算 1。正是这个性质让 RS 特别适合突发错误信道——突发再长只要它落在有限个符号内对 RS 译码器来说就跟几个孤立的错误符号没有区别。3. 用 Python 跑通 RS 编码从构造生成多项式到校验位3.1 一次性构造生成多项式 g(x)编码前的准备只有一件事按校验符号个数 nsym 生成 g(x)。因为 g(x) 的根必须包含 α^0 到 α^(2t-1)多项式构造本质上是把这 2t 个一次因式乘起来。GF(2^8) 里的加法和减法等价所以因式直接写成(x α^i)。def gf_poly_mul(p, q, exp, log): r [0] * (len(p) len(q) - 1) for i, a in enumerate(p): if a 0: continue for j, b in enumerate(q): if b 0: continue r[i j] ^ gf_mul(a, b, exp, log) return r def rs_generator_poly(nsym, exp, log): g [1] for i in range(nsym): # 每次乘一个 (x α^i)最后得到 g(x)∏(xα^i) g gf_poly_mul(g, [1, exp[i]], exp, log) return g这个实现里我用“数组下标代表 x 的幂次”的升序约定也就是poly[0]是 x^0 的系数poly[1]是 x^1 的系数。实际通信帧一般按“先发高次符号”处理所以后面对接协议时记得把符号顺序颠倒一下否则发出去的码字长度和位置都对不上。nsym传 32得到的g就是 33 个系数的数组最高次系数始终是 1。3.2 用循环码除法右移得到校验位让码字成为 g(x) 的倍式最经典的做法是系统码把消息多项式先乘以 x^nsym再用 g(x) 做多项式除法余数就是校验符号直接拼在信息符号后面。这套“除以生成多项式、余作校验”的思路是循环码的标准动作。def gf_poly_div(dividend, divisor, exp, log): out list(dividend) for i in range(len(dividend) - len(divisor) 1): coef out[i] if coef ! 0: for j in range(1, len(divisor)): if divisor[j] ! 0: # GF(2^m) 中减法和加法同为异或 out[i j] ^ gf_mul(divisor[j], coef, exp, log) rem_len len(divisor) - 1 return out[:-rem_len], out[-rem_len:] def rs_encode_msg(msg, nsym, exp, log): gen rs_generator_poly(nsym, exp, log) _, rem gf_poly_div(msg [0] * nsym, gen, exp, log) return msg rem注意gf_poly_div返回两个部分我只取余数。msg [0] * nsym相当于把消息乘上 x^nsym因为低次端补零后原有各项的幂次都没变而校验位最终落在低次端。编码完成后一个码字长度为 k nsym前 k 个符号是原始消息后 nsym 个是校验这正是系统码最友好的输出形式。fmt参数不漏写exp和log是第 2 章里gf_init返回的两张表任何编译出错十有八九是漏传了它们。3.3 自检编码器余数为零才算对编码器写完后不要急着调译码先验证一个基本事实任意编码结果除以 g(x)余数必须是全零。再顺手翻转 t 个以内的错误符号用第 4 章的译码器跑回来比对。exp, log gf_init() nsym 32 msg [0x0a] * 223 # 223 个信息符号全部填 0x0a codeword rs_encode_msg(msg, nsym, exp, log) # 自检一能被生成多项式整除 gen rs_generator_poly(nsym, exp, log) _, rem gf_poly_div(codeword, gen, exp, log) assert rem [0] * nsym, 编码结果没有通过整除校验 # 自检二原地改 3 个符号模拟信道破坏 corrupted list(codeword) corrupted[7] ^ exp[3] corrupted[40] ^ exp[9] corrupted[128] ^ exp[14]这 3 个符号溢出位随意不需要对齐任何规则。如果自检二最终能解回原消息说明 GF 运算、生成多项式、除法编码这条链路是闭环的。会在这里翻车的基本原因只有一个多项式系数的幂次顺序不一致导致编码器和后面的译码器对“位置”的理解南辕北辙。4. 从校验子到 Berlekamp-MasseyRS 译码的完整实现4.1 校验子计算先看码字有没有被污染RS 译码第一步永远是算校验子。发送端选了 2t 个根 α^0 到 α^(2t-1)接收端把收到的码字依次代入这些根如果所有结果都是 0说明收到的就是合法码字直接输出。否则结果非零的位置就构成了“哪里不对劲”的线索。译码器在这里还是一次通过开销几乎可以忽略。def rs_syndrome(codeword, nsym, exp, log): synd [] for i in range(nsym): # 在 x α^i 处计算收到的码字多项式 y codeword[0] for j in range(1, len(codeword)): y gf_mul(y, exp[i], exp, log) ^ codeword[j] synd.append(y) return synd这里用 Horner 方法计算多项式值省掉每次重新求 x 的幂。synd长度为 32如果全部是 0后面的 Berlekamp-Massey 和 Chien 搜索都不需要执行。实际产品里我会在这里先加一层快速路径校验子全零直接返回原消息省掉至少几微秒对高速链路是实打实的收益。4.2 Berlekamp-Massey 求解错误位置多项式当校验子非零时需要找到一个错误位置多项式 σ(x)满足 σ(x) 的根能精确定位出错的符号位置。Berlekamp-MasseyBM算法用迭代方式不断修正 σ(x)是 RS 译码最核心、也最容易被抄错的一段。下面这个实现以校验子数组和 nsym 为输入返回 σ(x) 的系数按升序排列。def rs_berlekamp_massey(synd, nsym, exp, log): # C 是当前的 σ(x)B 是上次更新的缓存 C [1] [0] * (nsym - 1) B [1] [0] * (nsym - 2) L 0 m 1 b 1 for n in range(nsym): # 计算当前步骤的不一致值 discrepancy d synd[n] for i in range(1, L 1): d ^ gf_mul(C[i], synd[n - i], exp, log) if d 0: m 1 continue T list(C) coef gf_mul(d, gf_inv(b, exp, log), exp, log) for j in range(nsym - m): C[j m] ^ gf_mul(coef, B[j], exp, log) if 2 * L n: L n 1 - L B T b d m 1 else: m 1 return C[:L 1]BM 的直观解释是每迭代一步它都在“当前已掌握的错误位置信息”上做最小修正把所有已知错误信息用一个次数尽量低的多项式表达。算法里的L是当前 σ(x) 的次数b是上一次非零的 discrepancym记录迭代距离。如果你把C[i]和synd[n-i]的下标理解成多项式的卷积关系整段代码就变得很好读每一步都在算预测值和真实值差了多少差了就补一个修正项。4.3 Chien 搜索定位Forney 公式求错误值σ(x) 的根对应的就是出错位置。Chien 搜索不是去解方程而是把 GF(2^8) 里的 255 个非零元素挨个代入 σ(x)谁让 σ(x) 为 0谁就是错误的符号位置。定位之后还需要知道该位置被改成了什么值这步叫 Forney 公式先由校验子和 σ(x) 求错误估计多项式 Ω(x)再代入各错误位置算出幅度。def gf_poly_eval(poly, x, exp, log): # Horner 法计算 poly(x)poly 下标对应 x 的幂次 y 0 for c in reversed(poly): y gf_mul(y, x, exp, log) ^ c return y def rs_find_errors(sigma, n, exp, log): locs [] for i in range(n): # σ(α^{-i}) 0 表示索引 i 处出错 x exp[(255 - (i % 255)) % 255] if gf_poly_eval(sigma, x, exp, log) 0: locs.append(i) return locs求幅值这一步容易写错。GF(2^m) 里的求导和实数域不同偶数次幂系数求导后变成 0奇数次幂系数保留下来原因是在特征为 2 的域里偶数项求导后系数乘以偶数等价于 0。所以 Forney 公式里的 σ(x) 只需要保留奇数下标项。def rs_calc_magnitudes(synd, sigma, locs, nsym, exp, log): # Ω(x) S(x) * σ(x) mod x^nsym omega [0] * nsym for i in range(nsym): z 0 for j in range(i 1): if i - j len(sigma): z ^ gf_mul(synd[j], sigma[i - j], exp, log) omega[i] z values [] for pos in locs: x_inv exp[(255 - (pos % 255)) % 255] omega_x gf_poly_eval(omega, x_inv, exp, log) deriv [sigma[l] if l % 2 1 else 0 for l in range(len(sigma))] deriv_x gf_poly_eval(deriv, x_inv, exp, log) values.append(gf_div(omega_x, deriv_x, exp, log)) return valuesrs_calc_magnitudes里gf_div的实现与前文gf_inv对应这里用到了域除法。如果求出来deriv_x是 0说明 σ(x) 有重根码字损坏到了不可恢复的程度直接判定译码失败比强修更安全。4.4 把步骤收进一个 rs_correct 主函数实际使用时不希望每一步都手动串联我会把它们收到一个函数里并加上两个失败保护的判断σ(x) 的次数超过 t或 Chien 找到的根数量与 σ 次数不一致都意味着错误数量超过纠错能力。出现这些情况时直接报错而不是输出一个“看似正常但实际错误”的码字。def rs_correct(received, nsym, exp, log): synd rs_syndrome(received, nsym, exp, log) if not any(synd): return received sigma rs_berlekamp_massey(synd, nsym, exp, log) if len(sigma) - 1 nsym // 2: raise ValueError(错误个数超过纠错能力 t) locs rs_find_errors(sigma, len(received), exp, log) if len(locs) ! len(sigma) - 1: raise ValueError(错误定位失败码字不可纠) vals rs_calc_magnitudes(synd, sigma, locs, nsym, exp, log) out list(received) for pos, val in zip(locs, vals): out[pos] ^ val return out这段代码里len(sigma)-1就是 σ(x) 的次数也等于它能解释的错误个数。只要 locs 个数对不上就别硬修。把失败显式抛出来比返回一个带残留错误的码字更安全。5. 实测两种错误模型RS 的边界条件与调参思路5.1 用突发错误模型跑一遍随机实验RS 的纠错能力在“错误均匀分散”和“错误连续成片”两种模型下差别很大。我一般先在本地做随机突发试验生成一帧编码随机选一个起点连续破坏若干个符号再送去译码统计成功修复的比例。import random def test_burst(nsym32, burst_len16, rounds1000): exp, log gf_init() ok 0 for _ in range(rounds): msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) start random.randint(0, 255 - burst_len) for i in range(start, start burst_len): bad[i] ^ random.randint(1, 255) try: fixed rs_correct(bad, nsym, exp, log) if fixed cw: ok 1 except ValueError: pass return ok / rounds当burst_len等于或小于 t16 时修复率接近 100%。一旦突发长度超过 16修复率不会立刻归零而是下降成零星的成功——因为某些“突发”虽然跨了很多符号但里面的符号可能受害不深或者错误位置恰好分散。这个实验能直观看到 RS 的硬上限不是概率性的超出的部分就是超了剩下的成功只是运气。5.2 突发超过 t 时的两种惯用补救一种常规做法是加交织器。把多帧码字按行写入、按列读出原本连续 32 个符号的突发就被摊到 4 个不同的 RS 码字里每个码字只损失 8 个符号回到了可纠范围内。另一种是缩短 RS 码把 RS(255,223) 缩短成 RS(200,168)帧长变短突发在帧内占的比例自然下降。这两种方法可以组合很多存储控制器就是这么做的先做一次小交织再用缩短的 RS 码兜底。配置nkt码率典型用途RS(255,223)255223160.874深空、DVB 外码RS(255,239)25523980.937光纤传输RS(255,251)25525120.984高速接口保护RS(32,28)322820.875光盘子码5.3 编码译码顺序错了会怎样通信系统里经常出现收发两端把“高次先发、低次先发”搞反的情况。RS 码本身对线性顺序不敏感但缩短码、交织器、以及前后级协议对符号顺序是敏感的。一个很典型的错误是发送端按升序排列信息符号接收端按降序解释位置结果 Chien 搜索找到的错误位置全部镜像错位译码器会纠结在一堆“解释不通”的校验子上。调试时先打印synd如果看到非零校验子的分布完全没规律先别急着调 BM回头核对帧格式的符号顺序。6. 验证 RS 译码器的三个硬指标误判率、残留错误和时间RS 译码器写完后比“能解通一个例子”更重要的是验证它的边界。我会固定跑三组测试超过 t 个错误时是否报错、修复后的码字是否真的能被 g(x) 整除、单次译码耗时是否稳定。def verify_decoder(): exp, log gf_init() nsym 32 # 第一组0 到 t 个错误必须全部修复 for err_count in range(1, 17): msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) for i in range(err_count): bad[i * 3] ^ 0xb3 assert rs_correct(bad, nsym, exp, log) cw # 第二组超过 t 个错误必须抛异常而不是静默输出 msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) for i in range(20): bad[i] ^ 0x1f try: rs_correct(bad, nsym, exp, log) raise AssertionError(超限错误未被拒绝) except ValueError: pass # 第三组修复后校验子必须全零 msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) bad[10] ^ 0x55 bad[200] ^ 0xaa fixed rs_correct(bad, nsym, exp, log) assert rs_syndrome(fixed, nsym, exp, log) [0] * nsym第二个用例里如果 BM 或 Chien 搜索写错常见的表现不是抛异常而是返回一个“看起来能解、实际残错”的结果最后在第三组测试里原形毕露。所以三组必须一起跑缺一不可。最后一类问题在工程里特别隐蔽GF 表是全局单例还是每次重建直接影响并发场景下的正确性。我的习惯是把exp、log作为参数传入所有函数或者封进一个RSDecoder类避免多线程下共享表被意外覆盖。时间指标上GF(2^8) 的 RS(255,223) 在普通 Python 里单次译码约在几十微秒量级主要开销在 Chien 搜索的 255 次多项式求值如果帧长变成 RS(4095, 3583)m12域表大小翻 16 倍Chien 搜索的循环次数也会线性增加。这时候用 C 扩展或者 SIMD 之前先检查一下是不是把校验子计算写成了 O(n^2) 的重复求幂。查表版本已经是 O(n)别再为它画蛇添足。本文还有配套的精品资源点击获取 RS码从数学到工程：GF(2^8)实现编码与突发错误纠错
简介RS码Reed-Solomon码作为应用广泛的信道编码方案在通信与数据存储领域承担着纠错重任。资源内提供了一套完整的RS码MATLAB实现面向通信工程、电子信息类学生及科研人员可辅助理解伽罗华域运算、RS编码、译码及AWGN信道下的仿真流程。压缩包共18个文件以15个.m函数脚本为主另有3个.asv自动备份文件整体仅7KB代码轻量清晰。文件覆盖GF(2^m)域生成、多项式乘加、BPSK调制映射、AWGN信道模拟、RS编码、Chien搜索与Forney译码等核心模块便于逐段研读与复用。目前已有275人学习下载适合需要结合理论进行仿真实验或完成课程设计的读者。通过运行demo可直观观察RS码在加性高斯白噪声环境下的纠错表现为后续深入通信系统研究打下基础。1. 突发误码与爆浆错误RS 码仍是信道编码的主心骨在真实信道里错误很少是均匀撒开的一两颗“椒盐”更多是连续一段符号整体被打翻。光盘上的划痕、闪存的坏块、无线信号被汽车点火干扰后的持续衰落这类突发错误让只能纠单个比特的码型立刻失去作用。Reed-Solomon 码RS 码恰恰是站在“符号”而不是“比特”的粒度上做编码译码的一个 8 比特符号整体被改掉在 RS 码眼里只是一个错误符号这与信道编码中常见的单比特纠错码是完全不同的纠错模型。也正因为如此RS 码从光盘、QR 码一路活到了 5G 的 eCPRI 接口和深空通信 CCSDS 标准里始终是信道编码里兜底的那一层。这篇文章不讨论某份现成源码而是从 GF(2^m) 的运算法则出发完整写好一个 RS(255,223) 的编码器与译码器再把它们放进突发错误模型里看参数边界在哪里。2. RS 码的数学底盘GF(2^m) 运算与 RS(n,k) 参数解读2.1 伽罗瓦域四个运算加法用 XOR乘法查表RS 码的“符号”不是普通整数而是有限域 GF(2^m) 里的元素。最常见的 m8 意味着一个符号占 8 bitGF(2^8) 里一共有 256 个合法元素每个元素都对应一个 8 bit 多项式。域上的加法就是按位异或这一点比整数加法简单得多乘法则要靠一个本原多项式做模约减。我一般会提前把“指数-元素”和“元素-指数”两张表算好之后所有乘法都变成查表跑 RS 译码时才不会在每次乘除上浪费周期。# 构建 GF(2^8) 的两张表exp 输出 α^ilog 给出元素的对数 def gf_init(prim0x11d, field256): exp [0] * (2 * field) log [0] * field x 1 for i in range(field - 1): exp[i] x log[x] i x 1 # 左移一位相当于乘以生成元 α if x field: # 溢出到第 9 位按本原多项式约减 x ^ prim # 把表扩到两倍后续 exp[log[a]log[b]] 就不必每次都取模 for i in range(field - 1, 2 * field - 1): exp[i] exp[i - (field - 1)] return exp, log def gf_mul(x, y, exp, log): if x 0 or y 0: return 0 return exp[log[x] log[y]] def gf_inv(x, exp, log): if x 0: raise ValueError(0 没有逆元) return exp[255 - log[x]]这里0x11d是 RS 常用的本原多项式 x^8 x^4 x^3 x^2 1。代码里x field判断左移后是否出现了 x^8 这一项有就异或本原多项式相当于做一次模 2 多项式除法。表建好之后gf_mul的代价只有一次数组寻址gf_inv用255 - log[x]是因为非零元素构成 255 阶循环群求逆即找指数互为相反数的元素。这套查表结构是后面生成多项式、校验子计算、Chien 搜索能跑得动的前提。2.2 最小距离与生成多项式的关系RS(n,k) 码的码字被设计成生成多项式 g(x) 的倍式。给定 2t 个校验符号生成多项式写作g(x) (x α^0)(x α^1)...(x α^(2t-1))只要发送端给出的码字是 g(x) 的倍式接收端用同样一组根去验证就能知道哪里不对。这也直接决定了码的最小距离 d_min n - k 1即码字之间最少可以相差多少个符号。RS 是一个最大距离可分码MDS 码这意味着同样冗余下它能达到的 d_min 已经顶到了理论天花板。信道编码选 RS 时看重的正是“已知 t 个错误符号必定能纠”的这个硬保证而不是那种靠概率才纠得动的软保证。2.3 RS(n,k) 参数表t 与码率的折中参数含义RS(255,223) 的取值m每个符号的比特数8GF(2^8)n码字包含的符号数255 2^8 - 1k一个码字里的信息符号数223n-k校验符号数32t可纠错误符号上限(n-k)/2 16d_min码最小距离n-k1 33码率 R有效信息占比k/n 0.8745这里最容易看走眼的是 t (n-k)/2RS 是按符号纠错的所以“t16”表示一帧里最多能容忍 16 个错误符号而不是 16 个错误比特。一个符号在信道里可能有好几个比特同时翻转但只要这 8 个比特同属一个错误符号计数仍然只算 1。正是这个性质让 RS 特别适合突发错误信道——突发再长只要它落在有限个符号内对 RS 译码器来说就跟几个孤立的错误符号没有区别。3. 用 Python 跑通 RS 编码从构造生成多项式到校验位3.1 一次性构造生成多项式 g(x)编码前的准备只有一件事按校验符号个数 nsym 生成 g(x)。因为 g(x) 的根必须包含 α^0 到 α^(2t-1)多项式构造本质上是把这 2t 个一次因式乘起来。GF(2^8) 里的加法和减法等价所以因式直接写成(x α^i)。def gf_poly_mul(p, q, exp, log): r [0] * (len(p) len(q) - 1) for i, a in enumerate(p): if a 0: continue for j, b in enumerate(q): if b 0: continue r[i j] ^ gf_mul(a, b, exp, log) return r def rs_generator_poly(nsym, exp, log): g [1] for i in range(nsym): # 每次乘一个 (x α^i)最后得到 g(x)∏(xα^i) g gf_poly_mul(g, [1, exp[i]], exp, log) return g这个实现里我用“数组下标代表 x 的幂次”的升序约定也就是poly[0]是 x^0 的系数poly[1]是 x^1 的系数。实际通信帧一般按“先发高次符号”处理所以后面对接协议时记得把符号顺序颠倒一下否则发出去的码字长度和位置都对不上。nsym传 32得到的g就是 33 个系数的数组最高次系数始终是 1。3.2 用循环码除法右移得到校验位让码字成为 g(x) 的倍式最经典的做法是系统码把消息多项式先乘以 x^nsym再用 g(x) 做多项式除法余数就是校验符号直接拼在信息符号后面。这套“除以生成多项式、余作校验”的思路是循环码的标准动作。def gf_poly_div(dividend, divisor, exp, log): out list(dividend) for i in range(len(dividend) - len(divisor) 1): coef out[i] if coef ! 0: for j in range(1, len(divisor)): if divisor[j] ! 0: # GF(2^m) 中减法和加法同为异或 out[i j] ^ gf_mul(divisor[j], coef, exp, log) rem_len len(divisor) - 1 return out[:-rem_len], out[-rem_len:] def rs_encode_msg(msg, nsym, exp, log): gen rs_generator_poly(nsym, exp, log) _, rem gf_poly_div(msg [0] * nsym, gen, exp, log) return msg rem注意gf_poly_div返回两个部分我只取余数。msg [0] * nsym相当于把消息乘上 x^nsym因为低次端补零后原有各项的幂次都没变而校验位最终落在低次端。编码完成后一个码字长度为 k nsym前 k 个符号是原始消息后 nsym 个是校验这正是系统码最友好的输出形式。fmt参数不漏写exp和log是第 2 章里gf_init返回的两张表任何编译出错十有八九是漏传了它们。3.3 自检编码器余数为零才算对编码器写完后不要急着调译码先验证一个基本事实任意编码结果除以 g(x)余数必须是全零。再顺手翻转 t 个以内的错误符号用第 4 章的译码器跑回来比对。exp, log gf_init() nsym 32 msg [0x0a] * 223 # 223 个信息符号全部填 0x0a codeword rs_encode_msg(msg, nsym, exp, log) # 自检一能被生成多项式整除 gen rs_generator_poly(nsym, exp, log) _, rem gf_poly_div(codeword, gen, exp, log) assert rem [0] * nsym, 编码结果没有通过整除校验 # 自检二原地改 3 个符号模拟信道破坏 corrupted list(codeword) corrupted[7] ^ exp[3] corrupted[40] ^ exp[9] corrupted[128] ^ exp[14]这 3 个符号溢出位随意不需要对齐任何规则。如果自检二最终能解回原消息说明 GF 运算、生成多项式、除法编码这条链路是闭环的。会在这里翻车的基本原因只有一个多项式系数的幂次顺序不一致导致编码器和后面的译码器对“位置”的理解南辕北辙。4. 从校验子到 Berlekamp-MasseyRS 译码的完整实现4.1 校验子计算先看码字有没有被污染RS 译码第一步永远是算校验子。发送端选了 2t 个根 α^0 到 α^(2t-1)接收端把收到的码字依次代入这些根如果所有结果都是 0说明收到的就是合法码字直接输出。否则结果非零的位置就构成了“哪里不对劲”的线索。译码器在这里还是一次通过开销几乎可以忽略。def rs_syndrome(codeword, nsym, exp, log): synd [] for i in range(nsym): # 在 x α^i 处计算收到的码字多项式 y codeword[0] for j in range(1, len(codeword)): y gf_mul(y, exp[i], exp, log) ^ codeword[j] synd.append(y) return synd这里用 Horner 方法计算多项式值省掉每次重新求 x 的幂。synd长度为 32如果全部是 0后面的 Berlekamp-Massey 和 Chien 搜索都不需要执行。实际产品里我会在这里先加一层快速路径校验子全零直接返回原消息省掉至少几微秒对高速链路是实打实的收益。4.2 Berlekamp-Massey 求解错误位置多项式当校验子非零时需要找到一个错误位置多项式 σ(x)满足 σ(x) 的根能精确定位出错的符号位置。Berlekamp-MasseyBM算法用迭代方式不断修正 σ(x)是 RS 译码最核心、也最容易被抄错的一段。下面这个实现以校验子数组和 nsym 为输入返回 σ(x) 的系数按升序排列。def rs_berlekamp_massey(synd, nsym, exp, log): # C 是当前的 σ(x)B 是上次更新的缓存 C [1] [0] * (nsym - 1) B [1] [0] * (nsym - 2) L 0 m 1 b 1 for n in range(nsym): # 计算当前步骤的不一致值 discrepancy d synd[n] for i in range(1, L 1): d ^ gf_mul(C[i], synd[n - i], exp, log) if d 0: m 1 continue T list(C) coef gf_mul(d, gf_inv(b, exp, log), exp, log) for j in range(nsym - m): C[j m] ^ gf_mul(coef, B[j], exp, log) if 2 * L n: L n 1 - L B T b d m 1 else: m 1 return C[:L 1]BM 的直观解释是每迭代一步它都在“当前已掌握的错误位置信息”上做最小修正把所有已知错误信息用一个次数尽量低的多项式表达。算法里的L是当前 σ(x) 的次数b是上一次非零的 discrepancym记录迭代距离。如果你把C[i]和synd[n-i]的下标理解成多项式的卷积关系整段代码就变得很好读每一步都在算预测值和真实值差了多少差了就补一个修正项。4.3 Chien 搜索定位Forney 公式求错误值σ(x) 的根对应的就是出错位置。Chien 搜索不是去解方程而是把 GF(2^8) 里的 255 个非零元素挨个代入 σ(x)谁让 σ(x) 为 0谁就是错误的符号位置。定位之后还需要知道该位置被改成了什么值这步叫 Forney 公式先由校验子和 σ(x) 求错误估计多项式 Ω(x)再代入各错误位置算出幅度。def gf_poly_eval(poly, x, exp, log): # Horner 法计算 poly(x)poly 下标对应 x 的幂次 y 0 for c in reversed(poly): y gf_mul(y, x, exp, log) ^ c return y def rs_find_errors(sigma, n, exp, log): locs [] for i in range(n): # σ(α^{-i}) 0 表示索引 i 处出错 x exp[(255 - (i % 255)) % 255] if gf_poly_eval(sigma, x, exp, log) 0: locs.append(i) return locs求幅值这一步容易写错。GF(2^m) 里的求导和实数域不同偶数次幂系数求导后变成 0奇数次幂系数保留下来原因是在特征为 2 的域里偶数项求导后系数乘以偶数等价于 0。所以 Forney 公式里的 σ(x) 只需要保留奇数下标项。def rs_calc_magnitudes(synd, sigma, locs, nsym, exp, log): # Ω(x) S(x) * σ(x) mod x^nsym omega [0] * nsym for i in range(nsym): z 0 for j in range(i 1): if i - j len(sigma): z ^ gf_mul(synd[j], sigma[i - j], exp, log) omega[i] z values [] for pos in locs: x_inv exp[(255 - (pos % 255)) % 255] omega_x gf_poly_eval(omega, x_inv, exp, log) deriv [sigma[l] if l % 2 1 else 0 for l in range(len(sigma))] deriv_x gf_poly_eval(deriv, x_inv, exp, log) values.append(gf_div(omega_x, deriv_x, exp, log)) return valuesrs_calc_magnitudes里gf_div的实现与前文gf_inv对应这里用到了域除法。如果求出来deriv_x是 0说明 σ(x) 有重根码字损坏到了不可恢复的程度直接判定译码失败比强修更安全。4.4 把步骤收进一个 rs_correct 主函数实际使用时不希望每一步都手动串联我会把它们收到一个函数里并加上两个失败保护的判断σ(x) 的次数超过 t或 Chien 找到的根数量与 σ 次数不一致都意味着错误数量超过纠错能力。出现这些情况时直接报错而不是输出一个“看似正常但实际错误”的码字。def rs_correct(received, nsym, exp, log): synd rs_syndrome(received, nsym, exp, log) if not any(synd): return received sigma rs_berlekamp_massey(synd, nsym, exp, log) if len(sigma) - 1 nsym // 2: raise ValueError(错误个数超过纠错能力 t) locs rs_find_errors(sigma, len(received), exp, log) if len(locs) ! len(sigma) - 1: raise ValueError(错误定位失败码字不可纠) vals rs_calc_magnitudes(synd, sigma, locs, nsym, exp, log) out list(received) for pos, val in zip(locs, vals): out[pos] ^ val return out这段代码里len(sigma)-1就是 σ(x) 的次数也等于它能解释的错误个数。只要 locs 个数对不上就别硬修。把失败显式抛出来比返回一个带残留错误的码字更安全。5. 实测两种错误模型RS 的边界条件与调参思路5.1 用突发错误模型跑一遍随机实验RS 的纠错能力在“错误均匀分散”和“错误连续成片”两种模型下差别很大。我一般先在本地做随机突发试验生成一帧编码随机选一个起点连续破坏若干个符号再送去译码统计成功修复的比例。import random def test_burst(nsym32, burst_len16, rounds1000): exp, log gf_init() ok 0 for _ in range(rounds): msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) start random.randint(0, 255 - burst_len) for i in range(start, start burst_len): bad[i] ^ random.randint(1, 255) try: fixed rs_correct(bad, nsym, exp, log) if fixed cw: ok 1 except ValueError: pass return ok / rounds当burst_len等于或小于 t16 时修复率接近 100%。一旦突发长度超过 16修复率不会立刻归零而是下降成零星的成功——因为某些“突发”虽然跨了很多符号但里面的符号可能受害不深或者错误位置恰好分散。这个实验能直观看到 RS 的硬上限不是概率性的超出的部分就是超了剩下的成功只是运气。5.2 突发超过 t 时的两种惯用补救一种常规做法是加交织器。把多帧码字按行写入、按列读出原本连续 32 个符号的突发就被摊到 4 个不同的 RS 码字里每个码字只损失 8 个符号回到了可纠范围内。另一种是缩短 RS 码把 RS(255,223) 缩短成 RS(200,168)帧长变短突发在帧内占的比例自然下降。这两种方法可以组合很多存储控制器就是这么做的先做一次小交织再用缩短的 RS 码兜底。配置nkt码率典型用途RS(255,223)255223160.874深空、DVB 外码RS(255,239)25523980.937光纤传输RS(255,251)25525120.984高速接口保护RS(32,28)322820.875光盘子码5.3 编码译码顺序错了会怎样通信系统里经常出现收发两端把“高次先发、低次先发”搞反的情况。RS 码本身对线性顺序不敏感但缩短码、交织器、以及前后级协议对符号顺序是敏感的。一个很典型的错误是发送端按升序排列信息符号接收端按降序解释位置结果 Chien 搜索找到的错误位置全部镜像错位译码器会纠结在一堆“解释不通”的校验子上。调试时先打印synd如果看到非零校验子的分布完全没规律先别急着调 BM回头核对帧格式的符号顺序。6. 验证 RS 译码器的三个硬指标误判率、残留错误和时间RS 译码器写完后比“能解通一个例子”更重要的是验证它的边界。我会固定跑三组测试超过 t 个错误时是否报错、修复后的码字是否真的能被 g(x) 整除、单次译码耗时是否稳定。def verify_decoder(): exp, log gf_init() nsym 32 # 第一组0 到 t 个错误必须全部修复 for err_count in range(1, 17): msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) for i in range(err_count): bad[i * 3] ^ 0xb3 assert rs_correct(bad, nsym, exp, log) cw # 第二组超过 t 个错误必须抛异常而不是静默输出 msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) for i in range(20): bad[i] ^ 0x1f try: rs_correct(bad, nsym, exp, log) raise AssertionError(超限错误未被拒绝) except ValueError: pass # 第三组修复后校验子必须全零 msg [random.randint(0, 255) for _ in range(223)] cw rs_encode_msg(msg, nsym, exp, log) bad list(cw) bad[10] ^ 0x55 bad[200] ^ 0xaa fixed rs_correct(bad, nsym, exp, log) assert rs_syndrome(fixed, nsym, exp, log) [0] * nsym第二个用例里如果 BM 或 Chien 搜索写错常见的表现不是抛异常而是返回一个“看起来能解、实际残错”的结果最后在第三组测试里原形毕露。所以三组必须一起跑缺一不可。最后一类问题在工程里特别隐蔽GF 表是全局单例还是每次重建直接影响并发场景下的正确性。我的习惯是把exp、log作为参数传入所有函数或者封进一个RSDecoder类避免多线程下共享表被意外覆盖。时间指标上GF(2^8) 的 RS(255,223) 在普通 Python 里单次译码约在几十微秒量级主要开销在 Chien 搜索的 255 次多项式求值如果帧长变成 RS(4095, 3583)m12域表大小翻 16 倍Chien 搜索的循环次数也会线性增加。这时候用 C 扩展或者 SIMD 之前先检查一下是不是把校验子计算写成了 O(n^2) 的重复求幂。查表版本已经是 O(n)别再为它画蛇添足。本文还有配套的精品资源点击获取