线性分组码编译码原理与Python仿真:从汉明距离到BER曲线

线性分组码编译码原理与Python仿真:从汉明距离到BER曲线 简介面向通信工程、电子信息及相关专业的学生、教师和科研初学者这份资源以MATLAB仿真形式完整呈现线性分组码的编译码链路旨在解决线性分组码原理理解困难、误码率性能评估缺少直观实验数据的问题既适合课堂辅助教学也适合课程设计或科研入门时的快速验证。压缩包内共8个文件其中包含7个功能独立的.m脚本与1个.mat参数文件脚本分工覆盖生成矩阵构建、信道模拟、译码算法以及误码率统计等关键环节参数文件提供预设仿真环境整体体积仅7KB结构紧凑且易于修改方便读者按需调用或重新组合仿真流程。该资源已有215人浏览学习虽然文件体量较小但代码逻辑完整能够支撑从编码到性能分析的完整实操同时脚本命名清晰、注释得当可按照功能模块独立学习也可以串联起来完成整体实验。通过运行这些脚本并调整信噪比等参数可以比较不同线性分组码在相同信道条件下的纠错表现观察码长、码率等因素对误码率曲线的影响从而加深对线性码代数结构及编解码策略的理解对于初学者可先借助参数文件运行基础仿真再逐步拓展到不同码型的性能对比降低入门难度。最终可输出直观的误码率统计结果为通信系统纠错编码方案的设计提供数据参考与排错思路。1. 线性分组码的纠错边界最小距离决定一切做信道编码仿真时最容易踩进去的直觉陷阱是码字越长纠错能力越强。实际上线性分组码能纠几个错误只由码字集合之间的最小汉明距离 dmin 决定dmin3 的 (7,4) 汉明码能纠 1 位错dmin5 的短码能纠 2 位错与码长并没有直接关系。下面这套方案把线性分组码编译码的完整闭环拆开来讲编码走生成矩阵译码走伴随式标准阵查表性能用 BPSKAWGN 信道下的 BER 曲线验证。读完你能跑通一套可复现的 Python 仿真弄懂 Eb/N0 到噪声方差的换算、错误帧统计和结果校准这些工程里最容易被忽略的环节。适合快速验证分组码方案的通信工程师也适合正在补信道编码实验的学生。2. 线性分组码编码实现生成矩阵、系统码与 H 矩阵校验2.1 编码方程与行空间约束线性分组码的编码方程只有一行c uG mod 2。k 位信息向量 u 左乘 k×n 的生成矩阵 G得到 n 位码字 c。线性二字体现在模 2 加法上任意两个码字相加仍是码字全部 2^k 个码字构成 n 维二元向量空间中的一个 k 维子空间。这意味着 G 的行向量必须线性无关否则不同信息向量会映射到相同码字接收端无法区分发送的是哪一条消息。生成矩阵的选取几乎决定了一个分组码的全部性能。对于同样的 (n,k)行空间可以有很多种取法每种取法对应不同的 dminBER 曲线高信噪比区的下降斜率也因此不同。手工构造时最常用的是循环移位和系统化拼接两种前者借用多项式环生成已知的好码后者直接控制编码结构。2.2 系统码把 G 写成 (I_k | P) 为什么省心系统形式 G(I_k | P) 的含义是码字前 k 位直接复制信息位后 rn-k 位是校验位。校验位由 P 决定parity uP mod 2。系统码不是理论必须却有实打实的工程便利接收端截取前 k 位即可预览信息调试时能直观看到某个校验方程由哪几位信息位参与形成。若你拿到的 G 不是系统形式对 G 做行初等变换可化为行最简形但若涉及列置换码字比特与物理信道位置的对应关系会改变查表译码时需要同步置换回原始坐标容易出错。我的做法是构造时就按系统码来绕开置换步骤。2.3 校验矩阵与 G 的一致性验证给定 G(I_k|P)校验矩阵取 H(P^T|I_r)。H 的每一行是一个奇偶校验方程核心性质是 G H^T 0模 2 意义下。接收端对任意码字 c 计算 c H^T结果必须为全零结果非零说明噪声翻转了某一位或某几位。用系统码构造时(7,4) 汉明码的经典参数为P [[1,1,0],[1,1,1],[1,0,1],[0,1,1]]对应 dmin3。选择 P 的总原则是让 G 行线性无关且 dmin 尽可能大对同参数对暴力搜索 P 是已知的组合爆炸问题实际工程中直接采用查到的经典码如汉明码、扩展汉明码更快。2.4 编码与校验的 numpy 最小实现下面的函数不依赖任何通信工具箱适合作为整个仿真链路的起点。import numpy as np def systematic_encode(u, P): 系统码编码c u (I_k | P) mod 2 参数: u: 形状 (k,) 或 (B, k) 的二进制信息向量 P: 形状 (k, r) 的校验分量矩阵 返回: 系统码码字形状 (k,) 或 (B, n)前 k 位为信息位 u np.asarray(u) single (u.ndim 1) if single: u u[None, :] parity np.mod(u P, 2) c np.concatenate([u, parity], axis1) return c[0] if single else c def compute_H(P): 由校验分量 P 构造校验矩阵 H (P^T | I_r) k, r P.shape H np.concatenate([P.T, np.eye(r, dtypeint)], axis0) return H.astype(int) def check_codeword(c, H): 返回全零表示 c 是合法码字 return np.mod(H c, 2)np.mod(u P, 2)完成二元域上的乘法与异或P 和 H 的 dtype 必须设成 int避免浮点矩阵乘法产生的取整误差。下面的使用示例可以直接跑P np.array([[1, 1, 0], [1, 1, 1], [1, 0, 1], [0, 1, 1]], dtypeint) H compute_H(P) u np.array([1, 0, 1, 1]) c systematic_encode(u, P) # [1 0 1 1 0 1 1] print(check_codeword(c, H)) # [0 0 0]参数说明P的每一行对应一个信息位每一列对应一个校验方程行数和列数任意设置但 G 的秩必须保持满秩。check_codeword返回三个零就是编码正确返回非零说明 c 不是码字问题多半出在 P 的选取导致行相关。常见误用是直接把浮点软信息喂给check_codeword模 2 会把非零软量折叠成 1结果伴随式失真。软信息必须先硬判决或走第 4 章讨论的软输入译码路线。下面列出几组常见短码的参考参数方便直接拿去对比仿真(n, k)rn-kdmin纠错能力 t典型场景(7,4) 汉明331信道编码入门仿真(8,4) 扩展汉明441可检 2 错纠检结合的小型帧(15,11) 汉明431存储系统 ECC 常见配置(31,26) 汉明531低速控制链路这些码的 dmin 相同瀑布区斜率接近交叉点位置则由 Rk/n 主导这个结论在第 4 章的曲线里能直接看到。3. 线性分组码伴随式硬判决译码标准阵与查表纠错3.1 伴随式到底在算什么接收向量 r_b 送入译码器之前先做硬判决。译码的目标是估计错误图案 e使 c_hat r_b xor e 是合法码字且 c_hat 与真实码字距离最近。伴随式 s H r_b^T把它展开若 r_bces H c^T H e^T H e^T可见 s 只由 e 决定与发送的具体码字无关。因此全系统只需一张从 s 到 e 的映射表2^r 个条目。这里的关键是最可能三个字怎么落地等概信息位下独立加性噪声导致重量较小的 e 概率更高所以查表策略等价于找最小重量错误图案。一个接收向量对应一个 s而所有与 s 相同伴随式的错误图案构成一个陪集每个陪集里重量最小的向量就是陪集首。标准阵正是这些陪集首与伴随式的映射表。3.2 用枚举构造标准阵标准阵的构造代码如下思路是按错误重量 0,1,2,... 依次枚举所有图案先出现者即最小重量代表。import itertools import numpy as np def build_syndrome_table(H): 构造伴随式到陪集首的查表 H: (r, n) 校验矩阵, r n-k 返回: dict, key 为伴随式 tuple, value 为陪集首错误图案 n H.shape[1] r H.shape[0] total 2 ** r table {} for weight in range(0, n 1): for idx in itertools.combinations(range(n), weight): e np.zeros(n, dtypeint) e[list(idx)] 1 s tuple(np.mod(H e, 2)) if s not in table: table[s] e if len(table) total: return table raise RuntimeError(伴随式空间未能填满检查 H 的列秩)逻辑说明外层循环按错误的二进制个数汉明重量从 0 到 n 遍历内层itertools.combinations产生该重量下所有位置组合。每个组合生成一个错误图案 e计算伴随式并写表。因为从小到大枚举后到的必然不是陪集首所以直接跳过已存在的 key。完整表必须达到 total2^r 项如果退出循环仍未填满说明 H 的列秩不足校验方程线性相关属于构造错误。以下表格给出常用短码的建表和内存参考(n,k)r伴随式数量 2^r查表内存需求约枚举到的最小终止重量(7,4)3856 B1(8,4)416128 B2(15,11)416240 B2(31,26)532992 B2对 (7,4) 码重量 0 的图案 1 个重量 1 的图案 7 个一共正好填满 8 个伴随式终止重量为 1。对 (8,4) 和更长的汉明码需要枚举到重量 2 才能凑齐全部陪集首。3.3 译码主循环与查表实现查表译码函数只有三步算伴随式、查表得到 e、异或恢复码字。def hard_decode(r_b, H, table): 硬判决伴随式译码 r_b: (n,) 0/1 向量 返回: (c_hat, e) s tuple(np.mod(H r_b, 2)) if s not in table: raise ValueError(伴随式未在表中找到请检查 H 与 table 是否来自同一码) e table[s] c_hat np.mod(r_b - e, 2) # 等价于 bitwise_xor return c_hat, enp.mod(r_b - e, 2)与np.bitwise_xor(r_b, e)完全等价前者写起来不容易把 dtype 搞乱。table 的 key 用 tuple是因为 numpy 数组无法直接作为 dict 的 key这是实现时最容易踩的一处。当 s 为全零时e 是全零向量c_hat 等于 r_b说明接收向量本身就是合法码字。使用示例P np.array([[1, 1, 0], [1, 1, 1], [1, 0, 1], [0, 1, 1]], dtypeint) H compute_H(P) table build_syndrome_table(H) r_b np.array([1, 0, 0, 1, 0, 1, 1]) # 在第3位人为翻转 c_hat, e hard_decode(r_b, H, table) print(c_hat) # [1 0 1 1 0 1 1]错误被纠正这个例子验证了单比特翻转可以被完整还原。如果要测试 2 比特翻转比如把第 2、第 5 位同时翻转你会看到 c_hat 已经不再等于原始码字这正好引向下一节的边界讨论。3.4 超过 t 个错误译码器何时给出错误答案当真实错误图案重量超过 t标准阵查表会把它误判为某个更轻的陪集首两者差一个非零码字输出的是另一个合法码字而非原始码字。此时错误比特数没有规律短码可能错 1 位也可能错多位而且帧头标志也无法发现——因为 c_hat 满足 H c_hat^T 0校验通过。这个现象说明了硬判决查表译码的本质它把纠错换成找最近码字码字集合的 dmin 直接决定半径。dmin3 时能保证半径 1 的球不交叠噪声把接收点推到球外时结果没有任何保证。低信噪比下多位错误概率升高查表反而放大错误BER 曲线会出现负增益区这是下一章仿真要重点观察的现象。4. 线性分组码 BER 性能仿真BPSK 信道下的蒙特卡罗闭环4.1 仿真链路的五步结构仿真链路按顺序做六件事随机信息生成、系统编码、BPSK 映射、加高斯噪声、硬判决解调、查表译码、统计错误。由于每个 SNR 点都要跑成千上万帧编码器和查表译码器都应当是纯 numpy 向量化实现避免逐比特 for 循环。4.2 蒙特卡罗仿真主循环下面的simulate_ber函数把整条链路封装在一个 while 循环里用最小差错帧数控制统计精度。def simulate_ber(P, table, ebno_db_list, min_frame_errors50, max_frames100000, seed42): 线性分组码 BPSKAWGN 硬判决查表译码 BER 仿真 P: (k, r) 校验分量 table: build_syndrome_table 的返回 min_frame_errors: 每个 SNR 点至少统计的差错帧数 rng np.random.default_rng(seed) k P.shape[0] n k P.shape[1] H compute_H(P) ber [] for ebno_db in ebno_db_list: ebno_lin 10 ** (ebno_db / 10.0) sigma 1.0 / np.sqrt(2.0 * (k / n) * ebno_lin) total_bits 0 total_err 0 frame_err 0 frames 0 while frame_err min_frame_errors and frames max_frames: u rng.integers(0, 2, sizek) c systematic_encode(u, P) x 1 - 2 * c # 0 - 1, 1 - -1 y x rng.normal(0.0, sigma, sizen) r_b y 0 # 硬判决 c_hat, _ hard_decode(r_b, H, table) err int(np.sum(c_hat ! c)) frames 1 total_bits k total_err err frame_err 1 if err 0 else 0 ber.append(total_err / total_bits) return np.array(ebno_db_list), np.array(ber)min_frame_errors和max_frames是统计可信度的两个旋钮前者保证每个点都不会因为恰好只出现几个错误而抖动后者防止高信噪比下循环跑不完。total_err只统计前 k 位信息位的错误个数这样 BER 的定义与理论公式口径一致。运行示例P np.array([[1, 1, 0], [1, 1, 1], [1, 0, 1], [0, 1, 1]], dtypeint) H compute_H(P) table build_syndrome_table(H) snr_pts np.arange(0.0, 8.0, 0.5) ebno_db_list, ber simulate_ber(P, table, snr_pts, min_frame_errors100) for db, b in zip(ebno_db_list, ber): print(f{db:.1f} dB BER{b:.2e})跑完的曲线如果是平滑的 S 形编码链路基本可信若曲线毛刺多先把min_frame_errors提高到 200 再试。4.3 Eb/N0 与 σ 换算漏乘 R 会让曲线凭空左移噪声方差的换算是最容易出错的一步。发送符号 x∈{1,-1}每个码符号能量 Es1。BPSK 每符号承载 1 个码比特但信息比特要靠码率摊算Es/N0 (k/n) · (Eb/N0)σ² N0/2于是噪声标准差为σ 1 / sqrt(2 · R · EbN0_linear)其中 Rk/nEbN0_linear 10^(EbN0_dB / 10)。漏乘 R 的直接后果是曲线整体左移约 10·log10(R) dB对 (7,4) 码是 -2.4 dB对 (15,11) 是 -1.3 dB。这个偏移不改变曲线形状但会让结果和未编码理论基准对不上。Eb/N0 (dB)(7,4) 码 R4/7 的 σ未编码 R1 的 σ3.00.3690.3144.00.3150.2675.00.2670.2257.00.1770.149工程上换码率时最容易忘记同步更新 R。如果你把 (7,4) 换成 (8,4)码率从 0.571 变成 0.5同一 Eb/N0 下 σ 要乘 sqrt(0.571/0.5)≈1.07看起来很小但在高信噪比区足够让曲线偏移出误差范围。4.4 曲线解读与负增益区仿真结果的典型三段式曲线可以用作正确性判据。Eb/N0 低于约 1dB 时编码线在未编码 BPSK 理论线之上这是负增益区噪声翻转多位概率大查表译码反而放大错误。中间 1.5dB 附近是交叉点。4dB 以后进入瀑布区BER 快速跌落下降斜率和 dmin 正相关。交叉点位置主要由码率决定。低码率码冗余多同样 Eb/N0 下每个码符号分到的能量更少负增益区更长但瀑布区的斜率更陡。做方案选型时如果信噪比预算低于交叉点用分组码不但没好处还会让时延变大。提示若曲线在低 SNR 处与未编码线完全重合或明显高于理论线优先检查 σ 换算和硬判决阈值 0若高 SNR 处出现平台优先怀疑 table 构造时 H 的列秩不足。4.5 软判决要不要做查表硬判决是先判决再译码等价于在 BPSK 星座上按符号距离最近判决。短码场景下可以做一个简单改进对 n 较小的码直接枚举 2^k 个码字计算接收向量 y 与每个码字映射符号向量的欧氏距离选最近者作为输出。def ml_decode(y, codebook): 最大似然软判决译码仅适合短码 codebook: (2^k, n) 的码字表0/1 y: 长度为 n 的软信息向量 返回 (索引, 估计码字) x_all 1 - 2 * codebook # 二进制码字转 ±1 符号 dist np.sum((y[None, :] - x_all) ** 2, axis1) idx np.argmin(dist) return idx, codebook[idx]这个函数能拿到约 2dB 左右的性能提升但枚举复杂度是 O(2^k·n)(7,4) 码只有 16 个码字秒出结果(15,11) 码 2048 个码字勉强能跑(20,10) 码 1024 个码字还好但 (24,12) 码是 4096 个——再往上就明显吃力。实际工程里短分组码软判决更多用分阶统计译码OSD这类近似方法那个属于另一套复杂度模型这里不展开。5. 线性分组码仿真结果校准三个验证技巧5.1 用未编码 BPSK 理论线做锚点未编码 BPSK 在 AWGN 下的误码率是分组码仿真的天然校准器。把仿真的收发链路中原封不动去掉编码和译码只保留 BPSK 调制解调跑出来的 BER 应当贴合 0.5·erfc(sqrt(EbN0_linear))。可以这样验证from math import erfc def uncoded_reference(ebno_db): ebno_lin 10 ** (ebno_db / 10.0) return 0.5 * erfc(np.sqrt(ebno_lin))如果未编码仿真与这条理论线的偏差超过 0.3dB优先检查噪声方差和硬判决方向未编码链路准了再套上编码器查表问题就只可能出在编码或译码环节。5.2 高 SNR 点的错误统计陷阱曲线末尾抖动不是随机性是错误帧太少。经验做法是让每个 SNR 点至少统计 50 个差错帧再退出循环配合max_frames做上限保护。要看到 10^-5 量级的 BER至少需要 10^6 比特量级的样本同理如果某个点只统计到 10 次错误就给出了 10^-4 的估计这个点直接不要画出曲线。固定随机种子同样重要。np.random.default_rng(seed)在调试时可以复现任何一条问题曲线正式跑多组对比时再换随机种子验证稳定性。固定种子不会掩盖 bug但能大幅降低排查时的变量维度。5.3 标准阵表溢出的判断线标准阵项数等于 2^r(15,11) 是 16 项(23,12) Golay 码是 2048 项还能接受r 超过 16 时2^16 × n 的查表构建时间已经明显膨胀再往上就没有工程意义了。粗估内存时按表项数 × n × 8 字节算超过几百 MB 就果断换译码思路轻量级可以用 BCH 类循环码的代数译码或者转到置信传播类软译码。对只需要快速出短码 BER 曲线的场景(31,26) 汉明码的 32 项表和 (8,4) 扩展汉明的 16 项表是最有性价比的两组基线参数既能体现分组码的纠错增益又不会在查表和枚举上浪费调试时间。本文还有配套的精品资源点击获取