% File: bch_example.m
% Description: Generates the Galois field GF(2^m) with an arbitrary
%              primitive polynomial and illustrates the decoding of
%              binary primitive BCH codes. Errors-only decoding performed
%              by direct solution of two sets of linear equations, or
%              Peterson-Gorenstein-Zierler (PGZ) decoder.
%
% Copyright (c) 2007. Robert Morelos-Zaragoza. All rights reserved.
%
clear all
p=2; m=4;                               % Parameters of GF(p^m)
n=15; k=7; d=5;                         % BCH code parameters
b=1;                                    % log(initial zero)

% -----------------------------------------------------------------------
%               Generate the Galois field, GF(2^m)
% -----------------------------------------------------------------------
% Primitive polynomial: x^4 + x + 1
pol=[1 1 0 0 1];
pold = 19;  % In decimal
% You can also find a primitive polynomial of degree m over GF(2) with
% pol=gfprimfd(m,'min',2)
% Compute the nonzero elements
n=2^m-1;
for i=1:n;
    field(i,:) = gftuple(i-1,pol,2);
end
% Print the table of GF(2^m)
fprintf('   - '); for j=1:m, fprintf('0'); end, fprintf('\n');
for i=1:n
    fprintf('%4d ',i-1);
    for j=1:m
        fprintf('%d',field(i,j));
    end
    fprintf('\n');
end
fprintf('\n\n');
% Primitive element in GF(2^m)
alpha = gf(2,m,pold);

% -----------------------------------------------------------------------
%                        Received polynomial
% -----------------------------------------------------------------------
rec_pol = gf([1 0 0 0 0 0 1 1 0 0 0 0 0 0 0],m,pold);

% -----------------------------------------------------------------------
%                       Compute the syndromes
% -----------------------------------------------------------------------
s = gf([zeros(1,n-k)],m,pold);
for i=b:b+d-2
    ev = gf([zeros(1,n)],m,pold);  val = alpha^i;
    for j=1:n
        ev(j) = val^(j-1);
    end
    s(i-b+1) =rec_pol*ev';
end
for i=1:d-1, fprintf('s(%d)=alpha^',i);
    if      s(i) ~=0,  fprintf('%d\n',log(s(i)));
    else    fprintf('-Inf\n'); end
end

% -----------------------------------------------------------------------
%                           PGZ decoder
% -----------------------------------------------------------------------
done = 0;
t = floor(d/2);

while (~done)
    % ===================================================================
    %   STEP 1: Check determinant of t-by-t matrix A in key equation
    % ===================================================================
    A = gf([zeros(t,t)],m,pold);
    for i1=1:t
        for i2=1:t
            A(i1,i2) = s(i1+i2-1);
        end
    end
    delta = det(A);
    if (det(A) ~= 0)
        done = 1;
        % ===============================================================
        % Step 2: If determinant is nonzero, t errors occurred. Find the
        % coefficients of sigma, the error locator polynomial, by solving
        % the matrix equation A x = B
        % ===============================================================
        B = gf([zeros(t,1)],m,pold);
        sigma = gf([zeros(1,t)],m,pold);
        for i=1:t
            B(i) = s(t+i);
        end
        sigma = (A \ B)'; % Transpose of solution of A * x = B
        % Constant coefficient of sigma, always = 1
        sigma(t+1) = 1;
        % ===============================================================
        %                   Find the error positions
        % ===============================================================
        % Roots of sigma
        r = roots(sigma);
        % Error positions are the inverses of the roots
        fprintf('\nError locations: \n')
        position = 1./r
        % ===============================================================
        %                    Find the error values
        % ===============================================================
        for i1=1:t, for i2=1:t
                A(i1,i2) = position(i2)^(i1+b-1);
            end, end
        for i=1:t, B(i)=s(i); end
        value = (A \ B);    % Solution of A * x = B
        fprintf('Error values: \n')
        value
        % ===============================================================
        %                Correct the received polynomial
        % ===============================================================
        est_error = gf([zeros(1,n)],m,pold);
        est_code_pol = gf([zeros(1,n)],m,pold);
        for i=1:t, for j=1:n
                val = alpha^(j-1);
                if val == position(i)
                    est_error(j) = est_error(j) + value(i);
                end
        end, end
        fprintf('Estimated codeword:\n');
        est_code_pol = rec_pol + est_error
    else
        % ===============================================================
        % Else: Rank is less than t. Decrease the value of t by one
        % ===============================================================
        t = t-1;
        if t == 0
            done = 1;
        end
    end
    if t == 0
        fprintf('Errors present but not correctable\n');
    end
end