### 6.02 Lab #5: Convolutional Codes

Due date: Wednesday, 3/18, at 11:59p

Goal

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:

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:

```G0: 111, i.e., p0[n] = 1*x[n] + 1*x[n-1] + 1*x[n-2]
G1: 110, i.e., p1[n] = 1*x[n] + 1*x[n-1] + 0*x[n-2]```

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 metric PM[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:

1. 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.

2. 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.

3. 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_β.

Continuing the example from step 2: assuming the received parity bits p_received = 00, then BM_α = hamming(11,00) = 2 and BM_β = hamming(01,00) = 1.

4. 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.

5. 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:

```if PM_α <= PM_β:
PM[s,n] = PM_α
Predecessor[s,n] = α
else
PM[s,n] = PM_β
Predecessor[s,n] = β```

Completing the example from step 4: PM[s,n]=4 and Predecessor[s,n]=β.

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

1. Initialize PM[…,0] as described above.

2. 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.

3. 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.

4. 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.

lab5_1.py is the template file for this task.

```# 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.

# 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()

# 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):

# 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
# figure out how many columns they'll be in the trellis

# 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
n += 1

# Fill in the next columns of PM, Predecessor based
# on info in the next r incoming parity bits

# 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))
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:

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.

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.

lab5_2.py is the template file for this task:

```# 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:

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

2. 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).

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

4. 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.

lab5_3.py is the template file for this task:

```# 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