Characterizing the structure of nonlinear systems using gradual wavelet reconstruction

In this paper, classical surrogate data methods for testing hypotheses concerning nonlinearity in time-series data are extended using a wavelet-based scheme. This gives a method for systematically exploring the properties of a signal relative to some metric or set of metrics. A signal continuum is defined from a linear variant of the original signal (same histogram and approximately the same Fourier spectrum) to the exact replication of the original signal. Surrogate data are generated along this continuum with the wavelet transform fixing in place an increasing proportion of the properties of the original signal. Eventually, chaotic or nonlinear behaviour will be preserved in the surrogates. The technique permits various research questions to be answered and examples covered in the paper include identifying a threshold level at which signals or models for those signals may be considered similar on some metric, analysing the complexity of the Lorenz attractor, characterising the differential sensitivity of metrics to the presence of multifractality for a turbulence time-series, and determining the amplitude of variability of the Hölder exponents in a multifractional Brownian motion that is detectable by a calculation method. Thus, a wide class of analyses of relevance to geophysics can be undertaken within this framework.


Introduction
Because of the wide ranging occurrence and varied nature of nonlinearity in geophysical time series (Johnson et al., 2005;Khan et al., 2005;Roux et al., 2009), gaining an understanding of the sources of any nonlinearity is an important topic.The presence of nonlinearity can be tested by ap-Correspondence to: C. J. Keylock (c.keylock@sheffield.ac.uk) plying some metric, such as time asymmetry (Schreiber and Schmitz, 1997) or maximal Lyapunov exponent to a time series of outputs from a system and comparing the response to surrogate data that are linear variants of the original signal.A test for any significant difference can be developed within this framework (Theiler et al., 1992;Schreiber and Schmitz, 1996).The intention of this paper is to pursue this matter in a new direction.The approach developed here moves away from the acceptance/rejection framework of the hypothesis test for nonlinearity to ask: How similar to the original data do the surrogates need to be to avoid rejection of the null hypothesis?From this, it is possible to develop new research questions within a surrogate data framework, such as: -Which parts of the time series need to be identical between the data and the surrogate in order to prevent the rejection of the null hypothesis (i.e.what are the most complicated parts of the original time series)?
-Does the range of values for a metric calculated for a set of surrogates that are not statistically different to the original data include the value of this metric for a model of that system (i.e. is the model validated)?
-Do different, but related nonlinear or chaotic time series exhibit differences in how similar their surrogates need to be to the original data to avoid rejecting the null hypothesis (i.e. are there differences in complexity of the series)?
-Do different measures applied to a nonlinear or chaotic time series exhibit differences in how similar the surrogates need to be to the original data to avoid rejecting the null hypothesis (i.e. are there differences in sensitivity of the measures used to characterise chaotic or nonlinear systems)?
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
Our approach, which we term Gradual Wavelet Reconstruction, permits these questions to be answered.This is illustrated by the examples presented in the final parts of the paper.Before this, we explain the gradual wavelet reconstruction approach, which requires us first to review briefly relevant literature on hypothesis testing for surrogate data in nonlinear science.

Hypothesis testing for nonlinear time series using surrogate data
The surrogate data methodology as proposed by Theiler et al. (1992) and enhanced by Schreiber and Schmitz (1996) is a common technique with applications in studies of river meandering (Frascati and Lanzoni, 2010), ice core data (Kwasniok and Lohmann, 2009), environmental turbulence (Basu et al., 2007;Keylock, 2009) and the magnetosphere (Pavlos et al., 1999), as well as other disciplines beyond geoscience, such as medicine (Mormann et al., 2005).Typically, one generates surrogate data that are stochastic realisations from a Gaussian linear system with the same values and (to some error tolerance) Fourier spectrum as the original data and employs a metric to see if the observed time series is significantly different to the surrogates.For a two-tailed hypothesis test at a significance level, a, if the value of the metric for the original data is less than or greater than that for all of the (2/a) − 1 surrogate datasets then the null hypothesis that the original data is a realisation of a Gaussian linear process will be rejected.
An effective method for producing surrogate data that preserve the values and, to some error level, the Fourier spectrum of the original data is the Iterated Amplitude Adjusted Fourier Transform (IAAFT) algorithm (Schreiber and Schmitz, 1996).Given a discrete time series g n , n = 1,...,N this algorithm proceeds as follows: 1. Store the squared amplitudes of the discrete Fourier transform of g n (i.e.G 2 f = | N 1 g n e i2π f n/N | 2 ); 2. perform a random shuffle of g n to give g (0) n ; 3. Subsequently, iterate a power spectrum step and a rankorder matching step on g (j ) n as follows: (a) Take the Fourier transform of g (j ) n and replace the squared amplitudes with G 2 f , while retaining the phases.Given the initial random sort, this means that the spectrum should be preserved but with random phases.Invert the transformation with the amplitudes replaced; (b) Replace the values in the new series g (j ) n by those in g n using a rank-order matching process.This preserves the set of original values in the dataset but deteriorates the quality of spectral matching, which explains why the Fourier amplitudes are only replicated approximately; 4. Repeat until a convergence criterion is fulfilled or any changes are too small to result in any re-ordering of the values.
The phase randomisation part of the algorithm will destroy temporal organisation in the original series that contributes to any nonlinearity, while the fact that the amplitudes of the spectrum are approximately preserved and the values of the original dataset are completely preserved mean that differences on some metric between the data and the surrogates cannot be attributed to these sources, which could be sources of difference between two linear time-series.Hence, a significant difference implies, either the presence of some form of nonlinearity in the original data or that these data are sampled from a non-Gaussian, linear process.Subsequently, this algorithm has been refined by groups such as Venema et al. (2006) who relaxed step 3(b), by imposing the values of g n more gradually in order to improve convergence.
The IAAFT algorithm was first implemented in the wavelet domain by Keylock (2006) using a Maximal Overlap Discrete Wavelet Transform (MODWT), which is described in the Appendix to this paper.Because a wavelet transform is a time-frequency decomposition (see A1), the use of a single IAAFT results in the constrained randomisation of a time-series of wavelet coefficients representing one particular frequency band (or scale).Thus, with J different scales, performing an IAAFT at each scale, results in a full-randomisation of the wavelet coefficients.Because the IAAFT algorithm does not alter the values for these coefficients, the wavelet power spectrum obtained from the MODWT (which is proportional to the variance of the wavelet coefficients) is unaffected by this transformation.The convergence of this method compared to the standard IAAFT and the enhanced method of Venema et al. (2006) was tested by Keylock (2008a), while the approach developed in Keylock (2006) has subsequently sparked interest in other new ways for describing stationarity of time series (Borgnat andFlandrin, 2009, 2010).Keylock (2007) presented a refinement to the earlier method, which still used the MODWT and the IAAFT, but fixed in place particular wavelet coefficients to provide a flexible means for designing surrogates.It is this algorithm that underpins gradual wavelet reconstruction, as is explained in the next section.

The algorithm
As established in the Appendix for the continuous wavelet transform and as stated for the MODWT, the square of the wavelet coefficients w 2 (j,k)/j 2 is the energy function of a time-series signal decomposed over different scales/frequencies, j , and positions along the time series, k (Vela- Arevalo and Marsden, 2004).For a signal of length N = 2 J there will be a total of k = 1,...,N wavelet coefficients at each scale, j produced by the MODWT.The total energy content of the real-valued wavelet transformed signal is proportional to and we define ρ as some chosen fraction of E. If the w 2 j,k are placed in a J × N length vector in descending rank order, the smallest number of w 2 j,k required to attain ρ can be determined by cumulating the squared wavelet coefficients until their sum as a fraction of E attains ρ.We term these the fixed wavelet coefficients.The other coefficients are randomised as explained below.Hence, ρ provides a measure that can be used to vary the degree of similarity between the surrogates and original data (Keylock, 2007).For ρ = 0 there are no fixed coefficients and the resulting surrogate will be similar to that obtained using the IAAFT.Trivially, for ρ = 1 all coefficients are fixed, no randomisation occurs and the surrogates and data are identical.
More formally, if W ∈ R + is the set of J × N coefficients, w 2 j,k , placed in descending rank order then, with 1 ≤ n ≤ (J ×N ) acting as an index for W , the set of fixed coefficients, F ⊆ W , is given by the first n elements of W that fulfils the condition n i=1 w 2 i E ≥ ρ.Hence, the {w 2 1 ,...,w 2 n } ∈ F are the smallest number of coefficients that fulfils the energy proportion, ρE.
The algorithm for generating a surrogate data series using this approach may now be stated.This is an improved version of the algorithm given by Keylock (2007).The wavelet used in this paper is a Daubechies (1993) least-asymmetric wavelet with 16 vanishing moments for effective frequency localisation.The centre frequency (i.e. the frequency that maximises the Fourier transform of the modulus of the wavelet) is 0.6774 and the relation between scale, j , and the negative logarithm of the pseudo frequencies has a proportionality constant of 0.693.We made use of MATLAB and software accompanying Percival and Walden (2000), written by Charlie Cornish and available from WMTSA (2006) to implement the MODWT.

Choose a value for ρ;
2. Perform a wavelet decomposition of the time series into a J × N array using the MODWT and determine the fixed coefficients for this ρ as explained, above; 3. For each wavelet scale, determine if any of the N wavelet coefficients are to be fixed.If they are not, apply the IAAFT algorithm to give a randomised realisation of the coefficients at this scale.If they are: 4. Invert the wavelet transform to produce a new time series of length N using the original approximation coefficients (which will be a constant for a stationary series if a full wavelet decomposition is undertaken, as is done throughout this paper) and the randomised detail coefficients (see Eq. A9-A12 for an explanation of MODWT approximation and detail coefficients).
5. Because of a loss of the original values in the dataset from this operation (and a subsequent loss of matching of the power spectrum when they are re-imposed), repeat stages (3) and (4) of the IAAFT algorithm to ensure convergence for the dataset as a whole.
Fixing in place more coefficients as ρ increases means that the surrogates become progressively more similar to the data.Applying the IAAFT algorithm to each scale of the wavelet transform ensures that the coefficients have the appropriate autocorrelation function and can be reconstructed appropriately because they are a feasible realisation of a MODWT.Note that by terminating the hierarchical MODWT algorithm at any stage (thereby increasing the number of frequencies contained within the approximation coefficients relative  to a full decomposition), this algorithm can be tailored to only operate at selected scales.However, in this paper, ρ pertains to the fraction of energy in a full decomposition of the time series.
The re-introduction of unfixed coefficients in our method provides a contrast with those techniques where some subset of the initial wavelet coefficients are used to reconstruct a process, such as the study by Venugopal and Foufoula-Georgiou (1996).The primary advantages of our approach are the preservation of the original values, the improved preservation of the Fourier spectrum, and the ability to determine the effect of the unfixed coefficients upon some metric applied to the data by consideration of the variability of the surrogates.However, in the examples considered in Sects.4 and 5, signal reconstructions based simply on the fixed coefficients in the style of Venugopal and Foufoula-Georgiou's work are also used.This helps illustrate the role played by the unfixed coefficients for the chosen signal metric.The black lines in Fig. 6b-d and the dotted lines in Fig. 8, are examples of this approach.
Given a set of stochastic surrogates at various choices for ρ, the research question stated in the introduction can be reexpressed as: At what choice of ρ is there no longer any difference between the value of our metric for the surrogates and for the original data?

Illustration and explanation of the method
Consider the 1.64 s (2 13 values) of the longitudinal, u 1 , and vertical, u 2 , components of turbulent velocity time series obtained at 5000 Hz in a 1 m cross-section wind tunnel in the wake of a 100 mm high fence, which are shown in Fig. 1.Data were obtained by the author at a Taylor Reynolds number for the far field of 150 and recorded 0.5 m downwind of the fence at a height of 55 mm above the base of the tunnel.Figure 2 shows the process of developing a surrogate for u 1 at ρ = 0.5, with just the operation of the algorithm at wavelet scale j = 8 shown, which was a local maximum for the wavelet spectrum and 23% of the coefficients w j =8,k were fixed for this ρ.As the IAAFT algorithm converges, the initially randomly located unfixed coefficients are adjusted to respect the stored Fourier amplitudes of the original set of coefficients.The fixed coefficients can be seen clearly as smooth regions in Fig. 2b, while Fig. 2e and f show that the parts of the surrogate time-series that differ most from the original data are where the coefficients are not fixed, as expected.

Surrogate representations of multifractal signals
Various types of geophysical data have been analysed in terms of their multifractal characteristics, including atmospheric processes (Tessier et al., 1993;Venugopal et al., 2006), topography (Gagnon et al., 2003), and seismicity (Nakaya and Hashimoto, 2002).The aim of this paper is not to replicate such characteristics explicitly in the surrogates, but to provide a means of generating surrogates that vary in their nature as a function of ρ.As ρ increases, any multifractality in the underlying dataset will be increasingly preserved.To see this, note that while it is possible to analyse the multifractal characteristics of a signal using windowed spectra (Pikovsky et al., 1995), it is more common to adopt a wavelet perspective (Muzy et al., 1991).It is well known (e.g., Mallat, 1999) that the multifractal characteristics of a signal can be approximated by calculating the wavelet transform modulus maxima, chaining together maxima across scales and then forming the partition function where q ∈ is a selected power that measures the scaling behaviour of Z(q,j ), ξ is a maximum of the wavelet transform modulus maxima, and Li indexes each of these maxima.Scaling exponents are calculated by and it has been shown by Bacry et al. (1993) and Jaffard (1997) that these scaling exponents can be related to the support of the multifractal distribution via a Legendre transform: where the Hölder/lipshitz exponents, α fall within the support of the multifractal spectrum D(α). Figure 3 shows that as ρ increases, the time series converges upon that in Fig. 1a.In addition, it illustrates how the resulting wavelet modulus maxima are preserved.For example, note that at scale ≈ 5 and t ≈ 10.9 s an energetic feature is fixed for ρ ≥ 0.3, but the feature at scale ≈ 5 and t ≈ 11.25 s is only fixed for ρ ≥ 0.5.
For the replication of the wavelet transform modulus maxima it is necessary for the surrogates to preserve the multifractal spectrum of the original data, and that this is only accomplished over all scales at high values for ρ.However, it is also the case that there is both an imprecision in the calculation of the α exponents due to the limitations of the resolution and length of the datasets, as well as imprecision in the signal itself owing to instrument noise etc.Hence, at a somewhat lower value for ρ there will be no significant difference between surrogates and a multifractal dataset, depending on the width of the support of the multifractal spectrum.This issue is examined in Sects.6.3 and 7.

Evaluating finite size effects on randomisation
Clearly, in the limit of ρ = 1 the surrogates and dataset are identical and no randomisation occurs.Hence, the issue of finite size effects is complex as it will be a function of the length and nature of the time series, and the chosen value for ρ.For example, the highest value (ρ = 0.999) used in this paper leaves, on average (calculated over 200 surrogates) 447 values in the same positions in Fig. 1 as they are in the surrogate (≈ 5% of N = 8192).For a very short segment of this time series of 256 velocity values and averaged over 500 surrogates, for ρ ∈ {0.5,0.9,0.99},2.7 (1%), 15.5 (6%) and 49.4 (19%) of values were fixed in place.These values may be compared to similar results for a time series also of 256 points, but of a very different structure -the sunspot data analysed in Keylock (2007).In that case, for ρ ∈ {0.5,0.9,0.99}, on average 2.8 (1%), 13.4 (5%) and 37.5 (15%) of values were fixed in place.Hence, finite size effects need to be considered for very high values for ρ when datasets are short because a perceived increase in ρ will not have altered the nature of the time series.However, the example applications of the technique in this paper retain sufficient degrees of freedom for sufficient randomisation to occur.
Please note that in all the examples presented in the rest of the paper, unless otherwise stated, 19 surrogates are generated.As explained in the introduction, this is sufficient for a one-tailed test at the 5% significance level (Sects.4, 6 and 7) or a two-tailed test at the 10% significance level (Sect.5).An increase in the number of surrogates can be used to either reduce the significance level or enhance the statistical power of the test but becomes computationally demanding when a number of surrogates must be generated for various choices of ρ.The value used here is not dissimilar to that used elsewhere in the literature (e.g., Pavlos et al., 1999) and can be increased for more accurate discrimination between choices of ρ.

An application to the time asymmetry behaviour of a laser intensity time series
Figure 4a shows 1024 values from the Santa Fe laser time series, x L , (Huebner et al., 1989).This is a well-known test data series in nonlinear science.Example surrogate series for different choices of ρ are given in Fig. 4b-f.The surrogates shown in Fig. 4 are those corresponding to the median value for the surrogate asymmetry in Fig. 5.While visually, a choice of ρ ≈ 0.5 qualitatively begins to resemble the original data, gradual wavelet reconstruction is used to study the behaviour of these data using the temporal asymmetry (or skewness) measure A (Schreiber and Schmitz, 1997): where the standard choice in testing for nonlinearity is to choose λ = 1 (Schreiber and Schmitz, 1997).There is a logic for choosing λ = 1 for this dataset because the autocorrelation function has crossed zero by λ = 2, going from 0.53 at λ = 1 to −0.19 at λ = 2. Adopting an ensemble of different choices for λ gives additional information on the nature of the nonlinearity within the dataset and provides further criteria that one could aim to replicate when attempting to model a time series.
Figure 5 gives values for A λ=1 for 19 surrogates, at 14 choices for ρ, as well as the value for the laser data, A λ=1 L , which is indicated by a dotted line.From Fig. 5, the null hypothesis is rejected until ρ = 0.97.Hence, although there are visual similarities between Fig. 4d and the original intensity data, a much higher choice for ρ is required to preserve the key elements of the signal with respect to asymmetry.Figure 6 shows the original Santa Fe laser series, x L , for reference (Fig. 6a) together with the difference, x / L , between x L and surrogate data series (in grey and displaced vertically by ±40 (Fig. 6b and c) or ±140 (Fig. 6d), as well as a series defined by the fixed wavelet coefficients (i.e. with no stochastic, unfixed coefficients added), which is shown in black, for the three choices of ρ.In every case, the surrogate series for x / L that is displaced downwards is that with the median value for A λ=1 in Fig. 5, and that displaced upwards is the series with the maximum.The key difference between the series whose value for A λ=1 exceeds that for A λ=1 L (the upper grey line in Fig. 6b) and all the other data shown (including the fixed part of the data series for ρ = 0.97) is that the unfixed coefficients have acted to remove the discontinuity that occurs after 70 samples, where there is a sudden transition in the behaviour of x L .Hence, it would appear that, in addition to the general saw-tooth nature of the laser intensity, representing this type of discontinuity correctly is essential if a model for this system is to replicate the asymmetry characteristics of the original data.
Studying A λ=4 and A λ=6 , which are the lags greater than zero with the minimum and maximum autocorrelations (R = −0.62 and R = 0.75, respectively), one finds that the null hypothesis is rejected until ρ = 0.97 and until ρ = 0.999, respectively.Thus, more rigorous model validation can be accomplished by deploying additional choices for λ.Going further, an ensemble of different metrics could also be employed, a topic that is considered in Sect.6. Huebner et al. (1989) propose two models for their laser data, one based on the Lorenz equations (Lorenz, 1963):  Huebner et al. (1989) found that these two models gave a reasonable fit to the original data in terms of their value for  6) and ( 7), which may be compared to the series, x L , in Fig. 4a and a surrogate with ρ = 0.97 in Fig. 4f.
the correlation dimension measure adopted in Sect.5, implying that Lorenz-type models are suitable for modelling such time-series.To test this hypothesis with respect to our skewness/asymmetry measure, we integrated both sets of equations using a time step of 0.01 and choosing the values: b = 0.25, R = 15, σ = 2, and δ = 0.05, as per Huebner et al. (1989).In both cases, the time series for y 2 1 was downsampled such that the time to the first zero crossing of the autocorrelation function matched that in the original dataset (2 samples) and A λ=1 was calculated for series of 1024 values.We obtained A λ=1 = 2.207 for the Lorenz model and A λ=1 = 2.624 for the complex Lorenz model.These results are higher than the value for the laser data of A λ=1 = 2.008, but it is not immediately clear if this difference is significant.
Using the results in Fig. 5 we see that at ρ = 0.97 there is no significant difference between the original data and the surrogates for A λ=1 .Both asymmetry values for the models are much greater than the largest value of A λ=1 = 2.017 found for the 19 surrogates at this choice of ρ.Going further, 200 surrogates were generated for ρ = 0.97 and a Ryan-Joiner test for normality showed the A λ=1 values to be normal at the 10% significance level.Based on the standard deviation of 0.0078, the asymmetry values for the models are 79 and 25 standard deviations from the value for the data.Hence, the probability of obtaining the models' asymmetry values based on surrogates at a value for ρ with the greatest intrinsic variability that preserve the asymmetry properties, is vanishingly small.The difference in the nature of the model signals is illustrated in Fig. 7.For the additional choices of A λ=4 and A λ=6 , marked differences are also evident, with only the results using Eq. ( 6) for λ = 4 anywhere close to those for the data (3.5 standard deviations too high).Hence, the gradual wavelet reconstruction approach to model validation would suggest that Lorenz-type models are not validated with respect to the asymmetry of the original data.

An application to the Lorenz system
The Lorenz equations are the paradigmatic chaotic system, although it is only relatively recently that a comprehensive study of all three parameters of this model has been undertaken (Barrio andSerrano, 2007, 2009).In this section of the paper we employ classical choices for b = 8/3 and σ = 10, but consider values for the Rayleigh number, R, that include the classical, globally attracting chaotic attractor (R = 28), a value (R = 24.29)that gives a chaotic attractor with a pair of stable attracting rest points, a value of R = 24.75 by which the stable points have been eliminated (see Kaplan and Yorke, 1979), and a value within the regime of an intermittent transition to chaos identified by Manneville and Pomeau (1979) and Pomeau and Manneville (1980 The Lorenz equations were solved with a time step of t = 0.01, with results recorded every tenth time step from t = 1000.The accuracy of our numerical method was checked by using the Gottwald and Melbourne (2005) test for chaos applied to y 1 .The transition to chaos at R = 166.0616found using this method was in agreement with the value found using the methods of Barrio and Serrano (2007) (personal communication from Roberto Barrio).
In this study we employed the correlation dimension, D c , as a means of characterising the attractor (Grassberger and Procaccia, 1983) based on a Gaussian kernel method and using the output for y 1 , with lags and Theiler windows defined based on the decay of the autocorrelation function and on false nearest neighbours, respectively.We tested our method for long (40 960 points) and short (4096 points) datasets.Grassberger and Procaccia (1983) quote a value of D c = 2.05 ± 0.01 for the correlation dimension at R = 28, which was matched successfully by both of our datasets (D c = 2.052 and D c = 2.055, for long and short datasets, respectively).Hence, we employed the shorter series in analysis.It is possible to obtain erroneous, finite correlation dimensions for stochastic processes (Schertzer et al., 2002).However, by working with a system that is known to exhibit chaos and by deploying high values for ρ that ensure the basic structure of the Lorenz attractor is fixed in place (much as it would be for data from a Lorenz attractor with noise),means that we have generated correlation dimensions for data that are approximating the original, chaotic attractors.
Figure 8 illustrates short time series for y 1 for the four choices of R examined here.Figure 9 gives the correlation dimension as a function of embedding dimension, D e , for our four choices of R. The embedding dimension is the dimension of the phase space into which the time series is embedded based on delayed versions of the original series (Takens, 1981).An accurate estimate for D c requires D e to be sufficienly large to capture the dimension of the attractor (i.e. at least 3 for the Lorenz system) but not so great as to introduce errors from finite size effects.For a time series, y 1 it is possible to form a set of N − (D e − 1)L vectors: each of which defines a point in this embedding space, where L is a lag and N is the number of values in the time series.Gradual wavelet reconstructions based on 19 surrogates are shown for two choices of ρ, 99% (blue error bars) and 99.9% (red error bars).In addition, reconstructions at these values for ρ are shown based purely on the fixed wavelet coefficients without using the IAAFT algorithm to re-introduce the unfixed coefficients (as described at the end of Sect.3.1).These are indicated by the dotted lines, with the circles showing ρ = 0.99 and the squares ρ = 0.999.The error bars are displaced a small horizontal distance from the integer value for D e and indicate the mean and ±2 standard deviations by horizontal lines.
The most similar plots are Fig.9b and c, both of which are within the same regime of behaviour for R according to Kaplan and Yorke (1979).In these cases, at ρ = 0.99 the dataset built from just the fixed coefficients clearly differs from the original data.The surrogate data (blue error bars) are generally even further from the original data on average, but within the ±2 standard deviation tolerance of both the original data and the fixed coefficient surrogate (particularly when R = 28.00).The higher choice for ρ results in a convergence of both types of surrogate to the original data.Thus, at ρ = 0.99, on average for these two cases, the addition of the unfixed w j,k to the surrogates results in greater error than a lack of precise preservation of the data histogram or wavelet spectrum from just using the fixed coefficients.However, by ρ = 0.999 this other error source is dominant and realisations built from just the fixed coefficients contain greater error.A more in-depth analysis could examine the precise values for ρ at which this transition in the dominance of different error sources occurred.The results for the Lorenz system are given by a black line, with those for series reconstructed from the fixed wavelet coefficients for ρ = 0.99 (dotted line with circles) and ρ = 0.999 (dotted line with squares) and surrogate series for ρ = 0.99 (blue) and ρ = 0.999 (red).The latter summarise results for 19 surrogates according to mean (central horizontal line) and ±2 standard deviations (vertical lines).They are translated slightly from the integer value for D e for clarity.The situation differs in Fig. 9a, where the error bars for the gradual wavelet reconstruction are all small and lie close to the original data.However, for ρ = 0.99 there is a clear difference for the surrogate built purely from the fixed coefficients, which sits outside the error bars for the surrogates.In this case, failure to preserve the wavelet spectrum and histogram accurately has had a significant effect on D c , while randomisation by the unfixed coefficients has a minimal effect.
The situation differs again in Fig. 9d, where this time at ρ = 0.99 it is the realisation from the fixed coefficients that lies significantly closer to the original data than the surrogates.Hence, the randomised coefficients generate greater error than the failure to preserve the histogram or spectrum.For this case, it is also notable that by ρ = 0.999, while both types of surrogates have converged upon one another, none have converged on the original data.Table 1 lists the proportion of coefficients fixed at the two chosen thresholds.For ρ = 0.99 and R = 167.0there is a small proportion of fixed coefficients (i.e.there is high energy in relatively few values).This means that the randomisation within the surrogates is causing a relatively weak convergence on the scaling behaviour of D c .However, by ρ = 0.999 more coefficients are preserved for this dataset than the others yet the surrogates still differ significantly (at the 10% level) from the original data.This shows that the attractor for the Pomeau and Manneville intermittency regime is more complex than for the other values for R used here in the sense that, the value for D c in the data can only be replicated by fixing in place a higher proportion of the wavelet energy and a greater proportion of the wavelet coefficients.In contrast, while the dimension of the attractor is higher when R = 24.29,surrogate series at ρ = 0.99 can adequately capture its behaviour.Thus, gradual wavelet reconstruction can be used to elucidate additional information on the nature of an attractor, here providing a classification of complexity ranging from the simpler case of Fig. 9a (lower ρ is sufficient to capture the behaviour) through the intermediate cases shown in Fig. 9bc, to the more complex case in Fig. 9d.That R = 24.75 and R = 28.00exhibit similar behaviour mimics the similar structure of their attractors.A wavelet method using a Daubechies wavelet with 8 vanishing moments is indicated by Dau8, while "O" indicates the oscillation method.The number following the "O" is the largest exponent used in the calculation for δ.Hence, "O8" means that the bins for δ ranged from 2 1 to 2 8 .Least-squares regression was used unless a suffix "L" (Lim Inf regression) or "P" (penalised least-squares regression) is included.See (FRACLAB, 2006) for more details.(c) Shows the sinusoidal function for α(t) used to generate a multifractional Brownian motion (blue line) and then estimated values for α(t) from that multifractional Brownian motion time series.The black line is the Dau8 algorithm, the green line is O5 and the red line is O10.
6 The Hölder characteristics of a turbulence time series have been a number of studies that have tried to characterise the multifractal characteristics of turbulence (e.g., Meneveau and Sreenivasan, 1987;She and Leveque, 1994) owing to the well-known intermittency characteristics (e.g., Frisch et al., 1978) that lead to a departure from Kolmogorov's 1 3 scaling as discussed by Kolmogorov (1962).However, the predictions of the latter's log-normal model differ from the log-Poisson model of She and Leveque (1994) and analyses based on universal multifractal scaling (Schertzer and Lovejoy, 1992;Schmitt et al., 1992) provide an alternative framework for classifiying these processes.
It has recently been proposed to make use of the pointwise roughness characteristics of a turbulence velocity time series, u(t), to isolate the periods of high activity in environmental turbulence data (Keylock, 2008b(Keylock, , 2009) ) using Hölder/Lipshitz exponents, α u (t).That is, the differentiability of a signal relative to polynomial approximations within the local domain of a specific point are used to derive α u (t).Hence, studying u(t) in a neighbourhood, δ, about a position, T , and taking t and T to be rescaled over the unit interval (ranging from 0.0 to 1.0), we obtain from a Taylor series expansion: where m is the number of times that u is differentiable in T ± δ. Defining the error in approximating u(t) at T by p T (t) as means that the order of differentiability of u(t) close to T gives an upper bound on T (t): This upper bound is then given by a non-integer Hölder/Lipshitz exponent, where a function u(t) has a pointwise Hölder exponent, α u ≥ 0 if a constant K > 0 and the polynomial p T (t) of degree m exists such that The Hölder regularity, α u (t), of u(t) at T is then given by the supremum of β that fulfil Eq. ( 12).
This approach was discussed by Kolwankar and Lévy Vehel (2002) and considered to be more accurate than wavelet based methods.For the comparison of the differential sensitivity of metrics to the presence of multifractality (Sect.6.2), using correlation, non-linear association by phase synchronisation, and on the variance of the estimated α u (t), as described below), we require a precise and consistent estimator for the Hölder regularity.
Figure 10 evaluates different methods for calculating Hölder exponents and supports the use of the oscillationbased method.Figure 10a shows the power spectral density for a fractional Brownian motion (N = 4096) and the fitted slope (red line) from this plot is shown as a Hölder exponent in Fig. 10b with a black dotted line.The DWT (see Appendix A) wavelet-based method using a Daubechies wavelet with 8 vanishing moments is clearly the least precise method.The accuracy and precision of the oscillation-based method increases as the size of the bins used to estimate δ increases.The results are much more sensitive to this than the particular method used to fit the log-log regression line.Figure 10c shows a sinusoidal curve in blue that was used to prescribe the variation of α(t) for a multifractional Brownian motion.This signal itself is not shown but attempts to back-estimate the α(t) from this signal are given in black (DWT-based method), green (least-squares oscillation-based method with the bins for δ ranging from 2 1 to 2 5 -O5) and red (least-squares oscillation-based method with the bins for δ ranging from 2 1 to 2 10 -O10).Our choice of the O10 algorithm is the most precise.These calculations were performed using the FRACLAB toolbox (FRACLAB, 2006).

The differential sensitivity of particular metrics to the presence of multifractality
For a choice of ρ = 0, the surrogate series will preserve the Fourier spectrum to some error level, but not the intermittency (multifractal characteristics).As ρ → 1, the variance of the Hölder series tends towards that for the original data.
Our gradual wavelet reconstruction of the turbulence data in Fig. 1 is based on 19 surrogates and compares the sensitivity of five metrics to the presence of (multifractal) nonlinearities: correlations between the velocity series, R(u 1 u 2 ), and Hölder series, R(α u 1 u 2 ), the phase synchronisation between these respective series, γ * S (u 1 u 2 ), and γ * S (α u 1 u 2 ), and the average standard deviation of the Hölder series across the two components: Phase synchronisation is a nonlinear method of association between data series and the procedure we adopted for its calculation is given in Appendix B. These results are shown in Fig. 11.The surrogate tests show that the two velocity series in this region of turbulent mixing appear independent as no significant difference between data and surrogates occurs for either their linear correlation or their phase synchronisation.This also shows that, as expected, linear correlation does not contain information on www.nonlin-processes-geophys.net/17/615/2010/ Nonlin.Processes Geophys., 17, 615-632, 2010 the multifractal characteristics of the signal.However, when the Hölder series are examined, both of these two measures show significant differences for ρ < 0.40.The null hypothesis is rejected at all our choices for ρ using the average standard deviation measure, indicating that this is the most sensitive to the multifractal characteristics of the two series.Our technique permits the relative sensitivity of different metrics to be determined empirically for particular data and it is interesting here that the sensitivity of the correlation and its nonlinear, phase synchronisation counterpart appears to be approximately the same.
From Fig. 11e it follows that intermittency in the surrogates has yet to converge on the data at the highest choices for ρ used in this paper.Part of the reason for this is the broad-band nature of the turbulence signal compared to the Lorenz attractor, for example.Table 1 has typical values of 57% and 68% of wavelet coefficients fixed for ρ = 0.99 and ρ = 0.999, respectively, while the equivalent percentages are 26.1% and 40.1% for u 1 , and 17.4% and 30.8% for u 2 .Increasing ρ to 0.9999 still only fixes 55.1% and 45.1% of the coefficients, respectively.This is why a sparse, wavelet representation of a turbulence signal is such an effective descriptor (half the coefficients contain 99.99% of the energy).
Figure 12 shows, for two wavelet scales, the modulus maxima of the wavelet coefficients for the data (black) and the difference between this series and that for surrogates generated at ρ = 0.99 (red) and ρ = 0.999 (blue).Note that while the values in Fig. 12d are roughly double those in Fig. 12a, the errors are an order of magnitude lower.Hence, for both choices for ρ there are a number of features at the finest wavelet scales whose energy is too small to be fixed, yet which contribute actively to the singularity structure of the time series, affecting the values for the Hölder series in the surrogates and the value for σ av for these data.
However, it is also the case that: 1. Determining the multifractal properties of a signal is a difficult task (Lux, 2004;Seuret, 2006); 2. The variance measure given by Eq. ( 14) is more dependent on the absolute accuracy of our method for evaluating Hölder exponents than the other metrics used in Fig. 11; and, 3. Different realisations of a stochastic multifractal process will lead to intrinsic variability in the estimated values for α(t).
Hence, the fact that Fig. 11e indicates that the multifractal characteristics are only preserved as ρ → 1 may be due to the nature of these particular data, or may be an artefact of the finite precision of the O10 algorithm.A comparison with another multifractal dataset helps to interpret these results appropriately.

The standard deviation of Hölder exponents for multifractal data
We employed a wavelet-based algorithm for generating multifractal data due to Benzi et al. (1993).This is a stochastic algorithm based on a discrete wavelet transform, whereby wavelet coefficients are assigned and then the inverse wavelet transform is used to construct the time series.we follow Benzi et al. (1993) and take an initial, arbitrary coefficient, χ 0,0 , representing the J + 1 wavelet scale, and then form the wavelet coefficients at scales j = J,...,1, hierarchically according to the recursion: where k = 1 2 k, j,k takes the values ±1 with equal probability, and the random variable η j,k here takes the values 2 −5/6 or 2 −1/2 with probabilities of 0.875 and 0.125, respectively.Fifty datasets were generated and two example realisations are shown in black in Fig. 13.
To reduce the variability in the data that would contribute to changes in the variance of the α benzi values for the 50 datasets, the stationarity of each realisation was improved by setting MODWT approximation and the detail coefficients at j = J − 2,...,J to zero (the red signals in Fig. 13).Thus, variability at scales that are affected by the finite length of the record (2 12 values) was removed.Furthermore, to eliminate any problems due to end effects, the α(t) values were calculated over the central 2 11 values.The degree of variability for σ (α benzi ) for all 50 realisations based on this procedure is given by the left-hand boxplot in Fig. 14.A gradual wavelet reconstruction for the dataset with the median value for σ (α benzi ) (shown in Fig. 13a) is then undertaken on the right-hand side of Fig. 14.
The first thing to note is that the multifractal properties of the time series are recovered some way before the limit of ρ → 1.Hence, the result from Fig. 11e is not generally true for all multifractal datasets.For the Benzi et al. (1993) process, IAAFT surrogates are able to match the values for σ (α benzi ) and, working from the right, the surrogates become significantly different to the dataset at ρ = 0.95.Given that we took care to minimise variability in the calculation of σ (α benzi ) resulting from finite size effects, it is also noteworthy that the variability for the 50 realisations on the left of Fig. 14 is greater than for even the IAAFT surrogates.Hence, the O10 algorithm used here would appear to be sufficiently precise both in absolute terms (Fig. 10c) and relative to the intrinsic variation of similarly generated, sotchastic, multifractal time series.Thus, the analysis of the turbulence dataset in Sect.6 shows that σ av is the measure most sensitive to the presence of multifractality, but the observation that even when ρ = 0.999 the surrogates differ from the data is not true in general for all multifractal series.I.e.gradual wavelet reconstruction can mimic relevant properties of multifractal data at values for ρ < 1.

The precision of our technique for evaluating Hölder exponents
The analysis in Sect.6.3 implies that we can use gradual wavelet reconstruction to study the precision of our technique for evaluating α(t).Multifractional Brownian motions were generated similar to those shown in Fig. 10c   where t ranges from 0 to 1 and contains N = 2048 values, and the amplitude, f ∈ {0.005,0.01,0.015}.Hence, the mean value for α u equates to that for inertial range turbulence.As f → 0 we should reach a value where imprecision in our algorithm means that we cannot discriminate between a truly (but weak) multifractional Brownian motion and a fractional Brownian motion with a Hurst exponent of 0.33.That is, the values for the metric applied to the IAAFT surrogates do not differ significantly than that for the data.The Hölder characteristics of the derived time series and their surrogates were evaluated using the O10 method.From Fig. 15, it is clear that as f increases, σ (α x ) for the data (the dotted blue lines) also increases, as expected.Moving leftwards from the right-hand side of these plots a significant difference emerges in Fig. 15a at ρ = 0.95, while it occurs at ρ = 0.97 in Fig. 15b and c.Both of these values are lower than seen in Fig. 11e, while the results in Fig. 15a are very similar to those in Fig. 14.The variability when ρ = 0 is relatively constant for varying f , meaning that, because σ (α x ) is lower in Fig. 15a, there is no significant difference between the observed value and those for IAAFT surrogates when f = 0.005.The complex nature of Fig. 15a (and Fig. 14) shows that no significant difference occurs when the surrogates are unconstrained, owing to the relatively large inherent variation in the time-series and the relatively weak expression of the multifractality.However, as ρ ∼ 0.9, sufficient energy has been fixed in place for the surrogates to be a trained upon the original data, but with such random variability that σ (α x ) is too low.It is only by ρ = 0.97 that sufficient energy is fixed for no significant difference to occur again.In contrast, the degree of multifractality in Fig. 15b and c is sufficient for the gradual wavelet reconstructions to be more simply structured, meaning that the precision of our O10 algorithm is sufficient to detect variability in Hölder exponents at a value for 0.005 < f < 0.01 and above.

Conclusions
This paper has presented a methodology for exploring properties of nonlinear time series through the systematic varying of an energy threshold and the construction of surrogate datasets that conform to this threshold using a wavelet transform.For a given threshold, either one realisation can be obtained based on the wavelet coefficients fixed at that threshold, or the unfixed coefficients can be added back to the wavelet template in an appropriately constrained, stochastic fashion using the IAAFT algorithm to give multiple realisations.Comparing the value of a metric with the values for the wavelet reconstructed series at multiple choices for ρ permits certain properties of the signal or the metrics to be elucidated including: 1. the parts of the signal that need to be preserved to give a value for the metric similar to the original data (laser data example in Sect.4); 2. assessing how closely model results match the value for the metric for the data (Sect.4); 3. classification of time series complexity (Lorenz equations example in Sect.5); 4. the sensitivity of different metrics (turbulence example in Sect.6); and, 5. the precision of a method for generating Hölder exponents (Sect.7).
With regards to multifractal characteristics of geophysical data, an alternative perspective on the methodology presented here would be to re-define the analysis such that the starting point is not a linear variant of the original signal, but a multifractal variant of the original signal.Recently, Palus (2008) has proposed a pertinent multifractal surrogate generation algorithm.This development would see the work presented herein take a new direction.In this paper, departures from monofractality of increasing strength will correspond to increases in the value for ρ at which significant departures are first detected (Sect.7).Imposing the multifractal structure at the start would be potentially of interest if one wished to see how a metric that was not directly related to the multifractality of the signal (e.g.Eq. 5) was conserved by the surrogates as a version of ρ increased but with the multifractal spectrum fixed in addition to the Fourier spectrum and histogram of values in the data, as is the case in this study.This alternative version of gradual wavelet reconstruction will be explored in the future.
account for this, data can be constructed from one of the time series before the phase differences are calculated.The mean value of γ for some finite number of phase-shuffled realisations (usually ∼ 10, which is the case in this paper), γS , can then be used to normalise the value of γ calculated from the data according to

Fig. 1 .
Fig.1.Two perpendicular components of turbulent velocity data (u 1 in black and u 2 in grey) obtained in a wind tunnel at 5000 Hz.

Fig. 2 .
Fig. 2. Illustration of the surrogate generating algorithm for u 1 with ρ = 0.5 just showing MODWT scale j = 8.The w j =8,k=1:N at this scale from the original MODWT are shown in (a), while (b) shows the fixed coefficients together with randomised values added to a cubic Hermitian polynomial that is used to interpolate between fixed values.The coefficients after one full iteration of the IAAFT algorithm has been applied to the coefficients in (b) are shown in (c) and the results after convergence (indicated by cv) are in (d).The difference between (a) and (d) is indicated by the prime and is given in (e), while (f) shows the differences for another realization of (b) after convergence of the IAAFT algorithm.

Fig. 3 .
Fig. 3. Surrogate time series for the u 1 dataset from Fig.1are shown in the left hand column for stated values of ρ, while the right hand column gives the corresponding absolute part of the continuous wavelet transform for each surrogate series using a Mexican Hat wavelet applied to the first 5 scales (formed using 81 voices).Differences would be less visible if more scales had been displayed because the higher energy at higher j means that a greater proportion of coefficients are fixed.

Fig. 4 .
Fig. 4. Example realisations of the Santa Fe laser intensity data, x L , for various choices of ρ. Figure 4a is the original data series.The other 5 time series are those for the surrogate with the median value for A λ=1 s at the appropriate value for ρ in Fig. 5.The values for ρ are 0.00 (3b); 0.20 (3c); 0.50 (3d); 0.70 (3e); 0.97 (3f).

Fig. 5 .
Fig. 5. Gradual wavelet reconstruction of the laser intensity data based on the one step time asymmetry, A λ=1 .The dotted line gives the value for A λ=1L .The boxplots are the values for the 19 surrogates, where the box indicates the upper and lower quartiles, the central line is the median and the whiskers extend up to 1.5 times the interquartile range, with outliers indicated by a +.Note the nonlinear scale on the abscissa.
Fig. 6. x L is shown in (a),while (b-d) illustrate the difference (indicated by a prime) between this time series and various surrogates series at 3 choices for ρ (0.97 in 6b, 0.95 in 6c, and 0.50 in 6d).The black line in these cases represents x L −x F , where x F is a data series produced solely from the fixed wavelet coefficients.The grey lines show x L x s where the upper surrogate series maximizes A λ=1 s at this ρ in Fig.5, and the lower series has the median value for A λ=1 s .

Fig. 7 .
Fig. 7.Output from the two models for the laser series, formed by squaring the y 1 output from Eqs. (6) and (7), which may be compared to the series, x L , in Fig.4aand a surrogate with ρ = 0.97 in Fig.4f.

Fig. 8 .
Fig. 8. Subsections for time series from the Lorenz equations for different choices of R.

Fig. 9 .
Fig. 9.The correlation dimension, D c , as a function of embedding dimension, D e , for the Lorenz system and surrogates with R = 24.29 (a), R = 24.75(b), R = 28.00(c), and R = 167.0(d).The results for the Lorenz system are given by a black line, with those for series reconstructed from the fixed wavelet coefficients for ρ = 0.99 (dotted line with circles) and ρ = 0.999 (dotted line with squares) and surrogate series for ρ = 0.99 (blue) and ρ = 0.999 (red).The latter summarise results for 19 surrogates according to mean (central horizontal line) and

Fig. 10 .
Fig.10.Testing the oscillation-based method and a wavelet-based algorithm for estimating Hölder exponents.(a) Shows the spectrum of a fractional Brownian motion and a fitted slope, which is translated into a value for α in (b) (dotted line).The box and whisker plots in (b) show the median, upper and lower quartiles (boxes) values extending up to 1.5 times the interquartile range (whiskers) and outliers (red crosses) for the set of 4096 estimated values α(t) as a function of the algorithm used.A wavelet method using a Daubechies wavelet with 8 vanishing moments is indicated by Dau8, while "O" indicates the oscillation method.The number following the "O" is the largest exponent used in the calculation for δ.Hence, "O8" means that the bins for δ ranged from 2 1 to 2 8 .Least-squares regression was used unless a suffix "L" (Lim Inf regression) or "P" (penalised least-squares regression) is included.See(FRACLAB, 2006)  for more details.(c) Shows the sinusoidal function for α(t) used to generate a multifractional Brownian motion (blue line) and then estimated values for α(t) from that multifractional Brownian motion time series.The black line is the Dau8 algorithm, the green line is O5 and the red line is O10.

Fig. 11 .
Fig. 11.Gradual wavelet reconstruction for the Hölder series of the two components of the velocity data from Fig. 1.The measures used are correlation between the velocity series (a) and Hölder series (c), the phase synchronisation between the velocity series (b) and the Hölder series (d), and the average standard deviation of the Hölder series (e).

Fig. 12 .Fig. 13 .
Fig. 12. Modulus maxima of wavelet coefficients for the time series for u 1 are shown in (a) and (d) at j = 1 and j = 4, respectively.The data in (b) and (c), as well as (e) and (f) show the difference (indicated by a prime) between |w j,k | and the values from a wavelet decomposition of a surrogate data series (for ρ = 0.99 in (b) and (e), and shown in red; for ρ = 0.999 in (c) and (f), and shown in blue).

NonlinFig. 14 .
Fig. 14.Gradual wavelet reconstruction of the dataset giving the median value for σ (α benzi ) based on fifty realisations of the multifractal process is shown in the right-hand figure.The left-hand figure is a boxplot giving the inherent variability of σ (α benzi ) for the 50 datasets.

Table 1 .
The proportion of wavelet coefficients fixed in place as a function of R for our two choices of ρ.