%bch decoder using Berlekamp-Massey and Chien search
%standard bch decoder simulation
%by: yumi
%10/21/2005
close all; clear all;

m=12;
n=2^m-1;
t=2;
k=2^m-1-m*t;
simplified = 1;

%the constants used in the algorithm
alpha = gf(2, m);
zero = gf(0, m);
one = gf(1, m);

%since all zero is always a codeword, only have to specify the error
%code = zeros(1, n);
code=gf(zeros(1, n), 1);

code(1)=gf(1, 1);
code(2)=gf(1, 1);

codev = code.x;
reccode=gf(codev, m);

%creating alpha array
%note that syndrome should be in the order [s3, s2, s1, s0]
alpha_tb=gf(zeros(1, 2*t), m);
for i=1:2*t,
    alpha_tb(i)=alpha^(2*t-i+1);
end;

%imba
lambda = gf([1, zeros(1, t)], m);
lambda0= lambda;
b=gf([0, 1, zeros(1, t)], m);
b2 = gf([0, 0, 1, zeros(1, t)], m);
k=0;
gamma = one;
delta = zero;
syndrome_array = gf(zeros(1, t+1), m);

if(simplified == 1) 
    for r=1:t,
        r1 = 2*t-2*r+2;
        r2 = min(r1+t, 2*t);
        num = r2-r1+1;
        syndrome_array(1: num) = syndrome(r1:r2);
        delta = syndrome_array*lambda';
        
        lambda0 = lambda;
        lambda = gamma*lambda-delta*b2(2:t+2);
        
        if((delta~= zero) && (k>=0))
            b2(3)=zero;
            b2(4:3+t) = lambda0(1:t);
            gamma = delta;
            k = -k;
        else
            b2(3:3+t) = b2(1:t+1);
            gamma = gamma;
            k=k+2;
        end
        joke=1;
    end
else
    for r=1:2*t,
        r1 = 2*t-r+1;
        r2 = min(r1+t, 2*t);
        num = r2-r1+1;
        syndrome_array(1:num) = syndrome(r1:r2);
        delta = syndrome_array*lambda';

        lambda0 = lambda;
        lambda = gamma*lambda-delta*b(1:t+1);

        if((delta ~= zero) && (k>=0))
            b(2:2+t)=lambda0;
            gamma = delta;
            k=-k-1;
        else
            b(2:2+t) = b(1:t+1);
            gamma = gamma;
            k=k+1;
        end   
        joke=1;
    end
end

%inverstable
inverse_tb = gf(zeros(1, t+1), m);
for i=1:t+1,
    inverse_tb(i) = alpha^(-i+1);
end;

%chien's search
lambda_v = zero;
accu_tb=gf(ones(1, t+1), m);
for i=1:n,
    lambda_v=lambda*accu_tb';
    accu_tb = accu_tb.*inverse_tb;
    if(lambda_v==zero)
        error(1,n-i+1)=1;
    else
        error(1,n-i+1)=0;
    end
end

found = find(error(1,:)~=0)