gffilter.m

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

function y = gffilter(b, a, x, p)
%GFFILTER GF(P) polynomial filter computation.
%       Y = GFFILTER(B, A, X) filters the data in vector X with the
%       filter described by vectors A and B to create the filtered
%       data Y in GF(2).
%
%       Y = GFFILTER(B, A, X, P) filters the data in GF(P).
%
%       The filter uses the standard difference equation:
%       y(n) = b_nb*x(n) + b_(nb-1)*x(n-1) + ... + b_0*x(n-nb)
%                  - a_(na-1)*y(n-1) - ... - a_0*y(n-na)
%
%       Different from all of the other GF functions in this toolbox, this
%       this function uses descending ordered polynomial, i.e.
%         A = [a_n, a_(n-1), ..., a_2, a_1, a_0] represents
%         A(X) = a_n X^n + ... + a_2 X^2 + a_1 X + a_0
%       a_i must be a non-negative integer less than P.
%
%       See also GFADD, GFMUL, GFDIV, GFTUPLE

%       Wes Wang 6/8/94, 10/7/95
%       Copyright (c) 1995-96 by The MathWorks, Inc.
%       $Revision: 1.2 $  $Date: 1996/04/25 15:57:10 $
%       One may directly use the built-in function filter.
%       However, the computation may not be accurate if
%       the computation order is high.

if nargin < 3
     error('Not enough input for GFFILTER')
elseif nargin < 4
 p = 2;
end

% variable assignment.
a = gftrunc(a(:)');
b = gftrunc(b(:)');
x = x(:)';
len_a = length(a);
len_b = length(b);
len_x = length(x);

% input check up
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)')
elseif ~isempty([find(x~=floor(x)), find(x<0), find(x>=p)])
 error('The input in Galois space must be in GF(P)');
end;

% make a(1) to be one.
if a(1) == 0
      error('First denominator filter coeficient must be non-zero');
elseif a(1) ~= 1
  t = a(1);
       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(1) == 1.
       a = rem(a * q, p)
end;

% computation length adjustment
xx = [zeros(1, len_b-1), x(:)'];
y = [zeros(1, len_a-1), zeros(1, len_x)];

% flip the computation variable
fa = fliplr(a(2:len_a));
fb = fliplr(b);

% difference equation iteration
%make fa = -fa:
fa = p - fa;
for i = 1:len_x,
    indx = i + len_a - 1;
    if len_a > 1
        y(indx) = fb * xx(i:len_b+i-1)' + fa * y(i:indx-1)';
    else
        y(indx) = fb * xx(i:len_b+i-1)';
    end;
    y(indx) = rem(y(indx), p);
end;

%remove the added zeros.
y = y(len_a:len_x+len_a-1);

%---end of GFFILTER---