Wavelet Transform Modulus Maxima
The following work has now been accepted by publication by Astrophysical Journal. It should be in the July 2007 edition. The paper is available on the ApJ website , the Max Millennnium eprint archive , as part of the package with the software, or by email request. Please note the following text was scripted at the preliminary stage, so the results in the paper may differ from those below. Please email me if you have any problems with implementing the code or if you have any questions or comments. I would also like to start a list of projects from people in different fields of study who are using the wtmm method. This could be a good way of disseminating knowledge and experience, and starting some discussion over problems. If you are interested in participating at any level, please let me know. If you can send a brief paragraph, weblink, or journal article I will be happy to put everything online.
Introduction
I found the following links and papers very useful in this work. Rather than me rewriting all the basics glance over these if you need a good introduction.
A Brief Overview of Multifractal Time Series is an excellent basic description of the reasoning behind the modulus maxima, and the results obtained from it. Also has C software which works well.
Stanley & Meakin, 1988 Nature review article
Chhabra et al., 1989 Physical Review A article, containing first description of Canonical methods, and why the D(q) approach and 'histogram' approach (which Lawrence et al use) fails. Also applies it to 1d turbulence and makes the link with thermodynamics.
Muzy et al. , 1991 Physical Review Letters article also uses the canonical method, and links into turbulence.
Muzy et al., 1993 Physical Review E article compares the wavelet transform modulus maxima to the structure function approach.
The Arneodo, et al., 1995, Physica A, 213, 232-275 article also belongs to this set of Muzy articles.
Degaudenzi., 1998 article describing the results of the WTMM as applied to pollen counts.
There are nice extensions in these two papers, which I have not coded yet.
These two Struzik papers ( struzik_a and struzik_b ) explain a method to extract the local effective Holder exponent (i.e., retain the local temporal information on the multifractality)
Osahi et al. , 1993 Physical Review E article describes the wtmm applied to positive and negative changes separately
The Reasoning of Modulus Maxima
The basic idea is to describe the partition function over only the modulus maxima of the wavelet transform of the signal. This allows for accurate calculation at negative q. Then by using the canonical method there is no need to Legendre transform the tau(q) to get the D(h), hence removing the errors associated with this. This also removes the problems associated with the histogram approach (no knowledge of prefactor) and the D(q) approach (poor log-log fit due to lacunarity, and scatter due to finite data). The WTMM also allows for the choice of mother wavelet which will determines the number of vanishing moments, important for removing any unwanted trends in the data. NB Despite this, the tau(q) curve is still best calculated using the standard log-log plots of scale against partition function.
The exact details of the algorithm are in the papers above, but in essence the idea is to calculate a normalised measure at each moment q, and scale, and then calculate the Hausdorff dimension of this support. This produces a plot of the D(h) spectrum, where D(h) is the Hausdorff dimension of the set of singularities of Holder exponent, h. For comparison, the F(alpha) spectrum is the Hausdorff dimension of singularities of strength, alpha. it can be shown these are the same. The D(h) works well because although the log-log plot still shows oscillations, these are periodic. Hence a linear fit is still accurate.
Interpreting Results
Here is the wavelet amplitude of a lightcurve (RHESSI 6-12KeV, below), with the Modulus maxima lines over plotted in black. The x axis is time (seconds from start of lightcurve) and the y axis is scale (both axis actually should be multiplied by the cadence of the data, 4 sec). A window of 4/5 of the entire lightcurve has been chosen, so as neglect any edge effects. The analysing wavelet is the third derivative of a Gaussian and hence polynomials of order two and below are removed. This figure alone exhibits a multitude of information. Its shows the branches of the detected singularities. At each scale, only the MM which can be traced from previous smaller scales are included, hence the (white) maxima at the right is omitted. The 'branch-like' nature of the modulus maxima is typical of self-organised systems in Nature. The largest power is in the decay phase of the flare, and extends to the largest scales (naturally, the flare takes longer to decay than to rise). Interestingly, the modulus maxima extend to shorter and shorter scales before the flare onset (i.e., tracking lines up from the bottom, they go up to smaller scales just before flare onset.
The D(h) spectrum is shown below. Loads of important parameters can be read directly from this plot. The D(h) spectrum is the Haussdorff dimension of set of singularities of Holder exponent, h. Hence the D_max(h) is representative of most frequent Holder exponent, and can be thought of as the Hurst exponent which best describes the entire system. A homogeneous system will have a delta function D(h), a multifractal system will show a range of D(h). This D_max(h) corresponds to the q=0, i.e., D(h(q=0)), and the D value will be approximately 1 for a 1-D system, i.e., the fractal dimension of the support is 1. Any overshoot or undershoot is a problem caused by the algorithm. Another test for algorithm problems is a discontinuity at D(h(q=0)) which is caused by a collection of insufficient data, poor mother wavelet choice (need to use a higher order mother wavelet) and the old negative q problem. This also shows as a strange blip in the tau(q) curves.
Various parameters can be used to describe this curve.
- The h value of D_max(h) (=1.40) is the most frequent Holder exponent and describes the memory of the system.
0 < h < 0.5 is an anti-persistent exponent; the system reverses itself more frequently and covers less distance than a random walk
0.5 < h < 1.0 is persistent exponent; the system exhibits long term memory effects and hence dependence on current conditions
h > 1.0 is a regular process, with long term memory
- The minimum value of h (=0.82) corresponds to the strongest singularity characterizing more the rarefied zone of most pronounced modulus maxima (i.e., enhances at q>0). Hence this is the Holder exponent of the reddest regions above, at the onset of the flare.
- The maximum value of h (=2.15) corresponds to the weakest singularity characterizing most dense zone of lower modulus maxima (i.e., survived by q<0). This is the Holder exponent of the blips in the decay of the flare.
- The symmetry information compares the inhomogeneity in q>0 and q<0
(h(max)-h(min))/2. = 1.4920934
h(q=0) = 1.4049253
So this data is more weighted toward the larger h (weaker singularities)
- The value of the generalized dimensions can also be read off this curve. As above,the Capacity Dimension D(h(q=0)) = 0.99. Furthermore the Information Dimension (the triangle) D(h(q=1)) = 0.85 (describing the scaling behaviour of information,; how much information is lost as the system evolved in time). The Correlation Dimension D(h(q=2)) = 0.64 (characterizing the chaotic attractor).
- The curve is multifractal - this has never been shown before. It exhibits the characteristic single hump in D(h) and single break in tau(q).
- The algorithm removes the low-order polynomial trend, the modulation of the grids (order of 4 seconds, again polynomial), and the errors due to the attenuator state changes (which I'm informed are small and will only show up at small scales anyway). Essentially the removal of instrumental effects is intrinsic to the algorithm
- The curves all undershoot slightly (I need to put on error bars) meaning D(h(q=0)) equals 1, or tau(q=0)=-1,
Future
First of all this is a pretty result, well worth putting out. Other groups are trying to use this method but no-one has published, yet. It uses very few parameters of an arbitrary nature - I've explored the result over a range of q, delta_q, order of derivative, removing edge effects by mirroring, removing edge effects by padding, taking different range of scales.
It would be nice to see the local Holder exponent (i.e. the struzik work above), to see if the onset of the flare (an extreme event) has a precursor. It would be good to see this at a number of flares. As Alex suggested, this could be a good way of characterizing different types of flares / other astronomical objects, e.g., GOES, BATSE data. Of particular interest could be a nice test of the Neupert effect - the DOG wavelet is a multiscale derivative.
This has to be applied to 2D data and compared to the improved box counting techniques Paul is working on. Is the MF spectrum of a magnetogram a precursor of a flare event (an extreme event). Interestingly, people are using a version of this technique for object detection .
The software and test data is available for download . Please email me if you have any problems with implementing the code or if you have any questions or comments.
