基于度分解和交替乘子法的LDPC码二次规划译码方法与流程

本发明属于通信技术领域，特别涉及一种对低密度奇偶校验ldpc码的译码方法，可用于磁存储、光纤通信和卫星数字视频。

背景技术

低密度奇偶校验码ldpc是目前一种能够逼近香农信道容量限的最佳编码方案之一，受到国内外学者研究的关注，并被广泛应用在各种通信领域中。通常采用置信传播bp算法对ldpc码进行译码。由于bp算法会受到码字结构图中存在的许多有害特性的影响，包括伪码字、瞬子、陷阱集合等，因而在高信噪比区域，往往受到错误平层的影响，误码率随着信噪比的增加几乎不再下降。受困于此，在误码性能要求较高的系统中，例如磁存储和光纤通信中，ldpc码的性能仍不足以满足系统的需求。

为了解决上述问题，随后feldman等人首次将最大似然ml译码模型松弛为线性规划lp译码模型，成功应用于二元线性分组码的译码，奠定了lp译码强大的数学理论支撑。同时，feldman还证明了lp译码具有ml特性，码字独立特性和全零假设等良好特性；lp译码可通过伪码字图、最小分式距离等分析译码性能；对于校验矩阵中存在短环的ldpc码，lp译码算法可通过增加冗余校验节点，消除短环影响，提高译码性能。lp译码优点很多，但是译码复杂度很高，求解困难，严重阻碍了它在实际场景中的应用。

为了解决这一难题，barman等人将传统的线性规划译码模型与交替方向乘子算法admm相结合，提出了一种迭代式投影译码算法，但是该算法的最大缺点是与传统的置信传播bp译码相比，在低信噪比时误码率性能较差。针对这个问题，liu等人提出了基于交替方向乘子算法admm的罚函数译码。此类译码是在barman提出模型的目标函数上增加不同的惩罚项来加大对非整数码字的惩罚，从而提高了误码性能，但其不足是由于每次迭代过程都需要奇偶多面体的投影导致译码速度较慢。针对此问题，xiaojiezhang和guoqiangzhang等人各自提出了优化算法去降低奇偶多面体的投影。除此之外，xiaopengjiao等人提出了一种减少投影的方法来降低交替方向乘子算法admm中投影的次数，但是该方法需要进行大量的奇偶多面体投影操作，消耗很大的译码时间。

上述所有方法中由于每次进行交替方向乘子算法admm迭代时仍然需要进行投影操作，延长了译码时间，而且这种操作在实际应用中也不易实现。虽然传统的bp译码相较于上述采用交替方向乘子算法admm的方法译码时间要短，但是在高信噪比时，bp译码比上述采用交替方向乘子算法admm的方法误码性能差，会出现错误平台。

技术实现要素：

本发明的目的在于针对上述现有技术的不足，提出一种基于度分解和交替乘子算法的ldpc码的二次规划译码方法，以在保证低误帧率的前提下，进一步提高ldpc码译码效率，减少译码时间，满足现代无线通信系统的需求。

本发明的技术方案是：将ldpc码的最大似然译码通过校验方程分解的方法等价转换为线性整数规划问题，再通过在目标函数上加惩罚项将等价的线性整数规划问题松弛为二次规划模型；利用此二次规划模型具有的稀疏性和正交性特性，采用基于分布式并行快速方法的交替方向乘子算法admm对基于度分解的二次规划模型进行求解。其实现步骤包括如下：

1、基于度分解和交替乘子法的ldpc码二次规划译码方法，包括：

(1)根据线性整数规划的定义，利用对数似然将最大似然ml译码模型转化为如下线性整数规划译码模型：

其中，γ表示对数似然比向量，xi表示发送的第i个码元，i＝1,2,...,n，n表示码元的总个数，j＝1,2,...,m表示第j个校验方程，m表示校验方程的总个数，x表示译码的码字，hji表示校验矩阵h＝(hji)m×n第j行第i列的数，表示校验方程，表示模2运算；

(2)分解校验方程:

当校验方程的元素个数即度等于3时，将校验方程变换为等价校验方程的不等式形式；

当校验方程的元素个数即度大于3时，先将校验方程分解成多个元素个数为3的子校验方程，再将这些子校验方程变换为等价校验方程的不等式形式；

(3)建立基于度分解的二次规划译码：

(3a)根据步骤(2)得到的等价校验方程的不等式形式，将<1>中的线性整数规划译码模型等价转化为如下基于度分解的线性整数规划译码模型：

其中，q表示扩展后对数似然比向量，t表示矩阵的转置，d表示扩展后的码字，a表示系数矩阵，b表示系数向量，n表示码字d的长度；

(3b)再将模型<2>中的目标函数上增加二次惩罚项，并将模型<2>中约束条件的0和1约束松弛到0到1区间，得到以下基于度分解的二次规划译码模型：

其中，表示增加的二次惩罚项，用于加大对分数解的惩罚，以更有利于获得整数解；α表示二次惩罚项系数；

(3c)对基于度分解的二次规划译码模型进行变形，即对模型<3>中的不等式约束条件增加辅助变量z，得到如下等式约束的二次规划译码模型：

(4)根据模型<4>求解期望的码字：

(4a)将模型<4>中的约束条件添加到目标函数中，得到如下建立的增广拉格朗日函数：

其中，lμ(d,z,λ)表示拉格朗日函数，λ表示拉格朗日对偶变量，μ表示惩罚参数，表示ad+z-b的2-范数平方；

(4b)定义第k+1次迭代后原始残差r^k+1＝ad^k+1+z^k+1-b，对偶残差s^k+1＝z^k+1-z^k，利用交替方向乘子算法对式<5>中进行循环迭代求解，当原始残差2-范数的平方和对偶残差2-范数的平方同时达到设定的阈值或者最大迭代次数时停止迭代，得到最优的扩展码字d^*，并从d^*中提取出译码的码字x^*。

本发明与现有技术相比具有以下优点：

1.降低了误码率。

常用的bp译码由于受到码字结构图中伪码字、陷阱集合等有害特性的影响，在高信噪比下会出现误码平台，误码率不再降低；本发明与bp译码相比，利用数学理论，使码字具有全零特性，减少了码字固有特性的影响，在低信噪比的情况下可以获得与bp相近的误码性能，在高信噪比下不会出现误码平台，保持着较好的瀑区性能。

2.提高了译码速度。

传统的基于交替方向乘子算法admm的惩罚译码，比常用的bp译码在高信噪比误码性能好，但是译码过程中需要调用投影算法，大大降低了译码速度；本发明相比调用投影算法的交替方向乘子算法admm译码，避免了复杂的投影操作，且充分利用了矩阵的稀疏性和正交性，因而在不降低误码性能的前提下，极大地提高了译码速度，并且易于实现和并行操作。

附图说明

图1是本发明的实现流程图；

图2是本发明中进行校验方程度分解原理图；

图3是用本发明与现有译码方法对常用[2640,1320]“margulis”ldpc码译码的误比特率比较图；

图4是用本发明与现有译码方法对常用[2640,1320]“margulis”ldpc码译码的平均译码时间比较图。

具体实施方式

以下结合附图对本发明的实施例及效果作进一步详细描述。

本实施实例是对规则的ldpc码进行信道译码。

参照图1，本实例的实现步骤如下：

步骤1：根据线性规划的标准型，利用对数似然比向量将最大似然ml译码模型转化为线性整数规划译码。

(1a)假设发送端发送的二进制低密度奇偶校验码ldpc码字为x＝{x1,…,xi,…,xn}，该码字对应的校验矩阵为h＝(hji)m×n，经过噪声为n＝{n1,…,ni,…,nn}的加性高斯白噪声awgn信道后，接收的码字为r＝{r1,…,ri,…,rn}，其中xi表示发送的第i个码元，ri表示接收的第i个码元，i＝1,2,...,n，n表示码元的总个数，j＝1,2,...,m表示第j个校验方程，m表示校验方程的总个数，(hji)m×n表示m×n的校验矩阵第j行第i列的数；

(1b)计算对数似然比向量：γ＝[γ1,...,γi,...,γn]^t：

其中第i个对数似然γi为：

由于在加性高斯白噪声awgn信道中，噪声是均值为0，方差为的高斯随机变量，服从正态分布，所以有：

得到第i个对数似然γi为：

最终得到对数似然比向量如下：

其中，t表示矩阵的转置，e表示指数，n0表示高斯白噪声功率；

(1c)将二进制低密度奇偶校验码ldpc的码字集合x∈{0,1}ⁿ代入到校验方程中，得到满足每一行校验方程的约束条件：

j＝1,2,...,m

其中，x∈{0,1}ⁿ表示n维向量x中的元素等于0或1，表示校验方程，表示模2运算；

(1d)根据(1b)中得到的γ和(1c)中得到的约束条件构建线性整数规划译码模型：

步骤2：分解校验方程。

对不同度的校验方程，采用不同的分解方法，由于本实例采用的分解方法是针对于校验方程度大于等于3的情况，因此将分解步骤分成度等于3和度大于3两种具体情况进行分解。

(2a)对于度等于3的校验方程，将校验方程变换为等价校验方程的不等式形式：令校验方程中的校验矩阵h第j行非零元素对应的索引为j1，j2和j3，且hj1＝hj2＝hj3＝1，并令与hj1，hj2和hj3相对应的参与校验方程的变量分别为xj1，xj2和xj3；将校验方程进行模2运算得到等价校验方程该等价校验方程的不等式形式如下：

(2b)对于度大于3的校验方程，先进行分解使每个子校验方程的度为3，再按(2a)将其变换为等价校验方程的不等式形式：

(2b1)假设第j个校验方程的度为λj＞3，则第j个约束条件为：

根据hji的取值为0或1，将第j个约束条件化简为：

其中，表示参与第j个校验方程的第λp个码元，λp＝1,2,...,λj，表示序列号；

(2b2)对于参与第j个校验方程的λj个码元，增加辅助变量以分解校验方程：

令分解之前码元的个数l⁽⁰⁾＝λj，对于第v次分解，令增加的辅助变量个数根据辅助变量l^(v)取值奇偶数的不同，为校验方程设置不同的满足条件：

当l^(v)为奇数时，校验方程满足如下条件：

当l^(v)为偶数时，校验方程满足如下条件：

根据校验方程满足的条件，得到每次分解后的子校验方程的度都为3，参与子校验方程的码元集合为其中，表示整数值上界，log2表示以2为底的对数，g表示序号；

根据因子图节点和边的关系，对度为λj≥3的校验方程，增加λj-3个辅助变量，通过分解得到λj-2个度为3的子校验方程。

例如，对于一个如图2所示的度为6的校验方程分解原理图，参与校验方程的原始码元分别为对校验方程的这6个原始码元进行分解时需要增加3个辅助码元，分别为

将原始码元和辅助码元按以下步骤进行结合：

将相邻的两个原始码元和与增加的辅助码元相结合，构成一个码元集合使子校验方程的度为3；

将相邻的两个原始码元和与增加的辅助码元相结合，构成一个码元集合使子校验方程的度为3；

将相邻的两个原始码元和与增加的辅助码元相结合，构成一个码元集合使子校验方程的度为3；

将增加的3个辅助码元相结合，构成一个码元集合使子校验方程的度为3；

通过以上分解，可将度为6的校验方程分解为4个度为3的子校验方程。

步骤3：建立基于度分解的二次规划译码模型并建立增广拉格朗日函数。

(3a)根据步骤2构造的等价的校验方程的不等式形式，建立基于度分解的线性整数规划译码模型：

(3a1)定义为辅助变量的总个数，为分解出的度为3的校验方程的总个数，将原始变量x和辅助变量合并扩展为将对数似然比向量扩展为则模型<1>中的目标函数转化为minq^td，其中，表示1行γa列的向量，表示1行γc列的向量全为0；

(3a2)对于第γc个子校验方程，假设与扩展后的码字d相连的码元变量为则定义其相对应的矩阵为根据步骤(2a)中等价校验方程的不等式，利用线性方程组的矩阵形式，令不等式右侧的值用向量t表示,即t＝[2,0,0,0]^t,不等式左侧的系数用矩阵f表示，即得到第γc个度为3的校验方程的矩阵形式为：

对于γc个度为3的校验方程，令系数矩阵系数向量则模型<1>中的约束条件转化为ad≤b，0≤d≤1，其中，表示第γc个子校验方程中与扩展后的码字d相连的第一个码元，表示第γc个子校验方程中与扩展后的码字d相连的第二个码元，表示第γc个子校验方程中与扩展后的码字d相连的第三个码元，矩阵中只有对应位置的元素为1，其他位置元素为零，表示笛卡尔积，表示长度为γc的元素全为1，表示第γc个度为3的校验方程的系数矩阵；

由于系数矩阵f中有12个非零元素且任意两列是相互正交的，且系数矩阵与系数矩阵f相比，只是把系数矩阵f的每一列分别对应的搬移到了第列的位置，不影响非零元素的个数，所以中只有12个非零元素且任意两列也是相互正交的；由于系数矩阵a是由γc个度为3的校验方程直接级联得到的，无需改变中任意两列的正交关系，所以系数矩阵a具有正交性；由于系数矩阵a有4γc×(n+γa)个元素，其中只有12γc个非零元素，所以系数矩阵a具有稀疏性；

(3a3)根据(3a1)得到的目标函数和(3a2)得到的约束条件，构建得到与ml译码等价的基于度分解的线性整数规划译码模型：

其中，q表示扩展后对数似然比向量，t表示矩阵的转置，d表示扩展后的码字，a表示系数矩阵，b表示系数向量，n表示码字d的长度；

(3b)在模型<2>的目标函数上增加二次惩罚项，并且将约束条件的0和1约束松弛为0到1区间，此时模型<2>可以松弛为以下基于度分解的二次规划译码模型：

其中，表示增加的二次惩罚项，其作用是为了加大对分数解的惩罚从而更有利于获得整数解，α表示二次惩罚项系数；

(3c)对基于度分解的二次规划译码模型进行变形，即对模型<3>的不等式约束条件增加辅助变量z将其转化成如下等式约束模型：

(3d)将模型<4>中的约束条件添加到目标函数中，得到如下建立的增广拉格朗日函数：

其中，lμ(d,z,λ)表示拉格朗日函数，λ表示拉格朗日对偶变量，μ表示惩罚参数，表示ad+z-b的2-范数平方。

步骤4：对扩展后的码字d，辅助变量z，拉格朗日对偶变量λ进行循环迭代求解。

循环迭代求解的方法包括增广拉格朗日乘子法和交替方向乘子算法admm，本发明采用交替方向乘子算法admm进行循环迭代求解，该admm算法是将对偶上升法的可分解性和增广拉格朗日乘子法的上界收敛性融合在一起的算法，其求解步骤如下：

(4a)利用交替方向乘子算法admm进行迭代更新，求解第k+1次迭代后的扩展后的码字d^k+1，辅助变量z^k+1，拉格朗日对偶变量λ^k+1：

(4a1)固定辅助变量z^k和拉格朗日对偶变量λ^k，对扩展后的码字d进行更新。因为矩阵a具有列正交性，所以a^ta是一个对角阵，对角元素为矩阵a每一列非零元素的个数。因而式<5>关于变量d是可拆分的。令矩阵a中第i列向量为ai和第i列的非零元素个数为ei，那么式<5>可以拆分成以下n+γa个子问题：

目标函数：

约束条件：0≤di≤1

其中，i∈{1,...,n+γa}。求解以上每个子问题的步骤为：对上式目标函数关于di进行求导，并令导数等于零，然后把求解的结果投影到[0,1]区间，得到更新后的第i码字

其中，π[0,1]表示在超立方体[0,1]上的投影操作，t表示矩阵的转置；

(4a2)对辅助变量z进行更新，即固定更新后的码字d^k+1和拉格朗日对偶变量λ^k，将使得式<5>为

再对该式lμ(d^k+1,z,λ^k)中的辅助变量z求导，并令导数等于零，得到更新后的辅助变量z^k+1：

其中，表示在区间上的投影操作；为了加快计算z^k+1，以上公式可以写成4γc个公式进行并行计算，即：

其中，j∈{1,...,4γc}，和表示矩阵a的第j行向量。

(4a3)将(4a1)更新后的码字d^k+1和(4a2)更新后的辅助变量z^k+1代入到admm算法中拉格朗日对偶变量的更新公式中，得到更新后的拉格朗日对偶变量λ^k+1：

(4b)定义第k+1次迭代后原始残差r^k+1＝ad^k+1+z^k+1-b，对偶残差s^k+1＝z^k+1-z^k，在迭代求解过程中，当原始残差2-范数的平方和对偶残差2-范数的平方同时小于设定的阈值或者达到最大迭代次数时停止迭代，得到最优的扩展码字d^*，并从中提取出译码的码字x^*。

本发明的效果通过以下仿真结果进一步说明：

仿真方法：本发明、现有基于交替方向乘子算法admm的lp译码方法、现有的基于高复杂度奇偶多面体交替方向乘子算法admm译码的对比方法1、现有基于有效奇偶多面体投影交替方向乘子算法admm译码的对比方法2和现有置信传播bp译码方法。

仿真1：用本发明、现有基于交替方向乘子算法admm的lp译码方法和现有置信传播bp译码方法分别对常用的[2640,1320]“margulis”规则ldpc码进行译码，比较其误帧率fer，其结果如图3所示；

由图3可见，将本发明和现有基于交替方向乘子算法admm的lp译码方法对[2640,1320]规则ldpc码的译码结果进行对比，本发明在低信噪比的误码性能比lp译码方法好很多，并且在高信噪比没有出现错误平层，这说明本发明在误码性能方面克服了lp译码在低信噪比性能差的缺点，并且保持了lp译码低错误平层的优点；

同样，将本发明和现有置信传播bp译码方法对[2640,1320]规则ldpc码的译码结果进行对比，在低信噪比区域，本发明达到了和bp译码相近的误码性能；在高信噪比区域，bp译码出现了错误平层，而本发明仍保持较好的瀑区性能，没有出现错误平层；说明本发明在误码性能上和传统的bp译码相比是有优越性的；

仿真2：用本发明、现有的基于高复杂度奇偶多面体交替方向乘子算法admm译码对比方法1、现有基于有效奇偶多面体投影交替方向乘子算法admm译码对比方法2和现有置信传播bp译码方法分别对常用的[2640,1320]“margulis”规则ldpc码进行译码，比较其平均译码时间，其结果如图4所示；

由图4可见，将本发明、现有的基于高复杂度奇偶多面体admm算法译码对比方法1、现有基于有效奇偶多面体投影admm算法译码对比方法2对[2640,1320]规则ldpc码的译码结果进行对比，得到本发明的平均译码时间比现有基于admm对比方法1和2所用的时间短，这说明本发明的平均译码速度比现有的基于admm对比方法1和2平均译码速度都快；

同样，将本发明和现有置信传播bp译码方法对[2640,1320]规则ldpc码的译码结果进行对比，得到本发明的平均译码时间比bp译码方法的平均译码时间短，说明本发明是快速的有效译码方法。