X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fcommunications-1.1.1%2Fcomms.m;fp=octave_packages%2Fcommunications-1.1.1%2Fcomms.m;h=8f8ecde21cc95d49ad4cebb775b3c5fb629b098f;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/communications-1.1.1/comms.m b/octave_packages/communications-1.1.1/comms.m new file mode 100644 index 0000000..8f8ecde --- /dev/null +++ b/octave_packages/communications-1.1.1/comms.m @@ -0,0 +1,665 @@ +## 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 . + +## -*- 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