Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / communications-1.1.1 / wgn.m

diff --git a/octave_packages/communications-1.1.1/wgn.m b/octave_packages/communications-1.1.1/wgn.m
new file mode 100644 (file)
index 0000000..fcacc29
--- /dev/null
+++ b/octave_packages/communications-1.1.1/wgn.m
@@ -0,0 +1,143 @@
+## Copyright (C) 2002 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} {@var{y} =} wgn (@var{m},@var{n},@var{p})
+## @deftypefnx {Function File} {@var{y} =} wgn (@var{m},@var{n},@var{p},@var{imp})
+## @deftypefnx {Function File} {@var{y} =} wgn (@var{m},@var{n},@var{p},@var{imp},@var{seed},)
+## @deftypefnx {Function File} {@var{y} =} wgn (@var{...},'@var{type}')
+## @deftypefnx {Function File} {@var{y} =} wgn (@var{...},'@var{output}')
+##
+## Returns a M-by-N matrix @var{y} of white Gaussian noise. @var{p} specifies
+## the power of the output noise, which is assumed to be referenced to an
+## impedance of 1 Ohm, unless @var{imp} explicitly defines the impedance. 
+##
+## If @var{seed} is defined then the randn function is seeded with this
+## value.
+##
+## The arguments @var{type} and @var{output} must follow the above numerial
+## arguments, but can be specified in any order. @var{type} specifies the
+## units of @var{p}, and can be 'dB', 'dBW', 'dBm' or 'linear'. 'dB' is
+## in fact the same as 'dBW' and is keep as a misnomer of Matlab. The
+## units of 'linear' are in Watts.
+##
+## The @var{output} variable should be either 'real' or 'complex'. If the 
+## output is complex then the power @var{p} is divided equally betwen the
+## real and imaginary parts.
+##
+## @seealso{randn,awgn}
+## @end deftypefn
+
+## 2003-01-28
+##   initial release
+
+function y = wgn (m, n, p, varargin)
+
+  if ((nargin < 3) || (nargin > 7))
+    error ("usage: wgn(m, n, p, imp, seed, type, output)");
+  endif
+    
+  if (!isscalar(m) || !isreal(m) || (m < 0) || !isscalar(n) || ...
+      !isreal(n) || (n<0))
+    error ("wgn: matrix dimension error");
+  endif
+  
+  type = "dBW";
+  out = "real";
+  imp = 1;
+  seed = [];
+  narg = 0;
+
+  for i=1:length(varargin)
+    arg = varargin{i};
+    if (ischar(arg))
+      if (strcmp(arg,"real"))
+        out = "real";  
+      elseif (strcmp(arg,"complex"))
+        out = "complex";  
+      elseif (strcmp(arg,"dB"))
+        type = "dBW";  
+      elseif (strcmp(arg,"dBW"))
+        type = "dBW";  
+      elseif (strcmp(arg,"dBm"))
+        type = "dBm";  
+      elseif (strcmp(arg,"linear"))
+        type = "linear";  
+      else
+        error ("wgn: invalid argument");
+      endif
+    else
+      narg++;
+      switch (narg)
+       case 1,
+         imp = arg;
+       case 2,
+         seed = arg;
+       otherwise
+         error ("wgn: too many arguments");
+      endswitch
+    endif
+  end
+
+  if (isempty(imp))
+    imp = 1;
+  elseif (!isscalar(imp) || !isreal(imp) || (imp < 0))
+    error ("wgn: impedance value illegal");
+  endif
+
+  if (!isempty(seed))
+    if (!isscalar(seed) || !isreal(seed) || (seed < 0) || 
+      ((seed-floor(seed)) != 0))
+      error ("wgn: random seed must be integer");
+    endif
+  endif
+    
+  if (!isscalar(p) || !isreal(p))
+    error("wgn: invalid power");
+  endif
+  if (strcmp(type,"linear") && (p < 0))
+    error("wgn: invalid power");
+  endif
+
+  if (strcmp(type,"dBW"))
+    np = 10 ^ (p/10);
+  elseif (strcmp(type,"dBm"))
+    np = 10 ^((p - 30)/10);
+  elseif (strcmp(type,"linear"))
+    np = p;
+  endif
+
+  if(!isempty(seed))
+    randn("state",seed);
+  endif
+
+  if (strcmp(out,"complex"))
+    y = (sqrt(imp*np/2))*(randn(m,n)+1i*randn(m,n));
+  else
+    y = (sqrt(imp*np))*randn(m,n);
+  endif
+  
+endfunction
+
+## Allow 30% error in standard deviation, due to randomness
+%!error wgn ();
+%!error wgn (1);
+%!error wgn (1,1);
+%!error wgn (1,1,1,1,1,1);
+%!assert (isreal(wgn(10,10,30,1,"dBm","real")));
+%!assert (iscomplex(wgn(10,10,30,1,"dBm","complex")));
+%!assert (abs(std(wgn(10000,1,30,1,"dBm")) - 1) < 0.3);
+%!assert (abs(std(wgn(10000,1,0,1,"dBW")) - 1) < 0.3);
+%!assert (abs(std(wgn(10000,1,1,1,"linear")) - 1) < 0.3);