RMDtraffic.m: generate Fractional Gaussian Noise

[Francesco Potorti`]

# Put in the public domain by:
# Francesco Potorti` <address@hidden> 1995

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

# Starting from seed, generate a random column vector of numbers of length
# 2^n representing a trace of the instantaneous traffic generated by a
# self similar source with average traffic rate M, burstiness a and Hurst
# parameter H during a time interval T.
#
# M, a and T should be positive numbers, 0<H<1, n > 0.  T is optional,
# defaulting to 2^n.  seed is optional, defaulting to the system time; a
# seed value of 0 also mean use system time as seed.
#
# M, a and H are only target values for the generation.  The trace will
# have generally different statistics.
#
# The whole output vector is considered to be of time length T, so the
#   values contained are the traffic intensities in time intervals of
#   length T/(2^n+1). 
# M is the value of the traffic rate.
# a is the value of the peakednes factor, that is the ratio of variance to
#   the mean of the traffic (M) in a unit time interval. 
#
# Lau, Erramilli, Wang, Willinger, "Self-similar Traffic Generation: The
# random Midpoint Displacement Algorithm and Its Properties", ICC '95,
# pp. 466-472.

  saveseed = "false";
  savedist = rand ("dist");

  if (nargin < 4 || nargin > 6)
    usage ("RNDtraffic (M, a, H, n, T, seed)");
  endif
  if (!is_scalar (M) || M < 0)
    error ("M should be a scalar not less than 0\n");
  endif
  if (!is_scalar (a) || a < 0)
    error ("a should a scalar be not less than 0\n");
  endif
  if (!is_scalar (H) || H < 0 || H > 1)  
    error ("H should be a scalar in the interval [0;1]\n");
  endif
  if (!is_scalar (n) || n <= 0)  
    error ("a should be a scalar greater than 0\n");
  endif
  if (nargin < 5)
    T = 2^n;
  else
    if (!is_scalar (T) || T <= 0)
      error ("T should be a scalar greater than 0\n");
    endif
    if (nargin == 6 && seed != 0)
      rand ("seed", seed);
    endif
  endif
  trace = (nargout > 1);                # tracing enabled

  N = 2^n;                              # length of the result vector
  rand ("normal");                      # create gaussian random numbers

  # Create a fractional brownian motion process approximation using the
  # random midpoint displacement algorithm.
  FBM = [0; rand (N, 1)] * T^H;         # create all the unscaled displacements
  scale = sqrt (1 - 2^(2*H-2));         # scaling factor for step 0
  if (trace) Trace(:,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 (trace) Trace(:,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 (trace) Trace(:,n+2) = [0;0;0;FGN]; endif
  FGN = M*T/N + sqrt (a*M) * FGN;
  FGN = max (0, FGN);
  
  if (is_scalar (saveseed)) rand ("seed", saveseed); endif
  rand (savedist);
endfunction