一、理论基础此部分参考资料LDPC码一种前向纠错码基础 译码算法 - 知乎1.1 概述置信传播Belief PropagationBP算法在编码理论中又常被称为和积算法Sum-Product AlgorithmSPA是一种在概率图模型上进行统计推断的消息传递机制其核心思想是通过在Tanner图的变量节点VN与校验节点CN之间迭代传递置信度信息逐步逼近最大后验概率MAP。BP算法由朱迪亚·珀尔Judea Pearl于1982年提出最初用于贝叶斯网络和马尔可夫随机场的概率推断后被引入通信领域并发展为LDPC码的核心译码方法。但BP译码算法不是只能用于LDPC码理论上可以用于任何拥有因子图Factor Graph表示的码型包括所有的线性分组码甚至卷积码Turbo码。Turbo码Turbo码的迭代译码本质上就是在两个分量码的因子图之间交换信息属于BP算法的一个特例。Polar码极化码Polar码最经典的就是SC串行消除和SCL串行消除列表译码但Polar码完全可以使用BP算法译码。Polar码的生成矩阵GN具有递归结构可以展开成一个递归蝶形因子图在此图上运行BP算法不仅可以实现并行译码降低时延还可以通过Soft-IC等手段进一步提升性能。PAC码极化调整卷积码PAC码在Polar码的极化编码之前加入了卷积预变换虽然无法再用递归蝶形因子图完整表示整个过程但其依然有确定的校验矩阵理论上也可以使用BP译码算法只是其校验矩阵的密度非常高直接运行BP算法的性能会很差。卷积码经典的BCJR算法其实就是定义在无环图Trellis图/网格图上的BP算法。RA码这是一类结构化的类LDPC码天然适合BP译码乘积码常用于光通信和存储其迭代硬判决或软判决译码也是BP思想的体现。1.2 图论基础与结构化表示前面提到BP算法是一个通用的框架理论上适用于所有可以用因子图表示的概率模型。在线性分组码中通过校验矩阵H构建的Tanner图就是因子图的一个特例。下面以LDPC码为例。1.2.1 线性分组码的矩阵定义LDPC码属于线性分组码的范畴。一个 (n,k) 线性分组码将 k 比特的信息序列映射为 n 比特的码字序列从而引入 n-k 个冗余比特用于纠错。该码由一个 m×n 的奇偶校验矩阵H唯一确定其中mn-k 。对于任意合法码字 x [x1,x2,...,xn]必须满足正交约束这意味着码字中的比特必须满足 m 个线性校验方程。LDPC码的核心特征在于校验矩阵H是“低密度”的即矩阵中绝大多数元素为“0”只有极少数元素为“1”。通常每行和每列中“1”的数量远小于 n 和 m 。这种稀疏性保证了BP算法的计算复杂度随码长线性增长O(n)而非指数增长。1.2.2 Tanner图二部图模型1981年Tanner提出用二部图Bipartite Graph来表示奇偶校验矩阵这种图后来被称为Tanner图。Tanner图由两个不同的节点集合组成变量节点Variable Node, VN和校验节点Check Node, CN这些节点通过边相互连接以图形化方式展示了校验矩阵H的结构。因此Tanner图实际上就是校验矩阵的可视化表示是定义BP算法消息传递路径的拓扑基础。例给定一个码长为7信息位长度为4的LDPC码的校验矩阵H的每一行对应一个校验节点也可以说对应一个校验方程。的每一列对应一个变量节点。表示第 j 个变量节点和第 i 个校验节点相连。根据该校验矩阵就可以画出对应的可视化Tanner图只要建立了这个Tanner图就可以在上面跑BP算法变量节点向校验节点发信息校验节点处理后再发回变量节点如此迭代。关于Tanner图还有几个比较重要的定义1边Tanner图中的边仅连接变量节点和校验节点不会连接同类节点。当且仅当校验矩阵元素时校验节点与变量节点之间存在一条边。2度分布Degree Distribution变量节点度与变量节点相连的边的数量等于对应列的汉明重。它表示该比特参与了多少个校验方程。校验节点度与校验节点相连的边的数量等于对应行的汉明重。它表示该方程约束了多少个比特。如果图中所有变量节点的度相同且所有校验节点的度也相同则该码称为规则LDPC码。如果度数不统一则称为非规则LDPC码。研究表明设计良好的非规则LDPC码通常具有比规则码更优异的瀑布区性能因为高度数的变量节点能更早收敛从而辅助低度数节点译码。2环Cycles与围长Girth环从一个节点出发经过若干条不同的边回到该节点的闭合路径。环的长度是指路径上边的数量。由于二部图的性质环长必然是偶数4、6、8...。围长图中长度最短的环的长度。在无环图树中BP算法可以精确收敛到最大后验概率。然而为了获得良好的距离特性实用的LDPC码必然包含环。短环特别是长为4的环会破坏BP算法的独立性假设导致信息错误在局部循环并被放大引发译码失效或震荡。因此构造LDPC码的一个重要目标就是消除短环最大化围长。1.3 消息传递机制在开始公式推导之前先回答两个问题更有利于对后面公式的理解在Tanner图上传播的消息具体是什么含义为什么BP算法是有效的1传播消息的具体含义在实际工程实现中为了简化乘法运算为加法运算通常在对数域进行计算传播的消息是对数似然比Log-Likelihood Ratio, LLR定义为如果 L(x) 0表示该比特更倾向为0。如果 L(x) 0表示该比特更倾向为1。| L(x) | 的大小表示置信度绝对值越大越确定。A. 变量节点 - 校验节点v - c含义根据信道给我的初始信息以及除了你该校验节点c之外的其他校验节点告诉我的信息我认为我自己是0或1的可能性是多少。B. 校验节点 - 变量节点c - v含义为了满足我的校验方程即所有连接到我的变量模2和必须为0根据除了你该变量节点v之外的其他变量节点目前的状态我推断你应该取值的可能性是多少。两个维度方向逻辑推断如果其他变量节点的模2和是0或1你该变量节点v也应该是0或1才能满足整个校验方程的模2和为0决定了LLR的正负号。强度置信度表示我推断你该变量节点v是0或1的这个消息有多确定决定了LLR绝对值的大小。校验节点传出的置信度受限于与之相连的绝对值最小最不可靠的那个变量节点这里用“木桶效应”比较好理解因为只要有一个变量节点的消息不靠谱例如L(x)0.01那么整个推断就是不靠谱的。2BP算法的有效性局部交互达成全局共识BP译码的核心思想是将一个复杂的全局概率计算问题分解为无数个简单的局部计算问题。每一个节点并不需要知道整个图的结构它只与它的“邻居”进行通信。变量节点负责收集关于该比特概率的信息校验节点负责检查这些概率是否符合奇偶校验约束。这种通信是双向迭代的。随着迭代次数增加信息在图中扩散原本不确定的节点受噪声严重干扰的比特会通过其他节点的强有力信息逐渐被纠正最终整个网络达成共识。利用了稀疏性BP算法之所以能成为LDPC码的核心译码方法就是因为LDPC码的校验矩阵是非常稀疏的这意味着每个变量只参与少数几个校验每个校验只关联少数几个变量围长大消息能够有效传播且这种稀疏性保证了复杂度随码长呈线性增长而非最大似然译码的指数级增长。外信息的交换这是BP算法的精髓当节点A向节点B传递消息时不会把B刚才告诉的它的消息再传回B这防止了信息在两个节点之间死循环保证了每次迭代传递的都是“新鲜”的独立证据。总而言之BP译码算法可以看成是一个迭代的“投票”过程。每个比特最初只有信道的模糊观测通过不断的变量节点汇总意见和校验节点施加约束之间的信息传播消除不确定性知道所有校验方程都被满足。二、详细公式推导此部分参考资料LDPC译码原理公式推导及其matlab代码实现超详细-CSDN博客高斯白噪声信道下对数似然比的推导--- QPSKBPSK - 哔哩哔哩高斯白噪声信道下后验概率的计算 - 哔哩哔哩LDPC码最小和译码优化及调度算法的研究 - 中国知网基于深度学习的LDPC码译码算法研究 - 中国知网卫星通信中LDPC译码器的研究与实现 - 中国知网2.1 引理1后面的推导需要用到引理1给定一个长为的相互独立的二进制序列第个比特为 1 的概率为则二进制序列 a 中包含偶数个 1 的概率为包含奇数个 1 的概率为证明如下采用构造函数和数学归纳法进行证明。首先构造关于的次多项式通过数学归纳法证明的系数为序列中包含个 1 的概率。1当时显然的系数为序列中包含 1 个 1 的概率的系数为序列中包含 0 个 1 的概率。2假设当时的次多项式为的系数为其序列中包含个 1 的概率。3当时的次多项式 为由上式可求得的系数为其中表示第个比特为 0 的概率表示的序列中有个 1 的概率表示第个比特为 1 的概率表示的序列中有个 1 的概率因此当时的系数也表示序列中包含个 1 的概率。综上通过数学归纳法证明了的系数为序列中包含个 1 的概率。在此基础上再求序列中包含偶数个 1 的概率。即需要求 m 次多项式中所有偶数次幂项的系数之和为此可以再构造一个函数则会消掉所有奇数次幂项得到所有偶数次幂项和的两倍再令得到序列中包含偶数个 1 的概率则序列中包含奇数个 1 的概率2.2 前置知识在介绍具体的译码算法之前还有一点前置知识需要先搞清楚整个通信链路。发送端的原始编码码字为经过调制得到发送信号再通过AWGN信道传输引入高斯白噪声接收端收到的信号其中噪声服从高斯分布。译码就是通过接收信号推断原始编码码字的过程推断的依据为原始编码码字是满足校验矩阵的如果译码结果也满足相同的校验矩阵即我们就认为找到了最可能的原始编码码字。2.3 概率域BP译码算法符号说明两个概念后验概率通信系统中我们关心的是在收到信号 y 的情况下判断发送方发送的比特是 0 或 1 的概率可以表示成和一般被称为后验概率。外信息原则不将节点A告诉我的信息回传给A以防止信息自我增强正反馈导致的过早收敛和偏差。1初始化对于每个变量节点根据信道接收值计算初始后验概率这是算法的输入代表了仅由信道观测提供的信息也是变量节点传递给校验节点的初始信息2校验节点更新校验节点传递给变量节点的外信息校验节点的任务是根据连接到它的变量节点除变量节点的信息计算它对变量节点的约束意见理解含义表示在第次迭代中校验节点推断变量节点为 0 的概率等价于除开该变量节点外的其他变量节点的模2和应该为0的概率这样才能满足连接到该校验节点的所有变量节点的模2和为0的校验约束。也就等价于求其他变量节点中包含偶数个 1 的概率即前面已经证明过的引理1。表示在第次迭代中除外的其他变量节点为 1 的概率由第次迭代时变量节点传递给校验节点的外信息确定。的含义同理可得。3变量节点更新变量节点传递给校验节点的外信息变量节点的任务是综合来自信道的初始信息和连接到它的校验节点除校验节点传来的外信息。其中的是校正因子使每次计算出的。理解含义表示在第次迭代中变量节点为 0 的概率这个概率包括初始为0的概率、除开该校验节点外的其他校验节点推测该变量节点为 0 的概率这些概率值相乘得到了第次迭代中的概率。的含义同理可得。4译码与终止条件在每次迭代结束时计算每个变量节点的总后验LLR这里包含了所有校验节点的信息其中的是校正因子目的是使。硬判决若则判否则判。校验计算。若说明所有校验方程满足译码成功提前终止。若且未达到最大迭代次数则回到步骤2继续下一轮迭代。若达到最大次数仍未满足校验条件则宣布译码失败。2.4 对数域BP译码算法前面介绍了概率域上的BP译码算法但在实际工程实现中并不会使用概率域BP译码算法因为其包含大量的连乘运算计算复杂度非常高。常见的做法是将概率信息用对数似然比表示得到对数域上的BP译码算法大量的连乘运算可以转化为加法运算。2.4.1 标准和积算法Sum-Product Algorithm, SPA和积算法是BP算法的标准实现不包含任何近似被公认为性能最优的软判决迭代算法。1初始化同样的第一步是初始化算法的输入。对于每个变量节点根据信道接收值计算初始信道LLR即变量节点发送给校验节点的初始信息。对于AWGN信道、BPSK调制下面对该公式进行详细推导对数似然比LLR的定义如前所述根据接收信号判断发送端发送的比特是 0 或 1 的概率可以表示为后验概率或。一般情况下我们更关心两者的比值因为如果大于我们就有有理由判断否则我们更倾向于判断为了简化运算一般还要再取对数这就是对数似然比的定义式判决准则相应地变为信道初始LLR的计算得到了LLR的定义式后却不能直接根据该定义式进行计算原因是后验概率只是一个离散分布即在接收端收到某一个具体的值后只可能是 0 或者 1我们无法写出它的函数解析式。在2.2节中提到过编码码字经过调制得到发送信号与是一一映射的在BPSK中0 一般映射为 1 1 一般映射为 -1 经过AWGN信道传输加入高斯白噪声得到接收信号因为噪声是服从高斯分布的所以接收信号也服从高斯分布。例如发送端发送了即也依然服从高斯分布因为给定的情况下是一个常数。噪声服从均值为 0 方差为的高斯分布可表示为满足以下概率密度函数代入这个式子可以写成进一步等价于所以式子最终可以写成这里的逻辑有点绕但其实很好理解因为和描述的是同一件事情。表示噪声取某个值的概率表示已知等价于已知因为它们是一一映射的的条件下取某个值的概率因为当已知时是一个固定的常数的值就完全由决定所以这两个表达式是完全等价的。还是前面的例子发送端发送了即这个时候求取某个值的概率也就是求时取某个值的概率。被称为似然概率根据上面的式子它是服从高斯分布的。还记得前面说的后验概率是一个离散分布无法写出它的函数解析式。因此我们可以利用贝叶斯公式将其转化为服从高斯分布的似然概率与先验概率的表达式到这里就可以对LLR进行计算了这里分成了两部分一是信道信息代入前面求的概率密度函数注意这里中代入的值需要看BPSK调制的具体映射关系这里是 0-11--1如果是更高阶的调制整个表达式会发生变化求解更加复杂。具体可以看这位博主的文章高斯白噪声信道下对数似然比的推导--- QPSKBPSK - 哔哩哔哩二是先验信息一般情况下发送比特是0或1的概率相等因此通常为0。因此信道初始LLR的计算式最终可以化简为这是对数域译码算法中信道初始化的常用式子基本上是作为一个结论使用的少有详细推导。其中是噪声方差它的求法也有一个推导过程这里就不赘述了以后有机会再补上。2校验节点更新校验节点传递给变量节点的外信息根据前面的定义对数域BP译码算法的校验节点更新公式为在数学中双曲正切函数有以下两个公式成立函数及其导数的图像如下且有以下几个关系式成立可以利用这几个公式对进行化简故校验节点更新的公式最终可以化简为这是SPA的标准公式。由于 tanh 的乘积计算复杂该标准公式还有另外一种表达形式可以使用函数将乘积转换为加法是符号部分。是幅度部分。符号部分表示校验节点推测变量节点是 0 或 1方向由输入消息的符号异或决定幅度部分表示校验节点推测变量节点是 0 或 1 的确定性强度。函数定义为是自逆函数即同时是单调递减函数且为正值当当时当时。3变量节点更新变量节点传递给校验节点的外信息变量节点更新的操作在LLR域即为简单的加法4译码与终止条件在每次迭代结束时计算每个变量节点的总后验LLR这里包含了所有校验节点的信息。硬判决若则判否则判。校验计算。若说明所有校验方程满足译码成功提前终止。若且未达到最大迭代次数则回到步骤2继续下一轮迭代。若达到最大次数仍未满足校验条件则宣布译码失败。2.4.2 最小和算法Min-Sum Algorithm, MSA虽然SPA提供了最优的误码率性能但 tanh函数或函数作为非线性运算在FPGA或ASIC实现中不仅占用大量面积还容易受到定点数量化误差的影响。为了解决这一问题工程界广泛采用基于”最小近似“的简化算法主要改进了校验节点更新的公式推导如下函数中当较大时会迅速衰减。因此在求和时主要的贡献来自于数值最大的项而是单调递减函数这意味着最大的值对应最小的输入幅度。因此我们可以用最小值来近似整个非线性运算MSA校验节点更新公式即符号部分不变对幅值部分进行了近似操作。这个近似结果是偏大的。如前所述是一个减函数且为正值因为都是正数所以一定大于等于其中最大的那一项又因为也是减函数因此。特性分析硬件优势复杂的双曲函数被简化为“比较”和“符号异或”操作。这极大地降低了逻辑门数量使得高速并行实现成为可能。性能损失由于实际上总是大于等于真实的值MSA 倾向于高估校验节点输出的可靠性LLR幅度偏大。这种过估计会导致性能下降通常在高斯信道下比SPA差约0.5dB到1dB。2.4.3 修正最小和算法1、归一化最小和算法Scaled Min-Sum Algorithm, NMSA为了修正MSA的幅度过估计问题NMSA引入了一个归一化因子缩放因子通常。NMSA的校验节点更新公式修正为参数选择的最优值取决于码率、度分布和信噪比通常通过密度进化理论计算或蒙特卡洛仿真搜索得到。对于大多数规则码是常用值。优势NMSA保留了MSA的线性特性仅增加了一个乘法可以通过移位加法实现却能恢复大部分因近似丢掉的性能。它是目前5G NR、DVB-T2/S2和Wi-Fi标准中LDPC译码器的主流选择。2、偏移最小和算法Offset Min-Sum Algorithm, OMSAOMSA采用减法偏移来降低LLR幅度修正过估计。OMSA的校验节点更新公式修正为其中是非负偏移量。优势在硬件中减法比乘法更容易实现尤其对于定点数。此外OMSA在某些高信噪比场景下对错误平层的抑制效果略优于NMSA。三、Matlab代码3.1 SPA基础版本严格对照译码流程逻辑清晰但包含大量for循环效率较低适合用来理解算法。function [decoded_bits, iter] SPA_BP(H, llr, max_iterations) % 实现SPA-BP译码算法 % 输入: % H - 校验矩阵 % llr - 接收到的LLR值 % max_iterations - 最大迭代次数 % 输出: % decoded_bits - 译码后的比特序列 % 获取校验矩阵尺寸 [m, n] size(H); % 初始化变量节点到校验节点的消息 % v2c[i,j]表示第i个变量节点到第j个校验节点的消息 v2c zeros(n, m); % 初始化校验节点到变量节点的消息 % c2v[j,i]表示第j个校验节点到第i个变量节点的消息 c2v zeros(m, n); % 初始化变量节点的后验LLR值 posterior_llr llr; % 创建校验节点和变量节点的连接列表避免遍历整个矩阵 check_nodes cell(m, 1); variable_nodes cell(n, 1); for i 1:m check_nodes{i} find(H(i, :) 1); end for j 1:n variable_nodes{j} find(H(:, j) 1); end % BP迭代过程 for iter 1:max_iterations % 校验节点到变量节点的消息更新 for i 1:m connected_vars check_nodes{i}; for j_idx 1:length(connected_vars) j connected_vars(j_idx); % 校验节点i到变量节点j的消息更新使用tanh-product规则 prod_tanh 1; for k_idx 1:length(connected_vars) k connected_vars(k_idx); if k ~ j % 防止数值不稳定 v2c_clamped max(min(v2c(k, i), 30), -30); prod_tanh prod_tanh * tanh(v2c_clamped / 2); end end % 防止数值不稳定避免log(0) if abs(prod_tanh) 0.9999 prod_tanh sign(prod_tanh) * 0.9999; end c2v(i, j) 2 * atanh(prod_tanh); % 限制消息的数值范围防止数值溢出 c2v(i, j) max(min(c2v(i, j), 30), -30); end end % 变量节点到校验节点的消息更新 for j 1:n connected_checks variable_nodes{j}; for i_idx 1:length(connected_checks) i connected_checks(i_idx); % 变量节点j到校验节点i的消息 % 原始LLR 所有其他连接的校验节点到该变量节点的消息 v2c(j, i) llr(j); for k_idx 1:length(connected_checks) k connected_checks(k_idx); if k ~ i v2c(j, i) v2c(j, i) c2v(k, j); end end end end % 计算每个变量节点的后验LLR值 for j 1:n posterior_llr(j) llr(j); connected_checks variable_nodes{j}; for i_idx 1:length(connected_checks) i connected_checks(i_idx); posterior_llr(j) posterior_llr(j) c2v(i, j); end end % 硬判决 decoded_bits double(posterior_llr 0); % 校验是否满足所有校验方程 syndromes mod(H * decoded_bits, 2); if sum(syndromes) 0 % 所有校验方程都满足提前结束迭代 break; end end end优化版本利用了Matlab擅长矩阵运算的特点去掉繁琐的for循环运行速度可以提升10倍左右function [decoded_bits, iter] SPA_BP_Optimized(H, llr, max_iterations) % SPA_BP_Optimized 实现矢量化BP译码算法 [m, n] size(H); % 预处理将H转换为逻辑矩阵用于快速掩膜(Masking) H_mask logical(H); % 初始化消息矩阵 (全部使用 m x n 维度) % V2C: 变量节点到校验节点 (初始为LLR值仅在H_mask位置有效) V2C repmat(llr, m, 1) .* H_mask; % C2V: 校验节点到变量节点 C2V zeros(m, n); % 数值稳定性常数 SMALL_VAL 1e-15; % 防止log(0) LARGE_VAL 30; % LLR截断阈值防止溢出 for iter 1:max_iterations % % 1. 校验节点更新 (Horizontal Step) % % 目标:计算 2 * atanh( prod( tanh(V2C/2) \ i ) ) % 优化策略: 使用 φ 函数技巧将乘积转为求和并将符号和数值分离处理 % A. 处理符号 (Sign) Sign_V2C sign(V2C); Sign_V2C(V2C 0) 1; % 避免0符号问题 Sign_V2C(~H_mask) 1; % 非连接处设为1不影响乘积 % 计算每行的总奇偶校验位 Total_Sign prod(Sign_V2C, 2); % 排除自身的符号 (通过除法或异或这里用乘法因为是1/-1) C2V_Sign Total_Sign .* Sign_V2C; % B. 处理数值 (Magnitude) - 使用 -log(tanh(|x|/2)) 变换 Abs_V2C abs(V2C); % 避免数值计算错误 log(tanh(0))Inf对极小值进行裁剪 Abs_V2C(Abs_V2C SMALL_VAL) SMALL_VAL; % 核心变换函数 phi(x) Phi_V2C -log(tanh(Abs_V2C / 2)); Phi_V2C(~H_mask) 0; % 非连接处设为0不影响求和 % 行求和 Total_Phi sum(Phi_V2C, 2); % 排除自身贡献 Phi_External Total_Phi - Phi_V2C; % 逆变换回去 (函数形式相同) % C2V_Mag -log(tanh(Phi_External / 2)); % 同样对极小值进行裁剪 Phi_External(Phi_External SMALL_VAL) SMALL_VAL; C2V_Mag -log(tanh(Phi_External / 2)); % C. 组合符号和数值并进行范围截断 C2V C2V_Sign .* C2V_Mag; C2V max(min(C2V, LARGE_VAL), -LARGE_VAL); % 应用掩膜确保非连接处为0 C2V C2V .* H_mask; % % 2. 变量节点更新 (Vertical Step) % % 目标: LLR_post LLR_intrinsic sum(C2V) % V2C_new LLR_post - C2V_old % 计算后验 LLR (列求和) Sum_C2V sum(C2V, 1); Posterior_LLR llr Sum_C2V; % 硬判决 decoded_bits double(Posterior_LLR 0); % 校验是否满足 H * x 0 % 注意Matlab中 dense matrix * dense vector 很快 syndrome mod(H * decoded_bits, 2); if all(syndrome 0) return; end % 更新 V2C 消息 (总和 - 来自该校验节点的消息) V2C repmat(Posterior_LLR, m, 1) - C2V; % 再次截断并掩膜 V2C max(min(V2C, LARGE_VAL), -LARGE_VAL); V2C V2C .* H_mask; end end3.2 MSA、NMSA、OMSAMSA基础版本将校验节点更新的部分相应改为% 校验节点到变量节点的消息更新 for i 1:m connected_vars check_nodes{i}; for j_idx 1:length(connected_vars) j connected_vars(j_idx); % 校验节点i到变量节点j的消息更新使用Min-sum规则 min_val Inf; sign_prod 1; for k_idx 1:length(connected_vars) k connected_vars(k_idx); if k ~ j sign_prod sign_prod * sign(v2c(k, i)); min_val min(min_val, abs(v2c(k, i))); % 遍历找到最小值 end end c2v(i, j) sign_prod * min_val; % 限制消息的数值范围防止数值溢出 c2v(i, j) max(min(c2v(i, j), 20), -20); end endNMSA基础版本将校验节点更新的部分相应改为% 校验节点到变量节点的消息更新 for i 1:m connected_vars check_nodes{i}; for j_idx 1:length(connected_vars) j connected_vars(j_idx); % 校验节点i到变量节点j的消息更新使用Min-sum规则 min_val Inf; sign_prod 1; for k_idx 1:length(connected_vars) k connected_vars(k_idx); if k ~ j sign_prod sign_prod * sign(v2c(k, i)); min_val min(min_val, abs(v2c(k, i))); end end % 应用归一化因子alpha c2v(i, j) sign_prod * alpha * min_val; % 限制消息的数值范围防止数值溢出 c2v(i, j) max(min(c2v(i, j), 20), -20); end endOMSA基础版本将校验节点更新的部分相应改为% 校验节点到变量节点的消息更新 for i 1:m connected_vars check_nodes{i}; for j_idx 1:length(connected_vars) j connected_vars(j_idx); % 校验节点i到变量节点j的消息更新使用Mins_sum规则 min_val Inf; sign_prod 1; for k_idx 1:length(connected_vars) k connected_vars(k_idx); if k ~ j sign_prod sign_prod * sign(v2c(k, i)); min_val min(min_val, abs(v2c(k, i))); end end % 应用偏移操作 offset_val max(min_val - beta, 0); c2v(i, j) sign_prod * offset_val; % 限制消息的数值范围防止数值溢出 c2v(i, j) max(min(c2v(i, j), 20), -20); end end优化版本将以上四种算法整合了一下直接一步到位function [decoded_bits, iter] BP_Decoder(H, llr, max_iterations, algo_type, param) % 实现统一的 BP 译码 % 输入: % H: 校验矩阵 % llr: 初始对数似然比 % max_iterations: 最大迭代次数 % algo_type: SPA, MSA, NMSA, OMSA % param: 对于 NMSA 是缩放因子 alpha (如 0.75); 对于 OMSA 是偏移量 beta (如 0.15) if nargin 4 algo_type SPA; end if nargin 5 % 默认参数设置 if strcmp(algo_type, NMSA) param 0.75; elseif strcmp(algo_type, OMSA) param 0.15; else param 0; end end [m, n] size(H); H_mask logical(H); % 初始化消息矩阵 V2C repmat(llr, m, 1) .* H_mask; C2V zeros(m, n); % 常数定义 LARGE_VAL 30; % 截断阈值 SMALL_VAL 1e-15; % 仅用于 SPA防止 log(0)Inf for iter 1:max_iterations % % 1. 校验节点更新 (Horizontal Step) % % --- A. 符号处理 (所有算法通用) --- Sign_V2C sign(V2C); Sign_V2C(V2C 0) 1; Sign_V2C(~H_mask) 1; % 计算总奇偶校验 (行乘积) Total_Sign prod(Sign_V2C, 2); C2V_Sign Total_Sign .* Sign_V2C; % 排除自身 % --- B. 数值处理 (根据算法区分) --- Abs_V2C abs(V2C); if strcmp(algo_type, SPA) % SPA 算法 Abs_V2C(Abs_V2C SMALL_VAL) SMALL_VAL; Phi_V2C -log(tanh(Abs_V2C / 2)); Phi_V2C(~H_mask) 0; Total_Phi sum(Phi_V2C, 2); Phi_External Total_Phi - Phi_V2C; % 防止数值问题 Phi_External(Phi_External SMALL_VAL) SMALL_VAL; C2V_Mag -log(tanh(Phi_External / 2)); else % MSA / NMSA / OMSA % 1. 预处理将非连接处(0)设为无穷大避免影响求最小值 Abs_V2C_Temp Abs_V2C; Abs_V2C_Temp(~H_mask) inf; % 2. 寻找每一行的 最小值(min1) 和 次小值(min2) % [val, idx] min(...) [min1_val, min1_idx] min(Abs_V2C_Temp, [], 2); % 为了找次小值先将最小值的位置设为无穷大 Abs_V2C_Temp_NoMin1 Abs_V2C_Temp; % 使用线性索引快速替换 linear_indices sub2ind([m, n], (1:m), min1_idx); Abs_V2C_Temp_NoMin1(linear_indices) inf; [min2_val, ~] min(Abs_V2C_Temp_NoMin1, [], 2); % 3. 构建 C2V_Mag 矩阵 % 逻辑如果所在位置正好是最小值 min1 则输出 min2否则是 min1 % 先假设所有位置输出的都是全局最小值 C2V_Mag repmat(min1_val, 1, n); % 将原本是 min1 的位置用 min2 代替 C2V_Mag(linear_indices) min2_val; % 4. 针对变种进行修正 if strcmp(algo_type, NMSA) % Normalized MSA: 乘以归一化因子 alpha (0 alpha 1) alpha param; C2V_Mag C2V_Mag * alpha; elseif strcmp(algo_type, OMSA) % Offset MSA: 减去偏移量 beta, 并保证不小于0 beta param; C2V_Mag max(C2V_Mag - beta, 0); end end % --- C. 组合与截断 (所有算法通用) --- C2V C2V_Sign .* C2V_Mag; C2V max(min(C2V, LARGE_VAL), -LARGE_VAL); C2V C2V .* H_mask; % % 2. 变量节点更新 (Vertical Step) % Sum_C2V sum(C2V, 1); Posterior_LLR llr Sum_C2V; decoded_bits double(Posterior_LLR 0); % 校验 H*x^T 0 syndrome mod(H * decoded_bits, 2); if all(syndrome 0) return; end V2C repmat(Posterior_LLR, m, 1) - C2V; V2C max(min(V2C, LARGE_VAL), -LARGE_VAL); V2C V2C .* H_mask; end end