using PyPlot
Homework #6 Programming Assignment Solution¶
Problem #1¶
In this homework, we examine the properties of binary convolutional codes and implement Viterbi decoders. Note: in the following, we always assume the input is a single sequence of bits.
Basic functions¶
We will represent each state by string. For example, all zero state with memory length 2 corresponds to "00".
The following basic operations will be useful.
str2int(x) = (x=="0"?0:1)
array2str(A) = string(A...)
str2array(parsetype, s) = [parse(parsetype, x) for x in split(s,"")]'
In Viterbi algorithm, when we find the best path, we need to track back the possible previous states.
prev(s, bit) = s[2:end]*bit
hammdist(s1,s2) = sum((str2array(Int, s1)+str2array(Int, s2))%2)
Functions related to channels:
bern(p) = rand(1)[1]<p ? 1 : 0
genRandBitVec(k) = [rand([0,1]) for x in 1:k]
function BSC(c, p)
# input: code vector c
# output: noisy vector y
y = str2array(Int, c)
for i in 1:length(c)
y[i] += bern(p)
end
return array2str(y % 2)
end
function AWGN(c, SNRdB)
# input: code vector c
# output: noisy vector y
y = 2*str2array(Float64, c)-1 # maps (0,1) to (-1,1)
N = 10^(-SNRdB/10)
y += sqrt(N)*randn(length(y))'
return y
end
(a)¶
Write a program for a function $\texttt{convEncode}$($\texttt{connections},$ $\texttt{inputbits}$) that takes an array of octal tap connections $\texttt{connections}$ and a sequence of bits $\texttt{inputbits}$ as inputs and outputs the encoded bit sequence. For example, for $\texttt{connections} = [\begin{matrix}4 & 5\end{matrix}]$ and $\texttt{inputbits} = 1001010001110,$ your function should output the bit sequence [ 11000111001000010011111001. ] Use this function to encode the bit sequence $111111111111111111111111111$ using
- the convolutional code with generator matrix $G(D) = [\begin{matrix}1+D^2 & 1+D+D^2 \end{matrix}],$ i.e., the code with octal tap connections $[\begin{matrix}5 & 7 \end{matrix}].$
- the convolutional code with generator matrix $G(D) = [\begin{matrix}1 & 1+D \end{matrix}],$ i.e., the code with octal tap connections $[\begin{matrix}2 & 3 \end{matrix}].$
First, given a tap-connections configuration, we first translate it into a matrix form.
octal2bin(oct) = *([bin(parse(Int, string(oct)[ii]), 3) for ii in 1:length(string(oct))]...)
function connmatrix(connections)
nummemory = Int(ceil(log2(maximum(connections))))-1
connections_matrix = zeros(Int64, length(connections), nummemory+1)
for ii in 1:length(connections)
connections_matrix[ii, :] = str2array(Int, octal2bin(connections[ii]))[end-nummemory:end]
end
return connections_matrix
end
For example,
connmatrix([12 13])
Using this function, we can easily implement the encoder.
function convEncode(connections, inputbits)
nummemory = Int(ceil(log2(maximum(connections))))-1
connections_matrix = connmatrix(connections)
outputbits = []
memory = zeros(Int64, nummemory+1)'
for j = 1:length(inputbits)
memory = [inputbits[j] memory[1:end-1]...]
output = (memory*connections_matrix')%2
append!(outputbits, output)
end
return array2str(outputbits)
end
connections = [4 5]
inputbits = str2array(Int, "1001010001110")
codebits = convEncode(connections, inputbits)
connections = [5 7]
inputbits = str2array(Int, "111111111111111111111111111")
codebits = convEncode(connections, inputbits)
connections = [2 3]
inputbits = str2array(Int, "111111111111111111111111111")
codebits = convEncode(connections, inputbits)
(c)¶
Write a program for a function $\texttt{ViterbiBSC}$($\texttt{connections},$ $\texttt{receivedsequence}$) that takes an array of octal tap connections $\texttt{connections}$ of a convolutional code, as well as the received sequence when the code is used over a BSC with $p<1/2,$ and outputs the most likely input sequence using Viterbi decoding. Verify that this function works as expected for the code with connections $[\begin{matrix}5 & 7 \end{matrix}].$ You should initialize your decoder at the all-zero state.
(e)¶
Write a program for a function $\texttt{ViterbiAWGN}$($\texttt{connections},$ $\texttt{receivedsequence}$) that takes an array of octal tap connections $\texttt{connections}$ of a convolutional code, as well as the received sequence when the code is used over an AWGN, and outputs the most likely input sequence using Viterbi decoding. Here, we assume that BPSK with unit power is used for transmission, i.e., the symbol $0$ is mapped to $-1$ and $1$ is mapped to $1$ while transmitting over the channel. Verify that this function works as expected for the code with connections $[\begin{matrix}5 & 7 \end{matrix}].$ You should initialize your decoder at the all-zero state.
First we define one-step codebook for a given encoder to ease the decoding process.
function convEncodeOneStep(connections, state, input)
s = [parse(Int, ss) for ss in split(state,"")]'
connections_matrix = connmatrix(connections)
memory = [input, s...]'
output = (memory*connections_matrix')%2
return array2str(output)
end
For example, given a tap-connection [5 7], if the state was "00" and the input bit was 1, then the output of the encoder is
convEncodeOneStep([5 7], "00", 1)
It would be helpful if we have a lookup table of the encoder.
function onestep(connections)
nummemory = Int(ceil(log2(maximum(connections))))-1
numoutput = length(connections)
state2idx = Dict(bin(i-1, nummemory) => i for i in 1:2^nummemory)
states = [bin(i-1,nummemory) for i in 1:2^nummemory]
onestep_codebook_index = Dict()
i = 0
onestep_codebook = zeros(Int64, 2^(nummemory+1), numoutput)
onestep_codebook = ["" for i in 1:2^(nummemory+1)]
for state in states
for input in [0 1]
i += 1
onestep_codebook_index[(state, input)] = i
onestep_codebook[i] = convEncodeOneStep(connections, state, input)
end
end
return (onestep_codebook, onestep_codebook_index)
end
For example, given a tap-connection [5 7],
onestep_codebook, onestep_codebook_index = onestep([5 7]);
onestep_codebook_index
onestep_codebook
You can read off the information as follows: for example, if the state was "01" and the input bit was 0, then the output of the encoder is stored as the third element of onestep_codebook.
Now, we implement the Viterbi algorithm. The only difference between BSC and AWGN is the metric we use.
function Viterbi(channel, connections, receivedsequence)
# Choose a metric.
if channel == "BSC"
dist = hammdist
elseif channel == "AWGN"
sqdist(code, received) = sum(((2*str2array(Float64, code)-1)'-received).^2)
dist = sqdist
else
print("ERROR: Channel\n")
return false
end
# Set parameters
nummemory = Int(ceil(log2(maximum(connections))))-1
numoutput = length(connections)
numinputbits = Int(length(receivedsequence)/numoutput)
state2idx = Dict(bin(i-1, nummemory) => i for i in 1:2^nummemory)
states = [bin(i-1,nummemory) for i in 1:2^nummemory]
# Initialize survivor metrics and paths
survivormetrics = zeros(2^(nummemory))
survivormetrics[2:end] = Inf
survivorpaths = ["" for i in 1:2^nummemory]
survivorinputpaths = ["" for i in 1:2^nummemory]
# One step codebook
onestep_codebook, onestep_codebook_index = onestep(connections)
# For each step, we do the following:
for ii in 1:numinputbits
received = receivedsequence[numoutput*(ii-1)+1:numoutput*ii]
distances = [dist(code, received) for code in onestep_codebook]
next_survivorpaths = copy(survivorpaths)
next_survivorinputpaths = copy(survivorinputpaths)
next_survivormetrics = copy(survivormetrics)
# For each state, we find the optimal path:
for state in states
state_idx = state2idx[state]
# Add
new_survmetric = [survivormetrics[state2idx[prev(state, s)]] + distances[onestep_codebook_index[(prev(state, s), parse(Int, state[1]))]] for s in ["0" "1"]]
# Compare & Select
optidx = new_survmetric[1]<new_survmetric[2]?1:2
opts = optidx==1?"0":"1"
optprev = prev(state, opts)
# Update
next_survivormetrics[state_idx] = new_survmetric[optidx]
opt_path = onestep_codebook[onestep_codebook_index[(optprev, parse(Int, state[1]))]]
next_survivorpaths[state_idx] = *(survivorpaths[state2idx[optprev]], opt_path)
next_survivorinputpaths[state_idx] = *(survivorinputpaths[state2idx[optprev]], string(state[1]))
end
# Update the survivor paths and metrics
survivorpaths = copy(next_survivorpaths)
survivorinputpaths = copy(next_survivorinputpaths)
survivormetrics = copy(next_survivormetrics)
# print(survivormetrics)
end
# print(survivormetrics, "\n")
# Pick the first minimum survivor metric
return survivorpaths[indmin(survivormetrics)], survivorinputpaths[indmin(survivormetrics)]
end
(d)¶
Use the functions in parts (a) and (c) to plot the bit error rate (BER) performance of the codes $[\begin{matrix}2 & 3\end{matrix}]$ and $[\begin{matrix}5 & 7\end{matrix}]$ over a BSC($p$) with $p\in [0\colon0.05\colon0.5]$. For each value of $p$, you should encode and decode a block of at least 2,000 symbols and calculate the average bit error rate. Please make sure that the encoder always start from the all-zero state. (No termination is required.)
iterN = 10
leninputbits = 2000
plist = 0:0.05:0.5
connectionslist = [[2 3], [5 7], [155, 117]]
BERlist = zeros(size(connectionslist,1),length(plist))
for ii in 1:size(connectionslist, 1)
connections = connectionslist[ii]
onestep_codebook, onestep_codebook_index = onestep(connections)
for jj in 1:length(plist)
print(jj, " ")
p = plist[jj]
BER = 0
for kk = 1:iterN
inputbits = genRandBitVec(leninputbits)
codebits = convEncode(connections, inputbits)
receivedbits = BSC(codebits, p)
decodedbits, decodedinputbits = Viterbi("BSC", connections, receivedbits)
BER += (hammdist(string(inputbits...), decodedinputbits))/length(inputbits)
end
BERlist[ii, jj] = BER/iterN
end
end
for ii = 1:length(connectionslist)
plot(plist, BERlist[ii,:], marker = "*", label = "Taps "*string(connectionslist[ii]))
end
plot(plist, plist)
grid("on")
title("BSC(p)")
xlabel("Bit flip probability p")
ylabel("Average bit error rate")
legend(loc = "upper left", fancybox = "true")
(f)¶
Use the functions in parts (a) and (e) to plot the bit error rate (BER) performance of the codes $[\begin{matrix}2 & 3\end{matrix}]$ and $[\begin{matrix}5 & 7\end{matrix}]$ over an AWGN with SNRdB ranging from $-1$ dB to $10$ dB. For each value of the SNRdB, you should encode and decode a block of at least 2,000 symbols and calculate the average bit error rate. Please make sure that the encoder always start from the all-zero state. (No termination is required.)
# Experiemnts: AWGN
iterN = 15
leninputbits = 2000
SNRdBlist = -1:1:10
connectionslist = [[2 3], [5 7], [155 117]]
BERlist = zeros(size(connectionslist,1),length(SNRdBlist))
for ii in 1:size(connectionslist, 1)
connections = connectionslist[ii]
for jj in 1:length(SNRdBlist)
print(jj, " ")
SNRdB = SNRdBlist[jj]
BER = 0
for kk = 1:iterN
inputbits = genRandBitVec(leninputbits)
codebits = convEncode(connections, inputbits)
receivedbits = AWGN(codebits, SNRdB)
decodedbits, decodedinputbits = Viterbi("AWGN", connections, receivedbits)
BER += (hammdist(string(inputbits...), decodedinputbits))/length(inputbits)
end
BERlist[ii, jj] = BER/iterN
end
end
for ii = 1:length(connectionslist)
plot(SNRdBlist, log10(BERlist[ii,:]), marker = "*", label = "Taps "*string(connectionslist[ii]))
end
grid("on")
title("AWGN(SNRdB)")
xlabel("SNR (dB)")
ylabel("log10(Average bit error rate)")
legend(loc = "upper right", fancybox = "true")
(b)¶
Write a program for a function $\texttt{freeDistance}$($\texttt{connections}$) that takes an array of octal tap connections $\texttt{connections}$ as input and outputs the free distance $d_{\textrm{free}}$ of the code. For example, for $\texttt{connections} = [\begin{matrix}4 & 5\end{matrix}],$ your function should output $$ d_{\textrm{free}} = 3. $$ Use this function to find $d_{\mathrm{free}}$ for the two codes in part (a), as well as for the code with octal tap connections $[\begin{matrix}155 & 117 \end{matrix}].$
We can easily modify the Viterbi algorithm.
function freeDistance(connections)
# Set parameters
nummemory = Int(ceil(log2(maximum(connections))))-1
numoutput = length(connections)
state2idx = Dict(bin(i-1, nummemory) => i for i in 1:2^nummemory)
states = [bin(i-1,nummemory) for i in 1:2^nummemory]
# Initialize survivor metrics and paths
survivormetrics = zeros(2^(nummemory))
survivormetrics[2:end] = Inf
survivorpaths = ["" for i in 1:2^nummemory]
# One step codebook
onestep_codebook, onestep_codebook_index = onestep(connections)
# Initialize the minimum distance and path.
mindist = Inf
minpath = ""
while true
received = ^("0", numoutput)
distances = [hammdist(code, received) for code in onestep_codebook]
next_survivorpaths = copy(survivorpaths)
next_survivormetrics = copy(survivormetrics)
for state in states
state_idx = state2idx[state]
if state == ^("0", nummemory) # For all-zero state, the optimal path is always from "00...01".
# Add
new_survmetric = [survivormetrics[state2idx[prev(state, "1")]] + distances[onestep_codebook_index[(prev(state, "1"), parse(Int, state[1]))]] for s in ["0" "1"]]
# Update
optidx = new_survmetric[1]<new_survmetric[2]?1:2
optprev = prev(state, "1")
else
new_survmetric = [survivormetrics[state2idx[prev(state, s)]] + distances[onestep_codebook_index[(prev(state, s), parse(Int, state[1]))]] for s in ["0" "1"]]
# Compare & Select
optidx = new_survmetric[1]<new_survmetric[2]?1:2
opts = optidx==1?"0":"1"
optprev = prev(state, opts)
end
# Update
next_survivormetrics[state_idx] = new_survmetric[optidx]
optpath = onestep_codebook[onestep_codebook_index[(optprev, parse(Int, state[1]))]]
next_survivorpaths[state_idx] = *(survivorpaths[state2idx[optprev]], optpath)
end
survivorpaths = copy(next_survivorpaths)
survivormetrics = copy(next_survivormetrics)
# Update the minimum distance
if mindist > survivormetrics[state2idx[^("0", nummemory)]]
mindist = survivormetrics[state2idx[^("0", nummemory)]]
minpath = survivorpaths[state2idx[^("0", nummemory)]]
end
if mindist < minimum(survivormetrics)
break
end
end
return Int(mindist), minpath
end
freeDistance([2 3])
freeDistance([5 7])
freeDistance([4 5])
freeDistance([155 117])