LDPC译码原理（公式推导）及其matlab代码实现（超详细）

博文更改记录

时间	版本号	更改内容
2022年4月12日 15:20	V1.0	新建博文
2022年4月20日 13:20	V1.1	完成初稿并发布该博文初稿内容包含相关背景概述及相关LDPC译码理论公式推导。
2022年4月21日 9:08	v1.2	补充2.4.2小节AWGN信道下的初始化的相关内容
2022年5月1日 14:25	v1.3	补充3.3节和3.4节关于归一化最小和译码算法和分层最小和译码算法等内容
2022年5月20日 11:26	V1.4	补充2.3.5小节关于改进的最小和译码算法增加LBP算法和NMSA收敛速度比较的Monte Carlo仿真matlab代码

一、背景概述

低密度奇偶校验码Low-Density Parity-Check Codes凭借其逼近Shannon限的纠错性能受到了信道编译码学者的广泛关注已成为DVB-S2、IEEE802.16e、CCSDS、5G等无线通信标准首选的信道编码方案。2001年S.Chung等人的研究结果表明码率1/2、码长10^7的非规则LDPC码在AWGN信道下采用置信传播算法进行迭代译码当错误概率为10 ^-6时距离Shnnon限仅相差0.0045dB。
LDPC码最迷人的地方在于译码算法本文主要关注基于置信传播的迭代译码算法博主在初学时在阅读论文和相关书籍过程中发现关于LDPC译码理论方面的讲解公式推导讲的不是很清楚需要阅读大量的书籍和文献才能对其有一个系统性的认识且关于LDPC译码matlab代码实现没有详尽的注释理解起来较为困难入门较难花费较多时间本文主要介绍基于置信传播的迭代译码算法及其改进算法对涉及的公式进行了详细的推导并附上自己的理解除此之外给出了matlab实现代码。因本文内容较多难免会出现错误还请各位读者批评指正希望这篇博文能对大家有所帮助和启发。

二、LDPC译码理论

2.1 LDPC码的表示方法

2.1.1LDPC码的矩阵表示

一个LDPC码$v$ 是一种$n$$k$线性分组码码长为$n$信息序列长度为$k$可由其校验矩阵H所唯一确定校验矩阵中1的数目远小于0的数目具有稀疏性。H的维数是$m\times n$每一行对应一个校验方程也称校验节点每一列对应码字的一位也称变量节点。每一行中非零元素的个数称为行重每一列中非零元素的个数称为列重。
下面是一个5 $\times$ 10的校验矩阵及其对应的校验方程
$$H=\left[\begin{array}{cccccccccc}1& 0& 0& 0& 0& 1& 1& 0& 0& 1\\ 0& 1& 1& 0& 0& 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 1& 0& 0& 0& 1& 1\\ 0& 0& 1& 1& 1& 0& 0& 1& 0& 0\\ 0& 1& 0& 1& 0& 0& 1& 1& 0& 0\end{array}\right]\to \{\begin{array}{c}{v}_0+v_5+v_6+v_9=0\\ v_1+v_2+v_5+v_8=0\\ v_0+v_4+v_8+v_9=0\\ v_2+v_3+v_4+v_7=0\\ v_1+v_3+v_6+v_7=0\end{array}$$

2.1.2 Tanner图表示

除了用校验矩阵表示LDPC码以外还可以用双向的图模型表示LDPC码其中Tanner图表示是比较方便的一种可以形象地刻画LDPC码的编译码特性。当H矩阵中对应位置不为0时在Tanner中有一条边将对应的校验节点和变量节点相连。校验矩阵的行重与列重与节点的度一致Tanner与校验矩阵一一对应。
下图所展示为5 $\times$ 10的校验矩阵对应的Tanner图

2.2符号说明

符号	含义
$v_i$	第 $i$个变量节点
$c_j$	第 $j$个校验节点
$r_{ji}^{(l)}(b)$	第$l$次迭代校验节点$j$传递给变量节点$i$的外信息,$b=0,1$
$q_{ij}^{(l)}(b)$	第$l$次迭代变量节点$i$传递给校验节点$j$的外信息,$b=0,1$
$C(i)$	与第$i$个变量节点相连的所有校验节点集合, C ( i ) = { i : h i j = 1 } C(i)=\{i:h_{ij}=1\} C(i)={i:hij​=1}
$V(j)$	与第$j$校验节点相连的所有变量节点集合, V ( j ) = { j : h i j = 1 } V(j)=\{j:h_{ij}=1\} V(j)={j:hij​=1}
$C\left(i\right)\backslash j$	除第$j$个校验节点外与第$i$个变量节点相连的其他校验节点集合, C ( i ) \ j = { k : h i k = 1 , k ≠ j } C\left( i \right)\backslash j=\{k:h_{ik}=1,k\ne j\} C(i)\j={k:hik​=1,k=j}
$V\left(j\right)\backslash i$	除第$i$个变量节点外与第$j$个校验节点相连的其他变量节点集合, V ( j ) \ i = { k : h k j = 1 , k ≠ i } V\left( j\right)\backslash i=\{k:h_{kj}=1,k\ne i\} V(j)\i={k:hkj​=1,k=i}
$P_i(b)$	接收端收到$y_i$后对应发送端码字$c_i=b$的后验概率,$b=0,1$
$q_i^{(l)}(b)$	第$l$次迭代第$i$个变量节点的后验概率信息,$b=0,1$

不同的文献对内信息和外信息的定义不一样有的文献把从信道接收来的信息称为外信息$externalmessages$把节点之间传递的信息称为内信息$internalmessages$;有的文献把从信道接收来的信息称为后验概率把节点之间传递的信息称为外信息本博文定义遵照后者。

2.3LDPC译码算法

LDPC码的译码方法可以分为两大类基于硬判决的译码和基于软判决的译码。基于硬判决的译码运算量小比较实用而软判决译码采用了后验概率信息并通过迭代运算使得LDPC码的性能得以逼近香农限。本博文主要关注基于软判决的译码及其改进算法。

2.3.1引理1及其证明

在学习LDPC译码相关的知识时都会涉及到引理1该引理关系具体LDPC译码算法但很少有文献给出关于引理的证明本博文关于引理的证明主要参考了博主stand_by_me123的博文。¹
引理1对于一个长为$m$的相互独立的二进制序列$a=[a_1,a_2,...,a_m]$第$h$个比特为1的概率$P(a_h)=p_h$那么二进制序列$a$中1的个数为偶数的概率为$\frac{1}{2}+\frac{1}{2}\prod _{h=1}^m(1-2p_h)$1的个数为奇数的概率为$\frac{1}{2}-\frac{1}{2}\prod _{h=1}^m(1-2p_h)$。
引理1的证明构造函数+数学归纳法
首先构造关于$t$的$m$次多项式$f(t)=\prod _{h=1}^m(1-p_h+p_{h}t)$通过数学归纳法证明$t^i$的系数为其序列中包含$i$个1的概率。
1当$m=1$$f(t)=1-p_h+p_{h}t$,显然$t$的一次幂及零次幂即为包含1个1和0个1的概率
2假设当$m=k$时$t^h$的系数$B_h$为$k$次多项式$f(t)=\prod _{h=1}^k(1-p_h+p_{h}t)=\sum _{h=0}^k{B}_h{t}^h$中包含$h$个1个概率
3当$m=k+1$时$$f(t)=\prod _{h=1}^{k+1}(1-p_h+p_{h}t)=(1-p_{k+1}+p_{k+1}t)\prod _{h=1}^k(1-p_h+p_{h}t)=(1-p_{k+1}+p_{k+1}t)\sum _{h=0}^k{B}_h{t}^h$$
由上式可求得$t^i$的系数为$(1-p_{k+1})B_i+p_{k+1}{B}_{i-1}$即为$k+1$次多项式中包含$i$个1的概率。
对于函数$g(t)=\prod _{h=1}^m(1-p_h-p_{h}t)$与$f(t)$的区别在于函数$g(t)$的奇次幂系数为负值。故可通过$f(t)+g(t)$来消去幂次为奇数的项并令$t=1$。
故可得到二进制序列$a$中1的个数为偶数的概率为$p_{even}=\frac{1+\prod _{h=1}^m(1-2p_h)}{2}$1的个数为奇数的概率为$p_{odd}=1-p_{even}=\frac{1-\prod _{h=1}^m(1-2p_h)}{2}$。

2.3.2和积译码算法概率域译码算法

1计算发送端发送的比特$v_i=1$或0的初始后验概率设定变量节点传递给校验节点的初始消息为发送端发送比特的初始后验概率
$$\{\begin{array}{c}{q}_{ij}^{(0)}(1)=P_i(1)=p(v_i=1|y_i)\\ q_{ij}^{(0)}(0)=P_i(0)=p(v_i=0|y_i)\end{array}$$
2计算校验节点传递给变量的外信息
$$\{\begin{array}{c}{r}_{ji}^{(l)}(0)=\frac{1}{2}+\frac{1}{2}\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(1-2q_{i^{{}^{\prime }}j}^{(l-1)}(1))\\ r_{ji}^{(l)}(1)=1-r_{ji}(0)=\frac{1}{2}-\frac{1}{2}\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(1-2q_{i^{{}^{\prime }}j}^{(l-1)}(1))\end{array}$$
理解首先理解校验节点传递给变量节点外信息的含义对于校验节点$c_j$即代表一个校验方程与其连接的一个变量节点$v_i$的值确定且除$v_i$外的其他变量节点为0或1的概率已知的条件下第$j$个校验方程满足的概率。
若变量节点$v_i$为0则第$j$个校验方程满足的条件为其他变量节点中有偶数个1。校验方程的计算法则为模2运算
若变量节点$v_i$为1则第$j$个校验方程满足的条件为其他变量节点中有奇数个1。
一个序列中有偶数个1或奇数1的概率计算公式由引理1给出。
3计算变量节点传递给校验节点的外信息
$$\{\begin{array}{c}{q}_{ij}^{(l)}(0)=K_{ij}{q}_{ij}^{(0)}(0)\prod _{j^{{}^{\prime }}\in C(i)\backslash j}{r}_{j^{{}^{\prime }}i}^{(l)}(0)\\ q_{ij}^{(l)}(1)=K_{ij}{q}_{ij}^{(0)}(1)\prod _{j^{{}^{\prime }}\in C(i)\backslash j}{r}_{j^{{}^{\prime }}i}^{(l)}(1)\end{array}$$
理解对于变量节点传递给校验节点外信息的含义当变量节点$v_i$连接的一个校验节点$c_j$的值确定即第$j$个校验方程确定为0或1与$v_i$相连的其他校验节点所在的校验方程能得到满足的概率已知此时变量节点$v_i=1$或$0$的概率。
$v_i=1$的概率包括初始为1的概率、为1时使得其他校验方程满足时的概率这些概率值相乘得到了$v_i=1$的概率$v_i=0$同理。
4译码判决计算经过迭代后变量节点的后验概率
$$\{\begin{array}{c}{q_{ij}}^{(l)}(0)=K_{ij}{q_{ij}}^{(0)}(0){\prod _{j^{{}^{\prime }}\in C(i)}{r}_{j^{{}^{\prime }}i}}^{(l)}(0)\\ {q_{ij}}^{(l)}(1)=K_{ij}{q_{ij}}^{(0)}(1){\prod _{j^{{}^{\prime }}\in C(i)}{r}_{j^{{}^{\prime }}i}}^{(l)}(1)\end{array}$$
0的概率大时将$v_i$判决为0否则判决为1。

2.3.3对数域BP译码算法LLR BP

概率域上的BP算法包含了大量的连乘运算硬件实现时具有较高的计算复杂度资源消耗大。将概率信息用似然比表示就得到对数域上的BP算法大量的乘法运算可以转化为加法运算。
1信道初始消息
$$L(P_i)=\ln \frac{P_i(0)}{P_i(1)}=\ln \frac{p(v_i=0|y_i)}{p(v_i=1|y_i)}$$
2校验节点传递给变量节点的外信息
$$r_{ji}=\ln \frac{r_{ji}(0)}{r_{ji}(1)}=\ln \frac{1+\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(1-2q_{i^{{}^{\prime }}j}(1))}{1-\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(1-2q_{i^{{}^{\prime }}j}(1))}$$
注意到${\tanh }^{-1}=\frac{1}{2}\ln \frac{1+x}{1-x}$$1-2q_{i^{{}^{\prime }}j}(1)=q_{i^{{}^{\prime }}j}(0)-q_{i^{{}^{\prime }}j}(1)$且$q_{i^{{}^{\prime }}j}(0)+q_{i^{{}^{\prime }}j}(1)=1$$q_{ij}=\ln \frac{q_{ij}(0)}{q_{ij}(1)}$$\tanh (x)=\frac{e^x-e^{-x}}{e^x+e^{-x}}$故可利用这几个公式对$r_{ji}$进行化简。
希望读者可以静下心仔细推不要惧怕公式。
$$r_{ji}=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(1-2q_{i^{{}^{\prime }}j}(1)))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(q_{i^{{}^{\prime }}j}(0)-q_{i^{{}^{\prime }}j}(1)))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(\frac{q_{i^{{}^{\prime }}j}(0)-q_{i^{{}^{\prime }}j}(1)}{q_{i^{{}^{\prime }}j}(0)+q_{i^{{}^{\prime }}j}(1)}))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(\frac{\frac{q_{i^{{}^{\prime }}j}(0)}{q_{i^{{}^{\prime }}j}(1)}-1}{\frac{q_{i^{{}^{\prime }}j}(0)}{q_{i^{{}^{\prime }}j}(1)}+1}))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(\frac{e^{q_{i^{{}^{\prime }}j}}-1}{e^{q_{i^{{}^{\prime }}j}}+1}))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(\frac{e^{\frac{q_{i^{{}^{\prime }}j}}{2}}(e^{\frac{q_{i^{{}^{\prime }}j}}{2}}-e^{-\frac{q_{i^{{}^{\prime }}j}}{2}})}{e^{\frac{q_{i^{{}^{\prime }}j}}{2}}(e^{\frac{q_{i^{{}^{\prime }}j}}{2}}+e^{-\frac{q_{i^{{}^{\prime }}j}}{2}})}))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(\frac{e^{\frac{q_{i^{{}^{\prime }}j}}{2}}-e^{-\frac{q_{i^{{}^{\prime }}j}}{2}}}{e^{\frac{q_{i^{{}^{\prime }}j}}{2}}+e^{-\frac{q_{i^{{}^{\prime }}j}}{2}}}))$$
$$=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}tanh\frac{q_{i^{{}^{\prime }}j}}{2})$$
故校验节点传递给变量节点的信息可以化简为
$r_{ji}=2tan{h}^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}tanh\frac{q_{i^{{}^{\prime }}j}}{2})$
3变量节点传递给校验节点的外信息
$$q_{ij}=\ln \frac{P_i(0)\prod _{j^{{}^{\prime }}\in C(i)\backslash j}{r}_{j^{{}^{\prime }}i}(0)}{P_i\left(1\right)\prod _{j^{{}^{\prime }}\in C(i)\backslash j}{r}_{j^{{}^{\prime }}i}(1)}=L(P_i)+\sum _{j^{{}^{\prime }}\in C(i)\backslash j}{r}_{j^{{}^{\prime }}i}$$
4译码判决经过迭代后变量节点的后验概率相应的修改为
$$q_{ij}=\ln \frac{P_i(0)\prod _{j^{{}^{\prime }}\in C(i)}{r}_{j^{{}^{\prime }}i}(0)}{P_i\left(1\right)\prod _{j^{{}^{\prime }}\in C(i)}{r}_{j^{{}^{\prime }}i}(1)}=L(P_i)+\sum _{j^{{}^{\prime }}\in C(i)}{r}_{j^{{}^{\prime }}i}$$
故后验概率大于0将$v_i$判决为0反之判决为1。对数值大于0说明真数大于1及分子对应的$v_i$为0的概率大于分母对应$v_i$为1的概率故将$v_i$判决为0

2.3.4最小和译码算法Min-Sum BP

对数域上的BP算法大大减少了乘法的运算次数但是引入了$tanh$函数在硬件实现时非线性运算需要消耗较多的硬件资源从实现的角度会采用查表方式性能优异的LDPC码通常码长较长计算量会随码长的增加而增大查表的引入让电路变得复杂带来较高的译码延迟和译码功耗。
为了降低计算复杂度简化BP译码算法提出了最小和译码算法最小和算法是对校验节点传递给变量节点外信息的改进将复杂的$tanh$函数和$tan{h}^{-1}$函数转化为基本的符号值计算和数值比较运算大大降低了算法的复杂度。
注意到$tanh$函数和$tan{h}^{-1}$函数是奇函数满足
$$tanh(x)=sgn(x)\cdot tanh(\left|x\right|)$$
$$tan{h}^{-1}(x)=sgn(x)\cdot tan{h}^{-1}(|x|)$$
$tanh(x)$函数的图像如下图所示可以看出$tanh(x)$为偶函数在$[0,1]$上单调递增。

故校验节点传递给变量节点的外信息可改写为
$$r_{ji}=2{\tanh }^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}\tanh \frac{q_{i^{{}^{\prime }}j}}{2})=2{\tanh }^{-1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}(sgn(\frac{q_{i^{{}^{\prime }}j}}{2})\cdot \tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|)))=2\prod _{i^{{}^{\prime }}\in V(j)\backslash i}sgn(\frac{q_{i^{{}^{\prime }}j}}{2})\cdot {\tanh }^{\text{-}1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}\tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|))$$
又
$$2{\tanh }^{\text{-}1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}\tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|))\approx 2{\tanh }^{-1}(\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|)))=2{\tanh }^{-1}(\tanh (\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|)))=\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)$$
故校验节点传递给变量节点的外信息最终可化简为
$$r_{ji}=\prod _{i^{{}^{\prime }}\in V(j)\backslash i}sgn(\frac{q_{i^{{}^{\prime }}j}}{2})\cdot \underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)=\prod _{i^{{}^{\prime }}\in V(j)\backslash i}sgn(q_{i^{{}^{\prime }}j})\cdot \underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)$$
在对数域BP算法中变量节点消息初始化时需要对信道参数进行估计才能获得参数${\delta }^2$但在最小和算法中数值计算转变成变量节点概率绝对值的比较运算且信道参数是其共同的系数可以省略不影响整个译码过程。在具体实现时不需要再额外添加信道参数估计模块进一步简化了电路提高了译码的时效性。
由图像可知$tanh(x)$函数小于1的相乘之后的结果肯定是越来越小所以只用最小值来代替连乘的结果但这样最小值肯定是大于连乘之后的结果的这时需要引入一个修正因子来补偿译码损失的性能。进而引出了归一化最小和译码算法和偏移最小和译码算法通过引入一个小于1的乘性因子或者减去一个值来补偿损失的译码性能。

2.3.5改进的最小和译码算法NMSA and OMSA

最小和译码算法将LLR BP的非线性乘法运算简化为符号乘法和最小值运算这种简化运算带来的精度的损失影响译码性能。注意到
$$2{\tanh }^{\text{-}1}(\prod _{i^{{}^{\prime }}\in V(j)\backslash i}\tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|))<2{\tanh }^{-1}(\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\tanh (\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|)))=2{\tanh }^{-1}(\tanh (\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|\frac{q_{i^{{}^{\prime }}j}}{2}\right|)))=\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)$$
即Min-Sum算法的校验节点更新值要大于LLR BP。为平衡两种译码算法校验节点更新值幅度上的差异两种改进最小和译码被提出归一化最小和算法NMSANormalized Min-Sum Algorithm和偏移最小和算法OMSAOffoset Min-Sum Algorithm算法即乘一个小于1的常数$\alpha$或者减去一个正值$\beta$。
NMSA
$$r_{ji}=\alpha \cdot \prod _{i^{{}^{\prime }}\in V(j)\backslash i}sgn(q_{i^{{}^{\prime }}j})\cdot \underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)$$
$\alpha$是一个乘性缩小因子也叫归一化因子$0\le \alpha \le 1$。在实际应用中$\alpha$通常在0.6~0.9之间选取一个固定值。NMSA虽然多了一次乘法运算对计算复杂度影响不大但是却可以很好的改善译码性能。
OMSA
$$r_{ji}=max(\underset{i^{{}^{\prime }}\in V(j)\backslash i}{\min }(\left|q_{i^{{}^{\prime }}j}\right|)-\beta ,0)\cdot \prod _{i^{{}^{\prime }}\in V(j)\backslash i}sgn(q_{i^{{}^{\prime }}j})$$
$\beta$为偏移因子其值大于0。同归一化因子$\alpha$一样最佳$\beta$值的选取受到LDPC码码长、码率、信噪比以及迭代次数的影响但在实际应用时$\beta$通常为常数。
由此可以看到NMSA算法和OMSA算法尽管处理方式不一样但核心思想都是补偿校验节点信息过大估计的问题。在实际应用中通过仿真确定最佳修正因子$\alpha$和$\beta$对于系统性能的提升是至关重要的。

三、LDPC译码算法的matlab实现

3.1BPSK调制体制下关于SNR和$E_b/N_0$与$\delta$之间的关系

1SNR在连续和离散系统中计算
假设连续时间高斯信道Signal Power = PNoise PSD噪声功率密度=$N_0/2$单边谱密度为$N_0$双边谱密度为$N_0/2$Bandwith = 2W。
可以推出$NoisePower=NoisePSD\times Bandwidth=N_{0}W$
$$SNR=P/N_{0}W$$
假设离散时间信道每个符号的能量$E_s=PT=P/2W$信号的能量乘以其持续时间噪声方差${\delta }^2$。
$$SNR=E_S/{\delta }^2=P/2W{\delta }^2$$
2计算BPSK的$SNR$
对于BPSK有$E_s=[(-1{)}^2+(+1{)}^2]/2=1$信道的高斯噪声Noise Power = ${\delta }^2$。
所以$SNR=1/{\delta }^2$。
编码后每个bit的能量$E_b=E_s/R=n{E}_s/k$$R=k/n$码率=待编码的bit数/编码后的bit数。能量一定需要增加因为编码之后的信息增加了冗余相当于多个bit只代表一个bit的信息。
所以可以推出$E_s=R\cdot E_b$。
又因为$N_0/2={\delta }^2$所以可以推出$SNR=E_s/{\delta }^2=R{E}_b/(N_0/2)=2R(E_b/N_0)$。
故$E_b/N_0=SNR/2R$对于BPSK调制来说$E_b/N_0=1/(2R{\delta }^2)$。

3.2AWGN信道下的初始化

下面介绍一下常用的AWGN信道下BPSK调制的LDPC码译码的消息初始化。²
在通信系统中采用二进制数字调制码字$c_i$映射为符号$x_i=(-1{)}^{c_i},i=1,2,...,n$接收符号$y_i=x_i+n_i$各个$n_i$为统计独立同分布的高斯随机变量双边功率谱密度为$N_0/2$可知$y_i$是均值为1方差为${\delta }^2$的高斯变量。在信源等概率分布时$p(x_i=+1)=p(x_i=-1)=\frac{1}{2}$此时概率BP译码的初始消息为
$$q_{ij}^{(0)}(1)=p(c_i=1|y_i)=p(x_i=-1|y_i)=\frac{1}{1+e^{2y_i/{\delta }^2}}$$
$$q_{ij}^{(0)}(0)=p(c_i=0|y_i)=p(x_i=+1|y_i)=\frac{1}{1+e^{-2y_i/{\delta }^2}}$$
证明
在AWGN信道中采用BPSK调制当功率谱密度为$N_0/2={\delta }^2$条件概率分布函数为
$$p(y_i|x_i=+1)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(y_i-1{)}^2}{2{\delta }^2}}$$
$$p(y_i|x_i=-1)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(y_i+1{)}^2}{2{\delta }^2}}$$
由贝叶斯公式可得
$$p(x_i=+1|y_i)=\frac{p(x_i=+1,y_i)}{p(y_i)}=\frac{p(y_i|x_i=+1)p(x_i=+1)}{p(y_i)}=\frac{p(y_i|x_i=+1)p(x_i=+1)}{p(y_i|x_i=+1)p(x_i=+1)+p(y_i|x_i=-1)p(x_i=-1)}$$
因为$p(x_i=+1)=p(x_i=-1)=\frac{1}{2}$所以有
$$p(x_i=+1|y_i)=\frac{\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(y_i-1{)}^2}{2{\delta }^2}}}{\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(y_i-1{)}^2}{2{\delta }^2}}+\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(y_i+1{)}^2}{2{\delta }^2}}}=\frac{1}{1+e^{\frac{-2y_i}{{\delta }^2}}}$$
$p(x_i=-1|y_i)$同理。

由此LLR BP译码的初始消息为
$$L(P_i)=\ln \frac{P_i(0)}{P_i(1)}=\ln \frac{p(v_i=0|y_i)}{p(v_i=1|y_i)}=\ln (\frac{1+e^{\frac{2y}{{\delta }^2}}}{1+e^{\frac{-2y}{{\delta }^2}}})$$
$$=ln(\frac{e^{\frac{y}{{\delta }^2}}(e^{\frac{y}{{\delta }^2}}+e^{\frac{-y}{{\delta }^2}})}{e^{\frac{-y}{{\delta }^2}}(e^{\frac{y}{{\delta }^2}}+e^{\frac{-y}{{\delta }^2}})})=ln(e^{\frac{2y}{{\delta }^2}})=\frac{2y}{{\delta }^2}$$

进行Monte Carlo仿真时假设信道传送符号的平均能量为$E_b$信道编码码率为$r$采用匹配滤波器的最佳接收解调器输入的信噪比$E_b/N_0$与噪声方差的关系为
$${\delta }^2=\frac{1}{2r(E_b/N_0)}$$

3.3基于归一化最小和译码算法Normalized min-sum algorithm,NMSA的matlab实现

该matlab代码主要参考了博主cea_wind的博客³。

function [ iter,decoderData ] = ldpcdecoderminsum(H,HRowNum,HColNum,receiveSignal,MAX_ITER_NUM,NORM_FACTOR)
%LDPCC decode algorithm , based on the Min-Sum algorithm
% H: check matrix
% HRowNum,HColNum: index generate from check matrix
% receiveSignal: received soft information from channel
% MAT_INTER_NUM: maximum iterations
% HRowNum,HColNum generating by the following codes
% 	[r_mark,c_mark] = find(H~=0);
% 	HColNum = sum(H);
% 	HRowNum = cell(1,size(H,1));
% 	for rowH = 1:size(H,1)
%		 HRowNum{rowH} = find(r_mark==rowH);
% 	end

[~,N] = size(H);    

vl = receiveSignal;
decoderData = zeros(1,N);

uml = zeros(1,sum(HColNum));
vml = uml;
ColStart = 1;
for L=1:length(HColNum)
    vml(ColStart:ColStart+HColNum(L)-1) = vl(L);
    ColStart = ColStart+HColNum(L);
end

for iter=1:MAX_ITER_NUM
    %check nodes information process
    for L_r=1:length(HRowNum)
        L_col = HRowNum{L_r};
        vmltemp = vml(L_col);
        vmlMark = ones(size(vmltemp));
        vmlMark(vmltemp<0) = -1;
        vmlMark = prod(vmlMark);
        minvml = sort(abs(vmltemp));
        for L_col_i = 1:length(L_col)
            if minvml(1)==abs(vmltemp(L_col_i))
                if vmltemp(L_col_i)<0
                    vmltemp(L_col_i) = -vmlMark*minvml(2);
                else
                    vmltemp(L_col_i) = vmlMark*minvml(2);
                end
            else
                if vmltemp(L_col_i)<0
                    vmltemp(L_col_i) = -vmlMark*minvml(1);
                else
                    vmltemp(L_col_i) = vmlMark*minvml(1);
                end
            end
        end
        uml(L_col) = NORM_FACTOR*vmltemp;
    end
    %variable nodes information process
    ColStart = 1;
    qn0_1 = ones(1,N);
    for L=1:length(HColNum)
        umltemp = uml(ColStart:ColStart+HColNum(L)-1);
        temp = sum(umltemp);
        qn0_1(L) = temp + vl(L); 
        umltemp = temp - umltemp;
        vml(ColStart:ColStart+HColNum(L)-1) = umltemp + vl(L);
        
        ColStart = ColStart+HColNum(L);
    end
    % decision decoding
    decoderData(qn0_1>=0) = 0;
    decoderData(qn0_1<0) = 1;
    if(mod(H*decoderData',2)==0)
        break;
    end
end

Monte Carlo仿真matlab代码本博文matlab代码在cea_wind的代码基础上有所修改。博主cea_wind将他的matlab代码上传到了GitHub上读者可以通过文末的超链接自行下载学习⁴。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%.M file for the simulation of BER performence.
%Create file date:2022-3-11 10:00;
%Designer:;
%H:parity check matrics which has M rows and N columns
%Qvc:
%Rcv:
%Qv:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear everything;
clear all;
close all;
clc;
tic;

%set the simulation parameter
EbN0_dB = 1:0.2:4;
Frames_Num =2000;          
MAX_Iter_Num =5;

%%  NMSA:'r-'
bitError_NMSA = zeros(1,length(EbN0_dB));
BER_NMSA = zeros(1,length(EbN0_dB));
iterNumTotal_NMSA = zeros(1,length(EbN0_dB));

INFO_LENGTH = 4096;
RATE = 4/5;
SIZE_M = 512;
NORM_FACTOR = 0.75;

H = ccsdscheckmatrix(SIZE_M,RATE);
G = ccsdsgeneratematrix(H,SIZE_M,RATE);
% %[H]=genH(128,256);
[M,N] = size(H);

[r_mark,c_mark] = find(H~=0);
HColNum = sum(H);
HRowNum = cell(1,size(H,1));
for rowH = 1:size(H,1)
    HRowNum{rowH} = find(r_mark==rowH);
end

for nEbN0 = 1:length(EbN0_dB) 
   
    fprintf('%2.2f',EbN0_dB(nEbN0));
    fprintf('dB,');
    fprintf('m = %d,codenum = ',MAX_Iter_Num);
    
    SNR_per_bit = 10^(EbN0_dB(nEbN0)/10); 
    N0 = 1/(SNR_per_bit*RATE);
    sigma = sqrt(N0/2);
    
    for nF = 1:Frames_Num
       
        fprintf('%d,',nF);
        
        %encode
        message = randi([0 1],1,INFO_LENGTH);
        encodeData = mod(message*G,2);
        transmitSignal = 1 - 2 * encodeData;            
        
        %AWGN Channel
        
        receiveSignal = transmitSignal + sigma*randn(1,length(transmitSignal));
        
%%decode
        [iter_Num_NMSA,recoverData_NMSA] =ldpcdecoderminsum(H,HRowNum,HColNum,receiveSignal,MAX_Iter_Num,NORM_FACTOR);
        
        %biterror,BER,iterNumToal caculation;  
        bitError_NMSA(nEbN0) = bitError_NMSA(nEbN0) + sum(abs(message - recoverData_NMSA(1:length(message))));
        iterNumTotal_NMSA(nEbN0) = iterNumTotal_NMSA(nEbN0) + iter_Num_NMSA;
        
        if (nF == Frames_Num)
            BER_NMSA(nEbN0)= bitError_NMSA(nEbN0)/(nF*length(message));
          
            break;
        end       
    end  %end for 'nF'  
    fprintf('\n');
end  % end for 'nEbN0'
semilogy(EbN0_dB,BER_NMSA,'*r-');
xlabel('Eb/N0(dB)'); ylabel('BER');
grid on;
toc;

下图所展示的为matlab仿真误结果总帧数为2000帧为获得更精确误码率曲线读者可自行修改参与仿真总帧数并可尝试不等间隔$E_b/N_0$

误码率曲线

3.4分层归一化最小和LDPC译码算法

从BP和NMSA译码算法变量节点更新过程可以看出总是需要更新完所有校验节点再更新所有变量节点这种更新策略称为泛洪Flooding调度机制使得本次迭代过程中产生的新信息只能在下次迭代中使用。⁵
LDPC码译码算法的收敛速度也会影响译码复杂度收敛速度越快所需迭代次数越少而迭代次数的减少既可以降低译码延时又可以降低总计算量。⁶
E.Sharon,J.Zhang等人提出了LBP⁷ 和SBP⁸ 静态策略串行译码算法该类算法不同于传统的BP算法的泛洪调度机制而是采用固定的顺序对节点消息进行更新每次迭代中校验节点的更新可利用本次迭代已更新的信息加快了译码收敛速度。

3.4.1LBP译码算法

LBP译码算法以校验节点为单位按校验矩阵行的顺序进行节点更新其算法流程如下(该小节图片取自参考文献6)结合3.4.2小节的matlab代码会加深对Flooding调度和分层串行调度算法的理解。

LBP译码算法消息传递过程如下图所示

3.4.2 分层最小和译码算法的matlab实现

function [n,recoverData] = LBP(receiveSignal, H,MAX_Iter_Num,NORM_FACTOR)


[M,N]=size(H);

Qv = receiveSignal;
Rcv = zeros(M,N);
Qvc = zeros(1,N);
recoverData = zeros(1,N);
for n = 1:MAX_Iter_Num
   for i = 1:M
       col = find(H(i,:));
       for k = 1:length(col)
           Qvc(1,col(k)) = Qv(col(k)) - Rcv(i,col(k));
       end    
       
       alpha = sign(Qvc);
       beta = abs(Qvc);
       signS = 1;
       min1 = 100000;
       min2 = 100000;
       index = 10000;
       
       for k = 1:length(col)
          signS = alpha(col(k))*signS;
          if beta(col(k)) < min1
             min2 = min1;
             min1 = beta(col(k));
             index = col(k);
          end
          if ((beta(col(k)) > min1) && (beta(col(k))<min2))
              min2 = beta(col(k));
          end
       end
       
       for k = 1:length(col)   %error 0315 k = 1:length(col(k))
          if col(k) == index
              Rcv(i,col(k)) = NORM_FACTOR*signS*alpha(col(k))*min2;
          else
              Rcv(i,col(k)) = NORM_FACTOR*signS*alpha(col(k))*min1;
          end
          Qv(col(k)) = Qvc(1,col(k)) + Rcv(i,col(k));
       end     
   end              %end for i    
   for k = 1:N
       if Qv(k) < 0
           recoverData(k) = 1;
       else
           recoverData(k) = 0;
       end
   end
   if rem(H*recoverData',2) == 0
       break;
   end
end                 %end for n 

Monte Carlo仿真matlab代码代码中ccsdscheckmatrix函数和ccsdsgeneratematrix函数为博主cea_wind在GitHub上传的代码读者可以通过文末的超链接自行下载

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%.M file for the simulation of BER performence.
%Create file date:2022-3-11 10:00;
%Designer:;
%H:parity check matrics which has M rows and N columns
%Qvc:
%Rcv:
%Qv:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear everything;
clear all;
close all;
clc;
tic;

%set the simulation parameter
EbN0_dB = 1:0.2:4;
Frames_Num =2000;          
MAX_Iter_Num =5;

%%  NMSA:'*r-'
bitError_NMSA = zeros(1,length(EbN0_dB));
BER_NMSA = zeros(1,length(EbN0_dB));
iterNumTotal_NMSA = zeros(1,length(EbN0_dB));

%%LBP:'+b-'
bitError_LBP = zeros(1,length(EbN0_dB));
BER_LBP = zeros(1,length(EbN0_dB));
iterNumTotal_LBP = zeros(1,length(EbN0_dB));

INFO_LENGTH = 4096;
RATE = 4/5;
SIZE_M = 512;
NORM_FACTOR = 0.75;

H = ccsdscheckmatrix(SIZE_M,RATE);
G = ccsdsgeneratematrix(H,SIZE_M,RATE);
% %[H]=genH(128,256);
[M,N] = size(H);

[r_mark,c_mark] = find(H~=0);
HColNum = sum(H);
HRowNum = cell(1,size(H,1));
for rowH = 1:size(H,1)
    HRowNum{rowH} = find(r_mark==rowH);
end

for nEbN0 = 1:length(EbN0_dB) 
   
    fprintf('%2.2f',EbN0_dB(nEbN0));
    fprintf('dB,');
    fprintf('m = %d,codenum = ',MAX_Iter_Num);
    
    SNR_per_bit = 10^(EbN0_dB(nEbN0)/10); 
    N0 = 1/(SNR_per_bit*RATE);
    sigma = sqrt(N0/2);
    
    for nF = 1:Frames_Num
       
        fprintf('%d,',nF);
        
        %encode
        message = randi([0 1],1,INFO_LENGTH);
        encodeData = mod(message*G,2);
        transmitSignal = 1 - 2 * encodeData;            
        
        %AWGN Channel
        
        receiveSignal = transmitSignal + sigma*randn(1,length(transmitSignal));
        
%%decode
        [iter_Num_NMSA,recoverData_NMSA] =ldpcdecoderminsum(H,HRowNum,HColNum,receiveSignal,MAX_Iter_Num,NORM_FACTOR);
        [iter_Num_LBP,recoverData_LBP] =LBP(receiveSignal,H,MAX_Iter_Num,NORM_FACTOR);
        
%biterror,BER,iterNumToal caculation;  
        bitError_NMSA(nEbN0) = bitError_NMSA(nEbN0) + sum(abs(message - recoverData_NMSA(1:length(message))));
        iterNumTotal_NMSA(nEbN0) = iterNumTotal_NMSA(nEbN0) + iter_Num_NMSA;
        bitError_LBP(nEbN0) = bitError_LBP(nEbN0) + sum(abs(message - recoverData_LBP(1:length(message))));
        iterNumTotal_LBP(nEbN0) = iterNumTotal_LBP(nEbN0) + iter_Num_LBP;
        if (nF == Frames_Num)
            BER_NMSA(nEbN0)= bitError_NMSA(nEbN0)/(nF*length(message));
          BER_LBP(nEbN0)= bitError_LBP(nEbN0)/(nF*length(message));
            break;
        end       
    end  %end for 'nF'  
    fprintf('\n');
end  % end for 'nEbN0'
semilogy(EbN0_dB,BER_NMSA,'*r-',EbN0_dB,BER_LBP,'+b-');
xlabel('Eb/N0(dB)'); ylabel('BER');
grid on;
toc;

3.4.3 LBP译码算法和NMSA译码算法收敛速度对比

仿真参数
%set the simulation parameter
EbN0_dB = 1:0.2:4;
Frames_Num =2000;
MAX_Iter_Num =5;
仿真结果

可以看出LBP具有更快的译码收敛速度。

四、LDPC译码的调度算法

---

1. LDPC码译码原理——概率公式推导 ↩︎

2. 袁东风张海刚等. LDPC码理论与应用[M]. 北京人民邮电出版社2008. ↩︎

3. CCSDS标准的LDPC编译码仿真 ↩︎

4. 博主cea_wind上传的CCSDS标准LDPC编译码仿真matlab代码 ↩︎

5. 王冰冰. 空间通信中LDPC译码算法研究与译码器设计[D].中国科学院大学(中国科学院国家空间科学中心),2021.DOI:10.27562/d.cnki.gkyyz.2021.000036. ↩︎

6. 康婧. 星地高速数传LDPC码编译码算法及高效实现技术研究[D].中国科学院大学(中国科学院国家空间科学中心),2021.DOI:10.27562/d.cnki.gkyyz.2021.000010. ↩︎

7. E. Sharon, S. Litsyn, J. Goldberger. An Efficient Message-Passing Schedule for LDPC Decoding[C]//23rd IEEE Convertion of Electrical and Electronics Engineer, Tel-Aviv, Israel, 2004: 223-226. ↩︎

8. J. Zhang, M. Fossorier. Shuffled Iterative Decoding[J]. IEEE Transactions on Communications, 2005, 53(2): 209-213. ↩︎