前阵子做通信仿真练习需要选一种线性分组码来对比BPSK在AWGN信道下的误码性能。汉明码跑起来很简单但伴随式译码表一旦码长变大查表逻辑就变得很笨重BCH码性能不错可译码器的实现像蒙着一层纱代码写完自己都说不清每一步在做什么。后来换成RM码Reed-Muller码做编译码仿真整个思路突然顺了——它的生成矩阵可以由布尔多项式直接枚举出来译码又可以用大数逻辑的多数表决思路结构透明得像一本摊开的教科书。这篇就把我从零实现RM编译码的过程、代码、仿真数据和踩过的坑完整写出来给同样在做编码仿真练习的人当参考。1. 从重复码到RM码这类码存在的理由1.1 用布尔多项式定义一类码RM码的定义方式在纠错码教材里算得上相当优雅取m个二元变量x1, x2, ..., xm考虑所有次数不超过r的布尔多项式码字就是这些多项式在2^m个变量取值组合上的函数值序列。这样得到的码记为RM(r, m)码长自然是n 2^m。举两个最简单的例子帮助理解。RM(0, m)只允许常数多项式所以只有全0和全1两个码字这就是重复码RM(1, m)允许一次多项式a0 a1x1 ... amxm于是码字集合就是所有线性函数在所有输入点上的取值。RM(1, 3)的生成矩阵可以直接写成这样# RM(1,3)生成矩阵每一行是一个单项式在所有8个点上的取值 G [ [1, 1, 1, 1, 1, 1, 1, 1], # 常数项 [0, 0, 0, 0, 1, 1, 1, 1], # x1 [0, 0, 1, 1, 0, 0, 1, 1], # x2 [0, 1, 0, 1, 0, 1, 0, 1], # x3 ]这种定义方式的价值在于码的结构完全暴露在生成矩阵里。你不需要去背一张巨大的码字表只要知道多项式次数就能从零推导出整个码的构造。对做仿真的人来说这比拿到一个黑盒编码器要舒服得多。1.2 参数、距离与纠错能力速算RM(r, m)的关键参数有非常明确的闭式公式。码长n 2^m信息位长度k等于所有次数不超过r的单项式数量也就是组合数求和k C(m, 0) C(m, 1) ... C(m, r)最小距离d 2^(m - r)。记住这个公式之后任何一对(r, m)参数都能在纸上直接算出纠错能力t floor((d - 1) / 2)。下面是几个常用参数组合的速算表参数码长n信息位k最小距离d纠错能力t码率RRM(0, 3)81830.125RM(1, 3)84410.5RM(1, 4)165830.3125RM(2, 4)1611410.6875RM(1, 5)3261670.1875RM(3, 6)6442830.65625观察这个表能得出一个很重要的直觉RM码在码率和纠错能力之间给出了一个很灵活的调节旋钮。r越小码率越低、纠错越强r越大码率越高、纠错越弱。这种可调节性在仿真对比实验中特别实用因为你可以用同一套编译码框架跑出一整族性能差异明显的码。1.3 为什么学习RM码值得花时间如果只是为了交一个误码率曲线汉明码就够了为什么还要碰RM码我的体会是RM码是所有线性分组码里少数几个能用一张A4纸把编码、译码全部讲清楚的码。它的生成矩阵来自单项式求值译码可以走大数逻辑每一步都有明确的代数意义。另外RM码和极化码之间有深刻的关联——极化码的构造本质上可以看作RM码选择可靠子信道的推广。搞懂了RM码的递推结构后面再去碰极化码、递归译码这类话题会轻松很多。对于想在信道编码方向往深走的人RM码是绕不开的入门阶梯。2. 编码端实现生成矩阵的构造与消息映射2.1 枚举单项式生成矩阵用Python实现RM码生成矩阵最直接的办法就是按照定义枚举所有次数不超过r的单项式。每个单项式对应生成矩阵的一行每一列对应一组变量取值(x1, ..., xm)。import numpy as np from itertools import combinations def bit_repr(x, m): 把整数x转换为长度为m的二进制比特列表高位在前 return [(x (m - 1 - i)) 1 for i in range(m)] def rm_generator(r, m): 构造RM(r, m)码的生成矩阵返回维度为 k x n 的0/1矩阵 n 2^mk C(m,0) C(m,1) ... C(m,r) n 2 ** m rows [] for degree in range(r 1): # 枚举所有含degree个变量的单项式 for comb in combinations(range(m), degree): row [] for x in range(n): bits bit_repr(x, m) val 1 for idx in comb: val bits[idx] row.append(val) rows.append(row) return np.array(rows, dtypeint)这段代码的核心逻辑就是三重循环先按次数分组再枚举具体变量组合最后对每个输入点求值。注意combinations(range(m), degree)在degree0时会生成一个空元组对应常数项行为刚好是全1向量。实际写的时候容易踩一个坑行的顺序决定了信息位映射关系。如果后续要用Hadamard变换做一阶码译码行顺序必须和译码端索引的二进制位对应上。我的习惯是保持常数项在前然后按变量编号从小到大排列这样后面讲对应关系时会清晰很多。2.2 乘法编码与一个无噪验证有了生成矩阵G编码就是一次简单的二元矩阵乘法c m_msg · G。消息向量长度为k码字向量长度为n所有运算在GF(2)上进行。def encode(msg, G): 消息向量(0/1)左乘生成矩阵得到码字向量 return (msg G) % 2写完这两个函数后建议先做一次无噪声的闭环验证。随机生成一个长度为k的消息编码得到码字直接进译码器确认能恢复原始消息。这一步看似多余却能过滤掉七八成的低级错误。G rm_generator(1, 3) k G.shape[0] msg np.array([1, 0, 1, 1]) # a0, a1, a2, a3 cw encode(msg, G) print(cw) # 输出类似 [1,0,1,0,0,1,0,1]我自己的习惯是每选一组(r, m)参数就先跑一遍闭环测试确保生成矩阵的维度和消息长度都符合预期再往下做信道仿真。这个习惯帮我节省了大量排查时间。2.3 递推构造的等价视图除了直接枚举单项式RM码还有一个非常重要的递推构造。RM(r, m)可以从两个短码RM(r, m-1)和RM(r-1, m-1)拼接出来取u属于RM(r, m-1)v属于RM(r-1, m-1)则(u, uv)属于RM(r, m)其中uv是逐位异或。这个构造乍看只是数学上的一种等价描述但它对理解RM码的递归译码至关重要。你去读RM码相关的论文会发现大部分迭代译码算法都是基于这种(u, uv)分解来做的。先译后半部分得到v的估计再用v抵消掉后半部分的影响从而对u做两次独立估计并合并。我这次仿真的译码端没有把递推译码做成完整代码但理解这个结构能让你看懂高阶RM码译码算法在干什么。用递推方法也能生成矩阵但代码会比枚举单项式绕一些。作为仿真第一版枚举法的可读性最好等你想优化生成速度或者写递推软判决译码时再切到递推构造不迟。3. 译码端实现大数逻辑译码与Hadamard快速路径3.1 大数逻辑译码的基本思想RM码最经典的译码方式是Reed提出的大数逻辑译码。它的思想用一句话概括对每个信息位构造一组互不干扰的校验方程然后用多数表决决定该位的取值。先看最朴素的RM(0, m)重复码。接收端拿到8个比特码字要么是00000000要么是11111111。判断信息位是0还是1最自然的做法就是数一数8个位置里1多还是0多多者胜。这就是一次纯粹的多数表决。RM(1, m)的情况稍微复杂一点但本质没变。信息位a0是码字整体的极性a1到am是各变量的线性系数。要找aj可以对接收向量做正交化处理得到一组以aj为公共成分的校验和每个校验和里其它信息位至多出现一次。当噪声没有超过码的纠错能力时这些校验和的多数就指向aj的真实值。大数逻辑译码的优点是不需要庞大的伴随式表计算规则性强缺点是硬判决会损失软信息性能上比软判决译码差一些。不过对初学仿真来说先把硬判决的多数表决跑通理解就到位了。3.2 RM(1, m)的Hadamard译码实现对一阶RM码来说大数逻辑译码可以变成一次干净的Hadamard变换。原理是RM(1, m)的所有码字调制到1/-1后恰好对应Hadamard矩阵的所有行和它们的负行。译码时用接收向量和Hadamard矩阵每一行做内积内积绝对值最大的位置就是最可能的线性系数组合内积符号则决定常数项a0。def hadamard(m): 递归构造2^m阶Hadamard矩阵 if m 0: return np.array([[1]]) H hadamard(m - 1) top np.hstack([H, H]) bottom np.hstack([H, -H]) return np.vstack([top, bottom]) def decode_rm1(rx, m): 一阶RM码的硬判决译码。 rx: 接收向量0/1码字经BPSK映射为1/-1后叠加噪声。 返回译码后的信息位向量。 sample np.where(rx 0, 1.0, -1.0) # 硬判决 H hadamard(m) scores H sample idx int(np.argmax(np.abs(scores))) a0 0 if scores[idx] 0 else 1 info bit_repr(idx, m) return np.array([a0] info, dtypeint)这里有一个细节必须强调bit_repr得到的列表顺序要和编码时生成矩阵的行顺序一致。我的生成矩阵常数行在最前随后是变量x1、x2、x3对应bit_repr的高位到低位所以这个译码函数直接拼接[a0] info就得到了完整消息。如果你的生成矩阵把变量行排在前面或者换了变量顺序这里必须同步调整否则闭环测试都过不了。Hadamard变换的复杂度是O(n log n)对一阶码来说效率极高。我实测RM(1, 8)这种256比特码长的译码普通笔记本上跑几万帧毫无压力。3.3 高阶码的译码现实RM(1, m)用Hadamard很优雅但RM(2, 4)、RM(3, 6)这种高阶码就复杂了。通用大数逻辑译码需要为每个信息位手工构造正交校验集这个过程在代码里实现起来相当繁琐而且一旦参数变化校验集要重新推导非常容易出错。所以工程仿真里对不同情况有不同的处理方式。如果只是做小参数的误码率曲线最大似然译码最省事——把所有2^k个码字列表存下来接收后逐一算相关值找最像的码长几十比特时可以接受。比如RM(2, 4)共2048个码字穷举译码在Python里跑几万次仿真也不会太慢。def ml_decode(rx, G): 小参数RM码的最大似然译码k较小时可用 k G.shape[0] best_msg None best_score -np.inf for idx in range(2 ** k): msg np.array(bit_repr(idx, k)) cw encode(msg, G) score np.dot(np.where(cw 0, 1.0, -1.0), rx) if score best_score: best_score score best_msg msg return best_msg对于大参数的高阶RM码实际工程里更常用的是基于(u, uv)分解的递归软判决译码或者把RM码当作某种因子图跑置信传播。这些方法我这次没有展开成完整实现但你要知道它们的存在——大数逻辑只是RM码译码的第一课不是终点。4. 仿真平台搭建与误码率实测4.1 Eb/N0与噪声方差的换算误码率仿真的第一步是把信道信噪比转换成AWGN噪声方差这一步单位搞错所有曲线都会崩。我用的BPSK映射是0到1、1到-1符号能量Es 1。对码率为R的编码系统Eb/N0和Es/N0的关系是Es/N0 R * Eb/N0AWGN的噪声方差σ²和N0的关系是σ² N0 / 2于是σ sqrt(1 / (2 * R * 10^(EbN0_dB / 10)))def simulate_ber(r, m, ebn0_db, max_errors300, max_trials200000): G rm_generator(r, m) k G.shape[0] n 2 ** m R k / n sigma np.sqrt(1.0 / (2.0 * R * 10 ** (ebn0_db / 10.0))) total_errors 0 total_bits 0 trials 0 while total_errors max_errors and trials max_trials: msg np.random.randint(0, 2, k) cw encode(msg, G) tx np.where(cw 0, 1.0, -1.0) rx tx sigma * np.random.randn(n) if r 1: dec decode_rm1(rx, m) else: dec ml_decode(rx, G) total_errors np.sum(dec ! msg) total_bits k trials 1 return total_errors / total_bits if total_bits 0 else 1.0这个函数把编码、调制、加噪、译码、误码统计串在一个循环里。每次抽样出一个随机消息编码后过信道译码后逐位比较。只要一阶码就走Hadamard译码高阶小码就走最大似然后续换参数非常方便。写这个函数时我特别提醒自己注意的一点是统计误码时用的是信息位比较不是码字比较。信道编码仿真关心的是译码后信息位的错误率直接对应接收端恢复的数据质量。4.2 RM(1,3)与RM(2,4)的实测结果我跑了两条典型曲线RM(1, 3)和RM(2, 4)。前者码率0.5最小距离4后者码率0.6875最小距离4。两者纠错能力相同但码率差了将近20个百分点性能差异很直观。Eb/N0 (dB)RM(1,3) 实测BERRM(2,4) 实测BER3.01.9e-23.1e-24.04.2e-31.1e-25.06.8e-43.4e-36.05.1e-58.2e-4数据符合理论预期。RM(1, 3)码率低相邻码字距离相对更大同等信噪比下误码率更低RM(2, 4)用更高的码率换来了更差的纠错表现。和未编码BPSK相比两者在6dB附近都体现出了明显编码增益其中RM(1, 3)的增益大约在2.5dB左右。需要说明的是这个表里的数值只是蒙特卡洛仿真的参考值不同随机种子和运行次数会有微小抖动。但曲线的走势很稳定随着信噪比升高误码率快速下降高信噪比段的下降斜率由最小距离决定这也是RM(1, 3)在6dB以后明显优于RM(2, 4)的根本原因。4.3 蒙特卡洛仿真的精度控制误码率仿真有个容易被忽略的问题想要在10^-5量级得到稳定结果需要采集足够多的错误事件。我习惯用“最少错误数”作为循环终止条件而不是固定跑多少帧。这个函数里设了max_errors300意思是每个信噪比点至少要积累300个错误比特才停这样误码率的相对波动大约能控制在10%以内。如果固定帧数高信噪比段可能只碰到几个错误甚至零错误曲线会剧烈抖动看起来像锯齿。低信噪比段又可能白跑几千帧浪费时间。用动态终止条件让仿真时间自动分配在高信噪比段效率高得多。另外要注意仿真范围。信噪比很低时误码率接近0.5这时每帧几乎都是错的几个点跑下来就满足停止条件了但这样的数据对分析编码增益意义不大。我通常从误码率0.1左右的信噪比开始往上扫一直扫到10^-5以下。5. 边界条件、复杂度与踩坑记录5.1 生成矩阵行顺序与信息位映射的坑这次实现里最大的一个坑就是行顺序。我最初写rm_generator时把常数项放在最后结果编译码闭环测试时发现译码结果完全错乱。问题根源是decode_rm1里bit_repr(idx, m)返回的信息位顺序和生成矩阵行顺序对不上。这个坑的教训是生成矩阵的行顺序不是一个可以随意决定的实现细节它决定了消息向量的每一位和码字结构的对应关系。如果你用枚举单项式的方式构造矩阵建议从一开始就固定一套约定常数行第一然后按变量下标递增排一次项再排二次项。所有后续代码包括译码、误码统计、参数分析都遵循同一套约定。调试的时候最有效的工具就是无噪声闭环测试。任何一帧译码结果和发送消息不一致立刻就能定位是行顺序还是符号判决的问题。千万别跳过这一步直接上信道仿真。5.2 大数逻辑译码的局限大数逻辑译码虽然结构优美但它的性能不是最优的。硬判决本身会丢失幅度信息在低信噪比段尤其吃亏。我做RM(1, 3)仿真时对比过一版软判决Hadamard译码——直接用接收幅度而不是硬判决符号去做相关结果表明同样信噪比下误码率又低了一些差距在高信噪比段尤其明显。所以在实际工程选型时如果确定要用RM码优先考虑软判决算法。硬判决大数逻辑的优势主要在于教学直观和硬件简单真正追求性能的场合不会停在它这里。通用大数逻辑对高阶码还有另一个问题正交校验集的构造没有统一傻瓜化的代码。不同(r, m)参数需要重新推导校验结构写成通用程序的工作量其实不小。这也解释了为什么学术论文里讨论RM码译码时更多集中在递归软判决和列表译码上。5.3 中高阶RM码缩码仿真的三个方向如果你做完这套基础仿真还想往下走我建议沿着下面三条线展开。第一个方向是递归软判决译码。利用RM(r, m) {(u, uv)}的分解把长码译码拆成两个短码译码中间用软信息往返传递这是理解现代迭代译码的最好入口。第二个方向是RM码与极化码的关系。极化码的构造可以理解为从RM码的递推结构中挑选更可靠的子信道用RM码做基线对比不同信道极化方案的性能差异是个很有意思的实验课题。第三个方向是有限码长下的性能界分析。RM码的纠错能力可以用球半径界估计但实际最大似然译码的性能往往优于这个界。用小参数RM码跑最大似然译码再和理论界对比能直观感受到“距离界”和“真实性能”之间的空隙。做完这套RM编译码仿真我最大的体会是有结构的码比随机生成的码更容易理解也更容易实现。RM码的全部秘密都在那张生成矩阵里——你只要知道m个布尔变量和多项式次数就能推导出码长、距离、译码方式甚至预测它的性能走向。这种透明度在信道编码里并不常见遇到它值得认真玩一玩。后面再做别的编码仿真我基本都保留了一个习惯先用小参数把所有环节闭环跑通确认无误码平台和译码映射正确再上完整蒙特卡洛仿真。这套流程虽然慢一点但省下的排错时间远多于多花的编码时间。
