bchpoly.m

经典通信系统仿真书籍《通信系统与 MATLAB (Proakis)》源代码

        end;
        if x_a > 2
            plot([3.5,3.5],[-1,hght+2],   'b-');
        end;

        lengt = length(setstr(n));
    end;

    ind = find(n_lst(:,1) < dim);
    if ~isempty(ind)
        base_i = sum(n_lst(ind, 3));
    else
        base_i = 0;
    end;

    for i = 1:n_cs-2
        if (nargout < 1)
             xx = num2str(n - sum( sum( ~isnan( cs(1 : i+1, :))))+1);
            xx = [setstr(ones(1, lengt-length(xx))*32) xx];
            if (x_a <= 1) | ((x_a <= 2) & (i <= hght)) | ((x_a <= 3) & (i<=hght))
                % the first column
                xx=text(0, hght-i, xx);
            elseif (x_a<=2) | ((x_a<=3) & (i<=2*hght))
                % the second colum
                xx=text(4-x_a, 2*hght-i, xx);
            else
                if (i <= 3*hght)
                    % the third column
                    xx=text(2, 3*hght-i, xx);
                elseif (i <= 4*hght)
                    % the forth column
                    xx=text(3, 4*hght-i, xx);
                elseif (i <= 5*hght)
                    % the fifth column
                    xx=text(4, 5*hght-i, xx);
                else
                    % the sixth column
                    xx=text(5, 6*hght-i, xx);
                end;
            end;
            set(xx,'Color',[1 0 1]);
            if dim <= 10
                xx = num2str(t_lst(i+base_i));
                xx = [setstr(ones(1, 4-length(xx))*32) xx];
                if (x_a<=1) | ((x_a<=2) & (i*2<=n_cs)) | ((x_a<=3) & (i*3<=n_cs))
                    % the first column
                    xx=text(1, hght-i, xx);
                elseif (x_a<=2) | ((x_a<=3) & (i*3<=2*n_cs))
                    % the second colum
                    xx=text(3, 2*hght-i, xx);
                else
                    % the third colum
                    xx=text(5, 3*hght-i, xx);
                end;
                set(xx,'Color',[0 1 1]);
            end;
            usda = [length(get(x, 'child')), [10, 10, 240*x_a, hght*12+100]];
            set(h, 'UserData', usda);
        else
            if dim <= 10
                pg = [pg;
                  [n, ...
                   n - sum( sum( ~isnan( cs(1 : i+1, :))))+1,...
                   t_lst(i+base_i)]];
            else
                pg = [pg;
                  [n, ...
                   n - sum( sum( ~isnan( cs(1 : i+1, :))))+1]];
            end;
        end;
               
    end;
    return;
end;

if (dim<1) | (k<1) | (floor(dim) ~= dim) | (floor(k) ~= k)
    error('Input variable for BCHPOLY is not valid')
end

if df_flag
    m = gfprimdf(dim);
end;

% compute the minmum polynomial
if k == n - dim
    % It is a trivial case, simply assign the primitive polynomial 
    % as BCH polynomial.
    pm = m;
    pg = m;
    t = 1;
    if (nargout > 2)
        cs = gfcosets(dim);
    end;
    n_i = 1;
else
    % list cosets and taking off the first row.
    cs = gfcosets(dim);
    [n_cs, m_cs] = size(cs);

    pl = gfminpol(cs(2:n_cs,1), m);

    n_terms = 0;
    n_i = 0;
    pm = [];
    while ((n_terms < (n - k)) & (n_i+1 < n_cs))
        n_i = n_i + 1;
        n_terms = n_terms + sum(~isnan( cs(n_i+1, :)));
        pm = [pm; pl(n_i, :)];
    end;

    if ((n_i+1 >= n_cs) & (n_terms < (n - k)))
        disp(['Warning: The maximum message length is ', num2str(n-n_terms)])
    end;

    pg = gftrunc(pm(1,:));
    [n_pm, m_pm] = size(pm);
    if n_pm > 1
        for i = 2 : n_pm
            pg = gfconv(pg, gftrunc(pm(i,:)));
        end;
    end;
end;

% construct the parity-check matrix
if nargout > 3
    tp = gftuple([-1:n-1]', m);
    h = tp(2:n+1,:)';
    if n_i > 1
        for i = 2 : n_i

            % build the rows leaded by new coset.
            ad = [];
            seed = cs(i+1,1);
            for j = 0:n-1
                ad = [ad tp(gfadd(-Inf, j*seed, tp)+2,:)'];
            end;

            %cleasr the linear related terms in the new rows.
            for j = dim : -1 : 1
                % clean the all zeros terms
                if max(ad(j, :)) == 0
                    ad(j,:) = [];
                else
                    % remove the repeat rows.
                    run_flag = 1;
                    jj = j - 1;
                    while run_flag & (jj > 0)
                        if max(rem(ad(j,:)+ad(jj, :), 2)) == 0
                            ad(jj,:)=[];
                            run_flag = 0;
                        end;
                        jj = jj - 1;
                    end;
                end;
            end;

            % add the new rows to the parity check matrix.
            h = [h; ad];
        end;
    end;
end;

if (~df_flag) & (dim <= 10)
    if (max(rem(m - gfprimdf(dim), 2)) == 0)
        df_flag == 0;
    end;
end;

% detect erorr-correction capability.
if df_flag & (dim <= 10)
    % look up the data base
    ind = find(n_lst(:,1) < dim);
    if ~isempty(ind)
        base_i = sum(n_lst(ind, 3));
    else
        base_i = 0;
    end;
    t = t_lst(base_i + n_i);
else
    % designed distance. By Corollary 9.2 in Peterson & Weldon
    %    t = floor((n - k) / dim);
    t = n_i;
    %    dis_min = 2 * n_i + 1;  % the minimum by theory.
    tt = 1;
    for i = 1:t+1
        tt = tt * i;
    end;

    % by Farr Theorem (Theorem 9.6) in Peterson & Weldon, that's the solution.
    if (dim > 1 + log(tt) / log(2)) & (t < dim/2)
        return;
    end;

    if dim > 3
        disp('It takes a long time to processs the minimum distance calculation');
    end;

    dis_min = 2 * t + 1;

    % The maximum distance will not more than the following value
    % by Theory 9.8 of Peterson and Weldon.
    dis = 2 * dis_min-1;
    h = h';
    [x1, x2] = size(h);
    dis = min(dis, x2);
    ind = [0 1 zeros(1, x1 - 2)];
    i = 3;
    px1 = 2^x1 - 1;
    tt = floor((dis - 1) / 2);
    fprintf(['The error-correction capability ',num2str(t),' and ',num2str(tt),'.\n']);
    fprintf([num2str(px1), ' number to be processed. \nPlease wait....']);
    count = 0;
    while (i < px1) & (tt > t)
        n_j = 1;
        i = i + 1;
        n_k = 1;
        while n_j
            if ind(n_k)
                ind(n_k) = 0;
                n_k = n_k + 1;
            else
                ind(n_k) = 1;
                n_j = 0;
            end;
        end;
        sum_ind = sum(ind);
        if (sum_ind >= dis_min)
            if (sum_ind < dis)
                count = count + 1;
                if count > 1000
                    fprintf('.');
                    count = 0;
                end;
                idx = find(ind == 1);
                if max(rem(sum(h(idx, :)), 2)) == 0
                    dis = sum_ind;
                    tt = floor((dis-1)/2);
                    if tt ~= t
                        fprintf(['\nThe error-correction capability is, now, between ', num2str(t), ' and ',num2str(tt), '\n']);
                    end;
                end;
            end;
        end;
    end;
    fprintf('..done.\n');
    fprintf(['The error-correction capability is found to be ',num2str(tt), '\n']);
    t = floor((dis - 1) / 2);
    h = h';
end;

%--end of bchpoly.m--