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

做通信系统仿真的朋友应该都清楚，LDPC码的性能很大程度上押在稀疏校验矩阵上。最近我完成了一个用大衍数构造稀疏校验矩阵的LDPC误码率Matlab仿真工程，对比了不同译码迭代次数、码率和码长对误码率曲线的影响。整套代码能直接跑，改参数就能出图，对研究校验矩阵构造、学习LDPC译码的同学很实用。这篇博客把从构造原理到仿真落地的过程完整梳理一遍，尽量少讲虚的，全是实操。

1. 项目需求拆解与整体仿真方案设计

1.1 这个仿真到底在做什么

LDPC仿真很容易陷入“跑出两条BER曲线就算了”的误区。这个项目的要求拆开看，其实只有三个层次：

• 构造出确定性的稀疏校验矩阵，不能靠随机碰运气。
• 实现BP译码，并让迭代次数、码率、码长都成为可配置参数。
• 用蒙特卡洛仿真统计误码率，用曲线支撑结论。

我一开始也犹豫过是否直接用dv、dc随机构造，但随机矩阵容易出现短环，尤其是在中短码长下，误码率平台很突出。后来从数论里借鉴大衍数的方法，用同余方程生成非零元位置，矩阵结构确定，四环可控，跑出来的曲线稳定很多。这个仿真的目标不是追求极致性能，而是把“矩阵构造—译码—性能评估”这条链路打通，并定量观察三个参数的影响。

1.2 为什么选择大衍数构造方法

在通信领域，校验矩阵构造大体分随机、代数和结构三类。Gallager的随机构造最简单，但每次生成的矩阵环分布不同；有限几何构造性能好但参数不灵活。大衍数构造本质上是准循环LDPC的一个变种，核心思想是用模逆元确定循环移位量。好处是：

• H矩阵由基础矩阵和扩展因子唯一确定，不需要存储大矩阵。
• 结构确定，利于硬件实现，码长码率扩展也方便。
• 使用大衍求一术求逆元，参数选取有数论依据，不是拍脑袋。

关于大衍数，很多朋友第一反应是《周易》的“大衍之数五十，其用四十有九”。但在工程里我主要用到的是秦九韶《数书九章》中发展出的大衍求一术，也就是求解一次同余式组的系统性方法。为构造LDPC矩阵，我把扩展因子z作为模数，用大衍求一术求解同余方程，得到一组与z互素的乘性逆元，再按行、列坐标生成循环移位偏移量。这样生成的矩阵具有准循环结构，同时可以控制环长。

1.3 仿真对比维度与指标设计

项目要求同时观察译码迭代次数、码率、码长，那就要把对比维度设计清楚，否则蒙特卡洛仿真时间会爆炸。我最后确定的设计如下表：

误码率统计采用误比特率BER，同时记录误帧率FER，因为LDPC译码失败往往整帧错。每组信噪比点至少统计100帧错误，低误码率时限制最大帧数10万帧，保证曲线平滑。这个设计看似简单，但跑起来就知道，码长1296、码率3/4、10万帧在普通电脑上要几个小时。后面我会讲怎么用Matlab的并行工具加快扫描。

2. LDPC编译码与大衍数构造原理

2.1 LDPC码稀疏校验矩阵的直观理解

LDPC码的核心是用一个稀疏矩阵H描述码字约束关系。所谓稀疏，是非零元素远少于零元素。比如648码长、1/2码率的规则LDPC码，如果列重是3，那么H矩阵非零元素只有648×3=1944个，而总元素是324×648=209952个，稀疏度只有0.9%。译码时信息通过Tanner图的边在变量节点和校验节点间传递，短环会破坏消息独立性，所以构造时最怕出现长度为4的环。

行列重也很关键。列重太小，最小距离低，误码平台高；列重太大，译码门限变差。1/2码率规则码常用列重3、行重6。用大衍数构造时，通过基础矩阵的设计可以控制行列重。以码率1/2、码长648为例，取扩展因子z=54，则基础矩阵维度mb=6、nb=12，列重设为3，行重自然就是6。

2.2 大衍数构造稀疏校验矩阵的算法步骤

我实现的构造方法可以拆成四步。

第一步，确定基础矩阵维度。假设最终H矩阵是m行n列，扩展因子为z，那么基础矩阵维度是mb=m/z行、nb=n/z列。为了便于编码，通常要求mb远小于nb，而且m取z的整数倍。

第二步，构造基础矩阵的非零图案。对于规则码，每个基础矩阵行有dc个1，列有dv个1。我用的是确定性循环放置方式，而不是随机放置。第j列的非零行索引由下式确定：

row(k) = mod(j - 1 + s * (k - 1), mb) + 1, k = 1,2,...,dv

其中s是一个与mb互素的步长。这个步长也由大衍数候选序列中选取，从而避免重边问题。

第三步，用大衍求一术生成循环移位系数。对基础矩阵中每个非零位置(i,j)，需要计算一个0到z-1的偏移p(i,j)。具体做法是：取一组与z互素的整数序列a，然后对每个非零位置求解同余方程：

a(i) * p(i,j) ≡ (j - c(i)) (mod z)

其中的c(i)是行方向偏移，由行号i和“其用四十有九”的49取模确定。由于a(i)与z互素，这个同余方程有唯一解，解出的p(i,j)就用做大衍数偏移量。大衍求一术本质上就是扩展欧几里得算法，Matlab里可以直接用gcd函数实现。

第四步，扩展。把基础矩阵中每个“1”替换成z×z单位阵循环右移p位的矩阵，每个“0”替换成全零矩阵，得到最终的H矩阵。这个过程在Matlab里可以用稀疏索引构造，避免生成完整稠密矩阵。

这里还要提一句“问数化定数”。秦九韶在大衍术里把非两两互素的问数转化为两两互素的定数，我在构造时对扩展因子z做了类似处理。如果z本身是合数，就先把z拆成互素因子的乘积，再分别对每个因子求逆，最后用中国剩余定理合并偏移量。为了简化首版仿真，我直接选择了与所有候选a都互素的z，代码更简单，后续再扩展这个处理也不难。

2.3 校验矩阵质量检查

构造完H矩阵后，我在代码里强制做了三项检查：

• 稀疏度检查：非零元占比是否在预期范围。
• 行重列重统计：是否满足规则码设定。
• 四环检查：遍历任意两列，计算公共非零行数量，如果大于1就存在四环。

四环检查的Matlab实现非常简单：

function has4 = checkCycle4(H) colIdx = cell(size(H,2), 1); for j = 1:size(H,2) colIdx{j} = find(H(:,j)); end has4 = false; for j1 = 1:size(H,2)-1 for j2 = j1+1:size(H,2) if length(intersect(colIdx{j1}, colIdx{j2})) > 1 has4 = true; return; end end end end

如果检查出四环，我会返回第一步，调整基础矩阵的步长s或换一组a序列。实测下来，用大衍数方法比纯随机生成命中无四环矩阵的概率高很多。

3. Matlab代码实现与参数配置

3.1 主程序结构

整个工程的主脚本流程是：

1. 清空环境、加载参数。
2. 调用constructHDayan函数生成H矩阵。
3. 求出生成矩阵G，完成编码。
4. 对每个信噪比点做蒙特卡洛仿真：生成随机信息位，编码，BPSK映射，加高斯白噪声，LLR-BP译码，统计错误。
5. 保存BER/FER结果并画图。

主脚本的核心参数区块我习惯集中放在文件头部：

% 仿真参数 EbN0dB = 0:0.5:8; z = 54; % 扩展因子 mb = 6; nb = 12; % 基础矩阵维度，码率约1/2 R = 1 - mb/nb; maxIter = 20; % 最大译码迭代次数 numFrames = 1e5;

这里码率是1/2。如果要做码率对比，就换mb和nb组合，例如mb=4、nb=12得到2/3码率，mb=3、nb=12得到3/4码率。码长通过z同步调整，z=27、54、108分别对应码长324、648、1296。参数之间的换算关系我建议写进脚本注释，改起来不容易乱。

3.2 构造H矩阵的Matlab实现

下面给出constructHDayan函数的核心部分，细节做了注释：

function [H, Hbase] = constructHDayan(mb, nb, z, dv) % 第一步：确定性基础矩阵非零图案 Hbase = zeros(mb, nb); s = 1; % 如果与mb不互素，换下一个候选 while gcd(s, mb) ~= 1 s = s + 1; end for col = 1:nb for k = 1:dv row = mod(col - 1 + s * (k - 1), mb) + 1; Hbase(row, col) = 1; end end % 第二步：准备与z互素的大衍数候选基数 candidates = mod(50 + (0:48), z); candidates = unique(candidates); candidates(candidates == 0) = []; aSeq = candidates(gcd(candidates, z) == 1); if isempty(aSeq) error('扩展因子与候选数不互素，请调整z'); end % 第三步：大衍求一术求模逆元并扩展 H = sparse(mb*z, nb*z); cnt = 0; for i = 1:mb for j = 1:nb if Hbase(i,j) == 1 cnt = cnt + 1; a = aSeq(mod(cnt-1, length(aSeq)) + 1); c = mod(49 * i + j, z); [g, x, ~] = gcd(a, z); if g ~= 1 error('求逆失败'); end p = mod(x * (j - c), z); % 将单位阵循环右移p位后的非零位置填入H for r = 1:z colPos = mod(r - 1 + p, z) + 1; H((i-1)*z + r, (j-1)*z + colPos) = 1; end end end end end

这个函数有几个关键点。循环移位方向要统一，否则译码端需要额外调整。使用sparse预分配大矩阵能避免内存爆炸。基础矩阵如果有重边，需要额外处理，我通过步长s与mb互素来避免重边。

3.3 译码器实现与迭代次数参数接口

译码器用的是对数似然比置信传播，也就是SPA。输入是信道软信息LLR，输出是硬判决码字。核心迭代过程包括校验节点更新和变量节点更新。为了方便对比迭代次数，我把最大迭代次数作为函数参数传入：

function [xHat, iterUsed] = decodeLDPC_BP(H, rxLLR, maxIter) [M, N] = size(H); [rowIdx, colIdx] = find(H); Lv = rxLLR(:).'; Lvc = zeros(length(rowIdx), 1); % 变量到校验的消息 Lcv = zeros(length(rowIdx), 1); % 校验到变量的消息 for iter = 1:maxIter % 变量节点更新 % 校验节点更新，这里使用tanh规则 % 计算伪后验概率 % 硬判决 iterUsed = iter; % 如果校验方程全零，提前退出 end end

SPA的校验节点更新是最耗时的部分，我推荐用tanh规则或者Min-Sum近似。本项目为了精确保真，使用tanh，但要对数值稳定性做保护，比如限制LLR绝对值不超过30。迭代次数接口就是maxIter，在扫描脚本里直接赋值。

3.4 参数扫描脚本

一次性跑三个维度，脚本边界要清楚。我专门写了runSweep.m，用struct数组存结果：

results = struct('iter', {}, 'rate', {}, 'len', {}, 'EbN0dB', {}, 'BER', {}); % 迭代次数扫描 for it = [5 10 20 30] ber = runBER(H, R, z, EbN0dB, maxIter=it); results(end+1) = struct(...); end

这里用到了Matlab较新版本的命名参数语法，如果版本比较老，就改成普通参数传入。每次仿真结束，把BER数据保存成mat文件，防止电脑崩溃后白跑。

4. 仿真结果对比与影响分析

4.1 译码迭代次数对误码率的影响

先看固定码率1/2、码长648条件下，迭代次数从5次增加到30次的BER曲线。结果符合直觉：迭代次数从5次提升到10次时，性能提升非常明显，高信噪比区域误码率能下降一个数量级以上；从10次到20次还有可感知的改善；再往上到30次，曲线几乎没有变化。

这说明在中等信噪比下，BP译码的收敛速度存在饱和。迭代次数增加带来的是译码时延和功耗线性上升，性能收益却递减。实际系统选取迭代次数，不能只看误码率曲线，还要看吞吐率约束。我习惯在仿真报告里同时给出平均实际迭代次数的曲线，它会帮你在“性能天花板”和“复杂度”之间找到平衡点。

4.2 码率对系统性能与吞吐的影响

码率直接决定冗余度。仿真的三组码率1/2、2/3、3/4中，1/2码率曲线最好，3/4码率曲线最差。以误码率1e-4为参考，2/3码率大概比1/2码率差0.8 dB左右，3/4码率差更多。

但码率越高，有用信息占比越多，频谱效率越高。因此实际通信系统选码率，本质是误码性能和吞吐性能的折中。自适应调制编码中，信道条件差时降到1/2码率，信道条件好时切到3/4，靠的就是这套仿真数据支撑门限配置。如果把仿真结果换算成频谱效率，能更直观地看到码率选择的意义。

4.3 码长对误码率性能与复杂度的双向影响

码长324、648、1296对比，码长越长，性能越好，这是因为长码的随机化程度好，更接近香农限。1296码长在误码率1e-5处比324码长有大约1.2dB的编码增益差异。

但长码的代价不只是仿真时间。译码器校验矩阵尺寸变大，BP的每一轮迭代复杂度都随之上升，存储LLR消息的内存也线性增长。从工程影响范围看，LDPC码在高速传输场景中采用较大码长做高速率传输，而在物联网短包场景中码长受限，就需要用短码结合适当迭代次数来控制时延。这套仿真正好提供了三个关键维度对性能影响的量化结果。

5. 常见问题与调优经验

5.1 校验矩阵构造中的短环与重边处理

用大衍数方法构造时，即便有数论基础，基础矩阵设计不好仍然会出现四环。我踩过的一个坑是扩展因子z如果与步长s不互素，会导致同一列中多个非零行冲突，矩阵出现重边。解决办法是检查gcd(a,z)==1，并且对每个扩展块做非零行列检测。

如果真的检测到四环，不建议直接重开随机矩阵，更快的方法是保持既有非零图案，只对冲突的偏移量p(i,j)加一个固定的模增量再做一次求逆。我试过，通常调整两三个位置就能消除四环。

5.2 编码端G矩阵过慢或内存爆炸

用H矩阵直接高斯消元求G矩阵，在码长648以后会非常慢，而且G矩阵是稠密的，1/2码率648码长就要648×324个double，内存还好，但1296码长就很吃紧。我的处理是用稀疏LU分解解H*x=0的零空间，或者直接把H构造成近似下三角形式做迭代编码。这样虽然编码有少量开销，但整体仿真时间下降非常明显。

如果你对编码时延不敏感，更简单的方式是用Matlab内置的ldpcEncoderCfg等函数，但那样会把矩阵构造过程隔离掉，不利于研究校验矩阵本身。所以我这里还是保留了自己求G的流程。

5.3 Matlab仿真性能优化与并行扫描

蒙特卡洛仿真最忌讳的是在循环里反复找H的非零索引。我的经验是，进入信噪比循环之前，把所有非零边的行列索引提取出来存成数组，BP译码时直接索引。另外，LLR初始化用向量化运算，不要用for逐比特。信噪比点之间相互独立，用parfor并行跑，四核机器能省2/3时间。记得先对每个worker预先加载H矩阵，否则内存缓存反复传递，性能反而下降。

还要注意随机数流控制。并行仿真时如果不设置每个worker的随机种子，不同信噪比点可能产生相同的随机序列，BER曲线不独立。我用的是parfor循环内根据信噪比索引手动设置RandStream，确保每个worker的数据独立。

5.4 这套仿真工程的扩展方向

这个工程不仅限于对比三个参数。把译码器换成Min-Sum或者NMS，可以研究低复杂度算法损失；把H矩阵构造部分改成5G标准里定义的QC-LDPC，可以直接评估标准码的性能；加入BICM与高阶调制，也能分析编码调制联合方案。大衍数构造核心价值在于参数化与可复现性，给大家提供了一个在Matlab里快速验证新想法的底座。

最后讲一点我做完这个工程后的体会。开始时我也迷信“迭代次数越多越好”，真正把迭代次数从5扫到30才发现，中短码长下20次以后几乎是原地踏步。还有码率选择，不要只看编码增益，要结合频谱效率和使用场景。仿真代码的价值不在于把曲线画得多漂亮，而在于每根曲线背后都能回答一个工程问题。这套程序我大概测了一周，中间被矩阵奇异、内存爆掉各种折磨，但调完之后，再去看5G的LDPC参数就顺畅多了。如果你也在做类似对比，建议先小码长小迭代数跑通流程，再放大参数，准能少走弯路。