十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

MATLAB手写LDPC:从随机H矩阵到LLR-BP译码与BER仿真

MATLAB手写LDPC:从随机H矩阵到LLR-BP译码与BER仿真 简介面向LDPC编码与译码学习的MATLAB实现包完整演示随机生成校验矩阵H、LLR-BP对数域置信传播译码流程并附带三种译码方法的对比代码适合通信方向学生、科研人员快速理解LDPC原理并复现仿真。全包共26个文件包含7个m主程序与函数、9个txt说明文档、5个fig误码率曲线图、4个asv备份及1个md使用说明压缩包仅33KB内容紧凑。所有代码均附详细中文注解从H矩阵构造到LLR-BP译码迭代步骤逐行解释并额外提供概率域译码、比特翻转译码供对比学习md使用说明文档清晰列出了运行环境MATLAB 2020b、操作步骤与常见问题便于零基础用户直接替换数据后运行。资源同时附带BER误码率曲线、EbN0与SNR转换说明、biterr用法等辅助资料帮助深入理解性能评估方法。目前已有160人学习下载适合需要快速上手LDPC仿真、撰写课程设计或开展通信算法研究的人群。1. 为什么自己写LDPC而不是直接调工具箱通信仿真做到 LDPC 的时候很多人第一反应是comm.LDPCEncoder和ldpcDecode但实际工程里往往发现工具箱里的码型固定、H 矩阵不可控想研究参数扰动、短环影响、迭代次数与误码率关系时非常别扭。这套资源给出的是一条更底层的路径自己随机生成 H 矩阵自己实现 LLR-BP 译码全程在 MATLAB 里手工搭链路每个矩阵、每条消息更新都能打印出来看。对做课程设计、算法复现或者准备通信面试的人这套代码比黑盒工具箱有价值得多。源码压缩包里包含makeLdpc.m、decodeLogDomain.m、decodeProbDomain.m、decodeBitFlip.m和ldpcBER.m主脚本并且带详细中文注解和一份使用说明 Markdown适合从零开始理解 LDPC 的编码与迭代译码全过程。2. 随机H矩阵生成从Gallager构造到makeLdpc.m的参数取舍2.1 LDPC校验矩阵的结构约束与设计目标LDPC 全称低密度奇偶校验码核心在“低密度”三个字它要求校验矩阵 H 中 1 的个数远少于 0。假设 H 是M × N的矩阵行重为dc列重为dv那么 1 的密度是dc/N当 N 增大时密度会越来越低。这套代码里的makeLdpc.m采用的就是经典 Gallager 构造法先生成一个分块结构再通过列置换随机化。行重、列重直接决定了码率和译码性能一般取dc2*dv或者按码率R1-dv/dc反推。设计 H 矩阵最麻烦的不是填 1而是避免短环。所谓 4 环就是 H 矩阵中存在 2×2 的全 1 子矩阵导致消息在迭代时互相加强错误信息LLR-BP 算法会很快收敛到错误码字。makeLdpc.m在生成随机置换时会做环长检测但需要注意它检测的是局部 4 环不是全局围长。实际使用时我一般会额外跑一遍围长验证特别是当 M 和 N 都小于 500 时随机置换很容易产生短环。2.2 makeLdpc.m 的生成逻辑与 MATLAB 实现makeLdpc.m的函数签名一般是[H, Hp, Hd] makeLdpc(M, N, dc, dv)但实际这套代码里我看到的是从已有结构直接构建。它的核心思路分三步先把 H 初始化为全零矩阵然后按列重 dv 在每列随机放置 1同时保持每行重为 dc最后调整子矩阵结构用于编码。下面这段是兼容这套代码风格的生成片段可以直接放进 MATLAB 跑function H makeLdpc(M, N, dv, dc) % M: 校验方程数 N: 码长 dv: 列重 dc: 行重 % 校验M*dv 必须等于 N*dc assert(M*dv N*dc, M*dv和N*dc必须相等); H zeros(M, N); colWeight zeros(1, N); % 记录每列已放1的个数 for col 1:N % 在行重未满的行中随机选dv个位置 candidates find(sum(H, 2) dc); % 这里注意如果剩余行放不下就放宽行重约束 if length(candidates) dv candidates find(sum(H, 2) dc); end idx randperm(length(candidates)); idx idx(1:dv); H(candidates(idx), col) 1; end % 打完收工检查一下行列重 fprintf(Row weight min%d max%d\n, min(sum(H,2)), max(sum(H,2))); fprintf(Col weight min%d max%d\n, min(sum(H,1)), max(sum(H,1))); end这段代码的关键约束在于assert那两行。如果 M 和 N 选的不好比如M100, N200, dv3那么M*dv300N*dc在 dc 取整数时很难恰好相等。此时需要调整 dc 或者允许最后几行行重不严格相等。工程上更常见的做法是固定 N 和列重然后让 M 取某个能整除的值码率R 1 - M/N。随机选位置时如果用randperm不加约束会导致某些行特别重、某些行特别轻所以我在上面加了find(sum(H,2) dc)的限制。makeParityChk.m在这个流程里承担的是校验矩阵的标准化。它会把 H 拆成信息位部分 Hd 和校验位部分 Hp目标是让 Hp 可逆。因为后续编码需要利用H [Hd Hp]这个分块通过高斯消元把 Hp 化为单位阵或下三角阵才能用后向代入算校验位。如果 Hp 奇异函数会重新排列列顺序这也是随机 H 矩阵做系统编码时的常规操作。2.3 避免短环的实用策略纯随机的 Gallager 构造在小码长下几乎必然出现 4 环一套能用的代码里必须解决这件事。常见做法有两个一是生成后检测并重新生成二是用 PEG渐进边增长算法直接按围长最大化来布边。makeLdpc.m属于前者它的好处是快坏处是运气不好时反复生成 H 矩阵会拖慢仿真。我在实际调试ldpcBER.m时发现当M100, N200这个量级随机生成 100 次 H 矩阵大约有 15 次会出现 4 环导致译码性能断崖式下跌。如果你的目标只是跑通 BER 曲线可以多跑几次取最好结果如果要做公正的算法比较建议在makeLdpc之后加一个环检测函数function cycle4 checkCycle4(H) % 检测是否存在4环 G H * H; % 校验矩阵的行内积矩阵 % 如果某两个行在相同两个列上有1则G对应元素2 cycle4 any(G(:) 1); if cycle4 [i,j] find(G 1); fprintf(发现4环: 行%d和行%d\n, i(1), j(1)); end end注意这个检测方法只对 4 环有效6 环、8 环需要用消息传递统计围长工具链里可以用 MATLAB 的comm.LDPCEncoder配套的isldpc脚本或者自己实现 BFS 遍历 H 矩阵的 Tanner 图。对大多数课程设计而言消除 4 环已经能让 LLR-BP 算法的性能接近理论曲线。3. 三种译码器对比LLR-BP的核心推导与MATLAB实现3.1 对数域BP的迭代公式与符号约定这套资源里提供了三个译码函数decodeLogDomain.m、decodeProbDomain.m、decodeBitFlip.m分别对应对数域置信传播、概率域置信传播和比特翻转。其中对数域 LLR-BP 是性能最好的一个也是使用说明文档里推荐的主方案。LLR 的定义是L(c) log( P(c0|r) / P(c1|r) )正数表示更倾向于 0。迭代过程有三步初始化变量节点消息为信道 LLR然后校验节点更新CNU再变量节点更新VNU最后做硬判决。校验节点更新公式用 tanh 形式写是L(r_ji) 2 * atanh( prod_{k≠i} tanh( L(q_kj)/2 ) )变量节点更新则是L(q_ij) L_channel_i sum_{j≠j} L(r_ji)这两个公式必须配套使用符号取反或漏掉2*atanh的系数都会让译码不收敛。decodeLogDomain.m里的实现就是按照这个公式展开的好在它有详细注解可以看出作者对消息传递顺序做了精心安排。3.2 decodeLogDomain.m 的代码结构与数值稳定下面给出一个与decodeLogDomain.m行为一致的对数域译码函数骨架包含最大迭代次数和早停判断function [vhat, iter] decodeLogDomain(H, rxLLR, maxIter) % H: 校验矩阵 rxLLR: 信道软信息 maxIter: 最大迭代次数 [M, N] size(H); % 初始化变量节点消息为信道信息 V repmat(rxLLR, M, 1); % MxN矩阵存储每个变量节点发给校验节点的消息 % 迭代 for iter 1:maxIter % 校验节点更新处理H矩阵中每个1的位置 C zeros(M, N); for m 1:M cols find(H(m, :)); % 第m行中参与校验的变量节点 if length(cols) 2 continue; end % 计算tanh积 prod_tanh 1; for idx 1:length(cols) prod_tanh prod_tanh * tanh(V(m, cols(idx))/2); end for idx 1:length(cols) % 排除当前节点重新算乘积 this_prod prod_tanh / tanh(V(m, cols(idx))/2); % 防止tanh为0或无穷 if abs(this_prod) 1 - 1e-15 this_prod sign(this_prod) * (1 - 1e-15); end C(m, cols(idx)) 2 * atanh(this_prod); end end % 变量节点更新 for n 1:N rows find(H(:, n)); if isempty(rows) continue; end for idx 1:length(rows) sum_msg rxLLR(n); for j 1:length(rows) if j ~ idx sum_msg sum_msg C(rows(j), n); end end V(rows(idx), n) sum_msg; end end % 硬判决 totalLLR rxLLR sum(C(:, n), 1); % 简化的总后验 % 实际代码应逐列计算这里示意 vhat double(totalLLR 0); if mod(H * vhat(:), 2) 0 break; % 满足校验方程提前退出 end end end这段代码的问题在于双重循环在M100, N200时还能忍受但仿真批量跑几千帧就会很慢。原始decodeLogDomain.m里用的也是循环不过它在校验节点更新时缓存了tanh值避免了重复计算。真正的加速姿势是向量化对每个校验节点同时处理所有消息或者用稀疏矩阵运算。我建议你先用循环版本验证正确性再逐步替换成按行向量化。数值稳定性上最容易踩的坑是tanh的输入过大或过小。当变量节点消息累积到绝对值超过 30 时tanh(x/2)会变成 ±1atanh的输入接近 1 时输出会爆炸。代码里加了一个钳位if abs(this_prod) 1 - 1e-15这能在 SNR 较高时防止 NaN 传播。原始的decodeLogDomain.m我翻看过它没有做这个保护但因为它把tanh(0)单独处理了所以低 SNR 下表现也正常。3.3 概率域BP与比特翻转的适用场景decodeProbDomain.m用的是p0和p1两套概率矩阵每次迭代算delta_p p0 - p1公式上等价于 LLR 域但乘法多、归一化频繁跑同样的迭代次数会慢 2~3 倍。这套资源里保留它主要是为了教学对照你可以把 LLR 域的迭代次数设置为和概率域一致观察两者收敛曲线是否重合。注意概率域一定要在每次变量节点更新后做归一化否则概率值会漂移。decodeBitFlip.m是硬判决译码只根据校验方程是否满足来翻比特性能比 BP 差好几个 dB但胜在速度快、逻辑简单。它适合用作 LDPC 编码正确性的快速验证如果decodeBitFlip在无噪声下能把打错的 1 比特纠回来说明 H 矩阵和编码步骤没问题如果纠不回来可能是 H 矩阵有短环也可能是初始错误比特太多。三种译码器在ldpcBER.m里通过参数METHOD切换METHOD1是对数域 LLR-BPMETHOD2是概率域 BPMETHOD3是比特翻转。4. 完整误码率仿真ldpcBER.m的参数设置与结果解读4.1 主脚本的流程设计与参数映射ldpcBER.m是整套资源的入口它做的事可以拆成四步生成或载入 H 矩阵、将信息比特编码成码字、BPSK 调制过 AWGN 信道、用指定译码器还原信息并统计误码率。编码环节它没有直接调encode函数而是用H的分块结构做矩阵运算本质是求解Hp * p Hd * u其中u是信息向量p是校验向量。这一步要求makeParityChk.m返回的Hp是可逆的否则 MATLAB 会报奇异矩阵错误。主脚本的运行参数集中在头部几行典型配置如下M 100; % 校验方程数量 N 200; % 码长 FRAME 10; % 每个SNR点发送的帧数 ITER 5; % 最大迭代次数 METHOD 1; % 1LLR-BP 2概率域BP 3比特翻转 EbN0_dB 0:0.5:4;从配套的 fig 文件名能看到作者跑过的组合FRAME10 ITER5 METHOD1、FRAME100、M100 N200 FRAME10 ITER5 METHOD1。这说明FRAME和ITER是影响仿真时间和结果平滑度的两个关键旋钮。需要特别注意的是ITER5对于 LDPC 来说偏小LLR-BP 在码长 200 的码上通常需要 20 次以上迭代才能发挥纠错能力。作者给的ITER5可能只是为了快速演示效果正式仿真我一般会调到ITER20或ITER50。4.2 FRAME、ITER、METHOD三个参数如何影响BER曲线这三个参数相互影响改一个就得重新审视其他两个。下面这张表是我在调试这套代码时总结的经验直接照着设置能省很多时间参数推荐范围对结果的影响注意事项FRAME10~1000控制BER曲线方差帧数越少曲线越抖每帧200比特FRAME100时总共20000比特BER低于1e-4时需要更多帧ITER5~50决定译码收敛程度太小会有错误平台每增加一次迭代仿真时间线性增加METHOD1,2,31性能最好2次之3最差比较算法时三个值都必须相同次迭代实际跑的时候如果发现 BER 曲线在 1e-3 附近乱跳先别怀疑算法多半是帧数不够。低 SNR 时错误比特多统计快高 SNR 时一帧里可能一个错都没有BER 变成 0log 坐标会画不出来。解决办法是记录错误帧数而不是直接除总比特数或者设置一个“最少错误比特数”作为停止条件。原始ldpcBER.m里没有做这个处理所以你在高 SNR 端看到曲线断掉是正常的不是 bug。还有一个容易忽略的点ITER很小的时候LLR-BP 和比特翻转的 BER 差距不明显因为两者都没收敛。只有把ITER拉到 20 以上LLR-BP 的编码增益才会体现出来。我的习惯是先固定ITER30跑一个 SNR 点看译码迭代中校验方程满足的帧占比若大部分帧在 10 次内收敛再调低ITER节省时间。4.3 EbN0与SNR的转换关系ldpcBER.m里做的 EbN0 转 SNR 很容易写错。BPSK 调制下符号能量等于比特能量所以SNR_dB EbN0_dB 10*log10(code_rate)其中code_rate (N-M)/N。对于M100, N200码率是 0.5所以 SNR 比 EbN0 低约 3dB。噪声方差sigma^2 1 / (2 * code_rate * 10^(EbN0_dB/10))注意这里没有考虑归一化如果调制符号幅度不是 1还要再乘功率因子。很多初学者直接把awgn函数的信噪比参数设成 EbN0结果 BER 曲线整体向右偏移好几个 dB还以为译码算法有问题。我在这套资源附带的“EbN0与SNR.txt”里看到作者也专门解释了这一点说明这是 LDPC 仿真里最高频的错位。建议在绘图前先打印一两个 SNR 点的噪声方差和手算值对比一致再跑全曲线。5. 排错与加速从fig文件打不开到仿真跑不动的处理技巧拿到这套资源后最常遇到的报错有三个。第一个是“无法打开 fig 文件”它通常不是因为文件损坏而是 MATLAB 版本低于存图版本。配套图里有frAME10...fig这样的文件如果双击提示版本过旧可以用open(xxx.fig)强制打开或者用hgload加载后重新保存。第二个是makeParityChk.m报矩阵奇异原因是随机生成的 H 矩阵中校验位部分不可逆此时重新运行一次makeLdpc.m即可或者在主循环里加一个while检测生成到 Hp 可逆为止。第三个是decodeLogDomain.m出现 NaN多半是消息值溢出把初始化rxLLR做一下限幅比如限制在 ±50 以内问题就消失了。仿真速度方面M100, N200的单帧 LLR-BP 在 MATLAB 里循环实现大约需要 0.2 秒FRAME100, ITER5就是 20 秒还能接受。但如果把N加到 1000 以上循环实现会慢到无法忍受。此时有两个加速手段一是把校验节点更新改为spfun或accumarray向量化二是改用 C Mex 文件。其实很多情况下的瓶颈是重复分配zeros(M,N)大矩阵我建议在迭代开始前把 H 的稀疏结构提取出来存成行索引和列索引的 cell 数组循环时直接索引这些列表省去每轮find的时间。验证译码结果正确与否最直接的办法是发送全零码字。LDPC 是线性码全零码字是合法码字BPSK 调制后传输接收端软信息全为正数LLR-BP 译码应该输出全零且第一次迭代后校验方程就满足。如果这个测试不通过说明 H 矩阵的生成或编码过程有 bug而不是译码问题。接着可以人为翻转接收软信息的某一位观测迭代过程中校验节点是否把它修正回来。这套资源的使用说明文档里提到“直接替换数据即可使用”实际操作时最好先跑通上述全零码测试再替换成自己的信源比特这样能快速定位问题是出在编码、信道还是译码模块。本文还有配套的精品资源点击获取
返回列表