## Copyright © 1995, 1997, 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} {} RMDinterp hurst (@var{A}, @var{H}, @var{N}, @var{n}, @var{sigma})
##
## Do fractal interpolation on a given vector
##
## Generates a row vector @var{B} of length @math{N*(lenght(A)-1)+1}.
## @var{B} interpolates the vector @var{A} with @math{2^n*(length(A)-1)+1}
## random numbers computed with the Random Midpoint Displacement algorithm,
## which represent one-dimensional fractional Brownian motion with Hurst
## parameter @var{H}.  @var{H} belongs in [0; 1].  The optional parameter
## @var{n} defaults to @math{ceil(log2(N))+3}. The optional parameter
## @var{sigma} is the standard deviation of the increments at distance
## @var{N}. If not given, it is computed as the standard deviation of the
## increments of @var{A}.
##
## For an example, try this several times:
## @code{A = [9.2 9.8 10.1 10.3 9.6 9.7 10.1 9.9 10 9.7 11 11.1 10.9 11];}
## @code{B = RMDinterp (A, 0.4, 50);}
## @code{plot ([0:601], B, "-", [0:13], A, "o");}
##
## @strong{Reference}
## Lau, Erramilli, Wang, Willinger, "Self-similar Traffic Generation: The
## random Midpoint Displacement Algorithm and Its Properties", ICC '95,
## pp. 466-472.
## @end deftypefn
##
## @seealso{hurst, fGn, RMDtraffic}

## Author: Francesco Potortì <Potorti@isti.cnr.it>
## 1.1
## Usage: B = RMDinterp (A, H, N, n, sigma)
## Description: Do fractal interpolation on a given vector

function B = RMDinterp (A, H, N, n, sigma)

  if (nargin < 3 || nargin > 5)
    usage ("RMDinterp (A, H, N, n, sigma)");
  endif
  if (!is_vector(A))
    error ("A should be a column vector");
  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)
    error ("N should be a scalar greater than 0");
  endif
  lA = length (A);
  if (nargin > 3)
    if (!is_scalar(n) || n <= 0 || floor(n) != n)
      error ("n should be an integer greater than 0");
    endif
  else
    n = ceil(log(N)/log(2))+3;
  endif
  if (nargin > 4)
    if (!is_scalar (sigma) || sigma <= 0)
      error ("sigma should be a scalar greater than 0");
    endif
  else
    sigma = std(A(2:lA) - A(1:lA-1));
  endif


  L = 2^n;				# length of the interpolation

  lB = (lA-1)*L+1;			# lenght of vector B
  given = 1:L:lB;			# indices of given points in B
  random = 1:lB;			# all points in B
  random(given) = [];			# points to generate randomly
  B(random) = randn (1, lB-lA);		# init B with unscaled displacements
  B(1:L:lB) = A;			# put in B the given points

  scale = sqrt (1 - 2^(2*H-2)) * sigma;	# scaling factor for step 0
  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:lB-d+1];		# indices of points to change in step i
    B(range) = 0.5 * (B(range-d) + B(range+d)) + B(range)*scale;
  endfor

  newlB = N*(length(A)-1)+1;		# length requested by the user
  newelements = ceil([0:newlB-1]*(lB-1)/(newlB-1)+0.5);
  B = B(newelements);
endfunction
