多进制LDPC的BP解码与调优实战)
简介一套基于BP信念传播算法的多进制LDPC码MATLAB仿真资源专注卫星导航场景下的纠错编译码研究适合通信工程、导航制导与控制专业学生及科研人员使用可帮助快速理解多进制LDPC码的Tanner图构造、BP消息传递机制以及在复杂信道下的完整译码实现。资源包共包含四个文件一个m脚本作为主译码程序实现基于信念传播的迭代译码三个txt文本文件分别保存该多进制LDPC码的指数信息、输入数据与码字元素三者配套驱动译码器运行压缩包整体仅4KB轻量易用便于部署到各类仿真环境。目前已有856人学习/下载特别适合课程设计、毕业设计或工程项目预研阶段使用。多进制LDPC码相比二进制LDPC码可获得更高编码增益在卫星导航这类易受电离层闪烁、多径衰落等影响的链路中能显著提升抗干扰能力、降低定位误差。运行代码可直接得到不同信噪比下的误码率性能曲线还可调整BP迭代次数与消息传递参数观察译码收敛速度及误码率变化从而优化译码策略为实际多进制LDPC编译码方案设计与参数选择提供可靠的参考依据。1. 多进制LDPC的BP解码为什么值得在MATLAB里重写一遍搜索“matlab_ldpc64_BP.rar”这类名字的人大多数想找一个能跑的二进制LDPC BP参考再改成多进制但多进制LDPC的BP和二进制BP不是同一个东西。变量节点消息从0/1的LLR变成GF(q)上的q维概率向量校验节点要做域上卷积代码里多出来的不只是循环而是整套消息模型。更反直觉的一点是在QPSK这类2比特调制下GF(4)符号可以直接映射到星座点短帧性能反而比同码率二进制LDPC更容易逼近香农界。这里用GF(4)作为最小例子拆开BP的变量更新、校验更新和参数调优给出一段能跑的MATLAB代码并说明导航LDPC场景怎么验收它。2. 多进制LDPC为什么要把BP消息从标量改成向量2.1 校验矩阵从GF(2)搬到GF(q)后码字和伴随式都变了二进制LDPC的校验矩阵是0/1矩阵码字满足H c^T0所有加法都在GF(2)上做。多进制LDPC先把符号域换成GF(2^m)也就是把原来的一位比特变成一个m比特符号。GF(4)是其中最小的例子符号集{0,1,2,3}加法为按位异或乘法按本原多项式x^2x1构造。这样一行校验方程不再是几个bit的异或而是GF(4)上系数与符号的线性组合。伴随式是否正确也不能只看奇偶性必须按GF(4)乘法表计算。换来的是码字符号与高阶调制星座的天然对应一个GF(4)符号带2bit正好放进QPSK的一个星座点GF(8)符号可匹配8PSKGF(16)可匹配16QAM。在导航这类信号帧短、功率受限的链路上短帧多进制码比二进制码更容易在同等复杂度下逼近香农限。真正变化的是约束的结构。二进制校验只回答“这几位加起来是不是偶数”GF(4)校验回答的是“这几位按域上的加权和是不是0”。生成矩阵和校验矩阵都不是0/1矩阵后编码也要相应换成GF(4)矩阵乘法解码端判断停止迭代时更不能用二进制的模2伴随式。这个区别在调试时最容易暴露用二进制LDPC的伴随式函数去验GF(4)码字几乎每一帧都会被误判为错帧。2.2 BP消息变成向量先验、V2C、C2V三种消息全是q维分布二进制BP里一条边缘消息是一个标量通常做成LLR。多进制BP里每个变量节点有q种可能取值所以任何经过这条边的消息都是一个长度为q的概率向量。记P0(x)P(c_ix)为信道先验变量节点发给校验节点的V2C消息为M_{v→c}(x) ∝ P0(x) * ∏_{c∈N(v){c}} M_{c→v}(x)校验节点发给变量节点的C2V消息则是在满足该校验方程的前提下对其他关联变量做GF(4)符号的“异或卷积”。一个度数为d_c的校验节点如果按穷举组合计算复杂度是O(d_c q^2)。GF(4)下q4穷举非常快GF(64)、GF(256)再这么做就会卡成幻灯片那时才需要FFT-QSPA或扩展最小和。导航上的LDPC很多是短帧GF(4)或GF(8)已经够用所以先不用急着上快速算法。消息从标量变成向量后归一化也变成向量归一化。每一轮变量更新结束必须把每条消息除以它自己的sum否则几十轮迭代后概率会下溢成NaN。这一点和二进制BP的LLR域做法很不一样很多人把二进制代码搬过来时漏掉的往往是这种看起来不起眼的归一化。2.3 先建好GF(4)的加法表和乘法表整个BP才有地基GF(4)加法等价于按位异或所以可以直接用MATLAB的bitxor生成。乘法表对应本原多项式x^2x1符号2是α符号3是α1。下面的代码先建表再用一个全1校验向量做自检。q 4; gf4_add zeros(q, q); for a 0:q-1 for b 0:q-1 gf4_add(a1, b1) bitxor(a, b) 1; % 存1索引方便后面查表 end end gf4_mul [0 0 0 0; 0 1 2 3; 0 2 3 1; 0 3 1 2]; % 自检全1向量 [1 2 3] 的GF(4)和应为0 acc 0; for c [1 2 3] acc bitxor(acc, c); end assert(acc 0);注意gf4_add存的是1索引a2、b3时bitxor(2,3)1存进矩阵的是2也就是符号1所对应的索引。gf4_conv这类函数直接拿这个表做输出索引能少一层“符号索引互换”的if。gf4_mul表则按0索引存符号查表时用gf4_mul(a1,b1)。两种索引不要混用这是多进制LDPC排错里最常见的暗坑。GF(4)乘法表如下行和列都是符号0到3乘012300000101232023130312这个表里223、332是特征2域的典型性质。写代码时如果这三个“乘法对角线”掉错了后面所有校验节点更新全会错而且不会报异常只是FER永远停在0.5附近。3. 在MATLAB里搭一个GF(4) LDPC BP解码器3.1 用三维数组存V2C和C2V避免cell的隐式开销这里采用m行n列乘q层的三维数组V2C(r, col, s)表示第r个校验节点到第n列变量节点之间、关于符号s的消息概率C2V方向相反。三维数组的好处是切片方便变量节点更新时可以直接用squeeze(C2V(r2, col, :))取出一行概率再用点乘做外信息乘积。虽然用cell数组也能写但每轮迭代都要cell2mat代码更长且速度更慢。初始化的原则是H为0的地方不参与任何更新所以V2C和C2V那里即使有数值也不会被读到。V2C只在find(H(:,col) ~ 0)的行填入P0C2V初始化为全1含义是“校验节点还未提供任何约束信息”。第一次变量更新时每条消息都等于信道先验P0和推导一致。3.2 完整的GF(4) BP主循环和GF(4)异或卷积下面这段代码能直接保存成gf4_bp_decode.m运行。它把变量节点更新、校验节点更新和伴随式检查放在同一个循环里。为了突出BP骨架这段代码的H暂时只处理非零元素为1的情况真实导航工程里的H会出现2、3通用化方法在下一章给。function [hard, iter] gf4_bp_decode(P0, H, max_iter, damp) % P0: n-by-4 信道先验概率P0(i,x1) P(c_i x) % H: m-by-n GF(4)校验矩阵本例只处理非零元素为1的情况 % max_iter: 最大迭代次数 % damp: 消息阻尼0~1建议取0.5 % hard: 译码符号0~3 q 4; [m, n] size(H); V2C zeros(m, n, q); C2V ones(m, n, q); gf4_add zeros(q, q); for a 0:q-1 for b 0:q-1 gf4_add(a1, b1) bitxor(a, b) 1; end end % 变量到校验消息初始化 for col 1:n rows find(H(:, col) ~ 0); for r rows V2C(r, col, :) P0(col, :); end end for it 1:max_iter % 变量节点更新信道先验乘以外来C2V消息 for col 1:n rows find(H(:, col) ~ 0); for r rows out P0(col, :); for r2 rows if r2 ~ r out out .* squeeze(C2V(r2, col, :)); end end out out / sum(out); % 阻尼旧消息占一部分新消息占一部分 V2C(r, col, :) (1 - damp) * squeeze(V2C(r, col, :)) damp * out; end end % 校验节点更新GF(4)上的异或卷积 for r 1:m cols find(H(r, :) ~ 0); if length(cols) 2 continue; end MsgIn zeros(length(cols), q); for idx 1:length(cols) MsgIn(idx, :) squeeze(V2C(r, cols(idx), :)); end MsgOut zeros(size(MsgIn)); for p 1:length(cols) others setdiff(1:length(cols), p); dist [1, 0, 0, 0]; % 单位元 for idx others dist gf4_conv_cyc(dist, MsgIn(idx, :), gf4_add); end MsgOut(p, :) dist; end for idx 1:length(cols) C2V(r, cols(idx), :) MsgOut(idx, :); end end % 后验概率与硬判决 post P0; for col 1:n rows find(H(:, col) ~ 0); for r rows post(col, :) post(col, :) .* squeeze(C2V(r, col, :)); end post(col, :) post(col, :) / sum(post(col, :)); end [~, hard] max(post, [], 2); hard hard - 1; % 索引1..4变成符号0..3 % 伴随式检查注意本例H非零元素全为1 synd_ok true; for r 1:m cols find(H(r, :) ~ 0); acc 0; for idx 1:length(cols) acc bitxor(acc, hard(cols(idx))); end if acc ~ 0 synd_ok false; break; end end if synd_ok iter it; return; end end iter max_iter; end function z gf4_conv_cyc(x, y, gf4_add) z zeros(1, length(x)); for a 1:length(x) for b 1:length(y) z(gf4_add(a, b)) z(gf4_add(a, b)) x(a) * y(b); end end end主函数里变量节点更新先取P0再乘上所有邻居C2V消息最后做一次归一化。这一步的物理含义是消息不能包含自己发出的信息所以要把被更新的那个邻居排除在外。阻尼系数damp在这里不是必需但对概率域BP很有用因为高信噪比时消息会在两个状态之间震荡阻尼能压低振荡。damp越大新消息占的比例越小迭代收敛越慢所以要配合max_iter一起调。校验节点更新用gf4_conv_cyc做GF(4)异或卷积。dist[1,0,0,0]相当于GF(4)离散卷积的单位元它本身不代表任何信道信息只让循环写法统一。对每一个目标边p它取其他边的消息做卷积得到的就是“在其他变量取值组合下当前符号必须取某个值”的概率分布。由于卷积保持概率和为1校验输出不需要再归一化如果你在工程里改了域大小建议仍加一次归一化避免浮点漂移。这个函数能直接跑通3×6这类小矩阵。测试时可以把P0设置为某个已知码字的one-hot概率硬判决应该一次就通过伴随式并且iter返回1。如果连无噪声回环都过不去问题多半在gf4_add表或伴随式检查而不是在迭代更新上。4. 多进制LDPC BP的必调参数与四个排错入口4.1 max_iter、damp、信道先验是三个旋钮max_iter的典型范围是5到20。GF(4)消息是4维向量信息在校验节点之间传播的速度比二进制标量慢迭代太少会让FER曲线出现平台加太多又只增加耗时。导航LDPC如果用IEEE 802.11ad风格的短帧码15到20次基本够如果是自己随便写的稀疏H先测5次和20次的FER差差得大就说明矩阵列重设计有问题。多进制码不是靠堆迭代救场的。damp的典型范围是0.3到0.7。概率域BP在SNR较高时消息会在两个候选符号之间反复横跳表现为FER在某个SNR点附近抖动。把damp从0.5提到0.7能压住抖动但每轮迭代的有效步长变小所以必须同步加几轮max_iter。如果damp超过0.8收敛会慢到让人误以为死循环。信道先验P0是最容易错的一环。导航接收机的软信息来自QPSK解映射接收复符号r和星座点x_s之间距离越近P0(s)越大。正确的做法是算欧氏距离后放进高斯表达式再做减去最大值的归一化。下面的代码片段是一般做法% rec: 1x2复数接收向量const: 1x4 QPSK星座点 % 星座点顺序必须和gf4符号0~3一一对应 const exp(1j * pi / 4 * [0 1 3 2]); % 按你的Gray映射重排 P_tmp zeros(1, 4); for s 1:4 d2 abs(rec - const(s)).^2; P_tmp(s) -d2 / N0; end P0 exp(P_tmp - max(P_tmp)); P0 P0 / sum(P0);exp(P_tmp - max(P_tmp))这种写法是为了防止N0很小时d2/N0直接溢出成-inf。先减最大值再exp得到的软概率和原始高斯比例一致只是整体缩放不影响BP后验。错误率最高的场景是const数组顺序和gf4符号编号不一致比如符号2在星座上是“01”你却在const里放了“10”那么先验概率即使再准BP也会把约束往错误方向推。4.2 一张表说清四个典型故障现象可能原因定位 / 调整FER一直0.5GF(4)乘法表或加法表写错手动查gf4_mul(2,2)3、gf4_mul(3,3)2迭代不收敛伴随式永远不满足H含2/3系数但校验更新没换元在卷积前把消息置换到yh*s域卷积后再反置换回来概率域出现NaN变量更新后没归一化检查每一处out/sum(out)C2V也建议加保护先验分布接近均匀FER反而不如硬判决星座映射和GF符号顺序错位用无噪声回环打印四个概率正确项应接近1H含非1系数是多进制LDPC和二进制实现最本质的分叉点。我的做法是在gf4_conv_cyc之前先把每条消息按h做符号置换如果校验方程是Σ h_i c_i0令y_ih_i c_i那么新消息的分布满足M(y)M(h^{-1}y)。置换后所有h都变成1就能继续用同一个异或卷积函数。卷积输出是关于y的分布要再按h^{-1}逆置换回c。这里只要忘掉逆置换伴随式就会在几轮迭代后爆出全零误判而错误帧率看起来却像模像样。4.3 把H的非1系数用换元方式接进现有循环如果你手里的校验矩阵不是全1千万不要去修改gf4_conv_cyc的内部循环。正确姿势是在进入校验节点更新前对当前行的每条消息做一个符号置换卷积完成后对输出消息再做逆置换。这个操作本质上是GF(4)上的排列不改变BP的消息总和也不引入额外信息。写成代码时需要一张gf4_mul表和一张gf4_inv表一般可以这样查% h为行稀疏矩阵H(r,col)的GF(4)系数 % MsgIn(idx,:)是符号0~3的概率向量 % 正向换元把符号s映射到h*s tmp zeros(1, 4); for s 0:3 d gf4_mul(h 1, s 1); % h*s tmp(d 1) tmp(d 1) MsgIn(idx, s 1); end MsgIn(idx, :) tmp;逆置换就是把上面的d换成h^{-1}*s。GF(4)里非零元素的逆是它自己的一部分1的逆是12的逆是33的逆是2所以逆置换表可以直接从gf4_mul反查。把正向和逆向写成两个小函数H里就算混入2和3主循环也不用再改。这个换元也解释了为什么很多现成的MATLAB二进制LDPC代码无法直接改成多进制二进制消息根本不需要做这种域上置换自然也没有人会提前留这个口子。5. 用无噪声回环和FER曲线验收多进制LDPC BP程序验证多进制BP程序按顺序做三步无噪声回环、单帧加噪、多帧FER统计。无噪声回环专门用来暴露代数语义错误比如gf4_add表的1索引写错、伴随式检查用了模2、星座点顺序和gf4符号不一致。做法是把已知码字转成one-hot先验直接喂给解码器cw random_gf4_codeword(); % 从生成矩阵取一个GF(4)码字 P0 zeros(length(cw), 4); for i 1:length(cw) P0(i, cw(i) 1) 1; % one-hot先验 end [hard, iter] gf4_bp_decode(P0, H, 20, 0.5); assert(isequal(hard(:), cw(:)), 无噪声回环失败);这一步如果失败不要急着看迭代过程先检查GF(4)加法表是不是从1索引开始、伴随式计算有没有漏掉h系数。无噪声回环通过后再丢进AWGN信道做FER统计。一个简单的验收循环如下fer zeros(1, 4); for snr_i 1:4 bad 0; for frame 1:200 cw random_gf4_codeword(); rx add_awgn(qpsk_map(cw), snr_dB(snr_i)); P0 symbol_prob(rx, N0, const); [hard, ~] gf4_bp_decode(P0, H, 20, 0.5); if any(hard(:) ~ cw(:)) bad bad 1; end end fer(snr_i) bad / 200; end观察FER曲线时有一条经验多进制BP程序如果某个SNR点FER比硬判决还高九成是P0的星座映射表错了而不是迭代次数不够。GNSS导航电文常用QPSK Gray映射GF(4)符号会按同相和正交拆成两个bit因此先验概率里“符号0对应哪个星座点”必须和发端bit序一致。我最常做的收尾验证是在解映射后打印第一个正确定位符号的四个概率值正确项应该明显超过0.5如果四个值都接近0.25先别调damp回去看const数组的顺序。导航LDPC还有一个比普通AWGN仿真多出的步骤接收信号往往带残余多普勒相位所以在把复符号送进symbol_prob之前先用导频或相位解卷绕把星座转正。这个补偿必须发生在星座映射之前不然BP算得再准也只会把相位误差当成噪声。最终确认过不了这一关时把GF(4)元素和QPSK星座点放进同一个由表驱动的函数里维护比在解码循环里反复试max_iter收敛得更快。本文还有配套的精品资源点击获取 MATLAB实现GF(4)多进制LDPC的BP解码与调优实战
简介一套基于BP信念传播算法的多进制LDPC码MATLAB仿真资源专注卫星导航场景下的纠错编译码研究适合通信工程、导航制导与控制专业学生及科研人员使用可帮助快速理解多进制LDPC码的Tanner图构造、BP消息传递机制以及在复杂信道下的完整译码实现。资源包共包含四个文件一个m脚本作为主译码程序实现基于信念传播的迭代译码三个txt文本文件分别保存该多进制LDPC码的指数信息、输入数据与码字元素三者配套驱动译码器运行压缩包整体仅4KB轻量易用便于部署到各类仿真环境。目前已有856人学习/下载特别适合课程设计、毕业设计或工程项目预研阶段使用。多进制LDPC码相比二进制LDPC码可获得更高编码增益在卫星导航这类易受电离层闪烁、多径衰落等影响的链路中能显著提升抗干扰能力、降低定位误差。运行代码可直接得到不同信噪比下的误码率性能曲线还可调整BP迭代次数与消息传递参数观察译码收敛速度及误码率变化从而优化译码策略为实际多进制LDPC编译码方案设计与参数选择提供可靠的参考依据。1. 多进制LDPC的BP解码为什么值得在MATLAB里重写一遍搜索“matlab_ldpc64_BP.rar”这类名字的人大多数想找一个能跑的二进制LDPC BP参考再改成多进制但多进制LDPC的BP和二进制BP不是同一个东西。变量节点消息从0/1的LLR变成GF(q)上的q维概率向量校验节点要做域上卷积代码里多出来的不只是循环而是整套消息模型。更反直觉的一点是在QPSK这类2比特调制下GF(4)符号可以直接映射到星座点短帧性能反而比同码率二进制LDPC更容易逼近香农界。这里用GF(4)作为最小例子拆开BP的变量更新、校验更新和参数调优给出一段能跑的MATLAB代码并说明导航LDPC场景怎么验收它。2. 多进制LDPC为什么要把BP消息从标量改成向量2.1 校验矩阵从GF(2)搬到GF(q)后码字和伴随式都变了二进制LDPC的校验矩阵是0/1矩阵码字满足H c^T0所有加法都在GF(2)上做。多进制LDPC先把符号域换成GF(2^m)也就是把原来的一位比特变成一个m比特符号。GF(4)是其中最小的例子符号集{0,1,2,3}加法为按位异或乘法按本原多项式x^2x1构造。这样一行校验方程不再是几个bit的异或而是GF(4)上系数与符号的线性组合。伴随式是否正确也不能只看奇偶性必须按GF(4)乘法表计算。换来的是码字符号与高阶调制星座的天然对应一个GF(4)符号带2bit正好放进QPSK的一个星座点GF(8)符号可匹配8PSKGF(16)可匹配16QAM。在导航这类信号帧短、功率受限的链路上短帧多进制码比二进制码更容易在同等复杂度下逼近香农限。真正变化的是约束的结构。二进制校验只回答“这几位加起来是不是偶数”GF(4)校验回答的是“这几位按域上的加权和是不是0”。生成矩阵和校验矩阵都不是0/1矩阵后编码也要相应换成GF(4)矩阵乘法解码端判断停止迭代时更不能用二进制的模2伴随式。这个区别在调试时最容易暴露用二进制LDPC的伴随式函数去验GF(4)码字几乎每一帧都会被误判为错帧。2.2 BP消息变成向量先验、V2C、C2V三种消息全是q维分布二进制BP里一条边缘消息是一个标量通常做成LLR。多进制BP里每个变量节点有q种可能取值所以任何经过这条边的消息都是一个长度为q的概率向量。记P0(x)P(c_ix)为信道先验变量节点发给校验节点的V2C消息为M_{v→c}(x) ∝ P0(x) * ∏_{c∈N(v){c}} M_{c→v}(x)校验节点发给变量节点的C2V消息则是在满足该校验方程的前提下对其他关联变量做GF(4)符号的“异或卷积”。一个度数为d_c的校验节点如果按穷举组合计算复杂度是O(d_c q^2)。GF(4)下q4穷举非常快GF(64)、GF(256)再这么做就会卡成幻灯片那时才需要FFT-QSPA或扩展最小和。导航上的LDPC很多是短帧GF(4)或GF(8)已经够用所以先不用急着上快速算法。消息从标量变成向量后归一化也变成向量归一化。每一轮变量更新结束必须把每条消息除以它自己的sum否则几十轮迭代后概率会下溢成NaN。这一点和二进制BP的LLR域做法很不一样很多人把二进制代码搬过来时漏掉的往往是这种看起来不起眼的归一化。2.3 先建好GF(4)的加法表和乘法表整个BP才有地基GF(4)加法等价于按位异或所以可以直接用MATLAB的bitxor生成。乘法表对应本原多项式x^2x1符号2是α符号3是α1。下面的代码先建表再用一个全1校验向量做自检。q 4; gf4_add zeros(q, q); for a 0:q-1 for b 0:q-1 gf4_add(a1, b1) bitxor(a, b) 1; % 存1索引方便后面查表 end end gf4_mul [0 0 0 0; 0 1 2 3; 0 2 3 1; 0 3 1 2]; % 自检全1向量 [1 2 3] 的GF(4)和应为0 acc 0; for c [1 2 3] acc bitxor(acc, c); end assert(acc 0);注意gf4_add存的是1索引a2、b3时bitxor(2,3)1存进矩阵的是2也就是符号1所对应的索引。gf4_conv这类函数直接拿这个表做输出索引能少一层“符号索引互换”的if。gf4_mul表则按0索引存符号查表时用gf4_mul(a1,b1)。两种索引不要混用这是多进制LDPC排错里最常见的暗坑。GF(4)乘法表如下行和列都是符号0到3乘012300000101232023130312这个表里223、332是特征2域的典型性质。写代码时如果这三个“乘法对角线”掉错了后面所有校验节点更新全会错而且不会报异常只是FER永远停在0.5附近。3. 在MATLAB里搭一个GF(4) LDPC BP解码器3.1 用三维数组存V2C和C2V避免cell的隐式开销这里采用m行n列乘q层的三维数组V2C(r, col, s)表示第r个校验节点到第n列变量节点之间、关于符号s的消息概率C2V方向相反。三维数组的好处是切片方便变量节点更新时可以直接用squeeze(C2V(r2, col, :))取出一行概率再用点乘做外信息乘积。虽然用cell数组也能写但每轮迭代都要cell2mat代码更长且速度更慢。初始化的原则是H为0的地方不参与任何更新所以V2C和C2V那里即使有数值也不会被读到。V2C只在find(H(:,col) ~ 0)的行填入P0C2V初始化为全1含义是“校验节点还未提供任何约束信息”。第一次变量更新时每条消息都等于信道先验P0和推导一致。3.2 完整的GF(4) BP主循环和GF(4)异或卷积下面这段代码能直接保存成gf4_bp_decode.m运行。它把变量节点更新、校验节点更新和伴随式检查放在同一个循环里。为了突出BP骨架这段代码的H暂时只处理非零元素为1的情况真实导航工程里的H会出现2、3通用化方法在下一章给。function [hard, iter] gf4_bp_decode(P0, H, max_iter, damp) % P0: n-by-4 信道先验概率P0(i,x1) P(c_i x) % H: m-by-n GF(4)校验矩阵本例只处理非零元素为1的情况 % max_iter: 最大迭代次数 % damp: 消息阻尼0~1建议取0.5 % hard: 译码符号0~3 q 4; [m, n] size(H); V2C zeros(m, n, q); C2V ones(m, n, q); gf4_add zeros(q, q); for a 0:q-1 for b 0:q-1 gf4_add(a1, b1) bitxor(a, b) 1; end end % 变量到校验消息初始化 for col 1:n rows find(H(:, col) ~ 0); for r rows V2C(r, col, :) P0(col, :); end end for it 1:max_iter % 变量节点更新信道先验乘以外来C2V消息 for col 1:n rows find(H(:, col) ~ 0); for r rows out P0(col, :); for r2 rows if r2 ~ r out out .* squeeze(C2V(r2, col, :)); end end out out / sum(out); % 阻尼旧消息占一部分新消息占一部分 V2C(r, col, :) (1 - damp) * squeeze(V2C(r, col, :)) damp * out; end end % 校验节点更新GF(4)上的异或卷积 for r 1:m cols find(H(r, :) ~ 0); if length(cols) 2 continue; end MsgIn zeros(length(cols), q); for idx 1:length(cols) MsgIn(idx, :) squeeze(V2C(r, cols(idx), :)); end MsgOut zeros(size(MsgIn)); for p 1:length(cols) others setdiff(1:length(cols), p); dist [1, 0, 0, 0]; % 单位元 for idx others dist gf4_conv_cyc(dist, MsgIn(idx, :), gf4_add); end MsgOut(p, :) dist; end for idx 1:length(cols) C2V(r, cols(idx), :) MsgOut(idx, :); end end % 后验概率与硬判决 post P0; for col 1:n rows find(H(:, col) ~ 0); for r rows post(col, :) post(col, :) .* squeeze(C2V(r, col, :)); end post(col, :) post(col, :) / sum(post(col, :)); end [~, hard] max(post, [], 2); hard hard - 1; % 索引1..4变成符号0..3 % 伴随式检查注意本例H非零元素全为1 synd_ok true; for r 1:m cols find(H(r, :) ~ 0); acc 0; for idx 1:length(cols) acc bitxor(acc, hard(cols(idx))); end if acc ~ 0 synd_ok false; break; end end if synd_ok iter it; return; end end iter max_iter; end function z gf4_conv_cyc(x, y, gf4_add) z zeros(1, length(x)); for a 1:length(x) for b 1:length(y) z(gf4_add(a, b)) z(gf4_add(a, b)) x(a) * y(b); end end end主函数里变量节点更新先取P0再乘上所有邻居C2V消息最后做一次归一化。这一步的物理含义是消息不能包含自己发出的信息所以要把被更新的那个邻居排除在外。阻尼系数damp在这里不是必需但对概率域BP很有用因为高信噪比时消息会在两个状态之间震荡阻尼能压低振荡。damp越大新消息占的比例越小迭代收敛越慢所以要配合max_iter一起调。校验节点更新用gf4_conv_cyc做GF(4)异或卷积。dist[1,0,0,0]相当于GF(4)离散卷积的单位元它本身不代表任何信道信息只让循环写法统一。对每一个目标边p它取其他边的消息做卷积得到的就是“在其他变量取值组合下当前符号必须取某个值”的概率分布。由于卷积保持概率和为1校验输出不需要再归一化如果你在工程里改了域大小建议仍加一次归一化避免浮点漂移。这个函数能直接跑通3×6这类小矩阵。测试时可以把P0设置为某个已知码字的one-hot概率硬判决应该一次就通过伴随式并且iter返回1。如果连无噪声回环都过不去问题多半在gf4_add表或伴随式检查而不是在迭代更新上。4. 多进制LDPC BP的必调参数与四个排错入口4.1 max_iter、damp、信道先验是三个旋钮max_iter的典型范围是5到20。GF(4)消息是4维向量信息在校验节点之间传播的速度比二进制标量慢迭代太少会让FER曲线出现平台加太多又只增加耗时。导航LDPC如果用IEEE 802.11ad风格的短帧码15到20次基本够如果是自己随便写的稀疏H先测5次和20次的FER差差得大就说明矩阵列重设计有问题。多进制码不是靠堆迭代救场的。damp的典型范围是0.3到0.7。概率域BP在SNR较高时消息会在两个候选符号之间反复横跳表现为FER在某个SNR点附近抖动。把damp从0.5提到0.7能压住抖动但每轮迭代的有效步长变小所以必须同步加几轮max_iter。如果damp超过0.8收敛会慢到让人误以为死循环。信道先验P0是最容易错的一环。导航接收机的软信息来自QPSK解映射接收复符号r和星座点x_s之间距离越近P0(s)越大。正确的做法是算欧氏距离后放进高斯表达式再做减去最大值的归一化。下面的代码片段是一般做法% rec: 1x2复数接收向量const: 1x4 QPSK星座点 % 星座点顺序必须和gf4符号0~3一一对应 const exp(1j * pi / 4 * [0 1 3 2]); % 按你的Gray映射重排 P_tmp zeros(1, 4); for s 1:4 d2 abs(rec - const(s)).^2; P_tmp(s) -d2 / N0; end P0 exp(P_tmp - max(P_tmp)); P0 P0 / sum(P0);exp(P_tmp - max(P_tmp))这种写法是为了防止N0很小时d2/N0直接溢出成-inf。先减最大值再exp得到的软概率和原始高斯比例一致只是整体缩放不影响BP后验。错误率最高的场景是const数组顺序和gf4符号编号不一致比如符号2在星座上是“01”你却在const里放了“10”那么先验概率即使再准BP也会把约束往错误方向推。4.2 一张表说清四个典型故障现象可能原因定位 / 调整FER一直0.5GF(4)乘法表或加法表写错手动查gf4_mul(2,2)3、gf4_mul(3,3)2迭代不收敛伴随式永远不满足H含2/3系数但校验更新没换元在卷积前把消息置换到yh*s域卷积后再反置换回来概率域出现NaN变量更新后没归一化检查每一处out/sum(out)C2V也建议加保护先验分布接近均匀FER反而不如硬判决星座映射和GF符号顺序错位用无噪声回环打印四个概率正确项应接近1H含非1系数是多进制LDPC和二进制实现最本质的分叉点。我的做法是在gf4_conv_cyc之前先把每条消息按h做符号置换如果校验方程是Σ h_i c_i0令y_ih_i c_i那么新消息的分布满足M(y)M(h^{-1}y)。置换后所有h都变成1就能继续用同一个异或卷积函数。卷积输出是关于y的分布要再按h^{-1}逆置换回c。这里只要忘掉逆置换伴随式就会在几轮迭代后爆出全零误判而错误帧率看起来却像模像样。4.3 把H的非1系数用换元方式接进现有循环如果你手里的校验矩阵不是全1千万不要去修改gf4_conv_cyc的内部循环。正确姿势是在进入校验节点更新前对当前行的每条消息做一个符号置换卷积完成后对输出消息再做逆置换。这个操作本质上是GF(4)上的排列不改变BP的消息总和也不引入额外信息。写成代码时需要一张gf4_mul表和一张gf4_inv表一般可以这样查% h为行稀疏矩阵H(r,col)的GF(4)系数 % MsgIn(idx,:)是符号0~3的概率向量 % 正向换元把符号s映射到h*s tmp zeros(1, 4); for s 0:3 d gf4_mul(h 1, s 1); % h*s tmp(d 1) tmp(d 1) MsgIn(idx, s 1); end MsgIn(idx, :) tmp;逆置换就是把上面的d换成h^{-1}*s。GF(4)里非零元素的逆是它自己的一部分1的逆是12的逆是33的逆是2所以逆置换表可以直接从gf4_mul反查。把正向和逆向写成两个小函数H里就算混入2和3主循环也不用再改。这个换元也解释了为什么很多现成的MATLAB二进制LDPC代码无法直接改成多进制二进制消息根本不需要做这种域上置换自然也没有人会提前留这个口子。5. 用无噪声回环和FER曲线验收多进制LDPC BP程序验证多进制BP程序按顺序做三步无噪声回环、单帧加噪、多帧FER统计。无噪声回环专门用来暴露代数语义错误比如gf4_add表的1索引写错、伴随式检查用了模2、星座点顺序和gf4符号不一致。做法是把已知码字转成one-hot先验直接喂给解码器cw random_gf4_codeword(); % 从生成矩阵取一个GF(4)码字 P0 zeros(length(cw), 4); for i 1:length(cw) P0(i, cw(i) 1) 1; % one-hot先验 end [hard, iter] gf4_bp_decode(P0, H, 20, 0.5); assert(isequal(hard(:), cw(:)), 无噪声回环失败);这一步如果失败不要急着看迭代过程先检查GF(4)加法表是不是从1索引开始、伴随式计算有没有漏掉h系数。无噪声回环通过后再丢进AWGN信道做FER统计。一个简单的验收循环如下fer zeros(1, 4); for snr_i 1:4 bad 0; for frame 1:200 cw random_gf4_codeword(); rx add_awgn(qpsk_map(cw), snr_dB(snr_i)); P0 symbol_prob(rx, N0, const); [hard, ~] gf4_bp_decode(P0, H, 20, 0.5); if any(hard(:) ~ cw(:)) bad bad 1; end end fer(snr_i) bad / 200; end观察FER曲线时有一条经验多进制BP程序如果某个SNR点FER比硬判决还高九成是P0的星座映射表错了而不是迭代次数不够。GNSS导航电文常用QPSK Gray映射GF(4)符号会按同相和正交拆成两个bit因此先验概率里“符号0对应哪个星座点”必须和发端bit序一致。我最常做的收尾验证是在解映射后打印第一个正确定位符号的四个概率值正确项应该明显超过0.5如果四个值都接近0.25先别调damp回去看const数组的顺序。导航LDPC还有一个比普通AWGN仿真多出的步骤接收信号往往带残余多普勒相位所以在把复符号送进symbol_prob之前先用导频或相位解卷绕把星座转正。这个补偿必须发生在星座映射之前不然BP算得再准也只会把相位误差当成噪声。最终确认过不了这一关时把GF(4)元素和QPSK星座点放进同一个由表驱动的函数里维护比在解码循环里反复试max_iter收敛得更快。本文还有配套的精品资源点击获取