A multifractal surrogate data generation algorithm that preserves pointwise H¨older regularity structure, with initial applications to turbulence C. J. Keylock Sheffield Fluid Mechanics Group and Department of Civil and Structural Engineering, University of Sheffield, Mappin Street, Sheffield, U.K., S1 3JD ∗ (Dated: January 4, 2017) An algorithm is described that can generate random variants of a time series or image while preservingtheprobabilitydistributionoforiginal valuesandthepointwiseH¨olderregularity. Thus, it preserves the multifractal properties of the data. Our algorithm is similar in principle to well- 7 known algorithms based on thepreservation of theFourieramplitude spectrum and original values 1 ofatimeseries. However,itisunderpinnedbyadual-treecomplexwavelettransformratherthana 0 Fouriertransform. Ourmethod,whichwetermtheIteratedAmplitudeAdjustedWaveletTransform 2 (IAAWT)method can beused togenerate bootstrapped versions of multifractal data and,because n itpreservesthepointwiseH¨olderregularitybutnotthelocalH¨olderregularity,itcanbeusedtotest a hypothesesconcerningthepresenceofoscillatingsingularities inatimeseries,animportantfeature J ofturbulenceandeconophysicsdata. Becausethelocations of thedatavaluesarerandomized with 3 respect to the multifractal structure, hypotheses about their mutual coupling can be tested, which is important for thevelocity-intermittency structure of turbulenceand self-regulating processes. ] n PACSnumbers: 05.45.Tp,05.45.Df,02.50.Fz,47.27.-i,07.05.Pj a - a t a d . s c i s y h p [ 1 v 9 7 5 0 0 . 1 0 7 1 : v i X r a ∗ c.keylock@sheffield.ac.uk 2 I. INTRODUCTION Givenameasurementofacomplex,andpotentiallynonlinearphenomenonintheformofanimageoratimeseries, it is often important and insightful to test the statistical significance of metrics that attempt to characterize and summarize the underlying data structure. Because a complex process does not necessarily have a Gaussian or even a known prior distribution function, conventional statistical analyses are not always readily applicable. Hence, the surrogatedataapproachwasdevelopedwithinnonlinearphysicsinthe1990sasameanstotreatthisclassofproblems [1–5]. This has proven to be a popular methodology, with applications spanning many disciplines, including medical physics [6–8], fluid mechanics [9, 10] and the geosciences [11, 12]. Perhaps the best known algorithm is the Iterated Amplitude Adjusted Fourier Transform (IAAFT) method [3], which gives surrogates constrained to the original values in the data, while also preserving the Fourier amplitude spectrum accurately. Consequently, differences between the data and surrogates will reside in a specific structure to the Fourier phases of the original data that cannot be replicated by chance in the synthetic surrogate data. Consequently, hypotheses concerning the possible nonlinear properties of an observed time series may be answered by comparing the values of some metric that characterises the nonlinearity (e.g. the correlation dimension, maximal Lyapunovexponent, derivative skewness,nonlinear predictionerror)[13] between the data and the surrogates. Thus, havingselectedastatisticalsignificancelevelαandassumingatwo-tailedstatisticaltest, sothat 2 1surrogatesare α− generated,if the value ofthe metric forthe datafalls outside the rangefor the surrogatesthena significantdifference exists. A development of this general approach is termed Gradual Wavelet Reconstruction (GWR) and envisages the originaldata to be atone endof acontinuum (ρ=1)andphase-randomised,IAAFT data to be atthe other(ρ=0). Surrogate data along this continuum are generated by fixing in place some of the phase information in the original time series [14]. This is achieved by adopting a wavelet transform and, indeed, wavelets and windowed Fourier transformmethodshavebeenusedinthepastasameanstogenerateIAAFT-likesurrogates[15]andforconstraining the randomization in other ways [16, 17]. With surrogates generated at different values for ρ it is then possible to generate bootstrapped confidence limits for data that preserve a sufficient degree of the nonlinearity present in the time series [18, 19], or to test relative degrees of complexity on a given metric by determining the value for ρ above which there is no longer a significant difference between data and surrogates [20]. Given that one way in which a time series may depart from its IAAFT surrogates is through the presence of multifractality, a different approach to GWR for this case is to construct surrogates that preserve the multifractal characteristics of a time series [21]. Such surrogates can be used to generate bootstrapped variants of multifractal phenomena such as human cardiac dynamics [8], earth surface topography [22], precipitation [23], stock market price fluctuations [24], or Quantum Hall systems [25]. They may also be used as a means for testing if the observeddegree of multifractality departs significantly from that of a constant pointwise Ho¨lder exponent. However, the existing multifractal algorithm of Palu˘s [21] does not impose the additional constraint that the values of the original dataset is preserved,making the algorithm more akin to the original surrogate data methods of Theiler and co-workers[1, 2] rather than the IAAFT algorithm, where differences have to reside in the phase information. This paper gives a new algorithm that not only incorporates this feature, but also approximates the ordering in time (or space for the image processing variant) of the Ho¨lder exponents for the original data. We demonstrate the utilityofthemethodbyapplyingittothedetectionofoscillatingsingularitiesintimeseries[26,27],whichisrelevant tofurtheringourunderstandingofturbulence[28,29]andtoanalysingturbulentsystemswherethereisadependence between the intermittency and the macroscalevelocity, as originally discussed by Kolmogorov[30], and subsequently shown to be of relevance in various contexts [31–33]. Before the new algorithm is explained and applied, we briefly review essential concepts of multifractality and Ho¨lder exponents, leading to methods for their estimation. We then review surrogate algorithms for IAAFT surrogate generation as well as the method introduced by Palu˘s [21]. After presentationofour method, we test its relativeeffectiveness andthen providesome example applicationsofrelevance to turbulence physics. II. HO¨LDER REGULARITY A. Definition of pointwise and local H¨older regularity In this paper we primarily focus on the pointwise Ho¨lder regularity, α (t), of a function, the variation in which is p approximatedbymethods forcharacterizingmultifractalitysuchasstructurefunctions[34,35],universalmultifractal exponents [36] or wavelet transform modulus maxima methods [37, 38]. However, in addition, one may define the local Ho¨lder regularity, α (t) to characterize oscillating singularities [39, 40], which are not distinguished from sharp L singularities when using α (t), and can go further and define a function of such exponents calculating all forms of p 3 regularityata point [40]. We makeuse ofα (t) when looking atthe statisticaldetection of anoscillatingsingularity, L as described below. Given a signal, u , a position, t , and assuming that α 0,...,1 then we seek values, β , for which t 0 p p ∈{ } u(t) u(t ) c t t βp (1) 0 p 0 | − |≤ | − | The value for α (t ) is then the supremum of this set of legitimate β values. Following Seuret and L´evy V´ehel [39], p 0 p to define α , we commence by considering a ball, B(t ,R), with a centre at t and with radius, R. If the constant, L 0 0 c , exists such that the following inequality holds: L u(t ) u(t ) c t t βL (2) i j L i j | − |≤ | − | where t ,t (i = j) fall within B(t ,R), then the choice for β is legitimate and α (t ,B(t ,R) is given by the i j 0 L L 0 0 6 supremum of these β . Finally, the local Ho¨lder exponent, α (t ) is given by L L 0 α (t )= lim α (t ,B(t ,R)) (3) L 0 L 0 0 R→0 and α α . L p ≤ This formulation permits us to contrast two functions - a cusp and a chirp - or in the terminology of Hunt and Vassilicos [28], a point vortex and a spiral vortex: Cusp: u(t)= t t α+ (4) 0 | − | 1 Chirp: u(t)= t t α+sin (5) | − 0| t t0 αω | − | In the former case, at t =0, α =α α . In the latter, α =α and α =(α /α ) 1. Hence, α measures the p L + p + ω p L p ≡ − density of singularities, but is insensitive to their type. On the other hand, α provides information on the topology L of these singularities. B. Combined estimation of α and α in a function space p L There are a number of techniques that exist for estimating α (t) [40] and α (t) [26, 27], but if both are to be p L determined simultaneously it is important that the methods used are consistent. This can be assured using a time domain, microlocalfrontier method [40, 41] that resolves both α (t ) and α (t ) simultaneously in the time domain. p 0 L 0 ′ (s,s ) AfunctionissaidtobelongtoK ,iffor0<δ <1/4,c >0,andforall(t ,t )where t t <δ and t t <δ, t0 k i j | i− 0| | j− 0| the following is true: ′ ′ u(t ) u(t ) c t t s+s (t t + t t )−s /2 i j K i j i j i 0 | − |≤ | − | | − | | − | ′ (t t + t t )−s /2 (6) i j j 0 × | − | | − | ′ (s,s ) ′ A comparison to (1) and (2) shows that K includes the definitions for both α and α . Setting s = 0 recovers t0 p L (2) and in the limit δ 0, one obtains s α (t ), meaning that the supremum is α (t ). Because the method is L 0 L 0 → ≤ ′ concernedwith first order differences (increment statistics), both s and s have a maximum of 1. (Note that higher − order singularities can be studied by replacing the left hand side of (6) with a higher order difference [41]). Hence, it follows that s+s′ = 1 is an upper bound on the meaningful behavior of K(s,s′). Forming the ratio between the left hand side of (6) as the numerator and the right hand side as the denominator, it is clear that this ratio tends to ′ ∞ if s+s >α (t ) for fixed t and increasing t . Given that α (t ) α (t ) it follows that the value for s that marks p j j i L 0 p 0 ′ ≤ the intersection between values on the convex frontier for K(s,s ) and s+s′ =0 gives α (t ). Similarly, the value for t0 p 0 ′ s where K(s,s ) intersects s′ = 0 gives α (t ) [40]. This is the method adopted here for the simultaneous estimation t0 L 0 of the pointwise and local exponents [42]. III. SURROGATE DATA ALGORITHMS Both the IAAFT method and the multifractal method of Palu˘s [21] are relevant to our approach. The IAAFT algorithm for a discrete time series x , t=1,...,N may be stated as [3]: t 4 1. Store the squared amplitudes of the discrete Fourier transform of x (i.e. X2 = Nx ei2πf(t/N) 2); t f | 1 t | P (j=0) 2. perform a random shuffle of x to give x ; t t 3. Then iterate a power spectrum step and a rank-order matching step on x(j) as follows: t (a) TaketheFouriertransformofx(j) andreplacethesquaredamplitudeswithX2,whileretainingthephases. t f Given the initial randomsort, this means that the spectrum should be preservedbut with randomphases. Invert the Fourier transform with the original amplitudes restored; (b) Replacethevaluesinthenewseriesx(j) bythoseinx usingarank-ordermatchingprocess. Thispreserves t t 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 from the previous iteration. In this way, the histogram of the original data is preserved precisely, and the Fourier spectrum is approximated to a given error tolerance. ThemultifractalapproachPalu˘s[21]isconceptuallysimilartoadiscretewavelettransform(DWT)basedalgorithm forrealisingmultifractaldata[43,44]. Syntheticwaveletcoefficientsaregeneratedonadyadictreethat,wheninverted yields a dataset with prescribed multifractal characteristics. If the chosen distribution function, P(η), is symmetric about zero then, given a value for the wavelet coefficient at the largest scale, w , the coefficients at scale, j, and 0,0 position, k, are found by: w =η w (7) j,k j,k j−1,1k 2 If P(η), is asymmetric about zero, then this is modified to w =ǫ η w (8) j,k j,k j,k j−1,1k 2 where, with equal probability, ǫ = 1. The Palu˘s [21] algorithm commences with a DWT of a dataset and, with j,k ± scales j 0,1 left fixed, then builds surrogate values for η for j >1, designated here by η˜ according to: j,k j,k ∈{ } w j,k η˜ = j >1 (9) j,k w j−1,1k 2 At each scale, j, the 2j multipliers η˜ are then permuted to give values µ˜ . Setting w˜ =w for j 0,1 , the j,k j,k j,k j,k ∈{ } new wavelet coefficients are constructed for j >1 by: w˜ =µ˜ w˜ (10) j,k j,k j−1,1k 2 Amplitude adjustment is then used to recover the original distribution of wavelet coefficients at each scale. That is, at a scale, c, the w˜ are replaced by the equivalent ranked values from w . In this way, the permutations of j=c,k j=c,k the DWT arerandomizedina manner thatpreservesthe originalhierachicalstructure. Fromthis randomizeddyadic structure of the w , the inverse DWT can be employed to generate a new data sequence. Note that the amplitude j,k adjustment step of this algorithm is applied to the w . Consequently, the histogram of the new data sequence will j,k not exactly replicate the values for the original time series. Figure 1showsanexamplemultifractionalBrownianmotion[45,46],m(t),generatedfromanunderlyingsinusoidal function for the pointwise Ho¨lder regularity, α (t), using the algorithm of Wood and Chan [47]: α (t) = 0.33+ p p 0.2sin(6πt) as the black solid line together with two example surrogates from the IAAFT algorithm (black and displaced upwards) and using the Palu˘s algorithm (black and displaced downwards). The IAAFT algorithm returns the correct average roughness because of the Fourier spectrum preservation. However, it clearly fails to capture the intermittency inthe testsignal. The Palu˘salgorithmdoes abetter jobin this respect,but itis clearvisually thatthe superior surrogates are those in gray, which are produced using the algorithm presented in this paper. In addition, while the data values in these latter surrogatesoccur at random locations, the Ho¨lder properties are clearly localised - the roughand smooth patches occur at the same positions in the data and surrogates. This is a desiredproperty of the new algorithm as shown in the example applications, below. 5 4 2 0 s e at g o urr−2 s + m(t) −4 −6 −8 0 0.2 0.4 0.6 0.8 1 t FIG. 1. A multifractional Brownian motion, m(t), is shown as a heavy black line at the origin and is generated from αp(t)= 0.33+0.2sin(6πt). TwoexampleIAAFTsurrogatesareshowninblackanddisplacedverticallyabovem(t),whiletwoexample Palu˘s surrogates are in black and displaced below m(t). Surrogates generated by the new algorithm are shown in gray and immediately below m(t). IV. THE IAAWT ALGORITHM A pair of dyadic wavelet trees may be designed to form a Hilbert transform pair [48] and such wavelets have been proposed for analysis of turbulence [49] and waveform encoding [50]. More specifically, we employ the dual-tree complex DWT [51, 52]. The analytic signal of some real signal, x(t) is given by x (t) = x(t)+ix (t), where x (t) a H H is the Hilbert transform of x(t), which is a convolution operator with a filter given by h(t)=1/(πt): ∞ x (t)= h(τ)x(t τ)dτ (11) H Z − −∞ Because the Fourier transform of h(t) lies completely in the imaginary plane, it follows that a Hilbert transform approach can be used to perform a complex-valued wavelet transform. Given two filters g(t) and h(t) and their FouriertransformsG(ω)andH(ω)thenitmaybeshownthatifG(ω)=H(ω)e−iω/2 for ω <π thentheirassociated | | wavelets form a Hilbert pair [53], which can be achieved for orthogonal wavelets by offsetting the scaling filters by one half sample. Hence, one approach would be to deploy two trees of linear phase filters, of even length in one tree and odd in the other [51]. However,suchfilters lack orthogonalityand the sub-sampling structure is not particularlysymmetric. Thus, Kingsbury [51] proposed the Q-shift dual tree where, below the coarsest scale, all filters are even length, but nolongerlinearinphase. Bydesigningthefiltersto haveadelayof 1 sampleandbyusingthetime reverseofoneset 4 of filters in the other tree, the required 1 sample delay canbe achieved. In this paper we use symmetric, biothogonal 2 filters with support widths of 13 and 19 values for the first level of the algorithm and Q-shift filters with a support of 14 values for all other levels on the dual tree (case C in Kingsbury [51]). The Q-shift dual tree approach retains properties that make undecimated transforms advantageous for use in surrogate generation, such as shift invariance [5], but at a computational cost that is merely double that for a standard discrete wavelet transform, O(N), rather than O(Nlog N), which is the case for an undecimated maximal overlap discrete wavelet transform [54]. 2 Our algorithmmay now be stated ina form thatparallelsthe IAAFT. Given a signaloflength N =2J we proceed as: 6 1. Performadual-treecomplexDWTandobtaintheamplitudesandphasesoverallJ scalesforthek =1,...,2j−1 complex valued w at each j; j,k 2. Randomly sort the original data and take the dual-tree complex DWT to produce randomised wavelet phases for each scale; 3. Produce new w by combining the original amplitudes with the randomised phases; j,k 4. Iterate the following steps until convergence is achieved: (a) Perform the inverse wavelet transform to give a new time series and then apply the same amplitude adjustment step used in the IAAFT algorithm; (b) Take the dual-tree complex DWT and obtain the new phases, combine these with the original amplitudes to give the latest w . j,k Aswiththe IAAFT algorithm,this newmethod, whichwetermthe IteratedAmplitude AdjustedWaveletTransform (IAAWT), preserves the probability distribution of original values in the dataset, with the multifractal structure preserved up to a convergence criterion. V. TESTING THE IAAWT ALGORITHM A. Tests using signals with prescribed regularity structure Figure 1showsthatthe IAAWT algorithmreplicatesthemultifractalcharacteristicsofasignalwhenthe α (t)are p givenby acontinuousanddifferentiablefunction oftime. Amore complexmultifractalsignalhas its α (t) prescribed p byacontinuousandnon-differentiablefunction,suchastheWeierstrassfunction[55],orbysomehierarchicalprocess. We consider the latter by generating a signal using the algorithm of Benzi et al. [43]. This provides a stringent test of the IAAWT algorithm because the similarities between the Benzi et al. [43] and Palu˘s [21] algorithms imply that the latter will yield better results than the IAAWT technique. With reference to eq. (7) and (8), we use an empirical, discrete distribution function for P(η) where η ∈ 0.4,0.5,0.6,0.7 and the corresponding probabilities are P(η) 0.45,0.1,0.1,0.55 , and with w selected ran- 0,0 { } ∈ { } domly from this distribution, we then apply (8) to derive a dyadic tree of wavelet coefficients over 16 scales, which we then invert with a least asymmetric Daubechies wavelet with 3 vanishing moments to produce a signal, B(t) containing 216 values. Nineteen further such signals were generated to characterise the inherent variability in the multifractal properties of stochastic signals with the same P(η). The multifractal properties of these data were de- termined using the structure function approach [34, 35] where, for a choice of moment exponent, n, the structure function is S (r) = [B(t+r) B(t)]n and the scaling exponent, ξ(n) describes the power-law relation between r n h − i and S (r) [34]. n Figure 2 showsanexample signalB in blackin (a)and(b) togetherwith surrogateseries, B/ in gray. Panel(b) is simply a magnification of a relatively rough regionof (a) showing that this section of the signal also contains regions of relatively high and low roughness, demonstrating the hierarchical nature of the multifractality of B. Note that both surrogate algorithms yield multifractal data but the enhanced ability of the IAAWT algorithm to preserve the multifractal properties of the original data is quantified in (c), where the variability for all 20 different variants of B areshownbythecircles(medianvalue)andheavyblacklines( 2standarddeviations). Thestructurefunctionresults ± for the data shown in Fig. 2a are given by a square and structure function estimates for 25 surrogates for this signal using the IAAWT algorithm are shown by the triangle and error bars, and using the Palu˘s algorithm, by a diamond with error bars. Not only is the median for the IAAWT surrogates very close to the value for the data (shown as the square), the error in estimating ξ using the IAAWT is much less than that inherent to different realisations of B. n The Palu˘s algorithm is both less accurate and less precise than our new method. B. Testing that characteristics other than multifractality are not preserved Tohighlightthatouralgorithmis onlypreservingthe multifractalstructure,weexaminethe abilitytopreservethe skewness of the increments and skewness of the absolute part of the velocity increments for a turbulence time series in 3. The data series consists of 217 samples of the longitudinal velocity component, u , measured 0.02 m ( 150 x wall units above the bed of a wind tunnel in a developed boundary layer with a mean flow of 6 m s−1, with∼more detailon the experimentaldesignin the originalreferences[56]). Results arepresentedfor choicesof the separationr 7 (a) 2.2 (c) 0.02 /B 2 nd 0 B a 1.8 −0.02 1.6 n 0 0.2 0.4 0.6 0.8 1ξ t 1.4 (b) −0.005 1.2 /B d −0.010 1 n a B 0.8 −0.015 2 3 4 5 6 0.1 0.12 0.14 0.16 0.18 0.2 t n FIG. 2. A multifractal signal, B, generated using the Benzi et al. [43] algorithm is shown in (a) as a black line. Panel (b) zooms in on a portion of B. Surrogate series generated using theIAAWTalgorithm are shown in gray and displaced upwards and those from the Palu˘s [21] algorithm are in gray and displaced downwards in these two figures. Panel (c) shows structure functionexponents,ξn asafunction ofmomentnumber,n,witherrorbarsindicating±2standarddeviations. Thedarklines with circles are for 20 different realisations of theprocess with constant P(η). The squares give the structurefunction scaling for the data in (a) and (c), while triangles are the results for 25 IAAWT surrogates of this data series, and diamonds are the results for 25 Palu˘s surrogates. Valuesare displaced horizontally from theinteger valueof n for clarity. z = −1.49 z = 4.69 r = 0.022 L 1 1 0.00069 0.00071 0.00073 0.00075 −0.00005 0 0.00005 0.0001 z = 1.42 z = 10.21 r = 0.10 L 1 1 0.0095 0.0096 0.0097 0.0098 0 0.001 0.002 0.003 z = 0.80 z = 10.87 r = 0.22 L 1 1 0.0262 0.0264 0.0266 0.0268 0 0.002 0.004 0.006 0.008 z = 0.62 z = 4.08 L 1 1 0.082 0.084 0.086 0.088 −0.005 0 0.005 0.01 〈 (|u − u |)3〉 〈 (u − u )3〉 x x−r x x−r FIG.3. Ananalysisofthethirdordermomentofthevelocitydifferencesforaturbulencetimeseriesofthelongitudinalvelocity component, ux. The values for r range from the Taylor scale ∼0.02L to the integral scale, L. The circles and squares are the values for the data, while the boxplots indicate the values for the surrogates: the median by a heavy vertical lines, the edges of the box by the lower and upper quartile, the whiskers extend up to ×1.5 the interquartile range, with outliers shown by crosses. that vary from the Taylor scale to the integral scale, L. Ten IAAWT surrogates were generatedto form the boxplots in 3. Because it is (u u )3 (circles), rather than (u u )3 (squares) that is the relevant expression in x x+r x x+r h | − | i h − i the structure function formulation of multifractality [34], it follows that surrogates and data should lie much closer to one another in the left hand column of Fig. 3. The z-scores included in each panel show this is the case - there is no significant difference between data and surrogates for any case on the left hand side of the figure at the 5% level, while all on the right are clearly significantly different at this level. Where difference are seen at the left-hand side, it is at the very smallest scales where there is inevitably some sensitivity to the random nature of the algorthm as differences are being taken over very small distances. 8 VI. EXAMPLE APPLICATIONS There are a great range of potential applications of this algorithm as it provides a means for assessing confidence limits for any algorithm applied to multifractal data and, indeed, to techniques used to determine the multifractal properties of observed data. Here we consider two contrasting applications, both of which are of potential relevance to turbulence physics. A. Oscillating singularities in time series Understanding the nature of any singularity or near-singularity (dissipation potentially smoothing discontinuities) in turbulence is an important area connecting topology to the Navier-Stokes equations [28, 57]. Related phenomena have been documented in studies of financial prices and rupture processes [58, 59]. We have already defined relevant cusp and chirp signals in terms of their pointwise Ho¨lder, α (t), and local Ho¨lder, α (t), exponents in (5), and p L describedthecoupledestimationoftheseexponents. Inaddition, notethatthesetofα (t)valuesisoftenrepresented p by its singularity spectrum, D(α ) [37], but various multifractal spectra can be defined: the Haussdorff multifractal p spectrum [60]; the Legendre transform spectrum [37]; and, the large deviation spectrum [61]. Prominent additional modes in the latter, resulting in a loss of convexity, are indicative that the α are an important part of the signal. L 1. Definition of artificial signals with a chirp Two test case signals are considered in this section: a fractional Brownian motion (fBm) with a chirp; and, a multifractionalBrownianmotion(mBm) witha chirp. They aredefined for 0 t 1.2andconsistof24 000discrete ≤ ≤ samples. The fBm has α = 0.33 (it is turbulence-like to the level of the second order structure function), while for p the mBm we have α (t) = 0.33+0.15sin16πt. Both are discretised into 20 000 values over the support 0 t 1. p The chirp in both cases has a support of 0.1 χ 0.1 discretised over 4000 values, and is given by χα+s≤in ≤1 , − ≤ ≤ | | |χ|αω with α =1.4 and α =0.33,using the notationin (5). Thus, α =0.33and α 0.14. The chirp was insertedinto ω + p L ∼ the fBm or mBm near the center of the signal at a position that eliminated introducing any additional discontinuity (that the detection might be sensitive to). Figure 4 shows these signals in black focusing on the region of the chirp(0.2 t 1). IAAFT surrogatesfor these ≤ ≤ twosignalsareshowningrayinpanels(b)and(d),anditisclearthattheoscillatingsingularityandtheintermittency have been destroyed. In contrast, Fig. 5 shows the signal from Fig. 4c together with three IAAWT surrogates. It is clearthattheIAAWTalgorithmpreservesthestructureofthemBm,althoughthedetailedbehavioroftheoscillation are lost (with only the broad characteristics giverned by α (t) preserved. p 2. Oscillating singularity detection ′ In addition to the Ks,s method discussed in the introduction and as noted above, the presence of oscillating singularities can be detected by a lack of convexity in the multifractal large deviation spectrum [60, 61]. Figure 6 showsthesingularityspectraforthefBm+chirp(toprow)andthemBm+chirp(bottomrow)signals. Inbothcases there is a lack of convexity to the multifractal spectrum and in the latter case the secondary mode is as large as the primary mode centered at α 0.33. The left hand figures are compared to 19 IAAFT surrogates of the respective p ∼ signals in gray, examples of which are given in Fig. 4. it is clear that convexity is recovered, indicating the removal of the oscillating singularity (the high values for α are a consequence of the edges of the chirp being smoother than p α =0.33in a pointwise sense). Incontrast,the IAAWT surrogatesinthe righthand panels are clearlyretaining the p general α characteristicsof the oscillation,with no significantdifference existing between the D(α for the data and p p surrogates. ′ Giventhatthe generalcharacteristicsofthe oscillationareretainedinthe IAAWT surrogates,wecanusethe Ks,s frontiers (section IIB) to undertake a more detailed analysis and extract the oscillating singularity to be detected with statistical confidence. We estimated values for α (t) and α (t) for the mBm+chirp signal and for 19 IAAFT p L and IAAWT surrogates from the intersection of the Ks,s′ frontier (estimated with 40 points) with s′ = 0 (α ) and L ′ with s+s = 0 (α ). This was undertaken for 1001 positions, centered on the middle of the chirp and traversing p double the supportofthe chirp(i.e. everyeighthvalue wascalculated). Givena position, t , we may derivea z-score 0 9 4 (a) ) (t 0 u −4 4 (b) (1)uS 0 −4 1 (c) ) (t 0 u −1 1 (d) (1)uS 0 −1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t FIG. 4. A fractional Brownian motion with an embedded chirp is shown in black in panel (a). An example IAAFT surrogate is shown in gray in (b). Panel (c) shows a multifractional Brownian motion with an embedded chirp is shown in black. An IAAFT surrogate for this is shown in (d). The panels show a smaller segment of the full dataset in each case, which has a support of 0≤t≤1.2. 1 (a) (t) 0 u −1 1 (b) (1)uS 0 −1 1 (c) 2) (uS 0 −1 1 (d) (3)uS 0 −1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t FIG. 5. A multifractional Brownian motion with an embedded chirp is shown in black in panel (a). Three example IAAWT surrogates are shown in gray in the lower panels. The preservation of the pointwise H¨older exponents and the location of the oscillating components is clear. The panels show a smaller segment of the full dataset in each case, which has a support of 0≤t≤1.2. estimate for the pointwise and local Ho¨lder exponents at this point. Using α (t ) as an example, p 0 α (t ) [α (t )] p 0 p 0 S z (t )= − (12) αp 0 σ[α (t )] p 0 S where the S subscript indicates the operation is taken over the 19 values from the surrogate data, the overbar is the arithmetic mean and σ is the standard deviation. Following Donoho and Johnstone [62], an expected upper bound on these z-scores (unit standard deviation by definition) is (2logN/√2)0.5. The results for z (t) (black) and z (t) (gray)are shown in Fig. 7, using IAAFT (a) and IAAWT (b) surrogates. αp αL The expected upper bound on these values is shown by dashed lines on each panel. It is clear from Fig. 4c that the 10 1 (a) 1 (b) ) p α ( 0.5 0.5 D 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 α α p p 1 (c) 1 (d) ) p α ( 0.5 0.5 D 0 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 α α p p FIG.6. Multifractalspectra,D(αp)estimatedusingthelargedeviationmethod(black)andcomparedtosurrogatedata(gray). Panels (a) and (b) are for a fractional Brownian motion with an embedded chirp, while (c) and (d) are for a multifractional Brownian motion with an embeddedchirp. Panels(a) and(c) showIAAFTsurrogate data, whileIAAWTsurrogates areused in(b)and(d). Themultimodalstructureofthespectrumindicatesthepresenceofoscillating singularities intheoriginal data in all cases, but this disappears in the IAAFTsurrogates (although is retained bythe IAAWTsurrogates). 20 (a) T) F A 10 A (IαL z 0 , αp z −10 0.4 0.45 0.5 0.55 0.6 0.65 0.7 T) 5 (b) W A A (IαL 0 z , αp z −5 0.4 0.45 0.5 0.55 0.6 0.65 0.7 t FIG. 7. Time series of z-score estimates based on the the original data and the mean and standard deivation for surrogates for αp (black) and αL (gray). The upper panel is the results for IAAFT surrogates and the lower for IAAWT surrogates. An almost sure threshold [62] is indicated by thedashed lines. center of the chirp lies at t 0.555. Because IAAFT surrogates neither preserve multifractal structure or oscillating ∼ singularities,departuresoutsidetheexpectedupperboundariseatavarietyoflocationsinFig. 7aforbothz (t)and αp z (t). Incontrast,apartfromonedatumatz (t 0.57),the z-scoresforthepointwiseHo¨lderexponentslie within αL αp ∼ the expected bound, which is to be expected for a multifractal preserving algorithm. In contrast, z remains within αL the bound, except at the center of the chirp, where a difference is observed across several consecutive values. Hence, the IAAWT algorithm preserves α (t) characteristics and, because it also fixes in place the position of particular p features,locates achirp-likefeature in the correctplace. However,it does notintrinsicallypreserveα , meaning that L the center of an oscillating singularity is extracted rather precisely.