errlocp.m

数字通信第四版原书的例程

function [sigma, err] = errlocp(syndrome, t, tp, pow_dim, err, type_flag)
%ERRLOCP Computes the error-location polynomial.
%       [SIGMA, ERR] = ERRLOCP(SYNDROME, T, TP, ERR, FLAG) computes the
%       error-location polynomial from length 2*T input vector SYNDROME.
%       T is the error-correction capability. Tp is the complete tuple list
%       of all elements in GF(2^M). POW_DIM = 2^M - 1, which is the dimension
%       of GF(2^M).  ERR contains the error information in the computation.
%       FLAG indicates the method to be used in the computation.
%       FLAG = 1, uses normal method.
%       FLAG = 0, use simplified method. The simplified method can be
%                      used for binary code only. 
%
%       Warning: This function does not provide error-check. Make sure
%                the parameters are signed correctly.

%       Wes Wang 10/5/95
%       Copyright (c) 1995-96 by The MathWorks, Inc.
%       $Revision: 1.1 $  $Date: 1996/04/01 17:57:28 $


[tp_n, tp_m] = size(tp);
tp_num = tp * 2.^[0 : tp_m-1]';
tp_inv(tp_num+1) = 0:pow_dim;

if type_flag 
    % use the berlekamp's algorithm
    t = 2 * t;
    mu = [-1, 0:t]';
    sigma_mu = [zeros(t+2,1), -ones(t+2, t)];
    d_mu = [0; syndrome(1); zeros(t, 1)];
    l_mu = [0 0 [1:t]]';
    mu_l_mu = mu - l_mu;

    % iteratiev start with row three. The first two rows are filled.
    for de_i = 3:t+2
        % no more effort to failed situation
        if (d_mu(de_i - 1) < 0) | err
            sigma_mu(de_i, :) = sigma_mu(de_i-1, :);
            l_mu(de_i) = l_mu(de_i - 1);
        else
            % find another row proceeding to row de_i -1
            % d_mu equals to zero
            % and mu - l_mu is the largest.
            indx = find(d_mu(1:de_i - 2) >= 0);
            rho  = find(mu_l_mu(indx) == max(mu_l_mu(indx)));
            rho = indx(rho(1));

            % by (6.25)
            % shifted = gfmul(d_mu(de_i - 1), pow_dim - d_mu(rho), tp);
            % shifted = gfmul(shifted, sigma_mu(rho, :), tp)';
            % multiply inreplace the above two lines.
            shifted = -ones(1, t+1);
            if (d_mu(de_i - 1) >= 0) & (pow_dim - d_mu(rho) >= 0)
                tmp = rem(pow_dim - d_mu(rho) + d_mu(de_i - 1), pow_dim);
                indx = find(sigma_mu(rho,:) >= 0);
                for de_k = 1 : length(indx)
                    shifted(indx(de_k)) = rem(tmp + sigma_mu(rho, indx(de_k)), pow_dim);
                end;            
            end;
            % end multiply
        
            shifting = (mu(de_i - 1) - mu(rho));

            % calculate new sigma_mu
            if ~isempty(find(shifted(t-shifting+2 : t+1) >= 0))
                % more than t errors, BCH code fails.
                err = 1;
            else
                % calculate the new sigma
                shifted = [-ones(1, shifting) shifted(1:t-shifting+1)];
                sigma_mu(de_i, :) = gfplus(sigma_mu(de_i-1,:), shifted, tp_num, tp_inv);
            end;
            l_mu(de_i) = max(l_mu(de_i-1), l_mu(rho) + (de_i - 1) - rho);
        end;
    
        % calculate d_mu. It is not necessary to do so if mu(de_i) == t
        if de_i < t+2
            % the constant term
            d_mu(de_i) = syndrome(mu(de_i) + 1);
            indx = find(sigma_mu(de_i, 2:t) >= 0);
    
            for de_j = 1 : length(indx)
                de_j_tmp = indx(de_j);
                % Before the "end", it is equivalent to 
                % d_mu(de_i) = gfadd(d_mu(de_i), ...
                %              gfmul(sigma_mu(de_i + 1, de_j_tmp+1), ...
                %              syndrome(mu(de_i) * 2 - de_j_tmp + 1), tp), tp);
                tmp = syndrome(mu(de_i) - de_j_tmp + 1);
                if (tmp < 0) | (sigma_mu(de_i, de_j_tmp + 1) < 0)
                    tmp = -1;
                else
                    tmp = rem(tmp + sigma_mu(de_i, de_j_tmp + 1), pow_dim);
                end;
                d_mu(de_i) = gfplus(d_mu(de_i), tmp, tp_num, tp_inv);
            end;
        end;

        % calculate 2*mu-l_mu
        mu_l_mu(de_i) = mu(de_i) - l_mu(de_i); 
    end;
else
    % use simplified algorithm
    mu = [-1/2, 0:t]';
    sigma_mu = [zeros(t+2,1), -ones(t+2, t)];
    d_mu = [0; syndrome(1); zeros(t, 1)];
    l_mu = [0 0 2*(1:t)]';
    mu2_l_mu = 2*mu - l_mu;
    % iteratiev start with row three. The first two rows are filled.
    for de_i = 3:t+2
        % no more effort to failed situation
        if (d_mu(de_i - 1) < 0) | err
            sigma_mu(de_i, :) = sigma_mu(de_i-1, :);
        else
            % find another row proceeding to row de_i -1
            % d_mu equals to zero
            % and 2*mu - l_mu is the largest.
            indx = find(d_mu(1:de_i - 2) >= 0);
            rho  = find(mu2_l_mu(indx) == max(mu2_l_mu(indx)));
            rho = indx(rho(length(rho)));

            % by (6.28)
            % shifted = gfmul(d_mu(de_i - 1), pow_dim - d_mu(rho), tp);
            % shifted = gfmul(shifted, sigma_mu(rho, :), tp)';
            % multiply inreplace the above two lines.
            shifted = -ones(1, t+1);
            if (d_mu(de_i - 1) >= 0) & (pow_dim - d_mu(rho) >= 0)
                tmp = rem(pow_dim - d_mu(rho) + d_mu(de_i - 1), pow_dim);
                indx = find(sigma_mu(rho,:) >= 0);
                for de_k = 1 : length(indx)
                    shifted(indx(de_k)) = rem(tmp + sigma_mu(rho, indx(de_k)), pow_dim);
                end;            
            end;
            % end multiply
        
            shifting = (mu(de_i - 1) - mu(rho)) * 2;

            % calculate new sigma_mu
            if ~isempty(find(sigma_mu(1:de_i-2,max(1,t-shifting+2):t+1) >=0))
                disp('potential error more than posible')
                disp(['de_i=',num2str(de_i), ' shifting=',num2str(shifting)]);
                disp(sigma_mu)
                disp(['shifted=',mat2str(shifted)]);
            end;
            if ~isempty(find(shifted(max(t-shifting+2, 1) : t+1) >= 0))
                % more than t errors, BCH code fails.
                err = 1;
            else
                % calculate the new sigma
                shifted = [-ones(1, shifting) shifted(1:t-shifting+1)];
                sigma_mu(de_i, :) = gfplus(sigma_mu(de_i-1,:), shifted, tp_num, tp_inv);
            end;
        end;
        l_mu(de_i) = max(find(sigma_mu(de_i,:) >= 0)) - 1;
    
        % calculate d_mu. It is not necessary to do so if mu(de_i) == t
        if de_i < t+2
            % the constant term
            d_mu(de_i) = syndrome(mu(de_i) * 2 + 1);
            indx = find(sigma_mu(de_i, 2:t) >= 0);
    
            for de_j = 1 : length(indx)
                de_j_tmp = indx(de_j);
                % Before the "end", it is equivalent to 
                % d_mu(de_i) = gfadd(d_mu(de_i), ...
                %              gfmul(sigma_mu(de_i, de_j_tmp+1), ...
                %              syndrome(mu(de_i) * 2 - de_j_tmp + 1), tp), tp);
                tmp = syndrome(mu(de_i) * 2 - de_j_tmp + 1);
                if (tmp < 0) | (sigma_mu(de_i, de_j_tmp + 1) < 0)
                    tmp = -1;
                else
                    tmp = rem(tmp + sigma_mu(de_i, de_j_tmp + 1), pow_dim);
                end;
                d_mu(de_i) = gfplus(d_mu(de_i), tmp, tp_num, tp_inv);
            end;
        end;

        % calculate 2*mu-l_mu
        mu2_l_mu(de_i) = mu(de_i) * 2 - l_mu(de_i); 
    end;
end;

% the error polynomial        
sigma = sigma_mu(t+2, :);
% truncate the redudancy
indx = find(sigma >= 0);
sigma = sigma(1:max(indx));
% completed constructing error polynomial

%-- end of errlocp --