octave_packages/communications-1.1.1/wgn.m

   1 ## Copyright (C) 2002 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} {@var{y} =} wgn (@var{m},@var{n},@var{p})
  18 ## @deftypefnx {Function File} {@var{y} =} wgn (@var{m},@var{n},@var{p},@var{imp})
  19 ## @deftypefnx {Function File} {@var{y} =} wgn (@var{m},@var{n},@var{p},@var{imp},@var{seed},)
  20 ## @deftypefnx {Function File} {@var{y} =} wgn (@var{...},'@var{type}')
  21 ## @deftypefnx {Function File} {@var{y} =} wgn (@var{...},'@var{output}')
  22 ##
  23 ## Returns a M-by-N matrix @var{y} of white Gaussian noise. @var{p} specifies
  24 ## the power of the output noise, which is assumed to be referenced to an
  25 ## impedance of 1 Ohm, unless @var{imp} explicitly defines the impedance.
  26 ##
  27 ## If @var{seed} is defined then the randn function is seeded with this
  28 ## value.
  29 ##
  30 ## The arguments @var{type} and @var{output} must follow the above numerial
  31 ## arguments, but can be specified in any order. @var{type} specifies the
  32 ## units of @var{p}, and can be 'dB', 'dBW', 'dBm' or 'linear'. 'dB' is
  33 ## in fact the same as 'dBW' and is keep as a misnomer of Matlab. The
  34 ## units of 'linear' are in Watts.
  35 ##
  36 ## The @var{output} variable should be either 'real' or 'complex'. If the
  37 ## output is complex then the power @var{p} is divided equally betwen the
  38 ## real and imaginary parts.
  39 ##
  40 ## @seealso{randn,awgn}
  41 ## @end deftypefn
  42
  43 ## 2003-01-28
  44 ##   initial release
  45
  46 function y = wgn (m, n, p, varargin)
  47
  48   if ((nargin < 3) || (nargin > 7))
  49     error ("usage: wgn(m, n, p, imp, seed, type, output)");
  50   endif
  51
  52   if (!isscalar(m) || !isreal(m) || (m < 0) || !isscalar(n) || ...
  53       !isreal(n) || (n<0))
  54     error ("wgn: matrix dimension error");
  55   endif
  56
  57   type = "dBW";
  58   out = "real";
  59   imp = 1;
  60   seed = [];
  61   narg = 0;
  62
  63   for i=1:length(varargin)
  64     arg = varargin{i};
  65     if (ischar(arg))
  66       if (strcmp(arg,"real"))
  67         out = "real";
  68       elseif (strcmp(arg,"complex"))
  69         out = "complex";
  70       elseif (strcmp(arg,"dB"))
  71         type = "dBW";
  72       elseif (strcmp(arg,"dBW"))
  73         type = "dBW";
  74       elseif (strcmp(arg,"dBm"))
  75         type = "dBm";
  76       elseif (strcmp(arg,"linear"))
  77         type = "linear";
  78       else
  79         error ("wgn: invalid argument");
  80       endif
  81     else
  82       narg++;
  83       switch (narg)
  84         case 1,
  85           imp = arg;
  86         case 2,
  87           seed = arg;
  88         otherwise
  89           error ("wgn: too many arguments");
  90       endswitch
  91     endif
  92   end
  93
  94   if (isempty(imp))
  95     imp = 1;
  96   elseif (!isscalar(imp) || !isreal(imp) || (imp < 0))
  97     error ("wgn: impedance value illegal");
  98   endif
  99
 100   if (!isempty(seed))
 101     if (!isscalar(seed) || !isreal(seed) || (seed < 0) ||
 102       ((seed-floor(seed)) != 0))
 103       error ("wgn: random seed must be integer");
 104     endif
 105   endif
 106
 107   if (!isscalar(p) || !isreal(p))
 108     error("wgn: invalid power");
 109   endif
 110   if (strcmp(type,"linear") && (p < 0))
 111     error("wgn: invalid power");
 112   endif
 113
 114   if (strcmp(type,"dBW"))
 115     np = 10 ^ (p/10);
 116   elseif (strcmp(type,"dBm"))
 117     np = 10 ^((p - 30)/10);
 118   elseif (strcmp(type,"linear"))
 119     np = p;
 120   endif
 121
 122   if(!isempty(seed))
 123     randn("state",seed);
 124   endif
 125
 126   if (strcmp(out,"complex"))
 127     y = (sqrt(imp*np/2))*(randn(m,n)+1i*randn(m,n));
 128   else
 129     y = (sqrt(imp*np))*randn(m,n);
 130   endif
 131
 132 endfunction
 133
 134 ## Allow 30% error in standard deviation, due to randomness
 135 %!error wgn ();
 136 %!error wgn (1);
 137 %!error wgn (1,1);
 138 %!error wgn (1,1,1,1,1,1);
 139 %!assert (isreal(wgn(10,10,30,1,"dBm","real")));
 140 %!assert (iscomplex(wgn(10,10,30,1,"dBm","complex")));
 141 %!assert (abs(std(wgn(10000,1,30,1,"dBm")) - 1) < 0.3);
 142 %!assert (abs(std(wgn(10000,1,0,1,"dBW")) - 1) < 0.3);
 143 %!assert (abs(std(wgn(10000,1,1,1,"linear")) - 1) < 0.3);