% 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 3];
trellis{2} = [1 3];
trellis{3} = [2 4];
trellis{4} = [2 4];

branchweight{1,1} = [0 0];   branchweight{1,2} = [1 1];
branchweight{2,1} = [1 1];   branchweight{2,2} = [0 0];
branchweight{3,1} = [1 0];   branchweight{3,2} = [0 1];
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 = 1e7 + 1;   % Number of message bits to generate
wD1 = 0;  wD2 = 0;   % Initial state of shift register
xvec = [];  yvec = [];

% Probability of bit error for BSC

if(0)
pvec = logspace(-3, 0, 8);    % Vector of p values
zi = find(pvec < 0.5);         % Keep only the values < 0.5
pvec = pvec(zi);
Np = length(pvec);             % Number of different p values
Np = Np+1;
pvec(Np) = 0.5;                % Add p=0.5 to end of vector
end

pvec = [.001, .003, .007, .02, .05, .1, .5];
Np = length(pvec);

% Matrix to decode states to message bits

sd = [ 0, -1,  1, -1;     % (i,j) element tells required input for
       0, -1,  1, -1;     % transition in state diagram from i to j.
      -1,  0, -1,  1; 
      -1,  0, -1,  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

dfree = 5;         % Free distance for this convolutional code
Npfine = length(pfine);
for ip=1:Npfine
p = pfine(ip);

done = 0; 
k = dfree; 
Pb(ip) = 0;
while (~done)
  Pk = 0;
  if (mod(k,2))    % k is odd
    for ee = (k+1)/2:k
      Pk = Pk + nchoosek(k,ee)*p^ee*(1-p)^(k-ee);
    end
  else             % k is even
    for ee = k/2:k
      Pk = Pk + nchoosek(k,ee)*p^ee*(1-p)^(k-ee);
    end
    Pk = Pk + 0.5*nchoosek(k,k/2)*p^(k/2)*(1-p)^(k/2);
  end
  ck = (k-4)*2^(k-5);
  dPb = ck*Pk;
  Pb(ip) = Pb(ip) + dPb;
% if (dPb < 1e-8)
  if (k >= 10)
    done = 1;
  else
    k = k+1;
  end    
end

end  % for loop ip

figure(1)
loglog(pvec, BER, 'bo:', pfine, Pb, 'b-', pfine, PeU, 'r-', ...
  pfine, PeL, 'g-', pfine, pfine, 'k-', pfine, PeRep3, 'c--');
legend('CONV. CODE SIMULATION', 'CONV. CODE U.B.', '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-8 1])
% print -depsc fig1_7.eps

figure(2)
loglog(pfine, Pb, 'b-', pfine, PeU, 'r-', ...
  pfine, PeL, 'g-', pfine, pfine, 'k-', pfine, PeRep3, 'c--');
legend('CONV. CODE U.B.', '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-8 1])
% print -depsc fig_theo.eps

save rjkvit2_out.mat

