多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

LDPC Decoding in MATLAB: From Belief Propagation to Min-Sum Optimization

LDPC Decoding in MATLAB: From Belief Propagation to Min-Sum Optimization 简介面向通信算法、FPGA开发与信道编码方向的研究人员和工程师提供一套基于MATLAB的LDPC译码算法实现。资源以码长960、码率0.5为参数覆盖校验矩阵构造、置信传播迭代译码、定点量化模拟等完整流程并按FPGA硬件实现的语言与时序特点进行仿真。压缩包共16个文件8个源码文件为核心涵盖主程序、节点更新、最小值查找、阈值控制等模块另有4个自动备份文件、2个文本说明文件和2个数据文件用于存放校验矩阵与中间变量整体仅21KB结构清晰。当前已有488人次查看学习。通过运行和修改代码可理解消息传递迭代机制、最小和近似处理以及浮点模型转定点模型时精度与资源之间的折中同时获得一套可直接参考的硬件化译码模型适合作为FPGA实现前的算法验证基础也可用于通信原理课程、研究生课题或毕业设计中的LDPC仿真与性能分析。1. LDPC 译码算法在 Matlab 仿真里比编码更值得花时间校验矩阵给定以后LDPC 编码就是一次稀疏线性变换ldpcEncode一行就能跑完真正决定误码率曲线形状、决定一次仿真要等多久的是译码器。5G NR、DVB-S2、Wi-Fi 6 的物理层都选了 LDPC在 Matlab 里做链路级仿真时译码算法选型、迭代次数上限、早停条件、最小和缩放因子每一项都同时影响 BER 曲线和仿真耗时。下面的内容面向要在 Matlab 里搭 LDPC 编译码链路、评估译码性能、做预研验证的工程师从置信传播的原理讲到ldpcDecode的配置再给出手写最小和译码、调参和排错的具体做法。2. 从置信传播到最小和Matlab 里先选对 LDPC 译码算法再写代码2.1 置信传播在译什么两类节点之间传的是外信息LDPC 译码的本质是在 Tanner 图上迭代传递消息。H 矩阵是 M×N 的每一行对应一个校验节点每一列对应一个变量节点编码比特。变量节点初始拿到信道 LLR每一轮迭代把「除了某个校验节点以外所有其他来源汇总的信息」发给那个校验节点。这种刻意排除对方自身信息的做法就是外信息extrinsic information——自己说过的话不能自己再听回去。两类节点的更新式是固定的。变量节点更新是一个求和结构v2c(m,n) L(n) Σ_{m∈N(n)\m} c2v(m,n)校验节点更新是置信传播的核心用 tanh 域表示c2v(m,n) 2·atanh( Π_{n∈N(m)\n} tanh(v2c(m,n)/2) )变量节点的更新可以优化成「总 LLR 减去对方上次发来的消息」每个变量节点只需要维护一个总 LLR这也是所有高效实现的基础。校验节点的问题在于 tanh/atanh 在 ±1 附近很陡直接算在数值上容易出问题常见的做法是换到 φ 对数域φ(x) −ln tanh(x/2)。换域之后乘法变成加法校验更新变成 φ⁻¹(Σφ(|m|))而这一步正是最小和近似的入口。2.2 最小和近似为什么快两个最小值就够φ 函数在正半轴衰减非常快Σφ(|m|) 这个求和里最大的那一项几乎就决定了结果。工程上的做法是直接忽略其余项把校验更新近似成c2v(m,n) ≈ (Π sign) × min_{n∈N(m)\n} |v2c(m,n)|这就是最小和Min-Sum。每个校验节点只需要在邻居消息里找出幅度最小的两个值回传给最小幅度那个变量节点时用第二小值回传给其他变量节点时用最小值。符号是全部符号的乘积再乘上自身符号——因为符号平方等于 1这就等价于刨掉自身的符号乘积。两次线性扫描就能完成一个校验节点的全部更新不再有任何超越函数。要注意最小和是「高估」消息幅度的精确值是 φ⁻¹(Σφ(|m|))它一定不大于 min|m|。消息幅度被高估意味着译码器过度自信表现就是高信噪比段出现错误平台。修正手段有两类归一化最小和把幅度乘一个 α常见 0.7~0.85偏移最小和把幅度减一个 β常见 0.3~0.6再截断到 0。四种算法的对比关系如下算法校验节点更新复杂度性能特点BP / 和积tanh 域精确计算最高基准瀑布区最优最小和符号 × 最小值低比 BP 差 0.2~0.5 dB高 SNR 有平台归一化最小和α × 最小值低α 标定后接近 BP硬件主流偏移最小和max(m−β, 0)实际仿真里我一般先用最小和把链路跑通确认逻辑没问题再决定要不要上修正项。对多数规则码归一化最小和跟 BP 的差距能压到 0.1 dB 左右换来的是接近一个数量级的加速。2.3 Matlab 里三条路线怎么选ldpcDecode、System object 与手写Matlab 里做 LDPC 译码有三条常见路线选哪条取决于你要标准码仿真、统计迭代次数还是改算法本身。R2021b 之后的新接口ldpcEncode/ldpcDecoderConfig/ldpcDecode配置对象风格直接支持 belief-propagation 和 min-sum支持伴随式早停处理 DVB-S2、5G NR 这种标准码最省事。一般做 BER 仿真默认走这条。老一代comm.LDPCEncoder/comm.LDPCDecoderSystem object旧代码兼容性好step的第二个返回值是每帧实际迭代次数想统计平均迭代次数时它比新接口方便。手写消息传递只有当你需要改调度方式分层 layered、shuffled、设计新的修正项、或者验证一个想法时才需要。我一般只在小矩阵上手写大码长直接用内置函数。三条路线不冲突用内置函数跑性能基准用手写版本验证算法想法两边共用同一份 H 和同一份 LLR。第 3 章先把内置接口这条主线讲透。3. 用 ldpcDecode 搭一条可运行的 LDPC 编译码链路3.1 校验矩阵从哪来dvbs2ldpc 优先中小矩阵用 ldpcQuasiCyclicMatrix译码器要工作第一件事是拿到校验矩阵 H。标准码直接调dvbs2ldpc它按 DVB-S.2 标准生成 64800 码长的稀疏校验矩阵H dvbs2ldpc(5/6); % DVB-S.2 高码率校验矩阵 [nParity, nBits] size(H); % nParity10800, nBits64800 kBits nBits - nParity; % 信息位长 54000dvbs2ldpc的入参是码率常见取值如下表注意返回矩阵的尺寸随码率变化很大码率 rH 尺寸 (行×列)信息位长 K适用场景1/232400×6480032400低码率瀑布区靠左3/416200×6480048600中等码率5/610800×6480054000高码率官方示例常用这个矩阵是 64800 码长的标准码单帧编码译码没问题但几千帧的 BER 扫描会明显变慢所以脚本里帧数要控制。想用中小矩阵练手可以自己构造准循环矩阵Z 27; % lifting 因子 B [ 1 0 -1 2 ; ... % base matrix-1 表示全零块 -1 1 3 -1 ; ... 0 -1 -1 1 ]; H ldpcQuasiCyclicMatrix(B, Z); % 81×108 的稀疏矩阵注意ldpcEncoderConfig要求 H 满行秩随机选的 shift 值有可能不满足报 rank 错误时换一组 B 里的值即可。3.2 译码器配置对象算法、迭代次数与早停参数ldpcDecode的输入除了 LLR 还有一个配置对象配置对象在译码之前一次性设定decCfg ldpcDecoderConfig(H); decCfg.Algorithm min-sum; % 或 belief-propagation decCfg.MaximumIterationCount 30; % 迭代预算 decCfg.IterationTerminationCondition parity-check; % 默认 max decCfg.MinSumScalingFactor 0.8; % 归一化最小和的缩放因子各参数的作用和调整方向整理如下配置项可选值默认值调参方向Algorithmbelief-propagation / min-sumbelief-propagation性能优先选 BP速度优先选 min-sumMaximumIterationCount正整数50小则延迟低大则压低错误平台IterationTerminationConditionmax / parity-checkmax伴随式满足即停省迭代MinSumScalingFactor(0, 1]1缩小到 0.7~0.85 修正高估DecisionTypehard / softhardsoft 输出 LLR方便级联外码MinSumScalingFactor只有在Algorithm为 min-sum 时才有意义而且它只修正幅度不碰符号。如果你用的版本里配置对象没有这个属性可以用comm.LDPCDecoder的同名属性或者直接用第 3.4 节的手写函数扫 α。3.3 完整 BER 仿真脚本编码、BPSK、AWGN 到 LLR下面是一段可以直接跑的 BER 仿真骨架编码、调制、加噪、LLR 计算、译码、误码统计都在%% LDPC 编译码 BER 仿真ldpcEncode ldpcDecode rng(42); % 固定随机种子保证结果可复现 r 5/6; H dvbs2ldpc(r); [nParity, nBits] size(H); kBits nBits - nParity; encCfg ldpcEncoderConfig(H); decCfg ldpcDecoderConfig(H); decCfg.Algorithm min-sum; decCfg.MaximumIterationCount 30; decCfg.IterationTerminationCondition parity-check; nframes 20; % 演示用小帧数跑通后自行加大 EbN0dB 1.5 : 0.5 : 3.0; ber zeros(size(EbN0dB)); for idx 1:numel(EbN0dB) errCnt 0; bitCnt 0; for f 1:nframes info randi([0 1], kBits, 1); cw ldpcEncode(info, encCfg); % 编码64800×1 tx 1 - 2*cw; % BPSK0→11→-1 ebn0 10^(EbN0dB(idx)/10); sigma2 1 / (2 * r * ebn0); % 单位能量实 AWGN 方差 rx tx sqrt(sigma2) * randn(nBits, 1); llr 2 * rx / sigma2; % LLR正数偏向比特 0 dec ldpcDecode(llr, decCfg, OutputType, info); errCnt errCnt sum(dec ~ info); bitCnt bitCnt kBits; end ber(idx) errCnt / bitCnt; fprintf(EbN0%4.2f dB BER%g\n, EbN0dB(idx), ber(idx)); end semilogy(EbN0dB, ber, o-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER);几个关键点说明。tx 1 - 2*cw把 0/1 映射到 1/−1接受符号rx tx n中 n 的方差是 σ²LLR 的推导是ln P(b0|rx)/P(b1|rx)展开后正好是2*rx/sigma2。σ² 由码率折算单位符号能量下 Es/N0 r·Eb/N0实 AWGN 中 N0 2σ²所以 σ² 1/(2·r·ebn0)。LLR 矩阵的约定是每一行对应一个编码比特位置每一列对应一帧列数就是帧数。提示第一次运行先把nframes改成 2、EbN0dB只留一个点确认程序能跑完且 BER 在 1e-2 附近再扩大规模否则调一次参数等半小时。3.4 手写归一化最小和两个循环看懂消息更新内置函数把复杂度都封装了想看懂 min-sum 到底在算什么手写一个泛洪式flooding版本最有帮助function [vhat, nIter] ldpc_minsum_decode(H, LLR, maxIter, alpha) % H: 稀疏校验矩阵 M×NLLR: N×1正数偏向比特 0 % alpha: 归一化缩放因子传 1 就是原始最小和 [M, N] size(H); vn cell(N,1); cn cell(M,1); for m 1:M, cn{m} find(H(m,:)); end for n 1:N, vn{n} find(H(:,n)); end L LLR(:); v2c zeros(M, N); % 变量节点 → 校验节点 c2v zeros(M, N); % 校验节点 → 变量节点 nIter maxIter; for it 1:maxIter % 变量节点更新总信息减去回传的那一条 for n 1:N nb vn{n}; for m nb others nb(nb ~ m); v2c(m,n) L(n) sum(c2v(others, n)); end end % 校验节点更新两个最小幅度 符号外积 for m 1:M nb cn{m}; mag abs(v2c(m, nb)); sgn sign(v2c(m, nb)); [min1, i1] min(mag); tmp mag; tmp(i1) inf; % 找第二小值 min2 min(tmp); for j 1:numel(nb) amp min1; if j i1, amp min2; end c2v(m, nb(j)) alpha * prod(sgn) * sgn(j) * amp; end end % 硬判决 伴随式校验 Ltot L sum(c2v, 1).; vhat double(Ltot 0); if all(mod(H * vhat, 2) 0) nIter it; return; end end end校验节点更新里prod(sgn) * sgn(j)是全部符号积乘自身符号等价于刨掉自身的符号积这就是外信息。幅度部分发给最小幅度那个邻居用min2发给其他邻居用min1。总 LLR 的符号决定硬判决Ltot 0判为比特 1与信道 LLR 的正负约定一致。这个版本用稠密矩阵存消息M×N 个 double几千比特的小码没问题64800 码长的 dvbs2 矩阵会直接撑爆内存只适合教学和小码验证。4. 调优最小和译码迭代次数、缩放因子与早停的取舍4.1 迭代次数曲线译到第几轮就不用再等了迭代次数和 BER 的关系是典型的「先陡后平」前 5 轮把大部分错误纠正掉10 轮以后每加一轮收益很小而每轮迭代的成本是线性的。扫一遍迭代预算就能看到拐点iterList [2 5 10 20 50]; berIt zeros(size(iterList)); for i 1:numel(iterList) decCfg.MaximumIterationCount iterList(i); % 固定一个 EbN0 点选 BER 在 1e-2~1e-3 的位置 % 把 3.3 节的内层循环搬进来统计该点的误码率 end选测试点有个经验太靠近瀑布区左边所有配置都是 0.5看不出差别太靠右全部是 0也看不出差别。要选在 BER 曲线上斜率最大的那段。对大多数规则码min-sum 配 10~20 轮、BP 配 30~50 轮就够超过 50 轮通常只是在压错误平台而那个平台往往不是迭代次数造成的而是消息幅度近似造成的该去调 α。4.2 归一化因子 α 怎么标定一个小循环扫出来α 的作用是抵消最小和对消息幅度的高估。α 太大会残留平台α 太小会让瀑布区右移两者都不行。常见的做法是固定一个中间信噪比点扫一遍 αalphas 0.55:0.05:0.95; for a alphas % 用 3.4 节手写函数在固定 EbN0 点上处理全部帧 for f 1:nframes [vhat, ~] ldpc_minsum_decode(H, llrList(:,f), 30, a); % 统计 vhat 与 infoList(:,f) 的不一致比特数 end end对规则 (3,6) 码α 通常落在 0.75~0.85我一般从 0.8 起步。不同码要重扫列重、行重变了最优值就变了不要跨码套用。α 标定完成后对全部 SNR 点都适用不需要每个点单独调。偏移最小和的 β 同理在 0.3~0.6 之间扫判断标准是瀑布区不明显右移的前提下高 SNR 平台压得最低。4.3 早停的三个注意点parity-check 不是免费的IterationTerminationCondition parity-check的含义是每一轮硬判决后算伴随式 H·vhat全零就停。它有三个容易忽略的坑。第一高 SNR 下收益明显帧基本在前 5 轮就满足校验低 SNR 下伴随式几乎不会全零实际还是跑满MaximumIterationCount别指望它降低低信噪比下的延迟。第二伴随式全零只能说明「译码结果是一个合法码字」不能证明它是对的。特别地如果 LLR 符号整体取反译码器会收敛到合法码字的按位取反——LDPC 是线性码补码仍然是合法码字伴随式照样全零。第三早停之后每次迭代多了一次稀疏矩阵乘法的开销帧长 64800 时这个开销不小帧很短时甚至可能抵消省下的迭代时间。5. 排错与验证用 BER 曲线确认 LDPC 译码器没有暗病5.1 三种典型异常BER 卡 0.5、BER 接近 1 和错误平台BER 曲线的形状能直接暴露译码链路的问题对照关系如下现象最可能的原因检查方法BER 始终在 0.5 附近不随 EbN0 下降帧同步错位信息位切片位置不对无噪声单帧自检比对 dec 与 info高 SNR 下 BER 接近 1LLR 符号约定反了给 llr 取负重跑一次瀑布区正常高 SNR 出现抬尾平台最小和未缩放或迭代不足加 α0.8迭代提到 50BER 接近 1 而不是 0.5是个很有诊断价值的现象符号整体反了之后译码器收敛到补码补码合法且伴随式全零所以错误率逼近 1。看到 BER 在 0.9 以上先别怀疑噪声模型直接检查 LLR 正负约定。5.2 无噪声单帧自检放在所有仿真之前不管改了什么参数先跑一遍无噪声自检它能一次性暴露维度、码率和符号约定三类问题info randi([0 1], kBits, 1); cw ldpcEncode(info, encCfg); llr0 1 - 2*cw; % 无噪声 BPSK 符号幅度不影响结果 dec ldpcDecode(llr0, decCfg, OutputType, whole); assert(all(mod(H * dec, 2) 0)); % 伴随式必须全零 assert(isequal(ldpcDecode(llr0, decCfg), info)); % 信息位必须逐位一致无噪声时 LLR 的幅度不重要符号对就行。这步过了说明 H 的维度、编码器和译码器的配置是自洽的再去加噪声才有意义。5.3 用 BP 做基准评估最小和修正是否到位最后一个小技巧把 BP 当作性能上限和归一化最小和在同一条链路上对比差距能控制在 0.1~0.2 dB 就说明 α 标定到位。cfgBp ldpcDecoderConfig(H); % 默认 belief-propagation cfgNm ldpcDecoderConfig(H); cfgNm.Algorithm min-sum; cfgNm.MinSumScalingFactor 0.8; % 用同一组 llr 帧分别译码统计两条 BER 曲线每轮调参时把rng种子、EbN0 点、帧数、α 一起写进文件名或 .mat同时记录每个 SNR 点的平均迭代次数。性能比较不只看 BER还要看延迟平均迭代次数就是最好的延迟指标。把无噪声自检和固定种子放在仿真脚本最前面后面画出来的每一条 BER 曲线才具备可比性。本文还有配套的精品资源点击获取
返回列表