--- /dev/null
+## Copyright (C) 2003 David Bateman
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} comms ('help')
+## @deftypefnx {Function File} {} comms ('info')
+## @deftypefnx {Function File} {} comms ('info', @var{mod})
+## @deftypefnx {Function File} {} comms ('test')
+## @deftypefnx {Function File} {} comms ('test', @var{mod})
+##
+## Manual and test code for the Octave Communications toolbox. There are
+## 5 possible ways to call this function.
+##
+## @table @code
+## @item comms ('help')
+## Display this help message. Called with no arguments, this function also
+## displays this help message
+## @item comms ('info')
+## Open the Commumications toolbox manual
+## @item comms ('info', @var{mod})
+## Open the Commumications toolbox manual at the section specified by
+## @var{mod}
+## @item comms ('test')
+## Run all of the test code for the Communications toolbox.
+## @item comms ('test', @var{mod})
+## Run only the test code for the Communications toolbox in the module
+## @var{mod}.
+## @end table
+##
+## Valid values for the varibale @var{mod} are
+##
+## @table @asis
+## @item 'all'
+## All of the toolbox
+## @item 'random'
+## The random signal generation and analysis package
+## @item 'source'
+## The source coding functions of the package
+## @item 'block'
+## The block coding functions
+## @item 'convol'
+## The convolution coding package
+## @item 'modulation'
+## The modulation package
+## @item 'special'
+## The special filter functions
+## @item 'galois'
+## The Galois fields package
+## @end table
+##
+## Please note that this function file should be used as an example of the
+## use of this toolbox.
+## @end deftypefn
+
+function retval = comms(typ, tests)
+
+ if (nargin < 1)
+ help ("comms");
+ elseif (nargin < 2)
+ tests = 'all';
+ endif
+
+
+ if strcmp(tests,"all")
+ nodename = "Top";
+ elseif strcmp(tests,"random")
+ nodename = "Random Signals";
+ elseif strcmp(tests,"source")
+ nodename = "Source Coding";
+ elseif strcmp(tests,"block")
+ nodename = "Block Coding";
+ elseif strcmp(tests,"convol")
+ nodename = "Convolutional Coding";
+ elseif strcmp(tests,"modulation")
+ nodename = "Modulations";
+ elseif strcmp(tests,"special")
+ nodename = "Special Fields";
+ elseif strcmp(tests,"galois")
+ nodename = "Galois Fields";
+ else
+ error ("comms: unrecognized package");
+ endif
+
+ if (strcmp(typ,"help"))
+ help ("comms");
+ elseif (strcmp(typ,"info"))
+ infopaths = ["."];
+ if (!isempty(char(strsplit (path, ":"))))
+ infopaths =[infopaths; char(strsplit (path, ":"))];
+ endif
+ if (!isempty(char(strsplit (DEFAULT_LOADPATH, ":"))))
+ infopaths =[infopaths; char(strsplit (DEFAULT_LOADPATH, ":"))];
+ endif
+ for i=1:size(infopaths,1)
+ infopath = deblank(infopaths(i,:));
+ len = length(infopath);
+ if (len)
+ if (len > 1 && strcmp(infopath([len-1, len]),"//"))
+ [status, showfile] = system(["find '", infopath(1:len-1), ...
+ "' -name ", infofile]);
+ else
+ [status, showfile] = system(["find '", infopath, "' -name ", ...
+ infofile, " -maxdepth 1"]);
+ endif
+ if (length(showfile))
+ break;
+ endif
+ endif
+ end
+ if (!exist("showfile") || !length(showfile))
+ error("comms: info file not found");
+ endif
+ if (showfile(length(showfile)) == "\n")
+ showfile = showfile(1:length(showfile)-1);
+ endif
+
+ if (exist("INFO_PROGAM"))
+ [testret, testout] = system(["'", INFO_PROGRAM, "' --version"]);
+ if (testret)
+ error("comms: info command not found");
+ else
+ system(["'", INFO_PROGRAM, "' --file '", showfile, "' --node '", ...
+ nodename, "'"]);
+ endif
+ else
+ [testret, testout] = system("info --version");
+ if (testret)
+ error("comms: info command not found");
+ else
+ system(["info --file '", showfile, "' --node '", nodename, "'"]);
+ endif
+ endif
+ elseif (strcmp(typ,"test"))
+ pso = page_screen_output();
+ unwind_protect
+ page_screen_output(0);
+
+ if (strcmp(tests,"random") || strcmp(tests,"all"))
+ fprintf("\n<< Random Signals Package >>\n");
+ fprintf(" Signal Creation: ");
+ n = 10;
+ m = 32;
+ x = randint(n,n,m);
+ if (any(size(x) != [10, 10]) || (max(x(:)) >= m) || (min(x(:) < 0)))
+ error ("FAILED");
+ endif
+ x = randsrc(n,n,[1, 1i, -1, -1i,]);
+ if (any(size(x) != [10, 10]) || ...
+ !all(all((x == 1) | (x == 1i) | (x == -1) | (x == -1i))))
+ error ("FAILED");
+ endif
+ x = randerr(n,n);
+ if (any(size(x) != [10, 10]) || any(sum(x') != ones(1,n)))
+ error ("FAILED");
+ endif
+
+ nse_30dBm_1Ohm = wgn(10000,1,30,1,"dBm");
+ nse_0dBW_1Ohm = wgn(10000,1,0,1,"dBW");
+ nse_1W_1Ohm = wgn(10000,1,1,1,"linear");
+ ## Standard deviations should be about 1... If it is greater than
+ ## some value flag an error
+ dev = [std(nse_30dBm_1Ohm), std(nse_0dBW_1Ohm), std(nse_1W_1Ohm)];
+ if (any(dev > 1.5))
+ error ("FAILED");
+ endif
+
+ x = [0:0.1:2*pi];
+ y = sin(x);
+ noisy = awgn(y, 10, "dB", "measured");
+ if (any(size(y) != size(noisy)))
+ error ("FAILED");
+ endif
+ ## This is a pretty arbitrary test, but should pick up gross errors
+ if (any(abs(y-noisy) > 1))
+ error ("FAILED");
+ endif
+ fprintf("PASSED\n");
+
+ fprintf(" Signal Analysis: ");
+ ## Protect!! Since bitxor might not be installed
+ try
+ n = 10;
+ m = 8;
+ msg = randint(n,n,2^m);
+ noisy = bitxor(msg,diag(3*ones(1,n)));
+ [berr, brate] = biterr(msg, noisy, m);
+ if ((berr != 2*n) || (brate != 2/(n*m)))
+ error ("FAILED");
+ endif
+ [serr, srate] = symerr(msg, noisy);
+ if ((serr != n) || (srate != 1/n))
+ error ("FAILED");
+ endif
+ catch
+ end
+ ## Can not easily test eyediagram, scatterplot!!
+ fprintf("PASSED\n");
+ endif
+ if (strcmp(tests,"source") || strcmp(tests,"all"))
+ fprintf("\n<< Source Coding Package >>\n");
+ fprintf(" PCM Functions: ");
+ fprintf("Not tested\n");
+ fprintf(" Quantization Functions: ");
+ x = [0:0.1:2*pi];
+ y = sin(x);
+ [tab, cod] = lloyds(y, 16);
+ [i, q, d] = quantiz(y, tab, cod);
+ if (abs(d) > 0.1)
+ error ("FAILED");
+ endif
+
+ mu = 0.1;
+ V = 1;
+ x = sin([0:0.1:2*pi]);
+ y = compand(x, mu, V, "mu/compressor");
+ z = compand(x, mu, V, "mu/expander");
+ ## Again this is a pretty arbitrary test
+ if (max(abs(x-z)) > 0.1)
+ error ("FAILED");
+ endif
+ fprintf("PASSED\n");
+ endif
+ if (strcmp(tests,"block") || strcmp(tests,"all"))
+ fprintf("\n<< Block Coding Package >>\n");
+ fprintf(" Cyclic Coding: ");
+ nsym = 100;
+ m = 4;
+ n = 2^m-1; # [15,11] Hamming code
+ k = n - m;
+ p = cyclpoly(n,k);
+ if (bi2de(p) != primpoly(m,"nodisplay"))
+ error("FAILED");
+ endif
+ [par, gen] = cyclgen(n,p);
+ if (any(any(gen2par(par) != gen)))
+ error("FAILED");
+ endif
+ if (gfweight(gen) != 3)
+ error("FAILED");
+ endif
+
+ msg = randint(nsym,k);
+ code = encode(msg,n,k,"cyclic");
+ noisy = mod(code + randerr(nsym,n), 2);
+ dec = decode(noisy,n,k,"cyclic");
+ if (any(any(dec != msg)))
+ error("FAILED");
+ endif
+ try # Protect! If bitshift isn't install!!
+ msg = randint(nsym,1,n);
+ code = encode(msg,n,k,"cyclic/decimal");
+ noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
+ dec = decode(noisy,n,k,"cyclic/decimal");
+ if (any(dec != msg))
+ error("FAILED");
+ endif
+ catch
+ end
+ fprintf("PASSED\n");
+
+ fprintf(" Hamming Coding: ");
+ nsym = 100;
+ m = 4;
+ [par, gen, n, k] = hammgen (m);
+ if (any(any(gen2par(par) != gen)))
+ error("FAILED");
+ endif
+ if (gfweight(gen) != 3)
+ error("FAILED");
+ endif
+ msg = randint(nsym,k);
+ code = encode(msg,n,k,"hamming");
+ noisy = mod(code + randerr(nsym,n), 2);
+ dec = decode(noisy,n,k,"hamming");
+ if (any(any(dec != msg)))
+ error("FAILED");
+ endif
+ try # Protect! If bitshift isn't install!!
+ msg = randint(nsym,1,n);
+ code = encode(msg,n,k,"hamming/decimal");
+ noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
+ dec = decode(noisy,n,k,"hamming/decimal");
+ if (any(dec != msg))
+ error("FAILED");
+ endif
+ catch
+ end
+ fprintf("PASSED\n");
+
+ fprintf(" BCH Coding: ");
+ ## Setup
+ m = 5;
+ nsym = 100;
+ p = bchpoly(2^m - 1);
+ ## Pick a BCH code from the list at random
+ l = ceil(size(p,1) * rand(1,1));
+ n = p(l,1);
+ k = p(l,2);
+ t = p(l,3);
+ ## Symbols represented by rows of binary matrix
+ msg = randint(nsym,k);
+ code = encode(msg,n,k,"bch");
+ noisy = mod(code + randerr(nsym,n), 2);
+ dec = decode(noisy,n,k,"bch");
+ if (any(any(dec != msg)))
+ error("FAILED");
+ endif
+ try # Protect! If bitshift isn't install!!
+ msg = randint(nsym,1,n);
+ code = encode(msg,n,k,"bch/decimal");
+ noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
+ dec = decode(noisy,n,k,"bch/decimal");
+ if (any(dec != msg))
+ error("FAILED");
+ endif
+ catch
+ end
+ fprintf("PASSED\n");
+
+ fprintf(" Reed-Solomon Coding: ");
+ ## Test for a CCSDS like coder, but without dual-basis translation
+ mrs = 8;
+ nrs = 2^mrs -1;
+ krs = nrs - 32;
+ prs = 391;
+ fcr = 112;
+ step = 11;
+
+ ## CCSDS generator polynomial
+ ggrs = rsgenpoly(nrs, krs, prs, fcr, step);
+
+ ## Code two blocks
+ msg = gf(floor(2^mrs*rand(2,krs)),mrs,prs);
+ cde = rsenc(msg,nrs,krs,ggrs);
+
+ ## Introduce errors
+ noisy = cde + [ 255,0,255,0,255,zeros(1,250); ...
+ 0,255,0,255,zeros(1,251)];
+
+ ## Decode (better to pass fcr and step rather than gg for speed)
+ dec = rsdec(noisy,nrs,krs,fcr,step);
+
+ if (any(dec != msg))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ endif
+ if (strcmp(tests,"convol") || strcmp(tests,"all"))
+ fprintf("\n<< Convolutional Coding Package >>\n");
+ fprintf(" Utility functions: ");
+ ## create a trellis, use poly2trellis and test with istrellis
+ fprintf("Not tested\n");
+ fprintf(" Coding: ");
+ ## use convenc, punturing??
+ fprintf("Not tested\n");
+ fprintf(" Viterbi: ");
+ ## use vitdec
+ fprintf("Not tested\n");
+ endif
+ if (strcmp(tests,"modulation") || strcmp(tests,"all"))
+ fprintf("\n<< Modulation Package >>\n");
+ fprintf(" Analog Modulation: ");
+ Fs = 100;
+ t = [0:1/Fs:2];
+ x = sin(2*pi*t);
+ xq = x + 1i * cos(2*pi*t);
+ ## Can not test FM as it doesn't work !!!
+ if ((max(abs(x - ademodce(amodce(x,Fs,"pm"),Fs,"pm"))) > 0.001) || ...
+ (max(abs(x - ademodce(amodce(x,Fs,"am"),Fs,"am"))) > 0.001) || ...
+ (max(abs(xq - ademodce(amodce(xq,Fs,"qam"),Fs,"qam/cmplx"))) > 0.001))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ fprintf(" Digital Mapping: ");
+ m = 32;
+ n = 100;
+ x = randint(n,n,32);
+ if ((x != demodmap(modmap(x,1,1,"ask",m),1,1,"ask",m)) || ...
+ (x != demodmap(modmap(x,1,1,"fsk",m),1,1,"fsk",m)) || ...
+ (x != demodmap(modmap(x,1,1,"msk"),1,1,"msk")) || ...
+ (x != demodmap(modmap(x,1,1,"psk",m),1,1,"psk",m)) || ...
+ (x != demodmap(modmap(x,1,1,"qask",m),1,1,"qask",m)) || ...
+ (x != demodmap(modmap(x,1,1,"qask/cir", ...
+ [floor(m/2), m - floor(m/2)]),1,1,"qask/cir", ...
+ [floor(m/2), m - floor(m/2)])))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ fprintf(" Digital Modulation: ");
+ fprintf("Not tested\n");
+ endif
+ if (strcmp(tests,"special") || strcmp(tests,"all"))
+ fprintf("\n<< Special Filters Package >>\n");
+ fprintf(" Hankel/Hilbert: ");
+ ## use hank2sys, hilbiir
+ fprintf("Not tested\n");
+ fprintf(" Raised Cosine: ");
+ ## use rcosflt, rcosiir rcosine, rcosfir
+ fprintf("Not tested\n");
+ endif
+ if (strcmp(tests,"galois") || strcmp(tests,"all"))
+ fprintf("\n<< Galois Fields Package >>\n");
+ ## Testing of the Galois Fields package
+ m = 3; ## must be greater than 2
+
+ fprintf(" Find primitive polynomials: ");
+ prims = primpoly(m,"all","nodisplay");
+ for i=2^m:2^(m+1)-1
+ if (find(prims == i))
+ if (!isprimitive(i))
+ error("Error in primitive polynomials");
+ endif
+ else
+ if (isprimitive(i))
+ error("Error in primitive polynomials");
+ endif
+ endif
+ end
+ fprintf("PASSED\n");
+ fprintf(" Create Galois variables: ");
+ n = 2^m-1;
+ gempty = gf([],m);
+ gzero = gf(0,m);
+ gone = gf(1,m);
+ gmax = gf(n,m);
+ grow = gf(0:n,m);
+ gcol = gf([0:n]',m);
+ matlen = ceil(sqrt(2^m));
+ gmat = gf(reshape(mod([0:matlen*matlen-1],2^m),matlen,matlen),m);
+ fprintf("PASSED\n");
+ fprintf(" Access Galois structures: ");
+ if (gcol.m != m || gcol.prim_poly != primpoly(m,"min", ...
+ "nodisplay"))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ fprintf(" Miscellaneous functions: ");
+ if (size(gmat) != [matlen, matlen])
+ error("FAILED");
+ endif
+ if (length(grow) != 2^m)
+ error("FAILED");
+ endif
+ if (!any(grow) || all(grow) || any(gzero) || !all(gone))
+ error("FAILED");
+ endif
+ if (isempty(gone) || !isempty(gempty))
+ error("FAILED");
+ endif
+ tmp = diag(grow);
+ if (size(tmp,1) != 2^m || size(tmp,2) != 2^m)
+ error("FAILED");
+ endif
+ for i=1:2^m
+ for j=1:2^m
+ if ((i == j) && (tmp(i,j) != grow(i)))
+ error("FAILED");
+ elseif ((i != j) && (tmp(i,j) != 0))
+ error("FAILED");
+ endif
+ end
+ end
+ tmp = diag(gmat);
+ if (length(tmp) != matlen)
+ error("FAILED");
+ endif
+ for i=1:matlen
+ if (gmat(i,i) != tmp(i))
+ error("FAILED");
+ endif
+ end
+ tmp = reshape(gmat,prod(size(gmat)),1);
+ if (length(tmp) != prod(size(gmat)))
+ error("FAILED");
+ endif
+ if (exp(log(gf([1:n],m))) != [1:n])
+ error("FAILED");
+ endif
+ tmp = sqrt(gmat);
+ if (tmp .* tmp != gmat)
+ error("FAILED");
+ endif
+
+ fprintf("PASSED\n");
+ fprintf(" Unary operators: ");
+ tmp = - grow;
+ if (tmp != grow)
+ error("FAILED");
+ endif
+ tmp = !grow;
+ if (tmp(1) != 1)
+ error("FAILED");
+ endif
+ if (any(tmp(2:length(tmp))))
+ error("FAILED");
+ endif
+ tmp = gmat';
+ for i=1:size(gmat,1)
+ for j=1:size(gmat,2)
+ if (gmat(i,j) != tmp(j,i))
+ error("FAILED");
+ endif
+ end
+ end
+ fprintf("PASSED\n");
+ fprintf(" Arithmetic operators: ");
+ if (any(gmat + gmat))
+ error("FAILED");
+ endif
+ multbl = gcol * grow;
+ elsqr = gcol .* gcol;
+ elsqr2 = gcol .^ gf(2,m);
+ for i=1:length(elsqr)
+ if (elsqr(i) != multbl(i,i))
+ error("FAILED");
+ endif
+ end
+ for i=1:length(elsqr)
+ if (elsqr(i) != elsqr2(i))
+ error("FAILED");
+ endif
+ end
+ tmp = grow(2:length(grow)) ./ gcol(2:length(gcol))';
+ if (length(tmp) != n || any(tmp != ones(1,length(grow)-1)))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ fprintf(" Logical operators: ");
+ if (grow(1) != gzero || grow(2) != gone || grow(2^m) != n)
+ error("FAILED");
+ endif
+ if (!(grow(1) == gzero) || any(grow != gcol'))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+
+ fprintf(" Polynomial manipulation: ");
+ poly1 = gf([2,4,5,1],3);
+ poly2 = gf([1,2],3);
+ sumpoly = poly1 + [0,0,poly2]; ## Already test "+"
+ mulpoly = conv(poly1, poly2); ## multiplication
+ poly3 = [poly,remd] = deconv(mulpoly, poly2);
+ if (!isequal(poly1,poly3))
+ error("FAILED");
+ endif
+ if (any(remd))
+ error("FAILED");
+ endif
+ x0 = gf([0,1,2,3],3);
+ y0 = polyval(poly1, x0);
+ alph = gf(2,3);
+ y1 = alph * x0.^3 + alph.^2 * x0.^2 + (alph.^2+1) *x0 + 1;
+ if (!isequal(y0,y1))
+ error("FAILED");
+ endif
+ roots1 = roots(poly1);
+ ck1 = polyval(poly1, roots1);
+ if (any(ck1))
+ error("FAILED");
+ endif
+ b = minpol(alph);
+ bpoly = bi2de(b.x,"left-msb");
+ if (bpoly != alph.prim_poly)
+ error("FAILED");
+ endif
+ c = cosets(3);
+ c2 = c{2};
+ mpol = minpol(c2(1));
+ for i=2:length(c2)
+ if (mpol != minpol(c2(i)))
+ error("FAILED");
+ endif
+ end
+ fprintf("PASSED\n");
+
+ fprintf(" Linear Algebra: ");
+ [l,u,p] = lu(gmat);
+ if (any(l*u-p*gmat))
+ error("FAILED");
+ endif
+ g1 = inv(gmat);
+ g2 = gmat ^ -1;
+ if (any(g1*gmat != eye(matlen)))
+ error("FAILED");
+ endif
+ if (any(g1 != g2))
+ error("FAILED");
+ endif
+ matdet = 0;
+ while (!matdet)
+ granmat = gf(floor(2^m*rand(matlen)),m);
+ matdet = det(granmat);
+ endwhile
+ matrank = rank(granmat);
+ smallcol = gf([0:matlen-1],m)';
+ sol1 = granmat \ smallcol;
+ sol2 = smallcol' / granmat;
+ if (any(granmat * sol1 != smallcol))
+ error("FAILED");
+ endif
+ if (any(sol2 * granmat != smallcol'))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+ fprintf(" Signal Processing functions: ");
+ b = gf([1,0,0,1,0,1,0,1],m);
+ a = gf([1,0,1,1],m);
+ x = gf([1,zeros(1,99)],m);
+ y0 = filter(b, a, x);
+ y1 = conv(grow+1, grow);
+ y2 = grow * convmtx(grow+1, length(grow));
+ if (any(y1 != y2))
+ error("FAILED");
+ endif
+ [y3,remd] = deconv(y2, grow+1);
+ if (any(y3 != grow))
+ error("FAILED");
+ endif
+ if (any(remd))
+ error("FAILED");
+ endif
+ alph = gf(2,m);
+ x = gf(floor(2^m*rand(n,1)),m);
+ fm = dftmtx(alph);
+ ifm = dftmtx(1/alph);
+ y0 = fft(x);
+ y1 = fm * x;
+ if (any(y0 != y1))
+ error("FAILED");
+ endif
+ z0 = ifft(y0);
+ z1 = ifm * y1;
+ if (any(z0 != x))
+ error("FAILED");
+ endif
+ if (any(z1 != x))
+ error("FAILED");
+ endif
+ fprintf("PASSED\n");
+
+ endif
+ fprintf("\n");
+ unwind_protect_cleanup
+ page_screen_output(pso);
+ end_unwind_protect
+ else
+ usage("comms: Unknown argument");
+ endif
+endfunction
+
+%!test
+%! try comms("test");
+%! catch disp(lasterr()); end
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff