Digital Link B¶
from IPython.core.display import HTML
import numpy as np
import matplotlib
import scipy
from scipy.stats import norm
from scipy.stats import binom
import pandas as pd
params = {'figure.figsize':(12,6), # These are plot parameters
'xtick.labelsize': 16,
'ytick.labelsize':16,
'axes.titlesize':18,
'axes.labelsize':18,
'lines.markersize':4,
'legend.fontsize': 20}
matplotlib.rcParams.update(params)
from matplotlib import pyplot as plt
import random
from ipywidgets import *
import numpy.linalg
from IPython.display import display
from IPython.core.display import HTML
from notebook.nbextensions import enable_nbextension
%matplotlib inline
print('The libraries loaded successfully')
The libraries loaded successfully
In this chapter, we continue our analysis of digital links. The LDPC codes are widely used in applications. Moreover, they are a great illustration of Bayes’ Rule.
LDPC Code¶
The LDPC code works as follows. One transmits an \(n\)-bit string \(\mathbf{x}\) that we represent as a column vector. To protect the transmission against errors, one appends an \(m\)-bit string \(\mathbf{y}\) that one computes as \(\mathbf{y} = H\mathbf{x}\) where \(H\) is an \(m\)-by-\(n\) matrix whose entries are \(0\) and \(1\) and the operations are modulo \(2\). Typically, most entries of \(H\) are zeros, hence the name low density parity codes. The discussion below is similar to that in the book. We include it here to be self-contained.
To explain the algorithm, consider the example illustrated in the figure below.
In this example, \(y_1 = x_1 + x_3 + x_4\). Also, \(\mathbf{x} = 0100101\) and \(\mathbf{y} = 001\). The transmitter sends \([\mathbf{x} | \mathbf{y}] = 0100101001\). There is one transmission error and the receiver gets \([\mathbf{x}' | \mathbf{y}'] = 0100101000\). Since \(\mathbf{y}' \neq H \mathbf{x}'\), the receiver knows that there was at least one transmission error. How can it locate the errors? Assuming that errors are unlikely, the receiver should find the minimum number of bits to flip in \([\mathbf{x}' | \mathbf{y}']\) to satisfy the equations. By looking exhaustively, the receiver can find that the equations are satisfied by flipping only the last bit. This takes \(n + m\) checks. However, if more than one bit has been flipped, this exhaustive exploration is too time-consuming to be feasible when \(m\) and \(n\) are large.
Instead of an exhaustive exploration, the algorithm proceeds iteratively as follows.
Step \(1\): The \(x\)-nodes \(\{x_1, \ldots, x_n\}\) on the left in the graph of the figure compute their estimates \(p_j\) of \(P(x_j = 1)\). For instance, since \(x_1^\prime = 0\) and the probability that an error flips the bit is \(\epsilon\), node \(x_1\) uses Bayes’ Rule to calculate
Each node \(x_j\) sends its estimate \(p_j\) to the \(y\)-nodes to which it is connected.
Node \(y_1\) updates its estimate \(q_1(1)\) of \(P(x_1 = 1)\) as follows. (The notation is \(q_i(j)\) for the estimate that node \(y_i\) maintains of \(P(x_j=1)\).) Node \(y_1\) has received the estimates \(p_1, p_3, p_4\) from the \(x\)-nodes and it knows that \(y_1 = x_1 + x_3 + x_4\) and \(y_1^\prime = 1\). Thus, \(x_1 = 1\) if \(y_1 = 1 + x_3 + x_4\), i.e., if an odd number of random variables among \(\{y_1, x_3, x_4\}\) are equal to one. Let us call that event \(E\). You can verify that \(z = [1 - (1 - 2y_1)(1 - 2x_3)(1 - 2x_4)]/2\) is equal to one if \(E\) occurs and to zero otherwise. Hence, assuming independence,
Given that \(y_1^\prime = 0\), we find as before, using Bayes’ Rule, that \(P[y_1 = 1|y_1^\prime = 0] = \epsilon\). Also, \(P(x_j = 1)= p_j\). Hence, we get
The other \(y\)-nodes update their estimates in the same way.
Step \(2\): The \(y\)-nodes send their estimates to the \(x\)-nodes to which they are connected. Consider node \(x_1\). It gets the estimates \(q_1(1), q_2(1)\) from nodes \(y_1\) and \(y_2\). Also, it remembers that \(x_1^\prime = 0\). Node \(x_1\) uses these observations to update \(p_1\) as follows:
To derive this formula, we use Bayes’ Rule again. Let \(Y_1\) and \(Y_2\) be the observations of nodes \(y_1\) and \(y_2\). These nodes provided their estimates \(q_1(1) = P[x_1 = 1 | Y_1]\) and \(q_2(1) = P[x_1 = 1|Y_2]\). Also, we know that \(x_1^\prime = 0\). We want to combine these values to compute \(P[x_1 = 1|x_1^\prime = 0, Y_1, Y_2]\). Bayes’ Rule gives
Now,
Moreover,
\end{eqnarray*}
Similarly, \(P[Y_1|x_1 = 0] = 2 P(Y_1)(1 - q_1(1))\). Putting together these formulas gives (6.1).
Step \(2k + 1\) similar to step \(1\).
Step \(2k + 2\) is similar to step \(2\).
Termination: The algorithm can run for a pre-determined number of steps or one may have a stopping rule based on the change in the estimates \(p_j\). When the algorithm stops, one decides that \(x_j = 1\) if \(p_j > 0.5\) and \(x_j = 0\) otherwise.
For a tutorial and very nice applications to image and sound coding, see https://pypi.python.org/pypi/pyldpc/0.7.4 Below, we propose a crude implementation, to illustrate the decoding algorithm.
Not surprisingly, the performance of the LDPC code depends strongly on the choice of the matrix \(H\).
The code below generates a suitable parity matrix \(H\) using a method propoed by Robert Gallager. It is from Hicham Janati: https://github.com/janatiH/pyldpc/blob/master/pyLDPC-Presentation.ipynb This code generates a regular \(H\) matrix with parameters \((n, a, b)\), defined as a matrix where each \(y_j\) is connected to the same number \(a\) of \(x_i\) and each \(x_i\) is connected to the same number \(b\) of \(y_j\). [Note: \(a < b\) and \(b\) divides \(n\).]
Here is our simple implementation of the LDPC decoding algorithm. The code generates a random vector \(\mathbf{x}\), chosen uniformly at random in \(\{0,1\}^n\). It then transmits \([\mathbf{x}, \mathbf{y} = H \mathbf{x}]\). The link flips each of those \(n+m\) bits independently with probability \(\epsilon\). The LDPC decoding algorithm runs for \(100\) steps. The code counts and print the number of decoding errors.
**Notes: **
(1) The code is not optimized and takes a few seconds for \(n = 36\).
(2) Recall that if \(- log_{10}(BER) = 2,\) then \(BER = 10^{2}\). Here, the value of \(- log_{10}(BER)\) is selected by the widget.
def dummy(berd,nd):
global ber, n
ber, n = float(berd), int(nd)
berd = widgets.Dropdown(options=['0.2', '0.4', '0.6', '0.8','1','1.2','1.4','1.6','1.8','2','2.2','2.4','2.6','2.8','3'],value='2',description='$BER$',disabled=False)
nd = widgets.Dropdown(options=['6', '8', '10', '12','14','18','20','22','24','26','28','30','32','34','36'],value='12',description='n',disabled=False)
z = widgets.interactive(dummy, berd = berd, nd = nd)
display(z)
def RegularH(n,d_v,d_c):
""" ------------------------------------------------------------------------------
Builds a regular Parity-Check Matrix H (n,d_v = a,d_c = b) following Gallager's algorithm :
----------------------------------------------------------------------------------
Parameters:
n: Number of columns (Same as number of coding bits)
d_v: number of ones per column (number of parity-check equations including a certain variable)
d_c: number of ones per row (number of variables participating in a certain parity-check equation);
----------------------------------------------------------------------------------
Errors:
The number of ones in the matrix is the same no matter how we calculate it (rows or columns), therefore, if m is
the number of rows in the matrix:
m*d_c = n*d_v with m < n (because H is a decoding matrix) => Parameters must verify:
0 - all integer parameters
1 - d_v < d_v
2 - d_c divides n
---------------------------------------------------------------------------------------
Returns: 2D-array (shape = (m,n))
"""
if n%d_c:
raise ValueError('d_c must divide n. Help(RegularH) for more info.')
if d_c <= d_v:
raise ValueError('d_c must be greater than d_v. Help(RegularH) for more info.')
m = (n*d_v)// d_c
a=m//d_v
Set=np.zeros((a,n),dtype=int)
# Filling the first set with consecutive ones in each row of the set
for i in range(a):
for j in range(i*d_c,(i+1)*d_c):
Set[i,j]=1
#Create list of Sets and append the first reference set
Sets=[]
Sets.append(Set.tolist())
#Create remaining sets by permutations of the first set's columns:
i=1
for i in range(1,d_v):
newSet = np.transpose(np.random.permutation(np.transpose(Set))).tolist()
Sets.append(newSet)
#Returns concatenated list of sest:
H = np.concatenate(Sets)
return H
def LDPC(BER, n):
# Print the regular H matrix
RegularH_test = 1
if RegularH_test:
H = RegularH(n,int(n/4), int(n/3))
print("Regular parity matrix (n, a, b) = ", (n, int(n/4),int(n/3)) , '\n \n H: \n', H)
N = 20 # steps in algorithm
eps = 10**(-BER) # BER
(m,n) = H.shape
x = np.arange(0.0, n)
xp = np.arange(0.0, n)
y = np.arange(0.0, m)
yp = np.arange(0.0, m)
p = np.arange(0.0, n) # guesses of P(x_i = 1) by xi
q = np.zeros((m,n)) # q(j,i) is guess about P(x_i = 1) by yj
err = 0 # number of transmission errors
s = 0
M = 100
for t in range(0,50): # simulations
for i in range(0,n):
x[i] = np.random.binomial(1,0.5) # transmitted string
w = np.random.binomial(1,eps)
err = err + w
xp[i] = np.remainder(x[i] + w,2) # received bits
v = np.dot(H,x)
for j in range(0,m):
y[j] = np.remainder(v[j],2) # transmitted parity check bits
w = np.random.binomial(1,eps)
err = err + w
yp[j] = np.remainder(y[j] + w,2) # received parity check
for i in range(0,n):
p[i] = eps + (1 - 2*eps)*xp[i] # initial guess
for t in range(0,N): # iterations
for j in range(0,m):
for i in range(0,n):
if H[j,i] == 1:
A = 1 - 2*(eps + (1 - 2*eps)*yp[j])
for k in range(0,n):
if (k != i) and H[j,k]==1:
A = A*(1 - 2*p[k])
q[j,i] = (1 - A)/2
for i in range(0,n):
A = eps + (1 - 2*eps)*xp[i]
B = 1 - A
for j in range(m):
if H[j,i] == 1:
A = A*q[j,i]
B = B*(1 - q[j,i])
p[i] = A/(A + B)
for i in range(0,n):
s = s + (x[i] != (p[i] > 0.5))
print('\n'"Number of bit transmission errors after 100 transmissions of n bits:",err)
print("Number of bit decoding errors after 100 transmissions of n bits:",s)
LDPC(ber, n)
Regular parity matrix (n, a, b) = (12, 3, 4)
H:
[[1 1 1 1 0 0 0 0 0 0 0 0]
[0 0 0 0 1 1 1 1 0 0 0 0]
[0 0 0 0 0 0 0 0 1 1 1 1]
[1 0 0 0 0 1 1 0 0 0 0 1]
[0 0 0 0 1 0 0 1 1 1 0 0]
[0 1 1 1 0 0 0 0 0 0 1 0]
[1 1 1 0 0 0 0 0 0 1 0 0]
[0 0 0 1 1 1 1 0 0 0 0 0]
[0 0 0 0 0 0 0 1 1 0 1 1]]
Number of bit transmission errors after 100 transmissions of n bits: 14
Number of bit decoding errors after 100 transmissions of n bits: 7