fractal functions

[Francesco Potorti`]

I just finished reviewing my old code for generating fractal noise and
computing the Hurst parameter of data series.  A presentation of the
functions is available at <http://fly.isti.cnr.it/software/#fractals>,
together with pointers to the function code.

I would gladly discuss improvements to these functions with whomever is
interested.

There are four Octave functions:

 fgn
This one generates fractional Gaussian noise (also known as Hurst noise)
using the spectral synthesis method.

 RMDtraffic
This one does the same, but using a quicker, approximate algorithm known
as random midpoint displacement.  It is written for easily producing
synthetic traces of long-range correlated network traffic for simulation
purposes.

 hurstVTP
This one computes the Hurst parameter of a data series, using the
variance-time plot method.  It computes a correlation coefficient for
estimating the reliability of the measurement.

 hurstIDC
This one does the same as above, using the index of dispersion for
counts (also known as peakedness).  The result computed for H is exactly
the same as the previous function, but the correlation coefficient is
generally different, and in most cases more reliable.  One should use
both this and the previous functions, and verify that both give a
correlation coefficient near 1.

The last two functions should maybe be bundled into one.

Similar functions could be easily written that use the rescaled range
statistics and the more accurate power spectrum analysis methods.  I
will do that if there is request for such features.

The "hurst" function currently distributed with Octave does the same
kind of computation, but in a rather crude way that does a rather
approximate job for values of H far from 0.5; also, it gives no
indication about the reliability of the computed result.  On the other
hand, it is very fast and works on matrices, rather than vectors, so it
is useful for data series whose Hurst parameter is not too far from 0.5
and whose fractal nature is known by other means.

For comparison, here are some results on how well the three H-meter
functions perform on some fractal and non-fractal data series.  The
column of results should be read like this:

column 1: the target H value
column 2: H computed using hurst (the one distributed with Octave)
column 3: H computed using hurstVTP
column 4: H computed using hurstIDC
column 5: correlation coefficient given by hurstVTP
column 6: correlation coefficient given by hurstIDC

This is how the functions behave with a fractional Gaussian noise of
2^20 samples.  The correlation coefficient is very near to 1 for both
VTP and IDC, indicating that the data seem to be fractal, apart from the
case H=0.5 in IDC, which is generally suspicious of data that do not
exhibit correlation between samples, which is the case for H=0.5.
Notice the second column, where the standard Octave hurst function does
not give accurate results for H far from 0.5.

octave> for H=0.1:0.1:0.9; a=fGn(H,20);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.10000  0.23751  0.10734  0.10734  0.95027  0.99601
  0.20000  0.30154  0.22623  0.22623  0.99047  0.99346
  0.30000  0.37557  0.31739  0.31739  0.99161  0.97528
  0.40000  0.42393  0.41069  0.41069  0.99906  0.98067
  0.50000  0.52204  0.47688  0.47688  0.99860  0.67502
  0.60000  0.56183  0.56030  0.56030  0.99738  0.82930
  0.70000  0.66740  0.69656  0.69656  0.99994  0.99930
  0.80000  0.77144  0.79488  0.79487  0.99992  0.99939
  0.90000  0.81954  0.85570  0.85570  0.99978  0.99875

In the previous cases results for H=0.9 are inaccurate, because high H
require a higher number of samples to be realiably measured.  Using 2^23
samples rather than 2^20 yields more precise results, but you need 1GB
memory to do this quickly:

octave> H=0.9; a=fGn(H,23);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.90000  0.83915  0.89033  0.89030  0.99999  0.99997

This is the same as the fGn case, but now the data are only an
approximation of fGn, because they are created with RMDtraffic:

octave> for H=0.1:0.1:0.9; a=RMDtraffic(10000,1,H,20);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.10000  0.25685  0.15169  0.15169  0.96829  0.99375
  0.20000  0.31538  0.20923  0.20923  0.96574  0.98181
  0.30000  0.39203  0.34153  0.34153  0.99428  0.97422
  0.40000  0.41432  0.38293  0.38293  0.99566  0.95630
  0.50000  0.52690  0.50207  0.50207  0.99912  0.09818
  0.60000  0.56717  0.59526  0.59525  0.99989  0.99570
  0.70000  0.70963  0.72213  0.72213  0.99970  0.99680
  0.80000  0.69264  0.73018  0.73018  0.99837  0.98391
  0.90000  0.79236  0.86371  0.86372  0.99991  0.99951

The following three cases are examples of what happens when non fractal
data are examined: in the first case the fGn function has been modified
like this:
  filter = (linspace(0, 5*beta*pi, 2*N));
while in the second case, it has been modified like this:
  filter = sin(linspace(0, 5*beta*pi, 2*N));
In all three cases, at least one of hurstIDC or hurstVTp gives an
indications that the results are unreliable, because of a correlation
coefficient that is far from 1.  The stock hurst function cannot
give any hint about this problem.

octave> for H=[0,1,2,3,4,6,7,8,9]/10; a=nonfGn1(H,20);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.10000  0.50306  0.50342  0.50342  0.99917  0.16481
  0.20000  0.51979  0.48926  0.48926  0.99893  0.42905
  0.30000  0.52251  0.51935  0.51935  0.99958  0.78882
  0.40000  0.52166  0.49980  0.49980  0.99915  0.00944
  0.50000  0.52413  0.50689  0.50689  0.99958  0.42480
  0.60000  0.53052  0.50865  0.50865  0.99958  0.50605
  0.70000  0.52084  0.49966  0.49966  0.99933  0.01836
  0.80000  0.51465  0.50755  0.50755  0.99926  0.36003
  0.90000  0.54396  0.50147  0.50147  0.99952  0.09520
octave> for H=[0,1,2,3,4,6,7,8,9]/10; a=nonfGn2(H,20);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.10000  0.19438  -0.01653  -0.01653  -0.35817   0.99653 
  0.20000  0.52547   0.48305   0.48305   0.99560   0.34940 
  0.30000  0.18170  -0.01286  -0.01286  -0.33167   0.99746 
  0.40000  0.51483   0.51304   0.51304   0.99966   0.69952 
  0.50000  0.13640   0.00246   0.00246   0.05014   0.99517 
  0.60000  0.14842  -0.00653  -0.00653  -0.28941   0.99909 
  0.70000  0.54319   0.49815   0.49815   0.99947   0.11376 
  0.80000  0.17482  -0.00743  -0.00743  -0.27003   0.99863 
  0.90000  0.20217   0.00277   0.00277   0.06424   0.99625 
octave> for H=0.1:0.1:0.9; a=max(0,-H+randn(2^20,1).^2);
>       [vh,vr]=hurstVTP(a); [ih,ir]=hurstIDC(a); r=hurst(a);
>       disp([H, r, vh, ih, vr, ir]); end;
  0.10000  0.49133  0.46585  0.46585  0.99507  0.59270
  0.20000  0.49312  0.49306  0.49306  0.99896  0.29497
  0.30000  0.52521  0.52937  0.52937  0.99838  0.69736
  0.40000  0.51067  0.49340  0.49340  0.99921  0.31923
  0.50000  0.49928  0.46205  0.46205  0.99753  0.75944
  0.60000  0.53441  0.51565  0.51565  0.99924  0.61504
  0.70000  0.49999  0.48260  0.48260  0.99949  0.74891
  0.80000  0.49707  0.45911  0.45911  0.99886  0.88108
  0.90000  0.50537  0.48215  0.48215  0.99870  0.58750

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 313 8091
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
Web: http://fly.isti.cnr.it/           Key:   fly.isti.cnr.it/public.key



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------