% Hamming(7,4) simulation
%
% Rich Kozick, October 2001

k = 4;    % Number of message bits per block
n = 7;    % Number of codeword bits per block

% Enter codeword table

xtable = [0, 0, 0, 0, 0, 0, 0; ...
          1, 1, 0, 1, 0, 0, 0; ...
          0, 1, 1, 0, 1, 0, 0; ...
          1, 0, 1, 1, 1, 0, 0; ...
          1, 1, 1, 0, 0, 1, 0; ...
          0, 0, 1, 1, 0, 1, 0; ...
          1, 0, 0, 0, 1, 1, 0; ...
          0, 1, 0, 1, 1, 1, 0; ...
          1, 0, 1, 0, 0, 0, 1; ...
          0, 1, 1, 1, 0, 0, 1; ...
          1, 1, 0, 0, 1, 0, 1; ...
          0, 0, 0, 1, 1, 0, 1; ...
          0, 1, 0, 0, 0, 1, 1; ...
          1, 0, 0, 1, 0, 1, 1; ...
          0, 0, 1, 0, 1, 1, 1; ...
          1, 1, 1, 1, 1, 1, 1];

% Probability of bit error for BSC

pvec = logspace(-3, 0, 12);    % 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

RUNs = 1e6;

for pind = 1:Np
p = pvec(pind);

bit_error_count = 0;

for run=1:RUNs

% Generate a block of 4 message bits
% (equally likely 0's and 1's)

z = unifrnd(0, 1, 1, 4);
w = round(z);

% Get the row index into xtable by binary-to-decimal conv.:

m = w(1) + w(2)*2 + w(3)*4 + w(4)*8;
x = xtable(m+1, :);

% Simulate errors of all 7 bits through BSC:

z = unifrnd(0, 1, 1, 7);
zi = find(z <= p);         % Positions where errors occur

% Error vector:

e = zeros(1, 7);
e(zi) = ones(size(zi));

% Received vector at output of BSC:

y = xor(x, e);

% Find codeword with minimum Hamming distance:

for mi=1:16
  dH = sum(xor(y,xtable(mi,:)));
  dHr(mi) = dH;
  if ( (dH == 0) | (dH == 1) )
    xhat = xtable(mi,:);
    break;
  end
end
what = xhat(4:7);
bit_error_count = bit_error_count + sum(xor(w,what));

end   % RUNs loop

BER(pind) = bit_error_count/(k*RUNs);

% fprintf('p = %.2f: %d bit errors in %d bits => BER = %f\n', p, ...
% bit_error_count, k*RUNs, BER(pind));

end   % p loop

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

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

figure(1)
loglog(pvec, BER, 'bo:', pfine, PeU, 'r-', pfine, PeL, 'y-', ...
pfine, pfine, 'k-')
legend('SIMULATION', 'UPPER BOUND', 'LOWER BOUND', 'UNCODED')
xlabel('ERROR PROBABILITY FOR BSC (p)')
ylabel('BER')
title('BER FOR HAMMING(7,4) CODE OVER BSC')
axis([1e-4 0.5 1e-7 1])

% Also add 3 times repetition code result
figure(2)
loglog(pvec, BER, 'go:', pfine, PeU, 'r-', pfine, PeL, 'b-', ...
pfine, pfine, 'k-', pfine, PeRep3, 'c--')
legend('SIMULATION', 'UPPER BOUND', 'LOWER BOUND', 'UNCODED', ...
'3 \times REPETITION CODE')
xlabel('ERROR PROBABILITY FOR BSC (p)')
ylabel('BER')
title('BER FOR HAMMING(7,4) CODE OVER BSC')
axis([1e-4 0.5 1e-7 1])

save ham74_out6.mat

