gfdeconv.m

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

function [q, r] = gfdeconv(b, a, p)
%GFDECONV GF(P) polynomial deconvolution or GF(P^M) elements dividing.
%       [Q, R] = GFDECONV(B, A) computes the quotient Q and remainder R of
%       deconvolution B by A in GF(2).
%
%       [Q, R] = BFDECONV(B, A, P) computes the quotient and remainder R of
%       deconvolution B by A in GF(P) when P is a scalar prime number.
%       When P is a matrix containing the tuple of all elements in GF(p^M),
%       this function takes A and B as indices (power number of the
%       exponential form) of GF(p^M) elements. The output Q is
%       alpha^C = alpha^A / alpha^B in GF(p^M). The computation is
%       element-by-element computation. You can generate the tuple of all
%       elements in GF(Q^M) by P = GFTUPLE([-1:Q^M-2]', M, Q). 
%   
%       In polynomial computation, A, B, and C are in ascending order, i.e.,
%       A = [a_0, a_1, a_2,..., a_(n-1), a_n] represents
%       A(X) = a_0 + a_1 X + a_2 X^2 +...+ a_(n-1) X^(n-1) + a_n X^n
%       a_i must be a element in GF(P).
%
%       In power representation form, [-Inf, 0, 1, 2, ...] represents 
%       [0, 1, alpha, alpha^2, ...] in GF(p^m).
%
%      See also GFADD, GFMUL, GFTUPLE, GFPRETTY

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

% input validation
if nargin < 2
        error('Not enough input for GFDECONV')
elseif nargin < 3
 p = 2;
end

a = gftrunc(a(:)', p);
b = gftrunc(b(:)', p);

len_b = length(b);
len_a = length(a);
len_p = length(p);
if len_a > len_b
    q = 0;
    if len_p > 1
        q = -Inf;
    end;
    r = b;
    return
end;

if length(p) == 1
    % computation in GF(p) field
    if ~isempty([find(a~=floor(a)), find(a<0), find(a>=p)])
      error('The polynomial coeficients must be in GF(P)')
    elseif ~isempty([find(b~=floor(b)), find(b<0), find(b>=p)])
     error('The polynomial coeficients must be in GF(P)')
    end;

    % make a(len_a) to be one.
    if a(len_a) ~= 1
          t = a(len_a);
           if isprime(p)
                   % by Fermat's theory (pp 96 in MacWilliams)
             % for any integer q less than p q^(p-1) = 1
             q = t^(p-2);
                    if  (rem(q * t, p) ~= 1)
                        disp('Warning: The computation may be not accurate because a large integer')
                            disp('is introduced in to the calculation.');
                   end;
            else
                    error('P must be a prime in GFFILTER')
          end
     % made a(len_a) == 1.
           a = rem(a * q, p)
    end;
    q = fliplr(gffilter(fliplr(b),fliplr(a),[1,zeros(1,length(b)-length(a))],p));
    r = gfadd(rem(p-b, p), gfconv(q, a, p), p, []);
else
    % computation in GF(p^m) field
    tp = p;
    p = max(max(tp)) + 1; % the prime number.
    [n_tp, m_tp] = size(tp);
    pow_dim = p^m_tp - 1;
    if pow_dim+1 ~= n_tp
        error('The given GF(P^M) M-tuple list is not a valid one');
    end;

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

    a = fliplr(a);
    b = fliplr(b);
    flag_a1 = 0;
    if a(1) ~= 0
        % have to make a(1) = 0
        factor = pow_dim - a(1);
        indxa = find(a<0);
        indxb = find(b<0);
        a = rem(a + factor, pow_dim);
        b = rem(b + factor, pow_dim);
        flag_a1 = 1;
        if ~isempty(indxa)
            a(indxa) = indxa - Inf;
        end;
        if ~isempty(indxb)
            b(indxb) = indxb - Inf;
        end;
    end;
    len_q = len_b - len_a + 1;
    for i = 1 : len_q
        if b(i) >= 0
            for k = 1 : len_a - 1
                % b(i+k) = gfadd(b(i+k), pow_dim - gfmul(b(i), a(k+1), tp), tp);
                % this line works for the binary tp only, if it is not binary tp
                % then, should be gfsub instead of gfadd.
                % tmp1 = gfmul(b(i), a(k+1), tp);
                % disp(['b(i)=', num2str(b(i)), ' a(k+1)=', num2str(a(k+1)),' tmp1=',num2str(tmp1)])
                %b(i+k) = gfplus(b(i+k), gfmul(b(i), a(k+1), tp), tp_num, tp_inv);
          if (b(i) < 0) | (a(k+1) < 0)
                    tmp = -1;
               else
                    tmp = rem(b(i) + a(k+1), pow_dim);
              end;
            b(i+k) = gfplus(b(i+k), tmp, tp_num, tp_inv);
            end;
        end;
    end;
    if len_a > 1
        r = gftrunc(fliplr(b(len_b-len_a+2 : len_b)), tp);
        if flag_a1
            r = rem(r + pow_dim-factor, pow_dim);
        end;
    else
        r = -Inf;
    end;
    q = fliplr(b(1:len_q));
end;
%--end of GFDIV--