## 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́
## 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