Next Previous Contents

2. Pointwise Hölder Exponent

Many techniques have been developed to estimate pointwise (and local) Hölder exponents, none of which give satisfactory results in all cases. Estimating a local irregularity index on discrete data without any a priori assumption is indeed a difficult task. One way to obtain robust results is to use a parametric approach. Such estimators have been developed mainly in the case of fractional Brownian motion and its extensions, such as multifractional Brownian motion. One such method is implemented in Fraclab. Non parametric methods, which are more numerous, generally give the correct estimation only when some technical conditions are met. Fraclab currently includes three such estimators, which are, in order of increasing quality : a standard method based on the continuous wavelet transform, with various improvements ; a second one relying on the use of generalized iterated functions systems. And a third one based on the evaluation of the oscillations. Note finally that more robust estimates of both the pointwise and local exponents, which work under broader assumptions, are available when one computes the 2-microlocal exponents. See the help corresponding to this topic below.

2.1 Parametric estimation of the exponent of an fBm or an mBm.

For an fBm or an mBm with smooth H(t), the local and pointwise exponent are the same. The parametric estimator implemented here is based on (4), and is proved to converge almost surely as soon as one is dealing with a discrete trace of an fBm or mBm.

As in most processing tools of Fraclab, the first line of the sub-menu window, entitled Input Signal, should display the signal currently analyzed. A known bug here is that the current signal of the Variables list sometimes does not appear by default. When this occurs, or if you wish to switch to another signal by highlighting another selection in the Variables list, press Refresh. There is only parameter that needs to be chosen before you launch the estimation, which is the window size. The optimal choice, as the name indicate, is the theoretical best size, i.e. the size which minimizes the risk. As is often the case in such situations, visually more pleasing results are obtained by using a slightly larger size, typically by multiplying the optimal size by a logarithmic factor: choosing enlarged just does this. In many situations, you'll find that it gives a clearer picture of the meaningful variations of H, because the variance has been decreased (at the expense of course of an increase in bias). When you suspect that your analyzed signal is an fBm, i.e. the exponent is the same all along the path, a more judicious choice is to use the whole signal as a window. This is the third choice, constant H. Once you have decided a window size, press Compute as usual. The output signal, denoted estim_Hölder_func_sig# (where sig is your input signal), is a 1D vector which gives the exponent at each point. The two other buttons are the customary Help and Close buttons. Note finally that computing times can get very large if your input signal has more than 1024 points.

A known bug in this sub-menu is that border effects have not been satisfactorily taken into account. As a result, the estimates on the first points are sometimes weird (e.g. you get a plateau which has no real meaning).

2.2 CWT-based estimation

The Continuous-Wavelet-Transform-based non parametric estimator of the Hölder function is probably the most well known method for estimating exponents, although it is not the most precise one. It is based on a remarkable property of wavelet coefficients: They are bounded from above by a quantity that depends on scale and on the Hölder exponent. Since this bound is optimal, it can serve as a basis for an estimation method. However, the bound is practical only of one assumes that the relevant coefficients for estimating the exponent at point t are located "nearly above" t. This is equivalent to assuming that the pointwise and local exponents coincide. This condition is hard to verify in practice, and somewhat restrictive.

When you launch this method, you get a window entitled "Local Hölder Exponent Estimation" (recall that this will also be the pointwise one if the procedure is applicable). In this simple form, you need just decide at which point you want to compute the exponent, by entering a value in the Estimation time box (the default is 1, which is a bit unfortunate, since the first point will always suffer from border effects). On hitting Compute, the value of the estimated exponent will appear in the Local Scaling Exponent box. If you need more control on the computations, hit Advanced compute. A new window appears, entitled Hölder Function or Local Exponent. Check first in the Input Signal box that you are dealing with the desired signal, otherwise Refresh it in the usual way. You will then specify various parameters related to the wavelet transform: fmin and fmax let you choose the minimum and maximum frequencies of analysis. The default values are the ones yielding maximal span compatible with the size of the signal. You may change the extreme frequencies either by typing values under fmin and fmax, or by using the predefined values on the menus to the right. The Voices parameters governs the number of intermediate frequencies between fmin and fmax at which the continuous wavelet coefficients will be computed. Be warned that giving an excessive number of voices may result in large computing times for long signals. Checking the Mirror item will deal the border effects by mirroring the signal at its extremities. Otherwise, zero-padding is used. Finally, you may choose the Size and Type of your analyzing wavelet: available wavelets are the Mexican Hat, and the real and analytic Morlet wavelet. The size may be any positive number (this parameter is not available for the Mexican hat). Once all the parameters that define the wavelet transform are chosen, hit Compute WT. The output signal is a matrix of size "number of voices" x "size of the original signal". It is called cwt_signal#, if "signal" is the name of your data, and where # is as usual an incremental parameter. It should appear in the Input CWT box just below the Compute WT button. You may want to view the continuous wavelet transform using the View menu. Note that Fraclab recognizes wavelet transforms, and display them differently form regular images. In particular, it uses a fixed aspect ratio (this is useful for instance if the number of voices is much smaller than the size of the signal), and the "jet" color-map, which often allows to highlight the important structures. If you want to view the transform as a normal image, or make other changes in the appearance, use the functionalities of the View menu described in the Overview help file.

The Refresh button to the left of the Input CWT box lets you load a wavelet transform which would already be present in the Variables list of the main window. This avoids computing several times the same transform. Once you are happy with your transform, move to the lower part of the window, which performs the actual computation of the exponents. Decide first if you wish to compute a Single Time Exponent (default), otherwise uncheck this box, and it will become Hölder Function. In the first case, you need to give the Time Instant at which the exponent will be estimated (with again the unfortunate default "1"). In the second case, this choice becomes unavailable, because you will compute the exponent at each point. In fact, this last statement is not quite true: Because the CWT-based estimator is a bit slow, it has been arbitrarily decided that only one point every four points would actually be processed. This is why the output signal is four times smaller than the input one. This fact may seem somewhat strange, because it means that the results of this menu cannot really be used in actual applications you may have. However, our justification is that the CWT-based estimation does not usually give good results. This window is thus mainly included because wavelet-based estimators are widely used, and for comparison purposes. Some of the reasons why the method often fails are explained at the end of this paragraph.

Recall that the exponent is obtained thanks to the fact that the coefficients are bounded from above by a certain quantity. As a consequence, in many cases, more relevant estimates are obtained if one chooses, at each scale, the largest coefficient among the ones which are located "nearly above" the current point t. This is the default in Fraclab, as is indicated by the fact that the Local Maxima box is checked. If you unmark this box by pressing the button to the left of Yes (which subsequently becomes No), then the program will use the coefficients which are in the column strictly above t. In case you want to use the local maxima, you need to tell Fraclab where to look for the largest coefficient. In other words, you must define precisely what is meant by "the coefficients located nearly above the point t". This is the purpose of the Radius and Scale Depth parameters. The default values of respectively 8 and 1 mean that, starting from the point t, the program will look for the largest coefficients at the scale immediately above (this is the "1"), and in a spatial neighbourhood of size 8. These two parameters are available for both the estimation at a single point and at all points. To the contrary, the Specify Regression Range choice is only meaningful for a single time exponent. When Specify Regression Range is chosen, you will be able to decide which range of frequencies is to be used for the estimation of the exponent (we explain how at the end of this paragraph). Otherwise, when Full Range Regression is selected, the estimation is performed using all coefficients between fmin and fmax. The exponents will be obtained through a regression of the logarithm of the wavelet coefficients modulus with respect to scale, and you may choose the type of regression: The default is Least Square Regression. Other choices are Weighted Least Square, Penalized Least Square, Maximum Likelihood and Lepskii Adaptive Procedure. Only the last one is not well known. See (6) for details.

When you hit Compute, three things may happen: if Hölder Function was checked, the program will output a signal called em/signal_Ht#/, with the usual naming conventions. If Single Time Exponent was checked instead, with the option Full Range Regression, the program will display a graphic window showing in abscissa the log-scale, and in ordinate the log of the modulus of the coefficients. This display helps to check if the regression is meaningful, i.e. if the points on this graph are well aligned. In that view, the regression line is displayed in red. In parallel, you get the value of the estimated exponent in the box Local Scaling Exponent. If Specify Regression Range was checked instead of Full Range Regression, you also get a graphic window, but this time you have additional information and control: to the left of this window, you see the part of the wavelet transform which lies in the neighbourhood of the selected time instant. Note that a grey levels color-map is used. You will also see two red curves, which show the region in which local maxima are searched for when this option is activated. Finally, you should see, at each scale, i.e. on each line of the wavelet transform, a circle which indicates which coefficient was chosen at this particular scale. If no local maxima were used, this circles will lie exactly above t, while otherwise they might be distributed anywhere between the two red curves. The graph to the right of the window is the same as in the previous case, i.e. it shows in abscissa the log-scale, and in ordinate the log of the modulus of the coefficients. Now Specify Regression Range means that, using the large cross that appears when you point to this window, you will be able to select a range of frequencies between which the regression will be performed. Of course, you want to choose a region where the points are reasonably well aligned. Sometimes there is no such region. In other cases, there might be two or more sub-parts in the graph where linear behaviours are observed. Since we are interested here in local exponents, you should in general prefer to choose the region containing the highest frequencies. To actually select the region, click first on the left end of the chosen interval, then click again on its right end, or vice-versa. Each time you have selected a frequency range, Fraclab will compute an estimated exponent and show its value in the caption above the graph to the right of the window. In addition, the regression line will be displayed in red. If you want to test another region, just click again its endpoints. Once you are satisfied with a result, hit Enter on your keyboard, and the cross will disappear. The exponent will then appear in the Hölder Function or Local Exponent window in front of Local Scaling Exponent. You may now close the graphic window showing the evolution of the wavelet coefficients.

To understand some of the reasons why the method is not so good, try it on a deterministic Weierstrass function synthesized with the defaults options. Launch the Advance compute procedure, and use again the defaults, except that you estimate at the middle point of your signal, i.e. at abscissa 128. Observe how the log-log plot display oscillations: By choosing a regression range on various parts of the frequency interval, you will be able to get essentially any exponent you want, including negative ones.

2.3 GIFS-based estimation

This method is applicable only when the pointwise and local exponents coincide. In contrast with the CWT-based one, it is very fast. The principle is the following: The extension of IFS called GIFS allows to construct signals which are able to approximate (in L2 and L-infinity norms) any L2 signal with arbitrary precision. GIFS are just IFS where the number of maps and the various parameters are allowed to change at each scale (see (2) for details). The first step in this estimation method is to compute a GIFS which approximates at best the original signal. In that view, one starts by computing the discrete wavelet transform of the signal. The parameters of the GIFS are then simply obtained as ratios of the wavelet coefficients. Once the GIFS is known, the estimation of the exponent is easy, because there is an analytical formula which gives the exponent at each point of a GIFS as a function of its parameters. Because the formula in only valid in the limit of infinite resolution, the obtained result is of course a finite size approximation.

From a practical point of view, you start as usual by checking and maybe updating the Input signal. You then choose which approximation procedure you want to use, i.e. the Limit type in a choice of regression and Cesaro (choosing one or the other option will in general hardly affect the result). Finally, you may choose which Wavelet to use, although in the current implementation of Fraclab only Daubechies 4 is available. On pressing Compute, you get an output signal, named alphagifs_signal#, which contains the estimated Hölder function. Note that this program assumes that the input data contains 2^N points. If this is not the case, only the first 2^N points will processed, where N is the largest integer such that the 2^N does not exceed the actual number of points in the signal. The remaining points will simply be ignored.

Test this method on a deterministic Weierstrass function synthesized with the defaults options.

2.4 Oscillation-based estimation

In contrast with previous ones, this method does not assumes that the local and pointwise exponents coincide. It truly tries to estimate the pointwise one. The principle is very simple: it just uses the definition of the exponent as a measure of how the oscillation of the signal in a neighbourhood of a given point t behaves with respect to the size of the neighbourhood. As usual, you first check the Input data name box and Refresh it if needed. You then need to choose what are the minimal and maximal sizes of the neighbourhood that will be used to investigate the behaviour of the oscillations. Enter the appropriate values in the Nmin and Nmax boxes, either directly or using the menus. Any integer will do, as long as Nmin is smaller than Nmax and Nmax is compatible with the size of the signal. These values should be understood as follows: Nmax = 8, for instance, means that the largest neighbourhood will be composed of 8 points to each side of the point where one wishes to perform the estimation. In other words, the neighbourhood will be a window of size 17 sample points centered at the point of interest. Increasing the values of both Nmin and Nmax yields smoother estimates, at the expense of precision. A value of Nmin larger than 1 roughly means that a high frequency cut-off is in effect. Larger values of Nmax let you include more low frequency information in the estimation.

Finally, you may choose a Regression Type from the usual choice of Least Square Regression, Weighted Least Square, Penalized Least Square, Maximum Likelihood and Lepskii Adaptive Procedure (see reference (6)). On hitting Compute, you get the output data pht_signal#, which contains the estimated exponent at each point. Note that refinements of this oscillation based method exist. They use a Bayesian framework that allows to minimize the effect of discretization. We hope to include these extensions in a future release of Fraclab.


Next Previous Contents