## [e45daf]: doc / comms.txi  Maximize  Restore  History

### 2184 lines (1728 with data), 75.3 kB

\input texinfo

@setfilename comms.info

@settitle Communications Toolbox for Octave

@titlepage
@title  Communications Toolbox for Octave
@subtitle March 2003
@author David Bateman
@author Laurent Mazet
@author Paul Kienzle
@page
@vskip 0pt plus 1filll

Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
are preserved on all copies.

Permission is granted to copy and distribute modified versions of this
manual under the conditions for verbatim copying, provided that the entire
notice identical to this one.

Permission is granted to copy and distribute translations of this manual
into another language, under the same conditions as for modified versions.
@end titlepage

@contents

@ifnottex
@node Top, Introduction
@top
@end ifnottex

* Introduction::
* Random Signals::
* Source Coding::
* Block Coding::
* Convolutional Coding::
* Modulations::
* Special Filters::
* Galois Fields::
* Function Reference::

@node Introduction, Random Signals, Top, Top
@chapter Introduction

This is the start of documentation for a Communications Toolbox for
Octave. As functions are written they should be documented here. In addition
many of the existing functions of Octave and Octave-Forge are important
in this Toolbox and their documentation should perhaps be repeated here.

This is preliminary documentation and you are invited to improve it and
submit patches.

@node Random Signals, Source Coding, Introduction, Top
@chapter Random Signals

The purpose of the functions described here is to create and add random
noise to a signal, to create random data and to analyze the eventually
errors in a received signal. The functions to perform these tasks can
be considered as either related to the creation or analysis of signals
and are treated separately below.

It should be noted that the examples below are based on the output of a
random number generator, and so the user can not expect to exactly recreate
the examples below.

* Signal Creation::
* Signal Analysis::

@node Signal Creation, Signal Analysis, , Random Signals
@section Signal Creation

The signal creation functions here fall into to two classes. Those that
treat discrete data and those that treat continuous data. The basic
function to create discrete data is @code{randint}, that creates a
random matrix of equi-probable integers in a desired range. For example

@example
octave:1> a = randint (3, 3, [-1, 1])
a =

0   1   0
-1  -1   1
0   1   1
@end example

creates a 3-by-3 matrix of random integers in the range -1 to 1. To allow
for repeated analysis with the same random data, the function @code{randint}
allows the seed-value of the random number generator to be set. For instance

@example
octave:1> a = randint (3, 3, [-1, 1], 1)
a =

0   1   1
0  -1   0
1  -1  -1
@end example

will always produce the same set of random data. The range of the integers
to produce can either be a two element vector or an integer. In the case
of a two element vector all elements within the defined range can be produced.
In the case of an integer range @var{M}, @code{randint} returns the equi-probable
integers in the range
@tex
$[0:2^m-1]$.
@end tex
@ifnottex
[0:2^@var{m}-1].
@end ifnottex

The function @code{randsrc} differs from @code{randint} in that it allows
a random set of symbols to be created with a given probability. The symbols
can be real, complex or even characters. However characters and scalars
can not be mixed. For example

@example
octave:1> a = randsrc (2, 2, "ab");
octave:2> b = randsrc (4, 4, [1, 1i, -1, -1i]);
@end example

are both legal, while

@example
octave:1> a = randsrc (2, 2, [1, "a"]);
@end example

is not legal. The alphabet from which the symbols are chosen can be either
a row vector or two row matrix. In the case of a row vector, all of the
elements of the alphabet are chosen with an equi-probability. In the case
of a two row matrix, the values in the second row define the probability
that each of the symbols are chosen. For example

@example
octave:1> a = randsrc (5, 5, [1, 1i, -1, -1i; 0.6 0.2 0.1 0.1])
a =

1 + 0i   0 + 1i   0 + 1i   0 + 1i   1 + 0i
1 + 0i   1 + 0i   0 + 1i   0 + 1i   1 + 0i
-0 - 1i   1 + 0i  -1 + 0i   1 + 0i   0 + 1i
1 + 0i   1 + 0i   1 + 0i   1 + 0i   1 + 0i
-1 + 0i  -1 + 0i   1 + 0i   1 + 0i   1 + 0i
@end example

defines that the symbol '1' has a 60% probability, the symbol '1i' has
a 20% probability and the remaining symbols have 10% probability each.
The sum of the probabilities must equal one. Like @code{randint},
@code{randsrc} accepts a fourth argument as the seed of the random
number generator allowing the same random set of data to be reproduced.

The function @code{randerr} allows a matrix of random bit errors to be
created, for binary encoded messages. By default, @code{randerr} creates
exactly one errors per row, flagged by a non-zero value in the returned
matrix. That is

@example
octave:1> a = randerr (5, 10)
a =

0   1   0   0   0   0   0   0   0   0
0   0   1   0   0   0   0   0   0   0
0   0   1   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   1
0   0   0   0   0   0   0   0   0   1
@end example

The number of errors per row can be specified as the third argument to
@code{randerr}. This argument can be either a scalar, a row vector or
a two row matrix. In the case of a scalar value, exactly this number of
errors will be created per row in the returned matrix. In the case of
a row vector, each element of the row vector gives a possible number of
equi-probable bit errors. The second row of a two row matrix defines
the probability of each number of errors occurring. For example

@example
octave:1> n = 15; k = 11; nsym = 100;
octave:2> msg = randint (nsym, k);       ## Binary vector of message
octave:3> code = encode (msg, n, k, "bch");
octave:4> berrs = randerr (nsym, n, [0, 1; 0.7, 0.3]);
octave:5> noisy = mod (code + berrs, 2)  ## Add errors to coded message
@end example

creates a vector @var{msg}, encodes it with a [15,11] BCH code, and then
add either none or one error per symbol with the chances of an error being
30%. As previously, @code{randerr} accepts a fourth argument as the seed of
the random number generator allowing the same random set of data to be
reproduced.

All of the above functions work on discrete random signals. The functions
@code{wgn} and @code{awgn} create and add white Gaussian noise to continuous
signals. The function @code{wgn} creates a matrix of white Gaussian noise
of a certain power. A typical call to @code{wgn} is then

@example
octave:1> nse = wgn (10, 10, 0);
@end example

Which creates a 10-by-10 matrix of noise with a root mean squared power
of 0dBW relative to an impedance of
@tex
$1\Omega$.
@end tex
@ifnottex
1 Ohm.
@end ifnottex

This effectively means that an equivalent result to the above can be
obtained with

@example
octave:1> nse = randn (10, 10);
@end example

The reference impedance and units of power to the function @code{wgn}
can however be modified, for example

@example
octave:1> nse_30dBm_50Ohm = wgn (10000, 1, 30, 50, "dBm");
octave:2> nse_0dBW_50Ohm = wgn (10000, 1, 0, 50, "dBW");
octave:3> nse_1W_50Ohm = wgn (10000, 1, 1, 50, "linear");
octave:4> [std(nse_30dBm_50Ohm), std(nse_0dBW_50Ohm), std(nse_1W_50Ohm)]
ans =

7.0805   7.1061   7.0730
@end example

All produce a 1W signal referenced to a
@tex
$50\Omega$.
@end tex
@ifnottex
50 Ohm
@end ifnottex
impedance. Matlab uses the misnomer "dB" for "dBW", and therefore "dB" is
an accepted type for @code{wgn} and will be treated as for "dBW".

In all cases, the returned matrix @var{v}, will be related to the input
power @var{p} and the impedance @var{Z} as

@tex
$$p = {\\sum_i \\sum_j v(i,j)^2 \\over Z} Watts$$
@end tex
@ifnottex
@var{p} = sum (@var{v}(:) .^ 2 ) / @var{imp} Watts
@end ifnottex

By default @code{wgn} produces real vectors of white noise. However, it can
produce both real and complex vectors like

@example
octave:1> rnse = wgn (10000, 1, 0, "dBm", "real");
octave:2> cnse = wgn (10000, 1, 0, "dBm", "complex");
octave:3> [std(rnse), std(real (cnse)), std(imag (cnse)), std(cnse)]
ans =

0.031615   0.022042   0.022241   0.031313
@end example

which shows that with a complex return value that the total power is the
same as a real vector, but that it is equally shared between the real and
imaginary parts. As previously, @code{wgn} accepts a fourth numerical argument
as the seed of the random number generator allowing the same random set of
data to be reproduced. That is

@example
octave:1> nse = wgn (10, 10, 0, 0);
@end example

will always produce the same set of data.

The final function to deal with the creation of random signals is
@code{awgn}, that adds noise at a certain level relative to a desired
signal. This function adds noise at a certain level to a desired
signal. An example call to @code{awgn} is

@example
octave:1> x = [0:0.1:2*pi];
octave:2> y = sin (x);
octave:3> noisy = awgn (y, 10, "measured")
@end example

@ifnotinfo
which produces a sine-wave with noise added as seen in Figure 1.

@center @image{awgn}

@center Figure 1: Sine-wave with 10dB signal-to-noise ratio
@end ifnotinfo

which adds noise with a 10dB signal-to-noise ratio to the measured power
in the desired signal. By default @code{awgn} assumes that the desired
signal is at 0dBW, and the noise is added relative to this assumed
power. This behavior can be modified by the third argument to @code{awgn}.
If the third argument is a numerical value, it is assumed to define the
power in the input signal, otherwise if the third argument is the string
"measured", as above, the power in the signal is measured prior to the

The final argument to @code{awgn} defines the definition of the power and
signal-to-noise ratio in a similar manner to @code{wgn}. This final
argument can be either "dB" or "linear". In the first case the numerical
value of the input power is assumed to be in dBW and the signal-to-noise
ratio in dB.  In the second case, the power is assumed to be in Watts
and the signal-to-noise ratio is expressed as a ratio.

The return value of @code{awgn} will be in the same form as the input
signal. In addition if the input signal is real, the additive noise will
be real. Otherwise the additive noise will also be complex and the noise
will be equally split between the real and imaginary parts.

As previously the seed to the random number generator can be specified
as the last argument to @code{awgn} to allow repetition of the same
scenario. That is

@example
octave:1> x = [0:0.1:2*pi];
octave:2> y = sin (x);
octave:3> noisy = awgn (y, 10, "dB", 0, "measured")
@end example

which uses the seed-value of 0 for the random number generator.

@node Signal Analysis, , Signal Creation, Random Signals
@section Signal Analysis

It is important to be able to evaluate the performance of a
communications system in terms of its bit-error and symbol-error
rates. Two functions @code{biterr} and @code{symerr} exist within this
package to calculate these values, both taking as arguments the
expected and the actually received data. The data takes the form
of matrices or vectors, with each element representing a single
symbol. They are compared in the following manner

@table @asis
@item Both matrices
In this case both matrices must be the same size and then by default the
return values are the overall number of errors and the overall error rate.
@item One column vector
In this case the column vector is used for comparison column-wise
with the matrix. The return values are row vectors containing the number
of errors and the error rate for each column-wise comparison. The number
of rows in the matrix must be the same as the length of the column vector.
@item One row vector
In this case the row vector is used for comparison row-wise
with the matrix. The return values are column vectors containing the number
of errors and the error rate for each row-wise comparison. The number
of columns in the matrix must be the same as the length of the row vector.
@end table

For the bit-error comparison, the size of the symbol is assumed to be the
minimum number of bits needed to represent the largest element in the
two matrices supplied. However, the number of bits per symbol can (and
in the case of random data should) be specified. As an example of the
use of @code{biterr} and @code{symerr}, consider the example

@example
octave:1> m = 8;
octave:2> msg = randint (10, 10, 2^m);
octave:3> noisy = mod (msg + diag (1:10), 2^m);
octave:4> [berr, brate] = biterr (msg, noisy, m)
berr = 32
brate = 0.040000
octave:5> [serr, srate] = symerr (msg, noisy)
serr = 10
srate = 0.10000
@end example

which creates a 10-by-10 matrix adds 10 symbols errors to the data and then
finds the bit and symbol error-rates.

Two other means of displaying the integrity of a signal are the
eye-diagram and the scatterplot. Although the functions
@code{eyediagram} and @code{scatterplot} have different appearance, the
information presented is similar and so are their inputs. The difference
between @code{eyediagram} and @code{scatterplot} is that @code{eyediagram}
segments the data into time intervals and plots the in-phase and
quadrature components of the signal against this time interval. While
@code{scatterplot} uses a parametric plot of quadrature versus in-phase
components.

Both functions can accept real or complex signals in the following
forms.

@table @asis
@item A real vector
In this case the signal is assumed to be real and represented by the vector
@var{x}.
@item A complex vector
In this case the in-phase and quadrature components of the signal are
assumed to be the real and imaginary parts of the signal.
@item A matrix with two columns
In this case the first column represents the in-phase and the second the
quadrature components of a complex signal.
@end table

An example of the use of the function @code{eyediagram} is

@example
octave:1> n = 50;
octave:2> ovsp = 50;
octave:3> x = 1:n;
octave:4> xi = 1:1/ovsp:n-0.1;
octave:5> y = randsrc (1, n, [1 + 1i, 1 - 1i, -1 - 1i, -1 + 1i]);
octave:6> yi = interp1 (x, y, xi);
octave:7> noisy = awgn (yi, 15, "measured");
octave:8> eyediagram (noisy, ovsp);
@end example

@ifnotinfo
which produces a eye-diagram of a noisy signal as seen in Figure 2. Similarly
an example of the use of the function @code{scatterplot} is

@center @image{eyediagram}

@center Figure 2: Eye-diagram of a QPSK like signal with 15dB signal-to-noise ratio
@end ifnotinfo
@ifinfo
Similarly an example of the use of the function @code{scatterplot} is
@end ifinfo

@example
octave:1> n = 200;
octave:2> ovsp = 5;
octave:3> x = 1:n;
octave:4> xi = 1:1/ovsp:n-0.1;
octave:5> y = randsrc (1, n, [1 + 1i, 1 - 1i, -1 - 1i, -1 + 1i]);
octave:6> yi = interp1 (x, y, xi);
octave:7> noisy = awgn (yi, 15, "measured");
octave:8> f = scatterplot (noisy, 1, 0, "b");
octave:9> hold on;
octave:10> scatterplot (noisy, ovsp, 0, "r+", f);
@end example

@ifnotinfo
which produces a scatterplot of a noisy signal as seen in Figure 3.

@center @image{scatterplot}

@center Figure 3: Scatterplot of a QPSK like signal with 15dB signal-to-noise ratio
@end ifnotinfo

@node Source Coding, Block Coding, Random Signals, Top
@chapter Source Coding

* Quantization::
* PCM Coding::
* Arithmetic Coding::
* Dynamic Range Compression::

@node Quantization, PCM Coding, , Source Coding
@section Quantization

An important aspect of converting an analog signal to the digital domain
is quantization. This is the process of mapping a continuous signal to a
set of defined values. Octave contains two functions to perform quantization,
@code{lloyds} creates an optimal mapping of the continuous signal to a fixed
number of levels and @code{quantiz} performs the actual quantization.

The set of quantization points to use is represented by a partitioning
table (@var{table}) of the data and the signal levels (@var{codes} to
which they are mapped. The partitioning @var{table} is monotonically
increasing and if x falls within the range given by two points of this
table then it is mapped to the corresponding code as seen in Table 1.

@center Table 1: Table quantization partitioning and coding

@multitable @columnfractions 0.1 0.4 0.4 0.1
@item @tab x < table(1)               @tab codes(1)   @tab
@item @tab table(1) <= x < table(2)   @tab codes(2)   @tab
@item @tab @dots{}                    @tab @dots{}    @tab
@item @tab table(i-1) <= x < table(i) @tab codes(i)   @tab
@item @tab @dots{}                    @tab @dots{}    @tab
@item @tab table(n-1) <= x < table(n) @tab codes(n)   @tab
@item @tab table(n-1) <= x            @tab codes(n+1) @tab
@end multitable

These partition and coding tables can either be created by the user of
using the function @code{lloyds}. For instance the use of a linear
mapping can be seen in the following example.

@example
octave:1> m = 8;
octave:2> n = 1024;
octave:3> table = 2*[0:m-1]/m - 1 + 1/m;
octave:4> codes = 2*[0:m]/m - 1;
octave:5> x = 4*pi*[0:(n-1)]/(n-1);
octave:6> y = cos (x);
octave:7> [i, z] = quantiz (y, table, codes);
@end example

If a training signal is known that well represents the expected signals,
the quantization levels can be optimized using the @code{lloyds} function.
For example the above example can be continued

@example
octave:8> [table2, codes2] = lloyds (y, table, codes);
octave:9> [i, z2] = quantiz (y, table2, codes2);
@end example

Which the mapping suggested by the function @code{lloyds}. It should be
noted that the mapping given by @code{lloyds} is highly dependent on the
training signal used. So if this signal does not represent a realistic
signal to be quantized, then the partitioning suggested by @code{lloyds}
will be sub-optimal.

@node PCM Coding, Arithmetic Coding, Quantization, Source Coding
@section PCM Coding

The DPCM function @code{dpcmenco}, @code{dpcmdeco} and @code{dpcmopt}
implement a form of predictive quantization, where the predictability
of the signal is used to further compress it. These functions are
not yet implemented.

@node Arithmetic Coding, Dynamic Range Compression, PCM Coding, Source Coding
@section Arithmetic Coding

The arithmetic coding functions @code{arithenco} and @code{arithdeco} are
not yet implemented.

@node Dynamic Range Compression, , Arithmetic Coding, Source Coding
@section Dynamic Range Compression

The final source coding function is @code{compand} which is used to
compress and expand the dynamic range of a signal. For instance
consider a logarithm quantized by a linear partitioning. Such a
partitioning is very poor for this large dynamic range. @code{compand}
can then be used to compress the signal prior to quantization, with
the signal being expanded afterwards. For example

@example
octave:1> mu = 1.95;
octave:2> x = [0.01:0.01:2];
octave:3> y = log (x);
octave:4> V = max (abs (y));
octave:5> [i, z, d] = quantiz (y, [-4.875:0.25:0.875], [-5:0.25:1]);
octave:6> c = compand (y, minmu, V, "mu/compressor");
octave:7> [i2, c2] = quantiz (c, [-4.875:0.25:0.875], [-5:0.25:1]);
octave:8> z2 = compand (c2, minmu, max (abs (c2)), "mu/expander");
octave:9> d2 = sumsq (y - z2) / length (y);
octave:10> [d, d2]
ans =

0.0053885   0.0029935
@end example

which demonstrates that the use of @code{compand} can significantly
reduce the distortion due to the quantization of signals with a large
dynamic range.

@node Block Coding, Convolutional Coding, Source Coding, Top
@chapter Block Coding

The error-correcting codes available in this toolbox are discussed here.
These codes work with blocks of data, with no relation between one block
and the next.  These codes create codewords based on the messages to
transmit that contain redundant information that allow the recovery of
the original message in the presence of errors.

* Data Formats::
* Binary Block Codes::
* BCH Codes::
* Reed-Solomon Codes::

@node Data Formats, Binary Block Codes, , Block Coding
@section Data Formats

All of the codes described in this section are binary and share
similar data formats. The exception is the Reed-Solomon coder which
has a significantly longer codeword length in general and therefore
using a different manner to efficiently pass data. The user should
reference to the section about the Reed-Solomon codes for the data
format for use with Reed-Solomon codes.

In general @var{k} bits of data are considered to represent a single
message symbol. These @var{k} bits are coded into @var{n} bits of
data representing the codeword. The data can therefore be grouped
in one of three manners, to emphasis this grouping into bits, messages
and codewords

@table @asis
@item A binary vector
Each element of the vector is either one or zero. If the data represents
an uncoded message the vector length should be an integer number of @var{k}
in length.
@item A binary matrix
In this case the data is ones and zeros grouped into rows, with each
representing a single message or codeword. The number of columns in the
matrix should be equal to @var{k} in the case of a uncoded message or
@var{n} in the case of a coded message.
@item A non-binary vector
In this case each element of the vector represents a message or codeword
in an integer format. The bits of the message or codeword are represented
by the bits of the vector elements with the least-significant bit
representing the first element in the message or codeword.
@end table

An example demonstrating the relationship between the three data
formats can be seen below.

@example
octave:1> k = 4;
octave:2> bin_vec = randint (k*10, 1);         # Binary vector format
octave:3> bin_mat = reshape (bin_vec, k, 10)'; # Binary matrix format
octave:4> dec_vec = bi2de (bin_mat);           # Decimal vector format
@end example

The functions within this toolbox will return data in the same format
to which it is given. It should be noted that internally the binary
matrix format is used, and thus if the message or codeword length is
large it is preferable to use the binary format to avoid internal
rounding errors.

@node Binary Block Codes, BCH Codes, Data Formats, Block Coding
@section Binary Block Codes

All of the codes presented here can be characterized by their

@table @asis
@item Generator Matrix
A @var{k}-by-@var{n} matrix @var{G} to generate the codewords @var{C} from
the messages @var{T} by the matrix multiplication
@tex
${\bf C} = {\bf T} {\bf G}$.
@end tex
@ifnottex
@var{C} = @var{T} * @var{G}.
@end ifnottex
@item Parity Check Matrix
A '@var{n}-@var{k}'-by-@var{n} matrix @var{H} to check the parity of the
@tex
${\bf H} {\bf R} = {\bf S} \ne 0$,
@end tex
@ifnottex
@var{H} * @var{R} = @var{S} != 0,
@end ifnottex
then an error has been detected. @var{S} can be used with the syndrome
table to correct this error
@item Syndrome Table
A 2^@var{k}-by-@var{n} matrix @var{ST} with the relationship of the error
vectors to the non-zero parities of the received symbols. That is, if
the received symbol is represented as
@tex
${\bf R} = ( {\bf T} + {\bf E} )\ mod\ 2$,
@end tex
@ifnottex
@var{R} = mod (@var{T} + @var{E}, 2),
@end ifnottex
then the error vector @var{E} is
@tex
${\bf ST}({\bf S})$.
@end tex
@ifnottex
@var{ST}(@var{S}).
@end ifnottex
@end table

It is assumed for most of the functions in this toolbox that the generator
matrix will be in a 'standard' form. That is the generator matrix can be
represented by

@tex
$${\bf G} = \left[\matrix{g_{11} & g_{12} & \ldots & g_{1k} & 1 & 0 & \ldots & 0 \cr g_{21} & g_{22} & & g_{2k} & 0 & 1 & & 0 \cr \vdots & & & \vdots & \vdots& & & \vdots \cr g_{k1} & g_{k2} & \ldots & g_{kk} & 0 & 0 & \ldots & 1}\right]$$
@end tex
@ifnottex
@example
@group
g(1,1)  g(1,2)  ...  g(1,k)  1  0  ...  0
g(2,1)  g(2,2)       g(2,k)  0  1  ...  0
.                    .     .          .
.                    .     .          .
.                    .     .          .
g(k,1)  g(k,2)  ...  g(k,k)  0  0  ...  1
@end group
@end example
@end ifnottex

or

@tex
$${\bf G} = \left[\matrix{1 & 0 & \ldots & 0 & g_{11} & g_{12} & \ldots & g_{1k} \cr 0 & 1 & & 0 & g_{21} & g_{22} & & g_{2k} \cr \vdots & & & \vdots & \vdots & & & \vdots \cr 0 & 0 & \ldots & 1 & g_{k1} & g_{k2} & \ldots & g_{kk}}\right]$$
@end tex
@ifnottex
@example
@group
1  0  ...  0  g(1,1)  g(1,2)  ...  g(1,k)
0  1  ...  0  g(2,1)  g(2,2)       g(2,k)
.          .     .                   .
.          .     .                   .
.          .     .                   .
0  0  ...  1  g(k,1)  g(k,2)  ...  g(k,k)
@end group
@end example
@end ifnottex

and similarly the parity check matrix can be represented by a combination
of an identity matrix and a square matrix.

Some of the codes can also have their representation in terms of a
generator polynomial that can be used to create the generator and parity
check matrices. In the case of BCH codes, this generator polynomial is
used directly in the encoding and decoding without ever explicitly forming
the generator or parity check matrix.

The user can create their own generator and parity check matrices, or
they can rely on the functions @code{hammgen}, @code{cyclgen} and
@code{cyclpoly}. The function @code{hammgen} creates parity check and
generator matrices for Hamming codes, while @code{cyclpoly} and
@code{cyclgen} create generator polynomials and matrices for generic
cyclic codes. An example of their use is

@example
octave:1> m = 3;
octave:2> n = 2^m - 1;
octave:2> k = 4;
octave:3> [par, gen] = hammgen (m);
octave:4> [par2, gen2] = cyclgen (n, cyclpoly (n, k));
@end example

which create identical parity check and generator matrices for the
[7,4] Hamming code.

The syndrome table of the codes can be created with the function
@code{syndtable}, in the following manner

@example
octave:1> [par, gen] = hammgen (3);
octave:2> st = syndtable (par);
@end example

There exists two auxiliary functions @code{gen2par} and @code{gfweight},
that convert between generator and parity check matrices and calculate
the Hamming distance of the codes. For instance

@example
octave:1> par = hammgen (3);
octave:2> gen = gen2par (par);
octave:3> gfweight (gen)
ans = 3
@end example

It should be noted that for large values of @var{n}, the generator,
parity check and syndrome table matrices are very large. There is
therefore an internal limitation on the size of the block codes that
can be created that limits the codeword length @var{n} to less than 64.
Which is still excessively large for the syndrome table, so use caution
with these codes. These limitations do not apply to the Reed-Solomon
or BCH codes.

The top-level encode and decode functions are @code{encode} and
@code{decode}, which can be used with all codes, except the Reed-Solomon
code. The basic call to both of these functions passes the message
to code/decode, the codeword length, the message length and the type
of coding to use. There are four basic types that are available with
these functions

@table @asis
@item "linear"
Generic linear block codes
@item "cyclic"
Cyclic linear block codes
@item "hamming"
Hamming codes
@item "bch"
Bose Chaudhuri Hocquenghem (BCH) block codes
@end table

It is not possible to distinguish between a binary vector and a decimal
vector coding of the messages that just happens to only have ones and
zeros. Therefore the functions @code{encode} and @code{decode} must be
told the format of the messages in the following manner.

@example
octave:1> m = 3;
octave:2> n = 7;
octave:3> k = 4;
octave:4> msg_bin = randint (10, k);
octave:5> cbin = encode (msg_bin, n, k, "hamming/binary");
octave:5> cdec = encode (bi2de (msg), n, k, "hamming/decimal");
@end example

which codes a binary matrix and a non-binary vector representation of a
message, returning the coded message in the same format. The functions
@code{encode} and @code{decode} by default accept binary coded
messages. Therefore "hamming" is equivalent to "hamming/binary".

Except for the BCH codes, the function @code{encode} and @code{decode}
internally create the generator, parity check and syndrome table
matrices.  Therefore if repeated calls to @code{encode} and @code{decode}
are made it will often be faster to create these matrices externally,
and pass them as an argument. For example

@example
n = 15;
k = 11;
[par, gen] = hammgen (4);
code1 = code2 = zeros (100, 15)
for i = 1:100
msg = get_msg (i);
code1(i,:) = encode (msg, n, k, "linear", gen);  # This is faster
code2(i,:) = encode (msg, n, k, "hamming");      # than this !!!
endfor
@end example

In the case of the BCH codes the low-level functions described in the
next section are used directly by the @code{encode} and @code{decode}
functions.

@node BCH Codes, Reed-Solomon Codes, Binary Block Codes, Block Coding
@section BCH Codes

The BCH coder used here is based on code written by Robert Morelos-Zaragoza
(r.morelos-zaragoza@@ieee.org). This code was originally written in C
and has been converted for use as an Octave oct-file.

@iftex
Called without arguments, @code{bchpoly} returns a table of valid BCH
error correcting codes and their error-correction capability
as seen in Table 1.

@center Table 2: Table of valid BCH codes with codeword length less than 511.

@multitable @columnfractions .083 .083 .083 .083 .083 .083 .083 .083 .083 .083 .083 .083
@item  N  @tab  K  @tab  T  @tab  N  @tab  K  @tab  T  @tab  N  @tab  K  @tab  T  @tab  N  @tab  K  @tab  T
@item   7 @tab   4 @tab   1 @tab 127 @tab  36 @tab  15 @tab 255 @tab  45 @tab  43 @tab 511 @tab 268 @tab  29
@item  15 @tab  11 @tab   1 @tab 127 @tab  29 @tab  21 @tab 255 @tab  37 @tab  45 @tab 511 @tab 259 @tab  30
@item  15 @tab   7 @tab   2 @tab 127 @tab  22 @tab  23 @tab 255 @tab  29 @tab  47 @tab 511 @tab 250 @tab  31
@item  15 @tab   5 @tab   3 @tab 127 @tab  15 @tab  27 @tab 255 @tab  21 @tab  55 @tab 511 @tab 241 @tab  36
@item  31 @tab  26 @tab   1 @tab 127 @tab   8 @tab  31 @tab 255 @tab  13 @tab  59 @tab 511 @tab 238 @tab  37
@item  31 @tab  21 @tab   2 @tab 255 @tab 247 @tab   1 @tab 255 @tab   9 @tab  63 @tab 511 @tab 229 @tab  38
@item  31 @tab  16 @tab   3 @tab 255 @tab 239 @tab   2 @tab 511 @tab 502 @tab   1 @tab 511 @tab 220 @tab  39
@item  31 @tab  11 @tab   5 @tab 255 @tab 231 @tab   3 @tab 511 @tab 493 @tab   2 @tab 511 @tab 211 @tab  41
@item  31 @tab   6 @tab   7 @tab 255 @tab 223 @tab   4 @tab 511 @tab 484 @tab   3 @tab 511 @tab 202 @tab  42
@item  63 @tab  57 @tab   1 @tab 255 @tab 215 @tab   5 @tab 511 @tab 475 @tab   4 @tab 511 @tab 193 @tab  43
@item  63 @tab  51 @tab   2 @tab 255 @tab 207 @tab   6 @tab 511 @tab 466 @tab   5 @tab 511 @tab 184 @tab  45
@item  63 @tab  45 @tab   3 @tab 255 @tab 199 @tab   7 @tab 511 @tab 457 @tab   6 @tab 511 @tab 175 @tab  46
@item  63 @tab  39 @tab   4 @tab 255 @tab 191 @tab   8 @tab 511 @tab 448 @tab   7 @tab 511 @tab 166 @tab  47
@item  63 @tab  36 @tab   5 @tab 255 @tab 187 @tab   9 @tab 511 @tab 439 @tab   8 @tab 511 @tab 157 @tab  51
@item  63 @tab  30 @tab   6 @tab 255 @tab 179 @tab  10 @tab 511 @tab 430 @tab   9 @tab 511 @tab 148 @tab  53
@item  63 @tab  24 @tab   7 @tab 255 @tab 171 @tab  11 @tab 511 @tab 421 @tab  10 @tab 511 @tab 139 @tab  54
@item  63 @tab  18 @tab  10 @tab 255 @tab 163 @tab  12 @tab 511 @tab 412 @tab  11 @tab 511 @tab 130 @tab  55
@item  63 @tab  16 @tab  11 @tab 255 @tab 155 @tab  13 @tab 511 @tab 403 @tab  12 @tab 511 @tab 121 @tab  58
@item  63 @tab  10 @tab  13 @tab 255 @tab 147 @tab  14 @tab 511 @tab 394 @tab  13 @tab 511 @tab 112 @tab  59
@item  63 @tab   7 @tab  15 @tab 255 @tab 139 @tab  15 @tab 511 @tab 385 @tab  14 @tab 511 @tab 103 @tab  61
@item 127 @tab 120 @tab   1 @tab 255 @tab 131 @tab  18 @tab 511 @tab 376 @tab  15 @tab 511 @tab  94 @tab  62
@item 127 @tab 113 @tab   2 @tab 255 @tab 123 @tab  19 @tab 511 @tab 367 @tab  17 @tab 511 @tab  85 @tab  63
@item 127 @tab 106 @tab   3 @tab 255 @tab 115 @tab  21 @tab 511 @tab 358 @tab  18 @tab 511 @tab  76 @tab  85
@item 127 @tab  99 @tab   4 @tab 255 @tab 107 @tab  22 @tab 511 @tab 349 @tab  19 @tab 511 @tab  67 @tab  87
@item 127 @tab  92 @tab   5 @tab 255 @tab  99 @tab  23 @tab 511 @tab 340 @tab  20 @tab 511 @tab  58 @tab  91
@item 127 @tab  85 @tab   6 @tab 255 @tab  91 @tab  25 @tab 511 @tab 331 @tab  21 @tab 511 @tab  49 @tab  93
@item 127 @tab  78 @tab   7 @tab 255 @tab  87 @tab  26 @tab 511 @tab 322 @tab  22 @tab 511 @tab  40 @tab  95
@item 127 @tab  71 @tab   9 @tab 255 @tab  79 @tab  27 @tab 511 @tab 313 @tab  23 @tab 511 @tab  31 @tab 109
@item 127 @tab  64 @tab  10 @tab 255 @tab  71 @tab  29 @tab 511 @tab 304 @tab  25 @tab 511 @tab  28 @tab 111
@item 127 @tab  57 @tab  11 @tab 255 @tab  63 @tab  30 @tab 511 @tab 295 @tab  26 @tab 511 @tab  19 @tab 119
@item 127 @tab  50 @tab  13 @tab 255 @tab  55 @tab  31 @tab 511 @tab 286 @tab  27 @tab 511 @tab  10 @tab 127
@item 127 @tab  43 @tab  14 @tab 255 @tab  47 @tab  42 @tab 511 @tab 277 @tab  28 @tab     @tab     @tab
@end multitable

@end iftex
@ifnottex
Called without arguments, @code{bchpoly} returns a table of valid BCH
error correcting codes and their error-correction capability.
@end ifnottex
The first returned column of @code{bchpoly} is the codeword length,
the second the message length and the third the error correction capability
of the code. Called with one argument, @code{bchpoly} returns similar
output, but only for the specified codeword length. In this manner codes
with codeword length greater than 511 can be found.

In general the codeword length is of the form @code{2^@var{m} - 1}, where
@var{m} is an integer. However if [@var{n},@var{k}] is a valid BCH
code, then it is also possible to use a shortened BCH form of the form
@code{[@var{n}-@var{x},@var{k}-@var{x}]}.

With two or more arguments, @code{bchpoly} is used to find the generator
polynomial of a valid BCH code. For instance

@example
octave:1> bchpoly (15, 7)
ans =

1   0   0   0   1   0   1   1   1

octave:2> bchpoly (14, 6)
ans =

1   0   0   0   1   0   1   1   1
@end example

show that the generator polynomial of a [15,7] BCH code with the default
primitive polynomial is

@tex
$$1 + x^4 + x^6 + x^7 + x^8$$
@end tex
@ifnottex
1 + @var{x} ^ 4 + @var{x} ^ 6 + @var{x} ^ 7 + @var{x} ^ 8
@end ifnottex

Using a different primitive polynomial to define the Galois Field over
which the BCH code is defined results in a different generator polynomial
as can be seen in the example.

@example
octave:1> bchpoly ([1 1 0 0 1], 7)
ans =

1   0   0   0   1   0   1   1   1

octave:2> bchpoly ([1 0 0 1 1], 7)
ans =

1   1   1   0   1   0   0   0   1
@end example

It is recommend not to convert the generator polynomials created by
@code{bchpoly} into generator and parity check matrices with the
BCH codes, as the underlying BCH software is faster than the generic
coding software and can treat significantly longer codes.

As well as using the @code{encode} and @code{decode} functions previously
discussed, the user can directly use the low-level BCH functions
@code{bchenco} and @code{bchdeco}. In this case the messages must be
in the format of a binary matrix with @var{k} columns

@example
octave:1> n = 31;
octave:2> pgs = bchpoly (n);
octave:3> pg = pgs(floor (rand () * (rows (pgs) + 1)),:);  # Pick a poly
octave:4> k = pg(2);
octave:5> t = pg(3);
octave:6> msg = randint (10, k);
octave:7> code = bchenco (msg, n, k);
octave:8> noisy = code + [ones(10, 1), zeros(10, n-1)];
octave:9> dec = bchdeco (code, k, t);
@end example

@node Reed-Solomon Codes, , BCH Codes, Block Coding
@section Reed-Solomon Codes

* Representation of Reed-Solomon Messages::
* Creating and Decoding Messages::
* Shortened Reed-Solomon Codes::

@node Representation of Reed-Solomon Messages, Creating and Decoding Messages, , Reed-Solomon Codes
@subsection Representation of Reed-Solomon Messages

The Reed-Solomon coder used in this package is based on code written by
Phil Karn (http://www.ka9q.net/code/fec). This code was originally written
in C and has been converted for use as an Octave oct-file.

Reed-Solomon codes are based on Galois Fields of even characteristics
GF(2^M). Many of the properties of Galois Fields are therefore important
when considering Reed-Solomon coders.

The representation of the symbols of the Reed-Solomon code differs from
the other block codes, in that the other block codes use a binary
representation, while the Reed-Solomon code represents each m-bit symbol
by an integer. The elements of the message and codeword must be elements
of the Galois Field corresponding to the Reed-Solomon code. Thus to
code a message with a [7,5] Reed-Solomon code an example is

@example
octave:1> m = 3;
octave:2> n = 7;
octave:3> k = 5;
octave:4> msg = gf (floor (2^m * rand (2, k)), m)
msg =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

5   0   6   3   2
4   1   3   1   2

octave:5> code = rsenc (msg, n, k)
code =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

5   0   6   3   2   3   5
4   1   3   1   2   6   3
@end example

The variable @var{n} is the codeword length of the Reed-Solomon coder,
while @var{k} is the message length. It should be noted that @var{k}
should be less than @var{n} and that @code{@var{n} - @var{k}} should be
even. The error correcting capability of the Reed-Solomon code is then
@code{(@var{n} - @var{k})/2} symbols. @var{m} is the number of bits per
symbol, and is related to @var{n} by @code{@var{n} = 2^@var{m} - 1}.
For a valid Reed-Solomon coder, @var{m} should be between 3 and 16.

@node Creating and Decoding Messages, Shortened Reed-Solomon Codes, Representation of Reed-Solomon Messages, Reed-Solomon Codes
@subsection Creating and Decoding Messages

The Reed-Solomon encoding function requires at least three arguments. The
first @var{msg} is the message in encodes, the second is @var{n} the codeword
length and @var{k} is the message length. Therefore @var{msg} must have
@var{k} columns and the output will have @var{n} columns of symbols.

The message itself is many up of elements of a Galois Field
GF(2^M). Normally, The order of the Galois Field (M), is related to the
codeword length by @code{@var{n} = 2^@var{m} - 1}. Another important
parameter when determining the behavior of the Reed-Solomon coder is
the primitive polynomial of the Galois Field (see @code{gf}). Thus
the messages

@example
octave:1> msg0 = gf ([0, 1, 2, 3], 3);
octave:2> msg1 = gf ([0, 1, 2, 3], 3, 13);
@end example

will not result in the same Reed-Solomon coding. Finally, the parity of
the Reed-Solomon code are generated with the use of a generator
polynomial. The parity symbols are then generated by treating the message
to encode as a polynomial and finding the remainder of the division of
this polynomial by the generator polynomial. Therefore the generator
polynomial must have as many roots as @code{@var{n} - @var{k}}. Whether
the parity symbols are placed before or afterwards the message will then
determine which end of the message is the most-significant term of the
polynomial representing the message. The parity symbols are therefore
different in these two cases. The position of the parity symbols can be
chosen by specifying "beginning" or "end" to @code{rsenc} and @code{rsdec}.
By default the parity symbols are placed after the message.

Valid generator polynomials can be constructed with the @code{rsgenpoly}
function. The roots of the generator polynomial are then defined by

@tex
$$g = (x - A^{bs}) (x - A^{(b+1)s}) \cdots (x - A ^{(b+2t-1)s}).$$
@end tex
@ifnottex

@example
@var{g} = (@var{x} - A^(@var{b}*@var{s})) * (@var{x} - A^((@var{b}+1)*@var{s})) * @dots{} * (@var{x} - A^((@var{b}+2*@var{t}-1)*@var{s})).
@end example
@end ifnottex

where @var{t} is @code{(@var{n} - @var{k})/2}, A is the primitive element
of the Galois Field, @var{b} is the first consecutive root, and @var{s}
is the step between roots. Generator polynomial of this form are constructed
by @code{rsgenpoly} and can be passed to both @code{rsenc} and @code{rsdec}.
It is also possible to pass the @var{b} and @var{s} values directly to
@code{rsenc} and @code{rsdec}. In the case of @code{rsdec} passing @var{b}
and @var{s} can make the decoding faster.

Consider the example below.

@example
octave:1> m = 8;
octave:2> n = 2^m - 1;
octave:3> k = 223;
octave:4> prim = 391;
octave:5> b = 112;
octave:6> s = 11;
octave:7> gg = rsgenpoly (n, k, prim, b, s);
octave:8> msg = gf (floor (2^m * rand (17, k)), m, prim);
octave:9> code = rsenc (msg, n, k, gg);
octave:10> noisy = code + [toeplitz([ones(1,17)], zeros(1,17)), zeros(17,238)];
octave:11> [dec, nerr] = rsdec (msg, n, k, b, s);
octave:12> nerr'
ans =

1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   -1

octave:13> any (msg' != dec')
ans =

0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1
@end example

This is an interesting example in that it demonstrates many of the
additional arguments of the Reed-Solomon functions. In particular
this example approximates the CCSDS standard Reed-Solomon coder,
lacking only the dual-basis lookup tables used in this standard.
The CCSDS uses non-default values to all of the basic functions
involved in the Reed-Solomon encoding, since it has a non-default
primitive polynomial, generator polynomial, etc.

The example creates 17 message blocks and adds between 1 and 17 error
symbols to these block. As can be seen @var{nerr} gives the number of
errors corrected. In the case of 17 introduced errors @var{nerr}
equals -1, indicating a decoding failure. This is normal as the
correction ability of this code is up to 16 error symbols. Comparing
the input message and the decoding it can be seen that as expected,
only the case of 17 errors has not been correctly decoded.

@node Shortened Reed-Solomon Codes, , Creating and Decoding Messages, Reed-Solomon Codes
@subsection Shortened Reed-Solomon Codes

In general the codeword length of the Reed-Solomon coder is chosen so
that it is related directly to the order of the Galois Field by the
formula @code{@var{n} = 2^@var{m} - 1}. Although, the underlying
Reed-Solomon coding must operate over valid codeword length, there
are sometimes reasons to assume that the codeword length will be shorter.
In this case the message is padded with zeros before coding, and the
zeros are stripped from the returned block. For example consider the
shortened [6,4] Reed-Solomon below

@example
octave:1> m = 3;
octave:2> n = 6;
octave:3> k = 4;
octave:4> msg = gf (floor (2^m * rand (2, k)), m)
msg =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

7   0   2   5
1   5   7   1

octave:5> code = rsenc (msg, n, k)
code =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

7   0   2   5   2   3
1   5   7   1   0   2
@end example

@node Convolutional Coding, Modulations, Block Coding, Top
@chapter Convolutional Coding

To be written.

@node Modulations, Special Filters, Convolutional Coding, Top
@chapter Modulations

To be written.

Currently have functions amodce, ademodce, apkconst, demodmap, modmap,

@node Special Filters, Galois Fields, Modulations, Top
@chapter Special Filters

To be written.

@node Galois Fields, Function Reference, Special Filters, Top
@chapter Galois Fields

* Galois Field Basics::
* Manipulating Galois Fields::

@node Galois Field Basics, Manipulating Galois Fields, , Galois Fields
@section Galois Field Basics

A Galois Field is a finite algebraic field. This package implements a
Galois Field type in Octave having 2^M members where M is an integer
between 1 and 16. Such fields are denoted as GF(2^M) and are used in
error correcting codes in communications systems. Galois Fields having
odd numbers of elements are not implemented.

The @emph{primitive element} of a Galois Field has the property that all
elements of the Galois Field can be represented as a power of this element.
The @emph{primitive polynomial} is the minimum polynomial of some primitive
element in GF(2^M) and is irreducible and of order M. This means that the
primitive element is a root of the primitive polynomial.

The elements of the Galois Field GF(2^M) are represented as the values
0 to 2^M -1 by Octave. The first two elements represent the zero and unity
values of the Galois Field and are unique in all fields. The element
represented by 2 is the primitive element of the field and all elements can
be represented as combinations of the primitive element and unity as follows

@multitable @columnfractions .33 .33 .33
@item Integer    @tab Binary     @tab Element of GF(2^M)
@item  0         @tab  000       @tab @code{0}
@item  1         @tab  001       @tab @code{1}
@item  2         @tab  010       @tab @code{A}
@item  3         @tab  011       @tab @code{A + 1}
@item  4         @tab  100       @tab @code{A^2}
@item  5         @tab  101       @tab @code{A^2 + 1}
@item  6         @tab  110       @tab @code{A^2 + A}
@item  7         @tab  111       @tab @code{A^2 + A + 1}
@end multitable

It should be noted that there is often more than a single primitive
polynomial of GF(2^M). Each Galois Field over a different primitive
polynomial represents a different realization of the Field. The
representations above however rest valid.

This code was written as a challenge by Paul Kienzle (octave forge) to
convert a Reed-Solomon coder I had in Octave to be compatible with
Matlab communications toolbox R13. This forced the need to have a complete
library of functions over the even Galois Fields. Although this code
was written to be compatible with the equivalent Matlab code, I did not
have access to a version of Matlab with R13 installed, and thus this code
is based on Matlab documentation only. No compatibility testing has been

* Creating Galois Fields::
* Primitive Polynomials::
* Accessing Internal Fields::
* Known Problems::

@node Creating Galois Fields, Primitive Polynomials, , Galois Field Basics
@subsection Creating Galois Fields

To work with a Galois Field GF(2^M) in Octave, you must first create a variable
that Octave recognizes as a Galois Field. This is done with the function
@code{gf (@var{a}, @var{m})} as follows.

@example
octave:1> a = [0:7];
octave:2> b = gf (a, 4)
b =
GF(2^4) array. Primitive Polynomial = D^4+D+1 (decimal 19)

Array elements =

0   1   2   3   4   5   6   7
@end example

This creates an array @var{b} with 8 elements that Octave recognizes as a
Galois Field. The field is created with the default primitive polynomial for
the field GF(2^4). It can be verified that a variable is in fact a Galois
Field with the functions @code{isgalois} or @code{whos}.

@example
octave:3> isgalois (a)
ans = 0
octave:4> isgalois (b)
ans = 1
octave:5> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
a           1x8                         24  double
b           1x8                         32  galois

Total is 16 elements using 56 bytes
@end example

It is also possible to create a Galois Field with an arbitrary primitive
polynomial. However, if the polynomial is not a primitive polynomial of
the field, and error message is returned. For instance.

@example
octave:1> a = [0:7];
octave:2> b = gf (a, 4, 25)
b =
GF(2^4) array. Primitive Polynomial = D^4+D^3+1 (decimal 25)

Array elements =

0   1   2   3   4   5   6   7

octave:3> c = gf (a, 4, 21)
error: gf: primitive polynomial (21) of Galois Field must be irreducible
@end example

The function @code{gftable} is included for compatibility with Matlab. In
Matlab this function is used to create the lookup tables used to accelerate
the computations over the Galois Field and store them to a file. However
Octave stores these parameters for all of the fields currently in use and
so this function is not required, although it is silently accepted.

@node Primitive Polynomials, Accessing Internal Fields, Creating Galois Fields, Galois Field Basics
@subsection Primitive Polynomials

The function @code{gf (@var{a}, @var{m})} creates a Galois Field using the default primitive
polynomial. However there exists many possible primitive polynomials for most
Galois Fields. Two functions exist for identifying primitive polynomials,
@code{isprimitive} and @code{primpoly}. @code{primpoly (@var{m}, @var{opt})} is
used to identify the primitive polynomials of the fields GF(2^M). For example

@example
octave:1> primpoly (4)

Primitive polynomial(s) =

D^4+D+1

ans = 19
@end example

identifies the default primitive polynomials of the field GF(2^M), which
is the same as @code{primpoly (4, "min")}. All of the primitive polynomials
of a field can be identified with the function @code{primpoly (@var{m}, "all")}.
For example

@example
octave:1> primpoly (4, "all")

Primitive polynomial(s) =

D^4+D+1
D^4+D^3+1

ans =

19   25
@end example

while @code{primpoly (@var{m}, "max")} returns the maximum primitive polynomial
of the field, which for the case above is 25. The function @code{primpoly}
can also be used to identify the primitive polynomials having only a
certain number of non-zero terms. For instance

@example
octave:1> primpoly (5, 3)

Primitive polynomial(s) =

D^5+D^2+1
D^5+D^3+1

ans =

37   41
@end example

identifies the polynomials with only three terms that can be used as
primitive polynomials of GF(2^5). If no primitive polynomials existing
having the requested number of terms then @code{primpoly} returns an
empty vector. That is

@example
octave:1> primpoly (5, 2)
warning: primpoly: No primitive polynomial satisfies the given constraints

ans = [](1x0)
@end example

As can be seen above, @code{primpoly} displays the polynomial forms the
the polynomials that it finds. This output can be suppressed with the
"nodisplay" option, while the returned value is left unchanged.

@example
octave:1> primpoly (4, "all", "nodisplay")
ans =

19   25
@end example

@code{isprimitive (@var{a})} identifies whether the elements of @var{a} can
be used as primitive polynomials of the Galois Fields GF(2^M). Consider
as an example the fields GF(2^4). The primitive polynomials of these fields
must have an order m and so their integer representation must be between
16 and 31. Therefore @code{isprimitive} can be used in a similar manner to
@code{primpoly} as follows

@example
octave:1> find (isprimitive (16:31)) + 15
ans =

19   25
@end example

which finds all of the primitive polynomials of GF(2^4).

@subsection Accessing Internal Fields

Once a variable has been defined as a Galois Field, the parameters of the
field of this structure can be obtained by adding a suffix to the variable.
Valid suffixes are '.m', '.prim_poly' and '.x', which return the order of the
Galois Field, its primitive polynomial and the data the variable contains
respectively. For instance

@example
octave:1> a = [0:7];
octave:2> b = gf (a, 4);
octave:3> b.m
ans = 4
octave:4> b.prim_poly
ans = 19
octave:5> c = b.x;
octave:6> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
a           1x8                         24  double
b           1x8                         32  galois
c           1x8                         64  double

Total is 24 elements using 120 bytes
@end example

@c Note that if code compiled with GALOIS_DISP_PRIVATES then '.n', '.alpha_to'
@c and '.index_of' are also available. These give 2^m-1, the lookup table
@c and its inverse respectively.

Please note that it is explicitly forbidden to modify the Galois field by
accessing these variables. For instance

@example
octave:1> a = gf ([0:7], 3);
octave:2> a.prim_poly = 13;
@end example

is explicitly forbidden. The result of this will be to replace the
Galois array @var{a} with a structure @var{a} with a single element
called '.prim_poly'. To modify the order or primitive polynomial of a
field, a new field must be created and the data copied. That is

@example
octave:1> a = gf ([0:7], 3);
octave:2> a = gf (a.x, a.m, 13);
@end example

An important consideration in the use of the Galois Field package is
that many of the internal functions of Octave, such as @code{roots}, can
not accept Galois Fields as an input. This package therefore uses Octave
classes to @emph{overload} the internal Octave functions with equivalent
functions that work with Galois Fields, so that the standard function names
can be used.

The version of the function that is chosen is determined by the first
argument of the function. This is a temporary situation until the
Galois Field class constructor can be rewritten to allow the use of the
@code{superiorto} function to define the galois class with a higher
precedence. So, considering the @code{filter} function,
if the first argument is a @emph{Matrix}, then the normal version of
the function is called regardless of whether the other arguments of the
function are Galois vectors or not.

Other Octave functions work correctly with Galois Fields and so overloaded
versions are not necessary. This include such functions as @code{size} and
@code{polyval}.

It is also useful to use the '.x' option discussed in the previous section,
to extract the raw data of the Galois field for use with some functions. An
example is

@example
octave:1> a = minpol (gf (14, 5));
octave:2> b = de2bi (a.x, [], "left-msb");
@end example

converts the polynomial form of the minimum polynomial of 14 in GF(2^5) into
an integer. Finally help for the Galois specific versions of the functions
must explicitly call the correct function as

@example
octave:1> help @@galois/conv
@end example

@subsection Known Problems

Please review the following list of known problems with the Galois type
before reporting a bug against this package.

@table @asis

Saving a Galois variable to a file is as simple as

@example
octave:1> a = gf (@dots{});
octave:2> save a.mat a
@end example

@noindent
where @var{a} is any Galois variable. Galois variables can be saved in the
Octave binary and ASCII formats, as well as the HDF5 format. To load a
Galois variable from a file, the Galois type must already be registered to
the Octave interpreter prior to the call to @code{load}. If no Galois
variables have been created yet, you will have to do something like

@example
octave:1> dummy = gf (1);
@end example

@item Logarithm of zero does not return NaN
The logarithm of zero in a Galois field is not defined. However, to avoid
segmentation faults in later calculations the logarithm of zero is defined
as @code{2^@var{m} - 1}, whose value is not the logarithm of any other value
in the Galois field. A warning is also shown to tell the user about the
problem. For example

@example
octave:1> m = 3;
octave:2> a = log (gf ([0:2^m-1], m))
warning: log of zero undefined in Galois field
a =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

7   0   1   3   2   6   4   5
@end example

To fix this problem would require a major rewrite of all code, adding
an exception for the case of NaN to all basic operators. These
exceptions will certainly slow the code down.

@item Speed
The code was written piecemeal with no attention to optimization. Some
operations may be slower than they could be. Contributions are welcome.

@end table

@node Manipulating Galois Fields, , Galois Field Basics, Galois Fields
@section Manipulating Galois Fields

* Expressions manipulation and assignment::
* Unary operations::
* Arithmetic operations::
* Comparison operations::
* Polynomial manipulations::
* Linear Algebra::
* Signal Processing::

@node Expressions manipulation and assignment, Unary operations, , Manipulating Galois Fields
@subsection Expressions, manipulation and assignment

Galois variables can be treated in similar manner to other variables within
Octave. For instance Galois fields can be accessed using index expressions
in a similar manner to all other Octave matrices. For example

@example
octave:1> a = gf ([[0:7]; [7:-1:0]], 3)
a =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

0   1   2   3   4   5   6   7
7   6   5   4   3   2   1   0

octave:2> b = a(1,:)
b =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

0   1   2   3   4   5   6   7
@end example

Galois arrays can equally use indexed assignments. That is, the data
in the array can be partially replaced, on the condition that the two
fields are identical. An example is

@example
octave:1> a = gf (ones (2, 8), 3);
octave:2> b = gf (zeros (1, 8), 3);
octave:3> a(1,:) = b
a =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

0   0   0   0   0   0   0   0
1   1   1   1   1   1   1   1
@end example

Implicit conversions between normal matrices and Galois arrays are possible.
For instance data can be directly copied from a Galois array to a real matrix
as follows.

@example
octave:1> a = gf (ones (2, 8), 3);
octave:2> b = zeros (2, 8);
octave:3> b(2,:) = a(2,:)
b =

0   0   0   0   0   0   0   0
1   1   1   1   1   1   1   1
@end example

The inverse is equally possible, with the proviso that the data in the matrix
is valid in the Galois field. For instance

@example
octave:1> a = gf ([0:7], 3);
octave:2> a(1) = 1;
@end example

is valid, while

@example
octave:1> a = gf ([0:7], 3);
octave:2> a(1) = 8;
@end example

is not, since 8 is not an element of GF(2^3). This is a basic rule of
manipulating Galois arrays. That is matrices and scalars can be used in
conjunction with a Galois array as long as they contain valid data
within the Galois field. In this case they will be assumed to be of the
same field.

Galois arrays can also be concatenated with real matrices or with other
Galois arrays in the same field. For example

@example
octave:1> a = [gf([0:7], 3); gf([7:-1:0], 3)];
octave:2> b = [a, a];
octave:3> c = [a, eye(2)];
octave:3> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
a           2x8                         64  galois
b           2x16                       128  galois
c           2x10                        80  galois

Total is 68 elements using 272 bytes
@end example

Other basic manipulations of Galois arrays are

@table @code
@item isempty
Returns true if the Galois array is empty.

@item size
Returns the number of rows and columns in the Galois array.

@item length
Returns the length of a Galois vector, or the maximum of rows or columns
of Galois arrays.

@item find
Find the indexes of the non-zero elements of a Galois array.

@item diag
Create a diagonal Galois array from a Galois vector, or extract a diagonal
from a Galois array.

@item reshape
Change the shape of the Galois array.

@end table

@node Unary operations, Arithmetic operations, Expressions manipulation and assignment, Manipulating Galois Fields
@subsection Unary operations

The same unary operators that are available for normal Octave matrices are
also available for Galois arrays. These operations are

@table @code
@item +@var{x}
Unary plus. This operator has no effect on the operand.

@item -@var{x}
Unary minus. Note that in a Galois Field this operator also has no effect
on the operand.

@item !@var{x}
Returns true for zero elements of Galois Array.

@item @var{x}'
Complex conjugate transpose. As the Galois Field only contains integer
values, this is equivalent to the transpose operator.

@item @var{x}.'
Transpose of the Galois array.
@end table

@node Arithmetic operations, Comparison operations, Unary operations, Manipulating Galois Fields
@subsection Arithmetic operations

The available arithmetic operations on Galois arrays are the same as on
other Octave matrices. It should be noted that both operands must be in
the same Galois Field. If one operand is a Galois array and the second is
a matrix or scalar, then the second operand is silently converted to the
same Galois Field. The element(s) of these matrix or scalar must however
be valid members of the Galois field. Thus

@example
octave:1> a = gf ([0:7], 3);
octave:2> b = a + [0:7];
@end example

is valid, while

@example
octave:1> a = gf ([0:7], 3);
octave:2> b = a + [1:8];
@end example

is not, since 8 is not a valid element of GF(2^3). The available arithmetic
operators are

@table @code
@item @var{x} + @var{y}
Addition. If both operands are Galois arrays or matrices, the number of rows
and columns must both agree.  If one operand is a is a Galois array with a
single element or a scalar, its value is added to all the elements of the
other operand. The @code{+} operator on a Galois Field is equivalent to an
exclusive-or on normal integers.

@item @var{x} .+ @var{y}
Element by element addition. This operator is equivalent to @code{+}.

@item @var{x} - @var{y}
As both @code{+} and @code{-} in a Galois Field are equivalent to an
exclusive-or for normal integers, @code{-} is equivalent to the @code{+}
operator

@item @var{x} .- @var{y}
Element by element subtraction. This operator is equivalent to @code{-}.

@item @var{x} * @var{y}
Matrix multiplication.  The number of columns of @var{x} must agree
with the number of rows of @var{y}.

@item @var{x} .* @var{y}
Element by element multiplication.  If both operands are matrices, the
number of rows and columns must both agree.

@item @var{x} / @var{y}
Right division.  This is conceptually equivalent to the expression

@example
(inverse (y') * x')'
@end example

@noindent
but it is computed without forming the inverse of @var{y'}.

If the matrix is singular then an error occurs. If the matrix is
under-determined, then a particular solution is found (but not minimum
norm). If the solution is over-determined, then an attempt is made
to find a solution, but this is not guaranteed to work.

@item @var{x} ./ @var{y}
Element by element right division.

@item @var{x} \ @var{y}
Left division.  This is conceptually equivalent to the expression

@example
inverse (x) * y
@end example

@noindent
but it is computed without forming the inverse of @var{x}.

If the matrix is singular then an error occurs. If the matrix is
under-determined, then a particular solution is found (but not minimum
norm). If the solution is over-determined, then an attempt is made
to find a solution, but this is not guaranteed to work.

@item @var{x} .\ @var{y}
Element by element left division.  Each element of @var{y} is divided
by each corresponding element of @var{x}.

@item @var{x} ^ @var{y}
@itemx @var{x} ** @var{y}
Power operator.  If @var{x} and @var{y} are both scalars, this operator
returns @var{x} raised to the power @var{y}.  Otherwise @var{x} must
be a square matrix raised to an integer power.

@item @var{x} .^ @var{y}
@item @var{x} .** @var{y}
Element by element power operator.  If both operands are matrices, the
number of rows and columns must both agree.

@end table

@node Comparison operations, Polynomial manipulations, Arithmetic operations, Manipulating Galois Fields
@subsection Comparison operations

Galois variables can be tested for equality in the usual manner. That is

@example
octave:1> a = gf ([0:7], 3);
octave:2> a == ones (1, 8)
ans =

0   1   0   0   0   0   0   0

octave:3> a != zeros (1, 8)
ans =

0   1   1   1   1   1   1   1
@end example

Likewise, Galois vectors can be tested against scalar values (whether they are
Galois or not). For instance

@example
octave:4> a == 1
ans =

0   1   0   0   0   0   0   0
@end example

To test if any or all of the values in a Galois array are non-zero, the
functions @code{any} and @code{all} can be used as normally.

In addition the comparison operators @code{>}, @code{>=}, @code{<} and
@code{<=} are available. As elements of the Galois Field are modulus
2^@var{m}, all elements of the field are both greater than and less than
all others at the same time.Thus these comparison operators don't make
that much sense and are only included for completeness. The comparison is
done relative to the integer value of the Galois Field elements.

@node Polynomial manipulations, Linear Algebra, Comparison operations, Manipulating Galois Fields
@subsection Polynomial manipulations

A polynomial in GF(2^M) can be expressed as a vector in GF(2^M). For instance
if @var{a} is the @emph{primitive element}, then the example

@example
octave:1> poly = gf ([2, 4, 5, 1], 3);
@end example

represents the polynomial

@tex
$$poly = a x^3 + a^2 x^2 + (a^2 + 1) x + 1$$
@end tex
@ifnottex
@example
poly = @var{a} * x^3 + @var{a}^2 * x^2 + (@var{a}^2 + 1) * x + 1
@end example
@end ifnottex

Arithmetic can then be performed on these vectors. For instance to add
to polynomials an example is

@example
octave:1> poly1 = gf ([2, 4, 5, 1], 3);
octave:2> poly2 = gf ([1, 2], 3);
octave:3> sumpoly = poly1 + [0, 0, poly2]
sumpoly =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

2   4   4   3
@end example

Note that @var{poly2} must be zero padded to the same length as poly1 to
allow the addition to take place.

Multiplication and division of Galois polynomials is equivalent to convolution
and de-convolution of vectors of Galois elements. Thus to multiply two
polynomials in GF(2^3).

@example
octave:4> mulpoly = conv (poly1, poly2)
mulpoly =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

2   0   6   0   2
@end example

Likewise the division of two polynomials uses the de-convolution function
as follows

@example
octave:5> [poly, remd] = deconv (mulpoly, poly2)
poly =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

2   4   5   1

remd =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

0   0   0   0   0
@end example

Note that the remainder of this division is zero, as we performed the inverse
operation to the multiplication.

To evaluate a polynomial for a certain value in GF(2^M), use the Octave
function @code{polyval}.

@example
octave:1> poly1 = gf ([2, 4, 5, 1], 3);  ## a*x^3+a^2*x^2+(a^2+1)*x+1
octave:2> x0 = gf ([0, 1, 2], 3);
octave:3> y0 = polyval (poly1, x0);
y0 =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

1   2   0

octave:4> a = gf (2, 3);               ## The primitive element
octave:5> y1 = a .* x0.^3 + a.^2 .* x0.^2 + (a.^2 + 1) .* x0 + 1
y1 =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

1   2   0
@end example

It is equally possible to find the roots of Galois polynomials with the
@code{roots} function. Using the polynomial above over GF(2^3), we can
find its roots in the following manner

@example
octave:1> poly1 = gf ([2, 4, 5, 1], 3);
octave:2> root1 = roots (poly1)
root1 =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

2
5
5
@end example

Thus the example polynomial has 3 roots in GF(2^3) with one root of
multiplicity 2. We can check this answer with the @code{polyval} function
as follows

@example
octave:3> check1 = polyval (poly1, root1)
check1 =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

0
0
0
@end example

which as expected gives a zero vector. It should be noted that both the
number of roots and their value, will depend on the chosen field. Thus
for instance

@example
octave:1> poly3 = gf ([2, 4, 5, 1], 3, 13);
octave:2> root3 = roots (poly3)
root3 =
GF(2^3) array. Primitive Polynomial = D^3+D^2+1 (decimal 13)

Array elements =

5
@end example

shows that in the field GF(2^3) with a different primitive polynomial,
has only one root exists.

The minimum polynomial of an element of GF(2^M) is the minimum degree
polynomial in GF(2), excluding the trivial zero polynomial, that has
that element as a root. The fact that the minimum polynomial is in GF(2)
means that its coefficients are one or zero only. The @code{minpol}
function can be used to find the minimum polynomial as follows

@example
octave:1> a = gf (2, 3);               ## The primitive element
octave:2> b = minpol (a)
b =
GF(2) array.

Array elements =

1   0   1   1
@end example

Note that the minimum polynomial of the primitive element is the primitive
polynomial. Elements of GF(2^M) sharing the same minimum polynomial form a
partitioning of the field. This partitioning can be found with the
@code{cosets} function as follows

@example
octave:1> c = cosets (3)
c =
@{
[1,1] =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

1

[1,2] =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

2   4   6

[1,3] =
GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)

Array elements =

3   5   7

@}
@end example

which returns a cell array containing all of the elements of the GF(2^3),
partitioned into groups sharing the same minimum polynomial. The function
@code{cosets} can equally accept a second argument defining the primitive
polynomial to use in its calculations (i.e. @code{cosets (@var{a}, @var{p})}).

@node Linear Algebra, Signal Processing, Polynomial manipulations, Manipulating Galois Fields
@subsection Linear Algebra

The basic linear algebra operation of this package is the LU factorization
of a Galois array. That is the Galois array @var{a} is factorized in the
following way

@example
octave:2> [l, u, p] = lu (a)
@end example

such that @code{@var{p} * @var{a} = @var{l} * @var{u}}. The matrix @var{p}
contains row permutations of @var{a}, such that @var{l} and @var{u} are
strictly upper and low triangular. The Galois array @var{a} can be
rectangular.

All other linear algebra operations within this package are based on this
LU factorization of a Galois array. An important consequence of this is that
no solution can be found for singular matrices, only a particular solution
will be found for under-determined systems of equation and the solution found
for over-determined systems is not always correct. This is identical to the
way Matlab R13 treats this.

For instance consider the under-determined linear equation

@example
octave:1> A = gf ([2, 0, 3, 3; 3, 1, 3, 1; 3, 1, 1, 0], 2);
octave:2> b = [0:2]';
octave:3> x = A \ b;
@end example

gives the solution @code{@var{x} = [2, 0, 3, 2]}. There are in fact 4
possible solutions to this linear system; @code{@var{x} = [3, 2, 2, 0]},
@code{@var{x} = [0, 3, 1, 1]}, @code{@var{x} = [2, 0, 3, 2]} and
@code{@var{x} = [1, 1, 0, 3]}. No particular selection criteria are
applied to the chosen solution.

In addition the fact that singular matrices are not treated, unless you
know the matrix is not singular, you should test the determinant of the
matrix prior to solving the linear system. For instance

@example
octave:1> A = gf (floor (2^m * rand (3)), 2);
octave:2> b = [0:2]';
octave:3> if (det (A) != 0); x = A \ b; y = b' / A; endif;
octave:4> r = rank (A);
@end example

solves the linear systems @code{@var{A} * @var{x} = @var{b}} and
@code{@var{y} * @var{A} = @var{b}}. Note that you do not need to take
into account rounding errors in the determinant, as the determinant can
only take values within the Galois Field. So if the determinant equals
zero, the array is singular.

@node Signal Processing, , Linear Algebra, Manipulating Galois Fields
@subsection Signal Processing with Galois Fields

Signal processing functions such as filtering, convolution, de-convolution
and Fourier transforms can be performed over Galois Fields. For instance
the @code{filter} function can be used with Galois vectors in the same
manner as usual. For instance

@example
octave:1> b = gf ([2, 0, 0, 1, 0, 2, 0, 1], 2);
octave:2> a = gf ([2, 0, 1, 1], 2);
octave:3> x = gf ([1, zeros(1, 20)], 2);
octave:4> y = filter (b, a, x)
y =
GF(2^2) array. Primitive Polynomial = D^2+D+1 (decimal 7)

Array elements =

1   0   3   0   2   3   1   0   1   3   3   1   0   1   3   3   1   0   1   3   3
@end example

gives the impulse response of the filter defined by @var{a} and @var{b}.

Two equivalent ways are given to perform the convolution of two Galois
vectors. Firstly the function @code{conv} can be used, or alternatively
the function @code{convmtx} can be used. The first of these function is
identical to the convolution function over real vectors, and has been
described in the section about multiplying two Galois polynomials.

In the case where many Galois vectors will be convolved with the same
vector, the second function @code{convmtx} offers an alternative method
to calculate the convolution. If @var{a} is a column vector and @var{x}
is a column vector of length @var{n}, then

@example
octave:1> m = 3;
octave:2> a = gf (floor (2^m * rand (4, 1)), m);
octave:3> b = gf (floor (2^m * rand (4, 1)), m);
octave:4> c0 = conv (a, b)';
octave:5> c1 = convmtx (a, length (b)) * b;
octave:6> check = all (c0 == c1)
check = 1
@end example

shows the equivalence of the two functions. The de-convolution function has
been previously described above.

The final signal processing function available in this package are the
functions to perform Fourier transforms over a Galois field. Three
functions are available, @code{fft}, @code{ifft} and @code{dftmtx}. The
first two functions use the third to perform their work. Given an element
@var{a} of the Galois field GF(2^M), @code{dftmtx} returns the @code{2^M - 1}
square matrix used in the Fourier transforms with respect to @var{a}. The
minimum polynomial of @var{a} must be primitive in GF(2^M). In the case of
the @code{fft} function @code{dftmtx} is called with the primitive element of
the Galois Field as an argument. As an example

@example
octave:1> m = 4;
octave:2> n = 2^m - 1;
octave:2> alph = gf (2, m);
octave:3> x = gf (floor (2^m * rand (n, 1)), m);
octave:4> y0 = fft (x);
octave:5> y1 = dftmtx (alph) * x;
octave:6> z0 = ifft (y0);
octave:7> z1 = dftmtx (1/alph) * y1;
octave:8> check = all (y0 == y1) & all (z0 == x) & all (z1 == x)
check = 1
@end example

In all cases, the length of the vector to be transformed must be
@code{2^M -1}. As the @code{dftmtx} creates a matrix representing the
Fourier transform, to limit the computational task only Fourier
transforms in GF(2^M), where M is less than or equal to 8, can be treated.

@node Function Reference, , Galois Fields, Top
@chapter Function Reference

@REFERENCE_SECTION(Communications)

@bye