A Study of Wavelet Analysis and Data Extraction from Second-Order Self-Similar Time Series

Abstract

Statistical analysis and synthesis of self-similar discrete time signals are presented. The analysis equation is formally defined through a special family of basis functions of which the simplest case matches the Haar wavelet. The original discrete time series is synthesized without loss by a linear combination of the basis functions after some scaling, displacement, and phase shift. The decomposition is then used to synthesize a new second-order self-similar signal with a different Hurst index than the original. The components are also used to describe the behavior of the estimated mean and variance of self-similar discrete time series. It is shown that the sample mean, although it is unbiased, provides less information about the process mean as its Hurst index is higher. It is also demonstrated that the classical variance estimator is biased and that the widely accepted aggregated variance-based estimator of the Hurst index results biased not due to its nature (which is being unbiased and has minimal variance) but to flaws in its implementation. Using the proposed decomposition, the correct estimation of the Variance Plot is described, as well as its close association with the popular Logscale Diagram.

1. Introduction

In the past decades, the self-similar processes and long-range dependence (LRD or long memory) have been applied to the study and modeling of many natural and man-made complex phenomena. These kinds of processes have been particularly attractive in the pursuit of optimal design and configuration of network communications.

The published work of Leland et al. in 1993 and 1994 [1, 2] demonstrated that Ethernet traffic is statistically self-similar and that the commonly used models are unable to capture that fractal behavior, highlighting that a burstiness and LRD are present when 𝐻>0.5. Since then, researchers have been studying extensively long memory processes and their impact on network performance, for example, Karagiannis et al. stated that the identification of LRD is not trivial and that not all scenarios in modern networks present LRD characteristics, for example, traffic in the Internet backbone is more likely to be Poisson type instead of LRD [3].

Many researchers have also addressed their studies to determine if network traffic is sufficiently modeled by self-similar processes or a more general model is needed, for example, one that considers multiscaling or multifractality [4–6]. The advantage of the capability to model complex systems with self-similar processes is that the correlation structure is defined by a single parameter: the Hurst index (𝐻).

Unlike other statistics, the Hurst index, although it is mathematically well defined, cannot be estimated unambiguously from real world samples. Several methods have been developed then in order to estimate it. Examples of classical estimators are those based on R/S statistic [7] (and its unbiased version [8]), detrended fluctuation analysis (DFA) [8, 9], maximum likelihood (ML) [10], aggregated variance (VAR) [7], wavelet analysis [11, 12], and so forth. In [13], Clegg developed an empirical comparison of estimators for data in raw form and corrupted. An important observation is that the estimation of the Hurst index may differ from one estimator to another, and the selection of the most adequate estimator is a difficult task. This selection depends greatly on how well the data sample meets the assumptions the estimator is based on. However, through analytical and empirical studies, it has been discovered that the estimators that have the best performance in bias and standard deviation, and, consequently, in mean squared error (MSE) are Whittle ML and the wavelet-based estimator proposed by Veitch and Abry in [11].

From these two estimators, the wavelet based is computationally simpler and faster [7, 11].

In addition to the Hurst index, other statistical characteristics are needed to describe the phenomenon under study. The most common are the first and second-order statistics, that is, mean, variance, and correlation. The classical estimators of these characteristics have been proposed decades ago, for example, Kenney and Keeping demonstrated in 1939 and 1951 that the classical variance estimator is unbiased for independent and identically distributed Gaussian observations [14, 15]. Confidence interval is also given for these estimations, for example, 𝑃[𝑋∈(−3/√𝑁,3/√𝑁)]≈99% (where 𝑋 is the sample mean estimated from a sample of size 𝑁) for a standardized white noise process. This confidence interval is narrower as the sample size increases.

It has been claimed that most processes satisfy those common assumptions [16]. As it has been expressed, however, other authors (Leland et al. [1], Taqqu et al. [5], Tsybakov and Georganas [17], Veitch and Abry [11], and many others) conclude that traffic characteristics present correlation and that the estimation of these statistics (including confidence interval) with the classical estimators (which do not consider correlation) may lead to estimation errors and, consequently, to wrong decisions or inaccurate models, especially when data presents accentuated LRD.

Several estimators of the Hurst index have been proposed, but many of them do not consider the effect of correlation on the estimation of first- and second-order statistics, thus applying incorrectly the classical formulae. Particularly, it has been claimed that the aggregated variance method can only be used as a heuristic method, and that the Variance Plot (also named Variance-time Plot; see Section 2.3) can only be used to check whether the time series is self-similar or not and, if so, to obtain a crude guess for the Hurst index [22, page 44]. This work clarifies this point, demonstrating that the Variance Plot can be estimated efficiently and that the estimation of the Hurst index from it is actually unbiased and has minimum variance (similarly to the wavelet-based estimator).

This work is motivated by the mentioned importance of the self-similar processes in many areas, especially in the analysis and modeling of Internet traffic, and by the fact that there are still some misunderstandings and bad practices that must be overcome.

This document is organized as follows. Section 2 defines the discrete time self-similarity and some of its statistical properties. Section 3 describes the proposed set of basis functions and the analysis and synthesis equations, which are used to generate self-similar samples whose Hurst index matches the estimator proposed by Veitch et al. and the corrected version of the variance-based estimator. Section 4 defines statistics for the sample mean and variance of self-similar discrete processes; Section 5 explains how the variance-based estimator has been misunderstood and defines the correct estimator, which coincides with the Haar wavelet based in a particular case. Section 6 presents some simulations and measurements and, finally, Section 7 concludes the work.

2. Self-Similarity

Self-similarity describes the phenomenon where certain properties are preserved irrespective of scaling in space or time. Deterministic self-similarity is clearly exemplified by popular figures as Sierpinski’s triangle or Koch’s snow flake. This form of self-similarity is named scale invariance and makes different scales of the same object undistinguishable. Stochastic self-similarity is not that obvious, it refers to how statistical properties of a stochastic process are preserved under time expansion. Stochastic self-similarity is defined for continuous and discrete time stochastic processes.

Self-similarity (either continuous or discrete time) is tightly related to short- and long-range dependencies (SRD and LRD, resp.). The degree of self-similarity is defined and measured through the so named Hurst index 𝐻 (0<𝐻<1). It is known that processes with 𝐻<0.5 are SRD, and processes with 𝐻>0.5 are LRD. If  𝐻=0.5, neither SRD nor LRD are present. For example, the commonly used white Gaussian noise (WGN) has always 𝐻=0.5 and does not present any time dependency.

Processes with LRD are also named long memory, as current and future realizations of these are strongly correlated. The dividing line between SRD and LRD processes is not ambiguous; for LRD processes, the autocovariance function is not absolutely convergent (i.e., the sum is not finite), while it is for SRD processes. This work refers only to stochastic discrete time self-similarity.

2.1. Discrete Time Self-Similarity

The definition of discrete time stochastic self-similarity is given in terms of the aggregated processes. Let {𝑋_𝑡;𝑡∈ℕ} be a discrete time series derived from a self-similar process with stationary increments and Hurst index 𝐻 (𝐻-SSSI). The aggregated time series, derived from 𝑋_𝑡 is the sequence given by [17] 𝑋(𝑚)=𝑋_𝑘^(𝑚);𝑘∈ℕ, (1) where each term 𝑋_𝑘^(𝑚) is defined as 𝑋_𝑘^(𝑚)=1𝑚^𝑘𝑚_{()𝑖=𝑘−1𝑚+1}𝑋_𝑖,𝑘∈ℕ, (2) where 𝑚 represents the aggregation level. That is, each new time series is obtained by partitioning the original time series into nonoverlapping blocks of size 𝑚 and then averaging each block to obtain its respective values.

Let 𝑋_𝑡 be a covariance stationary discrete time series with mean 𝜇_𝑋=0, variance 𝜎²_𝑋 and autocovariance function (ACvF) 𝛾_𝑋(𝑘) and 𝑋_𝑘^(𝑚) its aggregated series. Then it is said that 𝑋_𝑡 is self-similar (𝐻-SS), if the following holds [18]: 𝑋_𝑘^(𝑚)∼𝑚^𝐻−1𝑋_𝑡, (3) where ∼ means equality in distribution.

The many methods for generating artificial discrete time self-similar sequences are classified in sequential and fixed length. All these have particular advantages and shortcomings associated to accuracy, generation time, memory or processing resources, and so forth. Sequential generators are sometimes more appropriate for long duration simulations, but the level of approximation and several parameters are needed, while fixed length only depends on the desired Hurst index but may need to store part or all of the sequence in memory. For the purpose of the simulations described in Section 6, a fixed-length Davies-Harte fractional Gaussian noise (FGN) generator is used. FGN is a special type of noise where the autocovariance function has a special shape (expression (8)). Section 3.2 describes a theoretical method to synthesize self-similar discrete sequences using Haar wavelet-based decomposition. Although Daubechies wavelet- (DW-) based approximations are better in mean and variance, the Haar-based method generates a self-similar sequence of specified Hurst index from practically any given sequence, regardless of whether it is self-similar or not if some conditions are met [19]. For other applications, DW are preferred; see [20] as an example.

2.2. Properties of Self-Similar Discrete Time Series

The definition of discrete stochastic self-similarity (3) has important implications about the stochastic process 𝑋_𝑡; these implications include the following properties [21].(i)Zero-mean: 𝐸𝑋_𝑡=𝐸𝑋_𝑘^(𝑚)=0. (4)(ii)Power law of the 𝑞th order moments: 𝐸𝑋_𝑘^(𝑚)𝑞=𝑚^{𝑞(𝐻−1)}𝐸𝑋_𝑡^𝑞. (5)(iii)Power law of the 𝑞th order absolute moments: 𝐸𝑋_𝑘^(𝑚)𝑞=𝑚^{𝑞(𝐻−1)}𝐸𝑋_𝑡^𝑞. (6)

2.3. Second-Order Discrete Self-Similarity

The second-order definition of self-similarity is derived from (5) for 𝑞=2. The variance of the aggregated time series is defined by the following [17]: var𝑋_𝑘^(𝑚)=𝑚^2𝐻−2var𝑋_𝑡. (7)

An equivalent definition is 𝛾_𝑋^(𝑚)(𝑘)=𝜎²_𝑋2(𝑘+1)^2𝐻−2𝑘^2𝐻+(𝑘−1)^2𝐻,𝑘≥0, (8) where 𝛾_𝑋^(𝑚)(𝑘) is the autocovariance function of 𝑋_𝑘^(𝑚).

If a discrete time series 𝑋_𝑡 satisfies these conditions, it is called second-order self-similar with Hurst index 𝐻 (𝐻-SOSS). Note that the mean of an 𝐻-SOSS process is not necessarily zero.

The plot log[var(𝑋_𝑘^(𝑚))] versus log(𝑚) is known as Variance Plot. It is a straight line of slope 2𝐻−2 for self-similar processes. This plot is the basis of the variance-based estimator of the Hurst index. It has been “shown” in the literature that the variance-based estimator underestimates the Hurst index and that the variance-based estimator throws a coarse estimation of the true Hurst index. Section 5 demonstrates that this is a consequence of inadequate implementations of this estimator, that is, the aggregated variance is estimated with the classical formula (36), which is not correct if any correlation exists. An apparent solution is to use the proposed unbiased estimator, but that leads to an ill-conditioned problem: the Hurst index is needed to estimate the variance and vice versa. The solution for this situation is also described (see Section 5.1).

2.4. Wavelet Decomposition and the Logscale Diagram

The wavelet decomposition transforms a signal 𝑋_𝑡 into a sum of orthogonal components as follows: 𝑋_𝑡=^𝐽2_𝑗=𝑖^𝑗_𝑘=1𝑑_𝑋𝑗,𝑘𝜓_𝑗,𝑘(𝑡), (9) where each function 𝜓_𝑗,𝑘(𝑡) is derived from a basis function 𝜓₀(𝑡), namely, the mother wavelet, by scaling and displacement, that is, 𝜓_𝑗,𝑘(𝑡)=2^−𝑗/2𝜓₀2^−𝑗𝑡−𝑘, (10) and coefficients 𝑑_𝑋(𝑗,𝑘) is the value at time 𝑘 of scale 𝑗, computed as an inner product between the signal 𝑋_𝑡 and the wavelet function 𝜓_𝑗,𝑘(𝑡): 𝑑_𝑋𝑗,𝑘=𝑋(𝑡),𝜓_𝑗,𝑘(𝑡). (11)

The statistic 𝑠₂(𝑗) is then defined from these coefficients as 𝑆₂𝑗=𝐸𝑑_𝑋𝑗,⋅², (12) which, for an 𝐻-SOSS process, is related to the Hurst index as 𝑆₂𝑗=𝑐_𝑓𝐶2^{𝑗(2𝐻−1)}, (13) where the quantity 𝑐_𝑓𝐶, related to the power of the process, is considered a constant.

The plot log₂𝑆₂(𝑗) versus 𝑗 forms the widely known Logscale Diagram described by Veitch and Abry [11]. The Logscale Diagram of an 𝐻-SOSS process is a straight line of slope 2𝐻−1. To obtain an unbiased estimation of the Hurst index based on the Logscale Diagram, it is also necessary to subtract the bias that results of averaging the logarithms of the respective variance of real world time series, which is estimated as [11] 𝑔_𝑗=Ψ𝑛_𝑗/2ln(2)−log₂𝑛_𝑗2, (14) where 𝑛_𝑗 is the number of coefficients available at octave, that is 𝐸log₂var𝐶^𝑛,𝑖_𝑋,𝑡≈𝐸log₂var𝐶^𝑛,𝑖_𝑋,𝑡−𝑔_𝑗. (15)

If the Logscale Diagram of a time series cannot be adequately modeled with a linear model, possibly the scaling behavior needs to be described with more than one scaling parameter, that is, the Hurst parameter is not adequate (or insufficient) [22, 23]. However, even if the time series under study is not self-similar, the Logscale Diagram can show whether or not it presents LRD.

3. An Orthogonal Decomposition Performed by Subtracting Aggregated Series

The time series 𝑋_𝑡 can be decomposed into a set of time series, each one defined as 𝐶^𝑛,𝑗_𝑋,𝑡=𝑋^{(𝑛𝑗−1}_𝑡^𝐸)−𝑋^(𝑛𝑗_𝑡^𝐸),𝑛,𝑗∈ℕ, (16) where 𝑋^(𝑛𝑗_𝑡^𝐸) is the time series 𝑋_𝑡 after two operations, which are as follows.(1)Aggregation at level 𝑛^𝑗, as defined by (1) and (2), that is, 𝑚=𝑛^𝑗.(2)Expansion of level 𝑛^𝑗, which consists of  “repeat” each element of a time series 𝑛^𝑗 times, that is, 𝑋^(𝑛𝑗_𝑡^𝐸)=𝑋^(𝑛𝑗)_𝑘 for 𝑘=1+⌊(𝑡−1)/𝑛^𝑗⌋ and 𝑗∈ℕ.

These zero-mean components 𝐶^𝑛,𝑗_𝑋,𝑡 have three important properties.(i)They synthesize the original time series without loss (assuming zero-mean), that is, 𝑋_𝑡=_𝑗𝐶^𝑛,𝑗_𝑋,𝑡. (17)(ii)They are pair-wise orthogonal: 𝐶^𝑛,𝑗1_𝑋,𝑡,𝐶^𝑛,𝑗2_𝑋,𝑡=0,𝑗₁≠𝑗₂. (18)(iii)If 𝑋_𝑡 is exactly or at least second-order self-similar, then the variance of its components satisfies var𝐶^𝑛,𝑗_𝑋,𝑡=𝑟⋅var𝐶^{𝑛,𝑗−1}_𝑋,𝑡, (19) where 𝑟=𝑛^2𝐻−2. (20)

Another useful property relates the variance of the component to the variances of the aggregated series, that is, var𝐶^𝑛,𝑗_𝑋,𝑡=var𝑋^{(𝑛𝑗−1}_𝑡^𝐸)−var𝑋^(𝑛𝑗_𝑡^𝐸). (21)

It is easy to proof (21): from (16), it turns out that var[𝐶^𝑛,𝑗_𝑋,𝑡]=var[𝑋^{(𝑛𝑗−1}_𝑡^𝐸)]+var[𝑋^(𝑛𝑗_𝑡^𝐸)]−2cov[𝑋^{(𝑛𝑗−1}_𝑡^𝐸),𝑋^(𝑛𝑗_𝑡^𝐸)], but as 𝑋^(𝑛𝑗_𝑡^𝐸) is itself an aggregation of 𝑋^{(𝑛𝑗−1}_𝑡^𝐸) it turns out that cov[𝑋^{(𝑛𝑗−1}_𝑡^𝐸),𝑋^(𝑛𝑗_𝑡^𝐸)]=var[𝑋^(𝑛𝑗_𝑡^𝐸)], and (21) comes after a substitution.

Properties (i), (ii), and (iii) imply that 𝜎²_𝑋=_𝑖var𝐶^𝑛,𝑗_𝑋,𝑡,𝜎²_𝑋=1𝐶1−𝑟var^𝑛,1_𝑋,𝑡. (22)

Then, the variance of the 𝑗th component is related to the variance of 𝑋_𝑡 as var𝐶^𝑛,𝑗_𝑋,𝑡=(1−𝑟)𝑟^𝑗−1𝜎²_𝑋. (23)

It is easy to prove the following relation: var𝐶^𝑛,𝑗_𝑋,𝑡=𝑛^−𝑗𝑆₂𝑗. (24)

An immediate consequence of (24) is that the plot 𝑗+log_𝑛[var(𝐶^𝑛,𝑗_𝑋,𝑡)] versus 𝑗 is equivalent to the Logscale Diagram. It is a straight line for 𝐻-SOSS time series, and the slope is related to the Hurst index so that 𝑠=2𝐻−1.

For example, for 𝑛=2 the statistics var(𝐶^2,𝑗_𝑋,𝑡) and 𝑆₂(𝑗) are related as var𝐶^2,𝑗_𝑋,𝑡=2^−𝑗𝑆₂𝑗, (25) and in this case the basis function is the Haar wavelet, which is defined as [24, 25] 𝜓₀(𝑡)=1+1,0≤𝑡<2,1−1,2≤𝑡<1,0,otherwise. (26)

Figure 1 shows the components obtained from an 𝐻-SOSS sample of size 32 and 𝐻=0.9. The squared form of the components is due to the expansion of the aggregated series and disappears after the downsampling.

(a)

(b)

(c)

(d)

(e)

Figure 1 

Components of an 𝐻-SOSS sample of size 32 and 𝐻=0.9.

The authors of  [26] described the aggregation as an inner product with the signal and the Haar “father” wavelet, and then, the relation between wavelet coefficients and aggregation levels is obvious. At this point there are similarities between this section and that previous work, the most important is that the relation between aggregation levels and Haar wavelet is described by Abry et al. However, two differences must be highlighted: (1) the decomposition presented in this work only coincides with that definition for 𝑛=2 in (16); for higher values, the Haar wavelet is not sufficient to describe the components of (16); and (2) authors of [26] discard anyway the estimation of the Hurst index based on the so named “𝑎-aggregation.” In this work, it is clarified that the “𝑎-aggregation” is misunderstood leading to incorrect implementations (i.e., the “classical” variance estimator), and a generalization of the “𝑑-aggregation” is proposed.

3.1. General Waveform of the Basis Functions

The orthogonal decomposition defined by (16) can be expressed in terms of an inner product between the signal and a set of orthogonal wavelet-type functions. Let us describe the waveform for the general case.

The wavelet-based function is 𝜓_𝑛(𝑡)=11−𝑛1,0≤𝑡<𝑛,−1𝑛,1𝑛≤𝑡<1,0,otherwise, (27) and the wavelet functions derived from 𝜓_𝑛(𝑡) are obtained by three operations: scaling, displacement (similarly to (10)), and a phase shift, that is, 𝜓_{𝑛,𝑗,𝑘,𝜃}(𝑡)=1−𝑢𝑛^𝑗(𝑘−1)𝜓_{𝑛,𝑗,𝑘}𝑡−𝜃𝑛^𝑗−1+𝑢𝑛^𝑗(𝑘−1)𝜓_{𝑛,𝑗,𝑘}𝑡+𝑛^𝑗−𝜃𝑛^𝑗−1,𝜓_{𝑛,𝑗,𝑘}(𝑡)=𝑛^−𝑗/2𝜓_𝑛𝑛^−𝑗𝑡−𝑘, (28) for 𝑗=1,2,…,𝐽,𝑘=0,2,…,𝑛^𝐽−𝑗−1, and 𝜃=0,1,2,…,𝑛−1. Note that note that 𝜓_𝑛,0,0,0(𝑡)=𝜓_𝑛(𝑡).

The function defined by (27) is a generalized form of the Haar wavelet. It is always a rectangle-shaped function, but it is not symmetric about the horizontal axis except for 𝑛=2.

Figure 2 shows the basis function for 𝑛=2 without phase shift and with a phase shift of 1/2, respectively. Obviously, ⟨𝑋(𝑡),𝜓_{2,𝑗,𝑘,0}(𝑡)⟩=−⟨𝑋(𝑡),𝜓_{2,𝑗,𝑘,1}(𝑡)⟩, which means that both products give the same information. Redundant information can then be reduced by decimating the sequence of coefficients.

(a)

(b)

Figure 2 

Basis functions: (a) 𝜓2,0,0,0(𝑡) (no phase shift) and (b) 𝜓2,0,0,1(𝑡) (equal to 𝜓2,0,0,0(𝑡) with a phase shift of 1/2).

Figure 3 shows the basis function for 𝑛=3. Note that the phase shift moves the rectangle of height 2/3 from one-third to another. In this case, there exists also redundant information, as 𝜓_3,0,0,2(𝑡)=−𝜓_3,0,0,0(𝑡)−𝜓_3,0,0,1(𝑡) and in the general case 𝜓_{𝑛,𝑗,𝑘,𝜃}(𝑡)=−∑_𝑖≠𝜃𝜓_{𝑛,𝑗,𝑘,𝑖}(𝑡), which means that one of the coefficients, for example, the one obtained with the last phase shift, can be discarded. This means that the sequence of coefficients can be downsampled without loss of information. For a sample of length 𝑁, it is easy to verify that the number of observations that remains in all components (sequences of coefficients) after the downsampling is 𝑁−1, which can be complemented with the sample mean, in the case that this is not zero.

(a)

(b)

(c)

Figure 3 

Basis functions: (a) 𝜓3,0,0,0(𝑡) (no phase shift), (b) 𝜓3,0,0,1(𝑡) (𝜓3,0,0,0(𝑡) with a phase shift of 1/3), and (c) 𝜓3,0,0,2(𝑡) (𝜓3,0,0,0(𝑡) with a phase shift of 2/3).

3.2. Wavelet Synthesis of Self-Similar Time Series

A method to synthesize 𝐻₁-SOSS from practically any time series, regardless of whether it is or not self-similar or its marginal distribution, is proposed. This method consists of adjusting the sum expressed by (9) with a set of weights, that is, 𝑋_𝐻1(𝑡)=^𝐽_𝑗=𝑖𝑤_𝑗^2𝑗_𝑘=1𝑑_𝑋𝑗,𝑘𝜓_𝑗,𝑘(𝑡), (29) where these weights 𝑤_𝑗 are defined as 𝑤_𝑗=𝑐_𝑓𝐶⋅2^{𝑗(2𝐻1−1)}𝑆₂𝑗, (30) where 𝑆₂(𝑗) and 𝑐_𝑓𝐶 are the respective estimations of 𝑆₂(𝑗) and 𝑐_𝑓𝐶 (the associated power parameter [4]) from 𝑋_𝑡 and 𝐻₁ is the desired Hurst index of the new synthetic series. It is necessary that 𝑆₂(𝑗)>0 for all 𝑗=1,…,𝐽.

The weighted sum (29) can also be expressed in terms of the orthogonal components described in Section 3 and defined by (16) as follows: 𝑋_𝑡=^𝐿_𝑖=1𝑤_𝑖𝐶^𝑛,𝑖_𝑋,𝑡, (31) where the weights 𝑤_𝑖 are computed as 𝑤_𝑖=𝑟^𝑖−1⋅1−𝑟^𝐿⋅𝜎1−𝑟var𝐶^{2𝑋𝑛,𝑖}_𝑋,𝑡, (32) where 𝑟 is defined by (20). Note that the only restriction of (32), similarly to (30), is that the estimated variance 𝐶^𝑛,𝑖_𝑋,𝑡 must be nonzero. Even though this synthesis does not depend either on the marginal distribution or the correlation structure of the input signal, it is preferable that this is self-similar (e.g., FGN) and that its Hurst index is close to that of 𝑋_𝐻1(𝑡), that is, 𝐻₁≈𝐻.

Pathological behavior can be produced in the output series for some critical conditions, for example, noticeable steps, which may produce a nonstationary signal, result from transforming an SRD, or uncorrelated input signal to an LRD output with 𝐻 close to 1. Impulses of very large magnitude (outliers) can be also be produced when, for some 𝑗, 𝑆₂(𝑗) is close to zero.

A similar methodology was developed by Deléchelle et al. [27], to synthesize fractional Gaussian noise by performing a weighted sum of the intrinsic mode functions (IMF) of a white noise process. The advantages of the proposed synthesis compared to that of Deléchelle et al. are that the components defined by (16) are exactly orthogonal, the relation between the weights for the reconstruction sum are mathematically well defined, and the Hurst index of the synthesized time series matches perfectly the wavelet estimator proposed by Abry et al. [4], that is, the estimated 𝐻 of the synthesized series is unbiased (𝐸(𝐻)=𝐻) and has zero variance (var(𝐻)=0). The disadvantages are that the components are sequences of squared signals (because of the expansion described in Section 3) and not sinusoids, and noticeable steps arise when synthesizing a time series with high Hurst index, for example, close to 1, from an input that is SRD or weakly correlated. A solution for this problem is to apply interpolation (as in EMD) instead of expanding the series in order to produce softer components (sinusoids or polynomial) instead of square type, with the consequence that the Hurst index is no longer exact, but approximated.

4. Estimation of Mean and Variance Self-Similar Time Series

4.1. Sample Mean

The sample mean of a self-similar process is unbiased: its expected value is the process mean, that is, 𝐸(𝑋)=𝜇_𝑋, where 𝑋=1/𝑁∑^𝑁_𝑡=1𝑋_𝑡, regardless of the presence of correlation between observations. However, its variance does not depend only on the sample size (𝑁) but also on the degree of self-similarity (𝐻) of the process as follows: var𝑋=𝜎²_𝑋𝑁^2𝐻−2, (33) which becomes var𝑋=𝜎²_𝑋𝑁, (34) (classical estimator) for 𝐻=0.5 (uncorrelated observations). Figure 4 shows the probability distribution function (PDF) of the sample mean of standardized FGN.

Figure 4 

Distribution of the sample mean of standardized fractional Gaussian noise processes with 𝐻={0.50,0.70,0.90} and 𝑁=32.

To derive (33), consider that 𝑋, estimated from a sample of size 𝑁, behaves exactly the same as the stationary aggregated process 𝑋_𝑘^(𝑁), defined by (1), and its variance is determined by the definition of second-order self-similarity (7). Expression (33) can be also derived (for 𝐻>0.5) from the autocorrelation coefficient 𝜌(𝑘)=0.5[(𝑘+1)^2𝐻−2𝑘^2𝐻+(𝑘−1)^2𝐻] for 𝑘≥1 (𝜌(0)=1) and var(𝑋)=(𝜎²_𝑥/𝑁²)∑^𝑁_𝑖,𝑗=1𝜌(𝑘).

Important implications of (33) about the uncertainty of the mean are (1) that it increases with the Hurst index, for example, var(𝑋) tends to 𝜎²_𝑋 as 𝐻 tends to 1, which makes the sample mean worth a single observation and (2) that it cannot be zero for any case when estimated from a finite-size sample.

4.2. Sample Variance

For a high number of observations, sample variance is usually calculated as 𝜎²_𝑋=1𝑁^𝑁_𝑡=1𝑋_𝑡−𝜇_𝑋². (35) It is known that estimator (35) is biased, so for small samples, it is more adequate to use 𝜎²_𝑋=1𝑁−1^𝑁_𝑡=1𝑋_𝑡−𝑋². (36)

Particularly, if the observations 𝑋_𝑡 are independent and come from a normal distribution, 𝜎²_𝑋 is distributed as 𝜎²_𝑋∼(𝜎²_𝑋/𝑁)𝜒²_𝑛−1, as stated by Cochran’s theorem [28, 29]. Formula (36) is the most used estimator of the sample variance [14] but, as Beran indicates in [30], it is needed to know which assumptions this estimator is based on in order to apply it correctly; otherwise, it may be the source of errors that in practice cannot be negligible for all cases.

A self-similar process is uncorrelated only and only if the Hurst index is 0.5. In this particular case, the classical estimator of sample variance, defined by (36), is unbiased.

The expected value of the sample variance defined by (36) is 𝐸𝜎²_𝑋1=𝐸𝑁−1^𝑁_𝑡=1𝑋_𝑡−𝜇², (37) which can be expressed as 𝐸𝜎²_𝑋=𝑁𝐸𝑋𝑁−1^2𝑡−𝐸𝜇², (38) then, 𝐸𝜎²_𝑋=𝑁𝑋𝑁−1var_𝑡+𝐸𝑋_𝑡²−var𝜇+𝐸𝜇²,𝐸𝜎²_𝑋=𝑁𝜎𝑁−1^2𝑋1−𝑁^2𝐻−2. (39)

Expression (39) proves that the classical estimator (36) is biased (i.e., 𝐸(𝜎²_𝑋)≠𝜎²_𝑋) for 𝐻≠0.5. It is straightforward that the unbiased variance estimator for self-similar processes is then 𝜎²_𝑋=1𝑁−𝑁^{𝑁2𝐻−1}_𝑡=1𝑋_𝑡−𝑋², (40) which obviously becomes (36) for 𝐻=0.5.

A plot of log₁₀(𝜎²_𝑋−𝜎²_𝑋) versus 𝑁 is shown in Figure 5. Note that, for a fixed sample size, as 𝐻 increases the estimation of the variance by means of the classical estimator, (36) becomes less significant.

Figure 5 

Logarithm of the variance bias for different sample sizes.

Figure 5 exemplifies that as the sample size is greater, the classical variance estimator has less bias; however, as the Hurst index of the process is greater, the bias is also greater (in magnitude). Note that as 𝐻 approaches to 1, the variance is considerably underestimated, which makes the classical variance estimator useless.

The variance of the estimated variance of a self-similar time series can be approximated by applying the formula proposed by Yunhua in [31] for 𝑘=0, that is, var𝜎²_𝑋=2𝑁^4𝐻−3+8𝐻−72𝐻²(2𝐻−1)²𝑁()+14𝐻−3𝑁. (41)

This approximation is close to the variance of  𝜎²_𝑋, with the disadvantage that it has a discontinuity in 𝐻=0.75. Further work can be developed in order to verify this approximation and to quantify its error.

Let us mention that, although the proposed estimator of the sample variance is unbiased, its performance relies on the estimation of the Hurst index. This dependence is very noticeable as 𝐻 approaches to 1, as the statistic 1−𝑁²𝐻₋₂ is especially sensitive to the variation of  𝐻 under that condition. The derivative 𝑑(1−𝑁²𝐻₋₂)/𝑑𝐻=−2𝑁²𝐻₋₂ln(𝑁) versus 𝐻 is shown in Figure 6. Note that as the estimated 𝐻 increases, the estimation of the mean of the sample variance is more variable.

Figure 6 

Plot −2𝑁2𝐻−2ln(𝑁) versus 𝐻.

An immediate implication of this is that processes with Hurst index close to 1 must be carefully treated, as slight deviations of the Hurst index estimation derive in a nonnegligible error in the estimation of the process variance.

4.3. Statistics of the Aggregated Process

The aggregated process 𝑋_𝑘^(𝑚) (defined by (2)) derived from an 𝐻-SS process is also 𝐻-SS (self-similar with the same Hurst index). It is also true for the case of 𝐻-SOSS processes, and this aggregated process is, by definition, identically distributed to the sample mean obtained from a set of 𝑁=𝑚 observations (𝑚 is the aggregation level, as in (2)), that is, 𝑋_𝑘^(𝑚)∼𝜇. Also, the aggregated sample variance obtained with the classical estimator (36) is biased, as it is well known [32].

The variance of the aggregated series of an 𝐻-SOSS process must be then estimated with the unbiased formula (40) adapted to the number of observations in the sample, that is, var𝑋_𝑘^(𝑚)=1𝑁_𝑖−𝑁_𝑖^{𝑁2𝐻−1𝑖}_𝑘=1𝑋_𝑘^(𝑚)−𝑋_𝑘^(𝑚)2, (42) where 𝑁_𝑖 is the size of the series after aggregation (i.e., 𝑁_𝑖=𝑁/𝑚). Note that the estimation of the sample mean from the aggregated sample is also unbiased, that is, 𝐸[𝑋_𝑘^(𝑚)]=𝐸(𝑋_𝑡), and it is more reliable than the mean estimated from a sample of the same size, that is, var[𝑋_𝑘^(𝑚)]=𝑚^2𝐻−2var(𝑋) if the samples are of equal size.

4.4. Statistics of the Orthogonal Components

Let {𝑋_𝑡;𝑡=1,…,𝑁} be a finite-length self-similar time series such that 𝑁=𝑛^𝐽 (𝐽<∞) and 𝑛≥2 (i.e., 𝑁 is a power of 𝑛); then a set of nonzero 𝐽 components (𝐶^𝑛,𝑗_𝑋,𝑡; 𝑗=1,…,𝐽) can be obtained as expressed by the analysis (16). As the components are pair-wise orthogonal, the variance of 𝑋_𝑡 (𝜎²_𝑋) is the sum of a finite number of variances: 𝜎²_𝑋=^𝐽_𝑗=1var𝐶^𝑛,𝑗_𝑋,𝑡. (43)

Expressions (43) and (19) imply that the variance of the 𝑖th component (𝐶^𝑛,𝑖_𝑋,𝑡) (computed using formula (35)), which has also finite length, is var𝐶^𝑛,𝑖_𝑋,𝑡=1−𝑟1−𝑟^𝐽𝑟^𝑖−1𝜎²_𝑋, (44) and estimating var(𝐶^𝑛,𝑗_𝑋,𝑡) as in (23), it follows that var𝐶^𝑛,𝑖_𝑋,𝑡=var𝐶^𝑛,𝑗_𝑋,𝑡1−𝑁^2𝐻−21−𝑟^𝐽. (45)

Replacing 𝑟^𝐽 by 𝑁^2𝐻−2, it yields var𝐶^𝑛,𝑖_𝑋,𝑡=var𝐶^𝑛,𝑗_𝑋,𝑡. (46)

Expression (46) implies that the estimation of the variance of components is unbiased (a desirable property) and, as a consequence, so it is the estimation of the statistic 𝑆₂(𝑗). Another implication of (46) is that the estimations of the Hurst index and the power parameter 𝑐_𝑓𝐶 from the Logscale Diagram are unbiased, as have previously been proven by the authors of [11].

4.5. Correlation of the Wavelet Coefficients

Let us describe the autocovariance of the Haar wavelet coefficients, that is, the case for 𝑛=2. The structure of this 𝑗th component (𝐶^2,𝑗_𝑋,𝑡), as a function of the elements of 𝑋_𝑡, is 𝐶^2,𝑗_𝑋,𝜏=𝑥^(2𝑗−1)_2𝜏−1−𝑥^(2𝑗−1)_2𝜏2. (47) Note that 𝐶^2,𝑗_𝑋,𝜏 is a downsampled version of 𝐶^2,𝑗_𝑋,𝑡, that is, only the first observation of each 2^𝑗 of 𝐶^2,𝑗_𝑋,𝑡 remains.

Then, two consecutive coefficients are correlated as 𝛾_𝐶𝐽(𝑘)=𝐸𝐶^2,𝑗_𝑋,𝜏,𝐶^2,𝑗_{𝑋,𝜏+𝑘}𝑥=𝐸^{(2𝑗−1)2𝜏−1}−𝑥^(2𝑗−1)_2𝜏2𝑥^(2𝑗−1)_{2()𝜏+𝑘−1}−𝑥^(2𝑗−1)_2()𝜏+𝑘2. (48)

Assuming that 𝑋_𝑡 represents a zero-mean 𝐻-SOSS process, and according to the definition (8), 𝐸(𝑥^(2𝑗−1)_2𝜏−1𝑥^(2𝑗−1)_{2(𝜏+𝑘)−1})=𝐸(𝑥^(2𝑗−1)_2𝜏𝑥^(2𝑗−1)_2(𝜏+𝑘))=𝛾_𝑋(2𝑘), 𝐸(𝑥^(2𝑗−1)_2𝜏−1𝑥^(2𝑗−1)_2(𝜏+𝑘))=𝛾_𝑋(2𝑘+1) and 𝐸(𝑥^(2𝑗−1)_2𝜏𝑥^(2𝑗−1)_{2(𝜏+𝑘)−1})=𝛾_𝑋(2𝑘−1). Then, (48) is calculated as 𝛾_𝐶𝑗(𝑘)=2^2𝐻−2𝜎²_𝑋42𝜌_𝑋(2𝑘)−𝜌_𝑋(2𝑘+1)−𝜌_𝑋(2𝑘−1), (49) and the correlation coefficient of the 𝑗th component is then 𝜌_𝐶𝑗(𝑘)=𝜌_𝐶(𝑘)=𝜌_𝑋(2𝑘+1)−2𝜌_𝑋(2𝑘)+𝜌_𝑋(2𝑘−1)2^2𝐻.−4 (50)

As (50) shows, the correlation structure is the same for all components, that is, 𝜌_𝐶𝑗(𝑘) is independent of 𝑗. Other implications of (50) are that ^∞_𝑘=0𝜌_𝐶𝑗(𝑘)=1+−𝜌_𝑋(1)−𝜌_𝑋(2𝑘+1)+2∑^2𝐾+1_𝑖=1(−1)^𝑖−1𝜌_𝑋(𝑖)2^2𝐻,−4 (51) which can be approximated as ^∞_𝑘=0𝜌_𝐶𝑗(𝑘)≈−0.052𝐻²−0.311𝐻+1.168<∞, (52) that is, the coefficients are weakly correlated and the sum of correlations is finite [33]; furthermore, none of the components (sequences of wavelet coefficient) can be a self-similar time series except for 𝐻=0.5, for example, components of a white noise process are also white noise processes. A plot of 𝜌_𝐶𝑗(𝑘) versus 𝑘 is shown in Figure 7 for 𝐻={0.1,0.3,0.5,0.7,0.9}. The sum ∑^∞_𝑘=0𝜌_𝐶𝑗(𝑘) versus 𝐻 is shown in Figure 8.

Figure 7 

Correlation coefficient of the Haar wavelet components.

Figure 8 

Sum of correlation coefficients of the Haar wavelet components.

As the maximum magnitude of 𝛾_𝐶1(1) is 1/8 (when 𝐻→0), it is said that these coefficients are quasiuncorrelated. Note that 𝛾_𝐶1(1)=0 for 𝐻=0.5 and 𝐻→1.

Although the assumption that the estimation of the 1st component variance (var(𝐶^𝑛,1_𝑋,𝑡)) is unbiased is nearly accurate, it may not hold for components of greater order. As the wavelet coefficients 𝑑_𝑋(𝑗,𝑘) are almost uncorrelated, the estimation of the component variance (i.e., ∑^{𝑁/𝑛𝑖}_𝑡=1(𝐶^𝑛,𝑖_𝑋,𝑡)²) is approximately ((𝑁−𝑛^𝑖)/𝑁)var(𝐶^𝑛,𝑖_𝑋,𝑡).

5. Variance Plot-Based Estimation of the Hurst Index

As described in Section 2.3, the Variance Plot is a straight line for self-similar time series but, as many authors have claimed, it underestimates the Hurst index when working with real world data. This is a consequence of the inadequate estimation of the aggregated variance, caused by the application of the classical formula (36) regardless of whether the original process presents any type of correlation. The solution would be then to apply the unbiased formula (42), but it leads to an ill-conditioned problem: the Hurst index needs to be estimated and known at the same time. Then, it is not that the Variance Plot is not adequate to estimate the Hurst index, rather the flaw of those implementations is that the aggregated variance is underestimated.

Nevertheless, it is actually possible to estimate the Hurst index analytically from the Variance Plot: the key is to choose the aggregation levels so that they form a geometric series, as explained in Section 5.1. Note that a numerical method can also be applied to estimate simultaneously the Variance Plot and the Hurst index, but the proposed solution is computationally simpler and more efficient.

5.1. Analytical Solution to the Ill-Conditioned Problem

Let {𝑚_𝑖;𝑖=1,…,𝑀} be the set of aggregation levels such that 𝑚_𝑖=𝑎𝑚_𝑖−1=𝑎^𝑖−1𝑚₁, 𝑎,𝑚₁∈ℕ and 𝑎>1, that is, the levels of aggregation follow a geometric series, for example, {𝑚_𝑖}={2,4,8,…,2^𝑀} or {𝑚_𝑖}={10,100,1000,…,10^𝑀}, and let 𝜎²_{𝑋𝑖)(𝑚}be the variance of the aggregated series 𝑋^(𝑚𝑖)_𝑘 estimated with formula (35). Obviously, 𝜎²_{𝑋𝑖)(𝑚} is biased, as 𝜎²_{𝑋𝑖)(𝑚}=𝜎²_{𝑋𝑖)(𝑚}[1−(𝑁/𝑚_𝑖)^2𝐻−2]=𝜎²_𝑋𝑚_𝑖^2𝐻−2(1−𝑁^2𝐻−2), but then the difference between 𝜎²_{𝑋𝑖)(𝑚} and 𝜎²_{𝑋)(𝑚𝑖+1} is calculated as follows: Δ_𝑖=𝜎²_{𝑋𝑖)(𝑚}−𝜎²_{𝑋)(𝑚𝑖+1}=𝜎²_𝑋𝑚₁^2𝐻−2𝑎^{(𝑖−1)(2𝐻−2)}1−𝑎^2𝐻−2,𝑖=1,…,𝑀−1, (53) and its logarithm log_𝑎Δ_𝑖=log_𝑎𝜎²_𝑋𝑚₁^2𝐻−21−𝑎^2𝐻−2+(𝑖−1)(2𝐻−2). (54)

Finally, the slope (𝑠) of the plot log_𝑎Δ_𝑖 versus 𝑖 is obtained (e.g., with a weighted least square regression) and 𝐻 is estimated as 𝐻=𝑠/̂2+1. It can be easily proven as 𝜎²_{𝑋𝑖)(𝑚}=𝑋^{(𝑛𝑗−1}_𝑡^𝐸) and substituting it in (16) and (19).

The slope is computed by the following weighted formula: ̂∑𝑠=^{𝑀−1𝑖=1}𝑖Δ_𝑖𝑊_𝑖−∑^𝑀−1_𝑖=1𝑖𝑊_𝑖⋅∑^𝑀−1_𝑖=1Δ_𝑖𝑊_𝑖∑^𝑀−1_𝑖=1𝑖²𝑊_𝑖−∑^𝑀−1_𝑖=1𝑖𝑊_𝑖², (55) where the weights are such that ∑^𝑀−1_𝑖=1(𝑊_𝑖)=1 and they are adequate so that 𝑠̂ has minimal variance, for example, 𝑊_𝑖=𝑊_𝑖−1/𝑚₁.

Note that the Variance Plot can be estimated without bias, but the Hurt index is not estimated from it. Furthermore, if the aggregation levels are taken as 𝑚_𝑖=2^𝑖, the estimator is exactly the same than the one that uses Haar wavelet. The authors of [34] developed an empirical study of estimation of the Hurst index from series with the presence of trends. They conclude that a method named differenced-variance (a variation of the variance-bases estimator) should not be used for estimating the Hurst index. The proposed solution is also a differenced-variance type method, but it can be used to estimate the Hurst index without bias and with optimal variance. Evidently, when working with real world traces, the Variance Plot may differ from the straight line, and an additional bias results from the logarithm as 𝐸[log(⋅)]≠log[𝐸(⋅)]. This bias can be subtracted analogously to (15), that is, log_𝑎var𝑋^(𝑚𝑖)_𝑡=log_𝑎var𝑋^(𝑚𝑖)_𝑡−𝑔_𝑖, (56) where 𝑔_𝑖 is the bias defined in (14).

6. Simulation and Measurements

6.1. Estimation of the Sample Mean

In order to verify the equations that describe the mean and variance of the sample mean, a set of zero-mean, unitary variance, and FGN time series of size 𝑁_𝑝=10⁶ observations are generated using an implementation of the generator proposed by Davies and Harte in [35], each for a different Hurst index for 𝐻={0.30,0.50,0.70,0.90}. Then, the mean is estimated from blocks of size 𝑁=100 and the empirical PDF is obtained from the estimations and compared to the classical (34) and proposed (33) estimators. Figure 9(a) shows that the variance of the estimated mean does not fit the classical model when SRD or LRD is present (Figures 9(a), 9(c), and 9(d)), but only for the uncorrelated case (Figure 9(b)). Only proposed estimator (33) represents adequately this phenomenon for the four cases.

(a)

(b)

(c)

(d)

Figure 9 

Estimation of the sample mean with 𝑛=100 for four cases: (a) 𝐻=0.30, (b) 𝐻=0.50, (c) 𝐻=0.70, and (d) 𝐻=0.90.

6.2. Estimation of the Sample Variance

The followed procedure to verify the proposed estimator of the sample variance (𝜎²_𝑋) consists of the generation of a set of 100 FGN samples of size 1024 for each value of 𝐻={0.05,0.10,0.15,…,0.95} and the estimation of the variance using the classical formula (36) and the proposed estimator (40). The respective mean of  both estimations for each set was obtained. For the estimation of the Hurst index the wavelet estimator of Veitch and Abry [11] is used.

Figure 10 shows that the classical formula underestimates the variance noticeably for higher values of 𝐻. For 𝐻>0.95 the estimated variance is less than half of the process variance. The proposed formula (40) does not underestimate the variance, but for high values of 𝐻 the estimated variance is significantly different from one realization to another. This variation results from the estimation of the Hurst index, as the statistic 1−𝑁²𝐻₋₂ is very sensitive to the variations of 𝐻, which depends, in turn, on the efficiency of the sample generator. This is an indicator that the generator proposed by Davies and Harte may be less accurate as 𝐻 is closer to 1.

Figure 10 

Mean of the estimated variance for 𝑛=1024 and 𝐻={0.05,0.10,0.15,…,0.95}.

Figure 11 shows the variance of the estimated variance obtained with the proposed formula (40) compared to the approximation proposed by Yunhua in [31], as expected, the variance of the estimation is close to zero (lower than 0.001) for 𝐻<0.5, but as the Hurst index increases, it becomes noticeable when 𝐻>0.75. This verifies the observation of [31], which says that beyond 𝐻=0.75 the precision of the autocorrelation is about one order lower than when 𝐻<0.75.

Figure 11 

Variance of the estimated variance for 𝑛=1024 and 𝐻={0.05,0.10,0.15,…,0.95}.

6.3. Synthesis of 𝐻-SOSS Time Series

To exemplify the proposed wavelet-based synthesis (described in Section 3.2), four time series with respective Hurst index 0.3, 0.5, 0.7, and 0.9 were synthesized from an FGN sample of size 1024. The Logscale Diagram of the four new time series were obtained and compared to that of the original sample. The plot 𝑋_𝑡 versus 𝑡 for each one; the four synthesized series is shown in Figure 12. One can visually check the presence of positive correlation in Figures 12(c) and 12(d).

(a)

(b)

(c)

(d)

Figure 12 

Plot versus time of the four synthesized time series. (a) 𝐻=0.30, (b) 𝐻=0.50, (c) 𝐻=0.70, and (d) 𝐻=0.90.

The Logscale Diagram of these artificial series is shown in Figure 13. Note that the original Logscale Diagram of the source sample is not a straight line, but it so is for the synthesized series. Also, the estimated Hurst index of this series is 𝐻=0.56 and its estimated Logscale Diagram is not a straight line, but the estimated Hurst index of the four generated series is exactly the desired, for example, 𝐻=0.30 for the series shown in Figure 12(a) and the same for the others, and their respective Logscale Diagram is a straight line.

Figure 13 

Logscale Diagram of synthesized time series.

6.4. Voice over IP (VoIP) Measurements

The jitter behavior of Voice over Internet Protocol (VoIP) traffic by means of networks measurements is analyzed. As result of this analysis, detailed characterization and accurate modeling of this Quality of Service (QoS) parameter is provided. Previous studies have revealed that VoIP jitter can be modeled by self-similar processes with short-range dependence (SRD) or long-range dependence (LRD) [36]. The discovery of LRD (a kind of asymptotic fractal scaling) and weak self-similarity in the VoIP jitter data traces was followed by a further work that shows the evidence for multifractal behavior. The discovery of evidence for multifractal behavior is a richer form of scaling behavior associated with nonuniform local variability, which could lead to a complete and robust model of IP network traffic over all time scales.

Motivated by such concerns, the evidence for multifractal behavior of VoIP jitter data traces is reviewed. In order to accomplish this analysis, the time series of VoIP jitter into a set of time series or components 𝐶^2,𝑗_𝑋,𝜏 is decomposed as defined by (16). The behavior of these components is used to determine the kind of asymptotic fractal scaling. If the variance of the components of a time series is modeled by a straight line, the time series exhibits monofractal behavior, and a linear regression can be applied in order to estimate the Hurst parameter. On the other hand, if the variance of the components cannot be adequately modeled with a linear model, then the scaling behavior should be described with more than one scaling parameter, that is, the time series exhibits multifractal behavior [23]. In Figures 14 and 15, we show the components behavior of the collected VoIP jitter data traces.

Figure 14 

Components behavior of VoIP jitter data traces: monofractal behavior.

Figure 15 

Components behavior of VoIP jitter data traces: multifractal behavior.

Figure 14 shows the components behaviors of a VoIP jitter data trace that belongs to the data sets with SRD. It is observed that the variance of the components of this time series is modeled by a straight line; therefore, the time series exhibits monofractal behavior.

Figure 15 shows the components behaviors of a VoIP jitter data trace that belongs to the data sets with LRD. It is observed that the variance of the components of this time series cannot be adequately modeled with a linear model, and the scaling behavior should be described with multiple scaling parameters (biscaling); therefore, this time series exhibits multifractal behavior.

These results show that VoIP jitter with SRD or LRD exhibit monofractal or multifractal behavior, respectively. This phenomenon explains the behavior of the data traces with SRD and high degree of self-similarity (scale invariance), because the self-similarity is defined for a single scale parameter. On the other hand, the data traces with LRD exhibit weak self-similarity because they have associated nonuniform local variability (multifractal behavior).

The implication of this behavior for VoIP and other interactive multimedia services is that receiver dejitter buffer may not be large enough to mask the jitter with LRD and multifractal characteristics.

7. Conclusion

An orthogonal decomposition, that constitutes a powerful statistic tool to study discrete time series, is presented. The resulting components (𝐶^𝑛,𝑗_𝑋,𝑡) have zero-mean and three desirably properties: (1) they synthesize the source signal without loss, (2) they are pair wise orthogonal, and (3) their variances form a geometric series whose rate is related to the Hurst index.

This decomposition is firstly compared to Haar wavelet based. For a particular case these two coincide, but the proposed one is more general, as other levels of aggregation may occur. In this case, the Haar wavelet is not sufficient and a special class of basis functions are defined, and the wavelet coefficients are obtained by the inner products between the signal under study and a scaled, displaced, and phase shifted versions of the basis function, in contradistinction of other wavelet decomposition, that only apply scaling and displacement. For a fixed scaling and displacement, the last phase shift can be discarded (i.e., the components are downsampled) as it does not provide additional information.

The proposed decomposition can be used to estimate the Hurst index, as the plot 𝑗+log_𝑛[var(𝐶^𝑛,𝑗_𝑋,𝑡)] versus 𝑗 is equivalent to the Logscale Diagram proposed by Veitch and Abry in [11]. It is a straight line for 𝐻-SOSS time series, and the slope is related to the Hurst index so that 𝑠=2𝐻−1.

The components can be also used to synthesize 𝐻-SOSS time series by means of the weighted sum defined by (31) and (32), regardless of the distribution of the source signal and whether it is or not self-similar. The synthesis is exact, that is, the Hurst index of the synthesized series is exactly the desired, which is an advantage over the proposed synthesis, that uses IMF from an EMD decomposition, described in [27].

A study of the estimated mean and variance of self-similar time series is also presented. Both statistics, mean and variance, can be estimated without bias, by applying the classic sample mean (𝑋=1/𝑁∑^𝑁_𝑡=1𝑋_𝑡) and the proposed variance estimator (40). The variance of these estimators depends on the Hurst index of the process. In the case of the sample mean, its variance is an increasing function of 𝐻, as expressed by (33), such that var(𝑋)∈(𝜎²_𝑋𝑁⁻²,𝜎²_𝑋). Note that the sample mean becomes less significant as 𝐻 approaches to 1. For the variance of  𝜎²_𝑋, it can be observed that it has the best performance when 𝐻=0.5 (the variance of 𝜎²_𝑋 is minimal). For 𝐻<0.5, the variance increases as 𝐻 is lower, but it is still acceptable (e.g., var(𝜎²_𝑋)<0.001𝜎²_𝑋 when 𝐻→0. But when 𝐻 gets close to 1, the uncertainty of the sample variance increases rapidly, making the estimation of 𝜎²_𝑋 less significant.

It is demonstrated also (clarifying a popular misunderstanding) that the Variance Plot can be used to estimate efficiently the Hurst index. The claims of many researchers about the inefficiency of the Variance Plot are the result of an ill-conditioned problem: to estimate the Hurst index the variance of the aggregated series is needed and vice versa, leading to a vicious circle. It is shown how this problem is avoided by estimating the Variance Plot using a set of aggregation levels that follow a geometric series and calculating the differences between the variances of the corresponding aggregated series; then, a weighted linear regression is applied to estimate the slope (𝑠̂), and the estimated Hurst index is 𝐻=𝑠/̂2+1. This result gives more significance to the results published by Abry et al. in [26].

Acknowledgment

This work was supported by PROMEP Grant UQROO-PTC-110.