Multifractal formalism for fractal signals : The structure-function approach versus the wavelet-transform modulus-maxima method

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


I. INTRODUCTION
The multifractal formalism [1-9] has been established to account for the statistical scaling properties of singular measures arising in various physical situations [10][11][12][13][14][15][16][17][18].Notable examples include the invariant probability distribution on a strange attractor, the distribution of voltage drop across a random resistor network, the distribution of growth probabilities on the interface of a diffusionlimited aggregate, and the dissipation field in fully developed turbulent flows.This formalism lies upon the determination of the f (a) singularity spectrum [1] which associates the Hausdorff dimension f (a) to the subset of the support of the measure JL where the singularity strength is a: j(a)=dimH{xiJL<Bx(E))-Ea, for E~OJ , (1) where dimH denotes the Hausdorff dimension and Bx(E)  is an E-box centered at x.The so-called thermodynamical analogy provides a natural connection between the j(a) spectrum and an observable spectrum T(q) defined from the power-law behavior of a partition function [1-9] (in the limit E~O): Zq(E)= ~JL(B;(E))9-ET(q)' i (2) where the sum is taken over a partition of the support of the singular measure JL into boxes of size E.This T(q)   spectrum is directly related to the so-called "generalized" fractal dimensions [19-22] Dq =r(q)/(q -1).Using a standard steepest-descent argument to estimate the sum in Eq. (2), in the limit E~O, one gets r(q)=mina[qa-/(a)] .
(3) r(q) and j (a) are thus related by a Legendre transform.Actually, the variables q and r(q) play the same role as the inverse of temperature and the free energy in the thermodynamics [1][2][3][4][5][6][7][8][9] while the (inverse) Legendre transform j(a) =minq [qa -r(q)] (4) indicates that instead of the energy and the entropy, we have a and f (a) as the thermodynamical variables conju- gate to q and r(q).
In the context of the study of fully developed turbulence [2,10,[23][24][25], Parisi and Frisch [26] (PF) have proposed a similar description of the singular aspect of the longitudinal velocity signal v (x).They have called local singularity h (x 0 ) the exponent which characterizes the local scaling behavior of the velocity increment Bv (x 0 , /) around x 0 , in the limit of inertial separation I~o: (5) The singularity spectrum D (h) of the signal is then defined as the function that gives, for a fixed h, the Hausdorff dimension of the set of points x where the exponent h (x) is equal to h.In their original work, PP have actually pioneered the multifractal approach which consists in estimating the D(h) spectrum from the scaling exponents Sp of the p-order structure function [23,24,26], in the limit 1----+0: (6) Thus, by suitably inserting Eq. ( 5) into the above definition, one can bridge D (h) and Sp by a Legendre transform: (7) which is, a priori, the counterpart of the relationship (4) between the f(a) and r(q) spectra of singular measures.Hereafter, we will refer to this method, which consists in computing the D (h) singularity spectrum of a fractal signal by first estimating the exponents Sp [Eq.(6)] and then by Legendre transforming [Eq.(7)], as the structurefunction (SF) method.
The SF approach has raised a rebirth of experimental interest in the study of turbulent flows [27][28][29][30][31].The main feature of the data gathered in the last decade is that they all exhibit a nonlinear dependence of the exponents Sp versus p, the hall mark of multifractality, cor- roborating the results of earlier experimental investigations [32,33] of turbulent signals.But the elaboration of a multifractal formalism to characterize the singular nature of fractal signals is an important goal that is not restricted, in any case, to the context of fully developed turbulence.Fractal signals are commonly encountered in physics or other applied sciences.Well-known examples include all kinds of random walks (e.g., Brownian signals) used to mimic the noisy dynamical behavior observed in various experimental situations [10,[34][35][36][37][38][39][40][41][42], financial time series [36][37][38], geologic shapes [10,39], interfaces devel- oping in far from equilibrium growth processes [16][17][18]40], turbulent signals recorded in fractal growth phenomena (41], and DNA "walk" coding of nucleotide sequences [42].There have been some attempts [40,43] to extend the SF method to other fields than fully developed turbulence.But, as we will demonstrate in the following, even though the SF method can give some conspicuous information on the multifractality of self-affine signals, generally it does not allow a complete characterization of the distribution of singularities; moreover it has fundamental drawbacks which may introduce drastic bias in the estimate of the D (h) singularity spectrum.
In a previous work, we have elaborated on a method which turns out to be a very efficient tool to analyze fractal signals [44][45][46].This method, based on the wavelet- transform modulus-maxima representation of the signal, provides a practical way to determine the entire D (h) spectrum directly from any experimental signal.Moreover, it is likely to be a good candidate to achieve a unified statistical thermodynamic description of singular distributions (including measures and signals).Beyond the academic examples and the applications of our wavelet-based method given in our original work in Refs.[44,45], there was a need to compare its performance with that of the SF method.The aim of this paper is to point out using concrete examples the failures of the SF approach and to demonstrate that the wavelet-transform modulus-maxima (WTMM) method supplies to these drawbacks.
The paper is organized as follows.In Sec.11 we introduce the WTMM method within the mathematical framework of the wavelet transform and its modulusmaxima representation.In Sec.Ill we address the difficult issue of the estimate of the scaling exponents which are related to negative order moments in both the SF and the WTMM methods.In Sec.IV we show to what extent and we explain why the accessible range of singularities is intrinsically restricted when one uses the SF method.We argue theoretically and we demonstrate on specific examples that the WTMM method does not suffer from such limitations.In Sec.V we discuss the effects of the presence of highly regular parts in the signal on the estimate of the D (h) singularity spectrum.We conclude in Sec.VI with some perspectives for future research.

THE WAVELET-TRANSFORM MODULUS-MAXIMA METHOD
In this section, we briefly describe the WTMM method.This method has been introduced and tested numerically in Refs.(44,45].We refer the reader to Ref. [46] for a more rigorous and complete description.The wavelet transform (WT) of a function f consists in decomposing it into elementary space-scale contributions, associated to the so-called wavelets which are constructed from one single function, the analyzing wavelet 1/J, by means of translations and dilations [47][48][49][50][51].The WT of the function f is defined as where a ER+" is a scale parameter and b ER is a space parameter.The analyzing wavelet 1/J is generally chosen to be well localized in both space and frequency.Usually, 1/J is only required to be of zero mean but, for the particu- lar purpose of singularity tracking that is of interest here, we will further require 1/J to be orthogonal to some low- order polynomials [52][53][54]: A class of commonly used real-valued analyzing wavelets [55][56][57][58] which satisfies the above condition is given by the successive derivatives of the Gaussian function: (10) Let us recall the notion of local Holder exponent [52,54,59] which is more general than the one defined by Eq. ( 5).The Holder exponent h (x 0 ) of a function fat x 0 is defined as the largest exponent such that there exists a polynomial P n ( x) of order n that satisfies (11) for x in a neighborhood of x 0 • This definition provides a generalization of Eq. ( 5) to any exponent h characterizing a singularity in a higher derivative of f.For example, h (x 0 ) = 1. 5 implies that the function f is differentiable at x 0 , but its derivative is not: the singularity actually lies in the second derivative of f.This extension leads naturally to a generalization of the D (h) spectrum [Eq.(7)] introduced by PF.Henceforth we will denote D(h) the Hausdorff dimension of the set where the Holder exponent is equal to h: where h is no longer restricted to [0, 1) but a priori can take on positive as well as negative real values [46].
If one uses an analyzing wavelet 1/J that satisfies the condition (9), the local behavior of f in Eq. ( 11) is mir- rored by the wavelet transform which locally behaves like [52-60) in the limit a --+0, provided N satisfies N > h (x 0 ) in Eq. (9).The above equation mainly says that, when investigating the local scaling behavior of the wavelet coefficients computed with an analyzing wavelet whose N first moments vanish, one can generally detect (and estimate) all the Holder exponents off that are smaller than N.
The originality of the WTMM method consists in building a partition function from the modulus maxima of the wavelet transform.These maxima are defined [59,61), at each scale a, as the local maxima of IT, 1 ;[f)(x,a)l considered as a function of x.In Fig. l(b), we show the space-scale arrangement, in the (x,a) plane, of the modulus maxima of the wavelet transform of the fractional Brownian signal illustrated in Fig. l(a).The analyzing wavelet is the second derivative of the Gaussian function [Eq.(10)].These wavelet-transform modulus maxima are disposed on connected curves called maxima lines.Let us define [46) L(a 0 ) the set of all the maxima lines I that satisfy (x,a)El=a :Sa 0 , ( \;fa :Sa 0 , 3(x,a)E/ An important feature of these maxima lines is that, each time the analyzed signal has a local Holder exponent h (x 0 ) < N at the point x 0 , there is at least one maxima line pointing towards x 0 , along which Eq. (13) holds [46,59).In the case of fractal signals, which are typically characterized by a hierarchical distribution of singularities [52][53][54][55][56][57][58][59], we expect that the number of maxima lines will diverge in the limit a --+0.In fact, as emphasized in Refs.[44,45] the branching structure of the WT skeleton defined by the maxima line arrangement in the (x,a) half plane enlightens the hierarchical organization of the singularities.The WTMM method consists in taking advantage of the space-scale partitioning given by this skeleton to define the following partition function [44][45][46]: (c) r(q) vs q as computed with the WTMM method (e) and the SF method (6 ).(d) D(h) vs h from the Legend re transform [Eq.( 17)] of the r( q) spectrum computed with the WTMM method.The analyzing wavelet is the second derivative 1(!< 2 l of the Gaussian function. [ sup IT.p(f](x,a')ljq, (x,a')El (15) where q EIR.Z(a,q) plays a role similar to the partition function defined in Eq. (2) for singular measures.Indeed, a deep analogy [44][45][46]62] does exist between the classical partitions defined for measures and the ones provided by the wavelet-transform modulus-maxima representation.The analyzing wavelet 1/J can be seen as a box of a particular shape, the scale a being its size [ E in Eq. ( 2)], while the modulus maxima indicate how to position our special "boxes" to obtain a partition at the considered scale.[Let us note that the sup in Eq. ( 15) can be seen as a way to define a scale-adaptive partition which will prevent divergencies from showing up in the calculation of Z(a,q) for negative q values.]In the limit a--+0, one can again define the exponent r( q) from the power-law behavior of the partition function: Then, by using both the behavior of the wavelet coefficients along the maxima lines [Eq.(13)] and the extended definition of D (h) in Eq. (12), one can compute the D (h) singularity spectrum from the Legendre transform of r(q): D(h)=minq[qh --r(q)].(17) Let us point out that Eqs. ( 16) and ( 17) were rigorously derived in Ref. [46] for a class of signals whose singular part is related to the invariant measure of some affine dynamical systems (expanding Markov maps of the interval).
It is tempting to relate the exponents -r(q), defined from the wavelet-based partition function, to the exponents ~P of the structure functions [Eq.( 6)].A simple comparison of Eqs. ( 7) and ( 17) gives immediately (18) As we will discuss in Secs.Ill and IV, this relationship does not hold for every value of q.In fact, -r(q) turns out to be a more general spectrum than the ~q's, in the sense that .,-(q) is the exact Legendre transform of the D (h) singularity spectrum in most situations.

Ill. THE PROBLEM OF DIVERGENCIES IN THE COMPUTATION OF NEGATIVE q-ORDER MOMENTS
Let us first remark that in the original PF definition [Eq.( 6)], the structure functions are only defined for positive integer values of p.Even though one can extend this definition to real positive p values by replacing Bv (x, I) by I Bv (x, I) I in Eq. ( 6), the structure functions do not exist for negative p values.Indeed, there is no reason, a priori, that the probability density of Bv vanishes around Bv =0.Consequently, the Legendre transform ( 7) is not valid for p < 0 and only the part of the D (h) spectrum corresponding to the strongest singularities is amenable to the SF approach.
Recently, in the context of the study of growth processes resulting in self-affine interfaces, a slightly different "structure-function" method has been introduced in order to get rid of the divergencies problem encountered in the computation of negative-order moments [43].According to Barabasi and Vicsek, the multifractal properties of a self-affine signal f (x) can be investigated by calculating the q-order "height-height correlation function" defined as where the spatial average is performed by considering only the terms such that f(x)-f(x +l):¥=0.This restriction artificially cures the divergencies problems for q < 0 and the exponents ~q can be computed for any q.
Practically, this method consists in removing all the local increments that are below a fixed threshold (which could actually correspond to the numerical accuracy).The increment probability density functions are thus modified and set to zero over an interval whose length corresponds to the introduced threshold.However, if one chooses, as in Ref. [43], this artificial cutoff to be scale independent, one can expect to observe some self-similarity breaking when considering small increments [which are the ones that are dominating in Eq. ( 19) for negative values of q ].This phenomenon is likely to develop into a misleading phase transition [7,22,[63][64][65] in the ~q spectrum for some critical negative q value.This is the explanation of the very puzzling results obtained by Barabasi and Vicsek [43] when using their method to study fractional Brownian processes.These processes [34,66,67], generally indexed by a parameter H, are characterized by an increment probability density which is a Gaussian function whose variance behaves like z 2 H when considering increments over a distance l.If B H is a fractional Brownian process corresponding to the parameter H, since its increments BB H are Gaussian and stationary, one can compute spatial averages from ergodic formulas [68].In particular, the scaling behavior of the average in Eq. ( 19) can be derived analytically in the limit z_.. + oo (This limit is actually the one taken in Ref. [43]; it corresponds to the limit of infinite range of scales when a lower cutoff is imposed.Similar results can be derived in the limit l-+0): where the constant c is the cutoff introduced above.Then, using a straightforward change of variables, one gets One thus obtains the following exact formula for ~q: There is, thus, a singularity in the derivative of ~q for q = -1 which can be understood as a phase transition in the scaling properties of the signal.But one must emphasize again that this feature is only an artifact of the method.Even though one could try to extract from equations like Eq. ( 20) some relevant information about the singularity spectrum of the studied signals, this could be achieved only for particularly simple cases but not in more general (multifractal) situations [69].
This example clearly illustrates some intrinsic limitations in the SF method; this method cannot deal with the divergencies problems inherent to the computation of negative-order exponents without losing the natural Legendre transform bridge with the D (h) spectrum.The situation is much more transparent when one uses the WTMM method.According to the definition (15), the self-similarity properties of the signal are directly incorporated into the calculation of the partition function which never diverges whatever the value of q is.In Fig. 1 are reported the results of a statistical analysis of a fractional Brownian process ( H = +) using the WTMM method.The numerical signal was generated by filtering uniformly distributed pseudorandom noise in Fourier space in order to have the required k -513 spectral density [34].A BH=t/3 fractional Brownian signal is shown in Fig. l(a).The corresponding wavelet-transform modulus-maxima skeleton computed with the "Mexican hat" t/J< 2 l(x) [Eq.(10)] is illustrated in Fig. l(b).When plotted versus q, the exponent -r(q) extracted from the power-law behavior of Z(a,q) consistently fall on a line of slope h =0.33±0.01.Moreover, Fig. l(c) shows that the theoretical prediction 7(q) =q /3-1 (21) provides a remarkable fit of the data.No spurious phase transition is observed.By Legendre transforming 7(q) [Eq.( 17)], we obtain As expected theoretically, we find that the Brownian signal B 113 (x) is almost everywhere singular with a unique Holder exponent h =f.Throughout this simple example, we have illustrated the first crucial advantage of the WTMM method: this method is based on a discrete summation over well chosen points (scale-adaptive partition) rather than making a blind average over the whole sampling interval.Therefore divergencies are naturally removed without introducing any arbitrary artifice which might induce some meaningless phase transition phenomena in the observed spectra.

IV. LIMITATIONS IN THE RANGE OF ACCESSIBLE HOLDER EXPONENTS
In this section, we will focus our comparative analysis between the SF and the WTMM methods on the computation of both the partition and the structure functions for positive q values exclusively.Moreover, we will assume that the function f(x) under study is almost everywhere singular; this implies that the maximum value of the D (h) spectrum (which is the Hausdorff dimension of the support of the singularities) is 1.The case of signals which may be nonsingular on some set of finite Lebesgue measure will be discussed in Sec.V.  5) and (11), one can easily check that the exponent defined by the local behavior of the increments is equal to the Holder exponent if and only if the singularity lies in the first derivative of the function f.Indeed, in the definition (11), Pn(x 0 ) can be seen as the first n + 1 terms of the Taylor series off at x 0 ; it coincides with the constant term f (x 0 ) solely when f (x) is not differentiable in x 0 • When the function f is differentiable at least once, the behavior of the increment 8f(x 0 ,l) is generically dominated by the linear term j'(x 0 )/.Thus, whenever the function f has a local Holder exponent greater than l, the exponent identified by the increments will be h = 1.This is the demonstration that the SF method fails to detect the part of the D (h) spectrum [Eq.(12)] which lies beyond the value h 2:: l [even if this part corresponds to positive q values in Eq. ( 19)] and which corresponds to the weakest singularities.Furthermore, if the maximum value of D (h) is reached for some h 2:: l, the accessible range of singularities is truncated; the Legendre transform (7) becomes (23)

Therefore only the Holder exponents in the range [O,h * ],
where h* satisfies h*=l-[1-D(h*)]/D'(h*)<l, can be captured by the SF method.This limitation is also seen on the {;q spectrum which exhibits an unexpected abrupt change-of behavior for q smaller than some q *.From Eq. ( 23), one can easily check that for q >q*, one recovers the true "singular" spectrum (i.e.,

{;q=l+minh[qh -D(h)]
).But for q<q*, one observes the trivial behavior {;q =q, which actually corresponds to a singularity spectrum D (h)= h = 1.To overcome this difficulty, one is tempted to analyze the signal using a "higher-order" SF method based on the increments of the increments of the signal.This point of view will be discussed at the end of this section.
This problem in detecting singularities of Holder exponents h 2:: I does not show up when one uses the wavelet-transform analysis.According to Eq. ( 13), provided the analyzing wavelet t/J has enough vanishing moments, the Holder exponent can be correctly estimated from the local scaling behavior of the WT coefficients.From the orthogonality condition (9), one can see that any polynomial Pn(x -x 0 ) which may mask locally the singular behavior of the signal [Eq.( 11)], is canceled by the oscillations of the analyzing wavelet.Thus, within the WTMM framework, if the analyzing wavelet is well chosen, the estimate of the 7(q) spectrum [and in turn the D ( h ) spectrum] is not biased by the presence of regular behaviors {44-46].
In Fig. 2, we illustrate these considerations on a specific example.The considered signal is a recursive signal whose singularity spectrum can be computed analytically.Its construction rule involves mainly two steps.In a first step, we generate a singular measure that spreads over the whole sampling interval according to the following standard iterative rule: an interval is divided into four segments of equal length; to each of them we assign the "algebraic" weights p 1 =0.49, p 2 =0.21, p 3 = -0.21,and p 4 =0.09,respectively.In a second step, the soconstructed measure is integrated to obtain a singular signal.However, in order to obtain some Holder exponents greater than one, instead of simply integrating the measure, we perform a fractional integration of degree {J= 1. 7. The effect of a generalized integration of degree (3 (negative {J values actually correspond to generalized derivation) is to shift the D (h) singularity spectrum [70]

: D(h)--+D(h -(3).
The signal is plotted in Fig. 2(a).Some smooth parts are actually visible by a simple inspection of the signal; they correspond to Holder exponents that are greater than 1.The theoretical D (h) spectrum is represented by the continuous curve in Fig. 2(d); it extends over the interval [0.83,2.04].The results of the WTMM analysis are shown in Figs.2(b), 2(c), and 2(d).The computations were performed using the fourth derivative t/J< 4 > of the Gaussian function as analyzing wavelet.This wavelet, which has its four first moments equal to zero, is well adapted to explore the full range of Holder exponents present in the signal.In Fig. 2(b), the behavior of the partition function Z(a,q) [Eq.(15)] is displayed versus a in a "log-log" representation, for three values of q =-5, 0, and 5, respectively.From Eq. ( 16), the slope of these plots versus q defines the 7( q) spectrum.The overall results of our 7( q) measurements are reported in Fig. 2(c).The numerical data remarkably fall on the analytical 7( q) curve.This excellent agreement is confirmed when comparing, in Fig. 2(d), the numerical D (h) spectrum obtained by Legendre transforming the 7( q) data, with the theoretical prediction.The 7( q) spectrum obtained for q > 0 with the SF method [using Eq. ( 18)] is shown in Fig. 2(e) for comparison.Along the line of the above theoretical discussion [Eq.( 23)], the WTMM method provides a remarkable estimate of the 7(q) spectrum over the whole investigated range of q values, while below some critical value q = q *, the data corresponding to the SF method systematically deviate from the analytical 7(q) curve (solid line) and follow the trivial behavior 7( q) = q -1.In the range q < q *, the scaling beha vi or of the partition and the structure functions is dominated by the singularities of Holder exponents h ~ 1; the SF method is blind to these singularities which are misleadingly identified to h = 1.For q > q *, the contributions of the strongest singularities with h < 1 are dominating and the SF method gives the correct values 7(q)=minh[qh -D(h)]; in Fig. 2(e), the data from the WTMM and the SF methods converge to a unique numerical 7(q) spectrum, in good agreement with the 0.8 WTMM and SF analysis of a signal that possesses some singularities of HOlder exponents h ~ 1.(a) Graph of the signal s (x) constructed following an iterative deterministic rule (see text).(b) log2[Z(a,q)] vs log2(a) for different values of q [Eq.(15)].(c) r(q) spectrum obtained with the WTMM method when using the fourth derivative lfJ 14 l of the Gaussian function as analyzing wavelet (e).(d) D(h) spectrum from the Legendre transform of the numerical r(q) data in (c).(e) Comparison of the r(q) spectra obtained, respectively, with the WTMM (e) and the SF ( A. ) methods.The dashed curve corresponds to the straight line r(q)=q -1.In (c), (d), and (e), the solid curves represent the theoretical predictions.analytical curve.
Singularities of negative Holder exponent.A direct comparison can be performed between the wavelet analysis and the structure-function approach if one rewrites the definition of the increments in the following integral form: where a< 0 (x)=f>(x -1 )-f>(x).In this formulation, the local increment appears as a "wavelet coefficient" computed in x 0 at scale 1 with the analyzing wavelet a(!).
The above discussed inefficiencies of the increments to reveal Holder exponents greater than 1 come from the fact that only the first moment of a< 1 l is equal to zero.This special analyzing wavelet has been coined the "poor man wavelet" [66].This nickname is all the more appropriate as a 0 l suffers from another important drawback which appears when analyzing very strong singularities.To account for the strongest singularities which may exist in a signal, one needs to extend the definition (11) to tempered distributions.Roughly speaking, we will say that a distribution f has, at x 0 , a Holder exponent h (x 0 ), if its primitive (in the sense of distributions) has, at the same point, a local Holder exponent h (x 0 ) + 1 (we refer the reader to Ref. [59] for a more rigorous definition).According to this definition, the Holder exponent can take on negative values also [e.g., the Dirac distribution f>(x) can be considered as the derivative of the Heaviside function for which h (0)=0; thus, it has, at x =0, a Holder exponent h ( 0) = -1].The analyzing wavelet a 0 l is made of two Dirac distributions and therefore it cannot generally be integrated against a tempered distribution.Thus, when using Eq. ( 5) to study some distribution which has singularities of negative Holder exponent, one can expect to observe severe instabilities in the computation of the structure functions.The study of the scaling behavior of the increments of a signal that displays some discontinuities or other stronger singularities is a rather questionable procedure.
In Fig. 3 we show the results of both the WTMM and the SF analysis of a signal that possesses some singularities of negative Holder exponent.The signal is illustrated in Fig. 3(a).Its construction follows the same iterative rule as the one used to generate the multifractal recursive signal in Fig. 2(a).The model parameters are p 1 =0.36,Pz =0.24, p 3 = -0.24,p 4 =0.16; the generalized integration parameter is {3=0.4.In Fig. 3(a), one can observe several "jumps" in the signal: they correspond to singularities with negative Holder exponents.The analytical D (h) spectrum, represented by the solid line in Fig. 3(d), actually extends below the value h =0, over the range [ -0.1, 1. 36].The WTMM analysis is reported in Figs.3(b), 3(c), and 3(d).The analyzing wavelet 1/J(Zl is the second derivative of the Gaussian function.In Fig. 3(b) are shown some plots of log[ Z (a, q)] vs log( a), for different values of q; determining the exponents 7( q) just amounts to extracting the slope of these graphs from a least-squares linear regression fit.As illustrated in Fig. 3(c), the numerical data fall on a convex nonlinear curve which is particularly well fitted by the theoretical r(q) spectrum.Its Legendre transform D (h), in Fig. 3(d), is a single humped curve characteristic of the multifractality of the signal; it is also in remarkable agreement with the theoretical D (h) spectrum.The comparison of the SF and WTMM methods is reported in Fig. 3(e) where the data for r( q) obtained from both methods are compared to the analytical spectrum.As expected theoretically, the SF method gives correct values as long as q is smaller than some q value q * *, whereas for q > q * *, the SF numerical data systematically depart from the theoretical curve.The critical value q * * actually corresponds to the value of q where r(q) starts to decrease.This means that, for q > q **, negative Holder exponents are dominating the behavior of the partition function [from the properties of the Legendre transform, h corresponds to the derivative of r(q) with respect to q ].The analyzing function A.< 1 l is not suited to account for these very strong singularities and the SF method becomes unstable.On the other hand, the WTMM method, which involves smooth analyzing wavelets, does not possess this draw- ......... 3. WTMM and SF analysis of a signal that possesses singularities of negative Holder exponents.(a) Graph of the signal (see text for more details on its construction rule).(b) log2[Z(a,q)] vs log2(a) for different values of q.(c) r(q) spectrum obtained from the WTMM method using the second derivative ,p< 2 J of the Gaussian function as analyzing wavelet (e).(d) D(h) vs h from the Legendre transform of the numerical r(q) data in (c).(e) Comparison of the r(q) spectra computed, respectively, with the WTMM (e) and SF (.A.) methods.In (c), (d), and (e), the solid lines correspond to theoretical spectra.
back and provides reliable estimates on the entire range of q values.Comments on the use of higher-order SF methods.As we have just seen, the interval on which the D (h) spectrum can be determined from the SF method is not only limited in the weak singularities range (h < 1 ), but it does not extend below h =0 as well.The particular shape of A.< ll(x) restricts the investigation to the window hE [0, 1].In the same spirit, one could generalize this technique to higher-or lower-order SF approaches.For example, it can be shown that the square (or "box") function A. <Ol(x) is well appropriated to investigate the range of Holder exponents hE[ -1,0].Let us point out that the singular measures that fall under the scope of the multifractal formalism described by Eqs. ( 1) and ( 2) are precisely characterized by this range of Holder exponents.A "zero-order" SF method that uses A. <Ol as analyzing function corresponds exactly to the classical "box counting" algorithm [7,11] commonly used in the literature [22].In the same way, one could perform a "-1-order" SF method by using A.<-n, the primitive of A. <Ol; this would allow us to study the range h E [ -2, -1].Similarly, the "second-order" SF method would involve A.< 2 l(x)=B(x + 1)-2B(x +tHB(x) and would account for the range h E [ 1,2], and so on.One could even imagine a method which would consist in analyzing an experimental signal successively with A.< -ll, A.< 0 l, A. (I), A.< 2' 1, . . .; the overall spectrum is then reconstructed from the matching of all the associated spectra in their respective domains of validity.Besides the fact that such a method would be numerically cumbersome, it is not practically workable because of spurious boundary effects.For example, for A. (I), Eq. (23) shows that in the weakest singularities' direction, the D (h) spectrum is recovered only below some value h * < 1, which can be significantly smaller than one.In the opposite direction, the presence of singularities of negative Holder exponents might seriously disturb the estimate of D (h) for the strongest singularities of Holder exponent slightly larger than zero.These boundary mismatches of the successive spectra obtained using the A.<nl analyzing functions are rather difficult to control and they make the determination of the overall D (h) singularity spectrum rather uncertain.
In the previous discussions and numerical applications, we have pointed out that the excellent results obtained with the WTMM method are mainly a consequence of a good choice of the analyzing wavelet t/J.Indeed, one has essentially to choose the number of vanishing low-order moments of t/J greater than the greatest Holder exponent hmax which characterizes the weakest singularities present in the signal.But one may wonder what happens with this method when hmax = + oo, e.g., when the signal is nonsingular ( C 00 ) on a set of finite Lebesgue measure.
For the sake of simplicity, we will assume that the signal j(x)=s(x)+r(x) is a superposition of a singular part s ( x) living on a Cantor set [ s ( x) is assumed to be constant on each interval on which it is not singular] and a coo function r(x).Let us denote 7 8 (q) and D 8 (h) the multifractal spectra which characterize the singular signal s (x) alone.At each scale a 0 , the set of maxima lines .L 1 (a 0 ) off can be basically decomposed into two disjoint sets of maxima lines, L 8 (a 0 ) and .L,(a 0 ), corre- sponding to the lines created, respectively, by s (x) [and which are slightly perturbed by the presence of r (x)] and by the coo function r(x).It can be established [46] that along each line created by r(x) [E.L,(a 0 )] the wavelet coefficients behave like aN in the limit a -o, where N is the number of vanishing moments of the analyzing wavelet t/J [we suppose that N is chosen larger than the upper bound of the singularity range of s (x) ].Moreover, if the regular function r ( x) does not oscillate too much [in such a way that the number of maxima lines in .L,(a 0 ) be uniformly bounded], then the partition function defined in Eq. ( 15) splits into two parts: Z f(a,q)=Z 8 (a,q)+Zr(a,q)-a r,(ql +aqN, (25) where Z s and Z r correspond to summing over the maxi- ma lines in .Ls and .L,, respectively.From Eqs. ( 16) and (25), one can show [46] that there exists a critical value q crit < 0 such that q > qcrit =7(q)=7s(q) ' q <qcrit=7(q)=qN • (26) One thus predicts the existence of a singularity in the 7(q) spectrum.This nonanalyticity in the function 7(q) expresses the breaking of the self-similarity of the singular signal s(x) by the coo perturbation r(x).In the con- text of the thermodynamical analogy, this phenomenon defines a phase transition [7,22,[63][64][65].Below the critical value qcrit (which is the analog of the inverse of a transition temperature), one observes a regular phase whereas for q > qcrit• one switches to a multifractal phase.
Let us point out that, in contrast to the spurious phase transition previously obtained with the modified SF method [Eq.(20)], the phenomenon observed here contains fundamental information on the signal considered.Equation ( 26) indicates that the 7( q) spectrum in the "Coo phase" is governed by the number N of vanishing moments of the analyzing wavelet.Therefore the experimental observation of the variation of 7(q) according to Eq. ( 26) when N is changed, would be a decisive test of the presence of a highly regular part in the signal.
In Fig. 4, we illustrate our purpose for the specific example of a signal f (x) which combines a Coo component r(x)=R sin(81Tx) and a singular component s(x) which is the distribution function of a singular measure 11= s(x)= J~d/1 (this kind offunction is usually called "devil staircase").The measure 11 is actually a Bernoulli measure lying on the triadic Cantor set.It is obtained by iteratively assigning the weights p 1 =0.6 and P2 =0.4 to the two intervals left at each step of the Cantor set construction process.The signal is shown in Fig. 4(a).The spectra 7 8 (q)=7(q) and D 8 (h)=j(a=h) of the singular part s(x) are given by the 7(q) and /(a) spectra of the measure 11• as conventionally defined by Eqs. ( 1) and ( 2) in the context of the multifractal formalism for singular measures.The results of the WTMM measurement of the 7( q) and D ( h ) spectra of the signal f ( x) are shown, respectively, in Figs.4(b) and 4(c).Two different analyzing wavelets, namely, the first (t/J 0 l) and the second (t/!( 2 )) derivative of the Gaussian function [Eq.(10)] were used to compute the wavelet transform.For q > 0, both analyzing wavelets lead to numerically identical estimates for 7(q).On the contrary, for q <0, the numerical data obtained with t/J(I) and t/J( 2 ) separate from each other into two distinct straight lines of respective slope 1 and 2. We thus observe numerically the phase transition predicted by the above theoretical speculations.For positive q values, we recover the 7(q) spectrum of the underlying singular measure 11• while for q below some critical negative value qcrit> the shape of 7(q) is governed by the number N of vanishing moments of the analyzing wavelet: FIG. 4. WTMM analysis of a signal which is nonsingular on some intervals.(a) Graph of the signal f(x)=s(x)+r(x), with r (x)= R sin( 81rx) and s (x) is the distribution function of a Ber- noulli measure supported by the triadic Cantor set (see text).(b) r( q) vs q as obtained with the first [ ( o ) and ( .._ ) ] , the second [( o) and (e)], and the fourth [( o) and <•)] derivatives of the Gaussian function as analyzing wavelet; the solid lines correspond to the theoretical predictions [Eq.( 26)]; the dashed line is the part (q <qcrit) of the spectrum r,(q).(c) D(h) vs h from the Legendre transform of r(q); the symbols are the same as in (b).r(q}=Nq (N = 1 for t/J(l) and N =2 for t/J( 2 >).By Legendre transforming r(q) one gets the D (h) singularity spectrum of j(x).As shown in Fig. 4(c), as long as q > qcrit(N), the numerical results obtained with both analyzing wavelets remarkably fall on the theoretical D(h) curve (solid line).For q ~qcrit(N), however, the Legendre transform of the self-similarity breaking linear behavior of r( q) produces a linear behavior fall off of the D (h) curve towards the limiting value h = 1 for t/J( 1 > and h =2 for t/J( 2 l (actually h =N for t/J(Nl) where D (h) vanishes.This linear part is tangent to the theoretical D (h) spectrum (dashed line) and has a slope equal to qcrit(N).This is the signature of the phase transition phenomenon described above.
Remark.We have shown that a coo component superimposed on a signal that is not singular everywhere manifests in a phase transition phenomenon that masks the weakest singularities.However, since the wavelet coefficients behave like aN along the maxima lines created by the Coo function, by choosing N large enough and/or choosing a numerical threshold below which any local maximum is not considered, one can remove all the Coo maxima lines in .L,(a 0 ) [cf. the definition ( 14) of .L(a 0 ) in Sec.11] and thus "numerically restore" the self-similarity of s(x).The whole r 5 (q) and D 5 (h) spectra can then be estimated.
To show that this procedure is actually operational, we have reproduced the WTMM analysis on the same signal but with the fourth derivative of the Gaussian function [ t/J( 4 l(x)] as analyzing wavelet.The faster decrease of the wavelet coefficients along the maxima lines of .L, (T"'[f](.,a)l 1 e.L -a 4 ) makes more efficient the thresh-' old discrimination of the maxima lines emanating from the singular part s(x).The so-obtained r(q) spectrum is shown in Fig. 4(b).Now the theoretical spectrum of the singular measure is recovered and no phase transition phenomenon is observed.Let us point out that the choice of such a threshold (or analyzing wavelet) is somewhat uncertain and strongly depends on various parameters such as the number of sampling points, the relative amplitudes of r (x) and s (x) in the signal, etc.Indeed, a more reliably way to proceed consists in choosing N large enough [as compared to the largest Holder exponent of s (x)] so that the maxima lines induced by the regular part of the signal become easily distinguishable by the anomalously stiff decrease of the wavelet coefficients on the range of scale used to estimate the scaling exponents r(q) (supa~a 0 T1,b[f](.,a)lle.L, =O(a~)).

VI. CONCLUSION
Most previous characterizations of fractal signals have relied upon the measurement of the structure-function scaling exponents.In this paper, we have shown that these exponents reflect only a small part of the singular properties of these signals.We have reviewed the severe limitations to the increment-based technique that make the estimate of these exponents quite uncertain.We have proposed an alternative method, based on wavelettransform modulus-maxima tracking, that gives direct access to the D (h) singularity spectrum of any fractal signal.From the direct comparison of the WTMM method with the SF method for specific examples, we have demonstrated that the former does not introduce any bias in the estimate of the scaling exponents of some partition functions which are at the heart of these statistical analyses of multifractal signals.The efficiencies of the WTMM method originate from two main ingredients; on the one hand, the partition functions are based on discrete scaledependent summations: the skeleton of the wavelet trans- form defined by the local maxima of its modulus provides a practical guide to achieve a scale-adaptive partition; on the other hand, one uses sufficiently smooth and localized analyzing wavelets, with an arbitrary large enough number of vanishing moments, which makes the entire range of singularities accessible to this method, even in the presence of regular behavior in the signal.Preliminary results reported in Refs.[44,45] of an analysis of a fully developed turbulent velocity signal show that this method is readily applicable to experimental situations.Applications of the WTMM approach to turbulent dynamics generated by fractal growth phenomena, critical fluctuations in colloidal systems, and DNA "walks" nucleotide sequences are currently in progress.We believe that this method of determination of the singularity spectrum of fractal signals is likely to become as useful as the well-known phase-portrait reconstruction, Poincare section, and first-return map techniques for the analysis of chaotic time series.
FIG. l.WTMM and SF analysis of fractional Brownian motion.(a) Graph B 113 ( x) of a realization of a fractional Brownian process indexed by H = {.(b) Wavelet-transform maxima lines of B 113 (x).(c) r(q) vs q as computed with the WTMM method (e) and the SF method(6 ).(d) D(h) vs h from the Legend re transform [Eq.(17)] of the r( q) spectrum computed with the WTMM method.The analyzing wavelet is the second derivative 1(!< 2 l of the Gaussian function.

h > 1 .
If one compares Eqs. ( FIG.2.WTMM and SF analysis of a signal that possesses some singularities of HOlder exponents h ~ 1.(a) Graph of the signal s (x) constructed following an iterative deterministic rule (see text).(b) log2[Z(a,q)] vs log2(a) for different values of q [Eq.(15)].(c) r(q) spectrum obtained with the WTMM method FIG.3.WTMM and SF analysis of a signal that possesses singularities of negative Holder exponents.(a) Graph of the signal (see text for more details on its construction rule).(b) log2[Z(a,q)] vs log2(a) for different values of q.(c) r(q) spectrum obtained from the WTMM method using the second