% Name: MLdecoder_313code.m
% Maximum likelihood decoding of a binary (3,1,3) code.
% Copyright (c) 2007. Robert Morelos-Zaragoza. All rights reserved.
clear all

n = 3;                              % Length
k = 1;                              % Dimension
dmin = 3;                           % Minimum distance
rate = k/n;                         % Rate

G = [ 1 1 1 ];                      % Generator matrix

H = [ 1 1 0                         % Parity-check matrix
      1 0 1 ];                      % Integers 3 2 1

syntable = [ 0 0 0                  % Decoding look-up table
             0 0 1
             0 1 0
             1 0 0];

fprintf('   Simulation of HD and ML decoding of a binary (3,1,3) code\n');
fprintf('         with BPSK modulation over an AWGN channel\n');
fprintf('    (NOTE: This program may take a few minutes to finish)\n');

% Weight distribution W(C):
W = zeros(1,n+1);
W(1)=1; W(4)=1;

% The codewords and their modulated versions
for j=1:n
    v(1,j)=0; m(1,j)= 1;
    v(2,j)=1; m(2,j)=-1;
end

% Print the weight distribution
fprintf('\nWeight distribution: W(C)={1');
for i=2:n+1
    fprintf(',%d',W(i))
end
fprintf('}\n\n');

nsim = 10000;                          % Number of codewords per SNR value
init = 0;                               % Inital value of Eb/No (dB)
inc = 0.5;                              % Increment in Eb/No (dB)
final = 6;                              % Final value of Eb/No (dB)
num=0;

seed=123445678;
rand('state',seed);
randn('state',seed);

fprintf('Eb/No(dB) \tBER_HD   \tErr_HD \t    BER_ML  \tErr_ML \n');

for SNRdB = init:inc:final
    
    error_HD = 0;
    error_ML = 0;
    TotalNumError = 0;
    TotalN = 0;
    
    SNRdBs = SNRdB + 10*log(rate)/log(10);
    No = (10.^(-SNRdBs/10));
    Var = No/2;
    
    for Ns = 1:nsim
        
        msg = randint(1,k,[0,1]);       % Random binary message vector
        codeword = mod(msg*G,2);        % Compute the codeword in C
        x = (-1).^codeword;             % BPSK mapping
        rec = x + sqrt(Var)*randn(1,n); % Add AWGN to obtain received vector
        
        rh = (1-sign(rec))/2;           % Hard decision vector
        
        % Hard-decision decoding based on standard array (LUT)
        syn = bin2dec(num2str(mod(rh*H',2)));       % Syndrome as integer 
        codehard = mod(rh+syntable(syn+1,:),2);     % Decoded codeword
        
        % ML decoding
        [cor,I] = max(rec*m');          % I: Index of maximum correlation
        
        % Compute number of decoding information errors and update sum
        error_HD = error_HD + sum(xor(codeword(1:k),codehard(1:k)));
        error_ML = error_ML + sum(xor(v(I,1:k),codeword(1:k)));
        
    end % for Ns
    
    num = num + 1;
    ber_HD(num) = error_HD/(nsim*k);
    ber_ML(num) = error_ML/(nsim*k);
    snr(num) = SNRdB;
    
    fprintf('%5.2f\t%e %6.0f\t%e %6.0f\n', snr(num), ber_HD(num), ...
        error_HD, ber_ML(num), error_ML);

end % for SNRdB

% Compute the union bound for ML decoding
EbNo = 0:0.5:6; EbNoratio = 10.^(EbNo/10);
bound = 0;
for i=1:n
    bound = bound + i*W(i+1)/n * Q(sqrt(2*i*(k/n)*EbNoratio));
end
% and plot results
semilogy(snr,ber_HD,'-r^'), hold on, semilogy(snr,ber_ML,'-bs')
semilogy(EbNo,bound,'-bo'), axis tight, grid on
xlabel('E_b/N_0 (dB)'), ylabel('P_e(C)'), legend('HD sim','ML sim','Union bound')