6.02 Lab #5: Convolutional Codes

Write a decoder for convolutional codes based on the Viterbi
algorithm. Measure the bit error rate for both hard-decision
branch metrics and soft-decision branch metrics.

Instructions

See Lab #1 for a longer narrative about lab mechanics. The 1-sentence version:

Complete the tasks below, submit your task files on-line before the
deadline, and complete your check-off interview within a week of the
file submission deadline.

As always, help is available in 32-083 (see the Lab Hours web page).

Introduction

A convolutional encoder is characterized by two parameters: a
constraint length k and a set of r generator functions
{G0, G1, ...}. The encoder works through the message
one bit at a time, generating a set of r parity bits {p0, p1, ...}
by applying the generator functions to the current message bit,
x[n], and k-1 of the previous message bits,
x[n-1], x[n-2], ..., x[n-k-1]. The r
parity bits are then transmitted and the encoder moves on to the next
message bit. Since r parity bits are transmitted for
each message bit, the code rate is 1/r.

The operation of the encoder is best described as a state machine.
The figure below is the state transition diagram for a rate 1/2 encoder
with k=3 using the following two generator functions:

A generator function can described compactly by simply listing
its k coefficients as a k-bit binary sequence, or even more
compactly if we construe the k-bit sequence as an integer, e.g.,
G0: 7 and G1: 6.

In this diagram the states are labeled with the two previous message
bits, x[n-1] and x[n-2], in left-to-right order. The
arcs -- representing transitions between states as the encoder
completes the processing of the current message bit -- are labeled with
x[n]/p0p1.

You can read the transition diagram as follows: "If the encoder
is currently in state 11 and if the current message bit is a 0, transmit
the parity bits 01 and move to state 01 before processing the next
message bit." And so on, for each combination of current states and
message bits. The encoder starts in state 00.

The stream of parity bits arrives at the receiver; some of the bits
may be received with errors because of noise introduced by the
transmitter, channel and receiver. Based on the information in the
(possibly corrupted) parity bits, the decoder deduces the
sequence of states visited by the encoder and hence recovers the
transmitted message. Because of errors, the receiver cannot know
exactly what that sequence of states was, but it can determine
the most-likely sequence using the Viterbi algorithm.

The Viterbi algorithm works by determining a path
metricPM[s,n] for each state s and sample
time n. Consider all possible encoder state sequences that
leave the encoder in state s at time n. The
most-likely state sequence is the one that transmitted the
parity bit sequence that most closely matches the received parity bits, where
the closeness of the match is determined by the Hamming distance
between the transmitted and received parity bits: the smaller the
Hamming distance, the closer the match. Each increment in Hamming
distance corresponds to a bit error. The sequence with the smallest
Hamming distance between the transmitted and received parity bits
involves the fewest number of errors and hence is the most-likely
(fewer errors being more probable than more errors). PM[s,n]
records this smallest Hamming distance for each state at the specified
time.

The Viterbi algorithm computes PM[…,n] incrementally. Initially

PM[s,0] = 0 if s == starting_state else ∞

The algorithm uses the first set of r parity bits to
compute PM[…,1] from PM[…,0].
Then the next set of r parity bits is used to compute
PM[…,2] from PM[…,1] and so on until all the
received parity bits have been consumed.

Here's the algorithm for computing PM[…,n]
from PM[…,n-1] using the next set of r parity
bits to be processed:

For each state s:

Looking at the encoder's state transition diagram,
determine the two predecessor states α and β which
have transition arcs that arrive at state s.

Using the rate 1/2, k=3 encoder above,
if s is 10 then, in no particular order, α = 00 and
β = 01.

For the state transition α→s determine
the r parity bits the encoder would have transmitted;
call this r-bit sequence p_α. Similarly,
for the state transition β→s determine
the r parity bits the encoder would have transmitted;
call this r-bit sequence p_β.

Continuing the example from step 1: p_α = 11 and
p_β = 01.

Call the next set of r received parity bits
p_received. Compute the Hamming distance between
p_α and p_received. In the terminology
of the Viterbi algorithm this is branch metric for the
state transition α→s so we'll label it
BM_α. Similarly, compute the Hamming distance between
p_β and p_received, call it BM_β.

Compute two trial path metrics that correspond to the two
possible paths leading to state s:

PM_α = PM[α,n-1] + BM_α
PM_β = PM[β,n-1] + BM_β

PM_α is the Hamming distance between the transmitted
and received parity bits so far assuming we arrive at state s
via state α. Similarly, PM_β is the Hamming distance
between the transmitted and received parity bits so far assuming we
arrive at state s via state β.

Continuing the example from step 3: assuming
PM[α,n-1]=5 and PM[β,n-1]=3, then
PM_α=7 and PM_β=4.

Now compute PM[s,n] by picking the smaller of the two
trial path metrics. Also record which state we chose to be the
most-likely predecessor state for s:

We can use the following recipe at the receiver to determine
the transmitted message:

Initialize PM[…,0] as described above.

Use the Viterbi algorithm to compute PM[…,n] from
PM[…,n-1] and the next r parity bits;
repeat until all received parity bits have been consumed. Call the
last time point N.

Identify the most-likely ending state of the encoder by finding
the state s which has the mimimum value of PM[s,N].
If there are several states with the same minimum value, the end of
the message has been corrupted by errors, so decrement N and repeat
the search. Keep doing this until a unique s is found.

Determine message bit that caused the transition to
state s. (You can do this simply from knowing s;
you don't need to know the predecessor
state). Set s=Predecessor[s,N], decrement N and repeat this
step as long as N > 0, accumulating the message bits in reverse order.

Task 1: Implementing a Viterbi decoder? (6 points)

In this task we'll write the code for a Python class ViterbiDecoder.
One can make an instance of the class, supplying k and the parity generator functions,
and then use the instance to decode messages transmitted by the matching encoder.

The decoder will operate on a sequence of received voltage samples; the choice
of which sample to digitize to determine the message bit has already been made, so
there's one voltage sample for each bit. The transmitter has sent 0V for a "0"
and 1V for a "1" but those nominal voltages have been corrupted by noise described
by a Gaussian PDF with zero mean. 0's and 1's appear with equal probability in the
original message, so the decoder should use a 0.5V threshold when digitizing the voltages
into bits.

# This is the template file for Lab #5, Task #1
import numpy
import lab5
class ViterbiDecoder:
# given the constraint length and a list of parity generator
# functions, do the initial set up for the decoder. The
# following useful instance variables are created:
# self.k
# self.nstates
# self.r
# self.predecessor_states
# self.expected_parity
def __init__(self,k,glist):
self.k = k # constraint length
self.nstates = 2**(k-1) # number of states in state machine
# number of parity bits transmitted for each message bit
self.r = len(glist)
# States are named using (k-1)-bit integers in the range 0 to
# nstates-1. The bit representation of the integer corresponds
# to state label in the transition diagram. So state 10 is
# named with the integer 2, state 00 is named with the
# integer 0.
# for each state s, figure out the two states in the diagram
# that have transitions ending at state s. Record these two
# states as a two-element tuple.
self.predecessor_states = \
[((2*s+0) % self.nstates,(2*s+1) % self.nstates)
for s in xrange(self.nstates)]
# this is a 2D table implemented as a list of lists.
# self.expected_parity[s1][s2] returns the r-bit sequence
# of parity bits the encoder transmitted when make the
# state transition from s1 to s2.
self.expected_parity = \
[[lab5.expected_parity(s1,s2,k,glist) \
if s1 in self.predecessor_states[s2] else None
for s2 in xrange(self.nstates)]
for s1 in xrange(self.nstates)]
# expected is an r-element list of the expected parity bits
# (or you can also think of them as voltages given how we send
# bits down the channel). received is an r-element list of
# actual sampled voltages for the incoming parity bits.
# This is a hard-decision branch metric, so, as described in
# lab write up, digitize the received voltages to get bits and
# then compute the Hamming distance between the expected sequence
# and the received sequences, return that as the branch metric.
# Consider using lab5.hamming(seq1,seq2) which computes the
# Hamming distance between two binary sequences.
def branch_metric(self,expected,received):
pass # your code here...
# compute self.PM[...,n] from the batch of r parity bits and
# the path metrics for self.PM[...,n-1] computed on the previous
# iteration. Follow the algorithm described in the lab
# write up. In addition to making an entry for self.PM[n,s] for
# each state s, keep track of the most-likely predecessor
# for each state in the self.Predecessor array. You'll probably
# find the following instance variables and methods useful:
# self.predecessor_states
# self.expected_parity
# self.branch_metric()
def viterbi_step(self,n,received_voltages):
pass # your code here...
# Identify the most-likely ending state of the encoder by
# finding the state s which has the mimimum value of PM[s,n]
# where n points to the last column of the trellis. If there
# are several states with the same minimum value, the end of
# the message has been corrupted by errors, so decrement n
# and repeat the search. Keep doing this until a unique s is
# found. Return the tuple (s,n).
def most_likely_state(self,n):
pass # your code here...
# starting at state s at time n, use the Predecessor
# array to find all the states on the most-likely
# path. Each state contributes a message bit...
def traceback(self,s,n):
message = []
while n > 0:
# message bit that caused transition to
# state s is also the high-order bit of
# the state name
message.append(s >> (self.k-2))
# back to the next earlier state along the path
s = self.Predecessor[s,n]
n -= 1
message.reverse()
return message
# figure out what the transmitter sent from info in the
# received voltages
def decode(self,received_voltages,debug=False):
# figure out how many columns they'll be in the trellis
nreceived = len(received_voltages)
max_n = (nreceived/2) + 1
# this is the path metric trellis itself, organized as a
# 2D array: rows are the states, columns are the time points.
# PM[s,n] is the metric for the most-likely path through the
# trellis arriving at state s at time n.
self.PM = numpy.zeros((self.nstates,max_n),dtype=numpy.float)
# at time 0, the starting state is the most likely, the other
# states are "infinitely" worse.
self.PM[1:self.nstates,0] = 1000000
# a 2D array: rows are the states, columns are the time
# points, contents indicate the predecessor state for each
# current state.
self.Predecessor = numpy.zeros((self.nstates,max_n),
dtype=numpy.int)
# use the Viterbi algorithm to compute PM
# incrementally from the received parity bits.
n = 0
for i in xrange(0,nreceived,self.r):
n += 1
# Fill in the next columns of PM, Predecessor based
# on info in the next r incoming parity bits
self.viterbi_step(n,received_voltages[i:i+self.r])
# print out what was just added to the trellis state
if debug:
print self.PM[:,n],self.Predecessor[:,n]
# find the most-likely ending state from the last row
# of the trellis
s,n = self.most_likely_state(n)
# reconstruct message by tracing the most likely path
# back through the matrix using self.Predecessor.
return self.traceback(s,n)
# print out final path metrics
def dump_state(self):
print self.PM[:,-1]
if __name__=='__main__':
d = ViterbiDecoder(3,(7,6))
received = numpy.array([1,1,1,0,1,1,0,0,0,1,1,0,0,0])
message = d.decode(received,debug=True)
print "decoded message =",message

The template contains a start at the implementation -- you'll get to finish the job!
Here's a description of the functions you need to write:

number = branch_metric(self,expected,received)

expected is an r-element list of the expected
parity bits (or you can also think of them as voltages given how we
send bits down the channel). received is
an r-element list of actual sampled voltages for the incoming
parity bits. This is a hard-decision branch metric, so, as described
above, digitize the received voltages to get bits and then compute the
Hamming distance between the expected sequence and the received
sequences, return that as the branch metric.

Consider using lab5.hamming(seq1,seq2) which computes the Hamming
distance between two binary sequences.

viterbi_step(self,n,received_voltages)

compute self.PM[…,n] from the batch of r parity bits and
the path metrics for self.PM[…,n-1] computed on the previous
iteration. Follow the algorithm described in the lab
write up. In addition to making an entry for self.PM[s,n] for
each state s, keep track of the most-likely predecessor
for each state in the self.Predecessor array. You'll probably
find the following instance variables and methods useful:

self.predecessor_states
self.expected_parity
self.branch_metric()

s,n = most_likely_state(self,n)

Identify the most-likely ending state of the encoder by finding the
state s which has the mimimum value of PM[s,n] where
n points to the last column of the trellis. If there are several
states with the same minimum value, the end of the message has been
corrupted by errors, so decrement n and repeat the
search. Keep doing this until a unique s is found. Return
the info as a tuple (s,n).

The testing code at the end of the template makes an instance
of your decoder for the rate 1/2, k=3 code discussed in lecture.
The received sequence is identical to the example used
in Lecture #10, so the debugging printout showing the successive
calculation of the path metrics and predecessor states should match
what's shown on the slides.

Task 2: Error correction and BER (1 point)

In this task we'll run a simple experiment to measure the
improvement in bit-error rate (BER) achieved by using a convolutional
code. There's no code to write, just some results to analyze.

# This is the template file for Lab #5, Task #2
import lab5
from lab5_1 import ViterbiDecoder
if __name__=='__main__':
# encode a 500,000 bit message using the specified
# convolutional encoder, transmit it through a noisy
# channel and then decode it with an instance of the
# provided decoder class. Compare the BER with simply
# transmitting the raw message across the channel.
lab5.test_ber(ViterbiDecoder,k=3,glist=(7,6),nbits=500000)
# try it again with a stronger convolutional code
lab5.test_ber(ViterbiDecoder,k=4,glist=(0xE,0xD),nbits=500000)

Run the code (be patient, it takes a while). If it's taking
too long, try reducing nbits to 100000. Questions about the results:

How many actual errors were detected when sending the
message without coding?

How many uncorrected errors were detected when sending the
message with coding? Note that the encoded message takes twice
as many bits to transmit as the original message (i.e., 2*nbits
bits were transmitted across the channel).

Relate your answer to question 2 to the printout of
the trellis state at the end of the decoding.

What was the (most likely) last bit of the message? You can tell from
the final trellis state.

Task 3: Soft-decision branch metrics (3 points)

Let's change the branch metric to use a "soft" decision based on the
actually received voltages instead of just digitizing them directly to
bits. As pointed out in lecture, there's a difference in likelihood
between a "1" arriving as .999V and one arriving as .501V.

The soft metric we'll use is the square of the Euclidean distance
between the received vector (of dimension r) of analog voltage samples
and the expected r-bit vector of parities. Just treat them like a
pair of r-dimensional vectors and use the usual formula for computing
the square of the distance between the two vectors.

# This is the template file for Lab #5, Task #3
import lab5
from lab5_1 import ViterbiDecoder
class SoftViterbiDecoder(ViterbiDecoder):
# override the default branch metric with one based
# the square of the Euclidian distance between the
# expected and received voltages.
def branch_metric(self,expected,received):
pass # your code here...
if __name__=='__main__':
# try both decoders on exactly the same noisy
# received voltages
lab5.test_two(ViterbiDecoder,SoftViterbiDecoder,
k=3,glist=(7,6),nbits=100000)

Complete the coding of the branch_metric method for
the SoftViterbiDecoder class. Note that other than changing
the branch metric calculation, the decoders are identical.

The test_two test function tests two decoders on exactly
the same sequence of noisy received voltages. Since soft-decision
Viterbi decoding is more effective than hard-decision decoding, we've
increased the noise on the channel, so the results aren't directly
comparable to those of Task #2.

Try running your code several times. What were the results of the
bake-off? Should soft-decision decoders be preferred by better coders
nationwide?