LDPC译码算法MATLAB实现：软判决与硬判决对比

简介一套基于Matlab的LDPC编码与译码完整项目面向通信工程相关专业学生、科研人员及对信道编码有需求的开发者涵盖软判决与硬判决两大类译码算法便于对比不同算法下的误码率性能与实现复杂度。压缩包共8个文件以7个m脚本为主并含1个docx算法说明文档整体仅21KB轻量易用。核心功能包括稀疏校验矩阵生成、校验关系验证、概率域/对数域译码、比特翻转译码及简易对数域译码并提供BER仿真入口适合新手快速上手及有一定经验者做算法改造。目前已有319人学习下载代码经过实测校正若运行遇到问题可联系作者咨询调试性价比高。1. LDPC 译码算法的 MATLAB 实现先从软判决和硬判决的分岔说起LDPC 码的译码性能差距往往不在信道编码本身而在译码器对软信息的利用程度。同样是校验矩阵硬判决的比特翻转Bit Flipping在 BSC 信道下实现简单、延迟低但性能比软判决的对数域 BP 算法差一到两个数量级而概率域 BP 虽然数学形式直观直接实现却会遇到下溢和归一化开销。达摩老生这套 MATLAB 项目把三类译码路径放进同一份代码里还带上了 BSN 结构下的校验矩阵构建比较适合两类人一类是想把信道编码课上的公式落到可运行代码的学生另一类是需要在仿真链路里快速对比不同译码算法性能的工程师。我解压这份资源后注意到一个加分细节它没有把软判决和硬判决写进同一个函数再靠 flag 切换而是拆成了decodeBitFlip.m、decodeProbDomain.m、decodeLogDomain.m三个独立入口。这直接让误码率曲线的对比实验变成了「换函数名」级别的操作。接下来我按「构造校验矩阵 → 硬判决 → 软判决 → 性能仿真 → 进阶调试」的顺序拆解这份源码。2. makeLdpc 与校验矩阵构建LDPC 编码的起点是 H 矩阵而不是 G 矩阵LDPC 编码的常规误区是一上来就找生成矩阵 G。实际上makeLdpc.m这类工具的核心任务是先构造出满足稀疏性约束的校验矩阵 H再通过高斯消元把 H 变换成系统形式进而提取出 G。H 矩阵的行数对应校验方程个数列数对应码长稀疏性决定了 Tanner 图上短环的数量而短环直接制约 BP 译码的收敛质量。2.1 规则 LDPC 和非规则 LDPC 的参数约定常见的规则 LDPC 用 (n, j, k) 表示n 是码长j 是列重每列 1 的个数k 是行重每行 1 的个数且满足 j k。非规则 LDPC 则允许行列重量不一致性能通常更好但编码复杂度更高。这套代码里的makeLdpc.m我一般按以下方式调用% 构造一个码长 96、列重 3、行重 6 的规则 LDPC 校验矩阵 n 96; j 3; k 6; H makeLdpc(n, j, k);逻辑说明makeLdpc先建立一个全零矩阵然后对每一列随机放置 j 个 1再通过调整保证每一行 1 的数量等于 k。这个过程如果只靠随机填充大概率会出现行重不达标的情况所以函数内部通常会有一次重排修正。参数上列重 j 越小译码门槛越低但太小的 j比如 2会让最小距离退化行重 k 越大码率越高但校验约束更强错误平层也更高。2.2 从 H 矩阵生成 G 矩阵的系统化处理得到 H 之后编码需要 G。MATLAB 中常用mod(rref(H), 2)把 H 转化成行最简阶梯型再通过列置换换出单位阵。这里有一个关键坑LDPC 的 H 矩阵如果直接用rref会因为浮点舍入误差导致结果不是严格的 0/1 整数。在资源对应的源码中makeParityChk.m承担了从 H 提取校验位并检查是否满足mod(H * G, 2) 0的职责。它的输入就是makeLdpc产出的 H内部先做二值高斯消元再调整列顺序使得前若干列构成单位子矩阵。验证代码可以这样写% 验证 H 与 G 的正交性 G makeParityChk(H); % 返回生成矩阵系统形式 valid all(mod(H * G, 2) 0, all); fprintf(H*G 0: %d\n, valid);参数说明这里mod(..., 2)是必需的因为 H 和 G 的元素都是 0/1普通矩阵乘法会得到非零整数无法直接判断正交性。如果valid为 0优先怀疑makeParityChk里的高斯消元是否对列做了合理的置换而不是怀疑随机种子。参数含义典型值调参影响n码长96, 256, 1024码长增大BP 译码性能提升但仿真耗时线性增长j列重3, 4列重越小高信噪比平台越低但短环概率增大k行重6, 8决定码率近似 1 - j/kseed随机种子任意整数影响 H 矩阵的短环分布需固定才能复现仿真2.3 短环检测不查环的 LDPC 仿真没有意义这是最容易忽略的一步。随机构造的 H 矩阵中 4 环概率很高而 BP 译码在 4 环存在时外部信息会在两个变量节点之间来回传递导致误码率平台抬高。检查方法如下% 基于 Tanner 图检查 4 环 G_tanner H * H; % 如果存在 4 环非对角线上会出现大于 1 的元素 four_cycle_flag any(G_tanner(:) 1);注意这里的逻辑H 的任意两列如果同时在两行上都有 1就构成一个 4 环反映到 H * H 上就是行列重叠度大于 1。如果检测到 4 环建议重新调用makeLdpc换一个随机种子或者改用渐进边增长算法。但需要注意的是makeLdpc不一定自动做了环长优化所以这一检查必须由使用者自己完成我通常将它放到仿真脚本的第一段。3. 硬判决译码的比特翻转实现从校验子到迭代翻转硬判决译码的数学基础是校验子方程 s r * H其中 r 是接收序列硬判决后的 0/1 向量。s 的每个分量对应一个校验方程是否被破坏。比特翻转Bit Flipping, BF算法的核心思想是找出让不满足校验方程最多的那个比特把它翻转然后迭代重复。3.1 标准比特翻转的单次迭代流程decodeBitFlip.m里最朴素的版本流程如下计算校验子统计每个变量节点参与的失效校验方程数量翻转失败次数最大的那个比特更新校验子进入下一轮。代码层面的核心循环长这样function decoded decodeBitFlip(rx, H, max_iter) [M, N] size(H); r double(rx(:) 0.5); % 硬判决 for iter 1:max_iter s mod(r * H, 2); % 计算校验子 if all(s 0) break; % 所有校验方程满足 end % 统计每个变量节点涉及的失效校验数 fail_count sum((s * ones(1, N)) .* H, 1); [~, idx] max(fail_count); % 选择失效数最多的比特 r(idx) 1 - r(idx); % 翻转 end decoded r; end需要重点说明的是fail_count的计算方式。s * ones(1, N)把校验子向量 s 的每个分量复制到所有变量节点维度上再点乘 H这样每一列得到的就是该变量节点参与了多少个当前失效的校验方程。选max对应的变量节点翻转是最原始的贪心策略后续可以改成带权重的加权比特翻转Weighted BF用信道接收值的幅度作为权重这样可靠性低的比特更容易被翻。3.2 为什么硬判决在加性白高斯噪声信道下性能不佳硬判决在 BPSK AWGN 信道下的问题在于它把接收信号在 0 处一分为二彻底丢弃了软信息。一个幅度为 0.01 的接收值和幅度为 1.2 的接收值在硬判决后都被判定为同一符号但前者明显更不确定。这个不确定性正是软判决译码的增益来源。比特翻转的优势在于复杂度。每次迭代只涉及 0/1 矩阵乘法和整数计数不涉及浮点对数或指数运算在 FPGA 或嵌入式场景里非常友好。它的劣势在于收敛性依赖翻转门限的选取固定max策略在低信噪比区域容易出现振荡某个比特被翻转后另一个比特又变成最大失败数导致前后两次迭代把同一个码字来回改动。3.3 加权比特翻转的改进思路如果decodeBitFlip.m只实现了标准 BF我一般建议改成加权版本改动很小把信道接收值的绝对值当作可靠性权重与失效校验数相乘。比如w abs(rx); % 可靠性权重 score fail_count .* w; % 加权得分 [~, idx] max(score);这里w越大表示接收越可靠理论上得分应该更低才对。实际操作中需要做一次反向映射比如w 1 ./ (1 abs(rx))让不可靠比特得分更高。加权 BF 在高信噪比区间比标准 BF 有明显增益且不会显著增加复杂度。4. 软判决译码的对数域实现从概率域到对数似然比的工程选择概率域 BP 译码的公式推导很多教材都有但直接实现往往会遇到数值下溢——置信度传播里的概率乘积会随迭代指数级缩小在 MATLAB 的 double 精度下大约 20 次迭代后概率值就变成 0 了。decodeProbDomain.m存在的主要意义是教学演示而真正稳定高效的是decodeLogDomain.m它在对数域里把乘法转成加法把归一化消掉。4.1 概率域 BP 的更新公式回顾概率域 BP 里变量节点传给校验节点的消息是置信度比值校验节点更新需要计算 P(b_i 0) / P(b_i 1) 的乘积。其数值不稳定的根本原因是这些值在 0 到 1 之间连乘随迭代轮数增加指数级趋近 0。decodeProbDomain.m的典型结构如下% 概率域消息初始化 Llr 2 * rx / sigma^2; % 将接收值转为对数似然比 p0 1 ./ (1 exp(Llr)); % 转换成概率 p1 1 - p0; % 迭代更新校验节点消息更新需要重置为恒等元 for iter 1:max_iter % 每个校验节点重新计算发给变量节点的可靠性 for m 1:M idx find(H(m, :)); for n idx prod0 prod(p0(idx(idx ~ n))); prod1 prod(p1(idx(idx ~ n))); r0(m, n) prod0 / (prod0 prod1); r1(m, n) prod1 / (prod0 prod1); end end % 变量节点更新后归一化 end这段代码有三个问题值得指出。其一prod(...)对大量接近 0 的小数连乘结果会低于realmin直接变成 0其二循环嵌套导致矩阵规模稍大就慢得不可接受其三prod0 / (prod0 prod1)在两者都为 0 时产生 NaN。这些代码层面的现象解释了为什么工程实现几乎都转投对数域。4.2 对数域最小和算法decodeLogDomain 的工程化简对数域 BP 的关键是把校验节点更新中的双曲正切乘积转化为加法L(r_mn) 2 * atanh( ∏ tanh(L(q_mn)/2) )工程上继续做近似得到最小和Min-Sum形式L(r_mn) ≈ ∏ sign(L(q_mn)) * min |L(q_mn)|这就是decodeLogDomain.m和decodeLogDomainSimple.m之间的主要分野。decodeLogDomain实现了带校正因子的完整对数域 BP而decodeLogDomainSimple直接走 Min-Sum 近似省去了tanh和atanh的浮点计算。Min-Sum 的误码率性能会比完整 BP 差约 0.2~0.4 dB但计算量大幅降低。核心迭代代码可以简化成下面这个形式这也是我重构decodeLogDomainSimple时最常用的骨架for iter 1:max_iter % 校验节点更新 for m 1:M idx find(H(m, :)); for n idx others idx(idx ~ n); signs prod(sign(Lvar(others, m))); % 符号相乘 mini min(abs(Lvar(others, m))); % 取最小绝对值 Lcheck(m, n) signs * mini; % 最小和近似 end end % 变量节点更新 Lvar Llr sum(Lcheck, 2); % 累加所有校验消息 end在这个实现里Lvar是变量节点消息矩阵Llr是信道初始对数似然比。需要注意Lvar(others, m)的索引粒度它取的是第 m 个校验节点连接的所有变量节点中除当前 n 以外的消息集合。signs * mini是 Min-Sum 的全部秘密不需要计算双曲正切。算法外部信息计算复杂度性能损耗概率域 BP概率连乘 归一化O(N·j·k) 浮点乘除数值不稳定易下溢完整对数域 BPtanh/atanh高涉及双曲函数基准无近似Min-Sum符号 最小值极低适合硬件约 0.3 dB归一化 Min-Sum×0.75低只多一次乘法约 0.1~0.2 dB4.3 从概率域到对数域的代码等价转换一个容易被忽略的事实是概率域的p0/(p0p1)归一化和对数域的LLR累积是数学等价的只是换了个代数空间。转换公式是L log(p0/p1)如果你手上有decodeProbDomain.m的代码想验证它与decodeLogDomain.m的数值关系可以直接在两者末尾把软输出转回硬判决后比对码字。差异来源几乎只有两个概率域因下溢导致迭代发散或者两种实现使用了不同的最大迭代次数。5. 软判决译码在 ldpcBER.m 中的性能仿真与迭代参数调优ldpcBER.m是整个资源里最能直接看到效果的文件。它的功能是蒙特卡洛仿真画出不同信噪比下的误码率曲线。这类仿真脚本的框架不复杂但有几个关键细节会影响仿真结果的可靠性和速度。5.1 蒙特卡洛仿真的主循环结构一个标准的性能仿真循环包括数据生成、信道加噪、译码、错误统计四步。资源中的ldpcBER.m把这几个步骤组织如下for snr_idx 1:length(snr_list) num_err 0; num_bit 0; while num_err 100 total_frame max_frame % 1. 随机生成信息位并编码 info randi([0 1], 1, K); codeword mod(info * G, 2); % 2. BPSK 调制 AWGN 信道 x 1 - 2 * codeword; sigma2 1 / (2 * 10^(snr/10) * rate); rx x sqrt(sigma2) * randn(size(x)); % 3. 软判决译码 decoded decodeLogDomain(rx, H, max_iter); % 4. 统计错误比特 err_count sum(decoded ~ codeword); num_err num_err err_count; num_bit num_bit length(codeword); end ber(snr_idx) num_err / num_bit; end这个循环里最容易出错的是噪声方差的换算。BPSK 调制下若信息比特经过码率为rate K/N的编码则输入信噪比对应到符号级噪声方差要乘上码率。sigma2 1 / (2 * 10^(snr/10) * rate)的含义是snr是 Eb/N0每比特能量与噪声功率谱密度之比乘以码率后得到 Es/N0再换算成符号级噪声方差。如果这里忘记乘rate仿真的误码率曲线会比真实情况偏移码率相关的 dB 数。5.2 收敛阈值的选取准则decodeLogDomain.m里的最大迭代次数是一个需要反复权衡的参数。一般经验是在低信噪比区域增加迭代轮数对性能帮助明显在高信噪比区域20 到 30 轮之后误码率曲线逐渐平台化继续增加收益有限。我从这份资源的仿真脚本里验证了一个有用的规律信噪比区间推荐最大迭代原因1~2 dB50~100低信噪比下收敛慢需要更多轮数3~5 dB20~30高信噪比下多数码字 10 轮内收敛5 dB 以上10~15继续迭代容易进入错误平层耗时无收益判断是否达到错误平层的方法是观察最后一个信噪比点的误码率是否随着迭代次数增加而显著下降。如果 50 轮和 100 轮的结果几乎一致说明已经触及了该码字集合的译码极限。5.3 MATLAB 仿真提速的三个实际手段这类蒙特卡洛仿真在 MATLAB 里最大的痛点是速度。我常用的优化手段有三个按收益排序第一向量化校验节点更新。decodeLogDomain.m里如果使用了for m 1:M逐行更新改成对 H 矩阵的非零元整体操作速度能提升 5 到 10 倍。做法是把 H 的稀疏结构提取出来用find预先计算各校验节点的连接关系避免每次迭代都扫一遍 H。第二提前终止机制。每轮迭代后检查校验子是否全零如果是就跳出一轮迭代并更新统计。这个逻辑在decodeBitFlip.m和decodeLogDomain.m里都有体现但如果你自己写过容易忽略「校验子全零不代表译码正确」这一点因为 H 存在多个码字对应同一校验子的情况。第三减少不必要的软输出转换。如果只需要最终硬判决结果不必每轮迭代都计算s mod(r * H, 2)只在校验子更新后才需要用到。但要注意提前终止依赖校验子不能省。在这些提速手段之外还有一个容易踩的坑randi的随机数生成器在不同 MATLAB 版本下默认流不一致如果你需要精确复现别人的误码率曲线记得在脚本开头固定随机种子rng(2024); % 固定随机数生成器这会直接决定蒙特卡洛仿真的可重复性。导出图片时我通常用exportgraphics(gcf, ber_curve.png, Resolution, 300)避免手动缩放导致图例文字糊成一片。如果对比多条曲线记得把不同算法的Marker区分开黑白打印时不至于无法辨认。本文还有配套的精品资源点击获取