octave_packages/communications-1.1.1/comms.m

   1 ## Copyright (C) 2003 David Bateman
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {} comms ('help')
  18 ## @deftypefnx {Function File} {} comms ('info')
  19 ## @deftypefnx {Function File} {} comms ('info', @var{mod})
  20 ## @deftypefnx {Function File} {} comms ('test')
  21 ## @deftypefnx {Function File} {} comms ('test', @var{mod})
  22 ##
  23 ## Manual and test code for the Octave Communications toolbox. There are
  24 ## 5 possible ways to call this function.
  25 ##
  26 ## @table @code
  27 ## @item comms ('help')
  28 ## Display this help message. Called with no arguments, this function also
  29 ## displays this help message
  30 ## @item comms ('info')
  31 ## Open the Commumications toolbox manual
  32 ## @item comms ('info', @var{mod})
  33 ## Open the Commumications toolbox manual at the section specified by
  34 ## @var{mod}
  35 ## @item comms ('test')
  36 ## Run all of the test code for the Communications toolbox.
  37 ## @item comms ('test', @var{mod})
  38 ## Run only the test code for the Communications toolbox in the module
  39 ## @var{mod}.
  40 ## @end table
  41 ##
  42 ## Valid values for the varibale @var{mod} are
  43 ##
  44 ## @table @asis
  45 ## @item 'all'
  46 ## All of the toolbox
  47 ## @item 'random'
  48 ## The random signal generation and analysis package
  49 ## @item 'source'
  50 ## The source coding functions of the package
  51 ## @item 'block'
  52 ## The block coding functions
  53 ## @item 'convol'
  54 ## The convolution coding package
  55 ## @item 'modulation'
  56 ## The modulation package
  57 ## @item 'special'
  58 ## The special filter functions
  59 ## @item 'galois'
  60 ## The Galois fields package
  61 ## @end table
  62 ##
  63 ## Please note that this function file should be used as an example of the
  64 ## use of this toolbox.
  65 ## @end deftypefn
  66
  67 function retval = comms(typ, tests)
  68
  69   if (nargin < 1)
  70     help ("comms");
  71   elseif (nargin < 2)
  72     tests = 'all';
  73   endif
  74
  75
  76   if strcmp(tests,"all")
  77     nodename = "Top";
  78   elseif strcmp(tests,"random")
  79     nodename = "Random Signals";
  80   elseif  strcmp(tests,"source")
  81     nodename = "Source Coding";
  82   elseif  strcmp(tests,"block")
  83     nodename = "Block Coding";
  84   elseif  strcmp(tests,"convol")
  85     nodename = "Convolutional Coding";
  86   elseif  strcmp(tests,"modulation")
  87     nodename = "Modulations";
  88   elseif  strcmp(tests,"special")
  89     nodename = "Special Fields";
  90   elseif  strcmp(tests,"galois")
  91     nodename = "Galois Fields";
  92   else
  93     error ("comms: unrecognized package");
  94   endif
  95
  96   if (strcmp(typ,"help"))
  97     help ("comms");
  98   elseif (strcmp(typ,"info"))
  99     infopaths = ["."];
 100     if (!isempty(char(strsplit (path, ":"))))
 101       infopaths =[infopaths; char(strsplit (path, ":"))];
 102     endif
 103     if (!isempty(char(strsplit (DEFAULT_LOADPATH, ":"))))
 104       infopaths =[infopaths; char(strsplit (DEFAULT_LOADPATH, ":"))];
 105     endif
 106     for i=1:size(infopaths,1)
 107       infopath = deblank(infopaths(i,:));
 108       len = length(infopath);
 109       if (len)
 110         if (len > 1 && strcmp(infopath([len-1, len]),"//"))
 111           [status, showfile] = system(["find '", infopath(1:len-1), ...
 112                               "' -name ", infofile]);
 113         else
 114           [status, showfile] = system(["find '", infopath, "' -name ", ...
 115                               infofile, " -maxdepth 1"]);
 116         endif
 117         if (length(showfile))
 118           break;
 119         endif
 120       endif
 121     end
 122     if (!exist("showfile") || !length(showfile))
 123       error("comms: info file not found");
 124     endif
 125     if (showfile(length(showfile)) == "\n")
 126       showfile = showfile(1:length(showfile)-1);
 127     endif
 128
 129     if (exist("INFO_PROGAM"))
 130       [testret, testout] = system(["'", INFO_PROGRAM, "' --version"]);
 131       if (testret)
 132         error("comms: info command not found");
 133       else
 134         system(["'", INFO_PROGRAM, "' --file '", showfile, "' --node '", ...
 135                 nodename, "'"]);
 136       endif
 137     else
 138       [testret, testout] = system("info --version");
 139       if (testret)
 140         error("comms: info command not found");
 141       else
 142         system(["info --file '", showfile, "' --node '", nodename, "'"]);
 143       endif
 144     endif
 145   elseif (strcmp(typ,"test"))
 146     pso = page_screen_output();
 147     unwind_protect
 148       page_screen_output(0);
 149
 150       if (strcmp(tests,"random") || strcmp(tests,"all"))
 151         fprintf("\n<< Random Signals Package >>\n");
 152         fprintf("  Signal Creation:                          ");
 153         n = 10;
 154         m = 32;
 155         x = randint(n,n,m);
 156         if (any(size(x) != [10, 10]) || (max(x(:)) >= m) || (min(x(:) < 0)))
 157           error ("FAILED");
 158         endif
 159         x = randsrc(n,n,[1, 1i, -1, -1i,]);
 160         if (any(size(x) != [10, 10]) || ...
 161             !all(all((x == 1) | (x == 1i) | (x == -1) | (x == -1i))))
 162           error ("FAILED");
 163         endif
 164         x = randerr(n,n);
 165         if (any(size(x) != [10, 10]) || any(sum(x') != ones(1,n)))
 166           error ("FAILED");
 167         endif
 168
 169         nse_30dBm_1Ohm = wgn(10000,1,30,1,"dBm");
 170         nse_0dBW_1Ohm = wgn(10000,1,0,1,"dBW");
 171         nse_1W_1Ohm = wgn(10000,1,1,1,"linear");
 172         ## Standard deviations should be about 1... If it is greater than
 173         ## some value flag an error
 174         dev = [std(nse_30dBm_1Ohm), std(nse_0dBW_1Ohm), std(nse_1W_1Ohm)];
 175         if (any(dev > 1.5))
 176           error ("FAILED");
 177         endif
 178
 179         x = [0:0.1:2*pi];
 180         y = sin(x);
 181         noisy = awgn(y, 10, "dB", "measured");
 182         if (any(size(y) != size(noisy)))
 183           error ("FAILED");
 184         endif
 185         ## This is a pretty arbitrary test, but should pick up gross errors
 186         if (any(abs(y-noisy) > 1))
 187           error ("FAILED");
 188         endif
 189         fprintf("PASSED\n");
 190
 191         fprintf("  Signal Analysis:                          ");
 192         ## Protect!! Since bitxor might not be installed
 193         try
 194           n = 10;
 195           m = 8;
 196           msg = randint(n,n,2^m);
 197           noisy = bitxor(msg,diag(3*ones(1,n)));
 198           [berr, brate] = biterr(msg, noisy, m);
 199           if ((berr != 2*n) || (brate != 2/(n*m)))
 200             error ("FAILED");
 201           endif
 202           [serr, srate] = symerr(msg, noisy);
 203           if ((serr != n) || (srate != 1/n))
 204             error ("FAILED");
 205           endif
 206         catch
 207         end
 208         ## Can not easily test eyediagram, scatterplot!!
 209         fprintf("PASSED\n");
 210       endif
 211       if (strcmp(tests,"source") || strcmp(tests,"all"))
 212         fprintf("\n<< Source Coding Package >>\n");
 213         fprintf("  PCM Functions:                            ");
 214         fprintf("Not tested\n");
 215         fprintf("  Quantization Functions:                   ");
 216         x = [0:0.1:2*pi];
 217         y = sin(x);
 218         [tab, cod] = lloyds(y, 16);
 219         [i, q, d] = quantiz(y, tab, cod);
 220         if (abs(d) > 0.1)
 221           error ("FAILED");
 222         endif
 223
 224         mu = 0.1;
 225         V = 1;
 226         x = sin([0:0.1:2*pi]);
 227         y = compand(x, mu, V, "mu/compressor");
 228         z = compand(x, mu, V, "mu/expander");
 229         ## Again this is a pretty arbitrary test
 230         if (max(abs(x-z)) > 0.1)
 231           error ("FAILED");
 232         endif
 233         fprintf("PASSED\n");
 234       endif
 235       if (strcmp(tests,"block") || strcmp(tests,"all"))
 236         fprintf("\n<< Block Coding Package >>\n");
 237         fprintf("  Cyclic Coding:                            ");
 238         nsym = 100;
 239         m = 4;
 240         n = 2^m-1;                      # [15,11] Hamming code
 241         k = n - m;
 242         p = cyclpoly(n,k);
 243         if (bi2de(p) != primpoly(m,"nodisplay"))
 244           error("FAILED");
 245         endif
 246         [par, gen] = cyclgen(n,p);
 247         if (any(any(gen2par(par) != gen)))
 248           error("FAILED");
 249         endif
 250         if (gfweight(gen) != 3)
 251           error("FAILED");
 252         endif
 253
 254         msg = randint(nsym,k);
 255         code = encode(msg,n,k,"cyclic");
 256         noisy = mod(code + randerr(nsym,n), 2);
 257         dec = decode(noisy,n,k,"cyclic");
 258         if (any(any(dec != msg)))
 259           error("FAILED");
 260         endif
 261         try                     # Protect! If bitshift isn't install!!
 262           msg = randint(nsym,1,n);
 263           code = encode(msg,n,k,"cyclic/decimal");
 264           noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
 265           dec = decode(noisy,n,k,"cyclic/decimal");
 266           if (any(dec != msg))
 267             error("FAILED");
 268           endif
 269         catch
 270         end
 271         fprintf("PASSED\n");
 272
 273         fprintf("  Hamming Coding:                           ");
 274         nsym = 100;
 275         m = 4;
 276         [par, gen, n, k] = hammgen (m);
 277         if (any(any(gen2par(par) != gen)))
 278           error("FAILED");
 279         endif
 280         if (gfweight(gen) != 3)
 281           error("FAILED");
 282         endif
 283         msg = randint(nsym,k);
 284         code = encode(msg,n,k,"hamming");
 285         noisy = mod(code + randerr(nsym,n), 2);
 286         dec = decode(noisy,n,k,"hamming");
 287         if (any(any(dec != msg)))
 288           error("FAILED");
 289         endif
 290         try                     # Protect! If bitshift isn't install!!
 291           msg = randint(nsym,1,n);
 292           code = encode(msg,n,k,"hamming/decimal");
 293           noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
 294           dec = decode(noisy,n,k,"hamming/decimal");
 295           if (any(dec != msg))
 296             error("FAILED");
 297           endif
 298         catch
 299         end
 300         fprintf("PASSED\n");
 301
 302         fprintf("  BCH Coding:                               ");
 303         ## Setup
 304         m = 5;
 305         nsym = 100;
 306         p = bchpoly(2^m - 1);
 307         ## Pick a BCH code from the list at random
 308         l = ceil(size(p,1) * rand(1,1));
 309         n = p(l,1);
 310         k = p(l,2);
 311         t = p(l,3);
 312         ## Symbols represented by rows of binary matrix
 313         msg = randint(nsym,k);
 314         code = encode(msg,n,k,"bch");
 315         noisy = mod(code + randerr(nsym,n), 2);
 316         dec = decode(noisy,n,k,"bch");
 317         if (any(any(dec != msg)))
 318           error("FAILED");
 319         endif
 320         try                     # Protect! If bitshift isn't install!!
 321           msg = randint(nsym,1,n);
 322           code = encode(msg,n,k,"bch/decimal");
 323           noisy = mod(code + bitshift(1,randint(nsym,1,n)), n+1);
 324           dec = decode(noisy,n,k,"bch/decimal");
 325           if (any(dec != msg))
 326             error("FAILED");
 327           endif
 328         catch
 329         end
 330         fprintf("PASSED\n");
 331
 332         fprintf("  Reed-Solomon Coding:                      ");
 333         ## Test for a CCSDS like coder, but without dual-basis translation
 334         mrs = 8;
 335         nrs = 2^mrs -1;
 336         krs = nrs - 32;
 337         prs = 391;
 338         fcr = 112;
 339         step = 11;
 340
 341         ## CCSDS generator polynomial
 342         ggrs = rsgenpoly(nrs, krs, prs, fcr, step);
 343
 344         ## Code two blocks
 345         msg = gf(floor(2^mrs*rand(2,krs)),mrs,prs);
 346         cde = rsenc(msg,nrs,krs,ggrs);
 347
 348         ## Introduce errors
 349         noisy = cde + [ 255,0,255,0,255,zeros(1,250); ...
 350                        0,255,0,255,zeros(1,251)];
 351
 352         ## Decode (better to pass fcr and step rather than gg for speed)
 353         dec = rsdec(noisy,nrs,krs,fcr,step);
 354
 355         if (any(dec != msg))
 356           error("FAILED");
 357         endif
 358         fprintf("PASSED\n");
 359       endif
 360       if (strcmp(tests,"convol") || strcmp(tests,"all"))
 361         fprintf("\n<< Convolutional Coding Package >>\n");
 362         fprintf("  Utility functions:                        ");
 363         ## create a trellis, use poly2trellis and test with istrellis
 364         fprintf("Not tested\n");
 365         fprintf("  Coding:                                   ");
 366         ## use convenc, punturing??
 367         fprintf("Not tested\n");
 368         fprintf("  Viterbi:                                  ");
 369         ## use vitdec
 370         fprintf("Not tested\n");
 371       endif
 372       if (strcmp(tests,"modulation") || strcmp(tests,"all"))
 373         fprintf("\n<< Modulation Package >>\n");
 374         fprintf("  Analog Modulation:                        ");
 375         Fs = 100;
 376         t = [0:1/Fs:2];
 377         x = sin(2*pi*t);
 378         xq = x + 1i * cos(2*pi*t);
 379         ## Can not test FM as it doesn't work !!!
 380         if ((max(abs(x - ademodce(amodce(x,Fs,"pm"),Fs,"pm"))) > 0.001) || ...
 381             (max(abs(x - ademodce(amodce(x,Fs,"am"),Fs,"am"))) > 0.001) || ...
 382             (max(abs(xq - ademodce(amodce(xq,Fs,"qam"),Fs,"qam/cmplx"))) > 0.001))
 383           error("FAILED");
 384         endif
 385         fprintf("PASSED\n");
 386         fprintf("  Digital Mapping:                          ");
 387         m = 32;
 388         n = 100;
 389         x = randint(n,n,32);
 390         if ((x != demodmap(modmap(x,1,1,"ask",m),1,1,"ask",m)) || ...
 391             (x != demodmap(modmap(x,1,1,"fsk",m),1,1,"fsk",m)) || ...
 392             (x != demodmap(modmap(x,1,1,"msk"),1,1,"msk")) || ...
 393             (x != demodmap(modmap(x,1,1,"psk",m),1,1,"psk",m)) || ...
 394             (x != demodmap(modmap(x,1,1,"qask",m),1,1,"qask",m)) || ...
 395             (x != demodmap(modmap(x,1,1,"qask/cir", ...
 396                   [floor(m/2), m - floor(m/2)]),1,1,"qask/cir", ...
 397                   [floor(m/2), m - floor(m/2)])))
 398           error("FAILED");
 399         endif
 400         fprintf("PASSED\n");
 401         fprintf("  Digital Modulation:                       ");
 402         fprintf("Not tested\n");
 403       endif
 404       if (strcmp(tests,"special") || strcmp(tests,"all"))
 405         fprintf("\n<< Special Filters Package >>\n");
 406         fprintf("  Hankel/Hilbert:                           ");
 407         ## use hank2sys, hilbiir
 408         fprintf("Not tested\n");
 409         fprintf("  Raised Cosine:                            ");
 410         ## use rcosflt, rcosiir rcosine, rcosfir
 411         fprintf("Not tested\n");
 412       endif
 413       if (strcmp(tests,"galois") || strcmp(tests,"all"))
 414         fprintf("\n<< Galois Fields Package >>\n");
 415         ## Testing of the Galois Fields package
 416         m = 3;                  ## must be greater than 2
 417
 418         fprintf("  Find primitive polynomials:               ");
 419         prims = primpoly(m,"all","nodisplay");
 420         for i=2^m:2^(m+1)-1
 421           if (find(prims == i))
 422             if (!isprimitive(i))
 423               error("Error in primitive polynomials");
 424             endif
 425           else
 426             if (isprimitive(i))
 427               error("Error in primitive polynomials");
 428             endif
 429           endif
 430         end
 431         fprintf("PASSED\n");
 432         fprintf("  Create Galois variables:                  ");
 433         n = 2^m-1;
 434         gempty = gf([],m);
 435         gzero = gf(0,m);
 436         gone = gf(1,m);
 437         gmax = gf(n,m);
 438         grow = gf(0:n,m);
 439         gcol = gf([0:n]',m);
 440         matlen = ceil(sqrt(2^m));
 441         gmat = gf(reshape(mod([0:matlen*matlen-1],2^m),matlen,matlen),m);
 442         fprintf("PASSED\n");
 443         fprintf("  Access Galois structures:                 ");
 444         if (gcol.m != m || gcol.prim_poly != primpoly(m,"min", ...
 445                                                       "nodisplay"))
 446           error("FAILED");
 447         endif
 448         fprintf("PASSED\n");
 449         fprintf("  Miscellaneous functions:                  ");
 450         if (size(gmat) != [matlen, matlen])
 451           error("FAILED");
 452         endif
 453         if (length(grow) != 2^m)
 454           error("FAILED");
 455         endif
 456         if (!any(grow) || all(grow) || any(gzero) || !all(gone))
 457           error("FAILED");
 458         endif
 459         if (isempty(gone) || !isempty(gempty))
 460           error("FAILED");
 461         endif
 462         tmp = diag(grow);
 463         if (size(tmp,1) != 2^m || size(tmp,2) != 2^m)
 464           error("FAILED");
 465         endif
 466         for i=1:2^m
 467           for j=1:2^m
 468             if ((i == j) && (tmp(i,j) != grow(i)))
 469               error("FAILED");
 470             elseif ((i != j) && (tmp(i,j) != 0))
 471               error("FAILED");
 472             endif
 473           end
 474         end
 475         tmp = diag(gmat);
 476         if (length(tmp) != matlen)
 477           error("FAILED");
 478         endif
 479         for i=1:matlen
 480           if (gmat(i,i) != tmp(i))
 481             error("FAILED");
 482           endif
 483         end
 484         tmp = reshape(gmat,prod(size(gmat)),1);
 485         if (length(tmp) != prod(size(gmat)))
 486           error("FAILED");
 487         endif
 488         if (exp(log(gf([1:n],m))) != [1:n])
 489           error("FAILED");
 490         endif
 491         tmp = sqrt(gmat);
 492         if (tmp .* tmp != gmat)
 493           error("FAILED");
 494         endif
 495
 496         fprintf("PASSED\n");
 497         fprintf("  Unary operators:                          ");
 498         tmp = - grow;
 499         if (tmp != grow)
 500           error("FAILED");
 501         endif
 502         tmp = !grow;
 503         if (tmp(1) != 1)
 504           error("FAILED");
 505         endif
 506         if (any(tmp(2:length(tmp))))
 507           error("FAILED");
 508         endif
 509         tmp = gmat';
 510         for i=1:size(gmat,1)
 511           for j=1:size(gmat,2)
 512             if (gmat(i,j) != tmp(j,i))
 513               error("FAILED");
 514             endif
 515           end
 516         end
 517         fprintf("PASSED\n");
 518         fprintf("  Arithmetic operators:                     ");
 519         if (any(gmat + gmat))
 520           error("FAILED");
 521         endif
 522         multbl = gcol * grow;
 523         elsqr = gcol .* gcol;
 524         elsqr2 = gcol .^ gf(2,m);
 525         for i=1:length(elsqr)
 526           if (elsqr(i) != multbl(i,i))
 527             error("FAILED");
 528           endif
 529         end
 530         for i=1:length(elsqr)
 531           if (elsqr(i) != elsqr2(i))
 532             error("FAILED");
 533           endif
 534         end
 535         tmp = grow(2:length(grow)) ./ gcol(2:length(gcol))';
 536         if (length(tmp) != n || any(tmp != ones(1,length(grow)-1)))
 537           error("FAILED");
 538         endif
 539         fprintf("PASSED\n");
 540         fprintf("  Logical operators:                        ");
 541         if (grow(1) != gzero || grow(2) != gone || grow(2^m) != n)
 542           error("FAILED");
 543         endif
 544         if (!(grow(1) == gzero) || any(grow != gcol'))
 545           error("FAILED");
 546         endif
 547         fprintf("PASSED\n");
 548
 549         fprintf("  Polynomial manipulation:                  ");
 550         poly1 = gf([2,4,5,1],3);
 551         poly2 = gf([1,2],3);
 552         sumpoly = poly1 + [0,0,poly2];   ## Already test "+"
 553         mulpoly = conv(poly1, poly2);    ## multiplication
 554         poly3 = [poly,remd] = deconv(mulpoly, poly2);
 555         if (!isequal(poly1,poly3))
 556           error("FAILED");
 557         endif
 558         if (any(remd))
 559           error("FAILED");
 560         endif
 561         x0 = gf([0,1,2,3],3);
 562         y0 = polyval(poly1, x0);
 563         alph = gf(2,3);
 564         y1 = alph * x0.^3 + alph.^2 * x0.^2 + (alph.^2+1) *x0 + 1;
 565         if (!isequal(y0,y1))
 566           error("FAILED");
 567         endif
 568         roots1 = roots(poly1);
 569         ck1 = polyval(poly1, roots1);
 570         if (any(ck1))
 571           error("FAILED");
 572         endif
 573         b = minpol(alph);
 574         bpoly = bi2de(b.x,"left-msb");
 575         if (bpoly != alph.prim_poly)
 576           error("FAILED");
 577         endif
 578         c = cosets(3);
 579         c2 = c{2};
 580         mpol = minpol(c2(1));
 581         for i=2:length(c2)
 582           if (mpol != minpol(c2(i)))
 583             error("FAILED");
 584           endif
 585         end
 586         fprintf("PASSED\n");
 587
 588         fprintf("  Linear Algebra:                           ");
 589         [l,u,p] = lu(gmat);
 590         if (any(l*u-p*gmat))
 591           error("FAILED");
 592         endif
 593         g1 = inv(gmat);
 594         g2 = gmat ^ -1;
 595         if (any(g1*gmat != eye(matlen)))
 596           error("FAILED");
 597         endif
 598         if (any(g1 != g2))
 599           error("FAILED");
 600         endif
 601         matdet = 0;
 602         while (!matdet)
 603           granmat = gf(floor(2^m*rand(matlen)),m);
 604           matdet = det(granmat);
 605         endwhile
 606         matrank = rank(granmat);
 607         smallcol = gf([0:matlen-1],m)';
 608         sol1 = granmat \ smallcol;
 609         sol2 = smallcol' / granmat;
 610         if (any(granmat * sol1 != smallcol))
 611           error("FAILED");
 612         endif
 613         if (any(sol2 * granmat != smallcol'))
 614           error("FAILED");
 615         endif
 616         fprintf("PASSED\n");
 617         fprintf("  Signal Processing functions:              ");
 618         b = gf([1,0,0,1,0,1,0,1],m);
 619         a = gf([1,0,1,1],m);
 620         x = gf([1,zeros(1,99)],m);
 621         y0 = filter(b, a, x);
 622         y1 = conv(grow+1, grow);
 623         y2 = grow * convmtx(grow+1, length(grow));
 624         if (any(y1 != y2))
 625           error("FAILED");
 626         endif
 627         [y3,remd] = deconv(y2, grow+1);
 628         if (any(y3 != grow))
 629           error("FAILED");
 630         endif
 631         if (any(remd))
 632           error("FAILED");
 633         endif
 634         alph = gf(2,m);
 635         x = gf(floor(2^m*rand(n,1)),m);
 636         fm = dftmtx(alph);
 637         ifm = dftmtx(1/alph);
 638         y0 = fft(x);
 639         y1 = fm * x;
 640         if (any(y0 != y1))
 641           error("FAILED");
 642         endif
 643         z0 = ifft(y0);
 644         z1 = ifm * y1;
 645         if (any(z0 != x))
 646           error("FAILED");
 647         endif
 648         if (any(z1 != x))
 649           error("FAILED");
 650         endif
 651         fprintf("PASSED\n");
 652
 653       endif
 654       fprintf("\n");
 655     unwind_protect_cleanup
 656       page_screen_output(pso);
 657     end_unwind_protect
 658   else
 659     usage("comms: Unknown argument");
 660   endif
 661 endfunction
 662
 663 %!test
 664 %! try comms("test");
 665 %! catch disp(lasterr()); end