From help-request at octave dot org Mon Jan 17 11:17:36 2005 Subject: fractal functions From: Francesco Potorti` To: Octave users list Date: Mon, 17 Jan 2005 18:18:40 +0100 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 , 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: Potorti at isti dot cnr dot it 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 -------------------------------------------------------------