% Program to implement a 1/2-rate convolutional code operating over
% a binary symmetric channel (BSC) with bit error probability p 
% using hard decisions and Viterbi decoding at the output. The
% program runs the algorithm for T - 1 bits while keeping
% only WindowLen survivors. The program then calculates the number 
% of bits in error after the decoder and plots the resulting BER
% after channel decoding against the bit error probability in the
% channel and thus shows how much improvement in bit error rate
% the code provides.

% Code adopted from R. Kozick's rjkvit2.m. It employs functions
% initvit1.m, viterbi1.m, vitflush.m which are available from 
% the same course Web site as this code.

% The receiver trellis structure is defined by the structures
% trellis{}, branchweight{}, and sd. The structure trellis indicates 
% the successor nodes to every node in the trellis, for example 
% trellis{1} = [1 2]; means node 1's successor states are nodes 1 and
% 2. The structure branchweight indicates the labels on trellis 
% transitions for each "source" node. For example, branchweight{1,1}
% = [0 0] means the label on the first (upper) transition leaving node 
% 1 is 00, and branchweight{1,2} = [1 1] means the label on the second 
% (lower) transition leaving node 1 is 11. Note the presence of the 
% space in [0 0] and [1 1]. The structure sd describes which input 
% bits the transitions correspond to. Every row of the structure sd
% corresponds to one "source" node. In one row, each entry
% corresponds to a "destination" node. If there is no transition from
% the source node to that particular node, the entry is -1; if there 
% is a transition on input bit 0, the entry is 0; if there is a
% transition on input bit 1, the entry is 1. WindowLen is the duration
% the Viterbi algorithm keeps its search. It is generally accepted 
% that this number should be about 5 times the constraint length (or
% memory) of the code. The code below is for a convolutional code of 
% 2 and Windowlen is then 10. For convolutional codes with memory 3
% and 4, it should be increased to 15 and 20, respectively. 

% The structures above represent the trellis in the receiver. The 
% corresponding convolutional code is described by the variables x(2),
% x(1), and the description of the shift register contents wD2 = wD1, 
% wD1 = w. Here, wD2 is the rightmost content of the shift register,
% wD1 is the one left to the rightmost bit, and w represents the input 
% bit. Note wD1 and wD2 are initialized to 0 at the beginning of the 
% program. For more memory (larger constraint length), wD3 (and more 
% as needed) can be added. When this is done wD3 (and others) should
% also be initialized to zero at the beginning, and the shift register 
% shift wD3 = wD2 (and more, such as wD4 = wD3 if needed) should be 
% added.

% The variable T is the number of bits the algorithm runs over. For
% design, it can be kept at about T = 1e4 + 1. For reporting purposes,
% use a minimum of T = 1e5 + 1. A better number to use is T = 1e6 + 1,
% in which case the results for lower probabilities of error become
% meaningful but then running simulations for one convolutional code 
% can take of the order of an hour or more.

% You may wish to change the lcoation of the legend in the figure by
% an appropriate command. Upper left corner seems to be most suited.

% Convolutional codes are usually represented by the octal
% representations of their generator polynomials. This particular code
% corresponds to the octal polynomials 5,7. 

% rjkvit2.m -- Test Viterbi algorithm for decoding convolutional codes.
% Uses MATLAB functions initvit1.m, viterbi1.m, vitflush.m from the book
% T.K. Moon and W.C. Stirling, "Mathematical Methods and Algorithms for
% Signal Processing," Prentice-Hall, 2000.
%
% Rich Kozick, Nov. 2001

% ** This version stores only WindowLen past message bits

% Define the trellis for the convolutional code with rate 1/2 and
% g_1(D) = 1 + D + D^2 and g_2(D) = 1 + D^2.

trellis{1} = [1 2];
trellis{2} = [3 4];
trellis{3} = [1 2];
trellis{4} = [3 4];

branchweight{1,1} = [0 0];   branchweight{1,2} = [1 1];
branchweight{2,1} = [1 0];   branchweight{2,2} = [0 1];
branchweight{3,1} = [1 1];   branchweight{3,2} = [0 0];
branchweight{4,1} = [0 1];   branchweight{4,2} = [1 0];

WindowLen = 10;    % Window length for Viterbi algorithm
                   % (typically 5 times the memory in shift regs.)

% Initialize Viterbi algorithm
initvit1(trellis, branchweight, WindowLen, 'dH');
svec = [];

T = 1e5 + 1;   % Number of message bits to generate
wD1 = 0;  wD2 = 0;   % Initial state of shift register
xvec = [];  yvec = [];

% Probability of bit error for BSC

pvec = [.01, .05, .1];
Np = length(pvec);

% Matrix to decode states to message bits

sd = [ 0,  1, -1, -1; 
      -1, -1,  0,  1; 
       0,  1, -1, -1; 
      -1, -1,  0,  1];

for pind = 1:Np
p = pvec(pind);
errcnt = 0;
wvec = [];
% wrec = [];   whrec = [];

for t=1:T

% Generate message bit at time t
z = rand;
if (z <= 0.5)
  w = 0;
else
  w = 1;
end

if (t <= WindowLen+1)
  wvec = [wvec, w];
else
  wvec = [wvec(2:WindowLen+1), w];
end
% wrec(t) = w;

% Generate 2 bits from convolutional coder

x(2) = xor(w, wD2);
x(1) = xor(wD1, x(2));
wD2 = wD1;  wD1 = w;    % Update shift register

% Pass coded bits x through BSC to generate received bits y

y = x;
z = rand(1,2);
zi = find(z <= p);
y(zi) = not(x(zi));

% Apply Viterbi algorithm to y

s = viterbi1(y);  

if (t == WindowLen)
  sD1 = s;
elseif (t > WindowLen)
  wh = sd(sD1, s);
% whrec(t-WindowLen) = wh;
  errcnt = errcnt + (wvec(1) ~= wh);
  sD1 = s;
end

end  % Loop over time t=1:T

% Flush remaining paths

svec = vitflush(0); 
for t=1:WindowLen-1
  wh = sd(sD1, svec(t));
  sD1 = svec(t);
% whrec(T-WindowLen+t) = wh;
  errcnt = errcnt + (wvec(1+t) ~= wh);
end

fprintf('p = %f: %d errors in %d message bits\n', p, errcnt, T-1);
BER(pind) = errcnt/(T-1);

end   % Loop over p values

% Compute upper bound on probability of bit error

pfine = logspace(-4, 0, 200);
zi = find(pfine <= 0.5);
pfine = pfine(zi);
PeU = 1 - (1-pfine).^7 - 7*(pfine).*(1-pfine).^6;  % Hamming(7,4)
PeL = PeU/4;

PeRep3 = pfine.^3 + 3*(pfine.^2).*(1-pfine);       % 3x repetition

figure(1)
loglog(pvec, BER, 'bo:', pfine, PeU, 'r-', ...
  pfine, PeL, 'g-', pfine, pfine, 'k-', pfine, PeRep3, 'c--');
legend('CONV. CODE SIMULATION', 'HAMMING(7,4) U.B.', ...
'HAMMING(7,4) L.B.', 'UNCODED', '3 \times REPETITION CODE')
xlabel('ERROR PROBABILITY FOR BSC (p)')
ylabel('BER')
title('CODES OVER BSC')
axis([1e-4 0.5 1e-6 1])
% print -depsc fig1_7.eps

