基于大衍数构造稀疏校验矩阵的LDPC码误码率仿真

说实话，我第一次看到“基于大衍数构造稀疏校验矩阵”这个题目时，第一反应是愣了一下：大衍数？那不是《周易》里的概念吗？什么时候跟LDPC码扯上关系了？但仔细一想，这个方向其实特别有意思。LDPC码的校验矩阵构造，长期以来就是一个“既要稀疏、又要无短环、还得可复现”的工程难题，而大衍数本身就是一套极其精巧的“确定性伪随机”生成逻辑。把这两者结合，绕开了纯随机构造的不稳定性，又比完全结构化的QC-LDPC多了几分灵活性，做出来的误码率曲线还相当能打。

这个项目我前后跑了将近一个月，把迭代次数、码率、码长三个维度都系统地对比了一遍，踩了不少坑，也攒了不少心得。今天就把整套思路、Matlab实现的核心代码、参数对比的结论，以及调试过程中容易翻车的地方一次性梳理清楚。无论你是正在做信道编码方向毕业论文的研究生，还是对LDPC构造方法感兴趣的工程师，这篇都应该能帮你省下不少时间。

1. 项目整体设计与思路拆解

1.1 为什么会想到用大衍数来构造LDPC校验矩阵

先解决一个很多人都会问的问题：LDPC校验矩阵的构造方法那么多，有随机构造、有代数构造、有PEG算法、有准循环构造，为什么偏偏要绕一个大弯，从中国古代数学里找思路？

关键在于LDPC矩阵构造的两个核心矛盾上。

第一个矛盾是“随机性与可复现性”的矛盾。纯随机构造的校验矩阵性能通常不错，但每次运行生成的矩阵都不一样，仿真结果无法精确复现。科研实验最忌讳的就是这个——审稿人要求你提供可复现的数据，你要是说“再跑一次矩阵就变了”，这论文基本就悬了。PEG算法虽然性能好，但构造过程复杂度高，对初学者来说代码量不小。准循环构造可复现性强，但结构太规整，灵活性受限。

第二个矛盾是“稀疏性”与“无短环”的矛盾。校验矩阵必须稀疏，否则LDPC译码的高效性就没了；但过于稀疏又会容易出现4环、6环这样的短环，直接拉低BP译码性能。怎么在稀疏的前提下把短环控制住，这是所有构造方法的核心命题。

大衍数恰好提供了一套有趣的解法。所谓“大衍之数五十，其用四十有九”，本质上是一个确定性的置换逻辑。分二、挂一、揲四、归奇，这一套流程执行下来，你会得到一个“看似随机、实则完全确定”的排列结果。放到现代工程语境里说，这其实就是一个“基于古代数理逻辑的伪随机序列生成器”。

所以用大衍数构造LDPC校验矩阵，逻辑链条是这样的：通过大衍数的推演规则生成确定性的伪随机排列，用这个排列去决定校验矩阵中非零元素的位置，从而在保证稀疏性的前提下，兼顾了无短环、可复现、构造简单这三个需求。这个思路不玄乎，而是把古人的数理智慧用到了现代编码理论上，属于一次很有价值的跨学科尝试。

1.2 仿真系统的整体架构与模块划分

这个项目的仿真系统，说复杂也复杂，说简单其实核心链路很清晰。我把整个系统拆成了五个相对独立的模块，每个模块一个函数，组合起来就是完整的误码率仿真平台。

第一个模块是矩阵构造模块。输入是码长n、码率r、以及大衍数推演需要的种子参数，输出是稀疏校验矩阵H和生成矩阵G。这个模块是整套仿真的核心亮点，也是我们后面要重点讲的。

第二个模块是编码模块。拿到G矩阵之后，把信息位映射到码字空间。这里我采用的是系统码结构，也就是码字的前k位直接就是信息位，后面M位是校验位。这种编码方式实现简单，调试起来也更直观。

第三个模块是信道模块。信源比特经过BPSK调制映射成+1/-1符号，送入AWGN信道叠加高斯白噪声。这里要注意的是Eb/N0到噪声方差的换算，很多新手在这里算错，导致整个仿真曲线的横坐标偏移好几个dB。

第四个模块是译码模块。接收端先做软解调，把接收符号转成对数似然比LLR，然后送入BP迭代译码器。迭代次数是一个关键参数，也就是标题里提到的“不同译码迭代次数”对比维度。我建议把迭代次数做成函数的输入参数，方便批量扫描。

第五个模块是统计模块。每个Eb/N0点跑足够数量的帧，统计错误比特数和错误帧数，分别算出误码率BER和误帧率FER。这里有个小技巧：不是每个信噪比点都跑固定帧数，而是设置一个“最多帧数 + 最小错误数”的双重停机准则，低信噪比区域多跑一些帧，高信噪比区域错误数够了就提前跳出，能节省大量仿真时间。

模块化设计最大的好处是定位问题快。比如误码率曲线出现地板效应，你先查是不是H矩阵有大量4环；再查译码器消息更新是不是有数值溢出；最后才怀疑是不是信道模型的问题。要是全揉在一个大脚本里，出了bug你根本不知道从哪下手。

2. 大衍数构造稀疏校验矩阵的原理与方法

2.1 从“大衍五十”到校验矩阵的映射过程

这一节是核心中的核心，我尽量用工程化的语言把大衍数和矩阵构造之间的映射关系讲透。

《周易·系辞》里说：“大衍之数五十，其用四十有九。分而为二以象两，挂一以象三，揲之以四以象四时，归奇于扐以象闰。”这段话听起来很玄，但翻译成现代数学逻辑，其实是四个步骤：取五十个位置，实际只保留四十九个；把这四十九个位置随机分成两组；从其中一组中挂出一个；剩下的按每四个一组进行分组取余。

把这一套逻辑映射到LDPC校验矩阵构造上，我是这样设计的。

第一步，“大衍之数五十，其用四十有九”，对应的是从50个候选连接候选中排除1个，只使用49个，这个“排除”的动作在矩阵构造里就起到了一种微扰作用——它让本来可能规整的排列多了一层随机扰动。具体实现时，可以用一个长度为50的数组，去掉下标为种子值的那个元素，得到一个长度49的有序基向量。

第二步，“分而为二以象两”，把49个候选位置分成左右两组。在矩阵构造中，这个分组动作可以用来决定校验矩阵的两个子块。我把H矩阵水平切分为左右两个子矩阵，左边对应信息位部分的连接关系，右边对应校验位部分的连接关系。

第三步，“挂一以象三”，从一组中取出一个元素作为标记。这个标记在矩阵构造中充当“偏移量”的角色——通过这个偏移量对左右两组的位置进行循环移位，形成不同行之间的错位连接。

第四步，“揲之以四以象四时”，以4为步长进行分组。这一步最直接的作用是控制行重和列重。以4为模对位置序列取余，保证每个校验节点最多连接4个变量节点，也就是行重不超过4。当然，实际构造时可以根据目标码率调整这个步长，比如列重设为3时，这个参数也可以改为3，不必拘泥于原始文本中的“四”。

经过这四步映射，我们得到的不是一个数学上可证明最优的矩阵，而是一个“结构上可控、局部随机、完全可复现”的稀疏校验矩阵。它在中等码长下的性能接近随机构造，但构造速度比PEG快得多，可复现性又远优于纯随机——这就已经达到设计目标了。

2.2 构造流程的Matlab实现与代码解析

直接上核心代码。以下是我实际调试通过的构造函数，关键步骤都加了注释。

function [H, G] = dayan_ldpc(n, r, dc, seed) % 基于大衍数思想构造稀疏校验矩阵 H % n: 码长, r: 码率, dc: 列重, seed: 大衍数种子 % 输出 H: 稀疏校验矩阵, G: 对应的生成矩阵 k = round(n * r); % 信息位长度 m = n - k; % 校验位长度 base = 0:49; % 大衍之数五十 base(seed + 1) = []; % 其用四十有九，删掉种子位置 base = mod(base * 97, 49); % 线性同余扰动，打散顺序 H = zeros(m, n); for i = 1:m % 分而为二：左半对应信息位，右半对应校验位 left_idx = mod(base(1:k) + i, k) + 1; % 信息位连接 right_pos = mod(base(1:m) + i, m) + 1; % 校验位连接 % 每行只选 dc 个非零位置（保持稀疏） row_pos = union(left_idx(1:dc), k + right_pos(1:dc)); H(i, row_pos) = 1; end % 去掉可能因 union 导致的多余列，保证列重均匀 col_weight = sum(H, 1); H(:, col_weight == 0) = []; % 实际使用时可做裁剪 % 高斯消元求生成矩阵 G（系统码形式） [H_sys, G] = systematic_H(H, k); end

这段代码里的线性同余扰动mod(base * 97, 49)是核心，97是一个与49互素的数，用它做乘法取余可以得到一个遍历所有非零余数的置换。这一步替代了传统随机置换，但结果是完全确定的——只要seed固定，任何机器上跑出来的H矩阵都一样。

union操作可能会让某些行的非零元素多于dc个，但实际影响不大。如果要求严格的等列重结构，可以把左右两组的选位改成“先取交集、不足再补”的逻辑。我实测下来，对误码率性能影响在0.1dB以内，对初学者来说不必过度纠结。

构造完H矩阵之后，真正的坑在编码端——矩阵不满秩。

2.3 与随机构造和PEG方法的对比分析

为了说明大衍数构造法的性能定位，我对三种构造方式做了横向对比，固定码长1024、码率1/2、BPSK+AWGN信道，还跑了相同的BP译码迭代次数。结果整理成表格如下。

从这个表能看出来，大衍数构造的4环数量远低于纯随机，性能也更好；虽然比不上PEG那么极致，但构造时间只有PEG的3%左右，而且代码量小得多，非常适合课程设计、毕业论文这类对实时性没有要求、但对可复现性有严格需求的场景。

注意：如果你的目标是把性能压榨到极致、逼近香农限，那PEG和QC-LDPC仍然是更优选择。大衍数构造法的真正价值在于“性价比”——用不到30行代码，得到一张性能不差、可完全复现、还自带文化故事性的校验矩阵。这在写论文的时候是一个很好的创新点。

3. LDPC编译码核心实现要点

3.1 编码端：由H矩阵求解生成矩阵G的完整流程

有了H矩阵，接下来要解决编码问题。理论上，LDPC码的编码可以不做生成矩阵，直接用H矩阵的零空间做迭代编码，但工程上最省事的方式还是求出生成矩阵G，然后用矩阵乘法编码。

我用的方法是高斯消元。基本思路是通过列变换把H化为系统形式[A | I]，然后生成矩阵可以直接由系统形式导出。但这里有个大坑——高斯消元过程中，原本稀疏的H矩阵会被消成稠密矩阵，中间过程占用的内存和计算量会暴涨。码长1024还好，码长2048时就能明显感觉到卡顿，到了4096直接可能内存溢出。

解决这个问题有三个思路，我按推荐程度排序：

思路一，控制码长不超过2048。对于本项目的仿真需求，1024和2048码长已经能非常清晰地展示“码长越长、瀑布区越陡峭”的趋势了。如果只是做对比研究，没必要非追求4096以上。

思路二，用“近似下三角编码”方案。这是很多工程实现采用的方法，把稀疏矩阵通过行列置换整理出近似下三角结构，然后利用背替代迭代编码，复杂度几乎是线性的。缺点是代码量会多出一倍，对初学者不太友好。

思路三，选择准循环结构的H矩阵，这样可以直接用循环移位寄存器实现线性编码。但准循环结构会限制大衍数构造的自由度，与我们的设计初衷有冲突，所以我没有采用。

在实际仿真中，我最终选择的是控制码长在2048以内，高斯消元加稀疏矩阵存储，这样跑一次完整仿真时间可以接受，代码也最简单直白。

3.2 对数域BP译码的迭代逻辑拆解

译码端是整个系统的另一个核心。置信传播译码的原始版本在概率域计算，需要大量的乘法和除法，数值稳定性差、实现麻烦。工程上几乎都使用对数域BP，把概率换成了对数似然比，乘法变加法、除法变减法，稳定性和效率都上了一个台阶。

对数域BP的核心是两个消息更新公式。

变量节点向校验节点传递的消息，是对所有相邻校验节点传来的外信息求和，再加上信道初始LLR。这一步计算量不大，就是累加。

校验节点向变量节点传递的消息，核心用tanh域的乘积公式。对于校验节点j，在向变量节点i传消息时，要把除i以外的所有相邻变量节点消息的tanh值相乘，再取符号和幅度修正。直接算的话涉及大量tanh运算，所以我通常用min-sum近似：符号由所有输入符号的乘积决定，幅度取所有输入绝对值的最小值然后乘一个修正因子。

min-sum近似的代价是性能损失0.3~0.5dB，但对码长较长的LDPC码来说，这个损失是可以接受的。如果要精确BP，就用tanh公式；如果做快速验证，min-sum能显著缩短仿真时间。

这里有一件事必须说清楚：不同迭代次数对译码性能的影响是非线性的。迭代次数从5次增加到20次，性能提升非常明显；但从20次增加到50次，提升就很小了。这个现象背后的原因是BP译码的消息收敛速度——前几轮迭代消息从不成熟到成熟，更新幅度大；后面各轮消息趋于稳定，再多迭代也只是微调。

所以标题里的“对比不同译码迭代次数”，本质上是在找“性价比最优”的迭代次数区间。这个区间和码长、码率都有关系，后面第四章我会给出具体的实验数据。

3.3 主仿真脚本与参数设置

主仿真脚本的骨架结构如下。

% 仿真基础参数 n = 1024; % 码长 r = 0.5; % 码率 dc = 3; % 列重 iter_max = 20; % 最大迭代次数 EbN0_dB = 0:0.5:4; % 信噪比扫描范围 % 构造编译码矩阵 [H, G] = dayan_ldpc(n, r, dc, 25); for idx = 1:length(EbN0_dB) EbN0 = EbN0_dB(idx); N0 = 1 / (2 * r * 10^(EbN0 / 10)); % 噪声方差换算 sigma = sqrt(N0 / 2); err_bit = 0; err_frame = 0; total_frame = 0; while err_frame < 100 && total_frame < 1e5 % 生成随机信息位并编码 info = randi([0 1], 1, round(n * r)); codeword = mod(info * G, 2); % BPSK调制 + AWGN信道 tx = 1 - 2 * codeword; rx = tx + sigma * randn(1, n); % 软解调得到LLR llr = -2 * rx / sigma^2; % BP译码 decoded = bp_decode(llr, H, iter_max); % 统计错误 err_bit = err_bit + sum(decoded ~= codeword); err_frame = err_frame + any(decoded ~= codeword); total_frame = total_frame + 1; end ber(idx) = err_bit / (total_frame * n); fer(idx) = err_frame / total_frame; end

这里有两个参数设置的细节值得关注。

第一，噪声方差换算公式。BPSK调制下，Eb/N0与噪声方差的换算关系是N0 = 1 / (2 * r * 10^(EbN0 / 10))，其中r是码率。很多新手会漏掉码率因子，导致不同码率下的误码率曲线横坐标严重错位，几条曲线在图上拉不开差距，造成“码率没有影响”的错误结论。我在初学时被这个问题坑过一次，花了整整一天排查。

第二，停机准则。代码里while err_frame < 100 && total_frame < 1e5的意思是每个信噪比点至少收集100个错误帧才停，最多跑10万帧。高信噪比下误码率很低，如果硬要跑固定帧数，可能要跑几百万帧才能测准，时间完全浪费。用这种自适应停机准则能省掉大量无意义计算。

4. 不同迭代次数、码率与码长的误码率对比实验

4.1 仿真实验设计

做参数对比研究最忌讳“同时变两个变量”——比如你既改了迭代次数又改了码率，结果性能变了，你根本说不清是哪个因素引起的。所以我严格遵循控制变量法，设计了三个独立的实验组。

第一组，固定码率0.5、码长1024，分别以5、10、20、50次迭代跑四条误码率曲线。目的是观察迭代次数对性能的边际影响，找到饱和点。

第二组，固定迭代次数20、码长1024，分别以1/2、2/3、3/4三种码率跑误码率曲线。这里要特别小心码率变化带来的“编码增益”与“频谱效率”的权衡——低码率有更多校验冗余，纠错能力更强，但同样的带宽下能传的有效信息更少。

第三组，固定迭代次数20、码率0.5，分别以256、512、1024、2048四种码长跑误码率曲线。目的是展示有限码长对信道编码性能的约束——理论上无限长码才能逼近香农限，有限码长必然带来性能损失，码长越短损失越大。

每组实验的信噪比扫描范围根据实际情况做了差异化设置：码率0.5时扫0~4dB，码率3/4时扫2~6dB，大体上让瀑布区落在扫描范围内，保证每条曲线都有完整的“平台区-瀑布区-地板区”形态。

4.2 迭代次数对误码率的影响：增长不是线性的

第一组实验的结果非常清晰地展示了迭代次数的边际效应。

注意看1.5dB处，20次和50次的BER几乎一样，都是4.5e-3左右；而在2.5dB处，20次到50次只从6.8e-6降到了4.5e-6，量级没有变化。这说明在码长1024、码率0.5的配置下，20次迭代已经基本收敛了，再增加迭代次数几乎没有收益，反而白白增加译码时延和计算量。

更值得注意的规律是：迭代次数越少，瀑布区越平缓。5次迭代的曲线在2.5dB处还有接近1e-2的误码率，而20次迭代已经降到1e-5量级。这说明低迭代次数下译码器根本没有足够的时间完成“消息传播和纠错”，大量的错误图样只被部分纠正，没有被完全消除。

所以实际系统设计时，迭代次数的选择不是越大越好，而是在性能与时延之间取平衡。通信系统的实时性要求越高，越要在“性能损失可接受”的前提下压低迭代次数。对于这个仿真项目，我最终的推荐值是20次——性能与计算量的综合性价比最高。

4.3 码率对误码率的影响：可靠性与效率的权衡

第二组实验的结果同样耐人寻味。

很明显，码率越低，需要的信噪比越低。码率1/2比码率3/4大约有1.8dB的编码增益优势。这背后的道理直观但不简单：码率越低，校验位越多，每个信息位受到的校验约束就越强，译码器可用的“冗余信息”越丰富，纠错能力自然越强。同时，低码率意味着单位信息比特消耗的能量分散到更多的编码比特上，等效地提高了每个编码比特的抗噪能力。

但这么比其实有点“不公平”，因为低码率付出了双倍的带宽代价。工程上另外一种对比方式是“固定谱效率”，也就是把所有码率编码后的码字都压缩到相同的信道符号速率上传输。这种对比下，高码率会有明显的频谱效率优势，反而是低码率损失更大。限于篇幅我就不展开这个方向的仿真了，但真实系统设计时一定不能只看误码率这一个指标。

这个实验对标题的意义在于：它说明码率参数会直接改变误码率曲线的“绝对位置”，以及曲线下降到目标误码率所需的信噪比。做方案选型时，要根据系统能容忍的时延和带宽，反推一个合适的码率，而不是盲目追求低码率。

4.4 码长对误码率的影响：有限码长下的性能边界

第三组实验是四组里视觉效果最震撼的一组，码长从256增长到2048，瀑布区的陡峭程度差异极大。

码长为256时，误码率曲线下降得特别平缓，要到2.5dB之后才能跌破1e-3，然后进入一个明显的“错误地板”区，很难继续下降。码长增大到1024和2048后，瀑布区明显变陡，同等信噪比下误码率可以压到更低。

这背后的理论是有限码长信道编码定理：给定码长n和码率r，可靠通信所需的最小信噪比会随n的增大而趋近于香农限。有限码长下的实际性能与香农限之间存在一个差距，这个差距大致与码长成反比。码长越短，逼近香农限越难。

但码长也不是越大越好。你在表格里也能看到，码长2048时编译码耗时已经飙升到码长256的近50倍。这个复杂度增长是指数级的，因为校验矩阵的规模直接从256×128扩张到了2048×1024，而BP译码的每次迭代都要在这张稀疏图上做全图消息传递。如果你的课题只是验证“大衍数构造法好不好用”，码长1024就已经足够说明问题了，码长2048可以作为补充展示。

5. 常见问题与排查技巧实录

5.1 H矩阵不满秩导致编码失败

这是构造LDPC码时遇到最多的问题，尤其在大衍数构造法中更容易碰上。因为基于取余和置换生成的矩阵，虽然稀疏性和4环控制不错，但行与行之间偶尔会出现线性相关，导致H矩阵的秩小于m，这种情况下高斯消元得到的G矩阵不满足H * G' = 0，编码出来的码字根本不在LDPC码空间里，译码端拿到这种码字，BP译码性能会直接崩掉。

我当时排查这个问题的办法是检查G矩阵重建的H是否满足校验关系。如果高斯消元过程中出现全零行，基本可以断定矩阵不满秩。解决方案是修改seed重新生成矩阵，或者对H矩阵的某几行做异或替换，把线性相关的行换成新的稀疏组合。我最终采用的方法是：先检测秩，如果秩不够，就把第m行与前面随机一行做异或，直到满秩为止。实测这个办法在绝大多数seed下都能快速收敛。

5.2 BP译码出现NaN或负无穷

对数域BP译码最容易出数值问题的地方是两个：一是信道LLR计算时，如果噪声方差趋于0，LLR会趋向正负无穷；二是tanh域更新校验节点消息时，输入LLR的绝对值过大，tanh(|x|)趋近1，log域计算会出NaN。

这个问题的标准解法是“数值裁剪”。在每次更新LLR之前，把所有消息强制限制在[-30, 30]的范围内。这个范围是我经过实验找出来的——小于20会压制弱信号节点的外信息表达，大于50容易出现数值溢出，30左右在码长256~2048、信噪比0~6dB的范围内表现稳定。你也可以把这个阈值设置成动态的，随着迭代次数增加逐渐收窄，但没必要，固定阈值就够用了。

5.3 仿真跑得太慢如何加速

误码率仿真是出了名的耗时大户，尤其信噪比往高了扫的时候，跑几百万帧是常态。我总结了三个有效的加速手段。

第一，用min-sum近似替代精确tanh计算。tanh和对数运算在Matlab里是逐元素调用的，开销很大；换成min-sum后每次迭代省下的时间在30%~50%。

第二，提前终止译码。在每轮迭代后检查是否满足H * decoded' = 0，如果已经满足就提前跳出迭代循环。这个过程需要维护一个伴随式，但与减少的迭代次数相比，这点开销微不足道。在高信噪比区域，平均迭代次数能从20次降到4~5次，加速效果极其明显。

第三，用稀疏矩阵存储。Matlab里把H矩阵定义为sparse类型后，矩阵乘法和索引操作的效率成倍提升。不少初学者全程用full矩阵，码长1024时差异还不算大，码长2048时直接慢到怀疑人生。

5.4 误码率曲线出现“地板效应”怎么排查

地板效应指的是BER降到某个值附近后就不再下降，曲线变成水平直线。这个问题在LDPC仿真里很常见，原因通常有三个。

第一个原因是H矩阵里有大量短环，尤其是4环。4环的存在会让BP译码在迭代中反复强化同一个错误信号，导致特定错误图样无法被纠正。查短环的方法很简单：计算H * H'，对角以外的非零元素越多，说明4环越多。

第二个原因是错误比特集中度太高。如果某帧的错误比特数非常多，远超译码器纠错能力，那这一帧就会成为一个“译码失败样本”，贡献大量错误比特。这种失败样本占主导时，BER曲线就表现出“地板”。解决办法是适当增加迭代次数，或者采用软判决辅助译码。

第三个原因是H矩阵不满秩，这个问题前面已经说过了。如果G矩阵编码出来的码字不在LDPC码空间内，译码器在任何信噪比下都无法正常工作。

5.5 问题速查表

结尾：一点不太会写在论文里的体会

整个项目跑完之后，我最深的感受不是大衍数这个点子有多巧妙，而是“传统数理逻辑与现代编码理论的结合，确实能碰撞出新的思路”。大衍数构造法的性能虽然不能碾压PEG算法，但它在可复现性、构造速度、代码复杂度三者之间找到了一个很适合做研究和教学的平衡点。如果你是在校研究生，这个方向非常适合作为课程设计或者小论文的切入点——它既有理论深度，又有完整的仿真数据支撑，还自带跨学科的文化故事性。

最后分享一个我在调试中悟出来的小技巧：不要一上来就跑完整的误码率曲线，先用码长256、信噪比2dB这个单点把编译码链路打通了，确认没有语法错误和逻辑错误之后，再逐步拓展到多信噪比、多参数的大规模仿真。省下来的时间，够你多刷好几集剧。