## Copyright © 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} {} fGn (@var{H}, @var{n})
##
## generate fractional Gaussian noise using the spectral synthesis method
##
## Generates a fractional Gaussian noise data series with Hurst
## parameter @var{H} and length @math{2^n}, with null mean and unit
## variance by creating its spectrum and doing an inverse transform on it.
## @end deftypefn
##
## @seealso{hurst, RMDtraffic, RMDinterp}

## Author: Francesco Potortì <Potorti@isti.cnr.it>
## 1.3
## Usage: noise = fGn (H, n)
## Description: Generate fractional Gaussian noise

function noise = fGn(H, n)

  ##  N is half the number of frequency points used.
  N = 2^n;

  ## Create frequency white noise
  ##  First element is frequency 0, which we set equal to 0
  ##  Next N-1 elements have Gaussian modulus and uniform angle
  ##  Element N+1 is real
  ##  next N-1 elements are mirrored conjugates.
  FWnoise(1) = 0;
  FWnoise(2:N+1) = randn(1,N) .* exp(i*2*pi*rand(1,N));
  FWnoise(N+1) += conj(FWnoise(N+1));
  FWnoise(N+2:2*N) = conj(FWnoise(N:-1:2));

  ## Create filter with constant slope -beta
  beta = 2*H-1;
  filter = [0, 1:N, N-1:-1:1].^(-beta/2);
  filter(1) = 0;		# dummy value at frequency 0

  ## Fractional Gaussian noise
  noise = ifft(FWnoise .* filter);
  noise = real(noise(1:N));
  noise /= std(noise);		# this is about 2^-(1+beta*(n-2))

endfunction
