- Email: [email protected]
Singularity spectra of fractional Brownian motions as a multi-fractal
Chaos, Solitons and Fractals 19 (2004) 613–619 www.elsevier.com/locate/chaos
Singularity spectra of fractional Brownian motions as a multi-fractal T.S. Kim b
a,*
, S. Kim
b
a School of Electronic Engineering, Kyungpook National University, Taegu 702-701, Republic of Korea Department of Physics, Nonlinerar and Complex Systems Laboratory, Postech, Pohang 790-784, Republic of Korea
Accepted 16 April 2003 Communicated by M. Wadati
Abstract Fractional Brownian motion acts as a random process with statistical self-similarity in time and self-affinity in shape. From these properties, the complicated patterns can be suitably represented by it with a minimal parameter and less memory. By considering its statistical property through the power spectrum density we can see that this process is not stationary, even though its differential motion is stationary. So in this paper, by taking the wavelet transfom instead of Fourier transformation we investigate its multi-fractal spectrum as a multi-fractal model. 2003 Elsevier Ltd. All rights reserved.
1. Introduction Fractional Brownian motion (FBM) as a generalization of the usual Brownian motion can be used as a model to describe a large number of natural phenomena and shape such as the range of rivers, time series of economic action, terrain surfaces, ripples of water, etc. This process can be defined as a random function obtained by taking the moving average of usual Brownian motion with kernel ðt sÞH 1=2 , where H is the so-called Hurst index, between zero and one [3]. Statistical self-similarity in the scale of time and statistical self-affinity in its shape which are frequently observed in several natural phenomena are the most important properties of FBM. These advantages of FBM make it possible to describe the complicated natural scenes and patterns with a minimal parameter set and less memory. Thus recently there are many attempts to apply FBM to image processing and segmentation, speech wave form analysis, terrain surface modeling, medical imaging analysis and classification, flow modeling and physical data analysis, etc. But due to its complexity as a fractal, the utilization of FBM is not so easy in reality. Usually the behavior of fractal motions can be characterized by their fractal dimensions which represent how rough they are [2]. Even though several types of fractal dimensions can be defined, their calculation is not so easy [6]. Among them, the capacity dimension is usually used by its computational simplicity. In the 1-dimensional case, that is estimated from the relation LðrÞ r1D ; where LðrÞ is the length of the given curve measured by a ruler of size r and D is the capacity (fractal) dimension. For FBM, it can be seen that its fractal dimension D is related to the index H [3,6] by D ¼ 2 H:
*
Corresponding author. E-mail addresses: [email protected] (T.S. Kim), [email protected] (S. Kim).
0960-0779/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0960-0779(03)00187-5
614
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
2. Fractional Brownian motion and fractional Gaussian noise Consider a situation in which a particle moves on a line by jumping a step-length of either n or þn every s seconds. In modeling diffusion we consider n of some microscopic length, say the particle diameter, and s be a microscopic time, say the collision time. Let the step length n be given by a Gaussian or normal probability distribution: 1 n2 pðn; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi exp : ð1Þ 4Ds 4pDs pffiffiffiffiffiffiffiffi By replacing n with n= 2Ds, the Gaussian random process fni g is normalized so that the position of the particle is defined by X ðt ¼ nsÞ ¼
n X
ni :
i¼1
In practice, it is impossible to observe the Brownian motion with infinite resolution. Thus we will consider the particle position only at intervals bs, where b is some arbitrary number. First consider it only at every second time step, i.e. b ¼ 2. Then the increment n in the particle position is the sum of two independent increments n0 and n00 . Then the joint probability distribution pðn0 ; n00 ; sÞ is given by pðn0 ; n00 ; sÞ ¼ pðn0 ; sÞpðn00 ; sÞ as the product of two probability densities of each variable because the two increments are statistically independent. By integrating over all possible combinations of n0 and n00 , which add up to the total increment n, the probability density of n is given by Z 1 1 n2 0 0 0 : ð2Þ pðn; 2sÞ ¼ pðn n ; sÞpðn ; sÞ dn ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4D 2s 4pD 2s 1 This argument can be easily extended to the time interval bs between observations with the result 1 n2 pðn; bsÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4Dbs 4pDbs
ð3Þ
so that the increments constitute an independent Gaussian random process with hni ¼ 0 and hn2 i ¼ 2Dbs. By changing the time scale by a factor b and the length scale by a factor b1=2 in Eq. (1), the following scaling relation is given: pðn^ ¼ b1=2 n; s^ ¼ bsÞ ¼ b1=2 pðn; sÞ:
ð4Þ
From the above process, the probability distribution for the particle position X ðtÞ for the some reference time t0 is also defined by ! 1 ½X ðtÞ X ðt0 Þ2 P ðX ðtÞ X ð0ÞÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ð5Þ 4Djt t0 j 4pjt t0 j such that the following scaling relation P ðb1=2 ½X ðbtÞ X ðbt0 ÞÞ ¼ b1=2 P ðX ðtÞ X ðt0 ÞÞ
ð6Þ
holds. We see that the position X ðtÞ of a Brownian particle is a random function of time t such that hX ðtÞ X ðt0 Þi ¼ 0;
h½X ðtÞ X ðt0 Þ2 i ¼ 2Djt t0 j:
From Eq. (6) it follows that the Brownian paths are statistically self-similar in the sense that the paths ½X ðtÞ X ðt0 Þ and ½X ðbtÞ X ðbt0 Þ are indistinguishable, and the graphs are statistically self-affine in the sense that the scaling factor is different in its two directions, time t and location. The general Brownian motion XH with index H is defined for a normalized independent Gaussian random process fng so that the increment of the position of the Brownian particle must satisfy XH ðtÞ XH ðt0 Þ njt t0 jH
ðt P t0 Þ:
ð7Þ
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
615
In particular, when H ¼ 1=2, this process is an ordinary Brownian motion. In this sense, the FBM with index H has been defined by Z t 1 XH ðtÞ ¼ ðt t0 ÞH1=2 dBðt0 Þ ð8Þ CðH þ 0:5Þ 1 for the gamma function CðxÞ and an ordinary Brownian motion BðtÞ with 0 < H < 1 [3]. By choosing the time variable as an integer and dividing each time unit into n small time steps, it can be approximated as H 1=2 nt X 1 k t n1=2 nk : XH ðtÞ CðH þ 0:5Þ k¼1 n Since the integral in Eq. (8) does not exist as t0 ! 1, the definition can be modified by the modified kernel ðt t0 ÞH 1=2 ; 0 6 t0 6 t; Kðt t0 Þ ¼ 0 H 1=2 0 H 1=2 ðt t Þ ðt Þ ; t0 < 0
ð9Þ
so that XH ðtÞ XH ð0Þ ¼
1 CðH þ 0:5Þ
1 ¼ CðH þ 0:5Þ
Z
t
Kðt t0 Þ dBðt0 Þ
1
Z
0
½t s
H 1=2
ð sÞ
1
H1=2
dBðsÞ þ
Z
t
ðt sÞH1=2 dBðsÞ:
ð10Þ
0
Since Kðbt bt0 Þ ¼ bH 1=2 Kðt t0 Þ, we have XH ðbtÞ XH ð0Þ ¼ bH fXH ðtÞ XH ð0Þg:
ð11Þ
In particular, putting t ¼ 1 and Dt ¼ bt, we see that the increment in the position of the fractional Brownian particle is proportional to jDtjH : XH ðDtÞ XH ð0Þ ¼ jDtjH fXH ð1Þ XH ð0Þg jDtjH :
ð12Þ
For convenience, assume that XH ð0Þ ¼ 0 and use units so that s ¼ 1 and 2Ds ¼ 1. Then it can be verified [2] that the covariance function CXH ðs; tÞ of XH ðtÞ is given by VH 2H ð13Þ jsj þ jtj2H jt sj2H ; CXH ðs; tÞ ¼ E½XH ðsÞXH ðtÞ ¼ 2 where VH ¼ Var½XH ð1Þ ¼
Cð2 2H Þ cosðpH Þ : pH ð2H 1Þ
ð14Þ
Thus the correlation function of future increments XH ðtÞ with past increments XH ðtÞ may be written CðtÞ ¼
hXH ðtÞXH ðtÞi hXH ðtÞ2 i
¼ 22H 2:
In case that H ¼ 1=2 the correlation between past and future increments is zero and the given FBM is an independent random process. However, for H 6¼ 1=2 past increments are correlated with future increments. We can see that there exists an increment for H > 1=2 and a decrement for H < 1=2 on the average in future when there is a positive increment for some time in the past. Though FBM is not independent for H 6¼ 1=2, the increments of FBM are stationary with zero-mean and self-similar such that Var½XH ðtÞ XH ðsÞÞ ¼ VH ðt sÞ2H : From this equation we can estimate the value H from the sequence of observations corresponding to FBM. Now consider the discrete fractional Gaussian noise (DFGN) process NH ðnÞ ¼ XH ðnÞ XH ðn 1Þ
ð15Þ
defined from the sampled FBM or discrete FBM XH ðnÞ XH ðnTs Þ with the sampling period Ts . Then DFGN is a discrete-time, Gaussian-distributed, stationary procewss with zero mean at all samples and from Eq. (13) its autocorrelation function (ACF) satisfies
616
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
rH ðkÞ ¼ E½NH ðn þ kÞNH ðnÞ ¼
r2 jk þ 1j2H 2jkj2H þ jk 1j2H ; 2
ð16Þ
P 2 where r2 is the variance of NH ðnÞ. Here the ACF of DFGN decays very slowly even for large lags, so that 1 k¼0 jrH ðkÞj diverges. Since the sequence frH ðkÞg is not square summable, we can not find the power spectrum density (PSD) of DFGN by directly applying the Fourier transform onPits ACF [1]. In this sense DFGN can not be regarded as the signal consisting of purely white noise process. Let nðnÞ ¼ 1 NH ðnÞ into k¼0 bðkÞeðn kÞ be the orthonormal projection of the P1 the space generated by the basic normalized error process feðn kÞg. Then in spite of the divergence of jrH ðkÞj2 , i¼0 P1 2 2 i¼0 jb½ij ¼ E½jnðnÞj 6 rH ð0Þ, so there must exist the orthogonal complement signal gðnÞg. Now let us estimate the regular part nðnÞ of NH ðnÞ by LevinsonÕs recursion algorithm [4]. • Step 1: Take ¼ 104 and initialize by setting rH ð1Þ k ¼ 1; a1 ð1Þ ¼ ; q1 ¼ 1 ja1 j2 rH ð0Þ: rH ð0Þ • Step 2: Replace k by k þ 1 and evaluate the followings: P rH ðkÞ j¼k1 ak1 ðjÞrH ðk jÞ j¼1 ; ak ðkÞ ¼ qk1 ak ðiÞ ¼ ak1 ðiÞ ak ðkÞak1 ðk iÞ; i 6 i 6 k 1; qk ¼ 1 jak ðkÞj2 qk1 • Step 3: If
Pk1 i¼1
jak ðiÞ ak1 ðiÞj2 6 , quit; else go to Step 2.
P b H ;k ðnÞ ¼ k1 ak ðiÞNH ðn iÞ of NH ðnÞ is obtained in the Then from k recent past data, the one step ahead prediction N i¼1 b H;1 ðnÞ as k ! 1. Here the normalized error process feðnÞ ¼ xðnÞ=pffiffiffi least mean square sense and converges, say, to N qg b H ;1 ðnÞ is a white Gaussian noise determined by NH ðnÞ. Thus the of the prediction error process xðnÞ ¼ NH ðnÞ N coefficient bðkÞ of nðnÞ can be estimated by 1 X 1 aðiÞrH ðk þ iÞg: b^ðkÞ ¼ E½NH ðnÞeðn kÞ ¼ pffiffiffi frH ðkÞ q i¼1
ð17Þ
From this, the PSD Pn ðnÞ of nðnÞ can also be taken by 2 1 X j2pkf bðkÞ e Pn ðf Þ ¼ : k¼0 The remainder process fgðnÞg is then obtained by gðnÞ ¼ NH ðnÞ nðnÞ. 3. Multi-fractal decomposition of FBM through the wavelet transform As in the case of most random self-similar signals, FBM is not a simple fractal but a multi-fractal such that there exists a multi-fractal invariant measure with singularity distributed over various fractal sets with each local fractal dimension. So we can investigate the given process by investigating the singularity spectrum of its singular invariant measure and by investigating the local self-affinity through the local scaling difference exponent with index H as explained in the previous section. But even the rough computations of these processes are not so easy. The scaling properties of aforementioned process, suggests that the spectra of local self-similar measure is closely related with the coarse signal and self-affinity is also related with the detail signal of the wavelet transform of our FBM process. In this sense it is ideally suitable to use the orthonormal wavelet basis to represent and analyze the FBM. Moreover, this transformation can be implemented in a computationally efficiently manner. Let us cover the support of the locally selfsimilar measure with balls of size r so that lðBi Þ ¼ li denotes the probability of being in the ith ball. Then the scaling singularity index ai describing the local singularity strength is defined as li r ai
as r ! 0:
Let NðaÞ be the number of the balls with local singularity index between a and a þ da then the spectrum of singularity f ðaÞ can be defined as the fractal dimension of the set with a singularity strength a so that
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
617
NðaÞ rf ðaÞ : This formalism leads to a description, f ðaÞ, of a multi-fractal measure in terms of interwoven set with a singularity strength a and the Hausdorff dimension [5]. Also this spectrum of singularity can be related to a set of generalized dimension, Dq through the partition function, Hðr; qÞ, defined by X q Hðr; qÞ ¼ li rðq1ÞDq ¼ rsðqÞ ; ð18Þ i
or equivalently by Dq ¼
1 ln Hðr; qÞ lim : q 1 r!0 ln r
ð19Þ
Then the explicit relationship between the set of dimensions Dq and the singularity spectrum f ðaÞ is given by the following Legendre transform [7]: f ðaÞ ¼ q
d ½ðq 1ÞDq ðq 1ÞDq ¼ qaðqÞ sðqÞ: dq
ð20Þ
Let F be a metric space which supports a probability measure l. Then the wavelet transform of the measure l is defined by Z W ða; bÞ ¼ Wp;h ða; bÞ ¼ ap hððx bÞ=aÞ dlðxÞ; ð21Þ F
where a (>0) and b 2 R are the scale and translation parameters and hðxÞ is an analyzing wavelet which is well localized in both space and frequency with zero mean and p > 0. The scaling property of the wavelet transformation of our defined self-measure in FBM was proposed by W ðka; kb þ x0 Þ kaðx0 Þ W ða; b þ x0 Þ
ð22Þ
with k > 0, x0 ; b 2 R. Here aðx0 Þ is local exponent describing the scaling near a ball Bðx0 ; rÞ of radius r centered at x0 such that lðBðx0 ; rÞÞ raðx0 Þ as r ! 0. All the singularities of the signals can be detected by the local maxima of W ða; bÞ at a given scale a. It is found that the partition function, which is essential in computing the singularity spectrum of fractal signal in a thermodynamic formalism, scales as X Hða; qÞ ¼ jW ða; bi Þjq asðqÞ ; ð23Þ bi
where the sum is taken only over the wavelet transform modulus maxima bi at a given scale a. Since the computation of the modulus maxima by the above continuous wavelet transformation is difficult so we consider the discrete wavelet transform, called a multi-resolution decomposition. For a sequence of increasing resolution j, the approximate signal which approximates the original signal up to the resolution 2i is defined by the set of wavelet coefficients in the orthonormal basis expansion at each stage. While the detail at the stage j is defined as the difference of the information between its approximations at the stage j and j 1. Let Vj be the subspace of L2 ðRÞ generated by the orthonormal basis f/j;n ðtÞ ¼ 2j=2 /ð2j t nÞg and let Wj be the orthogonal complement of Vj in Vjþ1 which is generated by the orthogonal basis fwj;n ðtÞ ¼ 2j=2 wð2j t nÞg. Then the approximate signal of XH ðtÞ at resolution scale j as the orthogonal projection of XH ðtÞ on Vj is defined by Aj XH ðtÞ ¼
1 X
< XH ðtÞ; /j;n ðtÞ > /j;n ðtÞ ¼
1
X
aj;n /j;n ðtÞ
ð24Þ
n
and the detail signal Dj XH ðtÞ at resolution scale j as the orthogonal projection of XH ðtÞ on Wj is defined by Dj XH ðtÞ ¼
1 X 1
< XH ðtÞ; wj;n ðtÞ > wj;n ðtÞ ¼
X n
dj;n wj;n ðtÞ;
ð25Þ
618
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
R1
f ðuÞg ðuÞ du. Then the original signal XH ðtÞ is decomposed as j j X X X X XH ðtÞ ¼ A0 XH ðtÞ ¼ Aj XH ðtÞ þ Dk XH ðkÞ ¼ aj;n /j;n ðtÞ þ dk;n wk;n ðtÞ:
where hf ; gi ¼
1
k¼0
n
k¼0
ð26Þ
n
Since XH ðtÞ has a multi-fractal invariant measure as mentioned before, its approximation signal can only be represented at its support. Let gðtÞ be the projection of XH ðtÞ onto the support of the given self-similar measure. Then the scaling of g is given by gðktÞ ¼ ka gðtÞ:
ð27Þ
Then from Eqs. (24) and (27), XZ XZ Ajþk XH ðtÞ ¼ gðuÞ/jþk;n ðuÞ du /jþk;n ðtÞ ¼ gð2k uÞ/j;n ðuÞdð2k uÞ /jþk;n ðtÞ n
X
n ðaðnÞþ1Þk
2
hgðuÞ; /j;n ðuÞi/jþk;n ðtÞ:
ð28Þ
n
Therefore we have jajþk;n j ð2k ÞaðnÞ jaj;n j:
ð29Þ
By applying the Eqs. (22) and (23) to this cases, the partition function can be expressed in terms of approximate signals such that X jaj;i j q P ð30Þ Hðj; qÞ ¼ n jaj;n j j which scales as Hðj; qÞ ð2j ÞsðqÞ
as j ! 1:
ð31Þ
While the approximate signals are associated with the scaling of the singular invariant measure, we find that the detail signals provide the scaling of the H€ older continuity with exponent H . Note that the exponent H can be measured in terms of the detail signals because they characterize the differences of signals at a small intervals as in Eq. (12). Moreover this detail signals provide a more complicated and powerful tool for signal processing because of their compactness and function as band pass filters. Moreover, the relation jdj;n j ¼ jhXH ; wj;n ij 2jðH þ1=2Þ
ð32Þ
is derived as follows: Z X X XH ðDn Þwj;n ðun Þ DH wj;n ðun Þ ¼ 2jH 2j=2 ¼ 2jðHþ1=2Þ ; dj;n ¼ XH ðuÞwj;n ðuÞ du n
n
n
where Dn ¼ un un1 ¼ 2 . Since the FBM is a multi-affine fractal, H ¼ H ðnÞ must be a function of n. Then the qth order correlation function can be defined in terms of detail signals by cq ðxÞ ¼
N 1 X jdj;n jq N n¼1
ð33Þ
which scales as cq ðxÞ
N 1 X 2jðHþ1=2Þq ð2j ÞnðqÞ : N n¼1
ð34Þ
Therefore we can consider the singularity spectrum nðqÞ of the FBM in terms of detail signals. For the purpose of computational reduction, we define the multi-resolution spectrums of detailed signals and approximated signals such that the total singularity spectrum can be observed in the case of classifying fractional motions from Eqs. (33) and (34) by sq ¼
M X N 1 X jdj;n jq : MN j¼1 n¼1
ð35Þ
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
619
Here s1 and s2 s21 means the the average and the variance of absolute values of all wavelet coefficients, respectively. Since the mth moments sm for m P 3 are not computationally effective, we consider s1 and s2 only. Also by defining appropriate threshold in several levels for the spectrum cq as a multifractal spectrum, we can separate the region of n1 given data as a multifractal decomposition. For example in 2-dimensional cases, considering the low scaled data WLL in level n as an input image, its 2-dimensional wavelet transform can be attained fast by passing it to the low and high pass n n filters horizontally and vertically, respectively and then by downsampling so that we can find four subbands WLL ; WLH ; n n WHL ; and WHH . Its high level transformation can also be taken similarly so that the high level subbanded multiresolution spaces can be used to decompose their local regions through the computation of their moment spectrums.
Acknowledgement This work is supported by the Ministry of Science & Technology in Korea through the National Research Laboratory Project.
References [1] [2] [3] [4] [5] [6] [7]
Barton RJ, Poor HV. Signal detection in fractional Gaussian noise. IEEE Trans Inform Theory 1988;34:943–59. Falconer K. Fractal geometry. Brisbane: John Wiley & Sons; 1990. Feder J. Fractals. New York: Plenum Press; 1988. Haykin S. Adaptive filter theory. Englewood Cliffs, NJ: Prentice-Hall; 1996. Ghez JM, Vaienti S. On the wavelet analysis for multifractal sets. J Stat Phys 1989;57(1):415–20. Kim TS, Kim S. Relations between dimensions and differentiability of curves. Fract Calculus Appl Anal 2001;4(2):135–42. Uhm W, Kim S. Discrete wavelet analysis of multifractal measures and multi-affine signals. J Korean Phys Soc 1998;32(1):1–7.
Singularity spectra of fractional Brownian motions as a multi-fractal T.S. Kim b
a,*
, S. Kim
b
a School of Electronic Engineering, Kyungpook National University, Taegu 702-701, Republic of Korea Department of Physics, Nonlinerar and Complex Systems Laboratory, Postech, Pohang 790-784, Republic of Korea
Accepted 16 April 2003 Communicated by M. Wadati
Abstract Fractional Brownian motion acts as a random process with statistical self-similarity in time and self-affinity in shape. From these properties, the complicated patterns can be suitably represented by it with a minimal parameter and less memory. By considering its statistical property through the power spectrum density we can see that this process is not stationary, even though its differential motion is stationary. So in this paper, by taking the wavelet transfom instead of Fourier transformation we investigate its multi-fractal spectrum as a multi-fractal model. 2003 Elsevier Ltd. All rights reserved.
1. Introduction Fractional Brownian motion (FBM) as a generalization of the usual Brownian motion can be used as a model to describe a large number of natural phenomena and shape such as the range of rivers, time series of economic action, terrain surfaces, ripples of water, etc. This process can be defined as a random function obtained by taking the moving average of usual Brownian motion with kernel ðt sÞH 1=2 , where H is the so-called Hurst index, between zero and one [3]. Statistical self-similarity in the scale of time and statistical self-affinity in its shape which are frequently observed in several natural phenomena are the most important properties of FBM. These advantages of FBM make it possible to describe the complicated natural scenes and patterns with a minimal parameter set and less memory. Thus recently there are many attempts to apply FBM to image processing and segmentation, speech wave form analysis, terrain surface modeling, medical imaging analysis and classification, flow modeling and physical data analysis, etc. But due to its complexity as a fractal, the utilization of FBM is not so easy in reality. Usually the behavior of fractal motions can be characterized by their fractal dimensions which represent how rough they are [2]. Even though several types of fractal dimensions can be defined, their calculation is not so easy [6]. Among them, the capacity dimension is usually used by its computational simplicity. In the 1-dimensional case, that is estimated from the relation LðrÞ r1D ; where LðrÞ is the length of the given curve measured by a ruler of size r and D is the capacity (fractal) dimension. For FBM, it can be seen that its fractal dimension D is related to the index H [3,6] by D ¼ 2 H:
*
Corresponding author. E-mail addresses: [email protected] (T.S. Kim), [email protected] (S. Kim).
0960-0779/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0960-0779(03)00187-5
614
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
2. Fractional Brownian motion and fractional Gaussian noise Consider a situation in which a particle moves on a line by jumping a step-length of either n or þn every s seconds. In modeling diffusion we consider n of some microscopic length, say the particle diameter, and s be a microscopic time, say the collision time. Let the step length n be given by a Gaussian or normal probability distribution: 1 n2 pðn; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi exp : ð1Þ 4Ds 4pDs pffiffiffiffiffiffiffiffi By replacing n with n= 2Ds, the Gaussian random process fni g is normalized so that the position of the particle is defined by X ðt ¼ nsÞ ¼
n X
ni :
i¼1
In practice, it is impossible to observe the Brownian motion with infinite resolution. Thus we will consider the particle position only at intervals bs, where b is some arbitrary number. First consider it only at every second time step, i.e. b ¼ 2. Then the increment n in the particle position is the sum of two independent increments n0 and n00 . Then the joint probability distribution pðn0 ; n00 ; sÞ is given by pðn0 ; n00 ; sÞ ¼ pðn0 ; sÞpðn00 ; sÞ as the product of two probability densities of each variable because the two increments are statistically independent. By integrating over all possible combinations of n0 and n00 , which add up to the total increment n, the probability density of n is given by Z 1 1 n2 0 0 0 : ð2Þ pðn; 2sÞ ¼ pðn n ; sÞpðn ; sÞ dn ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4D 2s 4pD 2s 1 This argument can be easily extended to the time interval bs between observations with the result 1 n2 pðn; bsÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4Dbs 4pDbs
ð3Þ
so that the increments constitute an independent Gaussian random process with hni ¼ 0 and hn2 i ¼ 2Dbs. By changing the time scale by a factor b and the length scale by a factor b1=2 in Eq. (1), the following scaling relation is given: pðn^ ¼ b1=2 n; s^ ¼ bsÞ ¼ b1=2 pðn; sÞ:
ð4Þ
From the above process, the probability distribution for the particle position X ðtÞ for the some reference time t0 is also defined by ! 1 ½X ðtÞ X ðt0 Þ2 P ðX ðtÞ X ð0ÞÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ð5Þ 4Djt t0 j 4pjt t0 j such that the following scaling relation P ðb1=2 ½X ðbtÞ X ðbt0 ÞÞ ¼ b1=2 P ðX ðtÞ X ðt0 ÞÞ
ð6Þ
holds. We see that the position X ðtÞ of a Brownian particle is a random function of time t such that hX ðtÞ X ðt0 Þi ¼ 0;
h½X ðtÞ X ðt0 Þ2 i ¼ 2Djt t0 j:
From Eq. (6) it follows that the Brownian paths are statistically self-similar in the sense that the paths ½X ðtÞ X ðt0 Þ and ½X ðbtÞ X ðbt0 Þ are indistinguishable, and the graphs are statistically self-affine in the sense that the scaling factor is different in its two directions, time t and location. The general Brownian motion XH with index H is defined for a normalized independent Gaussian random process fng so that the increment of the position of the Brownian particle must satisfy XH ðtÞ XH ðt0 Þ njt t0 jH
ðt P t0 Þ:
ð7Þ
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
615
In particular, when H ¼ 1=2, this process is an ordinary Brownian motion. In this sense, the FBM with index H has been defined by Z t 1 XH ðtÞ ¼ ðt t0 ÞH1=2 dBðt0 Þ ð8Þ CðH þ 0:5Þ 1 for the gamma function CðxÞ and an ordinary Brownian motion BðtÞ with 0 < H < 1 [3]. By choosing the time variable as an integer and dividing each time unit into n small time steps, it can be approximated as H 1=2 nt X 1 k t n1=2 nk : XH ðtÞ CðH þ 0:5Þ k¼1 n Since the integral in Eq. (8) does not exist as t0 ! 1, the definition can be modified by the modified kernel ðt t0 ÞH 1=2 ; 0 6 t0 6 t; Kðt t0 Þ ¼ 0 H 1=2 0 H 1=2 ðt t Þ ðt Þ ; t0 < 0
ð9Þ
so that XH ðtÞ XH ð0Þ ¼
1 CðH þ 0:5Þ
1 ¼ CðH þ 0:5Þ
Z
t
Kðt t0 Þ dBðt0 Þ
1
Z
0
½t s
H 1=2
ð sÞ
1
H1=2
dBðsÞ þ
Z
t
ðt sÞH1=2 dBðsÞ:
ð10Þ
0
Since Kðbt bt0 Þ ¼ bH 1=2 Kðt t0 Þ, we have XH ðbtÞ XH ð0Þ ¼ bH fXH ðtÞ XH ð0Þg:
ð11Þ
In particular, putting t ¼ 1 and Dt ¼ bt, we see that the increment in the position of the fractional Brownian particle is proportional to jDtjH : XH ðDtÞ XH ð0Þ ¼ jDtjH fXH ð1Þ XH ð0Þg jDtjH :
ð12Þ
For convenience, assume that XH ð0Þ ¼ 0 and use units so that s ¼ 1 and 2Ds ¼ 1. Then it can be verified [2] that the covariance function CXH ðs; tÞ of XH ðtÞ is given by VH 2H ð13Þ jsj þ jtj2H jt sj2H ; CXH ðs; tÞ ¼ E½XH ðsÞXH ðtÞ ¼ 2 where VH ¼ Var½XH ð1Þ ¼
Cð2 2H Þ cosðpH Þ : pH ð2H 1Þ
ð14Þ
Thus the correlation function of future increments XH ðtÞ with past increments XH ðtÞ may be written CðtÞ ¼
hXH ðtÞXH ðtÞi hXH ðtÞ2 i
¼ 22H 2:
In case that H ¼ 1=2 the correlation between past and future increments is zero and the given FBM is an independent random process. However, for H 6¼ 1=2 past increments are correlated with future increments. We can see that there exists an increment for H > 1=2 and a decrement for H < 1=2 on the average in future when there is a positive increment for some time in the past. Though FBM is not independent for H 6¼ 1=2, the increments of FBM are stationary with zero-mean and self-similar such that Var½XH ðtÞ XH ðsÞÞ ¼ VH ðt sÞ2H : From this equation we can estimate the value H from the sequence of observations corresponding to FBM. Now consider the discrete fractional Gaussian noise (DFGN) process NH ðnÞ ¼ XH ðnÞ XH ðn 1Þ
ð15Þ
defined from the sampled FBM or discrete FBM XH ðnÞ XH ðnTs Þ with the sampling period Ts . Then DFGN is a discrete-time, Gaussian-distributed, stationary procewss with zero mean at all samples and from Eq. (13) its autocorrelation function (ACF) satisfies
616
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
rH ðkÞ ¼ E½NH ðn þ kÞNH ðnÞ ¼
r2 jk þ 1j2H 2jkj2H þ jk 1j2H ; 2
ð16Þ
P 2 where r2 is the variance of NH ðnÞ. Here the ACF of DFGN decays very slowly even for large lags, so that 1 k¼0 jrH ðkÞj diverges. Since the sequence frH ðkÞg is not square summable, we can not find the power spectrum density (PSD) of DFGN by directly applying the Fourier transform onPits ACF [1]. In this sense DFGN can not be regarded as the signal consisting of purely white noise process. Let nðnÞ ¼ 1 NH ðnÞ into k¼0 bðkÞeðn kÞ be the orthonormal projection of the P1 the space generated by the basic normalized error process feðn kÞg. Then in spite of the divergence of jrH ðkÞj2 , i¼0 P1 2 2 i¼0 jb½ij ¼ E½jnðnÞj 6 rH ð0Þ, so there must exist the orthogonal complement signal gðnÞg. Now let us estimate the regular part nðnÞ of NH ðnÞ by LevinsonÕs recursion algorithm [4]. • Step 1: Take ¼ 104 and initialize by setting rH ð1Þ k ¼ 1; a1 ð1Þ ¼ ; q1 ¼ 1 ja1 j2 rH ð0Þ: rH ð0Þ • Step 2: Replace k by k þ 1 and evaluate the followings: P rH ðkÞ j¼k1 ak1 ðjÞrH ðk jÞ j¼1 ; ak ðkÞ ¼ qk1 ak ðiÞ ¼ ak1 ðiÞ ak ðkÞak1 ðk iÞ; i 6 i 6 k 1; qk ¼ 1 jak ðkÞj2 qk1 • Step 3: If
Pk1 i¼1
jak ðiÞ ak1 ðiÞj2 6 , quit; else go to Step 2.
P b H ;k ðnÞ ¼ k1 ak ðiÞNH ðn iÞ of NH ðnÞ is obtained in the Then from k recent past data, the one step ahead prediction N i¼1 b H;1 ðnÞ as k ! 1. Here the normalized error process feðnÞ ¼ xðnÞ=pffiffiffi least mean square sense and converges, say, to N qg b H ;1 ðnÞ is a white Gaussian noise determined by NH ðnÞ. Thus the of the prediction error process xðnÞ ¼ NH ðnÞ N coefficient bðkÞ of nðnÞ can be estimated by 1 X 1 aðiÞrH ðk þ iÞg: b^ðkÞ ¼ E½NH ðnÞeðn kÞ ¼ pffiffiffi frH ðkÞ q i¼1
ð17Þ
From this, the PSD Pn ðnÞ of nðnÞ can also be taken by 2 1 X j2pkf bðkÞ e Pn ðf Þ ¼ : k¼0 The remainder process fgðnÞg is then obtained by gðnÞ ¼ NH ðnÞ nðnÞ. 3. Multi-fractal decomposition of FBM through the wavelet transform As in the case of most random self-similar signals, FBM is not a simple fractal but a multi-fractal such that there exists a multi-fractal invariant measure with singularity distributed over various fractal sets with each local fractal dimension. So we can investigate the given process by investigating the singularity spectrum of its singular invariant measure and by investigating the local self-affinity through the local scaling difference exponent with index H as explained in the previous section. But even the rough computations of these processes are not so easy. The scaling properties of aforementioned process, suggests that the spectra of local self-similar measure is closely related with the coarse signal and self-affinity is also related with the detail signal of the wavelet transform of our FBM process. In this sense it is ideally suitable to use the orthonormal wavelet basis to represent and analyze the FBM. Moreover, this transformation can be implemented in a computationally efficiently manner. Let us cover the support of the locally selfsimilar measure with balls of size r so that lðBi Þ ¼ li denotes the probability of being in the ith ball. Then the scaling singularity index ai describing the local singularity strength is defined as li r ai
as r ! 0:
Let NðaÞ be the number of the balls with local singularity index between a and a þ da then the spectrum of singularity f ðaÞ can be defined as the fractal dimension of the set with a singularity strength a so that
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
617
NðaÞ rf ðaÞ : This formalism leads to a description, f ðaÞ, of a multi-fractal measure in terms of interwoven set with a singularity strength a and the Hausdorff dimension [5]. Also this spectrum of singularity can be related to a set of generalized dimension, Dq through the partition function, Hðr; qÞ, defined by X q Hðr; qÞ ¼ li rðq1ÞDq ¼ rsðqÞ ; ð18Þ i
or equivalently by Dq ¼
1 ln Hðr; qÞ lim : q 1 r!0 ln r
ð19Þ
Then the explicit relationship between the set of dimensions Dq and the singularity spectrum f ðaÞ is given by the following Legendre transform [7]: f ðaÞ ¼ q
d ½ðq 1ÞDq ðq 1ÞDq ¼ qaðqÞ sðqÞ: dq
ð20Þ
Let F be a metric space which supports a probability measure l. Then the wavelet transform of the measure l is defined by Z W ða; bÞ ¼ Wp;h ða; bÞ ¼ ap hððx bÞ=aÞ dlðxÞ; ð21Þ F
where a (>0) and b 2 R are the scale and translation parameters and hðxÞ is an analyzing wavelet which is well localized in both space and frequency with zero mean and p > 0. The scaling property of the wavelet transformation of our defined self-measure in FBM was proposed by W ðka; kb þ x0 Þ kaðx0 Þ W ða; b þ x0 Þ
ð22Þ
with k > 0, x0 ; b 2 R. Here aðx0 Þ is local exponent describing the scaling near a ball Bðx0 ; rÞ of radius r centered at x0 such that lðBðx0 ; rÞÞ raðx0 Þ as r ! 0. All the singularities of the signals can be detected by the local maxima of W ða; bÞ at a given scale a. It is found that the partition function, which is essential in computing the singularity spectrum of fractal signal in a thermodynamic formalism, scales as X Hða; qÞ ¼ jW ða; bi Þjq asðqÞ ; ð23Þ bi
where the sum is taken only over the wavelet transform modulus maxima bi at a given scale a. Since the computation of the modulus maxima by the above continuous wavelet transformation is difficult so we consider the discrete wavelet transform, called a multi-resolution decomposition. For a sequence of increasing resolution j, the approximate signal which approximates the original signal up to the resolution 2i is defined by the set of wavelet coefficients in the orthonormal basis expansion at each stage. While the detail at the stage j is defined as the difference of the information between its approximations at the stage j and j 1. Let Vj be the subspace of L2 ðRÞ generated by the orthonormal basis f/j;n ðtÞ ¼ 2j=2 /ð2j t nÞg and let Wj be the orthogonal complement of Vj in Vjþ1 which is generated by the orthogonal basis fwj;n ðtÞ ¼ 2j=2 wð2j t nÞg. Then the approximate signal of XH ðtÞ at resolution scale j as the orthogonal projection of XH ðtÞ on Vj is defined by Aj XH ðtÞ ¼
1 X
< XH ðtÞ; /j;n ðtÞ > /j;n ðtÞ ¼
1
X
aj;n /j;n ðtÞ
ð24Þ
n
and the detail signal Dj XH ðtÞ at resolution scale j as the orthogonal projection of XH ðtÞ on Wj is defined by Dj XH ðtÞ ¼
1 X 1
< XH ðtÞ; wj;n ðtÞ > wj;n ðtÞ ¼
X n
dj;n wj;n ðtÞ;
ð25Þ
618
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
R1
f ðuÞg ðuÞ du. Then the original signal XH ðtÞ is decomposed as j j X X X X XH ðtÞ ¼ A0 XH ðtÞ ¼ Aj XH ðtÞ þ Dk XH ðkÞ ¼ aj;n /j;n ðtÞ þ dk;n wk;n ðtÞ:
where hf ; gi ¼
1
k¼0
n
k¼0
ð26Þ
n
Since XH ðtÞ has a multi-fractal invariant measure as mentioned before, its approximation signal can only be represented at its support. Let gðtÞ be the projection of XH ðtÞ onto the support of the given self-similar measure. Then the scaling of g is given by gðktÞ ¼ ka gðtÞ:
ð27Þ
Then from Eqs. (24) and (27), XZ XZ Ajþk XH ðtÞ ¼ gðuÞ/jþk;n ðuÞ du /jþk;n ðtÞ ¼ gð2k uÞ/j;n ðuÞdð2k uÞ /jþk;n ðtÞ n
X
n ðaðnÞþ1Þk
2
hgðuÞ; /j;n ðuÞi/jþk;n ðtÞ:
ð28Þ
n
Therefore we have jajþk;n j ð2k ÞaðnÞ jaj;n j:
ð29Þ
By applying the Eqs. (22) and (23) to this cases, the partition function can be expressed in terms of approximate signals such that X jaj;i j q P ð30Þ Hðj; qÞ ¼ n jaj;n j j which scales as Hðj; qÞ ð2j ÞsðqÞ
as j ! 1:
ð31Þ
While the approximate signals are associated with the scaling of the singular invariant measure, we find that the detail signals provide the scaling of the H€ older continuity with exponent H . Note that the exponent H can be measured in terms of the detail signals because they characterize the differences of signals at a small intervals as in Eq. (12). Moreover this detail signals provide a more complicated and powerful tool for signal processing because of their compactness and function as band pass filters. Moreover, the relation jdj;n j ¼ jhXH ; wj;n ij 2jðH þ1=2Þ
ð32Þ
is derived as follows: Z X X XH ðDn Þwj;n ðun Þ DH wj;n ðun Þ ¼ 2jH 2j=2 ¼ 2jðHþ1=2Þ ; dj;n ¼ XH ðuÞwj;n ðuÞ du n
n
n
where Dn ¼ un un1 ¼ 2 . Since the FBM is a multi-affine fractal, H ¼ H ðnÞ must be a function of n. Then the qth order correlation function can be defined in terms of detail signals by cq ðxÞ ¼
N 1 X jdj;n jq N n¼1
ð33Þ
which scales as cq ðxÞ
N 1 X 2jðHþ1=2Þq ð2j ÞnðqÞ : N n¼1
ð34Þ
Therefore we can consider the singularity spectrum nðqÞ of the FBM in terms of detail signals. For the purpose of computational reduction, we define the multi-resolution spectrums of detailed signals and approximated signals such that the total singularity spectrum can be observed in the case of classifying fractional motions from Eqs. (33) and (34) by sq ¼
M X N 1 X jdj;n jq : MN j¼1 n¼1
ð35Þ
T.S. Kim, S. Kim / Chaos, Solitons and Fractals 19 (2004) 613–619
619
Here s1 and s2 s21 means the the average and the variance of absolute values of all wavelet coefficients, respectively. Since the mth moments sm for m P 3 are not computationally effective, we consider s1 and s2 only. Also by defining appropriate threshold in several levels for the spectrum cq as a multifractal spectrum, we can separate the region of n1 given data as a multifractal decomposition. For example in 2-dimensional cases, considering the low scaled data WLL in level n as an input image, its 2-dimensional wavelet transform can be attained fast by passing it to the low and high pass n n filters horizontally and vertically, respectively and then by downsampling so that we can find four subbands WLL ; WLH ; n n WHL ; and WHH . Its high level transformation can also be taken similarly so that the high level subbanded multiresolution spaces can be used to decompose their local regions through the computation of their moment spectrums.
Acknowledgement This work is supported by the Ministry of Science & Technology in Korea through the National Research Laboratory Project.
References [1] [2] [3] [4] [5] [6] [7]
Barton RJ, Poor HV. Signal detection in fractional Gaussian noise. IEEE Trans Inform Theory 1988;34:943–59. Falconer K. Fractal geometry. Brisbane: John Wiley & Sons; 1990. Feder J. Fractals. New York: Plenum Press; 1988. Haykin S. Adaptive filter theory. Englewood Cliffs, NJ: Prentice-Hall; 1996. Ghez JM, Vaienti S. On the wavelet analysis for multifractal sets. J Stat Phys 1989;57(1):415–20. Kim TS, Kim S. Relations between dimensions and differentiability of curves. Fract Calculus Appl Anal 2001;4(2):135–42. Uhm W, Kim S. Discrete wavelet analysis of multifractal measures and multi-affine signals. J Korean Phys Soc 1998;32(1):1–7.