% softvit2.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, Dec. 2001

% ** This version stores only WindowLen past message bits, and does
% soft-decision decoding over AWGN channels.

% 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];

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

T = 1e7 + 1;   % Number of message bits to generate

% Define Eb/N0 for AWGN channel:

Eb_over_N0_dB = [4];
snrvec = 10.^(Eb_over_N0_dB/10);

pvec = Qf(sqrt(2*snrvec));
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);
snr = snrvec(pind);
errcnt = 0;

% Define weights for soft decision decoding:

branchweight{1,1} = [-1 -1]*snr;   branchweight{1,2} = [1 1]*snr;
branchweight{2,1} = [1 1]*snr;     branchweight{2,2} = [-1 -1]*snr;
branchweight{3,1} = [1 -1]*snr;    branchweight{3,2} = [-1 1]*snr;
branchweight{4,1} = [-1 1]*snr;    branchweight{4,2} = [1 -1]*snr;

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

wD1 = 0;  wD2 = 0;   % Initial state of shift register
wvec = [];

% Initialize vectors to store sequences
wrec = [];
xrec = [];
yrec = [];
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

% 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 AWGN channel,
% and generate matched filter outputs, y

y = normrnd( (x-0.5)*2*snr, sqrt([1, 1]*snr/2) );;

% Apply Viterbi algorithm to y

s = viterbi1(y);  

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

% wrec = [wrec, w];
% xrec = [xrec, x];
% yrec = [yrec, y];

end  % Loop over time t=1:T

% Flush remaining paths, but do NOT count errors at end of block:

svec = vitflush(0); 
for t=1:WindowLen-1
  wh = sd(sD1, svec(t));
  sD1 = svec(t);
% whrec= [whrec, 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);
fprintf('p = %f: %d errors in %d message bits\n', p, errcnt, T-WindowLen);
BER(pind) = errcnt/(T-WindowLen);

end   % Loop over p values

save softvit2_out.mat
return

% 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 soft.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 softvit2_out.mat

