octave_packages/communications-1.1.1/bchpoly.m

   1 ## Copyright (C) 2003 David Bateman
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{p} = } bchpoly ()
  18 ## @deftypefnx {Function File} {@var{p} = } bchpoly (@var{n})
  19 ## @deftypefnx {Function File} {@var{p} = } bchpoly (@var{n},@var{k})
  20 ## @deftypefnx {Function File} {@var{p} = } bchpoly (@var{prim},@var{k})
  21 ## @deftypefnx {Function File} {@var{p} = } bchpoly (@var{n},@var{k},@var{prim})
  22 ## @deftypefnx {Function File} {@var{p} = } bchpoly (@var{...},@var{probe})
  23 ## @deftypefnx {Function File} {[@var{p},@var{f}] = } bchpoly (@var{...})
  24 ## @deftypefnx {Function File} {[@var{p},@var{f},@var{c}] = } bchpoly (@var{...})
  25 ## @deftypefnx {Function File} {[@var{p},@var{f},@var{c},@var{par}] = } bchpoly (@var{...})
  26 ## @deftypefnx {Function File} {[@var{p},@var{f},@var{c},@var{par},@var{t}] = } bchpoly (@var{...})
  27 ##
  28 ## Calculates the generator polynomials for a BCH coder. Called with no input
  29 ## arguments @dfn{bchpoly} returns a list of all of the valid BCH codes for
  30 ## the codeword length 7, 15, 31, 63, 127, 255 and 511. A three column matrix
  31 ## is returned with each row representing a seperate valid BCH code. The first
  32 ## column is the codeword length, the second the message length and the third
  33 ## the error correction capability of the code.
  34 ##
  35 ## Called with a single input argument, @dfn{bchpoly} returns the valid BCH
  36 ## codes for the specified codeword length @var{n}. The output format is the
  37 ## same as above.
  38 ##
  39 ## When called with two or more arguments, @dfn{bchpoly} calculates the
  40 ## generator polynomial of a particular BCH code. The generator polynomial
  41 ## is returned in @var{p} as a vector representation of a polynomial in
  42 ## GF(2). The terms of the polynomial are listed least-significant term
  43 ## first.
  44 ##
  45 ## The desired BCH code can be specified by its codeword length @var{n}
  46 ## and its message length @var{k}. Alternatively, the primitive polynomial
  47 ## over which to calculate the polynomial can be specified as @var{prim}.
  48 ## If a vector representation of the primitive polynomial is given, then
  49 ## @var{prim} can be specified as the first argument of two arguments,
  50 ## or as the third argument. However, if an integer representation of the
  51 ## primitive polynomial is used, then the primitive polynomial must be
  52 ## specified as the third argument.
  53 ##
  54 ## When called with two or more arguments, @dfn{bchpoly} can also return the
  55 ## factors @var{f} of the generator polynomial @var{p}, the cyclotomic coset
  56 ## for the Galois field over which the BCH code is calculated, the parity
  57 ## check matrix @var{par} and the error correction capability @var{t}. It
  58 ## should be noted that the parity check matrix is calculated with
  59 ## @dfn{cyclgen} and limitations in this function means that the parity
  60 ## check matrix is only available for codeword length upto 63. For
  61 ## codeword length longer than this @var{par} returns an empty matrix.
  62 ##
  63 ## With a string argument @var{probe} defined, the action of @dfn{bchpoly}
  64 ## is to calculate the error correcting capability of the BCH code defined
  65 ## by @var{n}, @var{k} and @var{prim} and return it in @var{p}. This is
  66 ## similar to a call to @dfn{bchpoly} with zero or one argument, except that
  67 ## only a single code is checked. Any string value for @var{probe} will
  68 ## force this action.
  69 ##
  70 ## In general the codeword length @var{n} can be expressed as
  71 ## @code{2^@var{m}-1}, where @var{m} is an integer. However, if
  72 ## [@var{n},@var{k}] is a valid BCH code, then a shortened BCH code of
  73 ## the form [@var{n}-@var{x},@var{k}-@var{x}] can be created with the
  74 ## same generator polynomial
  75 ##
  76 ## @end deftypefn
  77 ## @seealso{cyclpoly,encode,decode,cosets}
  78
  79 function [p, f, c, par, t] = bchpoly(nn, k, varargin)
  80
  81   if ((nargin < 0) || (nargin > 4))
  82     error ("bchpoly: incorrect number of arguments");
  83   endif
  84
  85   probe = 0;
  86   prim = 0;    ## Set to zero to use default primitive polynomial
  87   if (nargin == 0)
  88     m = [3:9];
  89     n = 2.^m - 1;
  90     nn = n;
  91   elseif (isscalar(nn))
  92     m = ceil(log2(nn+1));
  93     n = 2.^m - 1;
  94     if ((n != floor(n)) || (n < 7) || (m != floor(m)) )
  95       error("bchpoly: n must be a integer greater than 3");
  96     endif
  97   else
  98     prim = bi2de(n);
  99     if (!isprimitive(prim))
 100       error ("bchpoly: prim must be a primitive polynomial of GF(2^m)");
 101     endif
 102     m = length(n) - 1;
 103     n = 2^m - 1;
 104   endif
 105
 106   if ((nargin > 1) && (!isscalar(k) || (floor(k) != k) || (k > n)))
 107     error ("bchpoly: message length must be less than codeword length");
 108   endif
 109
 110   for i=1:length(varargin)
 111     arg = varargin{i};
 112     if (ischar(arg))
 113       probe = 1;
 114       if (nargout > 1)
 115               error ("bchpoly: only one output argument allowed when probing valid codes");
 116       endif
 117     else
 118       if (prim != 0)
 119               error ("bchpoly: primitive polynomial already defined");
 120       endif
 121       prim = arg;
 122       if (!isscalar(prim))
 123               prim = bi2de(prim);
 124       endif
 125       if ((floor(prim) != prim) || (prim < 2^m) || (prim > 2^(m+1)) || ...
 126                 !isprimitive(prim))
 127               error ("bchpoly: prim must be a primitive polynomial of GF(2^m)");
 128       endif
 129     endif
 130   end
 131
 132   ## Am I using the right algo to calculate the correction capability?
 133   if (nargin < 2)
 134     if (nargout > 1)
 135       error ("bchpoly: only one output argument allowed when probing valid codes");
 136     endif
 137
 138     p = [];
 139     for ni=1:length(n)
 140       c = cosets(m(ni), prim);
 141       nc = length(c);
 142       fc = zeros(1,nc);
 143       f = [];
 144
 145       for t=1:floor(n(ni)/2)
 146               for i=1:nc
 147                 if (fc(i) != 1)
 148                   cl = log(c{i});
 149                   for j=2*(t-1)+1:2*t
 150                     if (find(cl == j))
 151                             f = [f, c{i}.x];
 152                             fc(i) = 1;
 153                             break;
 154                     endif
 155                   end
 156                 endif
 157               end
 158
 159               k = nn(ni) - length(f);
 160               if (k < 2)
 161                 break;
 162               endif
 163
 164               if (!isempty(p) && (k == p(size(p,1),2)))
 165                 p(size(p,1),:) = [nn(ni), k, t];
 166               else
 167                 p = [p; [nn(ni), k, t]];
 168               endif
 169       end
 170     end
 171   else
 172     c = cosets(m, prim);
 173     nc = length(c);
 174     fc = zeros(1,nc);
 175     f = [];
 176     fl = 0;
 177     f0 = [];
 178     f1 = [];
 179     t = 0;
 180     do
 181       t++;
 182       f0 = f1;
 183       for i=1:nc
 184               if (fc(i) != 1)
 185                 cl = log(c{i});
 186                 for j=2*(t-1)+1:2*t
 187                   if (find(cl == j))
 188                     f1 = [f1, c{i}.x];
 189                     fc(i) = 1;
 190                     ptmp = gf([c{i}(1), 1], m, prim);
 191                     for l=2:length(c{i})
 192                             ptmp = conv(ptmp, [c{i}(l), 1]);
 193                     end
 194                     f = [f; [ptmp.x, zeros(1,m-length(ptmp)+1)]];
 195                     fl = fl + length(ptmp);
 196                     break;
 197                   endif
 198                 end
 199               endif
 200       end
 201     until (length(f1) > nn - k)
 202     t--;
 203
 204     if (nn - length(f0) != k)
 205       error("bchpoly: can not find valid generator polynomial for parameters");
 206     endif
 207
 208     if (probe)
 209       p = [nn, k, t];
 210     else
 211
 212       ## Have to delete a line from the list of minimum polynomials
 213       ## since we've gone one past in calculating f1 above to be
 214       ## sure or the error correcting capability
 215       f = f(1:size(f,1)-1,:);
 216
 217       p = gf([f0(1), 1], m, prim);
 218       for i=2:length(f0)
 219               p = conv(p, [f0(i), 1]);
 220       end
 221       p = p.x;
 222
 223       if (nargout > 3)
 224               if (n > 64)
 225                 warning("bchpoly: can not create parity matrix\n");
 226                 par = [];
 227               else
 228                 par = cyclgen(n,p);
 229               endif
 230       endif
 231     endif
 232   endif
 233 endfunction