RM码编译码仿真从零实现：生成矩阵与Hadamard译码

前阵子做通信仿真练习，需要选一种线性分组码来对比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)。

下面是几个常用参数组合的速算表：

观察这个表能得出一个很重要的直觉：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^m，k = 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, dtype=int)

这段代码的核心逻辑就是三重循环：先按次数分组，再枚举具体变量组合，最后对每个输入点求值。注意combinations(range(m), degree)在degree=0时会生成一个空元组，对应常数项，行为刚好是全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, u+v)属于RM(r, m)，其中u+v是逐位异或。

这个构造乍看只是数学上的一种等价描述，但它对理解RM码的递归译码至关重要。你去读RM码相关的论文，会发现大部分迭代译码算法都是基于这种(u, u+v)分解来做的。先译后半部分得到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, dtype=int)

这里有一个细节必须强调：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, u+v)分解的递归软判决译码，或者把RM码当作某种因子图跑置信传播。这些方法我这次没有展开成完整实现，但你要知道它们的存在——大数逻辑只是RM码译码的第一课，不是终点。

4. 仿真平台搭建与误码率实测

4.1 Eb/N0与噪声方差的换算

误码率仿真的第一步是把信道信噪比转换成AWGN噪声方差，这一步单位搞错，所有曲线都会崩。我用的BPSK映射是0到+1、1到-1，符号能量Es = 1。对码率为R的编码系统，Eb/N0和Es/N0的关系是：

Es/N0 = R * Eb/N0

AWGN的噪声方差σ²和N0的关系是σ² = N0 / 2，于是：

σ = sqrt(1 / (2 * R * 10^(EbN0_dB / 10)))

def simulate_ber(r, m, ebn0_db, max_errors=300, max_trials=200000): 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个百分点，性能差异很直观。

数据符合理论预期。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_errors=300，意思是每个信噪比点至少要积累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, u+v)}的分解，把长码译码拆成两个短码译码，中间用软信息往返传递，这是理解现代迭代译码的最好入口。

第二个方向是RM码与极化码的关系。极化码的构造可以理解为从RM码的递推结构中挑选更可靠的子信道，用RM码做基线，对比不同信道极化方案的性能差异，是个很有意思的实验课题。

第三个方向是有限码长下的性能界分析。RM码的纠错能力可以用球半径界估计，但实际最大似然译码的性能往往优于这个界。用小参数RM码跑最大似然译码，再和理论界对比，能直观感受到“距离界”和“真实性能”之间的空隙。

做完这套RM编译码仿真，我最大的体会是：有结构的码比随机生成的码更容易理解，也更容易实现。RM码的全部秘密都在那张生成矩阵里——你只要知道m个布尔变量和多项式次数，就能推导出码长、距离、译码方式，甚至预测它的性能走向。这种透明度在信道编码里并不常见，遇到它值得认真玩一玩。后面再做别的编码仿真，我基本都保留了一个习惯：先用小参数把所有环节闭环跑通，确认无误码平台和译码映射正确，再上完整蒙特卡洛仿真。这套流程虽然慢一点，但省下的排错时间远多于多花的编码时间。