LDPC 7 - 解码简单例子

前言：

       前面我们知道NR 5G 里面BG1 46*68, BG2 42*52, 然后通过扩展因子生成一个特别大的奇偶校验矩阵H.这边用以（4,7）一个奇偶校验矩阵为例，阐述解码的过程。

解码过程很简单，SISO minsum原理后面会在补充章节单独讲。

主要参考NPTEL,台大的一些资料

目录：

     1： 总体流程

     2： 初始化

     3：行计算（SISO-minsum）

     4:  列计算（SISO-repetion）

     5：Decision

     6:  python code 

---

一 总体流程

      

收到数据

经过 Minsum 解码器，

得到码字

minsum 解码器主要工作：

         前面有讲过 SISO 的两个算法

         行运算： SISO minsum

        列运算：  SISO repetion

        然后迭代

         正常迭代8-10, 因为之后迭代error-rates就降低非常慢了。

---

二 初始化

     这里举个简单例子

    2.1  Storage Matrix L

         跟奇偶校验矩阵维度一致，通过这个矩阵更新每次计算后的结果，x是变量。

        非零值跟奇偶校验矩阵非零值的位置相对应，是一个in-place 操作

   

   

 2.2 初始化操作

     对照前面的奇偶校验矩阵，根据收到的信息r,  按列 做in-place  操作

      

      

例子

---

 三 行计算（SISO-minsum）

    3.1： min sum 操作

            min1 =  本行绝对值最小值

            min2 =  本行绝对值次小值

             min1 = 本行非零元素里面的绝对值最小值
             min2 = 本行非零元素里面的绝对值次小值
             设置所有值的大小= Min1
             设置最小值的大小 = min2

    2： sign 操作

             sign(parity) = 计算当前行元素的sign乘积
             sign( 新元素)= sign(原来元素)* sign(parity)

             如第一行，第一列元素：0.2 

             sign(parity)=sign(0.2)*sign(-0.3)*sign(1.2)*sign(0.8)=-1

             则经过operation 后

             sign(0.3)= sign(0.2)*sign(parity)=-1

             相当于做bit flip,奇偶校验失败就做Bit flip

    

   

行计算后结果：

    

---

四 列计算（SISO-repetion）

      sum of all entries in column j

      new entries = (old entry)

      本质上就是当前列某个码字更新为其它码字的和

---

五  Decision

    因为QPSK 是 0->1,1->-1

    所以最后要根据sumj的结果还原出原来的码字

    还原出来，后面可以根据 去检查一下decoder的结果

---

五 python code 

# -*- coding: utf-8 -*-
"""
Created on Mon May 23 16:14:00 2022

@author: chengxf2
LDPC 解码
"""
import numpy as np
from enum import Enum



class IndicsType (Enum):
      Row = 0 # 
      Column = 1

    


class LDPCDecoder():
    
    '''
    第一步初始化：
           生成L 矩阵存储Tanner图中的操作结果
    args
       rList: beliefs from channel
    return
       L: 存储操作结果
    '''
    def  init_para(self,rList):
        
        #奇偶校验矩阵
        H = [[1,1,1,0,1,0,0],
            [0,1,1,1,0,1,0],
            [1,1,0,1,0,0,1],
            [1,0,1,0,1,1,1]]
        
        self.H = np.array(H)
        m,n = np.shape(self.H)
        
        
        indics =np.where(self.H!=0)
        num = len(indics[0])
        
        self.rowIndics = {} #行
        self.colIndics = {} #列
        
        
        for i  in range(num):
            row = indics[0][i]
            col = indics[1][i] 
            if row not in self.rowIndics.keys():
                
                self.rowIndics[row]=[col]
            else:
                curItem = self.rowIndics[row]
                curItem.append(col)
                self.rowIndics[row]=curItem
                
        
        i = 0
        
        for i  in range(num):
            row = indics[0][i]
            col = indics[1][i]
            if col not in self.colIndics.keys():
                 
                self.colIndics[col]=[row]
            else:
                curItem = self.colIndics[col]
                curItem.append(row)
                self.colIndics[col]=curItem
    
                
        
        i = 0
        
        L = np.zeros((m,n))
        print(np.shape(L))
        
        
        for i  in range(m): #行
            for j in range(n): #列
                r=rList[j] #列
                h = self.H[i,j] #奇偶校验位
                belief = r*h
                #print("\n  i: %d j: %d"%(i,j),"\t item ",belief)
                L[i,j]=belief
        #print("\n第一步 根据channnel r 初始化 \n ",L)
        return L
    
    
    '''
    repetion code SISO decoding 算法
    args:
        r: data from channel,beliefs
        L: beliefs
    return
        L
    '''
    def column_operation(self, L,r):
        
        m,n = np.shape(L)
        
        sumjList =[]
        
        for col in range(n): #按列
        
            indics = self.get_indics(col, IndicsType.Column) #此列中非零元素的索引
            s =0.0
            for row in indics:
                s+=L[row,col]
            
            rj = r[col]
            sumj = np.round(s+rj,2) #此列和为更新后的beliefs
            sumjList.append(sumj)
            
            #print("\n column: %d:  sum:  %4.2f"%(col, sumj))
            for row in indics:
                old_entry =L[row,col]
                new_entry = sumj-old_entry
                L[row,col] = new_entry
        #print("\n column_operation \n",L)
        return L,sumjList
                
                
    '''
    根据收到的beliefs做出判断
    QPSK 是0 调制到幅度为1的信号
         1调制到幅度为-1的信息
    args
        rList: 收到的beliefs
    return
         hat_c: 解码得到码字
    '''       
    def decsion(self, rList):
        hat_c =[]
        for r in rList:
            if r>0:
                c =0
            else:
                c = 1
            hat_c.append(c)
        #print("\n rlist ",rList,"\n hat_c",hat_c)
        return np.array(hat_c).T
        
    '''
     找到此行或者此列非零元素的位置
     '''
    def get_indics(self,i,tp:IndicsType):

          if  tp is  IndicsType.Row: #找到此行中非零元素的下表
              indics = self.rowIndics[i]
          else:
              indics = self.colIndics[i]
          return indics
              
         
    
    '''
    SISO minsum 算法实现
    args
       L: tanner 图计算，bit node 把信息传递给check node, check node 执行
       minsum 计算，更新L 
    return
      L
    '''
    def row_operation(self, L):
        
        m,n = np.shape(L)
        
        for i in range(m): #执行奇偶校验,可以in parallely
            row = L[i,:]
            
            indics = self.get_indics(i, IndicsType.Row)
            
            min1,min2, min1_idx,min2_idx,sign = self.minsum(row,indics)
            
            
            #做in-place 操作#
            for index in indics:
                
                r = L[i,index]
                if index == min1_idx:
                    L[i,index]= min2*sign*self.sign(r)
                else:
                    L[i,index]=min1*sign*self.sign(r)
            

        #print("\n row_operation \n",L)
        return L
            
            
 
        
    '''
    获得最小值和次小值的位置
    args
       row: 行
       indics: 非零元素index的索引
    return
       min1: 最小值
       min2: 次小值
       min1_idx: 最小值索引
       min2_idx: 次小值索引
       
    '''
    def getMin(self,row ,indics):
        min1 = -1
        min2 = -1
    
        item = []
 
        for index  in indics:

            belief = row[index]
            item.append(abs(belief))
            
        item.sort()
        min1 = item[0]
        min2 = item[1]
        min1_idx =  np.where(np.abs(row) == min1)[0][0]
        min2_idx =  np.where(np.abs(row) == min2)[0][0]
        
        #print(min1,min2, min1_idx, min2_idx)
        
        return min1,min2, min1_idx,min2_idx
    
    '''
    符号函数
    args
       a
    return
      sign
    '''
    def sign(self, a):
        
        if a>0:
            return 1
        else:
            return -1
            
    '''
    sign: 符号函数
    '''
    def signProduct(self,lis,indics):
        
        s = 1
     
        
        for index in  indics:

            a = lis[index]
            sgn = self.sign(a)
           
            s =s*sgn
        #print("\n sign ",s)
        return s
    
    '''
    执行minsum 运算
    args
        找到此行非零元素，找出其中的最小值，次小值
        执行sign 连乘操作
    return
        查找结果  
    '''
    def minsum(self, row,indics):
        

        #print("\n row ",row, "\t indics ",indics)
        min1,min2, min1_idx,min2_idx= self.getMin(row, indics)
        sign= self.signProduct(row, indics)
        
        return min1,min2, min1_idx,min2_idx,sign
    


    '''
    CRC 校验结果
    args
       hat_c: 预测值
    '''
    def check_crc(self,hat_c):
        failnum = 0
        for row in self.H:
            a = np.multiply(row, hat_c)

            num = np.count_nonzero(a) #非零个数
            
            if num %2 == 1:
                failnum+=1
            
        
        
        #print("\n crc : ",failnum)
        return failnum
        
        
        
    
    def __init__(self,r):
        self.rList = r
        self.epoch = 5 #最大迭代次数
        self.rowIndics = None #非零元素行索引
        self.colIndics = None #非零元素列索引
        
        
    '''
    LDPC解码
    '''    
    def decoder(self):
        
        L= self.init_para(self.rList)
        m,n = np.shape(L)
        
       
        hat_c = self.decsion(self.rList)
        crc_fail = self.check_crc(hat_c)
        print("\n 开始CRC 校验失败的数目:  %d"%crc_fail)
        
        for step  in range(self.epoch):
        
           
            L= self.row_operation(L)
            L,rList= self.column_operation(L,self.rList)
            print("\n rList ",rList)
            
            #相当于repetion 过程
         
            hat_c = self.decsion(rList)
            crc_fail=self.check_crc(hat_c)
            print("\n step :",step,"\t 校验失败的数目: %d"%crc_fail)
      
#beliefs from channel        
r=[0.2,-0.3,1.2,-0.5,0.8,0.6,-1.1]
ldpc = LDPCDecoder(r)
ldpc.decoder()

    

开始CRC 校验失败的数目:  3

 step : 0      crc_fail number: 0

 step : 1      crc_fail number: 0

 step : 2      crc_fail number: 0

 step : 3      crc_fail number: 0

 step : 4      crc_fail number: 0