clear;

G = [1 0 1 0 0 1;
     0 1 1 0 1 0;
     1 1 0 1 0 0];
 
H = [1 0 0 1 0 1;
     0 1 0 1 1 0;
     0 0 1 0 1 1];

C = [3 2 5 1 6 4];
 
r = [0 0 1 1 1 0];

sigma = 0.1;
count = 1000;

x = [];
y = [];

for sigma = 0:0.05:1

    errors = 0;
  
    for i = 1:count
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Create the message
        n = floor(rand()*(2^3));
        m = bitget(n,3:-1:1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Encode the message
        U = (m*G == 1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Add Errors
        e = rand(1,6) <= sigma;
        r = (U+e == 1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Decode Message
        
        %Get syndrome
        s = (r*H' == 1);

        % Check syndrome
        S = sum(s);
        if S > 0 
            if S < 3
                % Correct errors
                c = C(sum(bitset(0,[3 2 1], s)));
                r(c) = (r(c) + 1 == 1);
            else
                % Too many errors, ignore
                errors = errors + 1;
                continue;
            end
        end
 
        % Decode
        d = r(6:-1:4);
       
        %%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Compare decoded to original
        if ~isequal(m,d)
            errors = errors + 1;
        end
    end

    x = [x sigma];
    y = [y errors/count];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot BER
loglog(x,y)
axis([0,1,0,1])
set(gca,'xtick',[0:0.1:1])
set(gca,'ytick',[0:0.1:1])
title('BER for (6,3) Hamming Code')
grid on