\input texinfo @setfilename comms.info @settitle @titlepage @title Communications Toolbox for Octave @subtitle March 2003 @author David Bateman @author Laurent Mazet @author Paul Kienzle @page @vskip 0pt plus 1filll Copyright @copyright{} 2003 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 resulting derived work is distributed under the terms of a permission 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 @ifinfo @node Top, Introduction @top @end ifinfo @menu * Introduction:: * Random Signals:: * Source Coding:: * Block Coding:: * Convolutional Coding:: * Modulations:: * Special Filters:: * Galois Fields:: * Function Reference:: @end menu @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. @menu * Signal Creation:: * Signal Analysis:: @end menu @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 @dfn{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 @dfn{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}, @dfn{randint} returns the equi-probable integers in the range @iftex @tex $[0:2^m-1]$. @end tex @end iftex @ifinfo [0:2^@var{m}]. @end ifinfo The function @dfn{randsrc} differs from @dfn{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 @dfn{randint}, @dfn{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 @dfn{randerr} allows a matrix of random bit errors to be created, for binary encoded messages. By default, @dfn{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 @dfn{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, @dfn{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 @dfn{wgn} and @dfn{awgn} create and add white Gaussian noise to continuous signals. The function @dfn{wgn} creates a matrix of white Gaussian noise of a certain power. A typical call to @dfn{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 @iftex @tex $1\Omega$. @end tex @end iftex @ifinfo 1 Ohm. @end ifinfo 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 @dfn{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 @iftex @tex $50\Omega$. @end tex @end iftex @ifinfo 50 Ohm @end ifinfo impedance. Matlab uses the misnomer "dB" for "dBW", and therefore "dB" is an accepted type for @dfn{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 @iftex @tex $$p = {\\sum_i \\sum_j v(i,j)^2 \\over Z} Watts$$ @end tex @end iftex @ifinfo @var{p} = sum(@var{v}(:) .^ 2 ) / @var{imp} Watts @end ifinfo By default @dfn{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, @dfn{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 @dfn{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 @dfn{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 @iftex 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 iftex which adds noise with a 10dB signal-to-noise ratio to the measured power in the desired signal. By default @dfn{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 @dfn{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 @dfn{awgn} defines the definition of the power and signal-to-noise ratio in a similar manner to @dfn{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 @dfn{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 @dfn{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 @dfn{biterr} and @dfn{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 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 @dfn{biterr} and @dfn{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 @dfn{eyediagram} and @dfn{scatterplot} have different appearance, the information presented is similar and so are their inputs. The difference between @dfn{eyediagram} and @dfn{scatterplot} is that @dfn{eyediagram} segments the data into time intervals and plots the in-phase and quadrature components of the signal against this time interval. While @dfn{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 @dfn{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 @iftex which produces a eye-diagram of a noisy signal as seen in Figure 2. Similarly an example of the use of the function @dfn{scatterplot} is @center @image{eyediagram, , } @center Figure 2: Eye-diagram of a QPSK like signal with 15dB signal-to-noise ratio @end iftex @ifinfo Similarly an example of the use of the function @dfn{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> hold off; octave:9> scatterplot(noisy,1,0,"b"); octave:10> hold on; octave:11> scatterplot(noisy,ovsp,0,"r+"); @end example @iftex 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 iftex @node Source Coding, Block Coding, Random Signals, Top @chapter Source Coding @menu * Quantization:: * PCM Coding:: * Arithmetic Coding:: * Dynamic Range Compression:: @end menu @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, @dfn{lloyds} creates an optimal mapping of the continous signal to a fixed number of levels and @dfn{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 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. @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 ... @tab ... @tab @item @tab table(i-1) <= x < table(i) @tab codes(i) @tab @item @tab ... @tab ... @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 @dfn{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 @dfn{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 @dfn{lloyds}. It should be noted that the mapping given by @dfn{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 @dfn{lloyds} will be sub-optimal. @node PCM Coding, Arithmetic Coding, Quantization, Source Coding @section PCM Coding The DPCM function @dfn{dpcmenco}, @dfn{dpcmdeco} and @dfn{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. @node Arithmetic Coding, Dynamic Range Compression, PCM Coding, Source Coding @section Arithmetic Coding The arithmetic coding functions @dfn{arithenco} and @dfn{arithdeco} are not yet implemented. @node Dynamic Range Compression, , Arithmetic Coding, Source Coding @section Dynamic Range Compression The final source coding function is @dfn{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. @dfn{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 @dfn{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. @menu * Data Formats:: * Binary Block Codes:: * BCH Codes:: * Reed-Solomon Codes:: @end menu @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 @iftex @tex $ {\bf C} = {\bf T} {\bf G}$. @end tex @end iftex @ifinfo @var{C} = @var{T} * @var{G}. @end ifinfo @item Parity Check Matrix A '@var{n}-@var{k}'-by-@var{n} matrix @var{H} to check the parity of the received symbols. If @iftex @tex $ {\bf H} {\bf R} = {\bf S} \ne 0$, @end tex @end iftex @ifinfo @var{H} * @var{R} = @var{S} != 0, @end ifinfo 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 @iftex @tex $ {\bf R} = ( {\bf T} + {\bf E} )\ mod\ 2$, @end tex @end iftex @ifinfo @var{R} = mod(@var{T} + @var{E}, 2), @end ifinfo then the error vector @var{E} is @iftex @tex ${\bf ST}({\bf S})$. @end tex @end iftex @ifinfo @var{ST}(@var{S}). @end ifinfo @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 @iftex @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 @end iftex @ifinfo @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 ifinfo or @iftex @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 @end iftex @ifinfo @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 ifinfo 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 @dfn{hammgen}, @dfn{cyclgen} and @dfn{cyclpoly}. The function @dfn{hammgen} creates parity check and generator matrices for Hamming codes, while @dfn{cyclpoly} and @dfn{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 @dfn{syndtable}, in the following manner @example octave:1> [par, gen] = hammgen(3); octave:2> st = syndtable(par); @end example There exists two auxiliary functions @dfn{gen2par} and @dfn{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 @dfn{encode} and @dfn{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 @dfn{encode} and @dfn{decode} must be told the format of the messages in the following manner. @example 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"); @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 @dfn{encode} and @dfn{decode} by default accept binary coded messages. Therefore 'hamming' is equivalent to 'hamming/binary'. Except for the BCH codes, the function @dfn{encode} and @dfn{decode} internally create the generator, parity check and syndrome table matrices. Therefore if repeated calls to @dfn{encode} and @dfn{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 !!! end @end example In the case of the BCH codes the low-level functions described in the next section are used directly by the @dfn{encode} and @dfn{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, @dfn{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 @ifinfo Called without arguments, @dfn{bchpoly} returns a table of valid BCH error correcting codes and their error-correction capability. @end ifinfo The first returned column of @dfn{bchpoly} is the codeword length, the second the message length and the third the error correction capability of the code. Called with one argument, @dfn{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, @dfn{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 @iftex @tex $$ 1 + x^4 + x^6 + x^7 + x^8 $$ @end tex @end iftex @ifinfo 1 + @var{x} ^ 4 + @var{x} ^ 6 + @var{x} ^ 7 + @var{x} ^ 8 @end ifinfo 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 @dfn{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 @dfn{encode} and @dfn{decode} functions previously discussed, the user can directly use the low-level BCH functions @dfn{bchenco} and @dfn{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(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); @end example @node Reed-Solomon Codes, , BCH Codes, Block Coding @section Reed-Solomon Codes @menu * Representation of Reed-Solomon Messages:: * Creating and Decoding Messages:: * Shortened Reed-Solomon Codes:: @end menu @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 @dfn{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 @dfn{rsenc} and @dfn{rsdec}. By default the parity symbols are placed after the message. Valid generator polynomials can be constructed with the @dfn{rsgenpoly} function. The roots of the generator polynomial are then defined by @iftex @tex $$ g = (x - A^{bs}) (x - A^{(b+1)s}) \cdots (x - A ^{(b+2t-1)s}). $$ @end tex @end iftex @ifinfo @example @var{g} = (@var{x} - A^(@var{b}*@var{s})) * (@var{x} - A^((@var{b}+1)*@var{s})) * ... * (@var{x} - A^((@var{b}+2*@var{t}-1)*@var{s})). @end example @end ifinfo 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 @dfn{rsgenpoly} and can be passed to both @dfn{rsenc} and @dfn{rsdec}. It is also possible to pass the @var{b} and @var{s} values directly to @dfn{rsenc} and @dfn{rsdec}. In the case of @dfn{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: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 @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 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 @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, qaskdeco and qaskenco. @node Special Filters, Galois Fields, Modulations, Top @chapter Special Filters To be written. @node Galois Fields, Function Reference, Special Filters, Top @chapter Galois Fields @menu * Galois Field Basics:: * Manipulating Galois Fields:: @end menu @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 performed and so I am most interested in comments about compatibility at the e-mail address dbateman@@free.fr. @menu * Creating Galois Fields:: * Primitive Polynomials:: * Accessing Internal Fields:: * Function Overloading:: * Known Problems:: @end menu @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 *** local user variables: prot type rows cols name ==== ==== ==== ==== ==== rwd matrix 1 8 a rwd galois 1 8 b @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: primitive polynomial (21) of Galois Field must be irreducible error: unable to initialize Galois Field error: evaluating assignment expression near line 3, column 3 @end example The function @dfn{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, @dfn{isprimitive} and @dfn{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 @dfn{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 @dfn{primpoly} returns an empty vector. That is @example octave:1> primpoly(5,2) primpoly: No primitive polynomial satisfies the given constraints ans = [](1x0) @end example As can be seen above, @dfn{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 @dfn{isprimitive} can be used in a similar manner to @dfn{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). @node Accessing Internal Fields, Function Overloading, Primitive Polynomials, Galois Field Basics @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 *** local user variables: prot type rows cols name ==== ==== ==== ==== ==== rwd matrix 1 8 a rwd galois 1 8 b rwd matrix 1 8 c @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 @node Function Overloading, Known Problems, Accessing Internal Fields, Galois Field Basics @subsection Function Overloading An important consideration in the use of the Galois Field package is that many of the internal functions of Octave, such as @dfn{roots}, can not accept Galois Fields as an input. This package therefore uses the @dfn{dispatch} function of Octave-Forge to @emph{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 @multitable @columnfractions .11 .036 .17 .122 .036 .17 .122 .036 .17 @item @code{conv} @tab - @tab @code{gconv}, @tab @code{convmtx} @tab - @tab @code{gconvmtx}, @tab @code{diag} @tab - @tab @code{gdiag}, @item @code{deconv} @tab - @tab @code{gdeconv}, @tab @code{det} @tab - @tab @code{gdet}, @tab @code{exp} @tab - @tab @code{gexp}, @item @code{filter} @tab - @tab @code{gfilter}, @tab @code{inv} @tab - @tab @code{ginv}, @tab @code{log} @tab - @tab @code{glog}, @item @code{lu} @tab - @tab @code{glu}, @tab @code{prod} @tab - @tab @code{gprod}, @tab @code{reshape} @tab - @tab @code{greshape}, @item @code{rank} @tab - @tab @code{grank}, @tab @code{roots} @tab - @tab @code{groots}, @tab @code{sum} @tab - @tab @code{gsum}, @item @code{sumsq} @tab - @tab @code{gsumsq}. @end multitable The version of the function that is chosen is determined by the first argument of the function. So, considering the @dfn{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 @dfn{size} and @dfn{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. @node Known Problems, , Function Overloading, Galois Field Basics @subsection Known Problems Before reporting a bug compare it to this list of known problems @table @asis @item Concatenation For versions of Octave prior to 2.1.58, the concatenation of Galois arrays returns a Matrix type. That is @code{[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 @code{gf([gf([1, 0],m) gf(1,m)],m)}. Since Octave version 2.1.58, @code{[gf([1, 0],m) gf(1, m)]} returns another Galois array as expected. @item 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 @example octave:2> save a.mat a @end example where @var{a} is a Galois variable. To reload the variable within octave, the Galois type must be installed prior to a call to @dfn{load}. That is @example octave:1> dummy = gf(1); octave:2> load a.mat @end example 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 @example octave:2> x = a.x; m = a.m; p = a.prim_poly; octave:3> save a.mat x m p; @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 however printed 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 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. @end table @node Manipulating Galois Fields, , Galois Field Basics, Galois Fields @section Manipulating Galois Fields @menu * Expressions manipulation and assignment:: * Unary operations:: * Arithmetic operations:: * Comparison operations:: * Polynomial manipulations:: * Linear Algebra:: * Signal Processing:: @end menu @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. 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. @example 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 @end example and for previous versions of Octave @example 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 @end example 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 @example octave:1> a = gf([0:7],3); octave:2> b = gf([a, a], a.m, a.prim_poly); @end example which explicitly reconverts @var{b} to the correct Galois Field. 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 @dfn{any} and @dfn{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 @iftex @tex $$ poly = a x^3 + a^2 x^2 + (a^2 + 1) x + 1 $$ @end tex @end iftex @ifinfo @example poly = @var{a} * x^3 + @var{a}^2 * x^2 + (@var{a}^2 + 1) * x + 1 @end example @end ifinfo 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 @dfn{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 @dfn{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 @dfn{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 @dfn{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 @dfn{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 [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 @} @end example 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 @dfn{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 the 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; end; 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 @dfn{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 @dfn{conv} can be used, or alternatively the function @dfn{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 @dfn{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, @dfn{fft}, @dfn{ifft} and @dfn{dftmtx}. The first two functions use the third to perform their work. Give an element @var{a} of the Galois field GF(2^M), @dfn{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 @dfn{fft} function @dfn{dftmtx} is called with the 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 @dfn{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. @node Function Reference, , Galois Fields, Top @chapter Function Reference @REFERENCE_SECTION(Communications) @bye