octave_packages/communications-1.1.1/awgn.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} =} awgn (@var{x},@var{snr})
  18 ## @deftypefnx {Function File} {@var{y} =} awgn (@var{x},@var{snr},@var{pwr})
  19 ## @deftypefnx {Function File} {@var{y} =} awgn (@var{x},@var{snr}, @var{pwr},@var{seed})
  20 ## @deftypefnx {Function File} {@var{y} =} awgn (@var{...}, '@var{type}')
  21 ##
  22 ## Add white Gaussian noise to a voltage signal.
  23 ##
  24 ## The input @var{x} is assumed to be a real or complex voltage  signal. The
  25 ## returned value @var{y} will be the same form and size as @var{x} but with
  26 ## Gaussian noise added. Unless the power is specified in @var{pwr}, the
  27 ## signal power is assumed to be 0dBW, and the noise of @var{snr} dB will be
  28 ## added with respect to this. If @var{pwr} is a numeric value then the signal
  29 ## @var{x} is assumed to be @var{pwr} dBW, otherwise if @var{pwr} is
  30 ## 'measured', then the power in the signal will be measured and the noise
  31 ## added relative to this measured power.
  32 ##
  33 ## If @var{seed} is specified, then the random number generator seed is
  34 ## initialized with this value
  35 ##
  36 ## By default the @var{snr} and @var{pwr} are assumed to be in dB and dBW
  37 ## respectively. This default behaviour can be chosen with @var{type}
  38 ## set to 'dB'. In the case where @var{type} is set to 'linear', @var{pwr}
  39 ## is assumed to be in Watts and @var{snr} is a ratio.
  40 ## @end deftypefn
  41 ## @seealso{randn,wgn}
  42
  43 ## 2003-01-28
  44 ##   initial release
  45
  46 function y = awgn (x, snr, varargin)
  47
  48   if ((nargin < 2) || (nargin > 5))
  49     error ("usage: awgn(x, snr, p, seed, type");
  50   endif
  51
  52   [m,n] = size(x);
  53   if (isreal(x))
  54     out = "real";
  55   else
  56     out = "complex";
  57   endif
  58
  59   p = 0;
  60   seed = [];
  61   type = "dB";
  62   meas = 0;
  63   narg = 0;
  64
  65   for i=1:length(varargin)
  66     arg = varargin{i};
  67     if (ischar(arg))
  68       if (strcmp(arg,"measured"))
  69         meas = 1;
  70       elseif (strcmp(arg,"dB"))
  71         type = "dB";
  72       elseif (strcmp(arg,"linear"))
  73         type = "linear";
  74       else
  75         error ("awgn: invalid argument");
  76       endif
  77     else
  78       narg++;
  79       switch (narg)
  80               case 1,
  81                 p = arg;
  82               case 2,
  83                 seed = arg;
  84               otherwise
  85                 error ("wgn: too many arguments");
  86       endswitch
  87     endif
  88   end
  89
  90   if (isempty(p))
  91     p = 0;
  92   endif
  93
  94   if (!isempty(seed))
  95     if (!isscalar(seed) || !isreal(seed) || (seed < 0) ||
  96         ((seed-floor(seed)) != 0))
  97       error ("awgn: random seed must be integer");
  98     endif
  99   endif
 100
 101   if (!isscalar(p) || !isreal(p))
 102     error("awgn: invalid power");
 103   endif
 104   if (strcmp(type,"linear") && (p < 0))
 105     error("awgn: invalid power");
 106   endif
 107
 108   if (!isscalar(snr) || !isreal(snr))
 109     error("awgn: invalid snr");
 110   endif
 111   if (strcmp(type,"linear") && (snr < 0))
 112     error("awgn: invalid snr");
 113   endif
 114
 115   if(!isempty(seed))
 116     randn("state",seed);
 117   endif
 118
 119   if (meas == 1)
 120     p = sum( abs( x(:)) .^ 2) / length(x(:));
 121     if (strcmp(type,"dB"))
 122       p = 10 * log10(p);
 123     endif
 124   endif
 125
 126   if (strcmp(type,"linear"))
 127     np = p / snr;
 128   else
 129     np = p - snr;
 130   endif
 131
 132   y = x + wgn (m, n, np, 1, seed, type, out);
 133
 134 endfunction
 135
 136                                 %!shared x, y, noisy
 137                                 %!       x = [0:0.01:2*pi]; y = sin (x);
 138                                 %!       noisy = awgn (y, 20, "dB", "measured");
 139
 140 ## Test of noisy is pretty arbitrary, but should pickup most errors
 141                                 %!error awgn ();
 142                                 %!error awgn (1);
 143                                 %!error awgn (1,1,1,1,1);
 144                                 %!assert (isreal(noisy));
 145                                 %!assert (iscomplex(awgn(y+1i,20,"dB","measured")));
 146                                 %!assert (size(y) == size(noisy))
 147                                 %!assert (abs(10*log10(mean(y.^2)/mean((y-noisy).^ 2)) - 20) < 1);