大衍数构造稀疏校验矩阵的LDPC码误码率仿真与参数对比

做信道编码方向的研究，绕不开LDPC码。我最近把一个基于大衍数构造稀疏校验矩阵的LDPC误码率仿真项目完整跑了一遍，对比了不同译码迭代次数、码率、码长对性能的影响，全程在Matlab里实现。这篇文章把整个项目的设计思路、核心原理、代码实现和避坑经验都过一遍，给正在做LDPC相关仿真或者准备开题的朋友做个参考。

项目名字看着长，拆开其实就三件事：第一，用大衍数这套数理序列去生成稀疏校验矩阵；第二，在Matlab里搭一条完整的误码率仿真链路，包括编码、调制、AWGN信道、置信传播译码和误码统计；第三，把译码迭代次数、码率、码长这三个变量分别拎出来做对比实验。整个过程跑完，你会对LDPC码的性能特征、矩阵构造对误码率的影响、仿真参数怎么设，都有一个非常直观的认识。

1. 项目全貌与整体设计思路

先说清楚这个项目到底在解决什么问题。LDPC码进入5G、卫星通信、光传输这些实际系统已经有相当年头了，它最大的特点就是校验矩阵极度稀疏，配合迭代译码可以逼近信道容量极限。但构造校验矩阵这件事，从学术到工程有很多种做法，其中"基于代数结构的确定性构造"是一大类，大衍数构造就属于这个方向。

1.1 一句话拆解项目目标

这个仿真项目的完整链条是：先生成稀疏校验矩阵H，然后把它用于编码、加噪、译码，最后统计误码率。对比维度有三个——译码迭代次数、码率、码长。

所谓"基于大衍数构造"，指的是用一组从大衍数衍生出来的序列，去填充准循环LDPC码的循环移位矩阵。这样可以得到一个确定性的、可复现的校验矩阵，而不是靠随机函数碰运气。这个矩阵满足稀疏性、无短环等基本要求，性能稳定，便于后续讨论参数变化带来的影响。

一句话概括：用大衍数构造结构化稀疏校验矩阵，在Matlab里跑完整误码率仿真链路，再用三组对照实验回答"迭代次数、码率、码长分别怎么影响LDPC性能"这个问题。

1.2 为什么用大衍数而不是随机构造

随机构造LDPC码的思路很直接：按行重、列重约束随机填0/1，然后检查并消去短环。这种方式的优点是灵活，缺点是每次跑出来的矩阵都不一样，性能有波动，而且校验矩阵需要完整存储。对仿真来说还能接受，但如果后续要做硬件实现，随机矩阵的存储和布线负担会非常大。

大衍数构造走的是"代数构造"路线。它用一组有明确数学规则的序列来生成循环移位值，再通过扩展得到校验矩阵。好处很明显：

• 确定性：同一套参数和序列，矩阵完全一致，实验可复现
• 存储友好：只需存基矩阵和扩展因子，不需要存整个大矩阵
• 易于控制结构：可以通过设计序列规避四环、六环等对迭代译码不利的结构
• 研究意义强：把古典数理序列与现代编码理论结合，是学术课题中比较讨巧的创新点

我个人的看法是，用大衍数构造本质上是提供了一种"看起来有规律、分布足够均匀"的循环移位值来源。它真正优秀的地方不在于神奇，而在于有确定性的排列逻辑，方便做性能分析和后续扩展研究。

1.3 Matlab在这个项目里的价值

选Matlab做这套仿真，核心原因是它把稀疏矩阵、逻辑索引、蒙特卡洛统计这些事都封装得很顺手。你不需要从零写一个稀疏矩阵库，也不需要手动管理大量循环，只要把算法思路表达清楚，剩下的交给矩阵运算。

而且Matlab的调试体验对这类研究型仿真特别友好。你可以逐步打印校验矩阵的稀疏结构，可以快速检查某个码字是否满足校验方程，可以用semilogy直接画出误码率曲线。对于需要反复调整参数、观察趋势的课题来说，这种交互式工作流非常重要。

当然，裸的Matlab代码在仿真大码长、高信噪比时会比较慢，这个我在后面的优化章节专门聊。但作为算法验证和课题研究工具，它是完全够格的。

2. 核心原理：大衍数如何变成稀疏校验矩阵

这一章是项目的理论地基。很多人在这一步卡住，就是没搞明白"大衍数"和"校验矩阵"之间是怎么搭上关系的。

2.1 大衍数的数学背景与序列化方法

"大衍之数五十，其用四十有九"出自古代典籍《周易·系辞传》。这个数理体系里最常被引用的数字包括：大衍之数50、其用数49、天数25（1+3+5+7+9）、地数30（2+4+6+8+10）、天地之数55。这些数字的特点是——离散、确定、相互之间有加减关系，天然适合当作构造序列的种子。

在LDPC构造里，大衍数不是直接拿来用的，而是先要转换成一个"循环移位值序列"。常见思路是：把上述经典数字作为种子集，通过某种扩张或交织规则生成长序列，再对扩展因子Z取模，映射到[0, Z-1]区间。

比如一个比较直观的序列化方案：

% 大衍数种子序列（示意） seed = [50, 49, 25, 30, 55, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; % 用线性同余方式扩张，取模后均匀分布 Z = 81; seqLen = 32; shiftSeq = zeros(1, seqLen); shiftSeq(1) = mod(seed(1), Z); for n = 2:seqLen shiftSeq(n) = mod(shiftSeq(n-1) + seed(mod(n-1, length(seed)) + 1), Z); end

这只是无数种实现中的一种，不同文献对种子选取、扩张方式都有自己的设计。但共同原则是：取模后的序列要尽量均匀覆盖[0, Z-1]，避免大量重复，这样展开后的校验矩阵才不容易出现短环。

注意：大衍数的文化背景了解一下即可，不要过度神话。它就是一组有确定分布特征的整数序列，在编码理论里等价于一种代数构造工具。

2.2 准循环基矩阵与单位矩阵扩展

拿到循环移位序列之后，下一步是排列成基矩阵，再做扩展。设基矩阵B的大小是J行K列，每个非零元素是一个移位值s，用Z×Z的单位阵循环右移s位替换；元素为-1或某个约定值时，用全零矩阵替换。于是得到大小为(J·Z)×(K·Z)的校验矩阵H，码长N=K·Z，校验位数量M=J·Z，设计码率R=(K-J)/K。

举个具体例子：取J=4、K=8、Z=81，那么H是324×648的矩阵，码长N=648，码率R=(8-4)/8=1/2。这个参数组合在LDPC研究里很常见，适合做迭代次数对比实验。

Matlab里扩展这一步的代码可以这样写：

function H = qc_ldpc_matrix(baseB, Z) [J, K] = size(baseB); H = zeros(J*Z, K*Z); for r = 1:J for c = 1:K s = baseB(r, c); if s >= 0 % 单位阵循环右移s位 I = circshift(eye(Z), s, 2); H((r-1)*Z+1 : r*Z, (c-1)*Z+1 : c*Z) = I; end end end H = logical(H); % 转成逻辑稀疏矩阵，节省内存 end

实际仿真时建议用稀疏矩阵存储，尤其码长上到几千之后，全零矩阵会占大量内存。Matlab里可以用sparse优化，但逻辑矩阵本身对仿真规模来说基本够用。

稀疏性怎么体现？每个校验方程只涉及少量变量节点。比如扩展前基矩阵每行非零元素个数是λ，扩展后每行权重仍然是λ，而码长是K·Z，λ远小于码长，这就是"低密度"的含义。正是这种稀疏性，让迭代译码成为可能。

2.3 避免短环的代数约束：四环检测与序列重排

短环是LDPC码性能的最大杀手之一。尤其是四环——两个变量节点连到同一对校验节点，会让消息在迭代中来回传递、自我加强，导致误码率出现平台，难以继续下降。

四环出现的代数条件可以精确表达：在准循环矩阵里，如果基矩阵存在两行r1、r2和两列c1、c2，满足：

(baseB(r1,c1) - baseB(r1,c2)) ≡ (baseB(r2,c1) - baseB(r2,c2)) (mod Z)

那么展开后的矩阵必定包含一个四环。所以构造基矩阵时，要保证任意两行、任意两列的移位差在模Z意义下不重复。

用大衍数序列构造时，我习惯先把候选序列填进基矩阵，然后跑一遍四环检测脚本，发现有冲突就调整序列排列顺序，或者更换种子扩张参数。这个过程是"构造-检测-微调"，不是一蹴而就的。

function hasCycle4 = check_cycle4(baseB, Z) [J, K] = size(baseB); hasCycle4 = false; for r1 = 1:J-1 for r2 = r1+1:J diffMap = containers.Map(); for c = 1:K s1 = baseB(r1, c); s2 = baseB(r2, c); if min(s1, s2) < 0 continue; end d = mod(s1 - s2, Z); key = num2str(d); if diffMap.isKey(key) hasCycle4 = true; return; end diffMap(key) = c; end end end end

这个脚本的重点是理解"模Z差值重复就出四环"这个逻辑。它在整个项目里属于质量门禁，每次改完基矩阵都要跑一遍。

3. Matlab仿真链路搭建与关键参数设计

理论讲完，进入实战。这一章给出完整的仿真链路结构，并且把最容易出错的能量换算、对数域译码实现讲透。

3.1 仿真链路总览：从编码统计到误码率

完整的BPSK-AWGN-误码率仿真链路长这样：

1. 生成随机信息比特序列m
2. 用生成矩阵G编码得到码字c
3. BPSK调制：0映射为+1，1映射为-1
4. 加高斯白噪声得到接收序列y
5. 计算对数似然比LLR，作为译码器初始输入
6. 运行对数域置信传播（BP）迭代译码
7. 迭代结束后的硬判决结果与原始信息比较，统计误码率

其中编码这一步，通常先把校验矩阵H通过高斯消元转换成系统形式，再求生成矩阵G。H本身由大衍数构造，往往不含完整的单位阵，需要做列变换。用Matlab的rref或者自写高斯消元都能完成，但要注意转换成系统形式后，矩阵不再保持稀疏，编码复杂度会上升。对码长在几百到几千的仿真来说，这个开销完全可以接受。

如果对编码效率要求更高，可以改用基于H的编码算法，把H化成近似下三角结构，用后向递推计算校验位。这个属于进阶优化，第一次跑通链路不建议上来就搞。

3.2 Eb/N0与噪声方差的换算：一个隐蔽的大坑

误码率仿真中最容易翻车的点，是信噪比换算。很多初学者直接把Eb/N0的dB值当作信噪比SNR用，结果曲线整体右移，性能看起来"差得离谱"。

正确关系是这样的。BPSK调制下，每个码字符号的能量Ec等于一个比特能量Eb乘以码率R的倒数：

Ec = Eb / R

因为一个信息比特对应1/R个码字比特。在Matlab里加噪声，用的是符号信噪比折算出来的噪声方差：

sigma2 = 1 / (2 * R * 10^(EbN0_dB / 10))

其中10^(EbN0_dB/10)是把dB转成线性值。以R=1/2、Eb/N0=2dB为例：

EbN0_dB = 2; R = 1/2; EbN0_lin = 10^(EbN0_dB / 10); sigma2 = 1 / (2 * R * EbN0_lin);

算出来sigma2约等于0.315。如果忘了除R，算出来是0.315的两倍即0.63，加噪功率变大，BER曲线就会往右移约3dB。

所以我在代码里永远单独写一个函数做换算，绝不手算，也绝不在循环里直接拼公式。

function noiseVar = ebno2noisevar(EbN0dB, R) noiseVar = 1 ./ (2 * R .* 10.^(EbN0dB / 10)); end

3.3 对数域置信传播译码的实现要点

译码是整个仿真链路的核心计算模块。我推荐用对数域BP，因为概率域BP在迭代中要反复做乘法，动态范围大，很容易数值下溢；对数域把所有消息换成对数似然比LLR，乘除变成加减，稳定性好得多。

对数域BP的初始化是：

% 接收信号y，噪声方差sigma2 Lch = 2 * y / sigma2; % 信道LLR，BPSK+AWGN

校验节点更新的标准公式是tanh规则：

L(r_mn) = 2 * atanh( prod( tanh( L(q_m'n) / 2 ) ) )

实际工程里，只要矩阵行列重不太小，完全可以用min-sum近似：

L(r_mn) ≈ prod( sign(L(q)) ) * min( |L(q)| )

这个近似省去大量tanh和atanh运算，性能损失通常不到0.2dB，在码长较大的时候尤其划算。

变量节点更新则相对直接：

L(q_mn) = Lch(n) + sum_{m' in M(n), m'≠m} L(r_m'')

每次迭代后计算后验LLR：

LQ_n = Lch(n) + sum_{m in M(n)} L(r_mn)

硬判决规则：LQ_n > 0判0，否则判1。

迭代终止条件有两个：当前硬判决码字满足所有校验方程，即H*c^T = 0（模2意义下）；或者达到预设的最大迭代次数。满足前者说明译码成功，可以提前退出，能省不少时间。

核心循环的伪代码结构如下：

for iter = 1:maxIter % 校验节点更新（min-sum近似） Lr = update_check_nodes(H, Lq, Lr); % 变量节点更新 Lq = update_var_nodes(H, Lch, Lr); % 后验LLR与硬判决 LQ = Lch + sum(Lr, 1); hard = double(LQ < 0); % 判为1 % 校验方程 syndrome = mod(H * hard.', 2); if all(syndrome == 0) break; end end

注意H如果存成稀疏矩阵，H * hard.'这一步要用mod(full(...))或者用稀疏矩阵乘法再mod避免类型问题。Matlab里稀疏矩阵乘逻辑向量是支持的，但结果有时是double类型，取mod没问题。

3.4 蒙特卡洛误码率统计的实操细节

误码率仿真的本质是蒙特卡洛实验：随机发大量帧，统计错误比特比例。但这里有个执行层面的讲究——怎么分配每帧试验次数，才能即跑得快又不浪费算力。

我的经验是：每个信噪比点，至少收集50~100个错误帧再停止；如果错误帧数不够，就继续跑直到最大帧数（比如5万帧或10万帧）。这样在高信噪比区域，曲线不会因为统计太少而剧烈抖动。低信噪比区域，错误非常频繁，几百帧就能达到统计量，跑得飞快。

我习惯在代码里加两层保护：

maxFrames = 50000; minErrFrames = 100; for snrIdx = 1:length(EbN0dB_list) errBits = 0; totalBits = 0; errFrames = 0; for frame = 1:maxFrames % 生成、编码、调制、加噪、译码 % 统计错误比特 errFrames = errFrames + (errBitsInFrame > 0); if errFrames >= minErrFrames && totalBits > 1e6 break; end end ber(snrIdx) = errBits / totalBits; end

这种动态截止策略，是我在实际仿真里反复调整出来的平衡点。既不会低噪底时磨洋工，也不会高噪底过早收工导致曲线尾端波动大。

4. 三组对比实验：迭代次数、码率、码长如何影响误码率

链路跑通之后，最核心的对比实验就来了。下面给出每组实验的配置、实测趋势和解读方式。数值是我在自己实现里观察到的典型量级，不同构造矩阵会略有出入，但趋势稳定。

4.1 实验一：不同译码迭代次数对误码率的影响

固定码长N=648、码率R=1/2、基矩阵和扩展因子不变，只把最大迭代次数设为1、5、10、20、50，分别统计BER曲线。

实测下来的典型结果是：在误码率10的-4次方量级附近，迭代1次和迭代5次之间差距非常大——1次迭代基本没有纠错效果，BER曲线斜率很小；5次迭代已经能看出瀑布区雏形；但从10次增加到20次，性能改善明显放缓；20次到50次，曲线几乎重合。

这个趋势背后的原因在于，LDPC迭代译码的收敛过程是指数速率衰减的。前面若干次迭代把大部分错误纠正掉，后续迭代只是在处理少量顽固错误。迭代次数增加带来的增益曲线必然是饱和的。

所以当你做这个对比实验时，解读重点应该放在"边际收益递减"上。实际系统里选择迭代次数不是越大越好，而要看硬件延时和功耗约束。理解这个trade-off，才是这个实验的真正价值。

4.2 实验二：不同码率对误码率的影响

固定码长N=648，让码率分别取1/2、2/3、3/4。注意改变码率意味着改变基矩阵的行列比，需要重新构造基矩阵，而不是简单删几行。构造时保证每种码率都无四环。

实测趋势是：码率越低，相同Eb/N0下误码率越低。这不是LDPC码的特例，而是编码理论的一般规律——冗余越多，纠错能力越强。在BER=10的-4次方处，典型量级大概是：R=1/2约需1.7dB，R=2/3约需2.4dB，R=3/4约需3.1dB，逐级差0.6~0.8dB。

这个实验的解读要强调"公平性"：比较不同码率要在相同Eb/N0下比，即每信息比特消耗相同能量的前提。如果换成符号信噪比Es/N0，低码率因为符号能量高，曲线优势会更大，这个换算必须说清楚。另外，高码率虽然有性能劣势，但频谱效率高，实际系统中往往需要结合调制阶数、吞吐量需求做联合选择。

4.3 实验三：不同码长对误码率的影响

固定码率R=1/2，码长N取162、648、2048三档。这里有个构造技巧：保持基矩阵结构不变，只改变扩展因子Z，可以快速得到不同码长且结构同源的矩阵。比如Z从27变到81再变到256。

实测现象最明显的是两件事：

第一，长码的瀑布区更陡峭，即误码率下降更快，越发接近信道容量极限。N=2048时，瀑布区从1.5dB开始急速下降；N=162时，瀑布区被明显拉宽，需要更高信噪比才能进入快速下降区。

第二，短码在误码率10的-3次方附近出现错误平台。这个平台是有限长效应和短环共同作用的结果，增加信噪比很难继续压低BER。长码虽然理论上平台更深，但实际在中等信噪比区间不会碰到。

这个实验的解释逻辑要落到"码长决定码字距离分布"上。长码有更大的最小距离和更好的距离谱，错误的错误图样种类更多但权重更高，译码器更容易区分。短码则受限于有限码长效应，性能离理论极限明显偏出。

4.4 三组实验的结果解读与曲线绘制

画误码率曲线，最标准的做法是用semilogy画BER对Eb/N0，也可以用loglog看双对数下瀑布区的斜率差异。每条曲线至少跑3个信噪比点，理想情况下瀑布区要覆盖2~3个数量级的BER下降。

看BER曲线时，我习惯把曲线分成三段来看：

• 高误码率区域（10的-1次方以上）：主要反映码字距离结构的整体特征
• 瀑布区（BER快速下降的中段）：这是迭代译码纠错能力最强的工作区域，斜率越陡说明码越好
• 错误平台（BER下降变缓甚至停滞）：说明译码陷入局部最优，短环或有限长效应开始主导

三组实验正好分别对应这三段现象：迭代次数影响整个曲线的饱和高度，码率决定瀑布区起点位置，码长决定瀑布区陡峭程度和平台位置。把这三张图放在一起看，LDPC码性能的主要规律就都齐了。

注意：画图时别忘了在图例和横轴标注里写清楚仿真参数。我见过太多人跑完仿真，半年后回看代码，死活想不起某个曲线用的是多少行重、多少次迭代。仿真代码里命名规范，图表标题写全，会省去大量重复劳动。

5. 仿真中的常见问题与排查技巧

这一章是从多次踩坑里总结出来的。每个跑LDPC仿真的人都可能遇到下面这些问题，有些问题排查起来相当隐蔽。

5.1 误码率曲线异常：错误平台与曲线"右移"

错误平台是LDPC仿真最经典的异常现象。如果BER曲线在10的-3次方左右突然变平，甚至反弹上去，优先怀疑三件事：校验矩阵是否有四环；迭代次数是否过小；数值实现是否有下溢或溢出。

四环可以通过2.3节的检测脚本排查。迭代次数如果低于5，出现高错误平台是正常的，这是译码不充分的体现。如果两个都排除了，再用高精度格式跑一遍看是否数值问题。

曲线整体右移超过1dB则要检查噪声方差换算。最常见的就是漏了码率因子，或者BPSK映射时把0映射成0、1映射成1、然后直接当QPSK用。LDPC仿真里0/1到±1的映射必须和LLR初始化公式严格一致。

5.2 仿真运行时间过长：从几天到几小时的优化手段

码长2048、信噪比3dB以上时，每帧跑50次迭代，一帧就要做几百万次消息更新。如果循环写得粗糙，仿真跑几天都正常。我一般按次序做三层优化：

第一层，对数域和min-sum近似。概率域BP的tanh计算是高耗时主因，换成min-sum直接提速5~10倍，性能损失在可接受范围内。

第二层，向量化校验节点更新。不要用双重循环逐行逐列更新校验消息，用矩阵切片整体计算。Matlab对矩阵运算有高度优化，同样的逻辑，向量化之后提速可以做到一个数量级以上。

第三层，并行化。不同信噪比点之间互不依赖，可以用parfor直接并行。我的做法是把每个Eb/N0点封装成一个函数，parfor里调用，然后用一个数组收集结果。注意随机数种子的设置要独立，避免多个worker生成完全相同的数据。

parfor snrIdx = 1:length(EbN0dB) ber(snrIdx) = ldpc_sim_one_point(ebno_list(snrIdx), params); end

5.3 矩阵构造有效性的快速自检清单

每次重新构造基矩阵，我都会按下面这个清单快速检查一遍，避免带病矩阵进仿真：

• H的维度是否正确：行数M = J·Z、列数N = K·Z
• H的行重是否接近基矩阵每行非零元素个数，列重是否一致
• H是否满行秩：用rank对较小的矩阵可以快速验证，满秩才能保证编码后码率正确
• 四环是否为零：跑cycle4检测
• 编码后是否满足H*c^T = 0：用10个随机码字抽样验证

这个清单10分钟就能跑完，但能挡掉90%的后续问题。

5.4 数值稳定性与调试技巧

对数域BP虽然规避了大范围乘法，但LLR本身在迭代中会快速增长。如果用的是单精度float，性能会明显下降。仿真中一律用double，并且在校验节点更新时对消息幅值做截断处理：

Lr = sign(Lr) .* min(abs(Lr), 50); % 截断阈值可调节

这种截断既防止LLR无限增长导致的NaN，也符合实际硬件中的消息量化思想，不会明显影响误码率。

调试阶段最有用的技巧是固定随机数种子。每次随机生成信息比特或噪声之前，统一用rng(2024, 'twister')之类的设置。这样一旦发现某个信噪比点出现不合理结果，可以重新运行到同一帧，逐条消息对比调试。没有这个习惯，数小时的仿真一出问题只能从头再来。

最后分享一点个人体会

这个大衍数LDPC仿真项目跑完，我最深的感觉是：LDPC的研究虽然门槛不低，但用Matlab完整复现一个"构造-编码-译码-分析"闭环，对整个编码理论的理解提升是巨大的。很多人只会在论文里看BER曲线，却不知道一条曲线背后要经过多少次蒙特卡洛帧的迭代、多少个参数之间的相互牵制。

如果让我给正在做类似课题的朋友一个建议，我会说：先别急着上大码长、高信噪比。用一个小矩阵，比如N=162、R=1/2，把全链路跑通、曲线画出来，再逐步放大规模。小矩阵出错容易排查，跑一轮也就十几秒，改参数试错的成本极低。链路稳了以后，再上N=2048去逼近理论极限，你会明显感受到大码长带来的性能飞跃。

至于扩展方向，这套东西的理论空间还很大，比如结合密度进化分析大衍数序列的译码门限、或者把min-sum替换成归一化min-sum看性能变化、再往下还可以讨论硬件量化位宽对误码率的影响。底层的构造和仿真框架不变，换的是更细粒度的问题视角。这就是我觉得这个项目最值得花时间的地方——它给后面一系列研究铺好了可复用的地基。