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.

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
===============

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 randint, that creates a random matrix of equi-probable integers in a desired range. For example

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

        0   1   0
       -1  -1   1
        0   1   1

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 randint allows the seed-value of the random number generator to be set. For instance

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

        0   1   1
        0  -1   0
        1  -1  -1

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 M, randint returns the equi-probable integers in the range [0:2^M].

The function randsrc differs from 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

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

are both legal, while

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

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

     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

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 randint, 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 randerr allows a matrix of random bit errors to be created, for binary encoded messages. By default, randerr creates exactly one errors per row, flagged by a non-zero value in the returned matrix. That is

     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

The number of errors per row can be specified as the third argument to 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

     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

creates a vector 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, 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 wgn and awgn create and add white Gaussian noise to continuous signals. The function wgn creates a matrix of white Gaussian noise of a certain power. A typical call to wgn is then

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

Which creates a 10-by-10 matrix of noise with a root mean squared power of 0dBW relative to an impedance of 1 Ohm.

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

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

The reference impedance and units of power to the function wgn can however be modified, for 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

All produce a 1W signal referenced to a 50 Ohm impedance. Matlab uses the misnomer "dB" for "dBW", and therefore "dB" is an accepted type for wgn and will be treated as for "dBW".

In all cases, the returned matrix V, will be related to the input power P and the impedance Z as

P = sum(V(:) .^ 2 ) / IMP Watts

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

     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

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, 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

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

will always produce the same set of data.

The final function to deal with the creation of random signals is 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 awgn is

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

which adds noise with a 10dB signal-to-noise ratio to the measured power in the desired signal. By default 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 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 addition of the noise.

The final argument to awgn defines the definition of the power and signal-to-noise ratio in a similar manner to 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 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 awgn to allow repetition of the same scenario. That is

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

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

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 biterr and 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

Both matrices
     In this case both matrices must be the same size and then by default the the return values are the overall number of errors and the overall error rate.

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.

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.

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 biterr and symerr, consider the 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

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 eyediagram and scatterplot have different appearance, the information presented is similar and so are their inputs. The difference between eyediagram and scatterplot is that eyediagram segments the data into time intervals and plots the in-phase and quadrature components of the signal against this time interval. While scatterplot uses a parametric plot of quadrature versus in-phase components.

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

A real vector
     In this case the signal is assumed to be real and represented by the vector X.

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.

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.

An example of the use of the function eyediagram is

     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);

Similarly an example of the use of the function scatterplot is

     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> hold off;
     octave:9> scatterplot(noisy,1,0,"b");
     octave:10> hold on;
     octave:11> scatterplot(noisy,ovsp,0,"r+");

Source Coding
*************

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, lloyds creates an optimal mapping of the continous signal to a fixed number of levels and quantiz performs the actual quantization.

The set of quantization points to use is represented by a partitioning table (TABLE) of the data and the signal levels (CODES to which they are mapped. The partitioning TABLE is monotonicly 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.

Table 1: Table quantization partitioning and coding

  x < table(1)                 codes(1)
  table(1) <= x < table(2)     codes(2)
  ...                          ...
  table(i-1) <= x < table(i)   codes(i)
  ...                          ...
  table(n-1) <= x < table(n)   codes(n)
  table(n-1) <= x              codes(n+1)

These partition and coding tables can either be created by the user of using the function lloyds. For instance the use of a linear mapping can be seen in the following 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);

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

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

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

PCM Coding
==========

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

Arithmetic Coding
=================

The arithmetic coding functions arithenco and arithdeco are not yet implemented.

Dynamic Range Compression
=========================

The final source coding function is 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. compand can then be used to compress the signal prior to quantization, with the signal being expanded afterwards. For 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

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

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
============

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 K bits of data are considered to represent a single message symbol. These K bits are coded into 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

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 K in length.

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 K in the case of a uncoded message or N in the case of a coded message.

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.

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

     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

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.

Binary Block Codes
==================

All of the codes presented here can be characterized by their

Generator Matrix
     A K-by-N matrix G to generate the codewords C from the messages T by the matrix multiplication C = T * G.

Parity Check Matrix
     A 'N-K'-by-N matrix H to check the parity of the received symbols. If H * R = S != 0, then an error has been detected. S can be used with the syndrome table to correct this error

Syndrome Table
     A 2^K-by-N matrix 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 R = mod(T + E, 2), then the error vector E is ST(S).

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

      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

or

       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)

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 hammgen, cyclgen and cyclpoly. The function hammgen creates parity check and generator matrices for Hamming codes, while cyclpoly and cyclgen create generator polynomials and matrices for generic cyclic codes. An example of their use is

     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));

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 syndtable, in the following manner

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

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

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

It should be noted that for large values of 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 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 encode and 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

'linear'
     Generic linear block codes

'cyclic'
     Cyclic linear block codes

'hamming'
     Hamming codes

'bch'
     Bose Chaudhuri Hocquenghem (BCH) block codes

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 encode and decode must be told the format of the messages in the following manner.

     octave:1> m = 3;
     octave:2> n = 7;
     ocatve: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");

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

Except for the BCH codes, the function encode and decode internally create the generator, parity check and syndrome table matrices.  Therefore if repeated calls to encode and decode are made it will often be faster to create these matrices externally, and pass them as an argument. For 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 !!!
     end

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

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.

Called without arguments, bchpoly returns a table of valid BCH error correcting codes and their error-correction capability. The first returned column of bchpoly is the codeword length, the second the message length and the third the error correction capability of the code. Called with one argument, 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 2^M-1, where M is an integer. However if [N,K] is a valid BCH code, then it is also possible to use a shortened BCH form of the form [N-X,K-X].

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

     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

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

1 + X ^ 4 + X ^ 6 + X ^ 7 + X ^ 8

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.

     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

It is recommend not to convert the generator polynomials created by 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 encode and decode functions previously discussed, the user can directly use the low-level BCH functions bchenco and bchdeco. In this case the messages must be in the format of a binary matrix with K columns

     octave:1> n = 31;
     octave:2> pgs = bchpoly(n);
     octave:3> pg = pgs(floor(rand(1,1)*(size(pgs,1) + 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);

Reed-Solomon Codes
==================

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

     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

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

Creating and Decoding Messages
------------------------------

The Reed-Solomon encoding function requires at least three arguments. The first MSG is the message in encodes, the second is N the codeword length and K is the message length. Therefore MSG must have K columns and the output will have 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 N = 2^M - 1. Another important parameter when determining the behavior of the Reed-Solomon coder is the primitive polynomial of the Galois Field (see gf). Thus the messages

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

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 N - 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 rsenc and rsdec. By default the parity symbols are placed after the message.

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

     G = (X - A^(B*S)) * (X - A^((B+1)*S)) * ... * (X - A^((B+2*T-1)*S)).

where T is (N-K)/2, A is the primitive element of the Galois Field, B is the first consecutive root, and S is the step between roots. Generator polynomial of this form are constructed by rsgenpoly and can be passed to both rsenc and rsdec. It is also possible to pass the B and S values directly to rsenc and rsdec. In the case of rsdec passing B and S can make the decoding faster.

Consider the example below.

     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:13> nerr'
     ans =

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

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

       0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1

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 NERR gives the number of errors corrected. In the case of 17 introduced errors 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.

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 N = 2^M = 1. Although, the underlying Reed-Solomon coding must operate over valid codeword length, there are sometimes reasons to assume the 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

     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

Convolutional Coding
********************

To be written.

Modulations
***********

To be written.

Currently have functions amodce, ademodce, apkconst, demodmap, modmap, qaskdeco and qaskenco.

Special Filters
***************

To be written.

Galois Fields
*************

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 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 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

Integer      Binary       Element of GF(2^M)
0            000         0
1            001         1
2            010         A
3            011         A + 1
4            100         A^2
5            101         A^2 + 1
6            110         A^2 + A
7            111         A^2 + A + 1

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 performed and so I am most interested in comments about compatibility at the e-mail address dbateman@free.fr.

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 gf(A,M) as follows.

     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

This creates an array 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 isgalois or whos.

     octave:3> isgalois(a)
     ans = 0
     octave:4> isgalois(b)
     ans = 1
     octave:5> whos

     *** local user variables:

     prot  type                       rows   cols  name
     ====  ====                       ====   ====  ====
      rwd  matrix                        1      8  a
      rwd  galois                        1      8  b

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.

     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: primitive polynomial (21) of Galois Field must be irreducible
     error: unable to initialize Galois Field
     error: evaluating assignment expression near line 3, column 3

The function 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.

Primitive Polynomials
---------------------

The function gf(A,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, isprimitive and primpoly. primpoly(M,OPT) is used to identify the primitive polynomials of the fields GF(2^M). For example

     octave:1> primpoly(4)

     Primitive polynomial(s) =

     D^4+D+1

     ans = 19

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

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

     Primitive polynomial(s) =

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

     ans =

       19  25

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

     octave:1> primpoly(5, 3)

     Primitive polynomial(s) =

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

     ans =

       37  41

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 primpoly returns an empty vector. That is

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

     ans = [](1x0)

As can be seen above, 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.

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

       19  25

isprimitive(A) identifies whether the elements of 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 isprimitive can be used in a similar manner to primpoly as follows

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

       19  25

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

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

     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

     *** local user variables:

     prot  type                       rows   cols  name
     ====  ====                       ====   ====  ====
      rwd  matrix                        1      8  a
      rwd  galois                        1      8  b
      rwd  matrix                        1      8  c

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

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

is explicitly forbidden. The result of this will be to replace the Galois array A with a structure 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

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

Function Overloading
--------------------

An important consideration in the use of the Galois Field package is that many of the internal functions of Octave, such as roots, can not accept Galois Fields as an input. This package therefore uses the dispatch function of Octave-Forge to overload the internal Octave functions with equivalent functions that work with Galois Fields, so that the standard function names can be used. However, at any time the Galois field specific version of the function can be used by explicitly calling its function name. The correspondence between the internal function names and the Galois Field versions is as follows

conv   -   gconv,
convmtx   -   gconvmtx,
diag   -   gdiag,
deconv   -   gdeconv,
det   -   gdet,
exp   -   gexp,
filter   -   gfilter,
inv    -   ginv,
log   -   glog,
lu    -   glu,
prod   -   gprod,
reshape   -   greshape,
rank    -   grank,
roots   -   groots,
sum   -   gsum,
sumsq   -   gsumsq.

The version of the function that is chosen is determined by the first argument of the function. So, considering the filter function, if the first argument is a 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 size and 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

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

converts the polynomial form of the minimum polynomial of 14 in GF(2^5) into an integer.

Known Problems
--------------

Before reporting a bug compare it to this list of known problems

Concatenation
     For versions of Octave prior to 2.1.58, the concatenation of Galois arrays returns a Matrix type. That is [gf([1, 0],m) gf(1, m)] returns a matrix went it should return another Galois array. The workaround is to explicitly convert the  returned value back to the correct Galois field using gf([gf([1, 0],m) gf(1,m)],m).

     Since Octave version 2.1.58, [gf([1, 0],m) gf(1, m)] returns another Galois array as expected.

Saving and loading Galois variables

     Saving of Galois variables is only implemented in versions of octave later than 2.1.53. If you are using a recent version of octave then saving a Galois variable is as simple as

          octave:2> save a.mat a

     where A is a Galois variable. To reload the variable within octave, the Galois type must be installed prior to a call to load. That is

          octave:1> dummy = gf(1);
          octave:2> load a.mat

     With versions of octave later than 2.1.53, Galois variables can be saved in the octave binary and ascii formats, as well as the HDF5 format. If you are using an earlier version of octave, you can not directly save a Galois variable. You can however save the information it contains and reconstruct the data afterwards by doing something like

          octave:2> x = a.x; m = a.m; p = a.prim_poly;
          octave:3> save a.mat x m p;

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 2^M-1, whose value is not the logarithm of any other value in the Galois field. A warning is however printed to tell the user about the problem. For 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

     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.

Speed
     The code was written piece-meal with no attention to optimum code. Now that I have something working I should probably go back and tidy the code up, optimizing it at the same time.

Manipulating Galois Fields
==========================

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

     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

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

     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

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.

     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

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

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

is valid, while

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

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.

As Octave supports concatenation of typed matrices only for version 2.1.58 and later, matrix concatenation will force the Galois array back to a normal matrix for earlier version. For instance for Octave 2.1.58 and later.

     octave:1> a = [gf([0:7],3); gf([7:-1:0],3)];
     octave:2> b = [a, a];
     octave:3> whos

     *** local user variables:

       Prot Name       Size  Bytes  Class
       ==== ====       ====  =====  =====
        rwd a          2x8      64  galois
        rwd b          2x16    128  galois

     Total is 49 elements using 192 bytes

and for  previous versions of Octave

     octave:1> a = [gf([0:7],3); gf([7:-1:0],3)];
     octave:2> b = [a, a];
     octave:3> whos

     *** local user variables:

     prot  type                       rows   cols  name
     ====  ====                       ====   ====  ====
      rwd  matrix                        2      8  a
      rwd  matrix                        2     16  b

This has the implication that many of the scripts included with Octave that should work with Galois fields, won't work correctly for versions earlier than 2.1.58. If you wish to concatenate Galois arrays with earlier versions, use the syntax

     octave:1> a = gf([0:7],3);
     octave:2> b = gf([a, a], a.m, a.prim_poly);

which explicitly reconverts B to the correct Galois Field. Other basic manipulations of Galois arrays are

isempty
     Returns true if the Galois array is empty.

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

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

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

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

reshape
     Change the shape of the Galois array.

Unary operations
----------------

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

+X
     Unary plus. This operator has no effect on the operand.

-X
     Unary minus. Note that in a Galois Field this operator also has no effect on the operand.

!X
     Returns true for zero elements of Galois Array.

X'
     Complex conjugate transpose. As the Galois Field only contains integer values, this is equivalent to the transpose operator.

X.'
     Transpose of the Galois array.

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

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

is valid, while

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

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

X + 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 + operator on a Galois Field is equivalent to an exclusive-or on normal integers.

X .+ Y
     Element by element addition. This operator is equivalent to +.

X - Y
     As both + and - in a Galois Field are equivalent to an exclusive-or for normal integers, - is equivalent to the + operator

X .- Y
     Element by element subtraction. This operator is equivalent to -.

X * Y
     Matrix multiplication.  The number of columns of X must agree with the number of rows of Y.

X .* Y
     Element by element multiplication.  If both operands are matrices, the number of rows and columns must both agree.

X / Y
     Right division.  This is conceptually equivalent to the expression

          (inverse (y') * x')'

     but it is computed without forming the inverse of 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.

X ./ Y
     Element by element right division.

X \ Y
     Left division.  This is conceptually equivalent to the expression

          inverse (x) * y

     but it is computed without forming the inverse of 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.

X .\ Y
     Element by element left division.  Each element of Y is divided by each corresponding element of X.

X ^ Y
X ** Y
     Power operator.  If X and Y are both scalars, this operator returns X raised to the power Y.  Otherwise X must be a square matrix raised to an integer power.

X .^ Y

X .** Y
     Element by element power operator.  If both operands are matrices, the number of rows and columns must both agree.

Comparison operations
---------------------

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

     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

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

     octave:4> a == 1
     ans =

       0  1  0  0  0  0  0  0

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

In addition the comparison operators >, >=, < and <= are available. As elements of the Galois Field are modulus 2^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.

Polynomial manipulations
------------------------

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

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

represents the polynomial

     poly = A * x^3 + A^2 * x^2 + (A^2 + 1) * x + 1

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

     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

Note that 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).

     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

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

     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

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 polyval.

     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

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

     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

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

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

     Array elements =

       0
       0
       0

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

     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

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 minpol function can be used to find the minimum polynomial as follows

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

     Array elements =

       1  0  1  1

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 cosets function as follows

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

       Array elements =

         1

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

       Array elements =

         2  4  6

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

       Array elements =

         3  5  7

     }

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

Linear Algebra
--------------

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

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

such that P * A = L * U. The matrix P contains row permutations of A, such that L and U are strictly upper and low triangular. The Galois array 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

     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;

gives the solution X = [2, 0, 3, 2]. There are in fact 4 possible solutions to this linear system; X = [3, 2, 2, 0], X = [0, 3, 1, 1], X = [2, 0, 3, 2] and 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

     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; end;
     octave:4> r = rank(A);

solves the linear systems A * X = B and Y * A = 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.

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 filter function can be used with Galois vectors in the same manner as usual. For instance

     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

gives the impulse response of the filter defined by A and B.

Two equivalent ways are given to perform the convolution of two Galois vectors. Firstly the function conv can be used, or alternatively the function 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 convmtx offers an alternative method to calculate the convolution. If A is a column vector and X is a column vector of length N, then

     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

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, fft, ifft and dftmtx. The first two functions use the third to perform their work. Give an element A of the Galois field GF(2^M), dftmtx returns the 2^M - 1 square matrix used in the Fourier transforms with respect to A. The minimum polynomial of A must be primitive in GF(2^M). In the case of the fft function dftmtx is called with the the primitive element of the Galois Field as an argument. As an 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

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

Function Reference
******************

_SECTION(Communications)