Kenneth J Falconer

Jacques Lévy

Horizons of fractional Brownian surfaces

Introduction and background

Fractional Brownian random surfaces have found several applications in recent years in such diverse elds as computer modeling of landscapes [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF][START_REF] Peitgen | The Science of Fractal Images[END_REF], modeling of rough surfaces in physics and chemistry and engineering and in medical imaging. The parameters for such random elds need to be chosen so as to obtain the best t to the data at hand. As an example, natural-looking landscapes typically have local H older exponent about 0.85 almost everywhere and fractal dimension 2.15. In some situations, an additional feature has to be taken into account: the `horizon' de ned by the surface may be relevant to modeling, and one needs to control its fractal characteristics [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF][START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF]. For instance, it is believed that some glasses may be well modeled by fractional Brownian surfaces. Several properties of the glass, relevant both for industrial and domestic use, depend on the roughness of its horizon, which governs its behaviour under skimming lighting. More generally, fractal electromagnetism studies how waves interact with `fractal' media (such as the sea in the case of radar imaging) (see [START_REF] Jaggard | Fractal Electrodynamics and Modeling[END_REF]. The horizon is then the appropriate geometrical object to consider for studying skimming waves. In imaging processes, especially medical imaging, one observes a two dimensional `horizon' of a three dimensional object, such as a lung, and tries to infer the properties of the original object [START_REF] Lundahl | Fractional Brownian motion : a maximum likelihood estimator and its application to image texture[END_REF]. Permeability elds of a porous medium provide another example (see [START_REF] Addison | Engineering applications for fractional Brownian motion: self-a ne and self-similar random processes[END_REF].

It is conjectured that, for a `su ciently random and homogeneous' surface, the H older exponent of a horizon equals that of the surface, and that the `fractal dimension' of the horizon is one less than that of the surface. In this paper we investigate this horizon conjecture for index-fractional Brownian surfaces. We show that the horizon conjecture holds when = 1 2 , and for other values of ; 0 < < 1 given a highly plausible assumption (the `maximum property') on the form near a maximum of index-fractional Brownian motion in one variable. We present direct and indirect computational evidence based on very large data sets to support the conjecture.

We rst consider a function x : 0; 1] d ! R, where d is a positive integer, with graph graphx f(t; x(t)) : t 2 0; 1] d g R d+1 . In particular, for d = 2 we think of graphx as a surface above the unit square S = 0; 1] 0; 1] and we de ne the horizon function z : 0; 1] ! R by z(t) = sup 0 u 1 x(t; u);

(1.1) where t = (t; u) in coordinate form. Thus graphz R 2 may be thought of as the horizon of the surface graphx R 3 . We are interested in the relationships between the fractality of the functions x and z and especially between their H older exponents and the dimensions of the surface, graphx, and its horizon, graphz. Whilst our main concern in this paper is with index-fractional Brownian surfaces, we rst gives some lemmas relating to functions of a general form. We de ne the (lower) H older exponent of x at t by H x (t) = lim inf v!t log jx(v) x(t)j log jv tj :

(1.2)

It is also useful to de ne the lower and upper approximate H older exponents by H ap x (t) = ap lim inf v!t log jx(v) x(t)j log jv tj ; H ap x (t) = ap lim sup v!t log jx(v) x(t)j log jv tj : (1.3) Recall that the lower approximate limit is de ned by ap lim inf v!t f(v) = a if there exists a set A of Lebesgue density 1 at t such that lim inf v!t;v2A f(v) = a, with a similar de nition involving limsup for the upper approximate limit. (The set A has Lebesgue density 1 at t if it is Lebesgue measurable and lim !0 L d (B(t; ) \ A) L d (B(t; )) = 1 where L d is d-dimensional Lebesgue measure on R d ). Clearly H x (t) H ap x (t) H ap x (t) If H ap x (t) = H ap x (t) we say that the approximate H older exponent exists and write H ap x (t)

for the common value. We write dim H ; dim B ; dim B for Hausdor , lower box and upper box dimension respectively (see [START_REF] Falconer | Fractal Geometry -Mathematical Foundations and Applications[END_REF] for their de nitions and basic properties.

There are some basic lemmas that are useful in calculating dimensions of graphs. The rst lemma provides an upper bound for the dimension of graphx, given a H older condition on x.

Lemma 1.1 Suppose that x : 0; 1] d ! R satis es jx(v) x(t)j cjv tj

(1.4) for all t; v 2 0; 1] d , where 0 < 1 and c > 0. Then H x (t) for all t 2 0; 1] d , and dim H graphx dim B graphx dim B graphx d + 1 .

Proof. See (Falconer 1990, Corollary 11.2). The next lemma is useful for obtaining a lower bound for the dimension of a graph.

Lemma 1.2 Let x : 0; 1] d ! R. Suppose that there is a number s (1 s 2) such that for L d -almost all t 2 0; 1] d there exists c > 0 such that L d fv : jv tj h and jx(v) x(t)j hg ch s for all su ciently small h. Then dim H graphx s.

Proof. We may de ne a probability measure on graphx by (A) = L d ft 2 (There is a rather weaker condition that leads to a similar lower bound for the boxdimensions, (see Falconer 1990 Corollary 11.2).) The following useful lemma shows that a H older condition on a surface implies a corresponding condition on its horizon, with corresponding bounds for the dimensions.

Lemma 1.3 Suppose that x : 0; 1] d ! R satis es jx(v) x(t)j cjv tj

(1.5) for all t; v 2 S, where 0 < 1 and c > 0. Then the horizon function z given by (1.1) satis es a H older condition jz(v) z(t)j cjv tj

(1.6) for all t; v 2 0; 1], so in particular dim H graphz dim B graphz dim B graphz 2 and H z (t) for all t 2 0; 1].

Proof. Fix t and v, and choose u such that z(t) = x(t; u) (since u is continuous the supremum is attained). Then z(v) x(v; u) x(t; u) cjv tj = z(t) cjv tj ; so z(v) z(t) cjv tj . Combining this with the symmetric inequality gives (1.6). The

Hausdor and box dimension estimates follow from Lemma 1.1, and the H older exponent from (1.2).

Lower bounds for dimensions of horizons are much more awkward to obtain, and we proceed to consider this problem for fractional Brownian surfaces.

(2.1) (see [START_REF] Adler | The Geometry of Random Fields[END_REF]. It follows easily that E(jX(t 1 ) X(t 2 )j 2 ) = jt 1 t 2 j 2 ;

(2.2)

and that E(X(t) 2 ) = jtj 2 , so in particular X(0) = 0 a.s. Thus X has isotropic stationary increments and is scale invariant, in the sense that for h 6 = 0 the Gaussian process jhj X(ht) has the same covariance functions and thus the same distributions as X(t).

It should be emphasised that the increments are not independent in any sense unless = 1 2 . Index-fBf has a continuous version (see [START_REF] Adler | The Geometry of Random Fields[END_REF]; indeed, given > 0, with probability one there is a uniform H older inequality jX(t 1 ) X(t 2 )j jt 1 t 2 j

(2.3) for all jt 1 t 2 j su ciently small. We will always work with this continuous version of the process. It is not di cult to show, using Lemmas 1.1 and 1.2, that dim H graphX = dim B graphX = dim B graphX = d + 1 a.s. [START_REF] Adler | The Geometry of Random Fields[END_REF][START_REF] Falconer | The horizon problem for random surfaces[END_REF].

There is an extremely useful way of regarding Gaussian elds in terms of Hilbert spaces. Let be the sample space underlying a Gaussian eld on 0; 1] d and write L 2 ( ) for the space of zero mean random variables X on with E(X 2 ) < 1. Then L 2 ( ) is a real Hilbert space under the inner product given by E(X 1 X 2 ) for X 1 ; X 2 2 L 2 ( ).

We write H 0 for the subspace spanfX(t) : t 2 0; 1] d g and H for the closure of H 0 in L 2 ( ). With this Hilbert space approach, the covariances E(X(t 1 )X(t 2 )) are just inner products. This provides a very convenient way of representing conditional variances of the X(t). For example, if A S, then var(X(t)jX(s) : s 2 A) is the Hilbert space distance from X(t) to the subspace K spanfX(s) : s 2 Ag of H, and there is a decomposition X(t) = X 1 (t) + X 2 (t), where X 1 (t) 2 K and X 2 (t) 2 K ? . To aid calculation there is a useful integral representation of the covariances of index-fBf on 0; 1] d . It is easy to check [START_REF] Geman | Occupation densities[END_REF]Kahane 1985) that (2.1) may be written as

E(X(t)X(v)) = k Z R d jaj d 2 (e it a 1
)(e iv a 1)da;

(2.4)

where k depends only on d and . In particular this gives that, for t 1 ; : : : ; t n 2 R d and 1 ; : : :

; n 2 R, E(j n X j=1 j X(t j )j 2 ) = k Z R d jaj d 2 j n X
j=1 j e it j a j 2 da;

(2.5) provided that P n j=1 j = 0. This integral form provides a very useful way of calculating conditional variances of increments of the X(t). For E an ellipsoid in R d and > 1, write E for the homothetic (similar and similarly situated) ellipsoid with the same centre as E obtained by scaling E by a factor about its centre. For T R d we write F T for the sigma eld (X(t 1 ) X(t 2 ) : t 1 ; t 2 2 T) generated by the di erences of X at pairs of points in T.

Proposition 2.1 Let E R d be an ellipsoid with semi-axis lengths 0 < s 1 : : : s d and let E be de ned as above for some > 1. If V E and T E c (the complement of E ), and v 2 V and t 2 T, then var var(X(v) X(t) j (F T ; F V )) cs d+2 1 (s 1 : : : s d ) 1 ;

(2.6) and hence, for all 1 < 2 , P( 1 < X(v) X(t) 2 j (F T ; F V )) c 0 j 2 1 j s d=2 1 (s 1 : : : s d ) 1=2 ;

(2.7)

where c and c 0 depend only on d; and .

Proof. Working in the Hilbert space setting, var is the square of the Hilbert space distance of X(v) X(t) from the subspace de ned by the conditioning; thus

var = inf E(jX(v) X(t) + m X j=1 j (X(v j ) X(v)) + m 0 X j=1 j (X(t j ) X(t)) : v j 2 V; t j 2 Tj 2 )
where the in mum is over all 1 ; : : : ; m ; 1 ; : : : ; m 0 2 R, v 1 ; : : : ; v m 2 V and t 1 ; : : : ; t m 0 2 V . By rede ning the j ; j by taking appropriate linear combinations in the sums, we get that var inf E(j X j j X(v j ) X j j X(t j ) : v j 2 V; t j 2 T; X j j = X j j = 1j 2 ):

Since P j j P j j = 0, the integral representation (2.5) gives var inf k ( Z j P j j e iv j a P j j e it j a j 2 jaj d+2 da : v j 2 V; t j 2 T; We take g to be a C 1 `bump function', with g(w) = 1 for w 2 B(0; 1) and g(w) = 0 for w 6 2 B(0; ). Let p 2 R d be the centre of the ellipsoid E and let be a non-singular self-adjoint linear operator on R d that maps the unit ball onto the translated ellipsoid E p; thus E = (B(0; 1)) + p and E = (B(0; )) + p. With P j j = P j j = 1, it follows, using the de nitions of and g, that

1 = X j j g( 1 (v j p)) X j j g( 1 (t j p)) = (2 ) d Z ( X j j e ib 1 (p v j ) X j j e ib 1 (p t j ) )ĝ(b)db = (2 ) d Z ( X j j e ia (p v j ) X
j j e ia (p t j ) )ĝ( (a))j det jda;

writing b = (a) and using the self-adjointness of . Thus by Cauchy's inequality, 1

Z je ia p j 2 j P j j e ia v j P j j e ia t j j 2 jaj d+2 da Z ĝ( (a)) 2 jaj d+2 j det j 2 da

(2 ) 2d Z j P j j e ia v j P j j e ia t j j 2 jaj d+2 da Z ĝ(b) 2 j 1 (b)j d+2 dbj det j

(2 ) 2d Z j P j j e ia v j P j j e ia t j j 2 jaj d+2 da Z ĝ(b) 2 jbj d+2 db s d 2 1 (s 1 : : : s d ); as det = s 1 : : : s d and k 1 k = s 1 1 . Since ĝ(b) is continuous and rapidly decreasing, the last integral is nite, so combining this estimate with (2.8) gives (2.6). Inequality (2.7) follows on noting that, if W is Gaussian with mean and variance 2 , then

P( 1 < W 2 ) = (2 ) 1=2 1 Z 2 1 exp( (s )=2 2 )ds (2 ) 1=2 1 j 2 1 j:
3 The maximum property for fBm

In this section we brie y consider indexfractional Brownian motion (fBm) X : 0; 1] ! R, that is the index-fBf de ned by (2.1) with d = 1. Then X has a continuous version which a.s. attains its maximum at a unique point of (0; 1). We say that index-fBm X on 0; 1] satis es the maximum property if, conditional on X attaining its maximum at t m 2 (0; 1), for all > 0 there exists a.s. > 0 such that X(t) X(t m ) jt m tj + for all t 2 0; 1] with jt m tj :

(3.1)
Of course, the opposite inequality, that for all > 0 X(t) X(t m ) jt m tj for all t close to t m , follows from the a.s. H older condition (2.3). Thus X has the maximum property if its sample paths fall away from their maximum at a rate with power law exponent close to . It should be emphasised that the maximum property is de ned in terms of fractional Brownian motion on the real interval and not in terms of fractional Brownian elds on R2 , although this is where our applications lie. Proposition 3.1 Usual Brownian motion, that is index-1

Proof. This fact, stated for example in [START_REF] Barlow | Levels at which every Brownian excursion is exceptional[END_REF], follows from the decomposition theorem [START_REF] Itô | Di usion Processes and their Sample Paths[END_REF][START_REF] Williams | Path decomposition and continuity of local times for one-dimensional di usions I[END_REF]) which identi es Brownian motion conditioned on a maximum at 0 with the 3-dimensional Bessel process (that is the radial component of Brownian motion in R 3 ).

The following maximum conjecture seems extremely likely, but it does not appear to have been established rigorously.

Conjecture 3.2 For all 0 < < 1, index-fBm has the maximum property. This conjecture holds in an `approximate' sense: conditional on X attaining its maximum at t, we have a.s. that for all > 0 ap lim sup t!t (X(t) X(t m ) + jt m tj + ) 0;

(see Geman & Horowitz 1980, Theorem 30.4), but this is not su cient for our arguments.

Probabilistic estimates

From now on, we take d = 2 and work with the index-fractional Brownian eld X : S ! R where S is the unit square 0; 1] 0; 1]. In this case we refer to graphX as the indexfractional Brownian surface. In this section we estimate some conditional variances and probabilities of X preparatory to our calculations on the fractality of the horizon of X.

We use the vector and coordinate forms t = (t; u); v = (v; w); t 1 = (t 1 ; u 1 ) etc. for points in R 2 . We introduce sigma-elds generated by the random variables X(t; u) on certain subsets of the square S. For t 2 0; 1] let F t = (X(t; u 1 ) X(t; u 2 ) : u 1 ; u 2 2 0; 1]), and for (v; u) 2 S and r > 0 let F (v;u);r = (X(v; w) X(v; u) : w 2 0; 1]; jw uj r). Thus F t determines X along a section with rst coordinate t, and F (v;u);r determines X along a section of length 2r with midpoint (v; u), in both cases to within a vertical displacement.

Proposition 4.1 For all r > 0 and all pairs of points (t; u); (v; u) 2 S with 0 < jv tj r, we have, almost surely, var(X(v; u) X(t; u) j (F t ; F (v;u);r )) cjv tj 1+2 r 1 ;

(4.1) and for all 1 < 2 ,

P( 1 < X(v; u) X(t; u) 2 j (F t ; F (v;u);r )) c 0 j 2 1 jjv tj 1=2 r 1=2 ; (4.2)
where c and c 0 depend only on d and . Proof. We apply Proposition 2.1, taking E to be the ellipse with centre (u; v) and semiaxes parallel to the coordinate axes, of lengths 1 3 jt vj and r, and taking = 2. Then E contains the interval (v; u r); (v; u+r)] and E is disjoint from the line f(t; w) : t 2 Rg, so the result is immediate from Proposition 2.1

The (random) horizon function of the index-fBf X(t; u) is given by Z

(t) = sup 0 u 1 X(t; u) (0 t 1):
Since X is almost surely continuous, for each t 2 0; 1] we can de ne the random variable U t by U t = inffu 2 0; 1] : X(u; t) = Z(t)g;

(4.3) so U t is F t -measurable. (In fact X(t; u) a.s. attains a unique maximum at u = U t 2 (0; 1), with Z(t) = X(t; U t ), see Section 5.) For r > 0 and t; v 2 0; 1], we de ne the random variable Z r t (v) = supfX(v; u) : u 2 0; 1] and ju U t j rg:

(4.4)

Then Z r t (v) gives the maximum of X(v; u) over a restricted interval of u close to U t .

Proposition 4.2 For all r > 0, 1 < 2 and t; v 2 0; 1] with 0 < jv tj r, we have

P( 1 < Z r t (v) Z(t) 2 ) c 0 j 2 1 jjv tj 1=2 r 1=2 ; (4.5)
where c 0 depends only on .

Proof. For (v; u) 2 S de ne the random variable Z(v; u) = supfX(v; w) : jw uj r and (v; w) 2 Sg: Since Z(t) X(t; u) and Z(v; u) X(v; u) are F t -measurable and F (v;u);r -measurable respectively, and therefore are both (F t ; F (v;u);r )-measurable, Proposition 4.1 implies that for all t 2 0; 1] and (v; u) 2 S, P( 1 < Z(v; u) Z(t) 2 j (F t ; F (v;u);r )) c 0 j 2 1 jjv tj 1=2 r 1=2 ; a.s.. For n = 1; 2; : : : de ne the random variables U n t = 1 n nU t ] where ] denotes `the integer part of'. Then U n t is F t -measurable and nite valued and gives a discrete approximation to U t (we need to invoke this discretisation for step (4.6) below). Since (U n t ) (F t ; F (v;u);r ), the tower property for conditional probabilities gives, for each i = 1; : : : ; n,

P( 1 < Z(v; i n ) Z(t) 2 j U n t ) c 0 j 2 1 jjv tj 1=2 r 1=2 : Then P( 1 < Z(v; U n t ) Z(t) 2 ) = n 1 X i=0 P( 1 < Z(v; i n ) Z(t) 2 j U n t = i n )P(U n t = i n ) c 0 j 2 1 jjv tj 1=2 r 1=2 : (4.6)
Since X is a.s. continuous, Z(v; U n t ) Z(t) converges a.s. to Z(v; U t ) Z(t) as n ! 1, and hence converges in distribution, so (4.5) follows from (4.6).

Fix > 0 and de ne, for t; v 2 0; 1], Z t (v) = supfX(v; u) : u 2 0; 1] and ju U t j jv tj 1 g: 4.4) taking r = jv tj 1 jv tj. Corollary 4.3 For all t 2 0; 1] and 0 < < 1 , we have a.s. that L(v 2 0; 1] : jv tj h and jZ t (v) Z(t)j h) h 2 (4.8) and ap lim sup v!t log jZ t (v) Z(t)j log jv tj (4.9) for all su ciently small h.

(4.7) Thus Z t (v) is just Z r t (v) of (
Proof. Setting r = jv tj 1 in Proposition 4.2 and replacing Z r t by Z t , noting (4.4) and (4.7), gives

P( 1 < Z t (v) Z(t) 2 ) c 0 j 2 1 jjv tj =2 : (4.10) For t 2 0; 1], E(L(v : jv tj h and jZ t (v) Z(t)j h)) Z t+h t h P(jZ t (v) Z(t)j h)dv 2c 0 h Z t+h t h jv tj =2 dv = 4c 0 (1 =2) 1 h 2 =2 : (4.11)
Thus for each t 2 0; 1], by applying (4.11) with h = 2 k (k = 1; 2; : : :) and using the Borel-Cantelli lemma, we get a.s. that (4.8) holds for all su ciently small h. Again using (4.10), we have E(L(v : jv tj h and jZ t (v) Z(t)j jv tj + ))

Z t+h t h P(jZ t (v) Z(t)j jv tj + )dv 2c 0 Z t+h t h jv tj =2 dv = 4c 0 (1 + =2) 1 h 1+ =2 ;
so another application of the Borel-Cantelli lemma gives 1 2h L(v 2 0; 1] : jv tj h and jZ t (v) Z(t)j jv tj + ) h =4 a.s. for su ciently small h, from which (4.9) follows using (1.3).

Horizons of fractional Brownian surfaces

Recall that the horizon of the index-fBf X on S = 0; 1] 0; 1] is the random function Z on 0; 1] given by Z(t) = sup 0 u 1 X(t; u):

(5.1) Since X satis es the H older condition (2.3) a.s., Lemma 1.3 implies that a.s. H Z (t) for all t 2 0; 1] and that dim H graphZ dim B graphZ dim B graphZ 2 .

Estimates in the other direction are more awkward. It was shown in [START_REF] Falconer | The horizon problem for random surfaces[END_REF] that H Z (t)

+ 1 2 and 1 1 2 dim H graphZ almost surely. Here we show that dim H graphZ = dim B graphZ = dim B graphZ = 2 almost surely if = 1 2 , and also for all 0 < < 1 if we assume the maximum Conjecture 3.2 for index-fractional Brownian motion on the line. Our strategy is to show that, for each t, we have Z(v) = Z t (v), for all v su ciently close to t, see (4.4), and then apply Lemma 1.2 using the estimate of Corollary 4.3 to obtain the dimensions and H older exponents. With U t the value of u at which the maximum sup 0 u 1 X(t; u) = Z(t) is attained (4.3), we consider how U t varies with t.

Proposition 5.1 Assume that index-fBm has the maximum property. Let > 0. Then for all t 2 0; 1], there exists a.s. > 0 such that, for all v 2 0; 1] with jv tj < , jU v U t j jv tj 1 :

(5.2)

Proof. Choose 0 < < such that 1 < ( )=( + ):

(5.3) Index-fBf satis es an a.s. uniform H older condition (2.3) which may be invoked in the form: there exists a.s. a number 1 > 0 such that jX(t; u) X(v; w)j 1 3 j(t; u) (v; w)j

(5.4) for all (t; u); (v; w) 2 S with j(t; u) (v; w)j 1 . Fix t 2 0; 1] for the remainder of the proof. The process X(t; ) X(t; 0) : 0; 1] ! R is index-fBm, being a section of index-

(t) = 1
2 for almost all t 2 0; 1] and dim H graphZ = dim B graphZ = dim B graphZ = 1 1 2 . Proof. This is immediate from Theorems 5.3 and 5.4, since by Proposition 3.1 index-1 2 fBm has the maximum property.

Synthesis of fractional Brownian surfaces

The remainder of this paper is devoted to numerical investigation of the maximum conjecture and the fractal characteristics of the horizon function (box dimension and H older exponent). To this end, we have synthesized a large number of traces of fBm-s and fBf-s. In this section we describe the methods used for generating sample paths for various values of , and present some graphs. Section 7 explains the procedures for investigating the maximum property. Sections 8 and 9 are devoted to the estimation of the box dimension and H older exponent of the horizon Z.

Several methods have been proposed for synthesizing fractional Brownian motion, and, more generally, fractional Brownian elds. They include Choleski decomposition, midpoint displacement and its various improvements [START_REF] Daniel | Modeling and estimation of fractional Brownian motion using multiresolution stochastic processes[END_REF], spectral synthesis [START_REF] Kaplan | An improved method for 2-D self-similar image synthesis[END_REF], wavelet-based methods, and synthesis based on di erential models [START_REF] Kesner | 1/f noises[END_REF]. The reason why so many algorithms have appeared is that synthesizing fBm-s is by no means an easy process, especially if one needs to build large traces. The problem lies mainly in the non-Markovian nature of fBm, the strong correlations that it displays (for < 1=2, each increment is negatively correlated with all the others, while for > 1=2 the spectral density of the increments is singular at the origin).

Though the Choleski method allows exact synthesis, a plain implementation su ers from two drawbacks: a large amount of computing time and a large amount of memory are needed to construct the signal. Under some minor restrictions on the generated traces, it is possible to use fast and e cient algorithms for Choleski decomposition, but various approximate methods using di erent approaches have also been designed that allow reasonable computer time/memory requirements.

We rst describe in some detail the principal steps involved in the classical Choleski decomposition method. Assume we wish to generate samples of an index-fBm X at N equidistant points of 0; 1]. Let DX denote the discrete increments of X, that is DX(k=N) = X(k=N) X((k 1)=N); k = 1; : : : N. These increments form a discrete (strict sense) stationary Gaussian process with zero mean, and the statistical properties of the vector DX N = (DX(1=N); DX(2=N); : : :; DX(1)) are entirely determinated by the autocovariance matrix A N = E(DX N (DX N ) T ), where U T denotes the transpose of U. From (2.1), E(DX(i=N)DX((i + k)=N)) = 1 2N 2 (jk + 1j 2 + jk 1j 2 2jkj 2 ). Since A N is positive de nite, it may be written using its Choleski decomposition as:

A N = L N L T N (6.1)

where L N is an invertible lower triangular matrix. Let DY N = (DY (1=N); DY (2=N); : : : ; DY (1)) be an N samples realization of a unit variance centered white Gaussian noise. It is easy to see that the autocovariance matrix of the random vector L N DY N is exactly A N . Indeed:

E(L N DY N (L N DY N ) T ) = L N E(DY N (DY N ) T )L T N = A N (6.2)
We may thus set DX N = L N DY N , and generate a realization of the index-fBm X as X(k=N) = P k p=1 DX(p=N). Since A N depends only on N and , it is entirely determined once we have xed the exponent and the number of points we wish to generate. The problem of synthetizing a sample of an fBm is thus reduced to that of computing L N from A N . Note that, so far, we have only used the fact that A N is a valid autocovariance matrix, so that it has a Choleski decomposition. Thus the procedure above may be applied for synthetizing any discrete Gaussian process. However, a direct method for general Choleski factorization has complexity O(N 3 ) and requires O(N 2 ) memory. This precludes the use of this approach for building large traces. Fortunately, when the process is stationary (this is why we work with DX rather than X) and in the common case where the samples are equi-spaced, the matrix considered is Toeplitz: one can then use fast algorithms, such as the Schur or Levinson algorithms, which have complexity O(N 2 ) and need O(N) memory. It is possible to do even better if one forces N to be a power of 2. Since such a requirement should not a ect the veri cation of the maximum property conjecture, we shall restrict ourselves to this case. Then, the doubling Schur algorithm [START_REF] Ammar | Superfast solution of real positive de nitive Toeplitz systems[END_REF] allows the complexity to be reduced to O(N(log 2 (N)) 2 ). This is very reasonable and permits the synthesis of quite large traces. We now turn to the oldest of the approximate algorithms, the midpoint displacement. In the case of Brownian motion, this is in fact the original construction by Paul L evy of the Wiener process. This method is very e ective, as its complexity is linear in the number of samples. When di ers from 1/2, however, it is no longer exact, and the resulting process, while still Gaussian, has second order properties that di er signi cantly from those of the fBm. In [START_REF] Daniel | Modeling and estimation of fractional Brownian motion using multiresolution stochastic processes[END_REF], an improvement of the classical scheme is proposed, which allows one to recover approximately the right covariance function with a low computational burden. The algorithm is based on the notion of a multiresolution process, and improves on the classical method by using statistical descriptions of both the interpolation and displacement steps.

Wavelet based methods rely mainly on the fact that the wavelet transform acts as a `whitening lter' on fractional Brownian motion. This allows easy synthesis of the wavelet coe cients of fBm. However, the problem of building the low frequency approximation of the signal remains, and, while this methods allow synthesis of large datasets, it does not appear to be exact enough for our purposes.

An interesting spectral method for generating an approximation to a two dimensional N N fBf is the incremental Fourier synthesis described in [START_REF] Kaplan | An improved method for 2-D self-similar image synthesis[END_REF]. The idea is to create rst a periodic random eld of size 2N 2N with statistics close to those of the increments of the fBf over half the spatial period. Periodic random elds are easy to generate because their Karhuenen-Lo eve transform is simply the two dimensional discrete Fourier transform. The approximate fBf is then obtained by adding up the increments. Thanks to the use of the fast Fourier transform, the complexity of this method is only O(N 2 log 2 (N)) for a N N image, with a memory requirement of O(N 2 ). An alternative spectral method that gives an approximation to fBf but has complexity O(N log N) is described in [START_REF] Yin | New Methods for Simulation of Fractional Brownian Motion[END_REF]. Finally, while methods based on di erential models (that is where the process is obtained as a solution of a stochastic di erential equation) are interesting because they have a physical meaning, they do not seem to yield correct approximations to fBm. In fact, these methods are mainly used to study generic `1/f' noises.

Since we wish to test numerically a conjecture related to a very ne property of fractional Brownian motion, it is important to ensure that the synthesized traces are as close as possible to real sample paths of fBm-s. Thus we have chosen to use the only exact method, i.e. the Choleski decomposition via the doubling Schur algorithm. This implementation allowed us to generate traces of fBm with 131072 sample points and 2048 2048 fBf. More speci cally, the following data were synthesized: veri cation of the maximum property: for each value of = 0:1; 0:2; : : : ; 0:9, 100 traces of fBm with 131072 sample points were generated. estimation of the box dimension and H older exponent of the horizon function: for each value of = 0:1; 0:2; : : : ; 0:9, 100 traces of fBf with 2048 2048 sample points were generated, from which horizons with 2048 sample points were obtained. In addition, two other kinds of traces were synthesized:

fBm with 131072 sample points and fBf with 2048 2048 sample points generated with the spectral synthesis method. fBm with 262144 sample points and fBf with 4096 4096 sample points generated again with the spectral synthesis method. This provided a test for possible discrepancies depending on the synthesis method and on the size of the samples. In all cases, to within statistical uctuations, the same results were obtained on the three kinds of traces. It is worth recalling that the maximum property and thus the fractal estimates are known to hold for the Brownian motion (with = 1 2 ), see Corollary 5.5, so index 1 2 fBf gives a way to assess the validity of our synthesis and estimation procedures. Figure 1 presents some realisations of fBf with = 0:1; 0:3; 0:5; 0:7; 0:9 and 2048 2048 sample points. (Note that, although X(0; 0) = 0, we have renormalised the colour scale in each case to make use of the whole colour range.) Superimposed is the function U t : the black dots are the points where X(t; u) attains a maximum w.r.t. to u. It is interesting to note how these points are distributed. Whilst Proposition 5.1 should be satis ed, U t has many of discontinuities, somewhat reminiscent of a L evy process. Figure 2 shows the corresponding horizon functions Z t .

Since we are always dealing with fBm (resp. fBf) de ned on 0; 1] (resp. 0; 1] 0; 1]), the number of samples must be interpreted as the inverse of the resolution. For instance, an fBm with 131072 sample points requires the following points to be generated: X(1=131072); X(2=131072); : : :; X(131071=131072); X(1). This is particularly important in testing the maximum property.

7 Numerical veri cation of the maximum conjecture

In this section we test the maximum Conjecture 3.2, that is whether for all xed 2 (0; 1) and all > 0, there exists a.s. = ( ; ) > 0 such that: X(t) X(t m ) jt m tj + for all t 2 0; 1] with jt m tj :

(7.1)
where X is index-fBm and t m is the point where X attains its maximum. As a qualitative test, we provide the following displays:

Figure 3 shows traces of fBm with = 0:1; 0:3; 0:5; 0:7; 0:9 and 131072 sample points, along with the corresponding curves Y (t) = X(t m ) jt m tj + , with = 0:05 (as the case of Brownian motion is of special interest, we show two realizations with = 0:5). Notice that, whilst in the majority of cases, the graph is only locally below the theoretical envelope (as expected), for several realizations, Y (t) is a tight upper bound for the whole path,

Figure 4 shows zooms of the traces of Figure 3 around the maxima in the cases of = 0:1; 0:5; 0:9, again with the Y (t) curves. In all cases, the path of X falls nicely below the envelope in the neighbourhood of t m , Figure 5 gives a zoom around the maximum of an fBm with = 0:5 along with both the upper envelope Y (t) and the lower theoretical envelope V (t) = X(t m ) jt m tj corresponding to the H older regularity of X(t). It may be seen that, for a given , the neighbourhood corresponding to the validity of the lower bound is smaller than the one for the upper bound. This property was observed on most traces.

In order to perform quantitative tests, we used the traces generated by the Choleski method described above. Thus, for each value of , we produced 100 sample paths with 131072 points. The di culty in obtaining a meaningful veri cation is that equation (7.1) does not tell us how and should be chosen. It could be that for a given , say 0.05, the corresponding is so small that traces much larger (i.e. with a much higher resolution) than the ones that were synthesised are needed. Conversely, if we choose large enough and small enough, even if X does not satisfy the maximum property, equation (7.1) is likely to hold for this particular choice of ( ; ). We adopted the following strategy. For each value of , and for each = 0:02; 0:04; : : :; 0:1, we computed ( ; ) as the largest number such that, in at least 95% of the cases (that is for at least 95 traces out of the 100 ones available for each ), equation (7.1) holds. In other words, the value of selected is such that, \almost surely", all points in a -neighbourhood of the maximum fall under the conjectured envelope, while this property fails for any larger neighbourhood.

Table 1 shows the results for the 500 traces corresponding to = 0:1; 0:3; 0:5; 0:7; 0:9. Note that is expressed in absolute units (not as a number of samples) obtained from the ratio (number of samples)=131072. The same tests were performed on the traces with the remaining values of and also on the sample paths obtained with the spectral synthesis method (both resolutions), yielding essentially the same behaviour.

In view of the results in Table 1, it is reasonable to conclude that the numerical simulations do con rm the maximum conjecture. Firstly, even for values of as small as 1=45 of (for = 0:9) and 1=5 of (for = 0:1), there exists a neighbourhood such that, for at least 95% of the realizations, the expected inequality holds. This neighbourhood has at least size 0.01 (corresponding roughly to 1300 samples points) or 1% of the size of the path, and in some case this size increases to 0:03. Secondly, Table 1 reveals a nice evolution of the size of the neighbourhood with respect to : for low values of , the \discontinuous" increase of matches the antipersistant property of the paths, whilst for large values of , the increase in is smooth paralleling the long range dependance in the realizations. Finally, recalling that the result holds true for = 0:5, we see that essentially the same things happen for all values of the exponent. Of particular interest is the fact that, for all values of except 0.08, the largest neighbourhood sizes are not attained in the case of usual (index-1

2 ) Brownian motion. Thus, according to the numerical experiments, there does not seem to be anything special about = 1 2 . Note that all the tests and comments above make a hidden assumption. Given a random element ! of the sample space, that is for a xed realization X t (!) of an fBm on 0; 1], we rst perform a sampling to get the points X t 1 (!); : : : ; X tn (!), and then nd the maximum (t m ; X tm ) over this nite set of samples. However, depending on n, it might well be that t m bears little relationship to the location of the \true" maximum t , that is the maximum of X t over 0; 1]. One way around this is as follows. Knowing the H older exponent of X t and its -H older norm, both of which are estimated by (2.3), it is possible to construct from X t 1 (!); : : : ; X tn (!) a continuous path X(t) that is an upper bound to X t for all t. If X tm is larger than X(t) for all t smaller than t m 1 and larger than t m+1 , then we are assured that the true maximum lies in t m 1 ; t m+1 ], with a value X(t ) \close" to X tm , and this is the best we can achieve. Otherwise, t might be located arbitrarily far from t m . This last situation may for instance occur in gure 13. In this case, the numerical veri cation of the maximum conjecture is meaningless, since the maximum is not well de ned on the sample path. In our experiments, we discarded such traces. For more on the problem of locating the maximum of a H older function from a sampling (see [START_REF] Lutton | H older functions and Deception of Genetic Algorithms[END_REF].

Estimation of the box dimension of the horizon

In this section we present numerical evidence to support Theorem 5.3, in particular that the box dimension of the horizon is equal to 2 for an index-fBf.

There are various methods for estimating the box dimension of a graph F. The crudest one, which simply consists in counting the number of boxes N( ) of a given size that intersect F and then making a linear regression of log(N( )) against log( ) su ers from many well-known defects, and does not in general give good results, so more re ned procedures are needed. The rst procedure we used is the variation method [START_REF] Tricot | Curves and Fractal Dimension[END_REF].

De ne the -variation of Z: The results obtained with the two methods on the 900 traces generated through Choleski decomposition are displayed in Table 2. For each value of and both estimators, the mean value and the standard deviation as measured on the 100 corresponding traces are shown. Once again, tests on horizon functions obtained through spectral synthesis yielded similar behaviour. Note that in all tests that consider the whole path, one hundred points at each end of the sample were discarded to avoid boundary problems.

Several conclusions may be drawn from Table 2. Firstly, both estimators agree, within statistical uctuations, for all values of . Since they are based on di erent procedures, this gives us con dence that the results are meaningful. Secondly, all the means are close to the predicted values (the largest error is 0.05 for the variation method and 0.04 for the regularization dimension), well within error bars given by the standard deviation. Thirdly, the values obtained for both the means and the standard deviations are of the order of what is usually reported in the litterature when testing the performances of new estimators of the box dimension on classical \fractal" curves such as fBm or Weierstrass functions. Recalling that such curves are self-similar while the horizon has no direct scale invariance properties gives even more con dence in the results. Finally, note that better results are obtained for medium range values of the box dimension, while low values are overestimated and large values are underestimated. Again, this is common in box dimension estimation, with similar behaviour observed in almost all situations. While methods exist to correct such bias, they usually make assumptions not applicable here. The estimate for = 0:5, where the theoretical value is known to be 1:5, is no better than for the neighbouring values = 0:3; 0:4; 0:6 and 0:7. This shows once again that there does not seem to be anything special about = 0:5. In summary, Table 2 clearly supports the parts of Theorem 5.3 pertaining to box dimension.

Estimation of H older exponents of the horizon

In this section we estimate the H older exponent H Z (t), given by (1.2), of the horizon function Z, and give numerical evidence to support the conclusion of Theorem 5.4. From the de nition, it is obvious that such a quantity is extremely di cult to estimate from numerical data, and we describe various steps that simplify the computation.

First, it is easy to see that (1.2) is equivalent to

H Z (t) = lim inf !0 + log osc Z (t; ) log : (9.1)
The form (9.1) leads to better estimates as it takes more data into account in the computation: instead of involving the values Z(u) Z(t) where only u varies in the neighbourhood of t, it depends on Z(u) Z(v) for all (u; v) in the neighbourhood of t. The second step is to eliminate the lower limit in the de nition of H Z (t). Indeed, estimating a lower limit from numerical data means that we are able to decide which neighbourhoods should be taken into account and which should be discarded. Without any further information on the behaviour of osc Z (t; ), this is a very hard task. Fortunately, we may invoke Theorem 5.4, which yields a ne measure theoretic property of H Z (t), much stronger that the basic de nition. Indeed, the proofs of Corollary 4.3 and Theorem 5.4 imply that, for almost all sequences of convergent to 0 at a geometric rate, we should obtain the right exponent. We may thus replace the lower limit by a plain limit, which we evaluate classically as a regression of log osc Z (t; ) vs. log . Thirdly, we shall use a quantity which is closely related to the pointwise H older exponent H Z (t), termed the local H older exponent, and de ned as H l Z (t) = lim inf !0 + log osc Z (t; ) log jt 1 t 2 j :

(9.2) where (t 1 ; t 2 ) is a pair of extreme points of the oscillation. Then H l Z is a more `stable' function of t that H Z , and is thus easier to estimate. This may be seen by noting that, while H l Z has to be a lower semi-continuous function, H Z need only be a lower limit of a sequence of continuous functions [START_REF] Guiheneuf | 2-Microlocal analysis and application in signal processing[END_REF]. We thus expect a less erratic behaviour for the local exponent. Further, though H Z and H l Z are not in general equal, they should coincide in our case: if the maximum property holds, H Z (t) must, a.s., equal for almost all t, implying that H l Z(t) will also a.s. equal for almost all t.

As before, we performed tests on the 900 traces generated by the Choleski method. The results are displayed in Table 2. For each value of , the mean value and the standard deviation of H Z as measured from 100 traces are given. Note that we are dealing here with two averaging procedures: one on all realizations and one along the path of each particular realization. Table 4 shows what happens if we consider only one particular sample path. For each , one horizon was chosen at random, and the mean value and standard deviation of the function H Z (t) was computed. Conversely, Table 3 displays the result obtained on averaging, for a xed t 0 , the values of H Z (t 0 ) on all 100 realizations. The conclusions of this numerical study parallel to a certain extent those of the preceeding section. However, due to its local nature, it is obvious that estimating a H older exponent is a much harder task than estimating a box dimension. In particular, the sampling has more in uence: in general, it entails two competing e ects that bias the estimation. Firstly, the oscillation as measured on the discrete data is always a lower bound for the \true" oscillation. This e ect is dominant when H Z (t) is small, because the discrepancy will be larger, and results in an overestimate of H Z (t). For large values of H Z (t), an \interference" phenomenon occurs, that is the oscillations around a \smooth" point t may be dominated by what happens in a smaller neighbourhood of an \irregular" point close to t. This translates into an underestimate of H Z (t). These heuristic expla- nations can be made precise under certain assumptions of scale invariance, these precise arguments are not available for the horizon, and we must be content with the expectation that H Z (t) will be overestimated for small values and underestimated for large values, a fact borne out in Tables 3, 4 and 5. The di culty in estimating H Z (t) is re ected in our tables by the relatively large values of standard deviation, especially when dealing with an average over one particular realization. This is not unexpected since this averaging procedure does not have a real statistical basis. Table 5 gives somewhat better resuls, although the biases are still large for extreme values of H Z (t), and the standard deviations are quite big. Note that in the three tables, the worst results are obtained for low values of H Z (t), possibly indicating that in our case, the discrepancy between the measured oscillation and the true one dominates the e ect of interference. Finally, the results for = 0:5 are again not very di erent from those of neighbouring values of the exponent, an argument in favour of the conjectured results. In conclusion, though the numerical estimates for the H older exponents are not as sharp as for the box dimensions, they do not discredit the conclusions of Theorem 5.4. Indeed, for = 0:4; 0:6 and 0:7, they match the theoretical predictions fairly closely.

Conclusion

The work in this paper gives strong support to the horizon conjecture, that the horizon of an index-fractional Brownian surface almost surely has the same H older exponents as the surface itself, with the box and Hausdor dimensions of the horison one less than those of the surface. We have proved this rigorously for the usual index-1 2 Brownian surface, and also for the index-Brownian surface for 0 < < 1 assuming the maximum property for index-Brownian functions. For all we have presented substantial computational evidence both for the horizon conjecture for index-Brownian surfaces and for the maximum conjecture for the index-Brownian functions.

0.02 0.04 0.06 0.08 0.1 = 0:1 0.03 0.05 0.05 0.05 0.05 = 0:3 0.03 0.05 0.05 0.07 0.08 = 0:5 0.02 0.03 0.05 0.08 0.08 = 0:7 0.01 0.02 0.03 0.04 0.05 = 0:9 0.02 0.025 0.03 0.033 0.04 0.15 0.18 0.25 0.17 0.33 0.18 0.39 0.12 0.49 0.1 0.57 0.15 0.72 0.13 0.83 0.15 0.87 0.17

Table 5: Estimated H older exponent at the horizon averaged at a given position over all realizations.

  Theorem 5.4 Assume that index-fBm has the maximum property. Writing Z(t) for the horizon function of the index-fBf X(t; u), we have that a.s.,

		H Z (t) = H ap Z (t) =	(5.11)
	for almost all t 2 0; 1].
	Proof. Lemma 1.3 and (2.3)imply that a.s.
		H Z (t) H ap Z (t) H ap Z (t)
	for all t. Proposition 5.2 along with (4.9) gives that for all t 2 0; 1], we have
		ap lim sup v!t	log jZ(v) Z(t)j log jv tj
	a.s., that is H ap Z (t)	a.s., and Fubini's theorem gives the conclusion.
	Corollary 5.5 Let Z be the horizon function of the usual (index-1 2 ) fBf. Then, with probability one, H Z (t) = H ap Z

  The main advantage of this procedure compared to the box method is that it is invariant under a change of scale. It also allows for a smooth evolution of the window size and leads to a fast implementation. Its well-known drawback is that oscillations are sensitive to noise, but this is not a problem in our setting.A second estimator was examined, based on the computation of the regularization dimension dim R , de ned as follows. Let (t) be a kernel function of Schwartz class S Thus, if the maximum property is true, we should have that dim R graph Z coincides with the common value of the box and Hausdor dimensions of the horizon Z. The advantage of using dim R is that it leads to a more robust estimator than other methods.

	where osc Z (t; ) is the oscillation of Z in the ball centered at t with radius ,
	osc Z (t; ) = sup t <u<t+	Z(u)	inf t <u<t+ Z(u)	(8.2)
	Then, assuming that the box dimension exists (that is dim B graphZ = dim B graphZ = dim B graphZ), we have
	dim B graphZ = max	( 1; lim !0 2 log(var Z ( )) log( )	!)	;	(8.3)
	see 5, Section 11.1], so dim B graphZ may be estimated as a regression of 1 2 log(var Z ( )) vs. j log j. such that Z = 1: (8.4)
	Let a (t) = 1 a ( t a ) be the dilation of at scale a > 0 and Z a be the convolution of Z with a . Since Z a 2 S, the length a of its graph is nite. The regularization dimension of graphZ is de ned as dim R graphZ = 1 + lim a!0 log( a ) log a (8.5)
	Thus dim R measures the speed at which the length of decreasingly regularised versions of a fractal graph tends to in nity. It is easy to check that, for a continuous curve C, one always has dim H C dim R C dim B C (8.6)
	var Z ( ) =	Z 1 0	osc Z (t; )dt	(8.1)

Table 1 :

 1 Evolution of the size of the neighbourhood where the maximum property is valid. .07 1.76 .09 1.71 .07 1.62 .05 1.47 .06 1.38 .07 1.29 .06 1.17 .07 1.15 .08 Regularization Dimension 1.88 .09 1.81 .06 1.69 .08 1.62 .06 1.49 .05 1.4 .06 1.31 .07 1.23 .07 1.14 .09

	: Variation Method	0.1 m std 1.87	0.2 m std	0.3 m std	0.4 m std	0.5 m std	0.6 m std	0.7 m std	0.8 m std	0.9 m std

Table 2 :

 2 Mean (m) and standard deviation (std) of the estimators for the box dimension at the horizon.

	: Estimated H older exponent	0.1 m std 0.16 0.12 0.26 0.15 0.34 0.13 0.42 0.11 0.52 0.1 0.59 0.13 0.72 0.12 0.84 0.15 0.86 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 m std m std m std m std m std m std m std m std

Table 3 :

 3 Estimated H older exponent at the horizon averaged over all positions and realizations.

	: Estimated H older exponent	0.1 m std 0.21 0.2 0.3 0.19 0.36 0.19 0.47 0.15 0.54 0.13 0.63 0.17 0.67 0.19 0.81 0.2 0.83 0.18 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 m std m std m std m std m std m std m std m std

Table 4 :

 4 Estimated H older exponent at the horizon averaged over one realization.

	: Estimated H older exponent	0.1 m std	0.2 m std	0.3 m std	0.4 m std	0.5 m std	0.6 m std	0.7 m std	0.8 m std	0.9 m std

fBm, has the maximum property.

fBf, so, using the maximum property (3.1), a.s. X(t; u) attains a unique maximum at u = U t 2 (0; 1), such that for some 2 > 0 X(t; u) X(t; U t ) ju U t j + if ju U t j 2 . Using the uniqueness of the maximum and a continuity argument, this implies that a.s. there is some (random) > 0 such that X(t; u) X(t; U t ) minfju U t j + ; g for all u 2 0; 1]:

(5.5) Set = minf 1 ; 1=( ) g. If jv tj < 1=( )

(5.6) then by (5.4) X(v; U t ) X(t; U t ) 1 3 jv tj :

(5.7) On the other hand, suppose (v; u) is such that jv tj < and ju U t j jv tj 1 :

(5.8) Then by (5.4),(5.5),(5.8),(5.6) and (5.3), X(v; u) X(t; u) + 1 3 jv tj X(t; U t ) minfju U t j + ; g + 1 3 jv tj X(t; U t ) minfjv tj ( + )(1 ) ; jv tj g + 1 3 jv tj X(t; U t ) 2 3 jv tj : Comparing with (5.7) it follows that if jv tj < the supremum sup 0 u 1 X(v; u) cannot be attained by u such that ju U t j jv tj 1 , and the result follows.

We can now conclude that Z(v) = Z t (v) if v is close to t. Proposition 5.2 Assume that index-fBm has the maximum property. Then for all t 2 0; 1], there exists a.s. > 0 such that, for all v 2 0; 1] with jv tj < , Z(v) = Z t (v):

(5.9) Proof. This is immediate from Proposition 5.1 and the de nitions (4.3) and (4.7) of U t and Z t (v).

We now combine our results to obtain the fractal properties of the horizon.

Theorem 5.3 Assume that index-fBm has the maximum property. Writing Z(t) for the horizon function of the index-fBf X(t; u), we have that dim H graphZ = dim B graphZ = dim B graphZ = 2 a.s.. Proof. By Corollory 4.3 and Proposition 5.2 we have for all > 0 and t 2 0; 1], that a.s. L(v : jv tj h and jZ(v) Z(t)j h) h 2

(5.10) for h su ciently small. Hence with probability one, for L-almost all t 2 0; 1], (5.10) holds for all 0 < h < h(t) for some h(t) > 0. Since > 0 is arbitrary, Lemma 1.2 gives dim H graphZ 2 a.s.. Since X satis es an almost sure uniform H older condition

(2.3), the opposite inequality follows from Lemma 1.3. We can obtain H older exponents in a similar way.


Jacques L Evy V Ehel
email: jacques.levyvehel@inria.fr

Evelyne Lutton
email: evelyne.lutton@inria.fr

Evolutionary signal enhancement based on H older regularity analysis

We present an approach for signal enhancement based on the analysis of the local H older regularity. The method does not make explicit assumptions on the type of noise or on the global smoothness of the original data, but rather supposes that signal enhancement is equivalent to increasing the H older regularity at each point. The problem of nding a signal with prescribed regularity that is as near as possible to the original signal does not admit a closed form solution in general. Attempts have been done previously on an analytical basis for simpli ed cases 11]. We address here the general problem with the help of an evolutionary algorithm. Our method is well adapted to the case where the signal to be recovered is itself very irregular, e.g. nowhere di erentiable with rapidly varying local regularity. In particular, we show an application to SAR image denoising where this technique yields good results compared to other algorithms. The implementation of the evolutionary algorithm has been done using the EASEA (EAsy speci cation of Evolutionary Algorithms) language.

Introduction

A large number of techniques have been proposed for signal enhancement. The basic frame is as follows. One observes a signal Y which is some combination F(X; B) of the signal of interest X and a \noise" B. Making various assumptions on the noise, the structure of X and the function F, one then tries to derive a method to obtain an estimate X of the original signal which is optimal in some sense. Most commonly, B is assumed to be independent of X, and, in the simplest case, is taken to be white, Gaussian and centered. F usually amounts to convoluting X with a low pass lter and adding the noise. Assumptions on X are almost always related to its regularity, e.g. X is supposed to be piecewise C n for some n 1. Techniques proposed in this setting resort to two domains: functional analysis and statistical theory. In particular, wavelet based approaches, developed in the last ten years, may be considered from both points of view [START_REF] Devore | Lucier Fast wavelet techniques for near optimal image processing[END_REF][START_REF] Donoho | De-noising by soft-thresholding[END_REF].

Our approach in this work is di erent from previous ones in several respects. First, we do not make explicit assumptions on the type of noise and the coupling between X and B through F. However, if some information of this type is available, it can readily be used in our method. Second, we do not require that X belong to a given global smoothness class but rather concentrate on its local regularity. More precisely, we view enhancement as equivalent to increasing the H older function Y (see next section for de nitions) of the observations. Indeed, it is generally true that the local regularity of the noisy observations is smaller than the one of the original signal, so that in any case, X should be greater than Y . We thus de ne our estimate X to be the signal \closest" to the observations which has the desired H older function. Note that since the H older exponent is a local notion, this procedure is naturally adapted for signals which have sudden changes in regularity, like discontinuities. From a broader perspective, such a scheme is appropriate when one tries to recover signals which are highly irregular and for which it is important that the restauration procedure yields the right regularity structure (i.e. preserves the evolution of X along the path). An example of this situation is when denoising is to be followed by image segmentation based on textural information: Suppose we wish to di erentiate highly textured zones (appearing for instance in MR or radar imaging) in a noisy image. Applying an enhancement technique which assumes that the original signal is, say, piecewise C 2 , will induce a loss of the information which is precisely the one needed for segmentation, since the denoised image will not contain much texture. The same di culty occurs in other situations such as change detection from noisy sequences of aerial images or automatic monitoring of the evolution of lung diseases from scintigraphic images: in such cases, the criterion for a change is often based on a variation of the irregularity in certain regions, and one needs to preserve this information.

In addition to the examples given above, other situations where such conditions occur are turbulence data analysis and characterization of non-voiced parts of voice signals. Since the method is highly non linear and quite complex, an analytic solution is not possible. We thus resort to a stochastic optimisation method based on evolutionary altgorithms.

The remaining of this paper is organized as follows. Section 2 recalls some basic facts about H older regularity analysis, which is the basis of our approach. The denoising method is explained in section 3. The evolutionary implementation based on EASEA is then detailed in section 4, nally numerical results on both 1D and 2D signals are displayed in section 5.

H older Regularity Analysis

A popular way to measure the local irregularity of signals is to consider H older spaces. We will focus in this paper on the H older characterizations of regularity. To simplify notations, we assume that our signals are nowhere di erentiable. Generalisation to other signals simply requires to introduce polynomials in the de nitions 7]. Let 2 (0; 1), R. One says that f 2 C l ( ) if: 9 C : 8x; y 2 : jf(x) f(y)j jx yj C

Let: l (f; x 0 ; ) = sup f : f 2 C l (B (x 0 ; ))g l (f; x 0 ; ) is non increasing as a function of .

We are now in position to give the de nition of the local H older exponent :

De nition 1 Let f be a continuous function. The local H older exponent of f at x 0 is the real number: l (f; x 0 ) = lim !0 l (f; x 0 ; ) Since l is de ned at each point, we may associate to f the function x ! l (x) which measures the evolution of its regularity. This regularity characterization is widely used in fractal analysis because it has direct interpretations both mathematically and in applications. It has been shown for instance that l indeed corresponds to the auditive perception of smoothness for voice signals. Similarly, simply computing the H older exponent at each point of an image already gives a good idea of its structure, as for instance its edges 10]. More generally, in many applications, it is desirable to model, synthesize or process signals which are highly irregular, and for which the relevant information lies in the singularities more than in the amplitude. In such cases, the study of the H older functions is of obvious interest.

In 12], a theoretical approach for signal denoising based on the use of a H older exponent and the associated multifractal spectrum was investigated. We develop here another enhancement technique that uses the information brought by the local H older function, which is simple from an algorithmic point of view, and yields good results on several kind of data.

Signal Enhancement

We adopt in this paper a functional analysis point of view. This means that we do not make any assumption about the noise structure, nor the way it interacts with the data. Rather, we seek a regularized version of the observed data that ful lls some constraints. A statistical approach, classically based on risk minimization, will be presented elsewhere 11].

Let X denote the original signal and Y the degraded observations. We seek a regularized version X of Y that meets the following constraints: a) X is close to Y in the L 2 sense, b) the (local) H older function of X is prescribed.

If X is known, we choose X = X . In some situations, X is not known but can be estimated from Y . Otherwise, we just set X = Y + , where is a userde ned positive function, so that the regularity of X will be everywhere larger than the one of the observations. We must solve two problems in order to obtain X. First, we need a procedure that estimates the local H older function of a signal from discrete observations. Second, we need to be able to manipulate the data so as to impose a speci c regularity. A third di culty arises from the following analysis: Assume the simplest case of an L 2 signal corrupted by independent white Gaussian noise. It is easy to check that almost surely Y = 1 2 everywhere, because a) 1 2 is the regularity of the noise, b) X 0 since X 2 L 2 , c) the regularity of the sum of two signals which have everywhere di erent H older exponents is the min of the two regularities. Thus Y does not depend on X, and one cannot go back from Y to X . This fact casts doubts on the e ciency on the whole approach, since the information it is based on is degenerate in this case.

All these problems are solved once one realizes that the mathematical notion of H older regularity is an abstraction that makes sense only asymptotically. One needs to analyze carefully how it should be adapted to a nite setting, much in the same way as what is done for abstract white noise. In particular, we are interested in a perceptual notion of regularity: If two 2D functions A and B are such that A < B , but an imaging of A and B at a given resolution yields the contrary visual impression that A is smoother than B, then of course our algorithm should go with the perceptual information. In other words, in practical applications we are not interested in the asymptotic behavior, but in the scales which are really in the image, and our estimate of should re ect this fact. This means precisely one thing: the estimation procedure should yield results in agreement with what is perceived, and not care for the \true" , which may or may not be accessible from the nite data. To go back to the example above, while it is true that at in nite resolution the sum \signal + white noise" would look much the same as white noise as far as regularity is concerned, this is not the case at nite resolution, where the in uence of the signal is still perceptible. In addition, since our procedure is di erential, i.e. we wish to impose X = Y + or X X = 0, for estimated X and Y , we do not care about constant bias.

We will use a wavelet based procedure for estimating and controlling the H older function. This is made possible by results in 7] and 6] which imply that: Proposition 1 Let f j;k g j;k be an orthonormal wavelet basis, where as usual j denotes scale and k position, and assumes that is regular enough and has su ciently many vanishing moments. Then, X has local H older exponent at t if and only if for all (j; k) such that t belongs to the support of j;k , jc j;k j C2 j( + 1 2 ) (1

)
where C is a constant and c j;k is the wavelet coe cient of X.

Although [START_REF] Devore | Lucier Fast wavelet techniques for near optimal image processing[END_REF] is only an inequality, it suggests that one may estimate X (t) by linear regression of log(jc j;k j) w.r.t. to the scale j (log denotes base 2 logarithm) considering those indices (j; k) such that the support of j;k is centered above t.

Of course this will be only approximate, but since (1) is a necessary and su cient condition, and if enough wavelet coe cients are \large", we may hope to obtain results su cient for our purpose. Two points are essential in this estimation procedure: { The estimation is obtained through a regression on a nite number of scales, de ned as a subset of the scales available on the discrete data. This avoids the pathologies described above concerning the regularity of the sum of two signals. In particular, it is possible to express the H older function of the noisy signal Y = X + Gaussian white noise as a function of X , and thus to estimate conversely X from Y 11].

{ The use of (orthonormal) wavelets allows to perform the reconstruction in a simple way: Starting from the coe cient (d j;k ) of the observations, we shall de ne a procedure that modi es them to obtain coe cients (c j;k ) that verify (1) with the desired , and then reconstruct X form the (c j;k ). We may now reformulate our program as follows: For a given set of observations Y = (Y 1 ; : : : ; Y 2 n ) and a target H older function , nd X such that jj X Y jj L 2 is minimum and the regression of the logarithm of the wavelet coe cients of X above any point i w.r.t. scale is ( (i) + 1 2 ). Note that we must adjust the wavelet coe cients in a global way. Indeed, each coe cient at scale j subsumes information about roughly 2 n j points. Thus we cannot consider each point i sequentially and modify the wavelet coe cients above it to obtain the right regularity, because point i + 1, which shares many coe cients with i, requires di erent modi cations. The right way to control the regularity is to write the regression contraints simultaneously for all points. This yields a system which is linear in the logarithm of the coe cients: L = A where is a (2 n ; 2 n+1 1) matrix of rank 2 n , and L = (log jc 1;1 j; log jc 2;1 j; log jc 2;2 j; : : : log jc n;2 n j); A = n(n 1)(n + 1)

12

(1) + 1 2 ; : : : ; (2 n ) + 1 2 Since we use an orthonormal wavelet basis, the requirements on the (c j;k ) may nally be written as: minimize X j;k (d j;k c j;k ) 2 subject to: 8 i = 1; : : : ; 2 n , n X j=1 s j log(jc j;E((i 1)2 j+1 n ) j) = M n ( (i) + 1 2 )

(2)

where E(x) denotes the integer part of x and the coe cients s j = j n+1 2 ,

M n = n(n 1)(n+1) 12 and equation ( 2) are deduced from the requirement that the linear regression of the wavelet coe cients of X above position i should equal ( (i) + 1 2 ). Finding the global solution to the above program is a di cult task ; in particular, it is not possible to nd a closed form formula for the c j;k . In the following, we show how this problem can be addressed with an evolutionary algorithm. [START_REF] Wall | GAlib home[END_REF] Evolutionary signal enhancement with EASEA An evolutionary technique seems to be appropriate for the optimisation problem described in equation ( 2): a large number of variables are involved, and the function to be optimised as well as the contraint are non linear. We describe in this section an implementation based on the EASEA 13] language and compiler.

EASEA (EAsy Speci cation of Evolutionary Algorithms) is a language dedicated to evolutionary algorithms. Its aim is to relieve the programmer of the task of learning how to use evolutionary libraries and object-oriented programming by using the contents of a .ez source le written by the user.

EASEA source les only need to contain the "interesting" parts of an evolutionary language, namely the tness function, a speci cation of the crossover, the mutation, the initialisation of a genome plus a set of parameters describing the run. With this information, the EASEA compiler creates a complete C++ source le containing function calls to an evolutionary algorithms library (either the GALIB or EO for EASEA v0.6). Therefore, the minimum requirement necessary to write evolutionary algorithms is the capability of creating non-object-oriented functions, speci c to the problem which needs to be solved.

In our case, the evolutionary optimisation involved to enhance a signal (1D or 2D) was implemented using a simple structure on which genetic operators were de ned. We used GALib 4] as the underlying evolutionary library.

We describe below the implementation for 1D signals, an implementation for image denoising was also produced based on the same principle, and results for 1D and 2D data are presented in the next section.

The Haar wavelet transform has been used to produce the d j;k associated to the observed signal Y . We also suppose that we know the desired H older exponents (i) (either (i) = Y (i)+ where the Y (i) are the H older exponents of Y and is a user de ned regularisation factor, or (i) is set a priori).

The problem is to nd some multiplication factor u j;k such that c j;k = u j;k d j;k ; j 2 0::n 1]; k 2 0::2 j 1]. As is usual in wavelet denoising, we let unchanged the rst l levels and seek for the remaining u j;k in 0; 1]. The genome is made of the u j;k coe cients, for j 2 l::n 1] and k 2 0::

2 j 1].
These coe cients are encoded as a real numbers vector of size SIZE MAX = 2 n 2 l , which can be written using EASEA syntax as :

GenomeClass { double U SIZE_MAX]; }
The EASEA Standard functions sections contain the speci c genetic operators, namely:

1. The initialisation function: Each u j;k coe cient is randomly set to a value in 0; 1]. Two initial solutions are also put in the initial population : u j;k = 1: and u j;k = 2 k .

2. The crossover function: a barycentric crossover has been easily de ned as follows : Let parent1 and parent2 be the two genomes out of which child1 and child2 must be generated, and let alpha be a random factor:

Fitness = X j;k ((1 u j;k ) c j;k ) 2 + W X i j u (i) (i)j
We use the GALib steady state genetic engine with replacement percentage of 60% and a selection by ranking. Crossover and mutation probabilities are xed respectively to 0.9 and 0.001. Genome size, Population size, and number of generations are xed for each experiment, see section 5.

Numerical Experiments

We rst show results of enhancement on synthetic 1D data. The original signal is a generalised Weierstrass function 2] with X (t) = t for all t, (i.e. X(t) = P 1 n=0 2 nt sin(2 nt ); t 2 0:1]) which has been corrupted by white Gaussian noise. Figure 1 shows the original signal, the noisy one, and the result of the enhancement procedure. For comparison, a denoising using a classical wavelet shrinkage is also displayed. For both procedures, the parameters were set so as to obtain the best t to the known original signal. It is seen that, for such irregular signals, the H older regularity based enhancement yields more satisfactory results. The constraint was to nd H older exponents that verify (t) = t for all t. Parameters of the evolutionary algorithm were as follows : Denoising using wavelet shrinkage method (left), and using the evolutionary scheme with prescribed (t) = t (right). Our second example deals with a synthetic aperture radar (SAR) image. A huge literature has been devoted to the di cult problem of enhancing these images, where the noise is non Gaussian, correlated and multiplicative. A ne analysis of the physics of the speckle suggests that it follows a K distribution 9]. Classical techniques speci cally designed for SAR image denoising include geometric ltering and Kuan ltering. Wavelet shrinkage methods have also been adapted to this case 3]. Figure 2 show an original SAR image, its denoising with the H older method and with soft thresholding. Notice how the river, with appears with a \ " shape in the middle of the image is nicely uncovered by the regularity based enhancement. As no a priori knowledge about H older exponents of the signal was available, the constraint was to nd H older exponents that verify denoised (t) = original (t) + for all t. Parameters of the evolutionary algorithm were as follows : We have shown in this paper how an evolutionary algorithm can be applied to a signal or image enhancement technique based on a \fractal" analysis. Good results have been obtained as well as on 1D or 2D signal in comparison to other techniques in an a ordable computation time.

Quality of the results is also very dependent from the quality of H older exponents estimation. Much better estimate of X (t) can be obtained by measuring the oscillations of X in balls centered at t and of radii k , and then regressing the logarithm of these oscillations w.r.t. the logarithm of the k . However, this procedure does lead to a more complex inverse problem, i.e. obtaining a signal with prescribed regularity. Future work will be devoted to an evolutionary formulation of this problem.

  \GenomeClass::crossover: double alpha = (double)random(0.,1.); if (&child1) for (int i=0; i<SIZE_MAX; i++) child1.U i] = alpha*parent1.U i] + (1.-alpha)*parent2.U i];if (&child2) for (int i=0; i<SIZE_MAX; i++) child2.U i] = alpha*parent2.U i] + (1.-alpha)*parent1.U i]; \end 3. The mutation function: Mutation is a random perturbation of radius SIGMA = 0:01, applied with probability PMut on each gene. \GenomeClass::mutator: // Must return the number of mutations as an int int NbMut=0; for (int i=0; i<SIZE_MAX; i++) if (tossCoin(PMut)){ NbMut++; Genome.U i]+=SIGMA*(double)random(-1.,1.); Genome.U i] = MIN(1.,Genome.U i]); Genome.U i] = MAX(0.,Genome.U i]);} if (NbMut==0) identicalGenome=true; // saves evaluation time return NbMut; \end 4. The evaluation function: The tness function has two aims: minimise P ((1 u j;k ) c j;k ) 2 , making sure constraint (2) is satis ed, i.e. the H older exponents are the ones we want. The constraint is integrated to the tness function using a high penalisation factor W :

Generalised

  Weierstrass funtion (left) + noise (right).

1

  Left: Estimated H older exponents of the original function (left), of the function + noise (middle) and of the reconstructed function (right).

Fig. 1 .

 1 Fig. 1. Results on a generalised Weierstrass function (t) = t

7 Fig. 2 .

 72 Fig. 2. Experiments on a 256x256 area of a SAR image (courtesy IRD)