☰
基于大衍数构造LDPC校验矩阵的误码率仿真与Matlab实现
2026/10/16 7:50:55 网站建设 项目流程

前段时间在做一个信道编码相关的仿真验证,需要对比不同译码迭代次数、码率和码长下LDPC码的误码率表现。一开始我图省事,直接用随机方法生成稀疏校验矩阵,结果短环问题、编码复杂度问题接踵而至,仿真出来的曲线怎么看都不对劲。后来偶然想到用“大衍数”这种确定性序列来构造校验矩阵,把非零元素的位置和循环移位量都交给序列控制,效果出奇地稳,短环能被有效抑制,码的构型也便于复现。这篇就基于这个思路,把整条LDPC误码率仿真链路完整复盘一遍,从矩阵构造原理到Matlab实现细节都讲清楚,希望能给正在做纠错码仿真或者LDPC性能研究的同学一些参考。

这套仿真重点做三件事:一是用大衍数构造稀疏校验矩阵H并完成LDPC编码,二是在AWGN信道下配合BPSK调制做蒙特卡洛误码率统计,三是分别扫描译码迭代次数、码率、码长三个维度,观察误码率曲线的变化规律。适合刚接触LDPC想跑通闭环仿真的人,也适合已经会用现成工具箱、但想尝试自定义矩阵构造方式的进阶读者。

1. 整体设计与思路拆解

1.1 为什么选择“大衍数”来构造校验矩阵

LDPC码的性能很大程度上被校验矩阵H的“长相”决定。H要足够稀疏,每行每列的非零元素不能太多,同时要尽量避免短环,尤其4环是译码性能的杀手。短环会让置信传播迭代时消息在环里打转,置信度互相强化,等于算法自己在骗自己。

传统做法是用完全随机的稀疏矩阵,靠天吃饭,不行就重新抽一版。Mackay的随机构造法和PEG算法也都在行业里用得很多,它们都能构造出性能不错的矩阵,但问题是:第一,随机性引入后每次复现结果不一致,调试起来很麻烦;第二,在Matlab里做大规模随机稀疏矩阵生成,还要做逐位短环检查,计算量不小。

大衍数这个思路有意思的地方在于,它绕开了“随机生成-检查-重来”的循环,改用一条确定性序列直接决定矩阵里非零元素的位置规律。“大衍之数五十”的说法出自古籍筮法,实际应用时我借用了“五十”这个模数概念,再配合“体”与“用”的映射思想,把序列映射到循环移位量上。这样做的好处非常直接:同一组参数生成的H矩阵永远一样,完全可复现,而且因为序列本身有规律,可以在构造阶段就尽量避开密度不均匀的坑。

用大衍数序列构造校验矩阵还有一种优势,就是天然能和准循环LDPC(QC-LDPC)结构结合。只要把序列当作基矩阵B中各个位置“1”的循环移位系数,再用大小为p的单位阵循环移位去扩展,就能得到一块块循环子矩阵组成的校验矩阵。准循环结构在编码端可以靠移位寄存器实现线性复杂度编码,这比普通稀疏矩阵做高斯消元编码要快得多,对长码长仿真来说尤其重要。

1.2 三个仿真维度的设计逻辑

我把迭代次数、码率、码长这三件事放在一起做对比,是因为它们分别对应LDPC译码的“计算深度”“冗余强度”和“结构尺度”,三者相互影响又各有侧重。

迭代次数直接决定置信传播算法的消息扩散范围。迭代次数太少,校验节点和变量节点之间的信息还没充分交换,一些本可以被纠正的比特只能带着错误信念继续传递,最终导致误码率居高不下。迭代次数太多也有问题,一是时间成本成倍增加,二是当置信度更新达到收敛后,多余迭代只会让消息在小范围环内反复震荡,有时甚至会因为数值偏离导致误码率抬升。所以我实际扫了5、10、20、50四个档位,就是为了观察性能和复杂度之间的平衡点在哪。

码率的设计则直接作用于系统的冗余程度。码率低,校验位占比高,纠错能力自然更强,但频谱效率也要付出代价。码率升高后,冗余减少,在同样的信道条件下误码率曲线会明显右移。仿真中我会通过固定扩展因子p、调整基矩阵行数来实现不同码率,这样得到的码字结构风格一致,对比起来更加公平。

码长的影响更偏向宏观。理论上无限长LDPC码才能逼近香农极限,有限码长下都会出现性能损耗,码长越短损耗越明显。但码长翻倍之后,H矩阵尺寸按平方增长,Matlab里无论是构造、编码还是译码,耗时都会快速增长。考虑到普通PC的性能限制,我只取了几档中等长度的码来做对比,比如576、1152、2304,这样既能看出长度趋势,又不至于跑一晚上只出两帧结果。

2. 核心细节解析与实操要点

2.1 大衍数构造稀疏校验矩阵的落地流程

把“大衍数”这个文化色彩很浓的概念变成可运行的算法,关键在于把抽象意象转化为数值操作。我采用的方案如下:

首先确定基矩阵的尺寸。设基矩阵B大小为m_b行、n_b列,其中m_b决定校验位的块数,n_b决定变量节点的块数。假设扩展因子为p,则最终校验矩阵H的大小是(m_bp)行、(n_bp)列。码率R可以由1减去(m_b/n_b)得到近似值,前提是去掉线性相关行,实际上构造时尽量让基矩阵行满秩。

其次生成大衍数序列。我取模数p=50,呼应“大衍之数五十”,序列形式可以采用一个简单的二次同余式:

p = 50; % 扩展因子,也是大衍模数 k = (0 : p - 1)'; a = 7; b = 23; c = 11; s = mod(a * k.^2 + b * k + c, p); % 生成大衍数序列 s = s + 1; % 映射到1..p,方便Matlab索引

这里选择二次式而不是简单线性式,是为了让序列在模p下分布更均匀,减少连续值过多造成的矩阵局部聚集。系数a、b、c可以微调,只要确保生成序列在0到p-1范围内差不多均匀覆盖就行。构造时不需要把所有50个值都用上,可以只挑序列的前面部分,对应“其用四十有九”的取用思想,这也给行重和列重的调节留出了余地。

有了序列s之后,基矩阵B中每个“1”的位置就对应一个循环移位量。比如B中第i行、第j列如果是1,那么H中对应分块就是一个p×p的单位阵循环右移s(t)位,t来自序列按顺序取出的索引。最终H可以按分块形式组装:

H = zeros(m_b * p, n_b * p); for i = 1 : m_b for j = 1 : n_b if B(i, j) == 1 shift = s(mod(t, p) + 1); % t为当前非零计数 block = circshift(eye(p), shift, 2); % 循环右移 rowIdx = (i - 1) * p + 1 : i * p; colIdx = (j - 1) * p + 1 : j * p; H(rowIdx, colIdx) = block; t = t + 1; end end end

这段代码只是构造阶段的核心骨架,实际使用中还需要根据目标码长和码率确定m_b、n_b和p的组合。比如总码长N = n_b * p,校验位M = m_b * p,码率就是1 - M/N。我一般先定N和R,再倒推基矩阵尺寸,最后调整p保证整除,避免出现非完整分块。

构造完成后一定要做两件检查:一是矩阵的秩,二是围长。秩可以通过rank(full(H))配合mod 2处理来检查,如果秩不等于m_b*p,说明存在行相关性,编码时生成矩阵会出问题。围长检测我写过一个小函数,用广度优先搜索从每个变量节点出发寻找最短返回路径,路径长度小于4就说明出现了4环,需要换一组大衍序列系数或调整基矩阵结构。

2.2 关键参数设计与译码算法选型

译码端我首选用对数似然比域上的置信传播算法,也就是LLR-BP,业内常称为和积算法SPA。核心思路是变量节点和校验节点之间来回传递消息,消息是对数似然比,正负号表示硬判决方向,绝对值大小表示置信程度。

初始化时给出信道观测的对数似然比,BPSK调制的映射通常选为0对应+1,1对应-1,于是LLR初始值可以写为:

L0 = 2 * receivedSignal ./ (noiseVar);

其中noiseVar是噪声方差。要注意这里默认发送符号是+1和-1,如果用户习惯0/1映射,符号关系要对齐。

校验节点更新是SPA里最容易写错的地方。LLR域下的标准形式包含双曲正切函数连乘,实际仿真中为了防止数值下溢,建议把双曲正切连乘改成基于符号和最小幅度近似的形式。单纯用最小和算法MSA实现起来更简单,但性能会比标准SPA差大约0.2到0.3dB。我在工程上通常取折中,使用归一化最小和算法,把最小幅度乘以一个0.75左右的归一化因子,既能保留大部分性能,又能避免Matlab里大量调用tanh带来的耗时。

迭代终止条件也在方案设计阶段就要想好。最常用的做法是每次迭代结束后对所有校验方程做一次硬判决校验,如果得到的码字乘以H转置模2全零,就提前退出,不用再继续迭代。这个early termination对长码长仿真特别重要,能省掉大量无用计算。我在代码里用mod(H * decodedHard, 2)是否等于零向量来判断,整个逻辑非常直接。

参数设计上需要保持各对比组之间的公平性。研究迭代次数时,固定码率和码长,只复制译码循环次数;研究码率时,固定码长和迭代次数,调整校验位块数;研究码长时,固定码率和迭代次数,按倍数增大p和基矩阵列数。这种控制变量思路能确保最终差异单纯来自目标参数,而不是混合因素。

3. 实操过程与核心环节实现

3.1 Matlab仿真环境与整体代码架构

整个仿真我在Matlab环境中跑通,不需要额外安装通信工具箱也可以写出核心代码,但如果有通信工具箱,部分调制和误码率统计函数可以省事一些。项目文件按功能拆成了五个模块,目录结构如下:

+ldpc_sim/ - configParams.m % 全局参数配置 - buildHMatrix.m % 基于大衍数构造校验矩阵 - ldpcEncode.m % 简单编码函数 - ldpcDecode.m % LLR-BP或归一化MSA译码 - runSimulation.m % 蒙特卡洛主循环

参数配置模块统一管理码长、码率、迭代次数、信噪比范围、仿真帧数上限,方便一次跑完所有对比。为了让每个信噪比点都拿到足够置信度的统计结果,我设了最低错误帧数阈值,比如累计正确帧数超过100万帧或者错误帧数超过200帧才停止当前信噪比的统计,否则继续仿真直到达到上限。

编码模块的设计在LDPC里通常会卡一下。若直接对H做高斯消元求生成矩阵G,Matlab里稀疏矩阵消元速度尚可,但消元过程中产生的列置换会让系统码字结构变得不直观,还可能出现某些列相关性导致G无法直接得到。我的做法是在构造基矩阵时就刻意让右侧区域接近单位阵形式,基矩阵右侧安排一个m_b×m_b的三角或对角块结构,这样H扩展后仍保住一个可逆的块结构,编码时只要做分块矩阵乘法和回代运算,不用经历全套通用消元。

如果用户不想那么在构造上花心思,也可以直接用高斯消元写一个通用编码函数,但一定要做模2运算,并记录列置换向量,否则编码出来的码字校验结果对不上。这个坑我踩过多次,后面会专门说。

3.2 校验矩阵构造与编译码核心代码逻辑

构造函数的核心是生成大衍数序列、组装基矩阵、再循环扩展。下面是完整度较高的一段逻辑:

function H = buildHMatrix(p, mBlocks, nBlocks, rowWeight, seqParams) k = (0 : p - 1)'; a = seqParams(1); b = seqParams(2); c = seqParams(3); s = mod(a * k.^2 + b * k + c, p) + 1; B = zeros(mBlocks, nBlocks); t = 1; for i = 1 : mBlocks % 每行放rowWeight个1,位置用序列s的索引决定 idx = mod(s(t : t + rowWeight - 1), nBlocks) + 1; B(i, idx) = 1; t = t + rowWeight; end H = sparse(mBlocks * p, nBlocks * p); for i = 1 : mBlocks for j = 1 : nBlocks if B(i, j) == 1 shift = s(mod(t + i * j, p) + 1); block = sparse(circshift(eye(p), shift, 2)); rowIdx = (i - 1) * p + 1 : i * p; colIdx = (j - 1) * p + 1 : j * p; H(rowIdx, colIdx) = block; end end end end

这里每一行放置rowWeight个1的位置也由大衍序列控制,而不是随机选,这样整体构造完全脱离随机数生成器。注意shift的计算为了避免序列索引重复,我加了一个i*j的扰动,保证不同分块的移位量存在差异。这种细节如果不处理,可能出现相邻两块循环移位量相同,导致H在局部产生规则性过强的结构,影响girth表现。

编码端基于系统形式实现,伪代码为:

function codeword = ldpcEncode(message, G) % G为k行n列的生成矩阵,message为1xk二进制向量 codeword = mod(double(message) * G, 2); end

如果H构造时保留了右侧块状结构,可以手动算出系统形式的校验部分,进而得到生成矩阵G。比如H可以写成[A B],其中B为可逆方阵,那么奇偶校验位p_vec可由message * A' * inv(B')的模2结果得到。为了避免在Matlab里对稀疏矩阵做逆运算,我一般用稀疏LU分解或者直接用B'\A'的方式求解,再取模2。

译码端是整套仿真最耗时的模块,核心逻辑如下:

function decodedBits = ldpcDecode(LLR0, H, maxIter, normFactor) [M, N] = size(H); [rowIdx, colIdx] = find(H); % 初始化变量节点消息 VMsg = LLR0; for iter = 1 : maxIter % 校验节点更新 for e = 1 : length(rowIdx) r = rowIdx(e); c = colIdx(e); neighbors = find(H(r, :)); others = neighbors(neighbors ~= c); Msg = VMsg(others); signProd = prod(sign(Msg)); minAbs = min(abs(Msg)); CMsg(e) = signProd * minAbs * normFactor; end % 变量节点更新 for e = 1 : length(rowIdx) r = rowIdx(e); c = colIdx(e); neighbors = find(H(:, c)); others = neighbors(neighbors ~= r); VMsg(c) = LLR0(c) + sum(CMsg(others)); end % 硬判决与提前终止 decodedBits = double(VMsg < 0); if mod(H * decodedBits', 2) == 0 break; end end end

这段代码只是教学级实现,实际跑大规模矩阵时用find和循环遍历会非常慢,尤其是列数达到2304甚至更高时。工程优化思路是预先计算每个节点的邻居索引表,把校验节点更新和变量节点更新改成矩阵批量操作,或者直接写MEX文件。但在仿真量不大时,纯Matlab循环加提前终止也能接受。

3.3 蒙特卡洛仿真流程与结果整理

主循环的核心结构是两层遍历,外层遍历信噪比点,内层逐帧仿真:

for snrIdx = 1 : length(snrAxis) errBits = 0; totalBits = 0; errFrames = 0; frameCnt = 0; while errFrames < minErrFrames && frameCnt < maxFrames info = randi([0 1], 1, k); code = ldpcEncode(info, G); symbols = 1 - 2 * code; % BPSK映射 noiseSigma = sqrt(noiseVar); r = symbols + noiseSigma * randn(size(symbols)); % LLR初始化 LLR0 = 2 * r / (noiseSigma^2); est = ldpcDecode(LLR0, H, maxIter, normFactor); errBits = errBits + sum(est ~= code); totalBits = totalBits + length(code); errFrames = errFrames + any(est ~= code); frameCnt = frameCnt + 1; end BER(snrIdx) = errBits / totalBits; end

每次运行前我建议用固定随机数种子设置随机流,保证不同参数组之间只有目标变量在变,随机噪声序列保持一致,这样出的曲线对比才公平。这一点很容易被忽略,但实验结果表明,如果每组参数都重新生成完全不同的随机序列,即使理论上性能相同,仿真出的曲线也可能相差0.2dB以上。

结果整理时我习惯把误码率数据存成结构体数组,再统一画图。在对比不同码率或码长时,横轴用Eb/N0而不是Es/N0,因为码率不同导致每个信息比特实际消耗的能量不同。换算公式为:

EbN0_dB = EsN0_dB + 10 * log10(1 / codeRate);

BPSK下EsN0就是归一化信噪比,这个换算在低码率时差别能达到3dB左右,很多人画出来的曲线位置对不上往往就是这里出了问题。

4. 不同维度下的误码率对比分析与解读

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

固定码率为0.5、码长为2304时,扫描迭代次数可以看到明显的性能分层。我取的档位是5、10、20、50,从结果趋势看,5次迭代时误码率曲线整体靠右,瀑布区的倾斜度也偏缓;提升到10次后曲线左移约0.5dB到0.8dB;继续提升到20次,增益只剩0.1dB左右;而50次迭代相比20次几乎没有肉眼可见的改善,反而在低信噪比区域偶尔出现误码率略抬升的现象。

这个现象背后就是消息扩散饱和与环结构双重作用的结果。迭代初期,消息沿着Tanner图向外扩散,每多迭代一次就多覆盖一层邻居节点,纠错能力快速增强;等扩散范围覆盖了整个连通区域,置信度更新进入收敛状态,额外迭代只会反复修正已经被确认的消息,增益趋于零。在存在短环的位置,消息经过环路返回自身时会带上过度自信的偏差,迭代次数越多这种偏差被放大得越明显,这就解释了为什么高迭代在部分信噪比点反而变差。

实际操作中的启示是,仿真时不要盲目把迭代次数拉到50甚至100,不仅耗时,结果也不一定更好。做初步扫参时用10次迭代足够看趋势,确定最优参数后如果想看到更精致的性能上限,再单独用20或30次迭代验证。

4.2 码率变化导致的性能差异

保持码长为2304、迭代次数为10,我对比了1/4、1/2、3/4三档码率。低码率1/4时冗余极多,纠错能力最强,误码率曲线在较低的信噪比就开始明显下降;码率1/2处于中间档,曲线位置比1/4右移了大约1.5dB到2dB,但频谱效率翻倍;码率3/4时曲线进一步右移,且瀑布区的斜率明显变缓,说明在高码率下有限码长LDPC的性能损失更加严重。

码率影响误码率的本质可以从两方面理解。第一是码字最小距离的变化,码率越高,码字集合在汉明空间中的分布越“拥挤”,最小距离相对较小,纠正随机错误的能力自然下降。第二是Tanner图密度,高码率意味着校验方程相对较少,变量节点从校验节点获得的外部信息更少,译码消息的约束强度不足,收敛速度变慢。

实际系统选码率时不能只盯住误码率曲线,还要把吞吐率和时延考虑进来。信道条件恶劣时1/4码率能带来明显增益,但有效信息速率只有符号速率四分之一;信道质量好时3/4码率显然更划算。仿真结果给我的参考是:1/2码率通常是最稳妥的中间选择,性能与效率兼顾。

4.3 码长对性能与仿真成本的双重影响

固定码率1/2、迭代次数10,对比576、1152、2304三个码长后发现,码长越长性能越逼近理想极限。2304码长的瀑布区位置比1152码长左移约0.3dB,1152又比576左移约0.4dB。这种增益在低信噪比区域更明显,越靠近香农限,长码的优势越突出。但代价也同样直观:2304码长的单帧仿真时间大约是576码长的15到20倍,因为矩阵尺寸变大让校验节点更新的消息规模成倍增长。

码长增加带来性能提升的原因在于有限长效应减弱。短码条件下,码字之间的最小距离受限,并且校验矩阵中的短环更容易形成,译码消息更容易陷入“自增强陷阱”。码长越长,H矩阵的围长理论上可以设计得越大,译码近似最优的概率就越高。

仿真成本会成为现实的瓶颈。576码长可以在几分钟内跑完一个信噪比的统计,2304码长就可能需要半小时以上。我实际跑的时候会把每个信噪比的最少错误帧数从200降到100,以换取时间;同时在迭代收敛后尽早终止,很多信噪比点实际只用2到3次迭代就满足校验约束,这样能大幅缩短总耗时。

4.4 综合对比与参数选择建议

把三维度放在一起看,可以从一组典型趋势中得到选择思路:

参数方向对性能的影响对仿真成本的影响实际选择建议
迭代次数增加先快速提升后饱和线性增长10~20次足够
码率降低性能提升明显校验位增加,编译码略增平衡信道的1/2左右
码长增加接近香农限平方级增长按时延预算选择

这些趋势看似简单,真在工程里做取舍时还要结合系统延时预算和硬件资源。我自己的经验是,先把迭代次数定在能力平衡点,再根据目标吞吐量选码率,最后根据可容忍的处理延迟选码长,这样比单纯追求最优误码率更符合实际落地需求。

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

5.1 问题速查表

仿真中几个高频问题整理成速查表,基本覆盖大多数人的卡点:

问题现象可能原因解决办法
编码后H乘以码字不为零H存在行相关性,生成矩阵G出错检查rank,确保秩等于校验位数;改用系统形式编码
误码率在低信噪比区域出现平台LLR初始化符号方向相反检查BPSK映射符号与硬判决阈值是否对齐
高信噪比下误码率不再下降矩阵存在4环或短环换大衍数序列系数,增大扩展因子或调整基矩阵行重
迭代次数提高后性能反而变差消息在短环中过度自信检查girth,增加缓存策略,限制最大迭代
仿真速度过慢纯循环实现,无提前终止增加校验终止判断,改用批量矩阵运算或MEX优化
不同SNR下曲线交叉严重随机噪声序列不一致固定随机种子,使用同一组噪声序列

5.2 独家避坑经验

第一个容易踩的坑是误比特率的统计口径。很多人译码输出后直接和原始信息比特比较,忽略了编码后码字和系统比特之间的位置对应关系。如果生成矩阵G经过了列置换,信息位不一定排在码字前k位,直接截断比较会得到完全错误的误码率。我后来统一改成编译码全程用同一套列索引变量,比较时也基于完整码字,先统计码字误码再换算信息位误码,才彻底解决这个问题。

第二个坑是噪声方差的单位换算。在BPSK+AWGN仿真里,噪声方差对应的是每符号噪声功率,但画BER曲线时横轴几乎都用Eb/N0。用户一旦把Eb/N0直接当作EsN0代入噪声生成公式,仿真曲线就会比理论位置偏移约1dB到3dB不等。我习惯先把Eb/N0换算成符号信噪比,再根据BPSK的符号能量为1推导噪声方差:

EsN0_dB = EbN0_dB + 10 * log10(codeRate); noiseVar = 10 ^ (-EsN0_dB / 10);

第三个经验是分组对比时务必保证“单变量”原则。我在同时扫码率和码长时发现,如果码长变了而码率只是近似相等,比如所谓的0.5码率一个取576/1152、另一个取600/1200,两者的校验方程数量和列重分布都不同,最终曲线差异根本分不清是码率影响还是码长影响。所以每次对比前先确认H矩阵的行重、列重、girth在同一风格下,再放心下手。

我还想提醒一点:大衍数构造的LDPC码不是万能的。它最大的优点是确定性强、方便复现、容易和准循环结构融合,但如果完全套用固定序列而不做任何girth优化,性能不一定能赢过精心调节的PEG矩阵。所以我把这套方法定位为“快速、稳定、可系统化对比的构造工具”,适合做参数扫描和教学验证,真正追求极限性能时还是要结合密度演进或围长优化来做针对性设计。

在使用大衍数构造矩阵时,最后再分享一个小心得:不要只把“五十”当成一个固定模数来用。我在不同码长下试过,扩展因子p取50到100之间的不同值,序列系数微调后,得到的矩阵性能会有细微差异。如果有耐心,可以在构造阶段多跑几组系数,用girth检测和短码仿真双重筛选,挑出最优的大衍参数组,这样整个LDPC仿真工程的根基会更扎实。对于刚上手的同学,建议先从固定参数跑通链路,再逐步加入自定义矩阵和性能对比,这样调试起来不会手忙脚乱。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询