% hardvit.m -- Hard decision decoding Viterbi algorithm.

% Describe the rate 1/2 convolutional code:
% g_1(D) = 1 + D + D^2  and  g_2(D) = 1 + D^2

% Define the paths in the trellis

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

% Define codewords along trellis branches

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

% Set the window length for the Viterbi algorithm:
% typically choose as 5 times the memory in the encoder

WindowLen = 10;

% Initialize the Viterbi algorithm.  dH.m is a function to compute the
% Hamming distance, which is the path metric that is minimized for
% hard decision decoding.

initvit1(trellis, cw, WindowLen, 'dH');
svec = [];     % Initialize state sequence through trellis to empty.
               % The output from the Viterbi algorithm will be svec.
               
T = 20;     % Number of message bits to send through channel
p = 4e-1;   % Probability of bit error for binary symmetric channel (BSC)
wD1 = 0;  wD2 = 0;   % Contents of shift register, initially 0.

% Simulate T message through BSC, and run Viterbi algorithm to decode.

for t = 1:T
    
    % Generate new message bit, w(t), at time t.
    
    z = rand;
    if (z <= 0.5)
        w(t) = 0;
    else
        w(t) = 1;
    end  
    
    % Find convolutionally encoded bits, x.
    
    %  x(1) = xor(xor(w(t), wD1), wD2);   % Dumb Kozick way
    x(2) = xor(w(t), wD2);
    x(1) = xor(x(2), wD1);
    
    % Update wD1 and wD2
    
     wD2 = wD1;
     wD1 = w(t);
    
    % Send x through the BSC, with output y.
    
    y = x;
    z = rand(1,2);
    zi = find(z < p);
    y(zi) = not(x(zi));
    
    % Feed y to the Viterbi algorithm for processing.
    
    s = viterbi1(y);    % returns state estimate WindowLen in the past,
    if (s)              % or s is empty if not enough bits through channel.
        svec = [svec, s];
    end
    
    disp([t, svec])
    
end  % Loop t=1:T

% Flush remaining paths through trellis

svec = [svec, vitflush(0)];

% Translate estimate of states, svec, to estimates of message bits
% state-to-message bit conversion matrix:

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

for t=1:T-1
    wh(t) = sd(svec(t), svec(t+1));
end

[w(1:T-1)', wh']
