- Email: [email protected]
Generating two-dimensional fractional Brownian motion using the fractional Gaussian process (FGp) algorithm
Physica A 311 (2002) 369 – 380
www.elsevier.com/locate/physa
Generating two-dimensional fractional Brownian motion using the fractional Gaussian process (FGp) algorithm Donald R. McGaugheya;∗ , G.J.M. Aitkenb a Royal
Military College of Canada, Box 17000, Station Forces, Kingston, Ont., Canada of Electrical and Computer Engineering, Queen’s University, Kingston, Ont., Canada
b Department
Received 9 June 2001
Abstract Fractional Brownian motion (FBM) is a random fractal that has been used to model many one-, two- and multi-dimensional natural phenomena. The increments process of FBM has a Gaussian distribution and a stationary correlation function. The fractional Gaussian process (FGp) algorithm is an exact algorithm to simulate Gaussian processes that have stationary correlation functions. The approximate second partial derivative of two-dimensional FBM, called 2D fractional Gaussian noise, is found to be a stationary isotropic Gaussian process. In this paper, the expected correlation function for 2D fractional Gaussian noise is derived. The 2D FGp algorithm is used to simulate the approximate second partial derivative of 2D FBM (FBM2) which is then numerically integrated to generate 2D fractional Brownian motion (FBM2). Ensemble averages of surfaces simulated by the FGp2 algorithm show that the correlation function and power spectral density have the desired properties of 2D fractional Brownian motion. Crown c 2002 Published by Elsevier Science B.V. All rights reserved. Copyright PACS: 65C20 models; Numerical methods; 60G18 self-similar processes; 68U20 simulation; 60H99 correlated noise processes Keywords: Fractional Brownian motion; Fractional Gaussian process; Stationary correlation function; Exact simulation
First appeared as conference paper D. McGaughey and G.J.M. Aitken, “Simulating phase screens”, SPIE Vol. 4007, pp. 705 –712, Munich, Germany, 2000. ∗
Corresponding author. Tel.: +1-613-541-6000; fax: +1-613-544-8107. E-mail addresses: [email protected] (D.R. McGaughey), [email protected] (G.J.M. Aitken).
c 2002 Published by Elsevier Science B.V. 0378-4371/02/$ - see front matter Crown Copyright All rights reserved. PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 0 7 7 8 - 1
370
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
1. Introduction Fractional Brownian motion (FBM) as introduced by Mandelbrot et al. [1,2] has Gaussian distributed increments and can be used to model many natural processes [2– 4]. Algorithms commonly used to generate FBM, such as the successive random additions (SRA) and spectrally synthesized FBM (SSFBM) [5,6] are known to be approximate methods [7–9]. The approximate Frst derivative of one-dimensional FBM is called fractional Gaussian noise (FGN) and is a stationary Gaussian process [1,9]. The fractional Gaussian process (FGp) algorithm generates stationary Gaussian processes which possess a prescribed correlation function [7]. It has been shown that this method can generate time series that are exact in principle. That is they are exact if there are no inaccuracies in computer precision and the random numbers generated are genuinely independent and random [7]. The FGp algorithm can be used to generate one-dimensional FGN which can then be numerically integrated to generate FBM. Caccia et al. [9] analyzed one-dimensional FBM simulated by the SRA, SSFBM and FGp algorithms. They concluded that of the algorithms tested, the FGp algorithm simulated one-dimensional FBM most accurately [9]. In this paper, 2D FBM (FBM2) is generated using the 2D FGp algorithm to simulate a stationary Gaussian derivative process, and then integrating this surface in two dimensions gives FBM2. It is shown that the approximate second partial derivative of 2D FBM is a stationary Gaussian process which will be called 2D fractional Gaussian noise (FGN2). This paper reviews the two-dimensional FGp algorithm as described by Wood and Chan [7]. One- and two-dimensional fractional Brownian motion processes are reviewed. A stationary Gaussian process called FGN2 which can be integrated to simulate FBM2 is then derived. Lastly, the correlation and power spectra of simulated two-dimensional surfaces generated with the FGp algorithm are compared to the desired correlation and power spectra, respectively. It is concluded that the FGp algorithm accurately simulates FBM2 with the desired correlation function and power spectrum.
2. Fractional Brownian motion (FBM) Mandelbrot and van Ness deFned FBM as the fractional derivation of Brownian motion [1]. Although FBM is not mean-square diJerentiable, the derivative process is approximated by smoothing the FBM over a Fnite time and taking a Fnite time diJerence. The derivative process is referred to as FGN. Just as white noise is the derivative of a Brownian motion process, FGN is the derivative process of FBM. The correlation of the FGN process is Cxx (n) = E{(BH (m) − BH (m − 1))(BH (m + n) − BH (m + n − 1))} =
2 (|n + 1|2H + |n − 1|2H − 2|n|2H ) 2
(1)
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
371
for 0 ¡ H ¡ 1, where H is the self-similarity parameter and n is the distance between the samples. The power spectral density (PSD) of FBM follows power-law scaling P(f) = cf− ;
(2)
where the exponent can be related to the self-similarity parameter H by [5] = 2H + 1 ;
(3)
where 1 ¡ ¡ 3 since 0 ¡ H ¡ 1. Similarly, the spectral exponent of FGN, , can be related to H by = 2H − 1 :
(4)
FBM can easily be extended to higher dimensions by changing the scalars to vectors and making the variance proportional to the vector modulus in Eq. (1) [5]. The properties of multi-dimensional FBM are [5] 1. the increments BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn ) are zero mean with a Gaussian distribution, 2. the variance of the increments BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn ) depends only on the distance between the points, 2H n E{|BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn )|2 } ˙ (xi − yi )2 ; (5) i=1
3. the SPSD given by
−(2H +n) n (i )2 S(1 ; 2 ; : : : ; n ) ˙
(6)
i=1
scales with power-law scaling with = 2H + n where ˜ = [1 ; 2 ; : : : ; n ] is the n-dimensional spatial frequency vector. FBM2 exhibits power-law scaling of the spatial power spectral density (SPSD) and the increments are isotropic. FBM2 has been used to model the distortion eJect of the earth’s atmosphere on images collected by large telescopes [4,8,10,11]. 3. Two-dimensional FGp The FGp algorithm has be derived for the multi-dimensional case [7] but is presented here for two-dimensions. The FGp algorithm has Fve steps. The procedure for the two-dimensional FGp algorithm (FGp2) is 1. Generate the expected correlation function for a surface at least twice the size desired for a symmetric matrix of distances. 2. Take the FFT2 or DFT2 to generate the desired SPSD. 3. If any of the elements of the SPSD are negative, create an approximate SPSD matrix from the SPSD matrix and use the approximate SPSD in the rest of the algorithm.
372
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
4. Take the square root of all the SPSD elements and multiply element-by-element by a symmetric, Gaussian, complex random matrix. 5. Perform the inverse FFT2 or DFT2 to generate the surface. In Step 1 a symmetric correlation is created. For a surface of dimension [M; N ], the correlations are calculated for all delays at least M points in the x direction and N points in the y direction. The distance along the x-axis is x=(0; 1; 2; : : : ; P; P−1; : : : ; 2; 1) for P even and x = (0; 1; 2; : : : ; P; P; P − 1; : : : ; 2; 1) for P odd where P ¿ M . Notice the symmetric vector of distances along the x-axis is at least twice as long as the desired dimension along the x-axis M . Similarly, the distance along the y-axis is y = (0; 1; 2; : : : ; Q; Q − 1; : : : ; 2; 1) for Q even and y = (0; 1; 2; : : : ; Q; Q; Q − 1; : : : ; 2; 1) for Q odd, where Q ¿ N . The desired correlation is calculated for a two-dimensional grid of shifts using the symmetric x and y vectors for the entry in the grid. The FFT2 or DFT2 of a symmetric, real correlation matrix gives the SPSD using the Wiener–Khinchine Theorem [12]. Thus, the correlation matrix was constructed so that the DFT2 of the desired correlation function gives the desired SPSD. Each entry of the SPSD should be positive. However, for some correlation functions the SPSD has negative entries except for very large surface dimensions [7]. For SPSD with negative values an approximate SPSD is used for the remainder of the FGp2 algorithm. The approximate SPSD (SPSDA ) is created by zeroing negative SPSD entries and scaling all the other entries to keep the power in the matrix the same. The SPSDA is given by SPSDA = 2 SPSDZ ;
(7)
where SPSDz is the SPSD matrix with all negative values zeroed. The power in a SPSD matrix can be found by summing all the SPSD entries. The scaling factor is a ratio of the power in the approximate matrix over the power in the SPSD and is given by [7]
(PSDz (k; l))
k; l : (8) = k; l (PSD(k; l)) The DFT2 is approximated by the square root of the SPSD or SPSDA [12]. The complex random scaling factor used in Step 4 must obey the symmetry properties of a two-dimensional DFT to ensure the surface is real valued. For a real function f(x; y) the DFT2 has the symmetry [13] F(n; m) = F ∗ (−n; −m) ;
(9)
∗
where indicates the complex conjugate. Also, the DFT2 is real where n = 0 or P and m=0 or Q at the same time. Thus the entries in the two-dimensional random matrix are 1. R(n; m) = Xn; m for n = [0; P]; m = [0; Q], 2. R(n; m) = √12 (Xn; m + jYn; m ) for n = [1; : : : ; P]; m = [1; : : : ; Q], 3. R(n; m) = R∗ (2P − n; 2Q − n) for n = [P + 1; : : : ; 2P − 1]; m = [Q; : : : ; 2Q − 1], where Xn; m and Yn; m are zero mean unit variance Gaussian variables. Notice the mean of the complex random variables is zero and the standard deviation is one for both the real and complex random variables.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
373
The square root of the SPSD or SPSDA is multiplied on an element-by-element basis with the random matrix R(n; m). Then the IDFT2 is taken to generate the surface. Analogous to the one dimensional case, only the Frst [N; M ] region of the surface is used so that the resulting surface will not have a symmetric correlation function. For a record length that is a power of two, the FFT and IFFT can be used by the FGp algorithm. The FGp algorithm is O(P log P) when the FFT is used where P = N × M is the number of points in the surface. The FFT is 2P points long and requires 4P memory locations since the FFT has complex valued parameters. 4. Correlation function for derivative process of FBM2 The FGp2 algorithm requires a stationary correlation function. However, the correlation function for FBM2 is nonstationary. A derivative process of FBM2, which has a stationary isotropic correlation function, is required for the FGp2 algorithm. FGp2 can generate a surface of the derivative process which can then be integrated appropriately to generate the FBM2. Many derivative processes exist for a two-dimensional surface; partial derivative in x; (@=@x)F(x; y), and y; (@=@y)F(x; y), the gradient [(@=@x)F(x; y); (@=@y)F(x; y)] and mixed second partial derivatives (@2 =@x@y)F(x; y). Considering Frst the correlation function for partial derivative in x process. The partial derivative in x for a sampled system can be approximated by x=0; BH (0; y) @ BH (x; y) = BH (x; y) = (10) @x BH (x; y) − BH (x − 1; y) x ¿ 0 : By deFning the derivative processes so that the Frst diJerence is the initial row or column value, the surface can be reconstructed from its partial derivative. If only the diJerences from adjacent samples are taken, the initial conditions are lost, the derivative surface has one less sample and the surface cannot be reconstructed from its derivative. The correlation function for the Frst partial derivative in x is derived in Appendix A and given by
2H c C(m; n) = − ( (m + 1)2 + n2 )2H − 2( m2 + n2 ) 2
(11) + ( (m − 1)2 + n2 )2H ; where c is a constant and (m; n) is the distance vector between the two samples. Using the formula the one point correlations √ in the x and y direction, respectively, are C(1; 0) = c[1 − 22H −1 ] and C(0; 1) = c[1 − ( 2)2H ]. Notice these two correlation are not equal. Thus the Frst partial derivative does not give an isotropic correlation function as required for FBM2. Thus, using the correlation function from a partial derivative in x in the FGp2 algorithm will not simulate FGN2. The correlation function of the Frst partial derivative in y can be similarly derived and will similar to the above expression except the ±1 factor will be added to the n variable. The gradient is the vector pair of the partial derivatives in the x and y directions, respectively. FGp2 can generate two independent surfaces for the x and y gradient.
374
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
However, the x and y gradient of FBM2 are not independent and the FGp2 algorithm cannot generate two dependent phase screens with the desired dependence. Thus the gradient is also eliminated as a derivative process to use. The second mixed partial derivative (@2 =@x@y)F(x; y) of FBM2 will have an isotropic correlation function and a Gaussian distribution. Thus, the mixed second partial derivative (@2 =@x@y)F(x; y) of FBM2 is used in the FGp2 algorithm. The mixed second partial derivative (@2 =@x@y)F(x; y) of FBM2 is two-dimensional fractional Gaussian noise (FGN2). The second partial derivative is approximated by @2 BH (x; y) = [BH (x; y) − BH (x − 1; y)] @x@y − [BH (x; y − 1) − BH (x − 1; y − 1)] ;
(12)
where BH (x; −1) = 0 and BH (−1; y) = 0. This can be found by taking the partial derivative with respect to x and then the partial derivative with respect to y on the (@=@x)BH (x; y) surface. The desired correlation of the FGN2 can be derived by deFning two points on the FGN2 surface. The two points are @2 BH (x1 ; y1 ) = BH (x1 ; y1 ) − BH (x1 − 1; y1 ) @x@y − BH (x1 ; y1 − 1) + BH (x1 − 1; y1 − 1)
(13)
@2 BH (x2 ; y2 ) = BH (x2 ; y2 ) − BH (x2 ; y2 − 1) @x@y − BH (x2 ; y2 ) + BH (x2 ; y2 − 1) :
(14)
and
The correlation can be found by taking the expected value of the product of these two points. The product of the terms can be expanded to the sum of squares. The variance of FBM2 depends on the distance between the points as seen in Eq. (5). @2 @2 The expected correlation E{( @x@y BH (x1 ; y1 ))( @x@y BH (x2 ; y2 ))} is [8]
C(m; n) = 12 [2( (m − 1)2 + n2 )2H + 2( (m + 1)2 + n2 )2H
+ 2( m2 + (n − 1)2 )2H + 2( m2 + (n + 1)2 )2H
− 4( m2 + n2 )2H
−( (m − 1)2 + (n − 1)2 )2H − ( (m − 1)2 + (n + 1)2 )2H
−( (m + 1)2 + (n − 1)2 )2H − ( (m + 1)2 + (n + 1)2 )2H ] ; (15) where m = x2 − x1 and n = y2 − y1 .
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
375
5. Testing FGp2 For FGN2 using the correlation derived above, the approximate SPSD (SPSDA ) was always used in our simulations. The surfaces generated by the FGp2 algorithm are not mathematically exact when the approximation is used [7]. Typically, ¿ 0:9998, indicating very little of the power was in the negative SPSD entries and the approximation is very good. The FGN2 has to be integrated twice, once along the x-axis and once along the y-axis, to generate the FBM2. If the FGN2 matrix is not zero mean, integrating will cause large ramp functions in the x or y direction, or in both. Thus, the mean is subtracted from the FBM2 before it is integrated in the y direction. An example of the 64 × 64 FBM2 generated using the FGp2 algorithm is shown in Fig. 1. 5.1. Two-dimensional correlation function The unbiased two-dimensional auto-correlation function is Cxx (m; n) = E{x(k; l)x(k + m; l + n)} N −|m|−1 1 1 = N − |m| M − |n|
M −|n|−1
k=0
x(k; l)x(k + m; l + n) ;
l=0
(16)
where x(k; l) is a real valued function. The two-dimensional, one-step correlation, Cxx (1; 0), is divided into four correlations, even–even, even–odd, odd–odd, and odd–even. The four correlations are based on the indices of the elements in the one delay product x(n; m)x(n + 1; m). A weighted sum of these four correlations, based on the number of points used in the discrete correlation, equals the one delay correlation Cxx (1; 0). If these four one-step correlations are not the same, the process is a nonstationary process. The one delay even–even (EE), even–odd (EO), odd–even (OE) and odd–odd (OO) correlations are N=2−1 1 CEE (1; 0) = (N − 1)=2 (M − 1)=2
CEO (1; 0) =
M=2−1
n=0
m=0
N=2−1 1 (N −1)=2 (M −1)=2
M=2−2
n=0
x(2n; 2m)x(2n + 1; 2m) ;
(17)
x(2n; 2m+1)x(2n+1; 2m+1) ;
m=0
(18) COO (1; 0) =
N=2−2 1
(N −1)=2 (M −1)=2 n=0
M=2−2
m=0
x(2n+1; 2m+1)x(2n+2; 2m+1) ; (19)
376
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
Fig. 1. An example of a 64 × 64 surface of FBM2 generated by integrating FGN2 created by the FGp2 algorithm.
COE (1; 0) =
N=2−2 1
(N − 1)=2 (M − 1)=2
M=2−1
n=0
x(2n + 1; 2m)x(2n + 2; 2m) ;
m=0
(20) where x and x indicate the next lowest and next highest integer, respectively. 5.2. Spatial power spectral density (SPSD) The two-dimensional discrete Fourier transform (DFT2) is given by [12] N −1 M −1 2( 2( n exp −j m x(n; m) exp −j X (k; l) = N M
(21)
n=0 m=0
and the two-dimensional periodogram estimate of the two-dimensional PSD (SPSD) is P(k; l) =
1 |X (k; l)|2 : NM
(22)
The slope of the SPSD on a log–log graph gives an estimate of the power-law scaling () of the power spectrum of the FBM2.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
377
Table 1 The one step auto-correlation Cxx (1; 0) of FBM2 with the correlation calculated in four parts for the FGp2 algorithm Correlation
FGp2
Even, even Even, odd Odd, odd Odd, even Average Expected
− − − − − −
0.1938 0.1964 0.1949 0.1959 0.1953 0.1953
6. Test results 6.1. Correlation The four one-step correlations, (even–even, even–odd, odd–odd and odd–even) were calculated for 100 instances of a 128×128 FBM2 with H = 56 using the FGp2 algorithm. Table 1 shows the four correlations, the average one step correlation and the expected one-step correlation. All four correlations (even–even, even–odd, odd–odd and odd–even) were close to the expected value for the FGp2 algorithms. This is strong evidence that the FBM2 generated by the FGp2 algorithm is stationary. The ensemble average of the auto-correlation function of the second partial derivative of a FBM2 was calculated for a 128 × 128 point surface with H = 56 for the FGp algorithm. The auto-correlation for the Frst Fve lags along the x-axis, along the y-axis and along the diagonal of the surface were plotted against the expected √ auto-correlation from Eq. (15). Note, along the diagonal each lag is a distance of 2 compared to the lag along the axes. The correlation function of FBM2 from the FGp2 algorithm agrees almost exactly with the expected correlation as seen in Fig. 2. 6.2. Spatial power spectral density (SPSD) The SPSD was calculated for 100 instances of a 128 × 128 FBM2 with H = 56 using the FGp2 algorithm. was estimated along the diagonal of the SPSD from each SPSD, from the average of all 100 √ SPSDs, and from √ SPSD averaged in groups of four between the digital frequencies of 2=128 and 0:425 2. Table 2 shows the estimates. The expected is =2H +2=3:667. The estimates from the FGp2 algorithm are very close to the desired values. 7. Conclusions The FGp2 algorithm is an eOcient algorithm for generating FBM2. The correlation function for FGN2 was derived and presented. The two-dimensional algorithm uses
378
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
1
1 AC expected
0.75 corr coef
corr coef
0.75 0.5 0.25
0.5 0.25
0 _0.25
0 0
1
2
3 lag
(a)
4
_0.25
5
(b)
0
1
4
5
2
3 lag
4
5
1 corr coef
0.75 0.5 0.25 0 _0.25
(c)
0
1
2
3 lag
Fig. 2. FGp2 correlation for Fve lags: (a) along the x direction, Cxx (k; 0), (b) along the y direction, Cxx (0; k), (c) along the diagonal, Cxx (k; k).
Table 2 estimates along the diagonal of 128 × 128 FBM2 created by the FGp2 algorithm FGP2 from each PSD from average of 4 PSD from average of 100 PSD
3.535 3.512 3.518
an approximate SPSD matrix in generating FGN2. Our testing indicates that the correlation of the FGN2 agrees almost exactly with the expected correlation even with this approximation. The FGp2 algorithm is recommended as an eOcient and accurate method of generating FBM2. Higher-dimensional FBM surfaces may be possible by generating the mixed partial derivative with one derivative in each dimension with the FGp algorithm and then integrating in each dimension.
Acknowledgements Financial support from NSERC (Canada) and from NATO CRG 951322 are gratefully acknowledged.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
379
Appendix A. Deriving correlation of FBM increments Consider a two-dimensional FBM2 BH (x; y). The approximate Frst partial derivative in x is deFned as @ x=0; BH (0; y) (A.1) BHx (x; y) = BH (x; y) = BH (x; y) − BH (x − 1; y) x ¿ 0 : @x The correlation function for two points of the partial derivative function that are separated by (m; n) samples is (x + m; y + n)} Cxx {m; n} = E{BHx (x; y)BHx
= E{[BH (x; y) − BH (x − 1; y)] [BH (x + m; y + n) − BH (x + m − 1; y + n)]} ;
(A.2)
where E denote the expected value operator. The product in the expected value is expanded term by term to give Cxx {m; n} = E{BH (x; y)BH (x + m; y + n) −BH (x − 1; y)BH (x + m; y + n) −BH (x; y)BH (x + m − 1; y + n) +BH (x − 1; y)BH (x + m − 1; y + n)} :
(A.3)
Each product term is then expanded as a square of diJerences plus two squared terms. The correlation function is then given as Cxx {m; n} =E{− 12 [BH (x; y) − BH (x + m; y + n)]2 + 12 BH2 (x; y) + 12 BH2 (x + m; y + n) + 12 [BH (x − 1; y) − BH (x + m; y + n)]2 − 12 BH2 (x − 1; y) − 12 BH2 (x + m; y + n) + 12 [BH (x; y) − BH (x + m − 1; y + n)]2 − 12 BH2 (x; y) − 12 BH2 (x + m − 1; y + n) − 12 [BH (x − 1; y) − BH (x + m − 1; y + n)]2 + 12 BH2 (x − 1; y) + 12 BH2 (x + m − 1; y + n)} :
(A.4)
Notice all the terms that are not a square of a diJerence cancel. Using Eq. (5), the correlation function for the partial derivative in x of FBM2 is c C(m; n) = − [( (m + 1)2 + n2 )2H 2
− 2( m2 + n2 )2H + ( (m − 1)2 + n2 )2H ] : (A.5)
380
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
References [1] B.B. Mandelbrot, J.V. Ness, Fractional brownian motions, fractional noises and applications, SIAM Rev. 10 (4) (1968) 422–437. [2] B.B. Mandelbrot, Fractal Geometry of Nature, Prentice-Hall, Englewood CliJs, NJ, 1984. [3] D.R. McGaughey, G.J.M. Aitken, Simulating phase screens, in: Astronomical Telescopes, Vol. 4007, SPIE, 2000, pp. 705 –712. [4] C. Schwartz, G. Baum, E.N. Ribak, Turbulence-degraded wavefronts as fractal surfaces, J. Opt. Soc. Am. A 11 (1) (1994) 444–451. [5] H-O. Peitgen, D. Saupe (Eds.), The Science of Fractal Images, Springer, New York, 1988. [6] H.A. Maske, S. Havlin, M. Schwartz, H.E. Stanley, Phys. Rev. E 53 (5) (1996) 5445–5449. [7] A.T.A. Wood, G. Chan, Simulation of stationary Gaussian processes in [0; 1]d , J. Comput. Graph. Stat. 3 (4) (1994) 409–432. [8] D. McGaughey, Spectral modelling and simulation of atmospherically distorted wavefront data, Ph.D. Thesis, Queen’s University, Kingston, Ont., Canada, 2000. [9] D.C. Caccia, D. Percival, M.J. Cannon, G. Raymond, J.B. Bassingthwaigthe, Analyzing exact fractal time series: evaluating dispersional analysis and rescaled range methods, Physica A 246 (1997) 609–632. [10] M.J. Korenberg, D. McGaughey, G.J.M. Aitken, Parallel cascade prediction of turbulence induced wavefront tilt, Electron. Lett. 32 (14) (1996) 1315–1316. [11] D. McGaughey, G.J.M. Aitken, Temporal analysis of stellar wavefront-tilt data, J. Opt. Soc. Am. A 14 (8) (1997) 1967–1974. [12] A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood CliJs, NJ, 1975. [13] J.S. Lim, A.V. Oppenheim, Advanced Topics in Signal Processing, Prentice-Hall, Englewood CliJs, NJ, 1988.
www.elsevier.com/locate/physa
Generating two-dimensional fractional Brownian motion using the fractional Gaussian process (FGp) algorithm Donald R. McGaugheya;∗ , G.J.M. Aitkenb a Royal
Military College of Canada, Box 17000, Station Forces, Kingston, Ont., Canada of Electrical and Computer Engineering, Queen’s University, Kingston, Ont., Canada
b Department
Received 9 June 2001
Abstract Fractional Brownian motion (FBM) is a random fractal that has been used to model many one-, two- and multi-dimensional natural phenomena. The increments process of FBM has a Gaussian distribution and a stationary correlation function. The fractional Gaussian process (FGp) algorithm is an exact algorithm to simulate Gaussian processes that have stationary correlation functions. The approximate second partial derivative of two-dimensional FBM, called 2D fractional Gaussian noise, is found to be a stationary isotropic Gaussian process. In this paper, the expected correlation function for 2D fractional Gaussian noise is derived. The 2D FGp algorithm is used to simulate the approximate second partial derivative of 2D FBM (FBM2) which is then numerically integrated to generate 2D fractional Brownian motion (FBM2). Ensemble averages of surfaces simulated by the FGp2 algorithm show that the correlation function and power spectral density have the desired properties of 2D fractional Brownian motion. Crown c 2002 Published by Elsevier Science B.V. All rights reserved. Copyright PACS: 65C20 models; Numerical methods; 60G18 self-similar processes; 68U20 simulation; 60H99 correlated noise processes Keywords: Fractional Brownian motion; Fractional Gaussian process; Stationary correlation function; Exact simulation
First appeared as conference paper D. McGaughey and G.J.M. Aitken, “Simulating phase screens”, SPIE Vol. 4007, pp. 705 –712, Munich, Germany, 2000. ∗
Corresponding author. Tel.: +1-613-541-6000; fax: +1-613-544-8107. E-mail addresses: [email protected] (D.R. McGaughey), [email protected] (G.J.M. Aitken).
c 2002 Published by Elsevier Science B.V. 0378-4371/02/$ - see front matter Crown Copyright All rights reserved. PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 0 7 7 8 - 1
370
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
1. Introduction Fractional Brownian motion (FBM) as introduced by Mandelbrot et al. [1,2] has Gaussian distributed increments and can be used to model many natural processes [2– 4]. Algorithms commonly used to generate FBM, such as the successive random additions (SRA) and spectrally synthesized FBM (SSFBM) [5,6] are known to be approximate methods [7–9]. The approximate Frst derivative of one-dimensional FBM is called fractional Gaussian noise (FGN) and is a stationary Gaussian process [1,9]. The fractional Gaussian process (FGp) algorithm generates stationary Gaussian processes which possess a prescribed correlation function [7]. It has been shown that this method can generate time series that are exact in principle. That is they are exact if there are no inaccuracies in computer precision and the random numbers generated are genuinely independent and random [7]. The FGp algorithm can be used to generate one-dimensional FGN which can then be numerically integrated to generate FBM. Caccia et al. [9] analyzed one-dimensional FBM simulated by the SRA, SSFBM and FGp algorithms. They concluded that of the algorithms tested, the FGp algorithm simulated one-dimensional FBM most accurately [9]. In this paper, 2D FBM (FBM2) is generated using the 2D FGp algorithm to simulate a stationary Gaussian derivative process, and then integrating this surface in two dimensions gives FBM2. It is shown that the approximate second partial derivative of 2D FBM is a stationary Gaussian process which will be called 2D fractional Gaussian noise (FGN2). This paper reviews the two-dimensional FGp algorithm as described by Wood and Chan [7]. One- and two-dimensional fractional Brownian motion processes are reviewed. A stationary Gaussian process called FGN2 which can be integrated to simulate FBM2 is then derived. Lastly, the correlation and power spectra of simulated two-dimensional surfaces generated with the FGp algorithm are compared to the desired correlation and power spectra, respectively. It is concluded that the FGp algorithm accurately simulates FBM2 with the desired correlation function and power spectrum.
2. Fractional Brownian motion (FBM) Mandelbrot and van Ness deFned FBM as the fractional derivation of Brownian motion [1]. Although FBM is not mean-square diJerentiable, the derivative process is approximated by smoothing the FBM over a Fnite time and taking a Fnite time diJerence. The derivative process is referred to as FGN. Just as white noise is the derivative of a Brownian motion process, FGN is the derivative process of FBM. The correlation of the FGN process is Cxx (n) = E{(BH (m) − BH (m − 1))(BH (m + n) − BH (m + n − 1))} =
2 (|n + 1|2H + |n − 1|2H − 2|n|2H ) 2
(1)
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
371
for 0 ¡ H ¡ 1, where H is the self-similarity parameter and n is the distance between the samples. The power spectral density (PSD) of FBM follows power-law scaling P(f) = cf− ;
(2)
where the exponent can be related to the self-similarity parameter H by [5] = 2H + 1 ;
(3)
where 1 ¡ ¡ 3 since 0 ¡ H ¡ 1. Similarly, the spectral exponent of FGN, , can be related to H by = 2H − 1 :
(4)
FBM can easily be extended to higher dimensions by changing the scalars to vectors and making the variance proportional to the vector modulus in Eq. (1) [5]. The properties of multi-dimensional FBM are [5] 1. the increments BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn ) are zero mean with a Gaussian distribution, 2. the variance of the increments BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn ) depends only on the distance between the points, 2H n E{|BH (x1 ; x2 ; : : : ; x n ) − BH (y1 ; y2 ; : : : ; yn )|2 } ˙ (xi − yi )2 ; (5) i=1
3. the SPSD given by
−(2H +n) n (i )2 S(1 ; 2 ; : : : ; n ) ˙
(6)
i=1
scales with power-law scaling with = 2H + n where ˜ = [1 ; 2 ; : : : ; n ] is the n-dimensional spatial frequency vector. FBM2 exhibits power-law scaling of the spatial power spectral density (SPSD) and the increments are isotropic. FBM2 has been used to model the distortion eJect of the earth’s atmosphere on images collected by large telescopes [4,8,10,11]. 3. Two-dimensional FGp The FGp algorithm has be derived for the multi-dimensional case [7] but is presented here for two-dimensions. The FGp algorithm has Fve steps. The procedure for the two-dimensional FGp algorithm (FGp2) is 1. Generate the expected correlation function for a surface at least twice the size desired for a symmetric matrix of distances. 2. Take the FFT2 or DFT2 to generate the desired SPSD. 3. If any of the elements of the SPSD are negative, create an approximate SPSD matrix from the SPSD matrix and use the approximate SPSD in the rest of the algorithm.
372
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
4. Take the square root of all the SPSD elements and multiply element-by-element by a symmetric, Gaussian, complex random matrix. 5. Perform the inverse FFT2 or DFT2 to generate the surface. In Step 1 a symmetric correlation is created. For a surface of dimension [M; N ], the correlations are calculated for all delays at least M points in the x direction and N points in the y direction. The distance along the x-axis is x=(0; 1; 2; : : : ; P; P−1; : : : ; 2; 1) for P even and x = (0; 1; 2; : : : ; P; P; P − 1; : : : ; 2; 1) for P odd where P ¿ M . Notice the symmetric vector of distances along the x-axis is at least twice as long as the desired dimension along the x-axis M . Similarly, the distance along the y-axis is y = (0; 1; 2; : : : ; Q; Q − 1; : : : ; 2; 1) for Q even and y = (0; 1; 2; : : : ; Q; Q; Q − 1; : : : ; 2; 1) for Q odd, where Q ¿ N . The desired correlation is calculated for a two-dimensional grid of shifts using the symmetric x and y vectors for the entry in the grid. The FFT2 or DFT2 of a symmetric, real correlation matrix gives the SPSD using the Wiener–Khinchine Theorem [12]. Thus, the correlation matrix was constructed so that the DFT2 of the desired correlation function gives the desired SPSD. Each entry of the SPSD should be positive. However, for some correlation functions the SPSD has negative entries except for very large surface dimensions [7]. For SPSD with negative values an approximate SPSD is used for the remainder of the FGp2 algorithm. The approximate SPSD (SPSDA ) is created by zeroing negative SPSD entries and scaling all the other entries to keep the power in the matrix the same. The SPSDA is given by SPSDA = 2 SPSDZ ;
(7)
where SPSDz is the SPSD matrix with all negative values zeroed. The power in a SPSD matrix can be found by summing all the SPSD entries. The scaling factor is a ratio of the power in the approximate matrix over the power in the SPSD and is given by [7]
(PSDz (k; l))
k; l : (8) = k; l (PSD(k; l)) The DFT2 is approximated by the square root of the SPSD or SPSDA [12]. The complex random scaling factor used in Step 4 must obey the symmetry properties of a two-dimensional DFT to ensure the surface is real valued. For a real function f(x; y) the DFT2 has the symmetry [13] F(n; m) = F ∗ (−n; −m) ;
(9)
∗
where indicates the complex conjugate. Also, the DFT2 is real where n = 0 or P and m=0 or Q at the same time. Thus the entries in the two-dimensional random matrix are 1. R(n; m) = Xn; m for n = [0; P]; m = [0; Q], 2. R(n; m) = √12 (Xn; m + jYn; m ) for n = [1; : : : ; P]; m = [1; : : : ; Q], 3. R(n; m) = R∗ (2P − n; 2Q − n) for n = [P + 1; : : : ; 2P − 1]; m = [Q; : : : ; 2Q − 1], where Xn; m and Yn; m are zero mean unit variance Gaussian variables. Notice the mean of the complex random variables is zero and the standard deviation is one for both the real and complex random variables.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
373
The square root of the SPSD or SPSDA is multiplied on an element-by-element basis with the random matrix R(n; m). Then the IDFT2 is taken to generate the surface. Analogous to the one dimensional case, only the Frst [N; M ] region of the surface is used so that the resulting surface will not have a symmetric correlation function. For a record length that is a power of two, the FFT and IFFT can be used by the FGp algorithm. The FGp algorithm is O(P log P) when the FFT is used where P = N × M is the number of points in the surface. The FFT is 2P points long and requires 4P memory locations since the FFT has complex valued parameters. 4. Correlation function for derivative process of FBM2 The FGp2 algorithm requires a stationary correlation function. However, the correlation function for FBM2 is nonstationary. A derivative process of FBM2, which has a stationary isotropic correlation function, is required for the FGp2 algorithm. FGp2 can generate a surface of the derivative process which can then be integrated appropriately to generate the FBM2. Many derivative processes exist for a two-dimensional surface; partial derivative in x; (@=@x)F(x; y), and y; (@=@y)F(x; y), the gradient [(@=@x)F(x; y); (@=@y)F(x; y)] and mixed second partial derivatives (@2 =@x@y)F(x; y). Considering Frst the correlation function for partial derivative in x process. The partial derivative in x for a sampled system can be approximated by x=0; BH (0; y) @ BH (x; y) = BH (x; y) = (10) @x BH (x; y) − BH (x − 1; y) x ¿ 0 : By deFning the derivative processes so that the Frst diJerence is the initial row or column value, the surface can be reconstructed from its partial derivative. If only the diJerences from adjacent samples are taken, the initial conditions are lost, the derivative surface has one less sample and the surface cannot be reconstructed from its derivative. The correlation function for the Frst partial derivative in x is derived in Appendix A and given by
2H c C(m; n) = − ( (m + 1)2 + n2 )2H − 2( m2 + n2 ) 2
(11) + ( (m − 1)2 + n2 )2H ; where c is a constant and (m; n) is the distance vector between the two samples. Using the formula the one point correlations √ in the x and y direction, respectively, are C(1; 0) = c[1 − 22H −1 ] and C(0; 1) = c[1 − ( 2)2H ]. Notice these two correlation are not equal. Thus the Frst partial derivative does not give an isotropic correlation function as required for FBM2. Thus, using the correlation function from a partial derivative in x in the FGp2 algorithm will not simulate FGN2. The correlation function of the Frst partial derivative in y can be similarly derived and will similar to the above expression except the ±1 factor will be added to the n variable. The gradient is the vector pair of the partial derivatives in the x and y directions, respectively. FGp2 can generate two independent surfaces for the x and y gradient.
374
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
However, the x and y gradient of FBM2 are not independent and the FGp2 algorithm cannot generate two dependent phase screens with the desired dependence. Thus the gradient is also eliminated as a derivative process to use. The second mixed partial derivative (@2 =@x@y)F(x; y) of FBM2 will have an isotropic correlation function and a Gaussian distribution. Thus, the mixed second partial derivative (@2 =@x@y)F(x; y) of FBM2 is used in the FGp2 algorithm. The mixed second partial derivative (@2 =@x@y)F(x; y) of FBM2 is two-dimensional fractional Gaussian noise (FGN2). The second partial derivative is approximated by @2 BH (x; y) = [BH (x; y) − BH (x − 1; y)] @x@y − [BH (x; y − 1) − BH (x − 1; y − 1)] ;
(12)
where BH (x; −1) = 0 and BH (−1; y) = 0. This can be found by taking the partial derivative with respect to x and then the partial derivative with respect to y on the (@=@x)BH (x; y) surface. The desired correlation of the FGN2 can be derived by deFning two points on the FGN2 surface. The two points are @2 BH (x1 ; y1 ) = BH (x1 ; y1 ) − BH (x1 − 1; y1 ) @x@y − BH (x1 ; y1 − 1) + BH (x1 − 1; y1 − 1)
(13)
@2 BH (x2 ; y2 ) = BH (x2 ; y2 ) − BH (x2 ; y2 − 1) @x@y − BH (x2 ; y2 ) + BH (x2 ; y2 − 1) :
(14)
and
The correlation can be found by taking the expected value of the product of these two points. The product of the terms can be expanded to the sum of squares. The variance of FBM2 depends on the distance between the points as seen in Eq. (5). @2 @2 The expected correlation E{( @x@y BH (x1 ; y1 ))( @x@y BH (x2 ; y2 ))} is [8]
C(m; n) = 12 [2( (m − 1)2 + n2 )2H + 2( (m + 1)2 + n2 )2H
+ 2( m2 + (n − 1)2 )2H + 2( m2 + (n + 1)2 )2H
− 4( m2 + n2 )2H
−( (m − 1)2 + (n − 1)2 )2H − ( (m − 1)2 + (n + 1)2 )2H
−( (m + 1)2 + (n − 1)2 )2H − ( (m + 1)2 + (n + 1)2 )2H ] ; (15) where m = x2 − x1 and n = y2 − y1 .
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
375
5. Testing FGp2 For FGN2 using the correlation derived above, the approximate SPSD (SPSDA ) was always used in our simulations. The surfaces generated by the FGp2 algorithm are not mathematically exact when the approximation is used [7]. Typically, ¿ 0:9998, indicating very little of the power was in the negative SPSD entries and the approximation is very good. The FGN2 has to be integrated twice, once along the x-axis and once along the y-axis, to generate the FBM2. If the FGN2 matrix is not zero mean, integrating will cause large ramp functions in the x or y direction, or in both. Thus, the mean is subtracted from the FBM2 before it is integrated in the y direction. An example of the 64 × 64 FBM2 generated using the FGp2 algorithm is shown in Fig. 1. 5.1. Two-dimensional correlation function The unbiased two-dimensional auto-correlation function is Cxx (m; n) = E{x(k; l)x(k + m; l + n)} N −|m|−1 1 1 = N − |m| M − |n|
M −|n|−1
k=0
x(k; l)x(k + m; l + n) ;
l=0
(16)
where x(k; l) is a real valued function. The two-dimensional, one-step correlation, Cxx (1; 0), is divided into four correlations, even–even, even–odd, odd–odd, and odd–even. The four correlations are based on the indices of the elements in the one delay product x(n; m)x(n + 1; m). A weighted sum of these four correlations, based on the number of points used in the discrete correlation, equals the one delay correlation Cxx (1; 0). If these four one-step correlations are not the same, the process is a nonstationary process. The one delay even–even (EE), even–odd (EO), odd–even (OE) and odd–odd (OO) correlations are N=2−1 1 CEE (1; 0) = (N − 1)=2 (M − 1)=2
CEO (1; 0) =
M=2−1
n=0
m=0
N=2−1 1 (N −1)=2 (M −1)=2
M=2−2
n=0
x(2n; 2m)x(2n + 1; 2m) ;
(17)
x(2n; 2m+1)x(2n+1; 2m+1) ;
m=0
(18) COO (1; 0) =
N=2−2 1
(N −1)=2 (M −1)=2 n=0
M=2−2
m=0
x(2n+1; 2m+1)x(2n+2; 2m+1) ; (19)
376
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
Fig. 1. An example of a 64 × 64 surface of FBM2 generated by integrating FGN2 created by the FGp2 algorithm.
COE (1; 0) =
N=2−2 1
(N − 1)=2 (M − 1)=2
M=2−1
n=0
x(2n + 1; 2m)x(2n + 2; 2m) ;
m=0
(20) where x and x indicate the next lowest and next highest integer, respectively. 5.2. Spatial power spectral density (SPSD) The two-dimensional discrete Fourier transform (DFT2) is given by [12] N −1 M −1 2( 2( n exp −j m x(n; m) exp −j X (k; l) = N M
(21)
n=0 m=0
and the two-dimensional periodogram estimate of the two-dimensional PSD (SPSD) is P(k; l) =
1 |X (k; l)|2 : NM
(22)
The slope of the SPSD on a log–log graph gives an estimate of the power-law scaling () of the power spectrum of the FBM2.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
377
Table 1 The one step auto-correlation Cxx (1; 0) of FBM2 with the correlation calculated in four parts for the FGp2 algorithm Correlation
FGp2
Even, even Even, odd Odd, odd Odd, even Average Expected
− − − − − −
0.1938 0.1964 0.1949 0.1959 0.1953 0.1953
6. Test results 6.1. Correlation The four one-step correlations, (even–even, even–odd, odd–odd and odd–even) were calculated for 100 instances of a 128×128 FBM2 with H = 56 using the FGp2 algorithm. Table 1 shows the four correlations, the average one step correlation and the expected one-step correlation. All four correlations (even–even, even–odd, odd–odd and odd–even) were close to the expected value for the FGp2 algorithms. This is strong evidence that the FBM2 generated by the FGp2 algorithm is stationary. The ensemble average of the auto-correlation function of the second partial derivative of a FBM2 was calculated for a 128 × 128 point surface with H = 56 for the FGp algorithm. The auto-correlation for the Frst Fve lags along the x-axis, along the y-axis and along the diagonal of the surface were plotted against the expected √ auto-correlation from Eq. (15). Note, along the diagonal each lag is a distance of 2 compared to the lag along the axes. The correlation function of FBM2 from the FGp2 algorithm agrees almost exactly with the expected correlation as seen in Fig. 2. 6.2. Spatial power spectral density (SPSD) The SPSD was calculated for 100 instances of a 128 × 128 FBM2 with H = 56 using the FGp2 algorithm. was estimated along the diagonal of the SPSD from each SPSD, from the average of all 100 √ SPSDs, and from √ SPSD averaged in groups of four between the digital frequencies of 2=128 and 0:425 2. Table 2 shows the estimates. The expected is =2H +2=3:667. The estimates from the FGp2 algorithm are very close to the desired values. 7. Conclusions The FGp2 algorithm is an eOcient algorithm for generating FBM2. The correlation function for FGN2 was derived and presented. The two-dimensional algorithm uses
378
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
1
1 AC expected
0.75 corr coef
corr coef
0.75 0.5 0.25
0.5 0.25
0 _0.25
0 0
1
2
3 lag
(a)
4
_0.25
5
(b)
0
1
4
5
2
3 lag
4
5
1 corr coef
0.75 0.5 0.25 0 _0.25
(c)
0
1
2
3 lag
Fig. 2. FGp2 correlation for Fve lags: (a) along the x direction, Cxx (k; 0), (b) along the y direction, Cxx (0; k), (c) along the diagonal, Cxx (k; k).
Table 2 estimates along the diagonal of 128 × 128 FBM2 created by the FGp2 algorithm FGP2 from each PSD from average of 4 PSD from average of 100 PSD
3.535 3.512 3.518
an approximate SPSD matrix in generating FGN2. Our testing indicates that the correlation of the FGN2 agrees almost exactly with the expected correlation even with this approximation. The FGp2 algorithm is recommended as an eOcient and accurate method of generating FBM2. Higher-dimensional FBM surfaces may be possible by generating the mixed partial derivative with one derivative in each dimension with the FGp algorithm and then integrating in each dimension.
Acknowledgements Financial support from NSERC (Canada) and from NATO CRG 951322 are gratefully acknowledged.
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
379
Appendix A. Deriving correlation of FBM increments Consider a two-dimensional FBM2 BH (x; y). The approximate Frst partial derivative in x is deFned as @ x=0; BH (0; y) (A.1) BHx (x; y) = BH (x; y) = BH (x; y) − BH (x − 1; y) x ¿ 0 : @x The correlation function for two points of the partial derivative function that are separated by (m; n) samples is (x + m; y + n)} Cxx {m; n} = E{BHx (x; y)BHx
= E{[BH (x; y) − BH (x − 1; y)] [BH (x + m; y + n) − BH (x + m − 1; y + n)]} ;
(A.2)
where E denote the expected value operator. The product in the expected value is expanded term by term to give Cxx {m; n} = E{BH (x; y)BH (x + m; y + n) −BH (x − 1; y)BH (x + m; y + n) −BH (x; y)BH (x + m − 1; y + n) +BH (x − 1; y)BH (x + m − 1; y + n)} :
(A.3)
Each product term is then expanded as a square of diJerences plus two squared terms. The correlation function is then given as Cxx {m; n} =E{− 12 [BH (x; y) − BH (x + m; y + n)]2 + 12 BH2 (x; y) + 12 BH2 (x + m; y + n) + 12 [BH (x − 1; y) − BH (x + m; y + n)]2 − 12 BH2 (x − 1; y) − 12 BH2 (x + m; y + n) + 12 [BH (x; y) − BH (x + m − 1; y + n)]2 − 12 BH2 (x; y) − 12 BH2 (x + m − 1; y + n) − 12 [BH (x − 1; y) − BH (x + m − 1; y + n)]2 + 12 BH2 (x − 1; y) + 12 BH2 (x + m − 1; y + n)} :
(A.4)
Notice all the terms that are not a square of a diJerence cancel. Using Eq. (5), the correlation function for the partial derivative in x of FBM2 is c C(m; n) = − [( (m + 1)2 + n2 )2H 2
− 2( m2 + n2 )2H + ( (m − 1)2 + n2 )2H ] : (A.5)
380
D.R. McGaughey, G.J.M. Aitken / Physica A 311 (2002) 369 – 380
References [1] B.B. Mandelbrot, J.V. Ness, Fractional brownian motions, fractional noises and applications, SIAM Rev. 10 (4) (1968) 422–437. [2] B.B. Mandelbrot, Fractal Geometry of Nature, Prentice-Hall, Englewood CliJs, NJ, 1984. [3] D.R. McGaughey, G.J.M. Aitken, Simulating phase screens, in: Astronomical Telescopes, Vol. 4007, SPIE, 2000, pp. 705 –712. [4] C. Schwartz, G. Baum, E.N. Ribak, Turbulence-degraded wavefronts as fractal surfaces, J. Opt. Soc. Am. A 11 (1) (1994) 444–451. [5] H-O. Peitgen, D. Saupe (Eds.), The Science of Fractal Images, Springer, New York, 1988. [6] H.A. Maske, S. Havlin, M. Schwartz, H.E. Stanley, Phys. Rev. E 53 (5) (1996) 5445–5449. [7] A.T.A. Wood, G. Chan, Simulation of stationary Gaussian processes in [0; 1]d , J. Comput. Graph. Stat. 3 (4) (1994) 409–432. [8] D. McGaughey, Spectral modelling and simulation of atmospherically distorted wavefront data, Ph.D. Thesis, Queen’s University, Kingston, Ont., Canada, 2000. [9] D.C. Caccia, D. Percival, M.J. Cannon, G. Raymond, J.B. Bassingthwaigthe, Analyzing exact fractal time series: evaluating dispersional analysis and rescaled range methods, Physica A 246 (1997) 609–632. [10] M.J. Korenberg, D. McGaughey, G.J.M. Aitken, Parallel cascade prediction of turbulence induced wavefront tilt, Electron. Lett. 32 (14) (1996) 1315–1316. [11] D. McGaughey, G.J.M. Aitken, Temporal analysis of stellar wavefront-tilt data, J. Opt. Soc. Am. A 14 (8) (1997) 1967–1974. [12] A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood CliJs, NJ, 1975. [13] J.S. Lim, A.V. Oppenheim, Advanced Topics in Signal Processing, Prentice-Hall, Englewood CliJs, NJ, 1988.