## Copyright © 1998, 2003, 2005  Francesco Potortì
##
## 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 2 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, write to the Free Software Foundation,
## Inc., 51 Franklin Street, 5th Floor, Boston, MA 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {} RMDtraffic (@var{M}, @var{a}, @var{H}, @var{n}, @var{T})
##
## generate approximate fractional Gaussian noise with positive values
##
## Generate a random row vector of numbers of length @math{2^n}
## representing a trace of the instantaneous traffic generated by a self
## similar source with mean traffic rate @var{M}, peakedness @var{a} and
## Hurst parameter @var{H} during a time interval @var{T}.
##
## @var{M}, @var{a} and @var{T} should be scalars: @math{M,a,T > 0},
## @math{0 < H < 1}, @math{n > 0}.  @var{T} is optional, defaulting to
## @math{2^n}. @var{M}, @var{a} and @var{H} are only target values for the
## generation. The trace will have generally different statistics, in
## particular the mean will generally be higher than @var{M} due to the
## truncation of the negative samples.
##
## The whole output vector is considered to be of time length @var{T}, so
##   the values contained are the traffic intensities in time intervals of
##   length @math{T*2^-n}.
## @var{M} is the mean value of the traffic rate.
## @var{a} is the peakedness factor, i.e. the ratio of the unit time
##   interval variance to the mean value of the traffic rate.
##
## @strong{Reference}
## Lau, Erramilli, Wang, Willinger, "Self-similar Traffic Generation: The
## random Midpoint Displacement Algorithm and Its Properties", ICC '95,
## pp. 466-472.
##
## @deftypefnx {Function File} {[@var{FGN}, @var{Hist}] =} RMDtraffic (@var{M}, @var{a}, @var{H}, @var{n}, @var{T})
## Return a matrix containing the results of the intermediate steps, for
## didactic or debugging purposes.
## @end deftypefn
##
## @seealso{hurst, fGn, RMDinterp}

## Author: Francesco Potortì <Potorti@isti.cnr.it>
## 1.8
## Usage: [FGN, Hist] = RMDtraffic (M, a, H, n, T)
## Description: Generate approximate fractional Gaussian noise with positive values

function [FGN, Hist] = RMDtraffic (M, a, H, n, T)

  oldrand = (exist("randn") == 0);	# no randn, Octave older than 2.0
  if (oldrand) savedist = rand ("dist"); endif

  if (nargin < 4 || nargin > 5)
    usage ("RMDtraffic (M, a, H, n, T)");
  endif
  if (!is_scalar (M) || M <= 0)
    error ("M should be a scalar greater than 0");
  endif
  if (!is_scalar (a) || a < 0)
    error ("a should be a scalar be not less than 0");
  endif
  if (!is_scalar (H) || H < 0 || H > 1)
    error ("H should be a scalar in the interval [0;1]");
  endif
  if (!is_scalar (n) || n <= 0 || round (n) != n)
    error ("n should be an integer greater than 0");
  endif
  if (nargin < 5)
    T = 2^n;
  else
    if (!is_scalar (T) || T <= 0)
      error ("T should be a scalar greater than 0");
    endif
  endif
  hist = (nargout > 1);			# history enabled

  N = 2^n;				# length of the result vector

  # Create a fractional Brownian motion process approximation using the
  # random midpoint displacement algorithm.
  FBM(1,1) = 0;				# init the FBM trace: first point is 0
  FBM(N+1,1) = 0;			# last point is 0
  if (oldrand)
    rand ("normal");			# create Gaussian random numbers
    FBM(2:N) = rand(N-1, 1);		# middle points: unscaled displacements
  else
    FBM(2:N) = randn(N-1, 1);		# middle points: unscaled displacements
  endif
  scale = sqrt (1 - 2^(2*H-2)) * T^H;	# scaling factor for step 0
  if (hist) Hist(:,1) = [2^n;scale;FBM]; endif

  for i = 1:n				# step n times
    d = 2^(n-i);			# distance between successive points
    scale = scale * 2^-H;		# scaling factor for step i
    range = [d+1:2*d:N-d+1];		# indices of points to change in step i
    FBM(range) = 0.5 * (FBM(range-d) + FBM(range+d)) + FBM(range)*scale;
    if (hist) Hist(:,i+1) = [d;scale;FBM]; endif
  endfor

  # Now differentiate to get the fractional gaussian noise process, scale it
  # to get the requested mean and peakedness, and truncate it to eliminate
  # the negative traffic spikes.
  FGN = FBM(2:N+1) - FBM(1:N);
  if (hist) Hist(:,n+2) = [0;0;0;FGN]; endif
  FGN = M*T/N + sqrt(a*M) * FGN;
  FGN = max (0, FGN');

  if (oldrand) rand (savedist); endif
endfunction
