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