Homework #5 Programming Assignment Solution¶
In [1]:
using PyPlot
Basic functions¶
In [2]:
function bern(p)
x = rand(1)[1]
return x<p ? 1 : 0
end
function BSC(c, p)
# input: code vector c
# output: noisy vector y
y = copy(c)
for i in 1:length(c)
y[i] += bern(p)
end
return y % 2
end
Out[2]:
Generate random input bits $u\in \{0,1\}^4$ and take $c=uG$.
In [3]:
function genRandBitVec(k)
# return a uniform binary (row) vector from [0:2^{k-1}]
u = [parse(Int, s) for s in bin(rand(0:2^(k-1)), k)]'
end
Out[3]:
Generate a codebook with a given generator matrix:
In [4]:
function genCodebook(G)
k, n = size(G)
codebook = Array{Int}(2^k, n)
for i in 0:2^k-1
u = [parse(Int, s) for s in bin(i, k)]' # input bit; row vector
codebook[i+1, :] = u*G % 2
end
return codebook
end
Out[4]:
Problem 1¶
In this homework, we examine the performance of (7,4) Hamming code compared to other (7,4) codes over BSC$(p)$, the binary symmetric channel with bit error probability $p.$ Throughout this problem, we consider the systematic generator matrix $G$ of the following form: $G = [I_k| A ],$ where $I_k$ is the $k\times k$ identity matrix and $A$ is a $k\times(n − k)$ matrix. For example, for the (7,4) Hamming code, we consider¶
$$ \begin{align*} G & = \begin{bmatrix} 1 & 0 & 0 & 0 & 1 & 1 & 0\\ 0 & 1 & 0 & 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 & 1 & 1 & 1\end{bmatrix}.\end{align*}$$
In this case, recall that the parity-check matrix $H$ is obtained easily by¶
$$\begin{align*} H & = [A^t | I_{n-k}].\end{align*}$$
Consider the systematic generator matrix $G = [I_{k}, A]$:
In [5]:
G = [1 0 0 0 1 1 0;
0 1 0 0 1 0 1;
0 0 1 0 0 1 1;
0 0 0 1 1 1 1]
Out[5]:
In [6]:
codebook = genCodebook(G)
Out[6]:
In [7]:
dist = sum(codebook, 2)'
Out[7]:
Then the parity check matrix can be easily obtained by $H=[A^T,I_{n-k}]$:
In [8]:
k, n = size(G)
H = Array{Int}(n-k, n)
H[:, 1:k] = G[:, k+1:end]'
H[:, k+1:end] = eye(n-k)
H
Out[8]:
Check orthogonality: $HG^T=0$
In [9]:
H*G' % 2
Out[9]:
a) Recall that the minimum distance decoding is maximum likelihood decoding for BSC$(p)$, where $p < 1/2.$ Write a minimum distance decoder function $\texttt{mindistdec(y,C)}.$ This function takes received bits $\texttt{y}$ and a codebook $\texttt{C}$ as inputs.
In [10]:
function mindistdec(y, codebook)
added = (codebook + repmat(y, size(codebook,1), 1)) % 2
dist = sum(added, 2)
mindist = minimum(dist)
minidx_list = find(x->x==mindist, dist)
return codebook[rand(minidx_list), :]' # if there is tie, pick one at random
end
Out[10]:
b) Write a syndrome decoder function $\texttt{syndromedec(y,H)},$ which takes received bits $\texttt{y}$ and a parity check matrix $\texttt{H}$ as inputs.¶
In [11]:
function syndromedec(y)
# input: received bits y
# output: corrected codeword chat
chat = copy(y)
r = H*y' % 2
erridx = 0
for j = 1:7
if norm(H[:, j] - r) < 1
erridx = j
break
end
end
if erridx != 0
chat[erridx] = (chat[erridx]+1) % 2
end
return chat
end
Out[11]:
c) For the $(7,4)$ Hamming code, verify that any one bit error in code bits can be corrected by both $\texttt{mindistdec}$ and $\texttt{syndromedec}$. Hence, for simplicity, we now use $\texttt{mindistdec}$ as our decoder for the rest of this problem.¶
For a sanity check of our decoders:
In [12]:
for i in 1:7
z = zeros(Int64, 7)'
z[i] = 1
corruptedcodebook = (codebook + repmat(z, size(codebook,1), 1)) %2
for j in 1:size(codebook, 1)
c = corruptedcodebook[j, :]' # row vector
chat1 = syndromedec(c)
chat2 = mindistdec(c, codebook)
if codebook[j,:]' != chat1
print("Hamming decoder is wrong!\n")
elseif codebook[j,:]' != chat2
print("Minimum distance decoder is wrong!\n")
end
# print(codebook[j,:]'==chat1, " ", codebook[j,:]'==chat2, "\n")
# print(codebook[j,:]', " ", chat1, " ", chat2, "\n")
end
end
d) For the (7, 4) Hamming code, find the formula of the probability of decoding to a wrong codeword when the code is used for BSC$(p)$.¶
The Hamming code can correct 1 error bit. Hence, decoding succeeds when if and only if less than 1 bit is flipped. Analytically therefore we have the probability of failure: $1-(1-p)^7-7(1-p)^6p=1-(1-p)^7(1+6p)$.
e) We want to run a simulation that verifies the result found in part (d). For $p = 0.05 : 0.05 : 0.5,$ do the following steps $1,000$ times to get a good estimate of the decoded word error rate (WER).¶
i. Generate 4 random binary digits.¶
ii. Encode them by calculating the 3 check digits, or using generator matrix G.¶
iii. Put the codeword through BSC(p).¶
iv. Decode the received 7-bit vector to a 4-bit word.¶
v. Compare the decoded word with the transmitted word and count word errors to find word error rate (WER). Also count the errors in the information positions of the code words to obtain an estimate of the decoded bit error probability (BER).¶
Plot your simulation results $p$ vs. BER$(p)$ and $p$ vs. WER$(p),$ along with the analytical expression found in part (d). Is the bit error probability larger or smaller than the word error probability?¶
Find error rate of a given code.
In [13]:
# Find BER and WER over BSC(p)
function findER(G, plist, nblock=10^4)
k, n = size(G)
codebook = genCodebook(G)
ber_list = zeros(Float64, 1:length(plist))
wer_list = zeros(Float64, 1:length(plist))
for jj in 1:length(plist)
p = plist[jj]
for ii = 1:nblock
u = genRandBitVec(k) # random input bits
c = u*G % 2 # code bits
y = BSC(c, p) # received bits
chat = mindistdec(y, codebook) # corrected bits
dist = sum((c[1:k]+chat[1:k])%2)
ber_list[jj] += dist/k
wer_list[jj] += dist>0?1:0
end
end
ber_list /= nblock
wer_list /= nblock
return ber_list, wer_list
end
Out[13]:
Now with (7,4) Hamming code:
In [25]:
plist = 0.05:0.05:0.5
ber_list_ham, wer_list_ham = findER(G, plist, 10^4)
plot(plist, ber_list_ham, label="BER", marker="o")
plot(plist, wer_list_ham, label="WER", marker="*")
# plot(plist, 4*ber_list_ham, label="4*BER", marker="*")
plot(plist, 1-(1-plist).^6.*(1+6*plist), label="WER: Theory")
title("(7,4) Hamming code")
xlabel("Error probability p")
ylabel("Error rate")
legend(loc = "upper left", fancybox = "true")
Out[25]:
f) Now we consider another systematic generator matrix¶
$$ \begin{align*} G & = \begin{bmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 1 & 1 & 1 & 1\end{bmatrix}.\end{align*} $$
Repeat the same experiment as in (e), and compare the results with BER and WER of the (7,4) Hamming code.¶
In [15]:
G2 = [1 0 0 0 1 0 0;
0 1 0 0 0 1 0;
0 0 1 0 0 0 1;
0 0 0 1 1 1 1]
Out[15]:
In [16]:
codebook2 = genCodebook(G2)
Out[16]:
The minimum distance of this code is 2.
In [17]:
weight = sum(codebook2, 2)'
Out[17]:
In [18]:
ber_list, wer_list = findER(G2, plist)
plot(plist, ber_list, label="BER (f)", marker="o")
plot(plist, wer_list, label="WER (f)", marker="o")
plot(plist, ber_list_ham, label="BER: (7,4) Hamming code")
plot(plist, wer_list_ham, label="WER: (7,4) Hamming code")
# plot(plist, 4*ber_list, label="4*BER", marker="*")
xlabel("Error probability p")
ylabel("Error rate")
legend(loc = "upper left", fancybox = "true")
Out[18]:
We see that Hamming code uniformly outperforms the second code. This is expected, since $d_{\text{min}}$ for the second code is only 2 (verify this!), whereas it is 3 for Hamming.
Assuming that if there are multiple codewords closest to the received vector, then the minimum distance decoder chooses one of the closest codewords uniformly at random, can you show that the word error probability for the second code over a BSC$(p)$ is given by $$ \begin{align*} P_e & = 3p(1-p)^6+18p^2(1-p)^5+35p^3(1-p)^3+21p^5(1-p)^2+7p^6(1-p)+p^7 ? \end{align*} $$
(Hint: you might want to use the following modified minimum distance decoder.)
In [19]:
function mindistdec_debug(y, codebook)
added = (codebook + repmat(y, size(codebook,1), 1)) % 2
dist = sum(added, 2)
mindist = minimum(dist)
minidx_list = find(x->x==mindist, dist)
return minidx_list # if there is tie, pick one at random
end
Out[19]:
Let us verify that theory matches experiments for the second code.
In [20]:
coeff = [0.0 3.0 18.0 35.0 35.0 21.0 7.0 1.0]
Out[20]:
In [21]:
function thwer(p, coeff)
s = zeros(length(p))
for w in 0:length(coeff)-1
s += coeff[w+1]*p.^w.*(1-p).^(7-w)
end
return s
end
Out[21]:
In [22]:
plot(plist, thwer(plist, coeff), label="WER: Theoretical")
plot(plist, wer_list, label="WER (f)", marker="o")
xlabel("Error probability p")
ylabel("Error rate")
legend(loc = "upper left", fancybox = "true")
Out[22]:
g) Finally we examine the performance of a randomly generated (7, 4) systematic generator matrix. Generate a random 4×3 binary matrix $A$ and use $G = [I_4 | A]$ as our generator matrix. For each $G,$ repeat the same experiment in (e), and take averages of BERs and WERs over the random matrices. Compare the results with BER and WER of the (7,4) Hamming code. Is there any $p$ where a randomly generated $G$ performs better than the Hamming (7,4) code on average?¶
Generate a random systematic generator matrix $G=[I_k, A]$. We decide $A$ as a $k\times(n-k)$ random matrix where each entry is drawn from Bernoulli(1/2) independently.
In [23]:
niter = 10^3
avgber_list = zeros(Float64, 1:length(plist))
avgwer_list = zeros(Float64, 1:length(plist))
n, k = 7, 4
tic();
for ll = 1:niter
A = reshape([bern(1/2) for i in 1:k*(n-k)], k, n-k)
Grand = [eye(Int64, k) A]
ber_list_new, wer_list_new = findER(Grand, plist, 10^2)
avgber_list += ber_list_new
avgwer_list += wer_list_new
end
toc();
avgber_list /= niter;
avgwer_list /= niter;
In [24]:
plot(plist, avgber_list, label="BER: Random code", marker="*")
plot(plist, avgwer_list, label="WER: Random code", marker="*")
plot(plist, ber_list, label="BER (f)", marker="o")
plot(plist, wer_list, label="WER (f)", marker="o")
plot(plist, ber_list_ham, label="BER: (7,4) Hamming code")
plot(plist, wer_list_ham, label="WER: (7,4) Hamming code")
title("Random (7,4) code")
xlabel("Error probability p")
ylabel("Error rate")
legend(loc = "upper left", fancybox = "true")
Out[24]:
We see that Hamming code uniformly outperforms both the random code and the second code, but the random codes perform marginally better than the second code.