【Python】LDPC码模拟实验

低密度奇偶校验码(LDPC)是在含噪声的通信信道中使用的纠错码,可降低错误率。

本文主要参考Ta, Tuan. “A Tutorial on Low Density Parity-Check Codes.” (2009)整理,并用Python实现LDPC解码过程和模拟信道传输。

构造校验矩阵

$m$行$n$列的校验矩阵,设每行有$w_r$个1,每列有$w_c$个1,且满足:

  • $w_r\ll m, w_c \ll n$
  • $w_r=w_c \frac{n}{m}$ (regular)

本文均假设$n=2m$,则$w_r=2w_c$。通过每次右移两位的方式(类似循环码),可以构造出符合要求的校验矩阵$H$。代码如下所示(用$p$代替$w_c$,最小为2)

def genH(m,n,p):
    H = np.zeros((m, n))
    for i in range(m):
        j=2*i
        for j in range(2*i,2*p+2*i):
            H[i,j%n]=1
            j=j+1
    return H

模拟信道传输

整个传输过程需要经过 编码->调制->传输->解调->解码 这些步骤。

为了省去生成矩阵的计算,我们直接使用全0码字(codeword)作为待传输消息,长度为$n$。

调制和解调

采用BPSK算法:

  • 调制(modulate):对于待发送的二元信息,将0编码为1,1编码为-1进行传输。
  • 解调(demodulate):对于接收到的信息,将大于0的值解调为0,小于0的值解调为1。

代码如下:

# BPSK(c for codeword)
def modulate(c):
    for i in range(np.size(c)):
        if c[i] == 0:
            c[i]=1
        else:
            c[i]=-1
    return c

def demodulate(y):
    for i in range(np.size(y)):
        if y[i] > 0:
            y[i]=0
        else:
            y[i]=1
    return y

噪声模拟

传输的过程会引入噪声,这里采用高斯噪声模拟,即消息的每个位置(entry)会受到一个满足高斯分布噪声信号的扰动。基础噪声分布为$N(0,1)$,并给予不同的信噪比(SNR),定义如下:

为了方便,取信号强度$p=1$,则噪声强度$\sigma=10^{-\frac{SNR}{20}}$。

这部分可以直接用Numpy包中的函数实现,注意使噪声信号长度和码字长度相等,便于直接相加。

sigma = 10**(-snr/20)
noise=np.random.normal(0,sigma,size=n)

LDPC解码

LDPC解码是个略为复杂的过程,这里实现其中较为简单的一种:Hard-decision Decoder

首先理解一下校验矩阵$H$的含义,$H$的每一行对应了一个奇偶校验方程,而每一列则分别对应了码字的每个位。$H[i,j]=1$表示码字的第$j$位需要由第$i$个方程进行校验。由这个思路,可将码字的每一位看作一个消息结点(message node),每一个方程看作一个校验节点(check node),由它们之间的联系可生成一幅信息传播图(也叫做Tanner Graph)。一个满足LDPC条件的$H$和对应的Tanner Graph如下图所示。

为了演示解码算法如何执行,假设我们有Codeword $c=[1,0,0,1,0,1,0,1]^T$,经过信道传输后收到消息$y=[1,1,0,1,0,1,0,1]$(第二位出现错误),下面通过信息节点和校验节点的相互通信来实现纠错(只是类比为通信,实际上没有发生信息传输),主要分为以下几步:

  • 信息节点将$y$的每一位依次发给与它相连的校验节点;
  • 校验节点根据每个信息节点发来的消息,构造出Receive表;
  • 校验节点对所有收到的信息求和,之后分别剔除每个信息节点的消息,作奇偶校验,以确定给每个节点发回的消息,构造出Send表;(Receive表和Send表共同组成下图中的Table 3)

  • 信息节点收到校验节点发回的消息后,与原有消息放入同一行,并依据每行1或0的个数作出Decision(少数服从多数,如第二行两个0和一个1,则decision为0,实现了纠错),组成Table4;

Table4的最后一行中$f1\to1$ 应该改为$f1\to 0$,属于paper原作者的笔误。

  • 最后输出解码后的值$y’$,过程结束。

在用Python实现的过程中需要注意索引的准确表达,特别是处理发送和接收时,可能需要多次遍历校验矩阵。个人实现的版本如下(这里同样假设n=2m):

# LDPC decode(hard decision)
def decode(H,y,m,n,p):
    fr = np.zeros((m, 2 * p)) # check nodes received
    fs = np.zeros((m, 2 * p)) # check nodes send
    sum = np.zeros(m) # check nodes received sum(for parity check)
    # message nodes table
    # 前p列为校验节点发来的消息,第p+1列为原始消息,第p+2列作为游标
    c=np.zeros((n,p+2)) 
    y1=np.zeros(n)

    # Fill the check nodes received table
    for i in range(m):
        count=0
        for j in range(n):
            if H[i][j] == 1:
                fr[i,count]=y[j]
                sum[i]=sum[i]+y[j]
                count = count+1

    # Calculate the check nodes send table
    for i in range(m):
        for j in range(2*p):
            fs[i,j]=(sum[i]-fr[i,j])%2

    # Fill the message node table
    for i in range(m):
        count=0
        for j in range(n):
            if H[i][j]==1:
                index=int(c[j,p+1])
                c[j,index]=fs[i,count]
                count = count+1
                c[j,p+1]+=1

    # Fill the last column with y
    for i in range(n):
        c[i, p] = y[i]

    # Decision
    for i in range(n):
        count=0
        for j in range(p+1):
            if c[i,j] == 1:
                count+=1
        if count > (p+1)/2:
            y1[i]=1
    return y1

用上面的消息$y=[1,1,0,1,0,1,0,1]$和对应的$H$测试,可以得到$y$被纠错后的结果$y1$。

m=4 # Number of rows
n=8 # Number of columns
p=2 # Number of 1s in a colomn

H=np.zeros((4,8))
H=[[0,1,0,1,1,0,0,1],[1,1,1,0,0,1,0,0],[0,0,1,0,0,1,1,1],[1,0,0,1,1,0,1,0]]

y=np.zeros(n)
y=[1,1,0,1,0,1,0,1]

y1 = decode(H,y,m,n,p)
print(y1)

错误率

最后计算一下不同密度和不同信噪比(即手动修改$p$和$SNR$,反复运行)的错误率:

def errorRate(c,y1):
    err=0
    total=np.size(c)
    i=0
    for i in range(total):
        if c[i] != y1[i]:
            err+=1
    r = err/total
    return r

可以描绘SNR和p对错误率影响的曲线:

结论

  • 根据上述曲线,可以初步判断出错误率随信噪比的提升而下降,前期总体呈线性趋势,最终缓慢收敛到0;

  • p的取值对错误率的影响较大。
  • 由于本文只进行了较为浅层的实验,构造的LDPC校验矩阵比较随意,且采用了Hard Decision作为解码和纠错标准,因此不能很好地体现出p的取值在LDPC中的作用。相关问题有更为细致的方法和深入的文献。

后记:LDPC码早在上世纪60年代就由Robert G. Gallager在博士论文中提出,当时并未引起太大轰动,近些年来随着软硬件的发展而变得应用广泛,因此LDPC码也称为Gallager码

ps. 本文的真实性质是《信息论与编码》课程的最后一个Groupwork,在写正式Report前给自己整理个思路。

晚饭加个🍗?