规则LDPC和不规则LDPC译码算法MATLAB仿真

Posted 2023-03-04 fpga和matlab

up目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

        LDPC码是根据低密度稀疏校验矩阵H来构造的。假设需要发送一组信息Tt_1,t_2,⋯,t_n，在发送前先使用生成矩阵G做线性变换，得到发送码字S=GTT，而校验矩阵H与生成矩阵G满足的关系为HGT=0，可以看到发送的码字是一组线性校验方程的解。

LDPC码与分组码的最大区别：在于译码方式的不同。LDPC码采用迭代译码，校验矩阵H使用Tanner图表示。

LDPC码的核心：围绕校验矩阵Hm×n的特性进行设计。

LDPC码本身的结构引入了交织(interleaving)。

不同LDPC码的构造方法都是为了实现如下的目标：① 增大Girth值；② 减小编码复杂度。

        根据H构造不同，分成了LDPC规则码和LDPC不规则码。直观地说，稀疏矩阵每行每列1的个数固定的LDPC码称为LDPC规则码。稀疏矩阵每行每列1的个数或者行和列的重量不固定的LDPC码称为LDPC不规则码。根据LDPC的码元是否为二进制，分成了二元LDPC码和高阶域LDPC码。高阶域LDPC码定义在GF(q)上，其中q= 2m。
       LDPC译码采用BP译码算法，该方法通过校验节点向变量节点传递信息以及变量节点向校验节点传递信息的方式，并经过迭代，得到近似ML的译码结果。
         这里说说自己的理解，LDPC编码的时候发送的原始比特会存在于多个校验方程中，一个校验方程中也会包含多个比特，这样就可以互相校验。对于某个校验方程，其他比特为1的个数为偶数，为使得校验结果为0，则本比特位置为0，反之则为1。而我们需要计算的是其他比特为偶数个1和奇数个1的概率。对于某个比特，由于其存在于多个校验方程中，因此，其他校验方程是可以给出本比特取值的概率的。在迭代过程中就是反复利用以上关系更新各节点的概率，最后得到译码结果的。

       根据H[c w]T=0以及上述矩阵结构，可以将w写成[p1 p2]两部分组成。其中p1的长度为4，p2的长度为42。那么式子H[s w]T=0可改写成，此处为便于区分用s指代以上的c。

可求得p1和p2的结果如下：

那么问题的关键是如何计算B-1。协议中Graph=1的情况下，B矩阵的形式如下： 

 

以上1表示对单位阵循环右移，0表示单位阵，X表示全0阵，那么根据逆矩阵的定义即可求得B-1，如下式：

      求取上述矩阵中的x1~x16即可得到逆矩阵。由于计算是定义在GF(2)上的，因此两个相同的矩阵相加为0，这样就能简便地计算出结果。另外，需要注意的是，在译码的时候，需要人为地在译码结果前面加上2Zc长度的0来补齐长度。 

二、核心程序

function x_hat = func_Dec(y,sigma,H,max_iter)

f1    = 1./(1+exp(-2*y/(sigma^2)));         
f0    = 1-f1;

[m,n] = size(H); 
if m>n
   H     = H'; 
   [m,n] = size(H); 
end

%如果不是稀疏矩阵，将H变为稀疏矩阵
if ~issparse(H) 
   [ii,jj,sH] = find(H);
   H = sparse(ii,jj,sH,m,n);
end

%初始化
[ii,jj] = find(H);        
%将非零元素标上
indx    = sub2ind(size(H),ii,jj); 
q0      = H * spdiags(f0(:),0,n,n);
sq0     = full(q0(indx)); 
sff0    = sq0;

q1 = H * spdiags(f1(:),0,n,n); 
sq1 = full(q1(indx));
sff1 = sq1;

k=0;
success = 0;

while ((success == 0) & (k < max_iter)),
   k = k+1;
   %若先验概率一样的话，令sdq为一个很小很小的数
   sdq = sq0 - sq1; 
   sdq(find(sdq==0)) = 1e-20; 
   %用概率差sdq构成m*n的稀疏矩阵
   dq = sparse(ii,jj,sdq,m,n);
   %将dq在对数域求和，即将dq相乘，构成数组
   Pdq_v = full(real(exp(sum(spfun('log',dq),2)))); 
   %更新了rij                                                                                                 
   Pdq = spdiags(Pdq_v(:),0,m,m) * H;
   sPdq = full(Pdq(indx));
   
   sr0 = (1+sPdq./sdq)./2; 
   sr0(find(abs(sr0) < 1e-20)) = 1e-20;%将每个元素的rij都归一化，更新rij。概率译码的第一步完成
   sr1 = (1-sPdq./sdq)./2; 
   sr1(find(abs(sr1) < 1e-20)) = 1e-20;%若比值非常小，则设为1e-20
   %将更新完的rij存入对应的位置，方便后面要用
   r0 = sparse(ii,jj,sr0,m,n);
   r1 = sparse(ii,jj,sr1,m,n);
   
   Pr0_v = full(real(exp(sum(spfun('log',r0),1))));
   Pr0 = H * spdiags(Pr0_v(:),0,n,n);
   sPr0 = full(Pr0(indx));
   Q0 = full(sum(sparse(ii,jj,sPr0.*sff0,m,n),1))';
   sq0 = sPr0.*sff0./sr0;
   
   Pr1_v = full(real(exp(sum(spfun('log',r1),1))));
   Pr1 = H * spdiags(Pr1_v(:),0,n,n);
   sPr1 = full(Pr1(indx)); 
   Q1 = full(sum(sparse(ii,jj,sPr1.*sff1,m,n),1))';
   sq1 = sPr1.*sff1./sr1;
   
   sqq = sq0+sq1;
   sq0 = sq0./sqq;
   sq1 = sq1./sqq;
   
   %伪后验概率
   QQ = Q0+Q1;
   Q0 = Q0./QQ;
   Q1 = Q1./QQ;
   
   %判决输出
    x_hat = (sign(Q1-Q0)+1)/2;
    if rem(H*x_hat,2) == 0, success = 1; end
end
up102

三、测试结果

在matlab2021a中仿真得到如下的效果：