New Approaches to Multiscale Signal Processing

Material Information

New Approaches to Multiscale Signal Processing
HU, JING ( Author, Primary )
Copyright Date:


Subjects / Keywords:
Brownian motion ( jstor )
Chaos ( jstor )
Chaos theory ( jstor )
Eggshells ( jstor )
Fractals ( jstor )
Radar ( jstor )
Seas ( jstor )
Signals ( jstor )
Time series ( jstor )
Trajectories ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Jing Hu. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
Resource Identifier:
659872455 ( OCLC )


This item has the following downloads:

hu_j ( .pdf )





















































































































Full Text






S2007 Jing Hu

To my family, teachers and friends


First of all, I would like to thank my advisor, Professor .Jianho Gao, for his invaluable

guidance and support. I was lucky to work with someone with so many original ideas

and such a sharp mind. Without his help, my Ph.D. degree and this dissertation would

have been impossible. During my years at the University of Florida, I was moved by his

magnificent personality. A famous Chinese adage ; .va, "One d .< as my advisor, entire life

as my father." Professor Gao highly qualifies.

I also want to thank my coninittee nienters (Professors Su-Shing C'I. in~ Yuguang

Fang, .Jose C. Principe, and K~eith D. White), for their interest in my work and valuable

-a-----Hul-My thanks of course also go to all the faculty, staff and my fellow students in

the Department of Electrical and Computer Engineering. I extend particular thanks to my

office mates .Jaentin Lee, Yi Z1n in_ and IUngsik K~in, for their helpful discussions.

I express my appreciation to all of my friends who offered encouragement. I specially

thank Qian Zhan, Xingyan Fan, Lingyu Gu and his wife (.Jing Zhang), Yannmin Xu, Yan

Liu, and Yana Liang. They gave me a lot of help during my early years in the University

of Florida. They taught me to cook and made my life more colorful.

Finally, I owe a great debt of thanks to my parents, grandparents, and my boyfriend,

.Jing Ai, for their pillar-like support and their unwavering relief and confidence in me

without this I do not think I would have made it this far!



ACK(NOWLEDGMENTS ......... . . 4

LIST OF FIGURES ......... .. . 7

ABSTRACT ......... ..... . 9


1 INTRODUCTION ......... . 10

1.1 Overview of Dynamical Systems . ..... .. 11
1.2 Examples of Multiscale Phenomena . .... .. 15
1.3 Brief Introduction to C'!s I... ......... .. 20
1.4 Brief Introduction to Fractal Geometry ... .. . . 21
1.5 Importance of Connecting ('!, I. .s and Random Fr-actal Theories .. .. 2:3
1.6 Importance of the Concept of Scale . ..... .. 24
1.7 Structure of this Dissertation ........ .. .. 24

CONDITIONS (PSIC) ......... ... 27

2.1 Dynamical Test for C'!s I... .... .. . .. .. 27
2.2 General Computational Framework for Power-law Sensitivity to Initial
Conditions (PSIC) .. . .. .. .. .. 3:3
2.3 C'll I l .terizing Edge of C'I!s I. by Power-law Sensitivity of Initial Conditions
(PSIC) .... ........ ............ :36


:3.1 C'll I l .terizing 1/ fP Processes by Power-law Sensitivity to Initial Conditions
(PSIC) .... ........ .. .. ...... ..... 41
:3.1.1 C'll I l .terizing Fractional Brownian Motion (fBm) Processes by
Power-law Sensitivity to Initial Conditions (PSIC) .. .. .. .. 4:3
:3.1.2 C'll I l .terizing ON/OFF Processes with Pareto distributed ON and
OFF Periods by Power-law Sensitivity to Initial Conditions (PSIC) 50
:3.2 C'll I l .terizing Levy Processes by Power-law Sensitivity to Initial Conditions
(PSIC) .... ........ ............ 54

INITIAL CONDITIONS (PSIC) ........ ... .. 60

4.1 Study of Transport Dynamics by Network Simulation .. .. .. .. 61
4.2 Complicated Dynamics of Internet Transport Protocols .. .. .. .. 69


5.1 Sea Clutter Data ........ ... .. 78
5.2 Non C'!s I..tic Behavior of Sea Clutter .... .. .. 79
5.3 Target Detection within Sea Clutter by Separating Scales .. .. .. .. 81
5.4 Target Detection within Sea Clutter by Power-law Sensitivity to Initial
Conditions (PSIC) ........ . .. 86

(SDLE) ........... ......... .. 89

6.1 Basic Theory
6.2 Classification of Complex Motions.
6.2.1 C'! I ..-~ Noisy (I I ..-~ and Noise-induced ChI ... .-
6.2.2 Processes Defined by Power-law Sensitivity to Initial
6.2.3 Complex Motions with Multiple Scaling Behaviors.
6.3 C'!I. .) ..terizing Hidden Frequencies


7 CONCLUSIONS ......... ... .. 104

REFERENCES .. .......... ........... 105



Figure page

1-1 Example of a trajectory in the phase space ..... .. .. 14

1-2 Giant ocean wave (tsunami). ......... .. 17

1-3 Example of sea clutter data. ......... .. 18

1-4 Example of Heart rate variability data for a normal subject. .. .. .. 19

2-1 Time-dependent exponent A(k) vs. evolution time k curves for Lorenz system. .32

2-2 Time-dependent exponent A(k) vs. evolution time k curves for the Mackey-Glass
system. ......... ..... 33

2-3 Time-dependent exponent A(t) vs. t curves for time series generated from logistic
map. ......... ............ .... 37

2-4 Time-dependent exponent A(t) vs. t curves for time series generated from Henon
map. ............ ............... 40

3-1 White noise and Brownian motion . ...... .. 46

3-2 Several fBm processes with different H. . ... .. 48

3-3 Time-dependent exponent A(k) vs. In k curves for fractional Brownian motion. 51

3-4 Example of ON/OFF processes. .... ... .. 51

3-5 Time-dependent exponent A(k) vs. In k curves for ON/OFF processes with Pareto
distributed ON and OFF periods. ... ... .. 53

3-6 Examples for Brownian motions and Levy motions. .. .. .. 57

3-7 Time-dependent exponent A(k) vs. In k curves for Levy processes. .. .. .. 58

4-1 Example of quasi-periodic congestion window size W(i) time series. .. .. .. 63

4-2 Examples of chaotic T(i) time series. . ..... .. .. 65

4-3 Complementary cumulative distribution functions (CCDFs) for the two chaotic
T(i) time series of Figs. 4-2(a,c). . ... ... .. 66

4-4 The time series T(i) extracted from a W(i) time series collected using net100
instruments over ORNL-LSU connection. ..... .. .. 68

4-5 Time series for the congestion window size cwnd for ORNL-GaTech connection. 72

4-6 Time-dependent exponent A(k) vs. k curves for cwnd data corresponding to
Figf.4-5. ......... ..... 73

4-7 Time-dependent exponent A(k) vs. In k curves for cwnd data corresponding to
Figf.4-5. ......... ..... 75

5-1 Collection of sea clutter data ......... ... 79

5-2 Typical sea clutter amplitude data. ........ .. .. 79

5-3 Two short segments of the amplitude sea clutter data severely affected by clipping. 80

5-4 Examples of the time-dependent exponent A(k) vs. k curves for the sea clutter
data ........ ... . .... 81

5-5 Variations of the Lyapunov exponent A estimated by conventional methods vs.
the 14 range bins for the 10 HH measurements. .... .. .. 84

5-6 Variations of the scale-dependent exponent corresponding to large scale vs. the
range bins for the 10 HH measurements. ..... .. .. 85

5-7 Examples of the A(k) vs. In k curves for the sea clutter data. .. .. .. .. 87

5-8 Variation of the P parameter with the 14 range bins. .. .. ... .. 87

5-9 Frequencies of the bins without targets and the bins with primary targets for
the HH datasets. .. ... .. .. 88

6-1 Scale-dependent Lyapunov exponent A(e) curves for the Lorenz system and logistic
map ........ ... . .... 93

6-2 Scale-dependent Lyapunov exponent A(e) curve for the Mackey-Glass system. .. 94

6-3 The functions F(x) and G(x). ......... ... .. 97

6-4 Scale-dependent Lyapunov exponent A(e) for the model described by Eq. 6-9. .97

6-5 Power-spectral density (PSD) for time series generated from Lorenz system. 100

6-6 Lorenz attractor. .. ... .. .. 101

6-7 Hidden frequency phenomenon of laser intensity data. ... .. .. .. 102

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



Jing Hu

May 2007

C'I I!r: Jianho Gao
Major: Electrical and Computer Engineering

Complex systems usually comprise multiple subsystems with highly nonlinear

deterministic and stochastic characteristics, and are regulated hierarchically. These

systems generate signals with complex characteristics such as sensitive dependence on

small disturbances, long nienory, extreme variations, and nonstationarity. Chaos theory

and random fractal theory are two of the most important theories developed for analyzing

complex signals. However, they have entirely different foundations, one being deterministic

and the other being random. To synergistically use these two theories, we developed a new

theoretical framework, power-law sensitivity to initial conditions (PSIC), to encompass

both chaos and random fractal theories as special cases. To show the power of this

framework, we applied it to study two challenging and important problems: Internet

dynamics and sea clutter radar returns. We showed that PSIC can readily characterize

the complicated Internet dynamics due to interplay of nonlinear AINID (additive increase

and niultiplicative decrease) operation of TCP and stochastic user behavior, and robustly

detect low observable targets from sea clutter radar returns with high accuracy.

We also developed a new niultiscale complexity measure, scale-dependent Lyapunov

exponent (SDLE), that can he computed from short noisy data. The measure readily

classified various types of complex motions, and simultaneously characterized behaviors of

complex signals on a wide range of scales, including complex irregular behaviors on small

scales, and orderly behaviors, such as oscillatory motions, on large scales.


A dynamical system is one that evolves in time. Dynamical systems can be stochastic

(in which case they evolve according to some random process such as the toss of a coin)

or deterministic (in which case the future is uniquely determined by the past according

to some rule or mathematical formula). Whether the system in question settles down to

equilibrium, keeps repeating in cycles, or does something more complicated, dynamics are

what we use to analyze the behavior of systems.

Complex dynamical systems usually comprise multiple subsystems with highly

nonlinear deterministic and stochastic characteristics, and are regulated hierarchically.

These systems generate signals with complex characteristics such as sensitive dependence

on small disturbances, long memory, extreme variations, and nonstationarity [1]. A

stock market, for example, is strongly influenced by multi- l-1~ we decisions made by

market makers, as well as influenced by interactions of heterogeneous traders, including

intradwl~i traders, short-period traders, and long-period traders, and thus gives rise to

highly irregular stock prices. The Internet, as another example, has been designed in a

fundamentally decentralized fashion and consists of a complex web of servers and routers

that cannot be effectively controlled or analyzed by traditional tools of queuing theory

or control theory, and give rise to highly bursty and multiscale traffic with extremely

high variance, as well as complex dynamics with both deterministic and stochastic

components. Similarly, biological systems, being heterogeneous, massively distributed, and

highly complicated, often generate nonstationary and multiscale signals. With the rapid

accumulation of complex data in life sciences, systems biology, nano-sciences, information

systems, and physical sciences, it has become increasingly important to be able to analyze

multiscale and nonstationary data.

Multiscale signals behave differently depending upon which scale the data are looked

at. In order to fully characterize such complex signals, the concept of scale has to be

explicitly considered, and existing theories need to be used synergistically, instead of

individually, since on different scale ranges different theories may be most pertinent.

This dissertation aints to build an effective arsenal for niultiscale signal processing by

synergistically integrating approaches based on chaos theory and random fractal theory,

and going beyond, to complement conventional approaches such as spectral analysis and

machine learning techniques. To make such an integration possible, two important efforts

have been made:

* Power-law sensitivity to initial conditions (PSIC): We developed a new
theoretical framework for signal processing, to create a coninon foundation for chaos
theory and random fractal theory, so that they can he better integrated.

* Scale-dependent Lyapunov exponent (SDLE): We developed an excellent
niultiscale measure, which is a variant of the finite size Lyapunov exponent (FSLE).
We proposed a highly efficient algorithm for calculating it, and showed that it can readily
classify different types of motions, including truly low-dintensional chaos, noisy chaos,
noise-induced chaos, random 1/ fP and cc-stable Levy processes, and complex motions with
chaotic behavior on small scales but diffusive behavior on large scales. The measure can
aptly characterize complex behaviors of real world niultiscale signals on a wide range of

In the rest of this chapter, we shall first present an overview of dynamics systems,

especially we shall point out why nonlinear systems are much harder to analyze than

linear ones. We then give a few examples of nmultiscale phenomena to show the difficulties

and excitement of niultiscale signal processing. After that, we shall briefly introduce some

background knowledge about chaos and fractal geometry. Then we discuss the importance

of connecting chaos theory and random fractal theory for characterizing the behaviors

of niultiscale signals on a wide range of scales. We also emphasize the importance of

explicitly incorporating the concept of scale in devising measures for characterizing

niultiscale signals. Finally, we outline the scope of the present study.

1.1 Overview of Dynamical Systems

There are two main types of dynamical systems: differential equations and iterated

maps (difference equations). Differential equations describe the evolution of systems

in continuous time, whereas iterated maps arise in problems where time is discrete. To

show how systems evolve, we take differential equations as examples. When confining our

attention to differential equations, the main distinction is between ordinary and partial

differential equations. For instance, the equation for a damped harmonic oscillator

m +1 b + kx = 0 (1-1)

is an ordinary differential equation, because it involves only ordinary derivatives dx/dt and

d2 ld2. That is, there is only one independent variable, the time t. In contrast, the heat

Bu 82U
dt 8x2

is a partial differential equation: it has both time t and sapce x as independent variables.

Our concern in this proposal is with purely temporal behavior, and so we deal with

ordinary differential equations almost exclusively.

A very general framework for ordinary differential equations is provided by the system


Here the overdots denote differentiation with respect to t. Thus xi = dxi/dt. The

variables xl, x, might represent concentrations of chemicals in a reactor, populations

of different species in an ecosystem, or the positions and velocities of the planets in the

solar system. The functions fl, f, are determined by the problem at hand.

For example, the damped oscillator (Eq. 1-1) can be rewritten in the form of

(Eq. 1-2), thanks to the following trick: we introduce new variables xl = x and x2 x

Then xl = Z2. Hence the equivalent system (Eq. 1-2) is

xi = x2
b k
xa2=--x 2- 1--i
m m

This system is said to be linear, because all the xi on the right-hand side appear to

the first power only. Otherwise the system would be nonlinear. Typical nonlinear terms

are products, powers, and functions of the Xi, such aS xlx2, 1 3), Or COSZ2.

For example, the swinging of a pendulum is governed by the equation

x + -sinx = 0,

where x is the angle of the pendulum from vertical, g is the acceleration due to gravity,

and L is the length of the pendulum. The equivalent system is nonlinear:

x z = x2

Z2 = L~smnxl.

Nonlinearity makes the pendulum equation very difficult to solve analytically. The

usual way around this is to fudge, by invoking the small angle approximation sinx e x for

x << 1. This converts the problem to a linear one, which can then be solved easily. But

by restricting to small x, we are throwing out some of the physics, like motions where the

pendulum whirls over the top. Is it really necessary to make such drastic approximation?

It turns out that the pendulum equation can be solved analytically, in terms of elliptic

functions. But there ought to be an easier way. After all, the motion of the pendulum is

simple: at low energy, it swings back and forth, and at high energy it whirls over the top.

There should be some way of extracting this information from the system directly -using

geometric methods.

Here is the rough idea. Suppose we happen to know a solution to the pendulum

system, for a particular initial condition. This solution would be a pair of functions xl(t)

and x2(t), representing the position and velocity of the pendulum. If we construct an

I (t), X 2(t))

(x (0), X2(0))

Figure 1-1. Example of a trajectory in the phase space

abstract space with coordinates (xl,,Z2), then the solution (xl(t),2 xa1)) COTTOSponds to

a point moving along a curve in this space, as shown Fig. 1-1. This curve is called a

trajectory, and the space is called the phase space for the system. The phase space is

completely filled with trajectories, since each point can serve as an initial condition.

Our goal is to run this construction in reverse: given the system, we want to draw

the trajectories, and thereby extract information about the solutions. In many cases,

geometric reasoning will allow us to draw the trajectories without actually solving the


The phase space for the general system (Eq. 1-2) is the space with coordinates

xl, x,. Because this space is n-dimensional, we will refer to Eq. 1-2 as an n-dimensional

system or an nth-order system. Thus a represents the dimension of the phase space.

As we have mentioned earlier, most nonlinear systems are impossible to solve

analytically. Why are nonlinear systems so much harder to analyze than linear ones?

The essential difference is that linear systems can be broken down into parts. Then each

part can be solved separately and finally recombined to get the answer. This idea allows a

fantastic simplification of complex problems, and underlies such methods as normal modes,

Laplace transforms, superposition arguments, and Fourier analysis. In this sense, a linear

system is precisely equal to the sum of its parts.

But many things in nature do not act this way. Whenever parts of a system interfere,

or cooperate, or compete, there are nonlinear interactions going on. Most of every w- life

is nonlinear, and the principle of superposition fails spectacularly. If you listen to your two

favorite songs at the same time, you will not get double the pleasure! Within the realm of

physics, nonlinearity is vital to the operation of a laser, the formation of turbulence in a

fluid, and the superconductivity of Josephson junctions.

1.2 Examples of Multiscale Phenomena

Multiscale phenomena are ubiquitous in nature and engfineeringf. Some of them are,

unfortunately, catastrophic. One example is the Tsunami that occurred in south Asia

around C!~!s-I it -~ 2004. Another example is the gigantic power outage occurred in North

America on August 14, 2003, which affected more than 4000 megawatts, was more than

300 times greater than mathematical models would have predicted, and cost between $4

billion and 71 billion, according to the U.S. Department of energy. Both events involved

scales that, on one hand, were so huge that our human being could not easily fathom,

and on the other hand, involved the very scales that were most relevant to individual life.

Below, we consider three more examples.

Multiscale phenomena in computer networks: Large-scale communications

networks, especially the Internet, are among the most complicated systems that man

has ever made, with multiscales in terms of a number of aspects. Intuitively speaking,

these multiscales come from the hierarchical design of protocol stack, the hierarchical

topological architecture, and the multipurpose and heterogeneous nature of the Internet.

More precisely, there are multiscales in

*Time, manifested by the prevailing fract al, multifract al, and long-range- dependent
properties in traffic,

* Space, essentially due to topology and geography and again manifested by scale-free

* State (e.g., queues, windows),

* Size (e.g., number of nodes, number of users).

Also, it has been observed that the failure of a single router may trigger routing

instability, which may be severe enough to instigate a route flap storm. Furthermore,

packets may be delivered out of order or even get dropped, and packet reordering is not

a pathological network behavior. As the next generation Internet applications such as

remote instrument control and computational steering are being developed, another facet

of complex niultiscale behavior is beginning to surface in terms of transport dynanxics.

The networking requirements for these next generation applications belong to (at least)

two broad classes involving vastly disparate time scales:

* High handwidths, typically multiples of 10Ghps, to support bulk data transfers,

* Stable handwidths, typically at much lower handwidths such as 10 to 100 1\bps, to
support interactive, steering and control operations.

Multiscale phenomena in sea clutter: Sea clutter is the radar backscatter

front a patch of ocean surface. The complexity of the signals contes front two sources,

the rough sea surface, sometimes oscillatory, sometimes turbulent, and the nmultipath

propagation of the radar backscatter. This can he appreciated by imagining how radar

pulses massively reflecting front the wavetip of Fig. 1-2. To be quantitative, in Fig. 1-3,

two 0.1 s duration sea clutter signals, sampled with a frequency of 1 K(Hz, are plotted

in Fig. 1-:3(a,h), a 2 s duration signal is plotted in Fig. 1-:3(c), and an even longer signal

(about 1:30 s) is plotted in Fig. 1-:3(d). It is clear that the signal is not purely random,

since the waveform can he fairly smooth on short time scales (Fig. 1-:3(a)). However, the

signal is highly nonstationary, since the frequency of the signal (Fig. 1-:3(a,b)) as well as

the randoniness of the signal (Fig. 1-:3(c,d)) change over time drastically. One thus can

perceive that naive Fourier analysis or deterministic chaotic analysis of sea clutter may

Figure 1-2. Giant ocean wave (tsunami). Suppose our field of observation includes the
wavetip of length scale of a few meters, it is then clear that the complexity of
sea clutter is mainly due to massive reflection of radar pulses from a wavy and
even turbulent ocean surface.

not e very,,,,,, useul From Fig 1-3(e,\ where Xt(m) is the non-overlapping running mean of

X over block size m, and X is the sea clutter amplitude data, it can be further concluded

that neither autoregressive (AR) models (As an example, we give the definition for the

first-order AR model, which is given by x, I = ax, + rl,, where the constant coefficient a

satisfies 0 / |a| < 1 and rlo is a white noise with zero mean.) nor textbook fractal models

can describe the data. This is because AR modeling requires exponentially decaying
autcorrelation (which,1 amounts, to Var(Xe(m)) ~ m-l, or a Hurst parameter of 1/2. See

later chapters for the definition of the Hurst parameter), while fractal modeling requires

the variation between Vanr(Xe"m)) and m to follow a power-law~. However, neither behavior

is observed in Fig. 1-3(e). Indeed, albeit extensive work has been done on sea clutter,

the nature of sea clutter is still poorly understood. As a result, the important problem of

80 80.02 80.04 80.06
Time (sec)

80.08 8C

Time (sec)

log2 m

U 2U 4U bU0 tU 1UU 12U
Time (sec)

Figure 1-3.

Example of sea clutter data. (a,b) two 0.1 s duration signal; (c) a 2 s duration
sea clutter signal; (d) the entire sea clutter signal (of about 130 s); and (e)
log2 m2Var(X(m))] vs. log2 m, Where
Xim) {Xe(m) : X(m) (Xtm-m+l + + Xtm)/m, t = 1, 2, } is the
non-overlapping running mean of X = {Xt : t = 1, 2, } over block size m,
and X is the sea clutter amplitude data. To better see the variation of
Vawr(X~) with mVa(e) is multiplied by m2. When the autocorrelation
of the data decays exponentially fast (such as modeled by an autoregressive
(AR) process), V'ar(Xe(m")) ~ m-1 Here Var(Xe(m)) decaty-s much faster. A
fractal process would have m2"Var(Xe"m)) ~ m ". However, th~is is not thle case.
Therefore, neither AR modeling nor ideal textbook fractal theory can be
readily applied here.

target detection within sea clutter remains a tremendous challenge. We shall return to sea

clutter later.

Multiscale and nonstationary phenomena in heart rate variability (HRV):

HRV is an important dynamical variable of the cardiovascular function. Its most salient

feature is the spontaneous fluctuation, even when the environmental parameters are



0 1 2 3 4 5 6 7
Index n4
x 104
100 110
(b) (c)


2.66 2.67 2.68 2.69 2.7 2.71 3.4 .9 396 .7 398 .9
Index n x 104 Index n x 104
(d) (e)
104 105

S102 0

10-2 10- 10-2 10-
f f

Figure 1-4. Example of Heart rate variability data for a normal subject. (a) the entire
signal, (hsc) the segments of signals indicated as A and B in (a); (d,e) power
spectral density for the signals shown in (b,c).

maintained constant and no perturbing influences can he identified. It has been observed

that HRV is related to various cardiovascular disorders. Therefore, analysis of HRV is

very important in medicine. However, this task is very difficult, since HRV data are highly

complicated. An example is shown in Fig. 1-4, for a normal young subject. Evidently,

the signal is highly nonstationary and niultiscaled, appearing oscillatory for some period

of time (Figs. 1-4(h,d)), and then varying as a power-law for another period of time

(Figs. 1-4(c,e)). The latter is an example of the so-called 1/ f processes, which will be

discussed in depth in later chapters. While the niultiscale nature of such signals cannot

he fully characterized by existing methods, the nonstationarity of the data is even more

troublesome, since it requires the data to be properly segmented before further analysis

hy methods from spectral analysis, chaos theory, or random fractal theory. However,

automated segmentation of complex biological signals to remove undesired components is

itself a significant open problem, since it is closely related to, for example, the challenging

task of accurately detecting transitions from normal to abnormal states in physiological


1.3 Brief Introduction to Chaos

Imagine we are observing an aperiodic, highly irregular time series. Can such a

signal arise from a deterministic system that can he characterized by only very few state

variables instead of a random system with infinite numbers of degrees of freedom? A

chaotic system is capable of just that! This discovery has such far reaching implications

in science and engineering that sometimes chaos theory is considered as one of the three

most revolutionary scientific theories of the twentieth-century, along with relativity and

quantum mechanics.

At the center of chaos theory is the concept of sensitive dependence on initial

conditions: a very minor disturbance in initial conditions leads to entirely different

outcomes. An often used metaphor illustrating this point is that a sunny weather in New

York could be replaced by a rainy one some time in the future after a butterfly flaps

its wings in Boston. Such a feature contrasts sharply with the traditional view, largely

based on our experience with linear systems, that small disturbances (or causes) can only

generate proportional effects, and that in order for the degree of randomness to increase,

the number of degrees of freedom has to be infinite.

No definition of the term chaos is universally accepted yet, but almost everyone would

agree on the three ingredients used in the following working definition:

C'!s I. .s is periodic long-term behavior in a deterministic system that exhibits sensitive

dependence on initial conditions.

*"Aperiodic long-term b. !! lit.-! I means that there are trajectories which do not settle
down to fixed points, periodic orbits, or quasi-periodic orbits as t oc. For practical
reasons, we should require that such trajectories are not too rare. For instance, we could

insist that there he an open set of initial conditions leading to periodic trajectories, or
perhaps that such trajectories should occur with nonzero probability, given a random
initial condition.

* "Deterministic" means that the system has no random or noisy inputs or parameters.
The irregular behavior arises from the system's nonlinearity, rather than from noisy
driving forces.

* "Sensitive dependence on initial conditions" means that nearby trajectories separate
exponentially fast (i.e., the system has a positive Lyapunov exponent).

C'!s I. .s is also commonly called a strange attractor. The term attractor is also difficult

to define in a rigorous way. We want a definition that is broad enough to include all

the natural candidates, but restrictive enough to exclude the imposters. There is still

disagreement about what the exact definition should he.

Loosely -lo' I1:;0s an attractor is a set to which all neighboring trajectories converge.

Stable fixed points and stable limit cycles are examples. 1\ore precisely, we define an

attractor to be a closed set ,4 with the following properties:

1. A is an invariant set: any trajectory x(t) that starts in ,4 stays in ,4 for all time.

2. A attracts an open set of initial conditions: there is an open set U containing 24 such
that if x(0) E U, then the distance from x(t) to A tends to zero as t 00. This
means that ,4 attracts all trajectories that start sufficiently close to it. The largest
such U is called the basin of attraction of ,4.

3. A is minimal: there is no proper subset of ,4 that satisfies conditions 1 and 2.

Now we can define a strange attractor to be an attractor that exhibits sensitive dependence

on initial conditions. Strange attractors were originally called strange because they are

often fractal sets. ?- ue- ul .vis this geometric property is regarded as less important than

the dynamical property of sensitivity dependence on initial conditions. The terms chaotic

attractor and fractal attractor are used when one wishes to emphasize one or the other of

those aspects.

1.4 Brief Introduction to Fractal Geometry

Euclidean geometry is about lines, planes, triangles, squares, cones, spheres, etc. The

common feature of these different objects is regularity: none of them is irregular. Now let

us ask a question: Are clouds spheres, mountains cones, and islands circles? The answer is

obviously no. Pursuing answers to such questions, Alandelbrot has created a new branch of

science -fractal geometry.

For now, we shall be satisfied with an intuitive definition of a fractal: a set that shows

irregular but self-similar features on many or all scales. Self-similarity means a part of

an object is similar to other parts or to the whole. That is, if we view an irregular object

with a microscope, whether we enlarge the object by 10 times, or by 100 times, or even

by 1000 times, we ah-- i-n find similar objects. To understand this better, let us imagine

we were observing a patch of white cloud drifting .l.- li- in the sky. Our eyes were rather

passive: we were staring more or less at the same direction. After a while, the part of the

cloud we saw drifted away, and we were viewing a different part of the cloud. Nevertheless,

our feeling remains more or less the same.

Mathematically, fractal is characterized by a power-law relation, which translates to

a linear relation in log-log scale. For now, let us again resort to imagination -we were

walking down a wild and very j I__- d mountain trail or coastline. We would like to know

the distance covered by our route. Suppose our ruler has a length of e which could be

our step-size, and different hikers may have different step-sizes -a person riding a horse

has a huge step-size, while a group of people with a little child must have a tiny step-size.

The length of our route is

L = N(e) (1-3)

where NV(e) is the number of intervals needed to cover our route. It is most remarkable

that typically NV(e) scales with e in a power-law manner,

N~e) E-D e 0(1-4)

with D being a non-integer, 1 < D < 2. Such a non-integral D is often called the fractal

dimension -to emphasize the fragmented and irregular characteristics of the object under


Let us now try to understand the meaning of the non-integral D. For this purpose, let

us consider how length, area, and volume are measured. A common method of measuring

a length, a surface area, or a volume is by covering it with intervals, squares, or cubes

whose length, area, or volume is taken as the unit of measurement. These unit intervals,

squares, and volumes are called unit boxes. Suppose, for instance, that we have a line

whose length is 1. We want to cover it by intervals (hoxes) whose length is e. It is clear

that we need NV(e) ~ e-1 hoxes to completely cover the line. Similarly, if we want to cover

an area or volume by boxes with linear length e, we would need NV(e) ~ E-2 to cover the

area, or NV(e) ~ e-3 hoxes for the volume. Such D is called the topological dimension, and

takes on a value of one for a line, two for an area, and three for a volume. For isolated

points, D is zero. That is why a point, a line, an area, and a volume are called an 0 D,

1 D, 2 D, and 3 D objects, respectively.

Now let us examine the consequence of 1 < D < 2 for a j I__- d mountain trail. It is

clear that the length of our route increases as e becomes smaller, i.e., when e 0 L o.

T> he more concrete, let us visualize a race between the Hare and the Tortoise on a fractal

trail with D = 1.25. Assume that the length of the average step taken by the Hare is 16

times that taken by the Tortoise. Then we have

LHare 2 ,Lartoise

That is, the Tortoise has to run twice the distance of the Hare! Put differently, if you

were walking along a wild mountain trail or coastline and tired, slowing down your pace,

shrinking your steps, then you were in trouble! It certainly would be worse if you also got


1.5 Importance of Connecting Chaos and Random Fractal Theories

Multiscale signals behave differently depending upon which scale the data are

looked at. How can the behaviors of such complex signals on a wide range of scales

he simultaneously characterized? One strategy we envision is to use existing theories

synergistically, instead of individually. To make this possible, appropriate scale ranges

where each theory is most pertinent need to be identified. This is a difficult task, however,

since different theories may have entirely different foundations. For example, chaos theory

is mainly concerned about apparently irregular behaviors in a complex system that are

generated by nonlinear deterministic interactions of only a few numbers of degrees of

freedom, where noise or intrinsic randomness does not pIIl w any important role. Random

fractal theory, on the other hand, assumes that the dynamics of the system are inherently

random. Therefore, to fully characterize the behaviors of multiscale signals on a wide

range of scales, different theories, such as chaos and random fractal theories, need to be

integrated and even generalized.

1.6 Importance of the Concept of Scale

Complex systems often generate highly nonstationary and multiscale signals because

of nonlinear and stochastic interactions among their component systems and hierarchical

regulations imposed by the operating environments. Rapid accumulation of such complex

data in life sciences, systems biology, nano-sciences, information systems, and physical

sciences, has made it increasingly important to develop complexity measures that

incorporate the concept of scale explicitly, so that different behaviors of signals on varying

scales can be simultaneously characterized by the same scale-dependent measure. The

most ideal scenario is that a scale-dependent measure can readily classify different types

of motions based on analysis of short noisy data. In this case, one can readily see that the

measure will not only be able to identify appropriate scale ranges where different theories,

including information theory, chaos theory, and random fractal theory, apply, but also be

able to automatically characterize the behaviors of the data on those scale ranges. The

importance of the concept of scale will be further illustrated in later chapters.

1.7 Structure of this Dissertation

We emphasized that to develop new approaches for multiscale signal processing, it is

most desirable to integrate and even generalize different theories. Recently, the defining

property for deterministic chaos, exponential sensitivity to initial conditions (ESIC) has

been generalized to power-law sensitivity to initial conditions (PSIC). In C'!s Ilter 2,

we first introduce the new concept of PSIC. Then we extend the concept of PSIC to

high-dimensional case. We shall present a general computational framework for assessing

PSIC from real world data. We also apply the proposed computational procedure to

study noise-free and noisy logistic and Henon maps at the edge of chaos. we show that

when noise is absent, PSIC is hard to detect from a scalar time series. However, when

there is dynamic noise, motions around the edge of chaos, he it simply regular or truly

chaotic when there is no noise, all collapse onto the PSIC attractor. Hence, dynamic

noise makes PSIC observable. In (I Ilpter 3, we demonstrate that the new framework of

PSIC can readily characterize random fractal processes, including 1/ fP processes with

long-range-correlations and Levy processes. From our study in C'!s Ilters 2 and 3, it is

clear that the concept of PSIC not just bridges standard chaos theory and random fractal

theory, but in fact provides a more general framework to encompass both theories as

special cases. To illustrate the power of the framework of PSIC, in OsI Ilpters 4 and 5, we

apply it to study two important but challenging problems: Internet dynamics and sea

clutter radar returns, respectively. Both are outstanding examples of complex dynamical

systems with both nonlinearity and randomness. The nonlinearity of Internet dynamics

comes from AllMD (additive increase and multiplicative decrease) operation of TCP

(Transmission Control protocol), while the stochastic component comes from the random

user behavior. Sea clutter is a nonlinear dynamical process with stochastic factors due

to interference of various wind and swell waves and to local atmospheric turbulence.

We shall show that the new theoretical framework of PSIC can effectively characterize

the complicated Internet dynamics as well as to robustly detect low observable targets

from sea clutter radar returns with high accuracy. In Chapter 6, we shall introduce

the basic theory of an excellent multiscale measure -the scale-dependent Lyapunov

exponent (SDLE), and develop an effective algorithm to compute the measure. To

understand the SDLE as well as appreciate its power, we apply it to classify various types

of complex motions, including truly low-dintensional chaos, noisy chaos, noise-induced

chaos, processes defined by PSIC, and complex motions with chaotic behavior on small

scales but diffusive behavior on large scales. We also discuss how the SDLE can help

detect hidden frequencies in large scale orderly motions. Finally we suninarize our work in


The fornialism of nonextensive statistical niechanics (NESM) [2] has found numerous

applications to the study of systems with long-range-interactions [3-6] and niultifractal

behavior [7, 8], and fully developed turbulence [8-10], among many others. In order to

characterize a type of motion whose complexity is neither regular nor fully chaotic/randons

recently the concept of exponential sensitivity to initial conditions (ESIC) of deterministic

chaos has been generalized to power-law sensitivity to initial conditions (PSIC) [7, 11].

Mathematically, the formulation of PSIC closely parallels that of NESM. PSIC has been

applied to the study of deterministic 1-D logisticlike maps and 2-D Henon niap at the

edge of chaos [7, 11-20], yielding considerable new insights. In this chapter, we first

briefly describe some background knowledge about chaos. In particular, we will discuss a

direct dynamical test for deterministic chaos developed by Gao and Zheng [21, 22]. This

method offers a more stringent criterion for detecting low-dintensional chaos, and can

simultaneously monitor motions in phase space at different scales. Then we introduce the

more generalized concept of PSIC, and extend the new concept to high-dintensional case.

Specifically, we present a general computational framework for assessing PSIC in a time

series. We also apply the proposed computational procedure to study two model systems:

noise-free and noisy logistic and Henon maps at the edge of chaos.

2.1 Dynamical Test for Chaos

C'!s I. .s is also coninonly called a strange attractor. Being an I11 .. 1 ..1-~", the

trajectories in the phase space are bounded. Being -lI", the nearby trajectories,

on the average, diverge exponentially fast. The latter property is characterized hv the

exponential sensitivity to initial conditions (ESIC), and can he niathentatically expressed

as follows. Let d(0) be the small separation between two arbitrary trajectories at time

0, and d(t) be the separation between them at time t. Then, for true low-dintensional

deterministic chaos, we have

d(t) ~ d(0)ex' (2-1)

where X1 is positive and called the largest Lyapunov exponent. Due to the boundedness

of the attractor and the exponential divergence between nearby trajectories, a strange

attractor typically is a fractal, characterized by a simple and elegant scaling law as defined

by Eq. 1-4. A non-integral fractal dimension of an attractor contrasts sharply with the

integer-valued topological dimension (which is 0 for finite number of isolated points, 1 for

a smooth curve, 2 for a smooth surface, and so on).

Conventionally, it has been assumed that a time series with an estimated positive

largest Lyapunov exponent and a non-integral fractal dimension is chaotic. However, it has

been found that this assumption may not he a sufficient indication of deterministic chaos.

One counter-example is the so-called 1/ fP processes. Such processes have spectral density

S(f ) f -(2H+1) (2-2)

where 0 < H < 1 is sometimes called the Hurst parameter. This type of processes

will be further discussed in much depth in ChI Ilpter 3. A trajectory formed by such a

process has a fractal dimension of 1/H. As we shall explain later, with most algorithms

of estimating the largest Lyapunov exponent, one obtains a positive number for the

"Lyapunov expei ill ~! Therefore, such processes could be interpreted as deterministic

chaos .

Due to the ubiquity of noise, experimental data are ahr-l- :- corrupted by noise to some

degree. mathematically speaking, a noisy system, no matter how weak the noise is, has

infinite dimensions. Experimentally speaking, one would be more interested in a certain

range of finite scales. If the noise is very weak, then its influence on the dynamics may be

limited to very small scales, leaving the dynamics on finite scales deterministic-like.

The above discussions motivate us to employ a more sophisticated method, the

direct dynamical test for deterministic chaos [21, 22], to determine whether the data

under investigation is truly chaotic or not. The method offers a more stringent criterion

for low-dimensional chaos, and can simultaneously monitor motions in phase space at

different scales. The method has found numerous applications in the study of the effects

of noise on dynamical systems [23-27], estimation of the strength of measurement noise

in experimental data [28, 29], pathological tremors [30], shear-thickening surfactant

solutions [31], dilute sheared aqueous solutions [32], and serrated plastic flows [33].

Now let us briefly describe the direct dynamical test for deterministic chaos [21, 22].

Given a scalar time series, x(1), x(2),..., x(NV), (.I--II11.11.- for convenience, that they have

been normalized to the unit interval [0,1]), one first constructs vectors of the following

form using the time delay embedding technique [34-36]:

1K = [x(i), x(i + L), ..., x(i + (m 1)L)] (2-3)

with m being the embedding dimension and L the delay time. For example, when m = 3

and L = 4, we have VI = [x(1), x(5), x(9)], VW = [x(2), x(6), x(10)], and so on. For the

analysis of purely chaotic signals, m and L has to be chosen properly. This is the issue of

optimal embedding (see [21, 22, 37] and references therein). After the scalar time series is

properly embedded, one then computes the time dependent exponent A(k) curves

A(k) = n +k-V (2-4)

with r < ||%~ Vy || < r + ar, where r and Ar are prescribed small distances. The angle

brackets denote ensemble averages of all possible pairs of (1M, Vj). The integer k, called

the evolution time, corresponds to time kbt, where 6t is the sampling time. Note that

geometrically (r, r + Ar) defines a shell, and a shell captures the notion of scale. The

computation is typically carried out for a sequence of shells. Comparing Eq. 2-4 with

Eq. 2-1, one notices that ||%~+k Vy, || phI i-4 the role of d(t), while ||%~ Vy1 || phi-4 the role

of d(0). Figure 2-1(a) shows the A(k) vs. k curves for the chaotic Lorenz system driven by

stochastic forcing:

=-16(x y) + Drll(t)

= -xz + 45.92x y + Drl2(t) (2-5)

= xy 4z + Dr/3(

with < rli(t) > 0,- < i6i,(')> b(t-t'), i, j = 1, 2, 3. Note that D2 is the variance of

the Gaussian noise terms, and D = 0 describes the clean Lorenz system. We observe from

Fig. 2-1(a) that the A(k) curves are composed of three parts. These curves are linearly

increasing for 0 < k < k,. They are still linearly increasing for k, < k < ky, but with a

slightly different slope. They are flat for k > k. Note that the slope of the second linearly

increasing part gives an estimate of the largest positive Lyapunov exponent [21, 22]. k,

is related to the time scale for a pair of nearby points (1M, Vyj) to evolve to the unstable

manifold of 1M or Vyj. It is on the order of the embedding window length, (m 1)L. kg

is the prediction time scale. It is longer for the A(k) curves that correspond to smaller

shells. The difference between the slopes of the first and second linearly increasing parts

is caused by the discrepancy between the direction defined by the pair of points (1K, Vyj)

and the unstable manifold of 1M or Vyj. This feature was first observed by Sato et al. [38]

and was used by them to improve the estimation of the Lyapunov exponent. The first

linearly increasing part can be made smaller or can even be eliminated by adjusting

the embedding parameters such as by using a larger value for m. Note that the second

linearly increasing parts of the A(k) curves collapse together to form an envelope. It is

this very feature that forms the direct dynamical test for deterministic chaos [21, 22].

This is because the A(k) curves for noisy data, such as white noise or 1/ f processes,

are only composed of two parts, an increasing (but not linear) part for k < (m 1)L

and a flat part [21, 22]. Furthermore, different A(k) curves for noisy data separate from

each other, hence an envelope is not defined. We also note, most other algorithms for

estimating the largest Lyapunov exponent is equivalent to compute A(k) for r < ro, where

ro is selected more or less arbitrarily, then obtain A(k)/kbt, for not too large k, as an

estimation of the largest Lyapunov exponent. With such algorithms, one can obtain a

'pm enLyapunov exponent for white noise and for 1/ f processes. However, the value of

the Lyapunov exponent estimated this way sensitively depends on the parameter ro chosen

in the computation, hence, typically is different for different researchers. We thus observe

a random element here! With these discussions, it should be clear why 1/ f processes may

be interpreted as chaotic.

One can expect that the behavior of the A(k) curves for a noisy chaotic system lies

in between that of the A(k) curves for a clean chaotic system and that of the A(k) curves

for white noise or for 1/ f processes. This is indeed so. An example is shown in Fig. 2-1(b)

for the noisy Lorenz system with D = 4. We observe from Fig. 2-1(b) that the A(k)

curves corresponding to different shells now separate. Therefore, an envelope is no longer

defined. This separation is larger between the A(k) curves corresponding to smaller shells,

indicating that the effect of noise on the sniall-scale dynamics is stronger than that on the

large-scale dynanxics. Also note that k, + k, is now on the order of the embedding window

length, and is almost the same for all the A(k) curves. With stronger noise (D > 4), the

A(k) curves will be more like those for white noise [21, 22].

Finally, we consider the Mackey-Glass delay differential system [39]. When a=

0.2, b = 0.1, c = 10, F = 30, it has two positive Lyapunov exponents, with the largest

Lyapunov exponent close to 0.007 [21].

axr(t +F
dxr/dt = bxr(t) (2-6)
1 + r(t +F)"

Having two positive Lyapunov exponents while the value of the largest Lyapunov exponent

of the system is not much greater than 0, one might he concerned that it may be difficult

to deal with this system. This is not the case. In fact, this system can he analyzed as

straightforwardly as other dynamical systems including the Logistic nmap, Henon niap and

the Rossler system. An example of the A(k) vs. k curves for the chaotic Alackey-Glass


O 50



Time-dependent exponent A(k) vs. evolution time k curves for (a) clean and
(b) noisy Lorenz system. Six curves, from bottom to top, correspond to shells
(2-(i+1)/2, 2-i/2) With i 8, 9, 10, 11, 12, and 13. The sampling time for the
system is 0.03 sec, and embedding parameters are m = 4, L = 3. 5000 points
are used in the computation.

Figure 2-1.

100 150

0 50 100 150
Evolution timne k

3.5 -




0 100 200 300 400 500 600
Evolution time k

Figure 2-2. Time-dependent exponent A(k) vs. evolution time k curves for the
Mackey-Glass system. Nine curves, from bottom to top, correspond to shells
(2-(i+1)/2, 2-i/2) With i 5, 6, and 13. The computation was done with
m = 5, L = 1, and 5000 points sampled with a time interval of 6.

delay differential system is shown in Fig. 2-2, where we have followed [21] and used

m = 5, L = 1, and 5000 points sampled with a time interval of 6. Clearly, we observe that

there exists a well defined common envelope to the A(k) curves. Actually, the slope of the

envelope estimates the largest Lyapunov exponent. This example illustrates that when

one works in the framework of the A(k) curves developed by Gao and Zheng [21, 22], one

does not need to be very concerned about non-uniform growth rate in high-dimensional

systems .

2.2 General Computational Framework for Power-law Sensitivity to Initial
Conditions (PSIC)

To understand the essence of power-law sensitivity to initial conditions (PSIC),

let us first consider the 1-D case and consider Eq. 2-7, where Ax(0) is the infinitesimal

discrepancy in the initial condition, Ax(t) is the discrepancy at time t > 0.

((t) =lim (2-7)

When the motion is chaotic, then

((t) = exlt (2-8)

((t) satisfies

Tsallis and co-workers [7, 11] have generalized Eq. 2-9 to

di = A,((t)V (2-10)

where q is called the entropic index, and X, is interpreted to be equal to K,, the

generalization of the K~olmogorov-Sinai entropy. Eq. 2-10 defines the PSIC in the 1-D
case. Obviously, PSIC reduces to ESIC when q 1 The solution to Eq. 2-10 is

s'(t) =[1 +(1 q)Agt] 'il) (2-11)

When t is large and q / 1, ((t) increases with t as a power-law,

((t) Ct 7 -')(2-12)

where C = [(1 q) A,] 1 (-V). For Eq. 2-12 to define an unstable motion with A, > 0, we

must have q < 1. Later we shall map different types of motions to different ranges of q.

To apply PSIC to the analysis of time series data, one can first construct a phase

space by constructing embedding vectors 1M as defined by Eq. 2-3. Eq. 2-11 can then be

generalized to high-dimensional case [40],

((t) = Ilimo~l =lvo~l lv t 1 + (1 q)X()t (2 13

where | |a V(0) | is the infinitesimal discrepancy between two orbits at time 0, | |a V(t) | is

the distance between the two orbits at time t > 0, q is Ithe entropic index, and A i) is the

first q-Lyapunov exponent, corresponding to the power-law increase of the first principle

axis of an infinlitesimral ball in thle phase space. AI1 mayi not be equal to K,. TIhisi is

understood by recalling that for chaotic systems, the K~olmogorov-Sinai entropy is the sum

of all the positive Lyapunov exponents. We conjecture a similar relation may hold between

the q-Lyapunov exponents and K. When there are multiple unstable directions, then in

general A I) maty- not be equal to Ki,. Whe~n t is large andc q / 1, Eq. 2-13 again gives a

power-law increase of ((t) with t. Note that under the above framework, the motion may

not be like fully developed chaos, thus not ergodic [41].

We now consider the general computational framework for PSIC. Given a finite

time series, the condition of ||av(0)|| 0 cannot be satisfied. In order to find the law

governing the divergence of nearby orbits, one thus has to examine how the neighboring

points, (1K, Vj), in the phase space, evolve with time, by forming suitable ensemble

averages. Notice that if (1MI, by) and (1%2, V/j2) arT tWO pairs of nearby points, when

| |%1 Vjll | |
|I |%2+t Vj2+tl CRllllOt be simply averaged to provide estimates for q anld Af ). Inl fact, it

would be most convenient to consider ensemble averages of pairs of points (1M, V4j) that

all fall within a very thin shell, rl < ||%~ g||1 < r2, Where rl and T2 arT ClOSe. These

arguments so-----~ -1 that the time-dependent exponent curves defined by Eq. 2-4 provide a

natural framework to assess PSIC from a time series. This is indeed so. In fact, we have

In I(t) a A(t) = In (||+ +t) (2-14)

Now, by the discussions in Sec. 2.1, it is clear that PSIC is a generalization of ESIC:

as long as the A(k) curves from different shells form a linear envelope, then q = 1 and the

motion is chaotic. The next question is: Does PSIC also include random fractals as special

cases? The answer is yes. It will be given in the next chapter. There, we will also gain

a better undlerstanding of the meaning of A l. The re~st of this chapter shows examples

of characterizing the edge of chaos by the new framework of PSIC. Specifically, we study

time series generated from noise-free and noisy logistic and Henon maps.

2.3 Characterizing Edge of Chaos by Power-law Sensitivity of Initial
Conditions (PSIC)

Let us first examine the logistic map around the edge of chaos. Considering that

dynamic noise is often an important component of a system (e.g., spontaneous emission

noise in semiconductor lasers [26], and physiological noise in biological systems [30]), we

study both the deterministic and the noisy logistic map:

X,+1 = ax,(1 x,) + arl, (2-15)

where a is the bifurcation parameter and rlo is a white Gaussian noise with mean zero

and variance 1. The parameter a characterizes the strength of noise. For the clean system

(a = 0), the edge of chaos occurs at the accumulation point, a, = 3.569945672 We

shall study three parameter values, al = a, 0.001, a,, and a2 = a, + 0.001. When

noise is absent, al corresponds to a periodic motion with period 25, while a2 COTTOSponds

to a truly chaotic motion. We shall only study transient-free time series. In Figs. 2-3(a-c),

we have plotted the A(t) vs. t curves for parameter values al, a,, and a2, TOSpectively.

We observe from Fig. 2-3(a) that the variation of A(t) with t is periodic (with period 16,

which is half of the period of the motion) when the motion is periodic. This is a generic

feature of the A(t) curves for discrete periodic attractors, when the radius of the shell is

larger than the smallest distance between two points on the attractor (when a periodic

attractor is continuous, A(t) can be arbitrarily close to 0). It has been found by Tsallis

and co-workers [11] that at the edge of chaos for the logistic map, ((t) is given by Eq. 2-11

with q m 0.2445. Surprisingly, we do not observe such a divergence in Fig. 2-3(b). In

fact, if one plots A(t) vs. In t, one only observes a curve that increases very slowly (similar

to that shown in Fig. 2-4(a)). The more interesting pattern is the one that is shown in

Fig. 2-3(c), where we observe a linearly increasingly A(t) vs. t curve. In fact, shown in

Fig. 2-3(c) are two such curves, corresponding to two different shells. Very interestingly,

the two curves collapse to form a common envelope in the linearly increasing part of the






0 1 23 456

Time-dependent exponent A(t) vs. t curves for time series generated from the
noise-free logistic map with (a) al a, 0.001, where the motion is periodic
with period 2s, (b) a, 3.569945672 and (c) a2 a, + 0.001, where the
motion is chaotic. Plotted in (d-f) are A(t) vs. In t curves for the noisy logistic
map with o- = 0.001. Very similar results were obtained when o- = 0.0001.
Shown in (c-f) are actually two curves, corresponding to two different shells.
104 pOintS Were used in the computation, with embedding parameters
m = 4, L = 1. However, so long as m > 1, the results are largely independent
of embedding. When m = 1, the A(t) curves are not smooth, and the
estimated 1/(1 q) value is much smaller than the theoretical value.

Figure 2-3.

curve. The slope of the envelope gives a good estimate of the largest positive Lyapunov

exponent. This is a generic feature of chaos [21, 22], as explained in Sec. 2.1. Since the

chaos studied here is close to the edge of chaos, the curves shown in Fig. 2-3(c) are less

smooth than those reported earlier [21-23, 26].

Why cannot the theoretical prediction of PSIC at the edge of chaos for the logistic

map be observed from a clean time series? In a recently published very interesting and

1.5 -


0 50 100 150 200 250
2.0 '


0.0 .
0 50 100 150 200 250

0 50 100 150 200 250

insightful paper, Beck [42] so__~- -;- that dynamic noise may be of fundamental importance

to the Tsallis NESM. May dynamic noise pl ai similarly significant role for PSIC? As is

shown below, the answer is yes. In Figs. 2-3(d-f), we have shown the A(t) vs. In t curves

for the three parameters considered, with noise strength al = 0.001. In fact, shown

in each figure are two curves, corresponding to two different shells. They parallel with

each other. The slopes of those curves are about 1.20, close to the theoretical value of

1/(1 0.2445) a 1.32. While it is very satisfactory to observe PSIC at the edge of chaos,

it is more thrilling to observe the collapse of regular as well as chaotic motions onto the

PSIC attractor around the edge of chaos. This signifies the stableness of PSIC when there

is dynamic noise. It is important to emphasize that the results shown in Figs. 2-3(d-f)

are largely independent of the noise strength, so long as noise is neither too weak nor too

strong. For example, very similar results have been observed with a2 = 1/10 = 0.0001.

Before we move on to discuss the Henon map near the edge of chaos, let us comment

on the difference between the ESIC and the PSIC. When the ESIC is the case, the A(t)

vs. t curves are straight for T, < t < T,, where T, is a prediction time scale, and T, is a

time scale for the initial separation to evolve to the most unstable direction. When A(t)

vs. t curves corresponding to different shells are plotted together, they collapse together

for a considerable range of t. When the PSIC is the case, the A(t) vs. In t curves are

fairly straight, and the A(t) vs. In t curves corresponding to different shells separate and

parallel with each other. This feature is the fundamental reason that ensemble averages

are most conveniently formed by requiring neighboring pairs of points to all fall within

a thin shell. For example, if two distinct thin shells, described by rl < ||Xi Xj| < T2 r

and T2 ||II Vy|| < T 3, arT joined together to form a single thicker shell, described by

rl < ||V~ Vy||1 < T3, and one averages the separation between Xi and Xj before taking

the logarithm, then the resulting slope between A(t) and In t, assumed to be still linear,

will have to be smaller than that estimated from the two thin shells separately. At this

point it is also worth noting that the linearly increasing part of the A(t) vs. In t curves

may only contain a small interval of t. Actually, for the simple logistic nmap, A(t) saturates

around i s 20, when the radius of the shell is about 10-4 and the length of the time series

is 104. For experimnertal tilie series, the t for A(t) to saturate may be even smaller (we

suspect it would be about 10 samples -irrespective of the sampling time, because of

the logarithmic timre scale). Tob get more accurate estimnates of the parametersi q and A ',

it may be better to fit ((t) using Eq. 2-13 when dealing with experimental time series.

In fact, we have found that if we do so, the estimated 1/(1 q) value from Figs. 1(d-f)

increases to about 1.27, much closer to the theoretical value of 1.32.

Next let us consider the Henon nmap, where ty., and ty, are white Gaussian noise, and

are uncorrelated with each other. The parameter a measures the strength of the noise.

r,z+l = 1 axr + Utz + wif, (77)

Un+l = br,z + aty,(iv)

Tirnakli [20] studied the deterministic version of this nmap at the edge of chaos, for

parameter values a,. = 1.40115518 b = O; a,. = 1.39966671 b = 0.001; and

a,. = 1.:386:37288 b = 0.01. We have studied both the noise-free and noisy nmap

for parameter values listed above, and found very similar results to those presented in

Fig. 2-:3. In Figs. 2-4(a,h), we have plotted A(t) vs. In t curves for the noise-free and noisy

nmap for a,. = 1.:386:37288 b = 0.01. As can he observed clearly, Fig. 2-4(b) is very

similar to Figs. 2-:3(d-f), while the slope for the curves in Fig. 2-4(a) is much smaller than

that in Fig. 2-4(b). Again, we have observed (but not shown by figures here, since they are

very similar to Figs. 2-:3(d-f) and Fig. 2-4(b)) that dynamic noise makes the regular and

chaotic motions to collapse onto the PSIC attractor for parameters around those definingf

the edge of chaos, and that such transitions are largely independent of the noise strength.

To suninarize, we have described an easily intplenientable procedure for computationally

examining PSIC front a time series. By studying two model systems: noise-free and noisy

logistic and Henon maps near the edge of chaos, we have found that when there is no

2.0 (a)''""'""''"" 6 (b

S0.5 -

0.0 -1 2

-0.5- -1 1
1.0O ~......... ......... ......... ......... .......... ............ O .
0123456 0123456
Int Int

Figure 2-4. Time-dependent exponent A(t) vs. t curves for time series generated from (a)
the noise-free Henon map with ac = 1.38637288 b = 0.01. Plotted in (b)
are two A(t) vs. In t curves (corresponding to two different shells) for the noisy
map with a 0.001. Similar results were obtained with a 0.0001. 104 pOintS
were used in the computation, with embedding parameters m = 4, L = 1.

noise, the PSIC attractor cannot be observed from a scalar time series. However, when

dynamic noise is present, motions around the edge of chaos, be it simply regular or

truly chaotic when there is no noise, all collapse onto the PSIC attractor. While these

examinations signify the ubiquity of PSIC, they also highlight the importance of dynamic

noise. The existence of the latter is perhaps the very reason that truly chaotic time series

can seldom be observed.


In C'!s Ilter 2, we have discussed that PSIC is a generalization of ESIC: as long as

the A(t) curves from different shells form a linear envelope, then q = 1 and the motion

is chaotic. We have also shown that edge of chaos can be well characterized by PSIC.

In this chapter, we shall study two 1!! I r~~ types of random fractals: 1 fP processes with

long-range-correlations and Levy processes. We shall show both analytically and through

numerical simulations that both types of processes can be readily characterized by PSIC.

3.1 Characterizing 1/ f Processes by Power-law Sensitivity to Initial
Conditions (PSIC)

A continuous-time stochastic process, X = {X(t), t > 0}, is said to be self-similar if

X(At) d AHX(t), t > 0 (3-1)

for A > 0, O < H < 1, where u~~denotes equality in distribution. H is called the

self-similarity parameter, or the Hurst parameter.

Before proceeding on, we note that, more rigorously speaking, processes defined by

Eq. 3-1 should be called self-affine processes [43] instead of self-similar processes, since X

and t have to be scaled differently to make the function look similar. A simple physical

explanation for this is that the units for X and t are different. In this dissertation,

however, we will continue to call such processes self-similar, to follow the convention in

some of the engineering disciplines.

The following three properties can be easily derived from Eq. 3-1:

E [X(t)] Mean (3-2)

Var [X(At)]
Var[X(t)]= 2HVariance (3-3)

R, (t, s) = 2HAutocorrelation (3-4)

Note that Var[X(t)] = R,(t, t); hence, Eq. 3-3 can be simply derived from Eq. 3-4; we

have listed it as a separate equation for convenience of future reference. If we consider

only second-order statistics, or if a process is Gaussian, Eq. 3-2 to Eq. 3-4 can be used

instead of Eq. 3-1 to define a self-similar process.

A very useful way of describing a self-similar process is by its spectral representation.

Strictly -p.' I1:;14 the Fourier transform of X(t) is undefined, due to the nonstationary

nature of X(t). One can, however, consider X(t) in a finite interval, ;?i, O < t < T:

X (t, T) = t)0< t< T (3-5)

and take the Fourier transform of X(t, T):

F(f, T) = lX(t)e-2xjftd (36)

|F(f, T)|12df is the contribution to the total energy of X(t, T) from those components with

frequencies between f and f + df. The (average) power spectral density (PSD) of X(t, T)
is then

S(f T) = FfT

and the spectral density of X is obtained in the limit as T oo

S( f) = lim S( f ,T)

Noting that the PSD for A-HX(At) is

A-2H limX(t-2ft _-2-S/)
T->oo AT$ I 'o X.X)ZiLt 2 X"-sfX

and that the PSD for A-HX(At) and X(t) are the same [44], we have

S(f) = S(f/A)A-2H-1

The solution to the above equation is

S( f) ~ f -P (3-7)


/7 = 2H + 1 (3-8)

Processes with power spectral densities as described by Eq. 3-7 are called 1/ f

processes. Typically, the power-law relationships that define these processes extend over

several decades of frequency. Such processes have been found in numerous areas of science

and engineering. Some of the older literatures on this subject can be found, for example,

in Press [45], Bak [46], and Wornell [47]. Some of the more recently discovered 1/ f

processes are in traffic engineering [48-50], DNA sequence [51-53], human cognition [54],

ambiguous visual perception [55, 56], coordination [57], posture [58], dynamic images [59,

60], and the distribution of prime numbers [61], among many others. It is further observed

that principle component analysis of such processes leads to power-law decaying eigenvalue

spectrum [62].

Two important prototypical models for 1/ f processes are the fractional Brownian

motion (fBm) processes and the ON/OFF intermittency with power-law distributed ON

and OFF periods. Below, we apply the concept of PSIC to both types of processes.

3.1.1 Characterizing Fractional Brownian Motion (fBm) Processes by Power-
law Sensitivity to Initial Conditions (PSIC)

We first briefly introduce Brownian motion (Bm). Brownian motion is defined as a

(nonstationary) stochastic process B(t) that satisfies the following criteria:

1. All Brownian paths start at the origin: B(0) = 0.

2. For 0 < 1 independent.

3. For all (s, t) > 0, the variable B(t + s) B(t) is a Gaussian variable with mean 0 and
variance s.

4. B(t) is a continuous function of t.
We may infer from the above definition that the probability distribution function of B

is given by:

Pt { B(t + ) B(t) ]
This function also satisfies the scaling property:

P{[B(t + s) B(t)] < x} = P{B(At + As) B(At) < X1/2X (310)

In other words, B(t) and A-1/2B(At) have the same distribution. Thus, we see from

Eq. 3-1 that Bm is a self-similar process with Hurst parameter 1/2.

Suppose we have measured B(tl), B(t2) 2~ 1 > 0. What can we ;?i about

B(s), 1 < s < 2? The answer is given by the Levy interpolation formula:

(S tl) (t2 s) 8 1l]~ 1/21
B(s) = B(tl) +[B(t2) 1 ~t
(t2 1l 2a 1

where W is a zero-mean and unit-variance Gaussian random variable. The first two terms

are simply a linear interpolation. The third term gives the correct variance for B(s)-B(t )

and B(t2) B 8), Which are s 1 and t2 8, TOSpectively.

A popular way of generating Bm is by the following random midpoint displacement

method, which is essentially an application of the Levy interpolation formula: suppose we

are given B(0) = 0 and B(1) as a sample of a Gaussian random variable with zero mean
and variance O.2. We can obtain B(1/2) from the following formula:

B(1/2) B(0) = n[B(1) B(0)] + D1

where D1 is a random variable with mean 0 and variance 2-2 2. D1 is simply the third

term of Eq. 3-11, and the coefficient 2-2 equals (t2 s) 8 1 2~t 1l), With tl = 0, t2 =

1, a = 1/2. Similarly, we have

B(1/4) B(0) = ,[B(1/2) B(0)] + D2

where D2 1S a random variable with zero mean and variance 2-3 2. We can apply the

same idea to obtain B(3/4). In general, the variance for D, is 2-("+l)a2.

We note that a Brownian path is not a differentiable function of time. Heuristically,

this can be understood in this way: consider the variable, B(t + s) B(t), with variance

s. Its standard deviation, which is a measure of its order of magnitude, is ~ 2. Thus, the

derivative of B at t behaves like the limit Z/S = S-1/2 aS s 0.

Although B(t) is almost surely not differentiable in t, symbolically one still often


B(t) = w~(T)dv (3-12)

where w(t) is stationary white Gaussian noise, and extends the above equation to t < 0

through the convention,

JO~ Jt
It should be understood that integrals with respect to the differential element w(t)dt

should be interpreted more precisely as integrals with respect to the differential element

dB(t), in the Riemann-Stieltjes integral sense. This has profound consequence when

numerically integrating Eq. 3-12.

To see how, let us partition [0, t] into a equally spaced intervals, at = t/n. Since w(t)

are Gaussian random variables with zero mean and unit variance, in order for B(t) to have

variance t, we should have

i= 1
That is, the coefficient is (At)1/2 inStead of at, as one might have guessed. Typically,

at at. Thus, if one incorrectly uses at instead of (At)1/2, One is

severely underestimating the strength of the noise.

(a) White noise

(b) Brownian motion

Figure 3-1. White noise (a) and Brownian motion (b). Axes are arbitrary.

When the time is genuinely discrete, or the units of time can be arbitrary, one can
take at to be 1 unit, and Eq. 3-13 becomes

B(t) = w~7(i) (3 -14)
i= 1

Eq. 3-14 provides perhaps the simplest method of generating a sample of Bm. An example
is shown in Fig. 3-1.

Eq. 3-14 is also known as a random walk process. A more sophisticated random
walk (or jump) process is given by summing up an infinite number of jump functions:

Je(t) = AiP(t ti), where P(t) is a unit-step function

P(t) = O < (3-15)

and Ai, ti are random variables with Gaussian and Poisson distributions, respectively.

Next we discuss fractional Brownian motion (fBm) processes. A normalized fBm

process, Z(t), is defined as follows:

Z(t) WiH (t > 0) (3-16)

where W is a Gaussian random variable of zero mean and unit variance, and H is the

Hurst parameter. When H = 1/2, the process reduces to the standard Brownian motion

(Bm). It is trivial to verify that this process satisfies the defining Eq. 3-1 for a self-similar
stochastic process.

Fr-actional Brownian motion BH 1) is a Gaussian process with mean 0, stationary
increments, variance

E [(BH )2] t2H n7

and covariance:

E[B 8)H 2H 2PH _2HJ I

where H is the Hurst parameter. Due to its Gaussian nature, according to Eq. 3-2

to Eq. 3-4, the above three properties completely determine its self-similar character.

Fig. 3-2 shows several fBm processes with different H.

Roughly, the distribution for an fBm process is [63]:

BH~l TTii~ H-/2 (319)

where F(t) is the gamma function:

0/OO"t-e~d 7

It is easy to verify that fBm satisfies the following scaling property:

P(BHt + ) BH~t I) = P 8H(X + ) BXH XHX} (3-20)

In other words, this process is invariant for transformations conserving the similarity

variable X/tH

(c) H=0.75

(d) H=0.90

Figure 3-2. Several fBm processes with different H.

The integral defined by Eq. 3-19 diverges. The more precise definition is

BH t) BH(0) =r( K1 -7dB0 (3-21)

where the kernel K(t -r) is given by [63]

K(t r) =l l (3-22)

The increment process of fBm, Xi = BH(( + )t BH i 1), i > i, Where at can

be considered a sampling time, is called fractional Gaussian noise (fGn) process. It is a

zero mean stationary Gaussian time series. Noting that

E(XXik) = E {[BH(Z + )t BH 8t~[ H ((i + 1 + k) at) BH ((i + k) at)]}

by Eq. 3-18, one obtains the autocovariance function y(k) for the fGn process:

q~~~~k)2` =. E(4ik/(X)= ( )2H- 2k2H + |k- 1|2H, k >0 (3 23)

Notice that the expression is independent of at. Therefore, without loss of generality, one

could take at = 1. In particular, we have

y(1) = 22-2

Let us first note a few interesting properties of y(k):

(i) When H = 1/2, y(k) = 0 for k / 0. This is the well-known property of white
Gaussian noise.

(ii) When 0 < H < 1/2, y(1) < 0.

(iii) When 1/2 < H < 1, y(1) > 0.

Properties (ii) and (iii) are often termed anti-persistent and persistent correlations [43],


Next, we consider the behavior of y(k) when k is large. Noting that when |x|

(1 + x)" a 1 + a~x + x

we have, when k > 1,

(k + 1)2H + |k~ 1|2H = k~2H[(1 + 1/k)2H + (1 1/k~)2H] = k~2H[2 + 2H(2H 1)k-2]


y(k) ~ H(2H 1)k2H-2, Sk 00(3-24)

when H / 1/2. When H = 1/2, y(k) = 0 for k > 1, the Xi's are simply white noise.

Now we are ready to characterize 1/ f processes in the framework of PSIC. Recalling

that the defining property for a 1 fP process is that its variance increases with t as t2H

Irrespective of which embedding dimension is chosen, this property can be translated to

I(t) = tH (3-25)

Therefore ,
dt)= HtH-1 (3-26)

Expressing t in terms of (, we have

= H( -$ H3-27)

Comparing with the defining equation of PSIC (Eq. 2-10), we find that

1 2
q =1- 1- (3-28)
H P-1

X( ) = H= (3-29)

Noticing 0 < H < 1, from Eq. 3-28, we have -oo < q < 0.

Computationally, the key is to demonstrate behaviors described by Eq. 3-25. This

can be readily done by calculating the time dependent exponent curves defined by

Eq. 2-14 (see also Eq. 2-4). In Fig. 3-3 we have shown a few examples for the fBm

processes with H = 0.25, 0.5, and 0.75. Clearly, the slopes of the straight lines correctly

estimate the H parameter used in the simulations.

3.1.2 Characterizing ON/OFF Processes with Pareto distributed ON and
OFF Periods by Power-law Sensitivity to Initial Conditions (PSIC)

ON/OFF intermittency is a ubiquitous and important phenomenon. For example,

a queueing system or a network can alternate between idle and busy periods, a fluid

flow can switch from a turbulent motion to a regular (called laminar) one. Fig. 3-4

shows an example for an ON/OFF traffic model with three independent traffic sources

Xi(t), i = 1, 3, where each Xi(t) is a 0/1 renewal process.

Let us denote an ON period by 1 and an OFF period by 0. We study ON/OFF trains

where ON and OFF periods are independent and both have heavy-tailed distribution.


Figure 3-4. ON/OFF sources with NV = 3, X (t), X2(t), X3(t), and their summation S3 -)


H = 0.25
H = 0.50
H = 0.75



2 2.5 3 3.5 4


Figure 3-3. Time-dependent exponent A(k) vs. In k curves for fractional Brownian motion.


X ,(t)



X 2t)



A heavy-tailed distribution is commonly expressed as

f (x) ~ x-"-1,- x 00

Equivalently, one may write:

P X > x] ~ x-". x 0 o 3-30)

This expression is often called the complementary cumulative distribution function

(CCDF) of a ]!. l.-i--r I.l:d distribution. In fact, it is more popular for expressing a

heavy-tailed distribution, since it emphasizes the tail of the distribution. It is easy to

prove that when p < 2, the variance and all higher than 2nd-order moments do not exist.

Furthermore, when p < 1, the mean also diverges.

When the power-law relation extends to the entire range of the allowable x, we have

the Pareto distribution:

P[X J > ~ I x]= >biU > 2> 2> (3-31)

where p and b are called the shape and the location parameters, respectively. It can be

proven [64] that when 1 < p < 2, the Hurst parameter of the ON/OFF processes with

Pareto distributed ON and OFF periods is given as

H = (3 p-) /2 (3-32)

Our purpose here is to check whether the ON/OFF processes with Pareto distributed

ON and OFF periods can be characterized by PSIC. In Fig. 3-5 we have shown a few

examples for p = 1.2, 1.6, and 2.0. Noticing Eq. 3-32, we see that the slopes of the

straight lines again correctly estimate the H parameter. Also note that when 0 < p < 1,

1 < H < 1.5. Therefore, Eqs. (3-25) and (3-28) are still correct when 1 < H < 1.5. The

entire range of H = (3 p)/2, with 0 < p I 2, determines that -1 I q < 1/3 for Pareto

distributed ON/OFF processes.

4L = 2.0
4 = 1.6
4 = 1.2



.53 3.5 4 4.5

Figure 3-5. Time-dependent exponent A(k) vs. In k curves for ON/OFF processes with
Pareto distributed ON and OFF periods.

Before ending this section, we would like to emphasize that the possibility of

misinterpreting 1/ f processes being deterministic chaos never occurs if one works

under the framework of PSIC. Under this framework, one monitors the evolution of two

nearby trajectories. If two nearby trajectories diverge exponentially fast, then the time

series is chaotic, if the divergence increases in a power-law manner, then the trajectories

belong to 1/ f processes. These ideas can be easily expressed precisely. Let | |% Vy | |

be the initial separation between two trajectories. This separation is assumed to be not

larger than some prescribed small distance r. After time t, the distance between the two

trajectories will be ||%~+t Vj+t||. For true chaotic systems,

|| K~+t Vj+t|| oc ||K 1,- Vj ||ext

where A > 0 is called the largest Lyapunov exponent. On the other hand, for 1/ f


The latter equation thus provides another simple means of estimating the Hurst


3.2 Characterizing Levy Processes by Power-law Sensitivity to Initial
Conditions (PSIC)

We now consider Levy processes, another important type of random fractal

models that have found numerous applications [65-70]. There are two types of Levy

processes [71]. One is Levy flights, which are random processes consisting of many

independent steps, each step being characterized by a stable law, and consuming a

unit time regardless of its length. The other is Levy walks, where each step takes time

proportional to its length. A Levy walk can be viewed as sampled from a Levy flight

with a uniform speed. The increment process of a Levy walk, obtained by differencing

the Levy walk, is very similar to an ON/OFF train with power-law distributed ON

and OFF periods. Therefore, in the following, we shall not be further concerned about

it. We shall focus on Levy flights. Note that Levy flights, having independent steps,

are memoryless processes characterized by H = 1/2, irrespective of the value of the

exponent a~ characterizing the stable laws [71]. In other words, methods such as detrended

fluctuation analysis (DFA) [72] cannot be used to estimate the a~ parameter.

Now we define Levy flights more precisely. A (standard) symmetric ac-stable Levy

process {L,(t), t > 0} is a stochastic process with the following properties [73]:

1. L,(t) is almost surely 0 at the origin t = 0;

2. L,(t) has independent increments; and

3. L,(t) L,(s) follows an ac-stable distribution with characteristic function e-(t-s) lo~

where 0 < s < t < 00.

A random variable Y is called (strictly) stable if the distribution for CE is the same as

that for niggY

i= 1
where Yi, Y2, are independent random variables, each having the same distribution as

Y. From Eq. 3-33, one then finds that nVarY = n2/"VarY. For the distribution to be

valid, O < a~ < 2. When a~ = 2, the distribution is Gaussian, and hence, the corresponding

Levy flight is just the standard Brownian motion. When 0 < a~ < 2, the distribution is

heavy-tailed, P[X > x] ~ x-", x 00o, and has infinite variance. Furthermore, when

0 < a~ < 1, the mean is also infinite.

At this point, it is worth spending a few paragraphs to explain how to simulate an

ac-stable Levy process. The key is to simulate its increment process, which follows an

ac-stable distribution. The first breakthrough for the simulation of stable distributions

was made by K~anter [74], who gave a direct method for simulating a stable distribution

with a~ < 1 and p = 1. The approach was later generalized to the general case by

C'I I1!.1wers et al. [75]. We first describe the algorithm for constructing a standard stable

random variable X ~ S(a~, p; 1).

Theorem Simulation of stable random variables: Let 8 and W be independent

with 8 uniformly distributed on (-xr/2, xr/2), W exponentially distributed with mean 1.

For 0 < a~ < 2, -1 < p < 1, and 8o = ard i I[Sltan(;ac/2)]/a,

Z (cosina(00+8) cos[aeo+(a-1)8] 1a/ 34

has a S(a~, p; 1) distribution.

For symmetric stable distribution with P = 0, the above expression can be greatly

simplified. To simulate an arbitrary stable distribution S(a~, P, y, 5), we can simply take

the following linear transform,

Y yZ + 6 a rJ 1 (3-35)
yZ+( S"qln n+6 a~= 1

where Z is given by Eq. 3-34.

We have simulated a number of symmetric ac-stable Levy motions with different

a~. Figure 3-6 shows two examples of Levy motions with a~ = 1.5 and 1. The difference

between the standard Brownian motions and Levy motions lies in that Brownian motions

appear everywhere similarly, while Levy motions are comprised of dense clusters connected

by occasional long jumps: the smaller the a~ is, the more frequently the long jumps appear.

Let us now pause for a while and think where Levy flights-like esoteric processes

could occur. One situation could be this: a mosquito head to a giant spider net and got

stuck; it struggled for a while, and luckily, with the help of a gust of wind, escaped. When

the mosquito was struggling, the steps it could take would be tiny; but the step leading

to its fortunate escape must be huge. As another example, let us consider (American)

football games. For most of the time during a game, the offense and defense would be

fairly balanced, and the offense team may only be able to proceed for a few yards. But

during an attack that leads to a touchdown, the hero getting the touchdown often "flies"

tens of yards -somewhat he has escaped the defense. While these two simple situations

are only meant as analogies, remembering them could be helpful when reading research

papers searching for Levy statistics in various types of problems such as animal foraging.

At this point, we also wish to mention that Levy-flights-like patterns have been used as a

type of screen saver for Linux operating systems.

The symmetric ac-stable Levy flight is 1/a~ self-similar. That is, for c > 0, the

processes {L,(ct), t > 0} and {cl/oL,(t), t > 0} have the same finite-dimensional

distributions. By this argument as well as Eq. 3-33, it is clear that the length of the

(a) (b)

(c) (d)

(e) (f)

Figure 3-6. Examples for Brownian motions and Levy motions. 1-D and 2-D Brownian
motions (a,b) vs. Levy motions with a~ = 1.5 and 1 for (c,d) and (e,f),

o a =1.0
-4.5 0 a =1.5




06~.5 1 1.5 2 2.5 3

Figure 3-7. Time-dependent exponent A(k) vs. In k curves for Levy processes.

motion in a time span of at, AL(At), is given by the following scaling:

AL(At) oc Atl/o (3-36)

This contrasts with the scaling for fractional Brownian motion:

ABHat oc OC (33H

We thus identify that 1/a ]I phi- the role of H. Therefore,

q = 1 a~ (3-38)

A 0 = 1/al (3-39)

Noticing 0 < a~ < 2, from Eq. 3-38, we have -1 < q < 1.

We have applied the concept of PSIC to Levy processes with different a~. Examples

for Levy processes with a~ = 1 and 1.5 are shown in Fig. 3-7. The slopes of the straight

lines correctly estimate the values of a~ used in the simulations.


Time series analysis is an important exercise in science and engineering. One of the

most important issues in time series analysis is to determine whether the data under

investigation is regular, deterministically chaotic, or simply random. Along this line, a

lot of work has been done in the past two decades. A hit surprisingly, quite often the

exponential sensitivity to initial conditions (ESIC) in the rigorous mathematical sense

cannot he observed in experimental data. While current consensus is to attribute this

fact to the noise in the data, we shall show that the more general concept of power-law

sensitivity to initial conditions (PSIC) provides an interesting framework for time series

analysis. We illustrate this by studying the complicated dynamics of Internet transport


The Internet is one of the most complicated systems that man has ever made.

In recent years, two types of fascinating multi-scale behaviors have been found in the

Internet. One is the temporal-domain fractal and multifractal properties of network traffic

(see [76, 77] and references therein), which typically represents .I_ _regation of individual

host streams. The other is the spatial-domain scale-free topology of the Internet (see [78]

and references therein). Also, it has been observed that the failure of a single router may

trigger routing instability, which may be severe enough to instigate a route flap storm [79].

Furthermore, packets may be delivered out of order or even get dropped, and packet

reordering is not a pathological network behavior [80]. As the next generation Internet

applications such as remote instrument control and computational steering are being

developed, it is becoming increasingly important to understand the transport dynamics

in order to sustain the needed control channels over wide-area connections. Typically, the

transport dynamics over the Internet connections are a result of the nonlinear dynamics

of the widely-deploi- II Transmission Control Protocol (TCP) interacting with the Internet

traffic, which is stochastic and often self-similar. This leads one to naturally expect the

transport dynamics to be highly complicated. To tackle the hard problem of studying

the transport dynamics, we first systematically carry out network simulation using a

coninonly used network simulator, the ns-2 simulator (we have chosen the Tahoe version

of TCP, see http://www.isi. edu/nsnam/ns/), and then carefully study the dynamics of

the TCP congestion window-size traces collected over the real Internet connections.

4.1 Study of Transport Dynamics by Network Simulation

By carrying out network simulation, we critically examine whether TCP dynamics

can he chaotic, and if yes, to identify a route to chaos. We shall also present Internet

measurements to complement the simulation results. Since the advent of the Internet,

it has been perceived that the dynamics of complicated coninunication networks may

be chaotic [81]. If transport dynamics are indeed chaotic, then steady-state analysis and

study of convergence to a steady-state will not he the sole focus in practice, particularly

for remote control and computational steering tasks. Rather, one should also focus on the

Il II-! 1!1" non-convergent transport dynanxics. Although recently a lot of research has

been carried out [82-88] to better understand this issue, a lot of fundamental problems

remain to be answered. For example, the observation of chaos in [84] cannot he repeated,

which leads to the suspicion that the chaotic motions reported in [84] may correspond to

periodic motions with very long periods [89]. Fundamentally, we shall show, by developing

a 1-D discrete nmap, that the network scenarios with only a few competing TCP flows

studied in [84] cannot generate chaotic dynanxics. However, with other parameter settings,

chaos is possible.

In order to examine the temporal evolution of a dynamical systent, one often

measures a scalar time series at a fixed point in the state space. Takens embedding

theorem [34-36 ensures that the collective behavior of the dynamics described by the

dynamical variables coupled to the measured one can he conveniently studied by the

measured scalar time series. To illustrate the qualitative aspects of TCP dynamics by only

measuring a scalar time series, it is most effective to consider a single bottleneck link, as

such a link must be tightly coupled with the rest of the network. For simplicity, we shall

follow the setup of Veres and Boda [84] consisting of long lived TCP flows on a single link.

This setup has four parameters, the link speed C in Mbps, the link propagation delay d in

ms, the buffer size B in unit of packets of 1000 bytes, and the number of competing TCP

flows NV. We have found that NV acts as a critical bifurcation parameter. For example,

when C = 0.1, d = 10 ms, B = 10 and NV is small, we have observed only periodic

and quasi-periodic motions (Figs. 4-1(a-d)). When NV is large, such as NV = 17, chaos is

observed (see Fig. 4-2(a,b)). We thus conclude that one route to chaos in TCP dynamics

is the well-known quasi-periodic route. In order to show more diversified dynamical

behavior such as quasi-periodic motions with more than two incommensurate (i.e.,

independent) frequencies, below, however, we shall vary all the parameters. In particular,

we shall show that when NV is small, chaos cannot happen.

For each of the NV competing TCP flows, one can record its congestion window

size data. Let us only record one of them, and denote it by W(i). An example of

quasi-periodic W(i) is shown in Fig. 4-1(a), where C = 0.1 Mbps, delay d = 10

ms, B = 10, NV = 3. Its power spectral density (PSD) is shown in Fig. 4-1(b). We

observe many discrete sharp peaks. Note that when only around 105 points are used

to compute the PSD, one does not observe ]rn lw: peaks in the spectrum. Rather one

observes a broad spectrum, and hence, is tempted to interpret the W(i) time series

to be chaotic. Therefore, the W(i) time series is not ideally suited for the study of

(quasi-)periodic motions with long periods. This motivates us to define a new time

series, T(i), which is the time interval between the onset of two successive exponential

backoff (i.e., multiplicative decrease) of TCP, as indicated in Fig. 4-1(a). The start of an

exponential backoff indicates triggering of a loss episode or heavy congestion. Hence, the

T(i) is related to the well-known round trip time (RTT). In fact, it can be considered

to be more or less equivalent to the time interval between two successive loss bursts.

The T(i) time series for the W(i) of Fig. 4-1(a) is shown in Fig. 4-1(c), together with its



0 105


0 0.05 0.1

10'" (d)

0 1000 2000 3000 4000 5000
index i

index i




0 0.1 0.2 0.3 0.4 0.5

Figure 4-1.

Example of quasi-periodic congestion window size IT(i) time series sampled by
an equal time interval of 10 ms (a) and its power spectral density (PSD) (b).
The parameters used are C = 0.1 Mbps, delay d = 10 ms, B = 10 packets,
where a packet is of size 1000 bytes, and the number of TCP streams NV = :3.
(c) The T(i) time series corresponding to (a) and (d) its PSD. (e) Another
T(i) time series corresponding to C = 0.5 Mbps, d = 10 ms, B = 10, NV = :3,
and (f) its PSD.

PSD in Fig. 4-1(d). We observe that this T(i) time series is a simple periodic sequence,

and hence, the IT(i) time series is quasi-periodic. However, T(i) time series can still be

quasi-periodic. An example is shown in Fig. 4-1(e), with C

0.5 Mbps, delay d

ms, B = 10 packets, NV = :3, together with its PSD is Fig. 4-1(f). Note that there are still

many discrete sharp peaks in Fig. 4-1(f). To appreciate how many independent frequencies



0 50 100 150
index i


may be contained in the T(i) time series of Fig. 4-1(e), we have further extracted two new

time series, Til) (i), which is the interval time series between the successive local maxima

of the T(i) time series, and T(2) i), Which is the interval time series between the successive

local maxima of the Til) (i) time series. We have found that the T(2) i) iS Simply periodic.

Hence, we conclude that the W(i) time series is quasi-periodic with four independent


Before we proceed to discuss the chaotic TCP dynamics, we note that simple periodic

or quasi-periodic W(i) time series with less or more than four independent frequencies

can all be observed. In fact, we have found the "chaotic" congestion window size data

studied in [84] to be quasi-periodic with two independent frequencies. Furthermore, the

probability distributions for the constant, periodic, and quasi-periodic T(i) time series are

only composed of a few discrete peaks. This -II__- -; that the observed quasi-periodic T(i)
time series constitutes a discrete torus.

Let us now turn to the discussion of chaotic TCP dynamics. Shown in Figs. 4-2(a,c)

are two irregular T(i) time series. The parameters C, d, and B used for Fig. 4-2(a) are the

same as those for Figs. 4-1(a-d). Noting from Fig. 4-2 that T(i) often can be on the order

of 104 to 10s, hence, a not too long W(i) time series is even more ill-suited for the study of

the underlying irregular dynamics.

To show that the T(i) time series in Figs. 4-2(a,c) is indeed chaotic, we employ

the direct dynamical test for deterministic chaos developed by Gao and Zheng [21, 22].

The A(k) curves for the T(i) time series of Figs. 4-2(a,c) are shown in Figs. 4-2(b,d),

where the four curves, from bottom to top, correspond to shells of sizes (2-(i+1)/2, 2-i/2)

i = 14, 15, 16, 17. we have simply chosen m = 6 and L = 1. Very similar curves have

been obtained for other choices of m and L. We observe a very well defined linear common

envelope at the lower left corner of Figs. 4-2(b,d). The existence of a common envelope

guarantees that a robust positive Lyapunov exponent will be obtained no matter which

shell is used in the computation. Hence, the two T(i) time series are indeed chaotic.





8 5

2 15

xr 101





0 0
0 2000 4000 6000 0 10 20 30 40 50

Figure 4-2. Two examples of chaotic T(i) time series obtained with parameters (a)
C = 0.1, delay d = 10 ms, B = 10, and the number of TCP streams NV = 17
and (c) C = 0.1, d = 10 ms, B = 12, and NV = 19. The A(k) curves for (a) and
(c) are shown in (b) and (d), respectively.

We have remarked that the probability distributions for the T(i) time series of regular

motions are comprised of a few discrete peaks. What are the distributions for the chaotic

T(i) time series such as shown in Figs. 4-2(a,c)? They are power-law-like, P(T < t) ~ t Y,
for almost two orders of magnitude in t, as shown in Figs. 4-3(a,b). The exponent y is not

a universal quantity, however.

We now develop a simple 1-D map to describe the operation of TCP. TCP relies

on two mechanisms to set its transmission rate: flow control and congestion control.

10' 10-'

-10 a 10-

10-3~ 10-3

1-41 1-4
102 103 104 105 102 103 104 105
t t

Figure 4-3. Complementary cumulative distribution functions (CCDFs) for the two chaotic
T(i) time series of Figs. 4-2(a,c).

Flow control ensures that the sender sends no more than the size of the receiver's last

advertised flow-control window based on its available buffer size while congestion control

ensures that the sender does not unfairly overrun the network's available bandwidth.

TCP implements these mechanisms via a flow-control window ( fwnd) advertised by the

receiver to the sender and a congestion-control window (cwnd) adapted based on inferring

the state of the network. Specifically, TCP calculates an effective window (ewnd), where

ewnd = min( fwnd, cwnd), and then sends data at a rate ewnd/RTT, where RTT is the

round trip time of the connection. Recently Rao and Chua [85] developed an analytic

model describing TCP dynamics. By assuming fwnd being fixed and the TCP's slow start

only contributing to transient behavior, Rao and Chua have developed the w-up~date map

M~: [1, Wmax]i [1, Wmax],

I, + 1/l, if It, E1 RI, r tE R1

M~re)= < It < n /2"" if It, E R1, I tE R2
I /2"" if n-, 6R2, I, 1 1 6

I, + 1/l, if no Ea R I, 1I R1 1

In region R1, there is no packet loss, as = 0, and the TCP is in additive increase mode. In

R2, there are ni > 0 packet losses, where ni may be a complicated function of time, and

the TCP is in multiplicative decrease model.

The above map (Eq. 4-1) can be simplified [90],

I / 2"' if n -, > Wmax

Since fwnd is not assumed to be a constant, the modified map is more general than

the original Rao and Chua's map. To apply the above map to the study of a TCP

flow competing with NV 1 other TCP flows, one may lump the effect of the NV 1

other TCP flows as a background traffic. Hence, Wmax is determined by cwnd and the

background traffic, and typically is a complicated function of time, noticing that the

available bandwidth of a bottleneck link can vary with time considerably due to the

dynamic nature of background traffic. We have carried out simulation studies of Eq. 4-2

by assuming Wmax to be periodic and quasi-periodic, and have indeed observed chaos.

We are now ready to understand why chaos cannot occur when the number of

competing TCP flows is small, such as 2. To see this, let us assume that I, and no +

al, as well as I,. It and I, t + Al, t are all smaller than Wmax. It is then easy to

see that 1An' ,I t| < An- |. That is, small disturbance decays. This means during the

additive increase phase of TCP, nearby trajectories contract instead of diverge. Hence,

the dynamics are stable. In order for the dynamics to be unstable so that the Lyapunov

exponent is positive, the transition from the additive increase phase to the multiplicative

decrease phase has to be fast. This can be true only when the number of competing TCP

streams is large. We thus conclude that the network scenarios considered in [84] cannot

generate chaotic TCP dynamics.

How relevant is our simulation result to the dynamics of the real Internet? Are the

irregular T(i) time series shown in Figs. 4-2(a,c) an artifact of Tahoe version of TCP



0 50 100 150 200 250 300 350 100 101 102 103

Figure 4-4. The time series T(i) extracted from a W(i) time series collected using net100
instruments over ORNL-LSU connection. (a) the time series T(i), (b) the
complementary cumulative distribution function (CCDF) for the T(i) time

we used or more generally of the ns-2 simulator? To find an answer, we collected W(i)

measurements between ORNL and Louisiana State University at millisecond resolution

using the net100 instruments and computed T(i) as shown in Fig. 4-4(a), whose profile is

qualitatively quite similar to Fig. 4-2(a). In fact, it also follows a (truncated) power-law,

as shown in Fig. 4-4(b). The exponent for the power-law part is even smaller than the

time series of Fig. 4-2. The truncation in the power-law could be due to the shortness

of the T(i) time series. As we have pointed out, the T(i) time series is related to RTT.

Cottrell and Bullot have been experimenting with many advanced versions of TCP,

and also observed RTT time series very similar to these. We also have observed from

RTT and loss data measured on geographically dispersed paths on the Internet (with

a resolution of 1 sec) by researchers at the San Diego Supercomputer Center that the

probability distributions for the time interval between successive loss bursts typically

follow a power-law-like distribution, with the exponent y also smaller than 1. Thus,

we have good reason to believe that the observed irregular T(i) time series is not an

artifact of the ns-2 simulator, but reflects reality to some degree. In fact, the complicated

dynamics described here is much simpler than the actual dynamics of the Internet,

considering that there are all kinds of randomness in the Internet, especially that typical

TCP flows are not long lived [91], but vary considerably in durations determined by

specific applications. In the next section, we will analyze the actual dynamics of the

Internet measurements.

4.2 Complicated Dynamics of Internet Transport Protocols

As we have discussed earlier, it is very difficult to understand the dynamics of a

transport protocol over Internet connections. The 1! in r~ difficulties lie in two aspects.

First, it has been technologically hard to collect high quality measurements of network

transport variables on hli,-' Internet connections. Secondly, it is very difficult to

analytically handle the interaction between the deterministic and stochastic parts of

transport dynamics. The deterministic dynamics are due to the highly non-linear behavior

of TCP's Additive Increase and Multiplicative Decrease (AIMD) method emploi-x & for

congestion control. The stochastic component is due to the often self-similar traffic [92]

with which TCP competes for bandwidth and router buffers. Roughly speaking, the

interaction between the two is in terms of the additions to TCP congestion window-size,

denoted by cwnd, in response to acknowledgments, and multiplicative decreases in

response to inferred losses (ignoring the initial slow-start phase). Any protocol, however

simple its dynamics are, is generally expected to exhibit apparently complicated

dynamics due to its interaction with the stochastic network traffic. TCP in particular

is shown (albeit in simulation) to exhibit chaotic or chaos-like behavior even when the

competing traffic is much simpler such as a single competing TCP [84] or User Datagram

Protocol (UDP) stream [85]. The difficulty is to understand the dynamics when both the

phenomena are in effect, and equally importantly to understand their impact on actual

Internet streams. Historically, the methods from standard chaos and stochastic theories

have been unable to offer much new perspective on the transport dynamics.

In this section, we utilize the recently developed net100 instruments (obtained front

http://www.net100. org) to collect high quality TCP cwnd traces for actual Internet

connections. We analyze the measurements using the concepts of time dependent exponent

curves [21, 22, 37]. Our major purpose is to elucidate how the deterministic and stochastic

components of transport dynamics interact with each other on the Internet. It is of

considerable interest to note that recently there have been several important works on

TCP dynamics with the purpose of improving congestion control [93-96]. In fact, there

has been considerable effort in developing new versions and/or alternatives to TCP so

that the network dynamics can he more stable. Conventional v- 0-<~ Of analyzing the

network dynamics are unable to readily determine whether newer methods result in stable

transport dynamics or help in designing such methods. Our analysis shows two important

features in this direction:

(a) Randoniness is an integral part of the transport dynamics and must he explicitly
handled. In particular, the step sizes utilized for adjusting the congestion parameters
must he suitably varied (e.g. using Robbins-1\onro conditions) to ensure the eventual
convergence [97]. This is not the case for TCP which utilizes fixed step sizes.

(b) The chaotic dynamics of the transport protocols do have a significant impact on
practical transfers, and the protocol design must he cognizant of its effects. In
particular, it might he worthwhile to investigate protocols that do not contain
dominant chaotic regimes, particularly for remote control and steering applications.

The traditional transport protocols are not designed to explicitly address the above two

issues, but justifiably so since their original purpose is data transport rather than finer

control of dynanxics.

In the rest of this section, we shall first briefly describe the ownd data studied here.

Then we shall briefly explain the analysis procedure and describe typical results of the


We have collected a number of cwnd traces using single and two competing TCP

streams. This data was collected on two different connections front Oak Ridge National

Laboratory (ORNL) to Georgia Institute of Technology (GaTech) and to Louisiana State

University (LSU). The first connection has high-bandwidth (OC192 at 10Gbps) with

relatively low backbone traffic and a round-trip time of about 6 milliseconds. Four traces

were collected for this connection, two with a single TCP stream and the others with two

competing TCP streams. The second connection has much lower bandwidth (10 Mbps)

with higher levels of traffic and a round trip-time of about 26 milliseconds. The sampling

time is approximately 1 millisecond, with an error of 10s of microseconds (due to a "busy

waitingt measurement loop). Eight traces were collected for the second connection. The

results based on these eight traces are qualitatively very similar to the ORNL-GaTech

traces, and we only discuss the results for the latter. To appreciate better what these data

look like, we have plotted in Figs. 4-5(a-d) segments of four datasets for ORNL-GaTech

connection. Power spectral analysis of these data does not show any dominant peaks,

and hence, the dynamics are not simply oscillatory. Since our data was measured on the

Internet with Int -' background traffic, it is apparently more complicated and realistic

than the traces obtained by network simulation.

Next let us analyze cwnd data, x(i), i = 1, n, using the concepts of time

dependent exponent A(k) curves [21, 22, 37]. As has been discussed earlier, we first

need to employ the embedding theorem [34-36, 98] to construct vectors of the form:

K = [x(i), x(i + L), ..., x(i + (m 1)L)], where m is the embedding dimension and L the

delay time. The embedding theorem [34-36] states that when the embedding dimension

m is larger than twice the box counting dimension of the attractor, then the dynamics

of the original system can be studied from a single scalar time series. We note that the

embedding dimension used in [84] is only two, which has to be considered not large enough

(this may call for a closer examination of their conclusions).

Before we move on, we point out a few interesting features of the A(k) curves:

(i) For noise, only for k up to the embedding window size (m 1)L will A(k) increase.
Thus, whenever A(k) increases for a much larger range of k, it is an indication of
non-trivial deterministic structure in the data.

(ii) For periodic signals, A(k) is essentially zero for any k.

x 104


(b) 1 TCP source





2000 4000 6000

2000 4000 6000

x 104

3~ i

X 104
41(d) 2 TCP sources
[Hdi n8 n, Ellrfli rfl

TCP sources


2000 4000
sample n


2000 4000
sample n


Figure 4-5. Time series for the congestion window size cwnd for ORNL-GaTech

(iii) For quasi-periodic signals, A(k) is periodic with an amplitude typically smaller than
0.1, hence, for practical purposes, A(k) can be considered very close to 0.

(iv) When a chaotic signal is corrupted by noise, then the A(k) curves break themselves
away from the common envelope. The stronger the noise is, the more the A(k)
curves break away till the envelope is not defined at all. This feature has actually
been used to estimate the amount of both measurement and dynamic noise [99].

Now we are ready to compute and understand the A(k) curves for cwnd traces. The

set of A(k) curves, corresponding to Fig. 4-5 are plotted in Figs. 4-6. In the computations,

3 x 104 pOintS are used, and m

10, L = 1. The seven curves, from the bottom to

(a) 1 TCP source



0 50 100 150 200 250

(c 12 TC urces




0 50 100 150 200 250

50 100 150 200
Evolution time k

50 100 150 200
Evolution time k


Figure 4-6.

Time-dependent exponent A(k) vs. k curves for cwnd data corresponding to
Fig. 4-5. In the computations, 3 x 104 pOintS are used, and m 10, L 1 .
Curves from the bottom to top correspond to shells with sizes (2-(i+1)/2, 2-i/2)
i = 2, 3, ---,8.

top, correspond to shells of sizes (2-(i+1)/2, 2-i/2), i

2, 3, 8. We make the following


(i) The dynamics are complicated and cannot be described as either periodic or
quasi-periodic motions, since A(k) is much larger than 0.

(ii) The dynamics cannot be characterized as pure deterministic chaos, since in no case
can we observe a well-defined linear envelope. Thus the random component of the
dynamics due to competing network traffic is evident and can not simply be ignored.

(iii) The data is not simply noisy, since otherwise we should have observed that A(k) is
almost flat when k > (m 1)L. Thus, the deterministic component of dynamics
which is due to the transport protocol phI i- an integral role and must be carefully
studied. The features (ii) and (iii) indicate that the Internet transport dynamics
contains both chaotic and stochastic components.

(iv) There are considerable differences between the data with only 1 TCP source and
with 2 competing TCP sources. In the latter case, the A(k) curves sharply rise when
k just exceeds the embedding window size, (m 1)L. On the other hand, A(k)
for Fig. 4-6(a) with only one TCP source increases much slower when k just exceeds
(m-1)L. Also important is that the A(k) curves in Figs. 4-6(c,d) are much smoother
than those in Figs. 4-6(a,b). Hence, we can ;?i that the deterministic component of
the dynamics is more visible when there are more than 1 competing TCP sources
(along the lines of [84]).

Since the increasing part of the A(k) curves are not very linear, let us next examine if

the cwnd traces can be better characterized by the more generalized concept of power-law

sensitivity to initial conditions (PSIC). Figure 4-7 shows the A(k) vs. In k curves for

ORNL-GaTech connection. Very interestingly, we have now indeed observed a bunch of

better defined linear lines, especially for small scales. In particular, the transport dynamics

with more than 1 competing TCP sources are better described by the concept of PSIC.

To summarize, by analyzing a number of high quality congestion window-size data

measured on the Internet, we have found that the transport dynamics are best described

by the concept of PSIC, especially for small scales. It is interesting to further examine

how one might suppress the stochasticity of the network by executing more controls on

the network when making measurements, such as using a large number of competing TCP

sources together with a constant UDP flow. Our analysis motivates that both chaotic and

stochastic aspects be paid a close attention to in designing Internet protocols that are

required to provide the desired and tractable dynamics.

02 4

Time-dependent exponent A(k) vs. In k curves for cwnd data corresponding to
Fig. 4-5. In the computations, 3 x 104 pOintS are used, and m 10, L 1 .
Curves from the bottom to top correspond to shells with sizes (2-(i+1)/2, 2-i/2)
i = 2, 3, ---,8.

Figure 4-7.


Understanding the nature of sea clutter is crucial to the successful modeling of

sea clutter as well as to facilitate target detection within sea clutter. To this end, an

important question to ask is whether sea clutter is stochastic or deterministic. Since

the complicated sea clutter signals are functions of complex (sonietintes turbulent) wave

motions on the sea surface, while wave motions on the sea surface clearly have their

own dynamical features that are not readily described by simple statistical features,

it is very appealing to understand sea clutter hv considering some of their dynamical

features. In the past decade, Haykin et al. have carried out analysis of some sea clutter

data using chaos theory [100, 101], and concluded that sea clutter was generated hv an

underlying chaotic process. Recently, their conclusion has been questioned by a number of

researchers [102-108]. In particular, Unsworth et al. [105, 106] have demonstrated that the

two main invariants used by Haykin et al. [100, 101], namely the in! I::!!Itsinl likelihood of

the correlation dimension alI lin II. and the "false nearest neighbors" are problematic in

the analysis of measured sea clutter data, since both invariants may interpret stochastic

processes as chaos. They have also tried an improved method, which is based on the

correlation integral of Grassherger and Procaccia [109] and has been found effective in

distinguishing stochastic processes front chaos. Still, no evidence of deternxinisni or chaos

has been found in sea clutter data.

To reconcile ever growing evidence of stochasticity in sea clutter with their chaos

hypothesis, recently, Haykin et al. [110] have -II_t-r-- -1.. that the non-chaotic feature of sea

clutter could be due to many types of noise sources in the data. To test this possibility,

1\kDonald and Dantini [111] have tried a series of low-pass filters to remove noise; but

again they have failed to find any chaotic features. Furthermore, they have found that

the coninonly used chaotic invariant measures of correlation dimension and Lyapunov

exponent, computed by conventional r- .--4, produce similar results for measured sea

clutter returns and simulated stochastic processes, while a nonlinear predictor shows little

intprovenient over linear prediction.

While these recent studies highly ell_-- -r that sea clutter is unlikely to be truly

chaotic, a number of fundamental questions are still unknown. For example, most of the

studies in [102-106] are conducted by comparing measured sea clutter data with simulated

stochastic processes. We can ask: can the non-chaotic nature of sea clutter he directly

demonstrated without resorting to simulated stochastic processes'? Recognizing that

simple low-pass filtering does not correspond to any definite scales in phase space, can we

design a more effective method to separate scales in phase space and to test whether sea

clutter can he decomposed as signals plus noise'? Finally, will studies along this line he of

any help for target detection within sea clutter'?

In this chapter, we employ the direct dynamical test for deterministic chaos discussed

in Sec. 2.1 to analyze 280 sea clutter data measured under various sea and weather

conditions. The method offers a more stringent criterion for detecting low-dintensional

chaos, and can simultaneously monitor motions in phase space at different scales.

However, no chaotic feature is observed front any of these data at all the different scales

examined. But interestingly, we show that scale-dependent exponent corresponding to

large scale appears to be useful for distinguishing sea clutter data with and without

targets. This -II_--- -is that sea clutter may contain interesting dynamic features and that

the scale-dependent exponent may be an important parameter for target detection within

sea clutter. Furthermore, we find that sea clutter can he conveniently characterized by the

new concept of power-law sensitivity to initial conditions (PSIC). We show that for the

purpose of detecting targets within sea clutter, PSIC offers a more effective framework.

Below, we shall first briefly describe the sea clutter data. Then we study the

sea clutter data by employing the direct dynamical test for low-dintensional chaos

developed by Gao and Zheng [21, 22]. And we show that the scale-dependent exponent

corresponding to large scale can he used to effectively detect low observable targets within

sea clutter, while the conventional Lyapunov exponent fails in such a task. Finally, we

apply the new concept of PSIC to detect targets within sea clutter.

5.1 Sea Clutter Data

14 sea clutter datasets were obtained from a website maintained by Professor Simon

Haykin: http: //soma.ece.mcmaster. ca/ipix/dartmouth/datasets.html. The measurement

was made using the McMaster IPIX radar at Dartmouth, Nova Scotia, Canada. The

radar was mounted in a fixed position on land 25-30 m above sea level, with the operating

(carrier) frequency 9.39 GHz (and hence a wavelength of about 3 cm). It was operated at

low grazing angles, with the antenna dwelling in a fixed direction, illuminating a patch

of ocean surface. The measurements were performed with the wave height in the ocean

varying from 0.8 m to 3.8 m (with peak height up to 5.5 m), and the wind conditions

varying from still to 60 km/hr (with gusts up to 90 km/hr). For each measurement, 14

areas, called antenna footprints or range bins, were scanned. Their centers were depicted

as B1, B2, B14 in Fig. 5-1. The distance between two .Il11 Il:ent range bins was 15 m.

One or a few range bins (;?-,, Bi_l, Bi and Bi 1) hit a target, which was a spherical block

of styrofoam of diameter 1 m, wrapped with wire mesh. The locations of the three targets

were specified by their azimuthal angle and distance to the radar. They were (128 degree,

2660 m), (130 degree, 5525 m), and (170 degree, 2655 m), respectively. The range bin

where the target is strongest is labeled as the primary target bin. Due to drift of the

target, bins .Il11 Il:ent to the primary target bin may also hit the target. They are called

secondary target bins. For each range bin, there were 217 complex numbers, sampled with

a frequency of 1000 Hz. Amplitude data of two polarizations, HH (horizontal transmission,

horizontal reception) and VV (vertical transmission, vertical reception) were an~ llh-. 1

Fig. 5-2 shows two examples of the typical sea clutter amplitude data without and

with target. However, careful examination of the amplitude data indicates that 4 datasets

are severely affected by clipping. This can be readily observed from Figs. 5-3(a, b). We

40 60 80 100 120 140
Time (sec)

11 11_1 11* I)

h : Antenna height
S: Grazing angle
R. : Range (distance from the radar)
B1 ~ B14 : Range bins


BB..B. Bi B.+ ... B,
1 2 1-1 1 111


Figure 5-1. Collection of sea clutter data



o 15

E 10



E 1

0 20


Figure 5-2. Typical sea clutter amplitude data (a) without and (b) with target.

discard those 4 datasets, and >.1, lli. .. the remaining 10 measurements, which contain 280

sea clutter time series.

5.2 Non Chaotic Behavior of Sea Clutter

In this section, we employ the direct dynamical test for deterministic chaos developed

by Gao and Zheng [21, 22], to 2.1, lli. .. sea clutter data. The method offers a more

stringent criterion for low-dimensional chaos, and can simultaneously monitor motions

in phase space at different scales. The method has found numerous applications in the

0 20 40 60 80 100 120 140
Time (sec)

2.5 2.5

02 02

E 1.5 E 1.5

1 1

0.5 0.5
4.00 4.01 4.02 4.03 4.04 5.22 5.23 5.24 5.25 5.26
Time (sec) Time (sec)

Figure 5-:3. Two short segments of the amplitude sea clutter data severely affected hv

study of the effects of noise on dynamical systems [2:3, 24, 26, 27], estimation of the

strength of measurement noise in experimental data [28, 29], pathological tremors [:30],

shear-thickening surfactant solutions [:31], dilute sheared aqueous solutions [:32], and

serrated plastic flows [:33]. In particular, this method was used by Gao et al. [107, 108] to

analyze one single set of sea clutter data. While chaos was not observed front that dataset,

no definite general conclusion could be reached, due to lack of large amount of data at

that time. Here we systematically study 280 sea clutter data measured under various sea

and weather conditions, and examine whether any chaotic features can he found front

these sea clutter data.

The explicit incorporation of scales in the Gao and Zheng's test [21, 22] enables us to

simultaneously monitor motions in phase space at different scales. We have systematically

analyzed 280 amplitude sea clutter time series measured under various sea and weather

conditions. However, no chaotic feature has been observed front any of these data at

all the scales examined. Typical examples of the A(k) vs. k curves for the sea clutter

amplitude data without targets are shown in Figs. 5-4(a, c, e) and the curves for the data

with targets shown in Figs. 5-4(h, d, f), respectively. We have simply chosen m = 6 and











5 10 15 20 25 30

5 10 15 20 25 30

5 10 15 20 25 30

5 10 15 20 25 30

Oi(e '" (f)c
0 5 10 15 20 25 30 0 5 10 15 20 25 30
k k

Figure 5-4. Examples of the time-dependent exponent A(k) vs. k curves for the sea clutter
data (a, c, e) without and (b, d, f) with the target. Six curves, from bottom to
top, correspond to shells (2-(i+1)/2, 2-i/2) With i 13, 14, 15, 16, 17, and 18.
The sampling time for the sea clutter data is 1 msee, and embedding
parameters are m 6, L 1 217 data points are used in the computation.

L = 1. Very similar curves have been obtained for other choices of m and L. We have

not observed a common envelope at any scales. In fact, the results of Fig. 5-4 are generic

among all the 280 sea clutter data on~ li. .1 here. Hence, we have to conclude that none of

the sea clutter data is chaotic.

5.3 Target Detection within Sea Clutter by Separating Scales

Robust detection of low observable targets within sea clutter is a very important issue

in remote sensing and radar signal processing applications, for a number of reasons: (i)

identifying objects within sea clutter such as submarine periscopes, low-flying aircraft,

and missiles can greatly improve our coastal and national security; (ii) identifying small

marine vessels, navigation buoys, small pieces of ice, patches of spilled oil, etc. can

significantly reduce the threat to the safety of ship navigation; (iii) monitoring and

policing of illegal fishing is an important activity in the environmental monitoring. Due

to the rough sea surface, which is a consequence of energy transfer from small scale to

large scale of sea surface wave, and the multipath propagation of the radar backscatter,

sea clutter is often highly non-Gaussian, even spiky, especially in heavy sea conditions.

Hence, sea clutter modelling is a very difficult problem, and a lot of effort has been made

to study sea clutter, both through the analysis of the distributions of sea clutter, including

Weibull [112], log-normal [113], K( [114, 115], and compound-Gaussian distributions [116],

as well as using chaos theory [100-104, 107, 108, 110, 111], wavelets [117], neural

networks [118, 119], wavelet-neural net combined approaches [120, 121], and the concept of

fractal dimension [122] and fractal error [123, 124]. However, no simple method has been

found to robustly detect low observable objects within sea clutter [125].

In this section, we examine whether the Lyapunov exponent estimated by conventional

methods and the scale-dependent exponent may be helpful for target detection within sea


We first check whether the Lyapunov exponent A estimated by conventional methods

can be used for distinguishing sea clutter data with and without targets. As pointed

out earlier, the conventional way of estimating the Lyapunov exponent is to compute

A(k)/kbt, where A(k) is defined as in Eq. 2-4, subject to the conditions that ||Xi-Xj|| < r

and ||Xi+k Xj,, || < R, where r and R are prescribed small distance scales. The condition

||Xi Xj || < r amounts to our smallest shell in computing the A(k) curves. The condition

||Xi+k Xj+k I < R is presumably to set the time scale k smaller than the prediction

time scale k. Our analysis finds that the Lyapunov exponent estimated this way fails to

detect targets from the sea clutter data on~ lli. here. This motivates us to try a modified

approach, as described below.

The modified approach for estimating the Lyapunov exponent A uses the least-squares

fit to the first few samples of the A(k) curve where the A(k) curve increases linearly or

quasi-linearly. This is done for all the 280 sea clutter time series. To better appreciate the

variation of the Lyapunov exponent A among the 14 range hins of the sea clutter data,

we have first subtracted the A parameter of each hin by the minimum of A values for that

measurement, and then normalized the obtained A values by its maximum. The variations

of the A parameters with the 14 range hins for the 10 HH measurements are shown in

Figs. 5-5(a)-(j), respectively, where open circles denote the range hins with the target,

and asteroids denote the hins without the target. The primary target hin is explicitly

indicated hv an arrow. Similar results have been obtained for the 10 VV measurements.

While Figs. 5-5(a) and (h) indicate that the primary target hin can he separated from bins

without the target, in general, we have to conclude that the Lyapunov exponent estimated

by methods equivalent or similar to conventional means cannot he used for distinguishing

sea clutter data with and without targets.

Next let us examine whether the scale-dependent exponent corresponding to large

scale may be useful for detecting targets within sea clutter data. By scale-dependent

exponent, we mean that if we use the least-squares fit to the linearly or quasi-linearly

increasing part of the A(k) curves of different shells, the slopes of those lines depend

on which shells are used for computing the A(k) curves. In other words, if we plot the

exponent estimated this way against the size of the shell, then the exponent varies

with the shell size or the scale. The scale-dependent exponent has been used in the

study of noise-induced chaos in an optically injected semiconductor laser model to

examine how noise affects different scales of dynamic systems [26]. We now focus on the

exponent corresponding to large shell or scale. Figures 5-6(a)-(j) show the variations of

the scale-dependent exponent corresponding to large scale with the 14 range hins for the

JTarget Bin

&z 05

0 5 10 15



0 5 10 15


0 5 10 15

(c) J

0 5 10 15

(e) +

0 5 10 15

5 10 15

00 5 10 15 00 5v 10 `15
Bln Number Bin Number

Figure 5-5. Variations of the Lyapunov exponent A estimated by conventional methods vs.
the 14 range bins for the 10 HHI- measurements. Open circles denote the range
bins with target, while denote the bins without target. The primary target
bin is explicitly indicated by an arrow.

same 10 HH measurements as studied in Fig. 5-5. We observe that the primary target bin

can be easily separated from the range bins without the target, since the scale-dependent

exponent for the primary target bin is much larger than those for the bins without the


It is thus clear that the scale-dependent exponent corresponding to large scale is

very useful for distinguishing sea clutter data with and without targets. This -II__- -;

that sea clutter may contain interesting dynamic features and that the scale-dependent

exponent may be an important parameter for target detection within sea clutter. This

(a) Primar (b)
Target Bin

0505 /q

0g 5 10 15 0 5 10 15

(c) (d)

4z 05 0

0 5 10 15 0 5 10 15

(e) ? (f) I

A 05 /05

0 5 10 15 0O 5 10 15

(g) / (h)

Az 05 05

0 0
0 5 10 15 0 5 10 15

4z 05 / 05

O' 0
0 5 10 15 0 5 10 15
Bln Number Bin Number

Figure 5-6. Variations of the scale-dependent exponent corresponding to large scale vs. the
range bins for the 10 HHI- measurements. Open circles denote the range bins
with target, while denote the bins without target. The primary target bin is
explicitly indicated by an arrow.

example clearly signifies the importance of incorporating the concept of scale in a measure.

In fact, the concept of scale is only incorporated in the time dependent exponent A(k)

curves [21, 22] in a static manner. In ('! .pter 6, we shall see that when a measure

dynamically incorporates the concept of scale, it will become much more powerful.

Note that one difficulty of using the scale-dependent exponent for target detection

within sea clutter is that we may need to choose a suitable scale for estimating the

exponent. This makes the method not easy to use, and it is hard to make the method


5.4 Target Detection within Sea Clutter by Power-law Sensitivity to Initial
Conditions (PSIC)

In this section, we show how the concept of PSIC can he applied to detect target

within sea clutter data. Examples of the A(k) vs. In k curves for the range hins without

target and those with target of one measurement are shown in Figs. 5-7(a) and (b),

respectively, where the curves denoted by asteroids are for data without the target,

while the curves denoted by open circles are for data with the target. We observe that

the curves are fairly linear for the first a few samples. Also notice that the slopes of

the curves for the data with the target are much larger than those for the data without

the target. For convenience, we denote the slope of the curve by the parameter /9. To

better appreciate the variation of the /9 parameter with the range hims, we normalize /9

of each hin by the nmaxinial /9 value of the 14 range hins within the single measurement.

Fig. 5-8 shows the variation of the /9 parameter with the 14 range hims. It is clear that

the /9 parameter can he used to distinguish sea clutter data with and without target.

Interestingly, the feature shown in Fig. 5-8 is generically true for all the measurements. It

is worth pointing out that usually the A(k) vs. In k curves for different scales are almost

parallel, thus the estimated /9 parameter is relatively less dependent on the scale. This is a

very nice feature of this method.

Let us examine if a robust detector for detecting targets within sea clutter can he

developed based on the /9 parameter. We have systematically studied 280 time series

of the sea clutter data measured under various sea and weather conditions. To better

appreciate the detection performance, we have first only focused on hins with primary

targets, but omitted those with secondary targets, since sometimes it is hard to determine

whether a hin with secondary target really hits a target or not. After omitting the range

hin data with secondary targets, the frequencies for the /9 parameter under the two

hypotheses (the hins without targets and those with primary targets) for HH datasets

0 0.5 1 1.5 2 2.5 3 3.5

In k

In k

Figure 5-7.

Examples of the A(k) vs. In k curves for range bins (a) without and (b) with
target. Open circles denote the range bins with target, while denote the bins
without target.

Target Bin


co. 0.85



0 2 4

Bin Number

10 12 14

Figure 5-8.

Variation of the parameter with the 14 range bins. Open circles denote the
range bins with target, while denote the bins without target.

are shown in Fig. 5-9. We observe that the frequencies completely separate for the HH

datasets. This means the detection accuracy can be 1011' .

no targets
12 primary targets




n mm

t0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 5-9. Frequencies of the bins without
the HH datasets.

targets and the bins with primary targets for


In C'!. Ilter 1, we have emphasized the importance of developing scale-dependent

measures to simultaneously characterize behaviors of complex multiscaled signals on a

wide range of scales. In this chapter, we shall develop an effective algorithm to compute

an excellent multiscale measure, the scale-dependent Lyapunov exponent (SDLE), and

show that the SDLE can readily classify various types of complex motions, including truly

low-dimensional chaos, noisy chaos, noise-induced chaos, processes defined by power-law

sensitivity to initial conditions (PSIC), and complex motions with chaotic behavior on

small scales but diffusive behavior on large scales. Finally, we shall discuss how the SDLE

can help detect hidden frequencies in large scale orderly motions.

6.1 Basic Theory

The scale-dependent Lyapunov exponent (SDLE) is a variant of FSLE [126]. The

algorithm for calculating the FSLE is very similar to the Wolf et al's algorithm [127].

It computes the average r-fold time by monitoring the divergence between a reference

trajectory and a perturbed trajectory. To do so, it needs to define 1.. neighbors",

as well as needs to perform, from time to time, a renormalization when the distance

between the reference and the perturbed trajectory becomes too large. Such a procedure

requires very long time series, and therefore, is not practical. To facilitate derivation of

a fast algorithm that works on short data, as well as to ease discussion of continuous but

non-differentiable stochastic processes, we cast the definition of the SDLE as follows.

Consider an ensemble of trajectories. Denote the initial separation between two

nearby trajectories by co, and their averagee separation at time t and t + at by et and et as,

respectively. Being defined in an average sense, et and et as can he readily computed from

any processes, even if they are non-differentiable. Next we examine the relation between et

and et+nt, where at is small. When at 0 we have,

et+nt = etex(tt)nt (6-1)

where A(et) is the SDLE. It is given by

In et+nt In et
A(et) = (6-2)

Given a time series data, the smallest at possible is the sampling time 7r.

The definition of the SDLE so_~- -;- Ma simple ensemble average based scheme to

compute it. A straightforward way would be to find all the pairs of vectors in the phase

space with their distance approximately e, and then calculate their average distance after a

time at. The first half of this description amounts to introducing a shell (indexed as k),

Ek | Vs Vy | < k Ek(6-3)

where 1%, Vj are reconstructed vectors, Ek, (the radius of the shell) and A~rk (the width

of the shell) are arbitrarily chosen small distances. Such a shell may be considered as

a differential element that would facilitate computation of conditional probability. To

expedite computation, it is advantageous to introduce a sequence of shells, k = 1, 2, 3, -

Note that this computational procedure is similar to that for computing the so-called

time-dependent exponent (TDE) curves [21, 22, 37].

With all these shells, we can then monitor the evolution of all of the pairs of vectors

(1M, Vj) within a shell and take average. When each shell is very thin, by assuming that

the order of averaging and taking logarithm in Eq. 6-2 can be interchanged, we have

where t and at are integers in unit of the sampling time, and the angle brackets denote

average within a shell. Note that contributions to the SDLE at a specific scale from

different shells can be combined, with the weight for each shell being determined by the

number of the pairs of vectors (1M, Vyj) in that shell. In the following, to see better how

each shell characterizes the dynamics of the data on different scales, we shall plot the A(e)

curves for different shells separately.

In the above formulation, it is implicitly assumed that the initial separation, ||% 1- Vy ||l,

aligns with the most unstable direction instantly. For high-dimensional systems, this is

not true, especially when the growth rate is non-uniform and/or the eigenvectors of the

Jacobian are non-normal. Fortunately, the problem is not as serious as one might be

concerned, since our shells are not infinitesimal. When computing the TDE curves [21,

22, 37], we have found that when difficulties arise, it is often sufficient to introduce an

additional condition,

I i> m 1 L6-5

when finding pairs of vectors within each shell. Such a scheme also works well when

computing the SDLE. This means that, after taking a time comparable to the embedding

window (m 1)L, it would be safe to assume that the initial separation has evolved to the

most unstable direction of the motion.

Before proceeding on, we wish to emphasize the 1!! ri ~ difference between our

algorithm and the standard method for calculating the FSLE. As we have pointed out,

to compute the FSLE, two trajectories, one as reference, another as perturbed, have to

be defined. This requires huge amounts of data. Our algorithm avoids this by employing

two critical operations to fully utilize information about the time evolution of the data: (i)

The reference and the perturbed trajectories are replaced by time evolution of all pairs of

vectors satisfying the inequality (6-5) and falling within a shell, and (ii) introduction of

a sequence of shells ensures that the number of pairs of vectors within the shells is large

while the ensemble average within each shell is well defined. Let the number of points

needed to compute the FSLE by standard methods be NV. These two operations imply

that the method described here only needs about z/Vpoints to compute the SDLE. In the

following, we shall illustrate the effectiveness of our algorithm by examining various types

of complex motions.

6.2 Classification of Complex Motions

T> understand the SDLE as well as appreciate its power, we apply it to classify

various types of complex motions.

6.2.1 Chaos, Noisy Chaos, and Noise-induced Chaos

Obviously, for truly low-dimensional chaos, A(e) equals the largest positive Lyapunov

exponent, and hence, must he independent of e over a wide range of scales. For noisy

chaos, we expect A(e) to depend on small e. T> illustrate both features, we consider the

chaotic Lorenz system with stochastic forcing described by Eq. 2-5. Figure 6-1(a) shows

five curves, for the cases of D = 0, 1, 2, 3, 4. The computations are done with 10000 points

and m = 4, L = 2. We observe a few interesting features:

For the clean chaotic signal, A(e) slightly fluctuates around a constant (which
numerically equals the largest positive Lyapunov exponent) when e is smaller than a
threshold value which is determined by the size of the chaotic attractor. The reason
for the small fluctuations in A(e) is that the divergence rate varies from one region of
the attractor to another.

When there is stochastic forcing, A(e) is no longer a constant when e is small, but
increases as -y In e when the scale e is decreased. The coefficient ]* does not seem to
depend on the strength of the noise. This feature -II__- -is that entropy generation is
infinite when the scale e approaches to zero. Note that the relation of A(e) ~ -y In e
has also been observed for the FSLE and the e-entropy. In fact, such a relation can
he readily proven for the e-entropy.

When the noise is increased, the part of the curve with A(e) ~ -y 1n e shifts to the
right. In fact, little chaotic signature can he identified when D is increased beyond
:3. When noise is not too -1i~~e.: this feature can he readily used to quantify the
strength of noise.

Next we consider noise-induced chaos. To illustrate the idea, we follow [2:3] and study

the noisy logistic map

r,w+1 = p~r, (1 r,w) + P,, O < r,z < 1 (6-6)

2.5 0.7
(a) LorenZ0. (b) Logistic

2 -05- (1): slope =0.28

S0.4 h (2): slope = 0

D. = o0.
D = 3 0.1

0.5 -v 0
-1.5 -1 -0.5 -2 -1.5 -1 -0.5 0
log ,E log ,E

Figure 6-1. Scale-dependent Lyapunov exponent A(e) curves for (a) the clean and the
noisy Lorenz system, and (b) the noise-induced chaos in the logistic map.
Curves from different shells are designated by different symbols.

where p is the bifurcation parameter, and P, is a Gaussian random variable with

zero mean and standard deviation o-. In [23], we reported that at p = 3.74 and

o- = 0.002, noise-induced chaos occurs, and thought that it may be difficult to distinguish

noise-induced chaos from clean chaos. In Fig. 6-1(b), we have plotted the A(et) for this

particular noise-induced chaos. The computation was done with m = 4, L = 1 and 10000

points. We observe that Fig. 6-1(b) is very similar to the curves of noisy chaos plotted in

Fig. 6-1(a). Hence, noise-induced chaos is similar to noisy chaos, but different from clean

chaos .

At this point, it is worth making two comments: (i) On very small scales, the effect

of measurement noise is similar to that of dynamic noise. (ii) The A(e) curves shown in

Fig. 6-1 are based on a fairly small shell. The curves computed based on larger shells

collapse on the right part of the curves shown in Fig. 6-1. Because of this, for chaotic

systems, one or a few small shells would be sufficient. If one wishes to know the behavior

of A on ever smaller scales, one has to use longer and longer time series.

1x 10-


10- 100

Figure 6-2. Scale-dependent Lyapunov exponent A(e) curve for the Mackey-Glass system.
The computation was done with ni = 5, L = 1, and 5000 points sampled with a
time interval of 6.

Finally, we consider the Mackey-Glass delay differential system [39], defined by

Eq. 2-6. The system has two positive Lyapunov exponents, with the largest Lyapunov

exponent close to 0.007 [21]. Having two positive Lyapunov exponents while the value of

the largest Lyapunov exponent of the system is not much greater than 0, one might he

concerned that it may be difficult to compute the SDLE of the system. This is not the

case. In fact, this system can he analyzed as straightforwardly as other dynamical systems

including the Henon map and the Rossler system. An example of the A(e) curve is shown

in Fig. 6-2, where we have followed [21] and used ni = 5, L = 1, and 5000 points sampled

with a time interval of 6. Clearly, we observe a well defined plateau, with its value close to

0.007. This example illustrates that when computing the SDLE, one does not need to be

very concerned about non-uniform growth rate in high-dimensional systems.

Up till now, we have focused on the positive portion of A(es). It turns out that when t

is large, A(e) becomes oscillatory, with mean about 0. Denote the corresponding scales by

ea, and call them the limiting or characteristic scales. They are the stationary portion of

et, and hence, they may still be a function of time. We have found that the limiting scales

capture the structured component of the data. This feature will be further illustrated

later, when we discuss identification of hidden frequencies. Therefore,the positive portion

of A(et) and the concept of limiting scale provide a comprehensive characterization of the

signals .

6.2.2 Processes Defined by Power-law Sensitivity to Initial Conditions (PSIC)

In C'!s Ilter 3, we have seen that power-law sensitivity to initial conditions (PSIC)

provides a common foundation for chaos theory and random fractal theory. Here, we

examine whether the SDLE is able to characterize processes defined by PSIC. Pleasingly,

this is indeed so. Below, we derive a simple equation relating the A, and q of PSIC to the


First, we recall that PSIC is defined by

( = lim = 1 ( ) )

Since et = ||av(0)||s,, it is now more convenient to express the SDLE as a function of (c.

Using Eq. 6-2, we find that
(()=In $t+nt In (a(67

When at 0, 1+ (1 q)Agt > (1 q)XAft. Simplifying Eq. 6-7, we obtain

A((s) = Aq tq-1 (6-8)

We now consider three cases:

*For chaotic motions, q = 1, therefore,

A((s) = Xq = COnStant

F'or 1/ fB noise, we, hae = 1_ (Eq. :328) and A," = H (Eq. :329). TIherefor~el

A((s) = Hg, H

For Le~vy fightsl weit have q = 1- a? (Eq. :3-38) and A,"l = 1/0~ (Eq. :3-39). Therefore,

6.2.3 Complex Motions with Multiple Scaling Behaviors

Some dynamical systems may exhibit multiple scaling behaviors, such as chaotic

behavior on small scales but diffusive behavior on large scales. To see how the SDLE can

characterize such systems, we follow Cencini et al. [128] and study the following map,

Xrn+l = [Xrn] + F(xrn [Xrnl) + aife (6-9)

where [r,z] denotes the integer part of r,t, r7< is a noise uniformly distributed in the interval

[-1, 1], o- is a parameter quantifying the strength of noise, and Fly) is given by

Fly) = (2 + A) y if y E [0, 1/2) (6-10)
(2 + A) y (1 + a) if y E (1/2, 1]

The map Fly) is shown in Fig. 6-3 as the dashed lines. It gives a chaotic dynamics with

a positive Lyapunov exponent In(2 + a). On the other hand, the term [:re] introduces a

random walk on integer grids.

It turns out this system is very easy to analyze. When a = 0.4, with only 5000 points

and m = 2, L = 1, we can resolve both the chaotic behavior on very small scales, and the

normal diffusive behavior (with slope -2) on large scales. See Fig. 6-4(a).

We now ask a question: Given a small dataset, which type of behavior, the chaotic

or the diffusive, is resolved first? To answer this, we have tried to compute the SDLE

with only 500 points. The result is shown in Fig. 6-4(b). It is interesting to observe that

the chaotic behavior can he well resolved by only a few hundred points. However, the





Figure 6-3.

The function F(x) (Eq. 6-10) for a = 0.4 is shown as the dashed lines. The
function G(x) (Eq. 6-11) is an approximation of F(x), obtained using 40
intervals of slope 0. In the case of noise-induced chaos discussed in the paper,
G(x) is obtained from F(x) using 104 HinerValS Of Slope 0.9.





clean data
noisy data

10-4 10-3 10-2 10 1

10-4 10-3 10-2

Scale-dependent Lyapunov exponent A(e) for
(a) 5000 points were used, for the noisy case,

Figure 6-4.

the model described by Eq. 6-9.
o- = 0.001. (b) 500 points were

diffusive behavior needs much more data to resolve. Intuitively, this makes sense, since

the diffusive behavior amounts to a Brownian motion on the integer grids and is of much

higher dimension than the small scale chaotic behavior. Therefore, more data are needed

to resolve it.

We have also studied the noisy map. The resulting SDLE for a = 0.001 is shown in

Fig. 6-4(a), as squares. We have again used 5000 points. While the behavior of the SDLE

-II- -_ -r ;noisy dynamics, with 5000 points, we are not able to well resolve the relation

A(e) ~ In e. This indicates that for the noisy map, on very small scales, the dimension is

very high.

Map (6-9) can be modified to give rise to an interesting system with noise-induced

chaos. This can be done by replacing the function Fly) in map (6-9) by G(y) to obtain

the following map [128],

xt+l = [Xt] + G(xt [xt]) + arlt (6-11)

where Orl is a noise uniformly distributed in the interval [-1, 1], a is a parameter

quantifying the strength of noise, and G(y) is a piecewise linear function which approximates

Fly) of Eq. 6-10. An example of G(y) is shown in Fig. 6-3. In our numerical simulations,

we have followed Cencini et al. [128] and used 104 interValS Of Slope 0.9 to obtain G(y).

With such a choice of G(y), in the absence of noise, the time evolution described by

the map (6-11) is nonchaotic, since the largest Lyapunov exponent In(0.9) is negative.

With appropriate noise level (e.g., a = 10-4 or 10-3), the SDLE for the system becomes

indistinguishable to the noisy SDLE shown in Fig. 6-4 for the map (6-9). Having a

diffusive regime on large scales, this is a more complicated noise-induced chaos than the

one we have found from the logistic map.

Before proceeding on, we make a comment on the computation of the SDLE from

deterministic systems with negative largest Lyapunov exponents, such as the map (6-11)

without noise. A transient-free time series from such systems is a constant time series.

Therefore, there is no need to compute the SDLE or other metrics. When the time series

contains transients, if the time series is sampled with high enough sampling frequency,

then the SDLE is negative. In the case of simple exponential decay to a fixed point, such

as expressible as e-xt, where A > 0, one can readily prove that the SDLE is -A. Since such

systems are not complex, we shall not be further concerned about them.

6.3 Characterizing Hidden Frequencies

Defining and characterizing large scale orderly motions is a significant issue in many

disciplines of science. One of the most important types of large scale orderly motions is

the oscillatory motions. An interesting type of oscillatory motions is associated with the

so-called hidden frequency phenomenon. That is, when the dynamics of a complicated

system is monitored through the temporal evolution of a variable x, Fourier analysis of

x(t) may not -II---- -1 any oscillatory motions. However, if the dynamics of the system is

monitored through the evolution of another variable, 11-, z, then the Fourier transform of

z(t) may contain a well-defined spectral peak, indicating oscillatory motions. An example

is the chaotic Lorenz system described by Eq. 2-5. In Figs. 6-5 (a-c), we have shown the

power-spectral density (PSD) of the x, y, z components of the system. We observe that the

PSD of x(t) and y(t) are simply broad. However, the PSD of z(t) shows up a very sharp

spectral peak. Recall that geometrically, the Lorenz attractor consists of two scrolls (see

Fig. 6-6). The sharp spectral peak in the PSD of z(t) of the Lorenz system is due to the

circular motions along either of the scrolls.

The above example illustrates that if the dynamics of a system contains a hidden

frequency that cannot be revealed by the Fourier transform of a measured variable

( 11-, x(t)), then in order to reveal the hidden frequency, one has to embed x(t) to a

suitable phase space. This idea has led to the development of two interesting methods

for identifying hidden frequencies. One method is proposed by Ortega [129, 130], by

computing the temporal evolution of density measures in the reconstructed phase space.

Another is proposed by C'I. I i. et al. [131], by taking singular value decomposition of local

neighbors. The Ortega's method has been applied to an experimental time series recorded

10 "


102 10-4~
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
10 Frequency f 14Frequency f
1(b) (e)


1021 10-4
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
109Frequency f 14Frequency f

106 100

1031 1 10-4
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
Frequency f Frequency f

F'igur~e 6-5. P~ower-spectral density (PSD>) for (a) x(t), (b) y(t), (c) z(t), (d) ei)(t), (e)

from a far-infrared laser in a chaotic state. The laser dataset can be downloaded from

the link http://www-psy ch. stanf ord .edu/~ andreas/Time- Ser ies /Sant aFe .htm1. It

contains 10000 points, sampled with a time interval of 80 ns. Fig. 6-7(a) shows the laser

dataset. The PSD of the data is shown in Fig. 6-7(b), where one observes a sharp peak

around 1.7 MHz. Fig. 6-7(c) shows the PSD of the density time series, where one notes

an additional spectral peak around 37 kHz. This peak is due to the envelope of chaotic

pulsations, which is discernable from Fig. 6-7(a).

Full Text








Firstofall,Iwouldliketothankmyadvisor,ProfessorJianboGao,forhisinvaluableguidanceandsupport.Iwasluckytoworkwithsomeonewithsomanyoriginalideasandsuchasharpmind.Withouthishelp,myPh.D.degreeandthisdissertationwouldhavebeenimpossible.DuringmyyearsattheUniversityofFlorida,Iwasmovedbyhismagnicentpersonality.AfamousChineseadagesays,\Onedayasmyadvisor,entirelifeasmyfather."ProfessorGaohighlyqualies.Ialsowanttothankmycommitteemembers(ProfessorsSu-ShingChen,YuguangFang,JoseC.Principe,andKeithD.White),fortheirinterestinmyworkandvaluablesuggestions.Mythanksofcoursealsogotoallthefaculty,staandmyfellowstudentsintheDepartmentofElectricalandComputerEngineering.IextendparticularthankstomyocematesJaeminLee,YiZheng,andUngsikKim,fortheirhelpfuldiscussions.Iexpressmyappreciationtoallofmyfriendswhooeredencouragement.IspeciallythankQianZhan,XingyanFan,LingyuGuandhiswife(JingZhang),YanminXu,YanLiu,andYanaLiang.TheygavemealotofhelpduringmyearlyyearsintheUniversityofFlorida.Theytaughtmetocookandmademylifemorecolorful.Finally,Ioweagreatdebtofthankstomyparents,grandparents,andmyboyfriend,JingAi,fortheirpillar-likesupportandtheirunwaveringbeliefandcondenceinme|withoutthisIdonotthinkIwouldhavemadeitthisfar! 4


page ACKNOWLEDGMENTS ................................. 4 LISTOFFIGURES .................................... 7 ABSTRACT ........................................ 9 CHAPTER 1INTRODUCTION .................................. 10 1.1OverviewofDynamicalSystems ........................ 11 1.2ExamplesofMultiscalePhenomena ...................... 15 1.3BriefIntroductiontoChaos .......................... 20 1.4BriefIntroductiontoFractalGeometry .................... 21 1.5ImportanceofConnectingChaosandRandomFractalTheories ...... 23 1.6ImportanceoftheConceptofScale ...................... 24 1.7StructureofthisDissertation ......................... 24 2GENERALIZATIONOFCHAOS:POWER-LAWSENSITIVITYTOINITIALCONDITIONS(PSIC) ................................ 27 2.1DynamicalTestforChaos ........................... 27 2.2GeneralComputationalFrameworkforPower-lawSensitivitytoInitialConditions(PSIC) ............................... 33 2.3CharacterizingEdgeofChaosbyPower-lawSensitivityofInitialConditions(PSIC) ...................................... 36 3CHARACTERIZINGRANDOMFRACTALSBYPOWER-LAWSENSITIVITYTOINITIALCONDITIONS(PSIC) ........................ 41 3.1Characterizing1=fProcessesbyPower-lawSensitivitytoInitialConditions(PSIC) ...................................... 41 3.1.1CharacterizingFractionalBrownianMotion(fBm)ProcessesbyPower-lawSensitivitytoInitialConditions(PSIC) .......... 43 3.1.2CharacterizingON/OFFProcesseswithParetodistributedONandOFFPeriodsbyPower-lawSensitivitytoInitialConditions(PSIC) 50 3.2CharacterizingLevyProcessesbyPower-lawSensitivitytoInitialConditions(PSIC) ...................................... 54 4STUDYOFINTERNETDYNAMICSBYPOWER-LAWSENSITIVITYTOINITIALCONDITIONS(PSIC) ........................... 60 4.1StudyofTransportDynamicsbyNetworkSimulation ............ 61 4.2ComplicatedDynamicsofInternetTransportProtocols ........... 69 5


............... 76 5.1SeaClutterData ................................ 78 5.2NonChaoticBehaviorofSeaClutter ..................... 79 5.3TargetDetectionwithinSeaClutterbySeparatingScales .......... 81 5.4TargetDetectionwithinSeaClutterbyPower-lawSensitivitytoInitialConditions(PSIC) ............................... 86 6MULTISCALEANALYSISBYSCALE-DEPENDENTLYAPUNOVEXPONENT(SDLE) ........................................ 89 6.1BasicTheory .................................. 89 6.2ClassicationofComplexMotions ....................... 92 6.2.1Chaos,NoisyChaos,andNoise-inducedChaos ............ 92 6.2.2ProcessesDenedbyPower-lawSensitivitytoInitialConditions(PSIC) .................................. 95 6.2.3ComplexMotionswithMultipleScalingBehaviors .......... 96 6.3CharacterizingHiddenFrequencies ...................... 99 7CONCLUSIONS ................................... 104 REFERENCES ....................................... 105 BIOGRAPHICALSKETCH ................................ 115 6


Figure page 1-1Exampleofatrajectoryinthephasespace ..................... 14 1-2Giantoceanwave(tsunami). ............................. 17 1-3Exampleofseaclutterdata. ............................. 18 1-4ExampleofHeartratevariabilitydataforanormalsubject. ........... 19 2-1Time-dependentexponent(k)vs.evolutiontimekcurvesforLorenzsystem. 32 2-2Time-dependentexponent(k)vs.evolutiontimekcurvesfortheMackey-Glasssystem. ........................................ 33 2-3Time-dependentexponent(t)vs.tcurvesfortimeseriesgeneratedfromlogisticmap. .......................................... 37 2-4Time-dependentexponent(t)vs.tcurvesfortimeseriesgeneratedfromHenonmap. .......................................... 40 3-1WhitenoiseandBrownianmotion ......................... 46 3-2SeveralfBmprocesseswithdierentH. ....................... 48 3-3Time-dependentexponent(k)vs.lnkcurvesforfractionalBrownianmotion. 51 3-4ExampleofON/OFFprocesses. ........................... 51 3-5Time-dependentexponent(k)vs.lnkcurvesforON/OFFprocesseswithParetodistributedONandOFFperiods. .......................... 53 3-6ExamplesforBrownianmotionsandLevymotions. ................ 57 3-7Time-dependentexponent(k)vs.lnkcurvesforLevyprocesses. ........ 58 4-1Exampleofquasi-periodiccongestionwindowsizeW(i)timeseries. ....... 63 4-2ExamplesofchaoticT(i)timeseries. ........................ 65 4-3Complementarycumulativedistributionfunctions(CCDFs)forthetwochaoticT(i)timeseriesofFigs. 4-2 (a,c). .......................... 66 4-4ThetimeseriesT(i)extractedfromaW(i)timeseriescollectedusingnet100instrumentsoverORNL-LSUconnection. ...................... 68 4-5TimeseriesforthecongestionwindowsizecwndforORNL-GaTechconnection. 72 4-6Time-dependentexponent(k)vs.kcurvesforcwnddatacorrespondingtoFig. 4-5 ........................................ 73 7


4-5 ........................................ 75 5-1Collectionofseaclutterdata ............................ 79 5-2Typicalseaclutteramplitudedata. ......................... 79 5-3Twoshortsegmentsoftheamplitudeseaclutterdataseverelyaectedbyclipping. 80 5-4Examplesofthetime-dependentexponent(k)vs.kcurvesfortheseaclutterdata. .......................................... 81 5-5VariationsoftheLyapunovexponentestimatedbyconventionalmethodsvs.the14rangebinsforthe10HHmeasurements. .................. 84 5-6Variationsofthescale-dependentexponentcorrespondingtolargescalevs.therangebinsforthe10HHmeasurements. ...................... 85 5-7Examplesofthe(k)vs.lnkcurvesfortheseaclutterdata. ........... 87 5-8Variationoftheparameterwiththe14rangebins. ............... 87 5-9FrequenciesofthebinswithouttargetsandthebinswithprimarytargetsfortheHHdatasets. ................................... 88 6-1Scale-dependentLyapunovexponent()curvesfortheLorenzsystemandlogisticmap. .......................................... 93 6-2Scale-dependentLyapunovexponent()curvefortheMackey-Glasssystem. .. 94 6-3ThefunctionsF(x)andG(x). ............................ 97 6-4Scale-dependentLyapunovexponent()forthemodeldescribedbyEq. 6{9 . 97 6-5Power-spectraldensity(PSD)fortimeseriesgeneratedfromLorenzsystem. .. 100 6-6Lorenzattractor. ................................... 101 6-7Hiddenfrequencyphenomenonoflaserintensitydata. ............... 102 8


Complexsystemsusuallycomprisemultiplesubsystemswithhighlynonlineardeterministicandstochasticcharacteristics,andareregulatedhierarchically.Thesesystemsgeneratesignalswithcomplexcharacteristicssuchassensitivedependenceonsmalldisturbances,longmemory,extremevariations,andnonstationarity.Chaostheoryandrandomfractaltheoryaretwoofthemostimportanttheoriesdevelopedforanalyzingcomplexsignals.However,theyhaveentirelydierentfoundations,onebeingdeterministicandtheotherbeingrandom.Tosynergisticallyusethesetwotheories,wedevelopedanewtheoreticalframework,power-lawsensitivitytoinitialconditions(PSIC),toencompassbothchaosandrandomfractaltheoriesasspecialcases.Toshowthepowerofthisframework,weappliedittostudytwochallengingandimportantproblems:Internetdynamicsandseaclutterradarreturns.WeshowedthatPSICcanreadilycharacterizethecomplicatedInternetdynamicsduetointerplayofnonlinearAIMD(additiveincreaseandmultiplicativedecrease)operationofTCPandstochasticuserbehavior,androbustlydetectlowobservabletargetsfromseaclutterradarreturnswithhighaccuracy. Wealsodevelopedanewmultiscalecomplexitymeasure,scale-dependentLyapunovexponent(SDLE),thatcanbecomputedfromshortnoisydata.Themeasurereadilyclassiedvarioustypesofcomplexmotions,andsimultaneouslycharacterizedbehaviorsofcomplexsignalsonawiderangeofscales,includingcomplexirregularbehaviorsonsmallscales,andorderlybehaviors,suchasoscillatorymotions,onlargescales. 9


Adynamicalsystemisonethatevolvesintime.Dynamicalsystemscanbestochastic(inwhichcasetheyevolveaccordingtosomerandomprocesssuchasthetossofacoin)ordeterministic(inwhichcasethefutureisuniquelydeterminedbythepastaccordingtosomeruleormathematicalformula).Whetherthesysteminquestionsettlesdowntoequilibrium,keepsrepeatingincycles,ordoessomethingmorecomplicated,dynamicsarewhatweusetoanalyzethebehaviorofsystems. Complexdynamicalsystemsusuallycomprisemultiplesubsystemswithhighlynonlineardeterministicandstochasticcharacteristics,andareregulatedhierarchically.Thesesystemsgeneratesignalswithcomplexcharacteristicssuchassensitivedependenceonsmalldisturbances,longmemory,extremevariations,andnonstationarity[ 1 ].Astockmarket,forexample,isstronglyinuencedbymulti-layereddecisionsmadebymarketmakers,aswellasinuencedbyinteractionsofheterogeneoustraders,includingintradaytraders,short-periodtraders,andlong-periodtraders,andthusgivesrisetohighlyirregularstockprices.TheInternet,asanotherexample,hasbeendesignedinafundamentallydecentralizedfashionandconsistsofacomplexwebofserversandroutersthatcannotbeeectivelycontrolledoranalyzedbytraditionaltoolsofqueuingtheoryorcontroltheory,andgiverisetohighlyburstyandmultiscaletracwithextremelyhighvariance,aswellascomplexdynamicswithbothdeterministicandstochasticcomponents.Similarly,biologicalsystems,beingheterogeneous,massivelydistributed,andhighlycomplicated,oftengeneratenonstationaryandmultiscalesignals.Withtherapidaccumulationofcomplexdatainlifesciences,systemsbiology,nano-sciences,informationsystems,andphysicalsciences,ithasbecomeincreasinglyimportanttobeabletoanalyzemultiscaleandnonstationarydata. Multiscalesignalsbehavedierentlydependinguponwhichscalethedataarelookedat.Inordertofullycharacterizesuchcomplexsignals,theconceptofscalehastobe 10


Thisdissertationaimstobuildaneectivearsenalformultiscalesignalprocessingbysynergisticallyintegratingapproachesbasedonchaostheoryandrandomfractaltheory,andgoingbeyond,tocomplementconventionalapproachessuchasspectralanalysisandmachinelearningtechniques.Tomakesuchanintegrationpossible,twoimportanteortshavebeenmade:Power-lawsensitivitytoinitialconditions(PSIC):Wedevelopedanewtheoreticalframeworkforsignalprocessing,tocreateacommonfoundationforchaostheoryandrandomfractaltheory,sothattheycanbebetterintegrated.Scale-dependentLyapunovexponent(SDLE):Wedevelopedanexcellentmultiscalemeasure,whichisavariantofthenitesizeLyapunovexponent(FSLE).Weproposedahighlyecientalgorithmforcalculatingit,andshowedthatitcanreadilyclassifydierenttypesofmotions,includingtrulylow-dimensionalchaos,noisychaos,noise-inducedchaos,random1=fand-stableLevyprocesses,andcomplexmotionswithchaoticbehavioronsmallscalesbutdiusivebehavioronlargescales.Themeasurecanaptlycharacterizecomplexbehaviorsofrealworldmultiscalesignalsonawiderangeofscales. Intherestofthischapter,weshallrstpresentanoverviewofdynamicssystems,especiallyweshallpointoutwhynonlinearsystemsaremuchhardertoanalyzethanlinearones.Wethengiveafewexamplesofmultiscalephenomenatoshowthedicultiesandexcitementofmultiscalesignalprocessing.Afterthat,weshallbrieyintroducesomebackgroundknowledgeaboutchaosandfractalgeometry.Thenwediscusstheimportanceofconnectingchaostheoryandrandomfractaltheoryforcharacterizingthebehaviorsofmultiscalesignalsonawiderangeofscales.Wealsoemphasizetheimportanceofexplicitlyincorporatingtheconceptofscaleindevisingmeasuresforcharacterizingmultiscalesignals.Finally,weoutlinethescopeofthepresentstudy. 11


dt2+bdx dt+kx=0(1{1) isanordinarydierentialequation,becauseitinvolvesonlyordinaryderivativesdx=dtandd2x=dt2.Thatis,thereisonlyoneindependentvariable,thetimet.Incontrast,theheatequation@u @t=@2u @x2 Averygeneralframeworkforordinarydierentialequationsisprovidedbythesystem _x1=f1(x1;;xn)..._xn=fn(x1;;xn):(1{2) Heretheoverdotsdenotedierentiationwithrespecttot.Thus_xi=dxi=dt.Thevariablesx1;;xnmightrepresentconcentrationsofchemicalsinareactor,populationsofdierentspeciesinanecosystem,orthepositionsandvelocitiesoftheplanetsinthesolarsystem.Thefunctionsf1;;fnaredeterminedbytheproblemathand. Forexample,thedampedoscillator(Eq. 1{1 )canberewrittenintheformof(Eq. 1{2 ),thankstothefollowingtrick:weintroducenewvariablesx1=xandx2=_x. 12


1{2 )is_x1=x2_x2=b mx2k mx1: Forexample,theswingingofapendulumisgovernedbytheequationx+g Lsinx=0; Lsinx1: Itturnsoutthatthependulumequationcanbesolvedanalytically,intermsofellipticfunctions.Butthereoughttobeaneasierway.Afterall,themotionofthependulumissimple:atlowenergy,itswingsbackandforth,andathighenergyitwhirlsoverthetop.Thereshouldbesomewayofextractingthisinformationfromthesystemdirectly|usinggeometricmethods. Hereistheroughidea.Supposewehappentoknowasolutiontothependulumsystem,foraparticularinitialcondition.Thissolutionwouldbeapairoffunctionsx1(t)andx2(t),representingthepositionandvelocityofthependulum.Ifweconstructan 13


Exampleofatrajectoryinthephasespace abstractspacewithcoordinates(x1;x2),thenthesolution(x1(t);x2(t))correspondstoapointmovingalongacurveinthisspace,asshownFig. 1-1 .Thiscurveiscalledatrajectory,andthespaceiscalledthephasespaceforthesystem.Thephasespaceiscompletelylledwithtrajectories,sinceeachpointcanserveasaninitialcondition. Ourgoalistorunthisconstructioninreverse:giventhesystem,wewanttodrawthetrajectories,andtherebyextractinformationaboutthesolutions.Inmanycases,geometricreasoningwillallowustodrawthetrajectorieswithoutactuallysolvingthesystem! Thephasespaceforthegeneralsystem(Eq. 1{2 )isthespacewithcoordinatesx1;;xn.Becausethisspaceisn-dimensional,wewillrefertoEq. 1{2 asann-dimensionalsystemorannth-ordersystem.Thusnrepresentsthedimensionofthephasespace. Aswehavementionedearlier,mostnonlinearsystemsareimpossibletosolveanalytically.Whyarenonlinearsystemssomuchhardertoanalyzethanlinearones?Theessentialdierenceisthatlinearsystemscanbebrokendownintoparts.Theneachpartcanbesolvedseparatelyandnallyrecombinedtogettheanswer.Thisideaallowsa 14


Butmanythingsinnaturedonotactthisway.Wheneverpartsofasysteminterfere,orcooperate,orcompete,therearenonlinearinteractionsgoingon.Mostofeverydaylifeisnonlinear,andtheprincipleofsuperpositionfailsspectacularly.Ifyoulistentoyourtwofavoritesongsatthesametime,youwillnotgetdoublethepleasure!Withintherealmofphysics,nonlinearityisvitaltotheoperationofalaser,theformationofturbulenceinauid,andthesuperconductivityofJosephsonjunctions. 15


Also,ithasbeenobservedthatthefailureofasingleroutermaytriggerroutinginstability,whichmaybesevereenoughtoinstigatearouteapstorm.Furthermore,packetsmaybedeliveredoutoforderorevengetdropped,andpacketreorderingisnotapathologicalnetworkbehavior.AsthenextgenerationInternetapplicationssuchasremoteinstrumentcontrolandcomputationalsteeringarebeingdeveloped,anotherfacetofcomplexmultiscalebehaviorisbeginningtosurfaceintermsoftransportdynamics.Thenetworkingrequirementsforthesenextgenerationapplicationsbelongto(atleast)twobroadclassesinvolvingvastlydisparatetimescales:Highbandwidths,typicallymultiplesof10Gbps,tosupportbulkdatatransfers,Stablebandwidths,typicallyatmuchlowerbandwidthssuchas10to100Mbps,tosupportinteractive,steeringandcontroloperations. 1-2 .Tobequantitative,inFig. 1-3 ,two0.1sdurationseacluttersignals,sampledwithafrequencyof1KHz,areplottedinFig. 1-3 (a,b),a2sdurationsignalisplottedinFig. 1-3 (c),andanevenlongersignal(about130s)isplottedinFig. 1-3 (d).Itisclearthatthesignalisnotpurelyrandom,sincethewaveformcanbefairlysmoothonshorttimescales(Fig. 1-3 (a)).However,thesignalishighlynonstationary,sincethefrequencyofthesignal(Fig. 1-3 (a,b))aswellastherandomnessofthesignal(Fig. 1-3 (c,d))changeovertimedrastically.OnethuscanperceivethatnaiveFourieranalysisordeterministicchaoticanalysisofseacluttermay 16


Giantoceanwave(tsunami).Supposeoureldofobservationincludesthewavetipoflengthscaleofafewmeters,itisthenclearthatthecomplexityofseaclutterismainlyduetomassivereectionofradarpulsesfromawavyandeventurbulentoceansurface. notbeveryuseful.FromFig. 1-3 (e),whereX(m)tisthenon-overlappingrunningmeanofXoverblocksizem,andXistheseaclutteramplitudedata,itcanbefurtherconcludedthatneitherautoregressive(AR)models(Asanexample,wegivethedenitionfortherst-orderARmodel,whichisgivenbyxn+1=axn+n,wheretheconstantcoecientasatises06=jaj<1andnisawhitenoisewithzeromean.)nortextbookfractalmodelscandescribethedata.ThisisbecauseARmodelingrequiresexponentiallydecayingautocorrelation(whichamountstoVar(X(m)t)m1,oraHurstparameterof1/2.SeelaterchaptersforthedenitionoftheHurstparameter),whilefractalmodelingrequiresthevariationbetweenVar(X(m)t)andmtofollowapower-law.However,neitherbehaviorisobservedinFig. 1-3 (e).Indeed,albeitextensiveworkhasbeendoneonseaclutter,thenatureofseaclutterisstillpoorlyunderstood.Asaresult,theimportantproblemof 17


Exampleofseaclutterdata.(a,b)two0.1sdurationsignal;(c)a2sdurationseacluttersignal;(d)theentireseacluttersignal(ofabout130s);and(e)log2[m2Var(X(m))]vs.log2m,whereX(m)=fX(m)t:X(m)t=(Xtmm+1++Xtm)=m;t=1;2;gisthenon-overlappingrunningmeanofX=fXt:t=1;2;goverblocksizem,andXistheseaclutteramplitudedata.TobetterseethevariationofVar(X(m)t)withm,Var(X(m)t)ismultipliedbym2.Whentheautocorrelationofthedatadecaysexponentiallyfast(suchasmodeledbyanautoregressive(AR)process),Var(X(m)t)m1.HereVar(X(m)t)decaysmuchfaster.Afractalprocesswouldhavem2Var(X(m)t)m.However,thisisnotthecase.Therefore,neitherARmodelingnoridealtextbookfractaltheorycanbereadilyappliedhere. targetdetectionwithinseaclutterremainsatremendouschallenge.Weshallreturntoseaclutterlater. 18


ExampleofHeartratevariabilitydataforanormalsubject.(a)theentiresignal,(b,c)thesegmentsofsignalsindicatedasAandBin(a);(d,e)powerspectraldensityforthesignalsshownin(b,c). maintainedconstantandnoperturbinginuencescanbeidentied.IthasbeenobservedthatHRVisrelatedtovariouscardiovasculardisorders.Therefore,analysisofHRVisveryimportantinmedicine.However,thistaskisverydicult,sinceHRVdataarehighlycomplicated.AnexampleisshowninFig. 1-4 ,foranormalyoungsubject.Evidently,thesignalishighlynonstationaryandmultiscaled,appearingoscillatoryforsomeperiodoftime(Figs. 1-4 (b,d)),andthenvaryingasapower-lawforanotherperiodoftime(Figs. 1-4 (c,e)).Thelatterisanexampleoftheso-called1=fprocesses,whichwillbediscussedindepthinlaterchapters.Whilethemultiscalenatureofsuchsignalscannotbefullycharacterizedbyexistingmethods,thenonstationarityofthedataisevenmoretroublesome,sinceitrequiresthedatatobeproperlysegmentedbeforefurtheranalysis 19


Atthecenterofchaostheoryistheconceptofsensitivedependenceoninitialconditions:averyminordisturbanceininitialconditionsleadstoentirelydierentoutcomes.AnoftenusedmetaphorillustratingthispointisthatasunnyweatherinNewYorkcouldbereplacedbyarainyonesometimeinthefutureafterabutteryapsitswingsinBoston.Suchafeaturecontrastssharplywiththetraditionalview,largelybasedonourexperiencewithlinearsystems,thatsmalldisturbances(orcauses)canonlygenerateproportionaleects,andthatinorderforthedegreeofrandomnesstoincrease,thenumberofdegreesoffreedomhastobeinnite. Nodenitionofthetermchaosisuniversallyacceptedyet,butalmosteveryonewouldagreeonthethreeingredientsusedinthefollowingworkingdenition: Chaosisaperiodiclong-termbehaviorinadeterministicsystemthatexhibitssensitivedependenceoninitialconditions.\Aperiodiclong-termbehavior"meansthattherearetrajectorieswhichdonotsettledowntoxedpoints,periodicorbits,orquasi-periodicorbitsast!1.Forpracticalreasons,weshouldrequirethatsuchtrajectoriesarenottoorare.Forinstance,wecould 20


Chaosisalsocommonlycalledastrangeattractor.Thetermattractorisalsodiculttodeneinarigorousway.Wewantadenitionthatisbroadenoughtoincludeallthenaturalcandidates,butrestrictiveenoughtoexcludetheimposters.Thereisstilldisagreementaboutwhattheexactdenitionshouldbe. Looselyspeaking,anattractorisasettowhichallneighboringtrajectoriesconverge.Stablexedpointsandstablelimitcyclesareexamples.Moreprecisely,wedeneanattractortobeaclosedsetAwiththefollowingproperties: 1. 2. 3. Nowwecandeneastrangeattractortobeanattractorthatexhibitssensitivedependenceoninitialconditions.Strangeattractorswereoriginallycalledstrangebecausetheyareoftenfractalsets.Nowadaysthisgeometricpropertyisregardedaslessimportantthanthedynamicalpropertyofsensitivitydependenceoninitialconditions.Thetermschaoticattractorandfractalattractorareusedwhenonewishestoemphasizeoneortheotherofthoseaspects. 21


Fornow,weshallbesatisedwithanintuitivedenitionofafractal:asetthatshowsirregularbutself-similarfeaturesonmanyorallscales.Self-similaritymeansapartofanobjectissimilartootherpartsortothewhole.Thatis,ifweviewanirregularobjectwithamicroscope,whetherweenlargetheobjectby10times,orby100times,orevenby1000times,wealwaysndsimilarobjects.Tounderstandthisbetter,letusimaginewewereobservingapatchofwhiteclouddriftingawayinthesky.Oureyeswereratherpassive:wewerestaringmoreorlessatthesamedirection.Afterawhile,thepartofthecloudwesawdriftedaway,andwewereviewingadierentpartofthecloud.Nevertheless,ourfeelingremainsmoreorlessthesame. Mathematically,fractalischaracterizedbyapower-lawrelation,whichtranslatestoalinearrelationinlog-logscale.Fornow,letusagainresorttoimagination|wewerewalkingdownawildandveryjaggedmountaintrailorcoastline.Wewouldliketoknowthedistancecoveredbyourroute.Supposeourrulerhasalengthof|whichcouldbeourstep-size,anddierenthikersmayhavedierentstep-sizes|apersonridingahorsehasahugestep-size,whileagroupofpeoplewithalittlechildmusthaveatinystep-size.Thelengthofourrouteis whereN()isthenumberofintervalsneededtocoverourroute.ItismostremarkablethattypicallyN()scaleswithinapower-lawmanner, withDbeinganon-integer,1




2 ,werstintroducethenewconceptofPSIC.ThenweextendtheconceptofPSICtohigh-dimensionalcase.WeshallpresentageneralcomputationalframeworkforassessingPSICfromrealworlddata.Wealsoapplytheproposedcomputationalproceduretostudynoise-freeandnoisylogisticandHenonmapsattheedgeofchaos.weshowthatwhennoiseisabsent,PSICishardtodetectfromascalartimeseries.However,whenthereisdynamicnoise,motionsaroundtheedgeofchaos,beitsimplyregularortrulychaoticwhenthereisnonoise,allcollapseontothePSICattractor.Hence,dynamicnoisemakesPSICobservable.InChapter 3 ,wedemonstratethatthenewframeworkofPSICcanreadilycharacterizerandomfractalprocesses,including1=fprocesseswithlong-range-correlationsandLevyprocesses.FromourstudyinChapters 2 and 3 ,itisclearthattheconceptofPSICnotjustbridgesstandardchaostheoryandrandomfractaltheory,butinfactprovidesamoregeneralframeworktoencompassboththeoriesasspecialcases.ToillustratethepoweroftheframeworkofPSIC,inChapters 4 and 5 ,weapplyittostudytwoimportantbutchallengingproblems:Internetdynamicsandseaclutterradarreturns,respectively.Bothareoutstandingexamplesofcomplexdynamicalsystemswithbothnonlinearityandrandomness.ThenonlinearityofInternetdynamicscomesfromAIMD(additiveincreaseandmultiplicativedecrease)operationofTCP(TransmissionControlprotocol),whilethestochasticcomponentcomesfromtherandomuserbehavior.Seaclutterisanonlineardynamicalprocesswithstochasticfactorsduetointerferenceofvariouswindandswellwavesandtolocalatmosphericturbulence.WeshallshowthatthenewtheoreticalframeworkofPSICcaneectivelycharacterizethecomplicatedInternetdynamicsaswellastorobustlydetectlowobservabletargetsfromseaclutterradarreturnswithhighaccuracy.InChapter 6 ,weshallintroducethebasictheoryofanexcellentmultiscalemeasure|thescale-dependentLyapunovexponent(SDLE),anddevelopaneectivealgorithmtocomputethemeasure.To 25


7 26


Theformalismofnonextensivestatisticalmechanics(NESM)[ 2 ]hasfoundnumerousapplicationstothestudyofsystemswithlong-range-interactions[ 3 { 6 ]andmultifractalbehavior[ 7 8 ],andfullydevelopedturbulence[ 8 { 10 ],amongmanyothers.Inordertocharacterizeatypeofmotionwhosecomplexityisneitherregularnorfullychaotic/random,recentlytheconceptofexponentialsensitivitytoinitialconditions(ESIC)ofdeterministicchaoshasbeengeneralizedtopower-lawsensitivitytoinitialconditions(PSIC)[ 7 11 ].Mathematically,theformulationofPSICcloselyparallelsthatofNESM.PSIChasbeenappliedtothestudyofdeterministic1-Dlogisticlikemapsand2-DHenonmapattheedgeofchaos[ 7 11 { 20 ],yieldingconsiderablenewinsights.Inthischapter,werstbrieydescribesomebackgroundknowledgeaboutchaos.Inparticular,wewilldiscussadirectdynamicaltestfordeterministicchaosdevelopedbyGaoandZheng[ 21 22 ].Thismethodoersamorestringentcriterionfordetectinglow-dimensionalchaos,andcansimultaneouslymonitormotionsinphasespaceatdierentscales.ThenweintroducethemoregeneralizedconceptofPSIC,andextendthenewconcepttohigh-dimensionalcase.Specically,wepresentageneralcomputationalframeworkforassessingPSICinatimeseries.Wealsoapplytheproposedcomputationalproceduretostudytwomodelsystems:noise-freeandnoisylogisticandHenonmapsattheedgeofchaos. 27


where1ispositiveandcalledthelargestLyapunovexponent.Duetotheboundednessoftheattractorandtheexponentialdivergencebetweennearbytrajectories,astrangeattractortypicallyisafractal,characterizedbyasimpleandelegantscalinglawasdenedbyEq. 1{4 .Anon-integralfractaldimensionofanattractorcontrastssharplywiththeinteger-valuedtopologicaldimension(whichis0fornitenumberofisolatedpoints,1forasmoothcurve,2forasmoothsurface,andsoon). Conventionally,ithasbeenassumedthatatimeserieswithanestimatedpositivelargestLyapunovexponentandanon-integralfractaldimensionischaotic.However,ithasbeenfoundthatthisassumptionmaynotbeasucientindicationofdeterministicchaos.Onecounter-exampleistheso-called1=fprocesses.Suchprocesseshavespectraldensity where0

23 { 27 ],estimationofthestrengthofmeasurementnoiseinexperimentaldata[ 28 29 ],pathologicaltremors[ 30 ],shear-thickeningsurfactantsolutions[ 31 ],diluteshearedaqueoussolutions[ 32 ],andserratedplasticows[ 33 ]. Nowletusbrieydescribethedirectdynamicaltestfordeterministicchaos[ 21 22 ].Givenascalartimeseries,x(1);x(2);:::;x(N),(assuming,forconvenience,thattheyhavebeennormalizedtotheunitinterval[0,1]),onerstconstructsvectorsofthefollowingformusingthetimedelayembeddingtechnique[ 34 { 36 ]: withmbeingtheembeddingdimensionandLthedelaytime.Forexample,whenm=3andL=4,wehaveV1=[x(1);x(5);x(9)],V2=[x(2);x(6);x(10)],andsoon.Fortheanalysisofpurelychaoticsignals,mandLhastobechosenproperly.Thisistheissueofoptimalembedding(see[ 21 22 37 ]andreferencestherein).Afterthescalartimeseriesisproperlyembedded,onethencomputesthetimedependentexponent(k)curves (k)=lnkVi+kVj+kk kViVjk(2{4) withrkViVjkr+r,whererandrareprescribedsmalldistances.Theanglebracketsdenoteensembleaveragesofallpossiblepairsof(Vi;Vj).Theintegerk,calledtheevolutiontime,correspondstotimekt,wheretisthesamplingtime.Notethatgeometrically(r;r+r)denesashell,andashellcapturesthenotionofscale.Thecomputationistypicallycarriedoutforasequenceofshells.ComparingEq. 2{4 withEq. 2{1 ,onenoticesthatkVi+kVj+kkplaystheroleofd(t),whilekViVjkplaystheroleofd(0).Figure 2-1 (a)showsthe(k)vs.kcurvesforthechaoticLorenzsystemdrivenby 29


dt=16(xy)+D1(t)dy dt=xz+45:92xy+D2(t)dz dt=xy4z+D3(t)(2{5) with=0,=ij(tt0),i;j=1;2;3.NotethatD2isthevarianceoftheGaussiannoiseterms,andD=0describesthecleanLorenzsystem.WeobservefromFig. 2-1 (a)thatthe(k)curvesarecomposedofthreeparts.Thesecurvesarelinearlyincreasingfor0kka.Theyarestilllinearlyincreasingforkakkp,butwithaslightlydierentslope.Theyareatforkkp.NotethattheslopeofthesecondlinearlyincreasingpartgivesanestimateofthelargestpositiveLyapunovexponent[ 21 22 ].kaisrelatedtothetimescaleforapairofnearbypoints(Vi;Vj)toevolvetotheunstablemanifoldofViorVj.Itisontheorderoftheembeddingwindowlength,(m1)L.kpisthepredictiontimescale.Itislongerforthe(k)curvesthatcorrespondtosmallershells.Thedierencebetweentheslopesoftherstandsecondlinearlyincreasingpartsiscausedbythediscrepancybetweenthedirectiondenedbythepairofpoints(Vi;Vj)andtheunstablemanifoldofViorVj.ThisfeaturewasrstobservedbySatoetal.[ 38 ]andwasusedbythemtoimprovetheestimationoftheLyapunovexponent.Therstlinearlyincreasingpartcanbemadesmallerorcanevenbeeliminatedbyadjustingtheembeddingparameterssuchasbyusingalargervalueform.Notethatthesecondlinearlyincreasingpartsofthe(k)curvescollapsetogethertoformanenvelope.Itisthisveryfeaturethatformsthedirectdynamicaltestfordeterministicchaos[ 21 22 ].Thisisbecausethe(k)curvesfornoisydata,suchaswhitenoiseor1=fprocesses,areonlycomposedoftwoparts,anincreasing(butnotlinear)partfork(m1)Landaatpart[ 21 22 ].Furthermore,dierent(k)curvesfornoisydataseparatefromeachother;henceanenvelopeisnotdened.Wealsonote,mostotheralgorithmsforestimatingthelargestLyapunovexponentisequivalenttocompute(k)forr

Onecanexpectthatthebehaviorofthe(k)curvesforanoisychaoticsystemliesinbetweenthatofthe(k)curvesforacleanchaoticsystemandthatofthe(k)curvesforwhitenoiseorfor1=fprocesses.Thisisindeedso.AnexampleisshowninFig. 2-1 (b)forthenoisyLorenzsystemwithD=4.WeobservefromFig. 2-1 (b)thatthe(k)curvescorrespondingtodierentshellsnowseparate.Therefore,anenvelopeisnolongerdened.Thisseparationislargerbetweenthe(k)curvescorrespondingtosmallershells,indicatingthattheeectofnoiseonthesmall-scaledynamicsisstrongerthanthatonthelarge-scaledynamics.Alsonotethatka+kpisnowontheorderoftheembeddingwindowlength,andisalmostthesameforallthe(k)curves.Withstrongernoise(D>4),the(k)curveswillbemorelikethoseforwhitenoise[ 21 22 ]. Finally,weconsidertheMackey-Glassdelaydierentialsystem[ 39 ].Whena=0:2;b=0:1;c=10;=30,ithastwopositiveLyapunovexponents,withthelargestLyapunovexponentcloseto0.007[ 21 ]. 1+x(t+)cbx(t)(2{6) HavingtwopositiveLyapunovexponentswhilethevalueofthelargestLyapunovexponentofthesystemisnotmuchgreaterthan0,onemightbeconcernedthatitmaybediculttodealwiththissystem.Thisisnotthecase.Infact,thissystemcanbeanalyzedasstraightforwardlyasotherdynamicalsystemsincludingtheLogisticmap,HenonmapandtheRosslersystem.Anexampleofthe(k)vs.kcurvesforthechaoticMackey-Glass 31


Time-dependentexponent(k)vs.evolutiontimekcurvesfor(a)cleanand(b)noisyLorenzsystem.Sixcurves,frombottomtotop,correspondtoshells(2(i+1)=2;2i=2)withi=8,9,10,11,12,and13.Thesamplingtimeforthesystemis0.03sec,andembeddingparametersarem=4,L=3.5000pointsareusedinthecomputation. 32


Time-dependentexponent(k)vs.evolutiontimekcurvesfortheMackey-Glasssystem.Ninecurves,frombottomtotop,correspondtoshells(2(i+1)=2;2i=2)withi=5;6;,and13.Thecomputationwasdonewithm=5;L=1,and5000pointssampledwithatimeintervalof6. delaydierentialsystemisshowninFig. 2-2 ,wherewehavefollowed[ 21 ]andusedm=5;L=1,and5000pointssampledwithatimeintervalof6.Clearly,weobservethatthereexistsawelldenedcommonenvelopetothe(k)curves.Actually,theslopeoftheenvelopeestimatesthelargestLyapunovexponent.Thisexampleillustratesthatwhenoneworksintheframeworkofthe(k)curvesdevelopedbyGaoandZheng[ 21 22 ],onedoesnotneedtobeveryconcernedaboutnon-uniformgrowthrateinhigh-dimensionalsystems. 2{7 ,wherex(0)istheinnitesimaldiscrepancyintheinitialcondition,x(t)isthediscrepancyattimet>0. x(0)(2{7) 33


Tsallisandco-workers[ 7 11 ]havegeneralizedEq. 2{9 to whereqiscalledtheentropicindex,andqisinterpretedtobeequaltoKq,thegeneralizationoftheKolmogorov-Sinaientropy.Eq. 2{10 denesthePSICinthe1-Dcase.Obviously,PSICreducestoESICwhenq!1.ThesolutiontoEq. 2{10 is Whentislargeandq6=1,(t)increaseswithtasapower-law, whereC=[(1q)q]1=(1q).ForEq. 2{12 todeneanunstablemotionwithq>0,wemusthaveq1.Laterweshallmapdierenttypesofmotionstodierentrangesofq. ToapplyPSICtotheanalysisoftimeseriesdata,onecanrstconstructaphasespacebyconstructingembeddingvectorsViasdenedbyEq. 2{3 .Eq. 2{11 canthenbegeneralizedtohigh-dimensionalcase[ 40 ], jjV(0)jj=1+(1q)(1)qt1=(1q)(2{13) wherejjV(0)jjistheinnitesimaldiscrepancybetweentwoorbitsattime0,jjV(t)jjisthedistancebetweenthetwoorbitsattimet>0,qistheentropicindex,and(1)qistherstq-Lyapunovexponent,correspondingtothepower-lawincreaseoftherstprincipleaxisofaninnitesimalballinthephasespace.(1)qmaynotbeequaltoKq.Thisisunderstoodbyrecallingthatforchaoticsystems,theKolmogorov-Sinaientropyisthesum 34


2{13 againgivesapower-lawincreaseof(t)witht.Notethatundertheaboveframework,themotionmaynotbelikefullydevelopedchaos,thusnotergodic[ 41 ]. WenowconsiderthegeneralcomputationalframeworkforPSIC.Givenanitetimeseries,theconditionofjjV(0)jj!0cannotbesatised.Inordertondthelawgoverningthedivergenceofnearbyorbits,onethushastoexaminehowtheneighboringpoints,(Vi;Vj),inthephasespace,evolvewithtime,byformingsuitableensembleaverages.Noticethatif(Vi1;Vj1)and(Vi2;Vj2)aretwopairsofnearbypoints,whenjjVi1Vj1jjjjVi2Vj2jj1,thentheseparationssuchasjjVi1+tVj1+tjjandjjVi2+tVj2+tjjcannotbesimplyaveragedtoprovideestimatesforqand(1)q.Infact,itwouldbemostconvenienttoconsiderensembleaveragesofpairsofpoints(Vi;Vj)thatallfallwithinaverythinshell,r1jjViVjjjr2,wherer1andr2areclose.Theseargumentssuggestthatthetime-dependentexponentcurvesdenedbyEq. 2{4 provideanaturalframeworktoassessPSICfromatimeseries.Thisisindeedso.Infact,wehave ln(t)(t)=lnkVi+tVj+tk kViVjk(2{14) Now,bythediscussionsinSec. 2.1 ,itisclearthatPSICisageneralizationofESIC:aslongasthe(k)curvesfromdierentshellsformalinearenvelope,thenq=1andthemotionischaotic.Thenextquestionis:DoesPSICalsoincluderandomfractalsasspecialcases?Theanswerisyes.Itwillbegiveninthenextchapter.There,wewillalsogainabetterunderstandingofthemeaningof(1)q.TherestofthischaptershowsexamplesofcharacterizingtheedgeofchaosbythenewframeworkofPSIC.Specically,westudytimeseriesgeneratedfromnoise-freeandnoisylogisticandHenonmaps. 35


26 ],andphysiologicalnoiseinbiologicalsystems[ 30 ]),westudyboththedeterministicandthenoisylogisticmap: whereaisthebifurcationparameterandnisawhiteGaussiannoisewithmeanzeroandvariance1.Theparametercharacterizesthestrengthofnoise.Forthecleansystem(=0),theedgeofchaosoccursattheaccumulationpoint,a1=3:569945672.Weshallstudythreeparametervalues,a1=a10:001,a1,anda2=a1+0:001.Whennoiseisabsent,a1correspondstoaperiodicmotionwithperiod25,whilea2correspondstoatrulychaoticmotion.Weshallonlystudytransient-freetimeseries.InFigs. 2-3 (a-c),wehaveplottedthe(t)vs.tcurvesforparametervaluesa1,a1,anda2,respectively.WeobservefromFig. 2-3 (a)thatthevariationof(t)withtisperiodic(withperiod16,whichishalfoftheperiodofthemotion)whenthemotionisperiodic.Thisisagenericfeatureofthe(t)curvesfordiscreteperiodicattractors,whentheradiusoftheshellislargerthanthesmallestdistancebetweentwopointsontheattractor(whenaperiodicattractoriscontinuous,(t)canbearbitrarilycloseto0).IthasbeenfoundbyTsallisandco-workers[ 11 ]thatattheedgeofchaosforthelogisticmap,(t)isgivenbyEq. 2{11 withq0:2445.Surprisingly,wedonotobservesuchadivergenceinFig. 2-3 (b).Infact,ifoneplots(t)vs.lnt,oneonlyobservesacurvethatincreasesveryslowly(similartothatshowninFig. 2-4 (a)).ThemoreinterestingpatternistheonethatisshowninFig. 2-3 (c),whereweobservealinearlyincreasingly(t)vs.tcurve.Infact,showninFig. 2-3 (c)aretwosuchcurves,correspondingtotwodierentshells.Veryinterestingly,thetwocurvescollapsetoformacommonenvelopeinthelinearlyincreasingpartofthe 36


Time-dependentexponent(t)vs.tcurvesfortimeseriesgeneratedfromthenoise-freelogisticmapwith(a)a1=a10:001,wherethemotionisperiodicwithperiod25,(b)a1=3:569945672,and(c)a2=a1+0:001,wherethemotionischaotic.Plottedin(d-f)are(t)vs.lntcurvesforthenoisylogisticmapwith=0:001.Verysimilarresultswereobtainedwhen=0:0001.Shownin(c-f)areactuallytwocurves,correspondingtotwodierentshells.104pointswereusedinthecomputation,withembeddingparametersm=4;L=1.However,solongasm>1,theresultsarelargelyindependentofembedding.Whenm=1,the(t)curvesarenotsmooth,andtheestimated1=(1q)valueismuchsmallerthanthetheoreticalvalue. curve.TheslopeoftheenvelopegivesagoodestimateofthelargestpositiveLyapunovexponent.Thisisagenericfeatureofchaos[ 21 22 ],asexplainedinSec. 2.1 .Sincethechaosstudiedhereisclosetotheedgeofchaos,thecurvesshowninFig. 2-3 (c)arelesssmooththanthosereportedearlier[ 21 { 23 26 ]. WhycannotthetheoreticalpredictionofPSICattheedgeofchaosforthelogisticmapbeobservedfromacleantimeseries?Inarecentlypublishedveryinterestingand 37


42 ]suggeststhatdynamicnoisemaybeoffundamentalimportancetotheTsallisNESM.MaydynamicnoiseplaysimilarlysignicantroleforPSIC?Asisshownbelow,theanswerisyes.InFigs. 2-3 (d-f),wehaveshownthe(t)vs.lntcurvesforthethreeparametersconsidered,withnoisestrength1=0:001.Infact,shownineachgurearetwocurves,correspondingtotwodierentshells.Theyparallelwitheachother.Theslopesofthosecurvesareabout1.20,closetothetheoreticalvalueof1=(10:2445)1:32.WhileitisverysatisfactorytoobservePSICattheedgeofchaos,itismorethrillingtoobservethecollapseofregularaswellaschaoticmotionsontothePSICattractoraroundtheedgeofchaos.ThissigniesthestablenessofPSICwhenthereisdynamicnoise.ItisimportanttoemphasizethattheresultsshowninFigs. 2-3 (d-f)arelargelyindependentofthenoisestrength,solongasnoiseisneithertooweaknortoostrong.Forexample,verysimilarresultshavebeenobservedwith2=1=10=0:0001. BeforewemoveontodiscusstheHenonmapneartheedgeofchaos,letuscommentonthedierencebetweentheESICandthePSIC.WhentheESICisthecase,the(t)vs.tcurvesarestraightforTatTp,whereTpisapredictiontimescale,andTaisatimescalefortheinitialseparationtoevolvetothemostunstabledirection.When(t)vs.tcurvescorrespondingtodierentshellsareplottedtogether,theycollapsetogetherforaconsiderablerangeoft.WhenthePSICisthecase,the(t)vs.lntcurvesarefairlystraight,andthe(t)vs.lntcurvescorrespondingtodierentshellsseparateandparallelwitheachother.Thisfeatureisthefundamentalreasonthatensembleaveragesaremostconvenientlyformedbyrequiringneighboringpairsofpointstoallfallwithinathinshell.Forexample,iftwodistinctthinshells,describedbyr1jjXiXjjjr2andr2jjViVjjjr3,arejoinedtogethertoformasinglethickershell,describedbyr1jjViVjjjr3,andoneaveragestheseparationbetweenXiandXjbeforetakingthelogarithm,thentheresultingslopebetween(t)andlnt,assumedtobestilllinear,willhavetobesmallerthanthatestimatedfromthetwothinshellsseparately.Atthispointitisalsoworthnotingthatthelinearlyincreasingpartofthe(t)vs.lntcurves 38


2{13 whendealingwithexperimentaltimeseries.Infact,wehavefoundthatifwedoso,theestimated1=(1q)valuefromFigs.1(d-f)increasestoabout1.27,muchclosertothetheoreticalvalueof1.32. NextletusconsidertheHenonmap,wherexandyarewhiteGaussiannoise,andareuncorrelatedwitheachother.Theparametermeasuresthestrengthofthenoise. Tirnakli[ 20 ]studiedthedeterministicversionofthismapattheedgeofchaos,forparametervaluesac=1:40115518;b=0;ac=1:39966671;b=0:001;andac=1:38637288;b=0:01.Wehavestudiedboththenoise-freeandnoisymapforparametervalueslistedabove,andfoundverysimilarresultstothosepresentedinFig. 2-3 .InFigs. 2-4 (a,b),wehaveplotted(t)vs.lntcurvesforthenoise-freeandnoisymapforac=1:38637288;b=0:01.Ascanbeobservedclearly,Fig. 2-4 (b)isverysimilartoFigs. 2-3 (d-f),whiletheslopeforthecurvesinFig. 2-4 (a)ismuchsmallerthanthatinFig. 2-4 (b).Again,wehaveobserved(butnotshownbygureshere,sincetheyareverysimilartoFigs. 2-3 (d-f)andFig. 2-4 (b))thatdynamicnoisemakestheregularandchaoticmotionstocollapseontothePSICattractorforparametersaroundthosedeningtheedgeofchaos,andthatsuchtransitionsarelargelyindependentofthenoisestrength. Tosummarize,wehavedescribedaneasilyimplementableprocedureforcomputationallyexaminingPSICfromatimeseries.Bystudyingtwomodelsystems:noise-freeandnoisylogisticandHenonmapsneartheedgeofchaos,wehavefoundthatwhenthereisno 39


Time-dependentexponent(t)vs.tcurvesfortimeseriesgeneratedfrom(a)thenoise-freeHenonmapwithac=1:38637288;b=0:01.Plottedin(b)aretwo(t)vs.lntcurves(correspondingtotwodierentshells)forthenoisymapwith=0:001.Similarresultswereobtainedwith=0:0001.104pointswereusedinthecomputation,withembeddingparametersm=4;L=1. noise,thePSICattractorcannotbeobservedfromascalartimeseries.However,whendynamicnoiseispresent,motionsaroundtheedgeofchaos,beitsimplyregularortrulychaoticwhenthereisnonoise,allcollapseontothePSICattractor.WhiletheseexaminationssignifytheubiquityofPSIC,theyalsohighlighttheimportanceofdynamicnoise.Theexistenceofthelatterisperhapstheveryreasonthattrulychaotictimeseriescanseldombeobserved. 40


InChapter 2 ,wehavediscussedthatPSICisageneralizationofESIC:aslongasthe(t)curvesfromdierentshellsformalinearenvelope,thenq=1andthemotionischaotic.WehavealsoshownthatedgeofchaoscanbewellcharacterizedbyPSIC.Inthischapter,weshallstudytwomajortypesofrandomfractals:1=fprocesseswithlong-range-correlationsandLevyprocesses.WeshallshowbothanalyticallyandthroughnumericalsimulationsthatbothtypesofprocessescanbereadilycharacterizedbyPSIC. for>0,0

3{3 canbesimplyderivedfromEq. 3{4 ;wehavelisteditasaseparateequationforconvenienceoffuturereference.Ifweconsideronlysecond-orderstatistics,orifaprocessisGaussian,Eq. 3{2 toEq. 3{4 canbeusedinsteadofEq. 3{1 todeneaself-similarprocess. Averyusefulwayofdescribingaself-similarprocessisbyitsspectralrepresentation.Strictlyspeaking,theFouriertransformofX(t)isundened,duetothenonstationarynatureofX(t).Onecan,however,considerX(t)inaniteinterval,say,0

with ProcesseswithpowerspectraldensitiesasdescribedbyEq. 3{7 arecalled1=fprocesses.Typically,thepower-lawrelationshipsthatdenetheseprocessesextendoverseveraldecadesoffrequency.Suchprocesseshavebeenfoundinnumerousareasofscienceandengineering.Someoftheolderliteraturesonthissubjectcanbefound,forexample,inPress[ 45 ],Bak[ 46 ],andWornell[ 47 ].Someofthemorerecentlydiscovered1=fprocessesareintracengineering[ 48 { 50 ],DNAsequence[ 51 { 53 ],humancognition[ 54 ],ambiguousvisualperception[ 55 56 ],coordination[ 57 ],posture[ 58 ],dynamicimages[ 59 60 ],andthedistributionofprimenumbers[ 61 ],amongmanyothers.Itisfurtherobservedthatprinciplecomponentanalysisofsuchprocessesleadstopower-lawdecayingeigenvaluespectrum[ 62 ]. Twoimportantprototypicalmodelsfor1=fprocessesarethefractionalBrownianmotion(fBm)processesandtheON/OFFintermittencywithpower-lawdistributedONandOFFperiods.Below,weapplytheconceptofPSICtobothtypesofprocesses. 1. AllBrownianpathsstartattheorigin:B(0)=0. 2. For0

WemayinferfromtheabovedenitionthattheprobabilitydistributionfunctionofBisgivenby: 2sdu(3{9) Thisfunctionalsosatisesthescalingproperty: Inotherwords,B(t)and1=2B(t)havethesamedistribution.Thus,weseefromEq. 3{1 thatBmisaself-similarprocesswithHurstparameter1/2. SupposewehavemeasuredB(t1);B(t2);t2t1>0.WhatcanwesayaboutB(s);t1

WenotethataBrownianpathisnotadierentiablefunctionoftime.Heuristically,thiscanbeunderstoodinthisway:considerthevariable,B(t+s)B(t),withvariances.Itsstandarddeviation,whichisameasureofitsorderofmagnitude,isp AlthoughB(t)isalmostsurelynotdierentiableint,symbolicallyonestilloftenwrites wherew(t)isstationarywhiteGaussiannoise,andextendstheaboveequationtot<0throughtheconvention,Zt0Z0t 3{12 Toseehow,letuspartition[0;t]intonequallyspacedintervals,t=t=n.Sincew(t)areGaussianrandomvariableswithzeromeanandunitvariance,inorderforB(t)tohavevariancet,weshouldhave Thatis,thecoecientis(t)1=2insteadoft,asonemighthaveguessed.Typically,t1;hence,(t)1=2t.Thus,ifoneincorrectlyusestinsteadof(t)1=2,oneisseverelyunderestimatingthestrengthofthenoise. 45


Whitenoise(a)andBrownianmotion(b).Axesarearbitrary. Whenthetimeisgenuinelydiscrete,ortheunitsoftimecanbearbitrary,onecantakettobe1unit,andEq. 3{13 becomes Eq. 3{14 providesperhapsthesimplestmethodofgeneratingasampleofBm.AnexampleisshowninFig. 3-1 Eq. 3{14 isalsoknownasarandomwalkprocess.Amoresophisticatedrandomwalk(orjump)processisgivenbysummingupaninnitenumberofjumpfunctions:Ji(t)=Ai(tti),where(t)isaunit-stepfunction andAi;tiarerandomvariableswithGaussianandPoissondistributions,respectively. 46


whereWisaGaussianrandomvariableofzeromeanandunitvariance,andHistheHurstparameter.WhenH=1=2,theprocessreducestothestandardBrownianmotion(Bm).ItistrivialtoverifythatthisprocesssatisesthedeningEq. 3{1 foraself-similarstochasticprocess. FractionalBrownianmotionBH(t)isaGaussianprocesswithmean0,stationaryincrements,variance andcovariance: 2ns2H+t2Hjstj2Ho(3{18) whereHistheHurstparameter.DuetoitsGaussiannature,accordingtoEq. 3{2 toEq. 3{4 ,theabovethreepropertiescompletelydetermineitsself-similarcharacter.Fig. 3-2 showsseveralfBmprocesseswithdierentH. Roughly,thedistributionforanfBmprocessis[ 63 ]: (H+1=2)Zt(t)H1=2dB()(3{19) where(t)isthegammafunction:(t)=Z10yt1eydy Inotherwords,thisprocessisinvariantfortransformationsconservingthesimilarityvariablex=tH. 47


SeveralfBmprocesseswithdierentH. TheintegraldenedbyEq. 3{19 diverges.Themoreprecisedenitionis (H+1=2)ZtK(t)dB()(3{21) wherethekernelK(t)isgivenby[ 63 ] TheincrementprocessoffBm,Xi=BH((i+1)t)BH(it);i1,wheretcanbeconsideredasamplingtime,iscalledfractionalGaussiannoise(fGn)process.ItisazeromeanstationaryGaussiantimeseries.NotingthatE(XiXi+k)=Ef[BH((i+1)t)BH(it)][BH((i+1+k)t)BH((i+k)t)]g


3{18 ,oneobtainstheautocovariancefunction(k)forthefGnprocess: 2n(k+1)2H2k2H+jk1j2Ho;k0(3{23) Noticethattheexpressionisindependentoft.Therefore,withoutlossofgenerality,onecouldtaket=1.Inparticular,wehave(1)=1 222H2 (i) WhenH=1=2,(k)=0fork6=0.Thisisthewell-knownpropertyofwhiteGaussiannoise. (ii) When00. Properties(ii)and(iii)areoftentermedanti-persistentandpersistentcorrelations[ 43 ],respectively. Next,weconsiderthebehaviorof(k)whenkislarge.Notingthatwhenjxj1,(1+x)1+x+(1) 2x2 hence, whenH6=1=2.WhenH=1=2,(k)=0fork1,theXi'saresimplywhitenoise. Nowwearereadytocharacterize1=fprocessesintheframeworkofPSIC.Recallingthatthedeningpropertyfora1=fprocessisthatitsvarianceincreaseswithtast2H. 49


Therefore, Expressingtintermsof,wehave ComparingwiththedeningequationofPSIC(Eq. 2{10 ),wendthat 2(3{29) Noticing0

Time-dependentexponent(k)vs.lnkcurvesforfractionalBrownianmotion. ON/OFFsourceswithN=3,X1(t),X2(t),X3(t),andtheirsummationS3(t). 51


Thisexpressionisoftencalledthecomplementarycumulativedistributionfunction(CCDF)ofaheavy-taileddistribution.Infact,itismorepopularforexpressingaheavy-taileddistribution,sinceitemphasizesthetailofthedistribution.Itiseasytoprovethatwhen<2,thevarianceandallhigherthan2nd-ordermomentsdonotexist.Furthermore,when1,themeanalsodiverges. Whenthepower-lawrelationextendstotheentirerangeoftheallowablex,wehavetheParetodistribution: x;xb>0;2>0(3{31) whereandbarecalledtheshapeandthelocationparameters,respectively.Itcanbeproven[ 64 ]thatwhen12,theHurstparameteroftheON/OFFprocesseswithParetodistributedONandOFFperiodsisgivenas OurpurposehereistocheckwhethertheON/OFFprocesseswithParetodistributedONandOFFperiodscanbecharacterizedbyPSIC.InFig. 3-5 wehaveshownafewexamplesfor=1:2;1:6,and2.0.NoticingEq. 3{32 ,weseethattheslopesofthestraightlinesagaincorrectlyestimatetheHparameter.Alsonotethatwhen0<<1,1

Time-dependentexponent(k)vs.lnkcurvesforON/OFFprocesseswithParetodistributedONandOFFperiods. Beforeendingthissection,wewouldliketoemphasizethatthepossibilityofmisinterpreting1=fprocessesbeingdeterministicchaosneveroccursifoneworksundertheframeworkofPSIC.Underthisframework,onemonitorstheevolutionoftwonearbytrajectories.Iftwonearbytrajectoriesdivergeexponentiallyfast,thenthetimeseriesischaotic;ifthedivergenceincreasesinapower-lawmanner,thenthetrajectoriesbelongto1=fprocesses.Theseideascanbeeasilyexpressedprecisely.LetjjViVjjjbetheinitialseparationbetweentwotrajectories.Thisseparationisassumedtobenotlargerthansomeprescribedsmalldistancer.Aftertimet,thedistancebetweenthetwotrajectorieswillbejjVi+tVj+tjj.Fortruechaoticsystems,jjVi+tVj+tjj/jjViVjjjet


65 { 70 ].TherearetwotypesofLevyprocesses[ 71 ].OneisLevyights,whicharerandomprocessesconsistingofmanyindependentsteps,eachstepbeingcharacterizedbyastablelaw,andconsumingaunittimeregardlessofitslength.TheotherisLevywalks,whereeachsteptakestimeproportionaltoitslength.ALevywalkcanbeviewedassampledfromaLevyightwithauniformspeed.TheincrementprocessofaLevywalk,obtainedbydierencingtheLevywalk,isverysimilartoanON/OFFtrainwithpower-lawdistributedONandOFFperiods.Therefore,inthefollowing,weshallnotbefurtherconcernedaboutit.WeshallfocusonLevyights.NotethatLevyights,havingindependentsteps,arememorylessprocessescharacterizedbyH=1=2,irrespectiveofthevalueoftheexponentcharacterizingthestablelaws[ 71 ].Inotherwords,methodssuchasdetrendeductuationanalysis(DFA)[ 72 ]cannotbeusedtoestimatetheparameter. NowwedeneLevyightsmoreprecisely.A(standard)symmetric-stableLevyprocessfL(t);t0gisastochasticprocesswiththefollowingproperties[ 73 ]: 1. 2. 3. 54


whereY1;Y2;areindependentrandomvariables,eachhavingthesamedistributionasY.FromEq. 3{33 ,onethenndsthatnVarY=n2=VarY.Forthedistributiontobevalid,0<2.When=2,thedistributionisGaussian,andhence,thecorrespondingLevyightisjustthestandardBrownianmotion.When0<<2,thedistributionisheavy-tailed,P[Xx]x;x!1,andhasinnitevariance.Furthermore,when0<1,themeanisalsoinnite. Atthispoint,itisworthspendingafewparagraphstoexplainhowtosimulatean-stableLevyprocess.Thekeyistosimulateitsincrementprocess,whichfollowsanstabledistribution.TherstbreakthroughforthesimulationofstabledistributionswasmadebyKanter[ 74 ],whogaveadirectmethodforsimulatingastabledistributionwith<1and=1.TheapproachwaslatergeneralizedtothegeneralcasebyChambersetal.[ 75 ].WerstdescribethealgorithmforconstructingastandardstablerandomvariableXS(;;1). (cos0cos)1=ncos[0+(1)] hasaS(;;1)distribution. Forsymmetricstabledistributionwith=0,theaboveexpressioncanbegreatlysimplied.TosimulateanarbitrarystabledistributionS(;;;),wecansimplytake 55


whereZisgivenbyEq. 3{34 WehavesimulatedanumberofsymmetricstableLevymotionswithdierent.Figure 3-6 showstwoexamplesofLevymotionswith=1:5and1.ThedierencebetweenthestandardBrownianmotionsandLevymotionsliesinthatBrownianmotionsappeareverywheresimilarly,whileLevymotionsarecomprisedofdenseclustersconnectedbyoccasionallongjumps:thesmallertheis,themorefrequentlythelongjumpsappear. LetusnowpauseforawhileandthinkwhereLevyights-likeesotericprocessescouldoccur.Onesituationcouldbethis:amosquitoheadtoagiantspidernetandgotstuck;itstruggledforawhile,andluckily,withthehelpofagustofwind,escaped.Whenthemosquitowasstruggling,thestepsitcouldtakewouldbetiny;butthestepleadingtoitsfortunateescapemustbehuge.Asanotherexample,letusconsider(American)footballgames.Formostofthetimeduringagame,theoenseanddefensewouldbefairlybalanced,andtheoenseteammayonlybeabletoproceedforafewyards.Butduringanattackthatleadstoatouchdown,theherogettingthetouchdownoften\ies"tensofyards|somewhathehasescapedthedefense.Whilethesetwosimplesituationsareonlymeantasanalogies,rememberingthemcouldbehelpfulwhenreadingresearchpaperssearchingforLevystatisticsinvarioustypesofproblemssuchasanimalforaging.Atthispoint,wealsowishtomentionthatLevy-ights-likepatternshavebeenusedasatypeofscreensaverforLinuxoperatingsystems. ThesymmetricstableLevyightis1=self-similar.Thatis,forc>0,theprocessesfL(ct);t0gandfc1=L(t);t0ghavethesamenite-dimensionaldistributions.BythisargumentaswellasEq. 3{33 ,itisclearthatthelengthofthe 56


ExamplesforBrownianmotionsandLevymotions.1-Dand2-DBrownianmotions(a,b)vs.Levymotionswith=1:5and1for(c,d)and(e,f),respectively. 57


Time-dependentexponent(k)vs.lnkcurvesforLevyprocesses. motioninatimespanoft,L(t),isgivenbythefollowingscaling: L(t)/t1=(3{36) ThiscontrastswiththescalingforfractionalBrownianmotion: BH(t)/tH(3{37) Wethusidentifythat1=playstheroleofH.Therefore, Noticing0<2,fromEq. 3{38 ,wehave1q<1. 58


3-7 .Theslopesofthestraightlinescorrectlyestimatethevaluesofusedinthesimulations. 59


Timeseriesanalysisisanimportantexerciseinscienceandengineering.Oneofthemostimportantissuesintimeseriesanalysisistodeterminewhetherthedataunderinvestigationisregular,deterministicallychaotic,orsimplyrandom.Alongthisline,alotofworkhasbeendoneinthepasttwodecades.Abitsurprisingly,quiteoftentheexponentialsensitivitytoinitialconditions(ESIC)intherigorousmathematicalsensecannotbeobservedinexperimentaldata.Whilecurrentconsensusistoattributethisfacttothenoiseinthedata,weshallshowthatthemoregeneralconceptofpower-lawsensitivitytoinitialconditions(PSIC)providesaninterestingframeworkfortimeseriesanalysis.WeillustratethisbystudyingthecomplicateddynamicsofInternettransportprotocols. TheInternetisoneofthemostcomplicatedsystemsthatmanhasevermade.Inrecentyears,twotypesoffascinatingmulti-scalebehaviorshavebeenfoundintheInternet.Oneisthetemporal-domainfractalandmultifractalpropertiesofnetworktrac(see[ 76 77 ]andreferencestherein),whichtypicallyrepresentsaggregationofindividualhoststreams.Theotheristhespatial-domainscale-freetopologyoftheInternet(see[ 78 ]andreferencestherein).Also,ithasbeenobservedthatthefailureofasingleroutermaytriggerroutinginstability,whichmaybesevereenoughtoinstigatearouteapstorm[ 79 ].Furthermore,packetsmaybedeliveredoutoforderorevengetdropped,andpacketreorderingisnotapathologicalnetworkbehavior[ 80 ].AsthenextgenerationInternetapplicationssuchasremoteinstrumentcontrolandcomputationalsteeringarebeingdeveloped,itisbecomingincreasinglyimportanttounderstandthetransportdynamicsinordertosustaintheneededcontrolchannelsoverwide-areaconnections.Typically,thetransportdynamicsovertheInternetconnectionsarearesultofthenonlineardynamicsofthewidely-deployedTransmissionControlProtocol(TCP)interactingwiththeInternettrac,whichisstochasticandoftenself-similar.Thisleadsonetonaturallyexpectthe 60


81 ].Iftransportdynamicsareindeedchaotic,thensteady-stateanalysisandstudyofconvergencetoasteady-statewillnotbethesolefocusinpractice,particularlyforremotecontrolandcomputationalsteeringtasks.Rather,oneshouldalsofocusonthe\transient"non-convergenttransportdynamics.Althoughrecentlyalotofresearchhasbeencarriedout[ 82 { 88 ]tobetterunderstandthisissue,alotoffundamentalproblemsremaintobeanswered.Forexample,theobservationofchaosin[ 84 ]cannotberepeated,whichleadstothesuspicionthatthechaoticmotionsreportedin[ 84 ]maycorrespondtoperiodicmotionswithverylongperiods[ 89 ].Fundamentally,weshallshow,bydevelopinga1-Ddiscretemap,thatthenetworkscenarioswithonlyafewcompetingTCPowsstudiedin[ 84 ]cannotgeneratechaoticdynamics.However,withotherparametersettings,chaosispossible. Inordertoexaminethetemporalevolutionofadynamicalsystem,oneoftenmeasuresascalartimeseriesataxedpointinthestatespace.Takensembeddingtheorem[ 34 { 36 ]ensuresthatthecollectivebehaviorofthedynamicsdescribedbythedynamicalvariablescoupledtothemeasuredonecanbeconvenientlystudiedbythemeasuredscalartimeseries.ToillustratethequalitativeaspectsofTCPdynamicsbyonlymeasuringascalartimeseries,itismosteectivetoconsiderasinglebottlenecklink,as 61


84 ]consistingoflonglivedTCPowsonasinglelink.Thissetuphasfourparameters,thelinkspeedCinMbps,thelinkpropagationdelaydinms,thebuersizeBinunitofpacketsof1000bytes,andthenumberofcompetingTCPowsN.WehavefoundthatNactsasacriticalbifurcationparameter.Forexample,whenC=0:1,d=10ms,B=10andNissmall,wehaveobservedonlyperiodicandquasi-periodicmotions(Figs. 4-1 (a-d)).WhenNislarge,suchasN=17,chaosisobserved(seeFig. 4-2 (a,b)).WethusconcludethatoneroutetochaosinTCPdynamicsisthewell-knownquasi-periodicroute.Inordertoshowmorediversieddynamicalbehaviorsuchasquasi-periodicmotionswithmorethantwoincommensurate(i.e.,independent)frequencies,below,however,weshallvaryalltheparameters.Inparticular,weshallshowthatwhenNissmall,chaoscannothappen. ForeachoftheNcompetingTCPows,onecanrecorditscongestionwindowsizedata.Letusonlyrecordoneofthem,anddenoteitbyW(i).Anexampleofquasi-periodicW(i)isshowninFig. 4-1 (a),whereC=0:1Mbps,delayd=10ms,B=10,N=3.Itspowerspectraldensity(PSD)isshowninFig. 4-1 (b).Weobservemanydiscretesharppeaks.Notethatwhenonlyaround105pointsareusedtocomputethePSD,onedoesnotobservemanypeaksinthespectrum.Ratheroneobservesabroadspectrum,andhence,istemptedtointerprettheW(i)timeseriestobechaotic.Therefore,theW(i)timeseriesisnotideallysuitedforthestudyof(quasi-)periodicmotionswithlongperiods.Thismotivatesustodeneanewtimeseries,T(i),whichisthetimeintervalbetweentheonsetoftwosuccessiveexponentialbacko(i.e.,multiplicativedecrease)ofTCP,asindicatedinFig. 4-1 (a).Thestartofanexponentialbackoindicatestriggeringofalossepisodeorheavycongestion.Hence,theT(i)isrelatedtothewell-knownroundtriptime(RTT).Infact,itcanbeconsideredtobemoreorlessequivalenttothetimeintervalbetweentwosuccessivelossbursts.TheT(i)timeseriesfortheW(i)ofFig. 4-1 (a)isshowninFig. 4-1 (c),togetherwithits 62


Exampleofquasi-periodiccongestionwindowsizeW(i)timeseriessampledbyanequaltimeintervalof10ms(a)anditspowerspectraldensity(PSD)(b).TheparametersusedareC=0:1Mbps,delayd=10ms,B=10packets,whereapacketisofsize1000bytes,andthenumberofTCPstreamsN=3.(c)TheT(i)timeseriescorrespondingto(a)and(d)itsPSD.(e)AnotherT(i)timeseriescorrespondingtoC=0:5Mbps,d=10ms,B=10,N=3,and(f)itsPSD. PSDinFig. 4-1 (d).WeobservethatthisT(i)timeseriesisasimpleperiodicsequence,andhence,theW(i)timeseriesisquasi-periodic.However,T(i)timeseriescanstillbequasi-periodic.AnexampleisshowninFig. 4-1 (e),withC=0:5Mbps,delayd=10ms,B=10packets,N=3,togetherwithitsPSDisFig. 4-1 (f).NotethattherearestillmanydiscretesharppeaksinFig. 4-1 (f).Toappreciatehowmanyindependentfrequencies 63


4-1 (e),wehavefurtherextractedtwonewtimeseries,T(1)(i),whichistheintervaltimeseriesbetweenthesuccessivelocalmaximaoftheT(i)timeseries,andT(2)(i),whichistheintervaltimeseriesbetweenthesuccessivelocalmaximaoftheT(1)(i)timeseries.WehavefoundthattheT(2)(i)issimplyperiodic.Hence,weconcludethattheW(i)timeseriesisquasi-periodicwithfourindependentfrequencies. BeforeweproceedtodiscussthechaoticTCPdynamics,wenotethatsimpleperiodicorquasi-periodicW(i)timeserieswithlessormorethanfourindependentfrequenciescanallbeobserved.Infact,wehavefoundthe\chaotic"congestionwindowsizedatastudiedin[ 84 ]tobequasi-periodicwithtwoindependentfrequencies.Furthermore,theprobabilitydistributionsfortheconstant,periodic,andquasi-periodicT(i)timeseriesareonlycomposedofafewdiscretepeaks.Thissuggeststhattheobservedquasi-periodicT(i)timeseriesconstitutesadiscretetorus. LetusnowturntothediscussionofchaoticTCPdynamics.ShowninFigs. 4-2 (a,c)aretwoirregularT(i)timeseries.TheparametersC,d,andBusedforFig. 4-2 (a)arethesameasthoseforFigs. 4-1 (a-d).NotingfromFig. 4-2 thatT(i)oftencanbeontheorderof104to105,hence,anottoolongW(i)timeseriesisevenmoreill-suitedforthestudyoftheunderlyingirregulardynamics. ToshowthattheT(i)timeseriesinFigs. 4-2 (a,c)isindeedchaotic,weemploythedirectdynamicaltestfordeterministicchaosdevelopedbyGaoandZheng[ 21 22 ].The(k)curvesfortheT(i)timeseriesofFigs. 4-2 (a,c)areshowninFigs. 4-2 (b,d),wherethefourcurves,frombottomtotop,correspondtoshellsofsizes(2(i+1)=2;2i=2),i=14;15;16;17.wehavesimplychosenm=6andL=1.VerysimilarcurveshavebeenobtainedforotherchoicesofmandL.WeobserveaverywelldenedlinearcommonenvelopeatthelowerleftcornerofFigs. 4-2 (b,d).TheexistenceofacommonenvelopeguaranteesthatarobustpositiveLyapunovexponentwillbeobtainednomatterwhichshellisusedinthecomputation.Hence,thetwoT(i)timeseriesareindeedchaotic. 64


TwoexamplesofchaoticT(i)timeseriesobtainedwithparameters(a)C=0:1,delayd=10ms,B=10,andthenumberofTCPstreamsN=17and(c)C=0:1,d=10ms,B=12,andN=19.The(k)curvesfor(a)and(c)areshownin(b)and(d),respectively. WehaveremarkedthattheprobabilitydistributionsfortheT(i)timeseriesofregularmotionsarecomprisedofafewdiscretepeaks.WhatarethedistributionsforthechaoticT(i)timeseriessuchasshowninFigs. 4-2 (a,c)?Theyarepower-law-like,P(Tt)t,foralmosttwoordersofmagnitudeint,asshowninFigs. 4-3 (a,b).Theexponentisnotauniversalquantity,however. Wenowdevelopasimple1-DmaptodescribetheoperationofTCP.TCPreliesontwomechanismstosetitstransmissionrate:owcontrolandcongestioncontrol. 65


Complementarycumulativedistributionfunctions(CCDFs)forthetwochaoticT(i)timeseriesofFigs. 4-2 (a,c). Flowcontrolensuresthatthesendersendsnomorethanthesizeofthereceiver'slastadvertisedow-controlwindowbasedonitsavailablebuersizewhilecongestioncontrolensuresthatthesenderdoesnotunfairlyoverrunthenetwork'savailablebandwidth.TCPimplementsthesemechanismsviaaow-controlwindow(fwnd)advertisedbythereceivertothesenderandacongestion-controlwindow(cwnd)adaptedbasedoninferringthestateofthenetwork.Specically,TCPcalculatesaneectivewindow(ewnd),whereewnd=min(fwnd;cwnd),andthensendsdataatarateewnd=RTT,whereRTTistheroundtriptimeoftheconnection.RecentlyRaoandChua[ 85 ]developedananalyticmodeldescribingTCPdynamics.ByassumingfwndbeingxedandtheTCP'sslowstartonlycontributingtotransientbehavior,RaoandChuahavedevelopedthew-updatemapM:[1;Wmax]![1;Wmax], 66


Theabovemap(Eq. 4{1 )canbesimplied[ 90 ], Sincefwndisnotassumedtobeaconstant,themodiedmapismoregeneralthantheoriginalRaoandChua'smap.ToapplytheabovemaptothestudyofaTCPowcompetingwithN1otherTCPows,onemaylumptheeectoftheN1otherTCPowsasabackgroundtrac.Hence,Wmaxisdeterminedbycwndandthebackgroundtrac,andtypicallyisacomplicatedfunctionoftime,noticingthattheavailablebandwidthofabottlenecklinkcanvarywithtimeconsiderablyduetothedynamicnatureofbackgroundtrac.WehavecarriedoutsimulationstudiesofEq. 4{2 byassumingWmaxtobeperiodicandquasi-periodic,andhaveindeedobservedchaos. WearenowreadytounderstandwhychaoscannotoccurwhenthenumberofcompetingTCPowsissmall,suchas2.Toseethis,letusassumethatwiandwi+wiaswellaswi+1andwi+1+wi+1areallsmallerthanWmax.Itistheneasytoseethatjwi+1jjwij.Thatis,smalldisturbancedecays.ThismeansduringtheadditiveincreasephaseofTCP,nearbytrajectoriescontractinsteadofdiverge.Hence,thedynamicsarestable.InorderforthedynamicstobeunstablesothattheLyapunovexponentispositive,thetransitionfromtheadditiveincreasephasetothemultiplicativedecreasephasehastobefast.ThiscanbetrueonlywhenthenumberofcompetingTCPstreamsislarge.Wethusconcludethatthenetworkscenariosconsideredin[ 84 ]cannotgeneratechaoticTCPdynamics. HowrelevantisoursimulationresulttothedynamicsoftherealInternet?AretheirregularT(i)timeseriesshowninFigs. 4-2 (a,c)anartifactofTahoeversionofTCP 67


ThetimeseriesT(i)extractedfromaW(i)timeseriescollectedusingnet100instrumentsoverORNL-LSUconnection.(a)thetimeseriesT(i),(b)thecomplementarycumulativedistributionfunction(CCDF)fortheT(i)timeseries. weusedormoregenerallyofthens-2simulator?Tondananswer,wecollectedW(i)measurementsbetweenORNLandLouisianaStateUniversityatmillisecondresolutionusingthenet100instrumentsandcomputedT(i)asshowninFig. 4-4 (a),whoseproleisqualitativelyquitesimilartoFig. 4-2 (a).Infact,italsofollowsa(truncated)power-law,asshowninFig. 4-4 (b).Theexponentforthepower-lawpartisevensmallerthanthetimeseriesofFig. 4-2 .Thetruncationinthepower-lawcouldbeduetotheshortnessoftheT(i)timeseries.Aswehavepointedout,theT(i)timeseriesisrelatedtoRTT.CottrellandBullothavebeenexperimentingwithmanyadvancedversionsofTCP,andalsoobservedRTTtimeseriesverysimilartothese.WealsohaveobservedfromRTTandlossdatameasuredongeographicallydispersedpathsontheInternet(witharesolutionof1sec)byresearchersattheSanDiegoSupercomputerCenterthattheprobabilitydistributionsforthetimeintervalbetweensuccessivelossburststypicallyfollowapower-law-likedistribution,withtheexponentalsosmallerthan1.Thus,wehavegoodreasontobelievethattheobservedirregularT(i)timeseriesisnotan 68


91 ],butvaryconsiderablyindurationsdeterminedbyspecicapplications.Inthenextsection,wewillanalyzetheactualdynamicsoftheInternetmeasurements. 92 ]withwhichTCPcompetesforbandwidthandrouterbuers.Roughlyspeaking,theinteractionbetweenthetwoisintermsoftheadditionstoTCPcongestionwindow-size,denotedbycwnd,inresponsetoacknowledgments,andmultiplicativedecreasesinresponsetoinferredlosses(ignoringtheinitialslow-startphase).Anyprotocol,howeversimpleitsdynamicsare,isgenerallyexpectedtoexhibitapparentlycomplicateddynamicsduetoitsinteractionwiththestochasticnetworktrac.TCPinparticularisshown(albeitinsimulation)toexhibitchaoticorchaos-likebehaviorevenwhenthecompetingtracismuchsimplersuchasasinglecompetingTCP[ 84 ]orUserDatagramProtocol(UDP)stream[ 85 ].Thedicultyistounderstandthedynamicswhenboththephenomenaareineect,andequallyimportantlytounderstandtheirimpactonactualInternetstreams.Historically,themethodsfromstandardchaosandstochastictheorieshavebeenunabletooermuchnewperspectiveonthetransportdynamics. 69


21 22 37 ].OurmajorpurposeistoelucidatehowthedeterministicandstochasticcomponentsoftransportdynamicsinteractwitheachotherontheInternet.ItisofconsiderableinteresttonotethatrecentlytherehavebeenseveralimportantworksonTCPdynamicswiththepurposeofimprovingcongestioncontrol[ 93 { 96 ].Infact,therehasbeenconsiderableeortindevelopingnewversionsand/oralternativestoTCPsothatthenetworkdynamicscanbemorestable.Conventionalwaysofanalyzingthenetworkdynamicsareunabletoreadilydeterminewhethernewermethodsresultinstabletransportdynamicsorhelpindesigningsuchmethods.Ouranalysisshowstwoimportantfeaturesinthisdirection: (a) Randomnessisanintegralpartofthetransportdynamicsandmustbeexplicitlyhandled.Inparticular,thestepsizesutilizedforadjustingthecongestionparametersmustbesuitablyvaried(e.g.usingRobbins-Monroconditions)toensuretheeventualconvergence[ 97 ].ThisisnotthecaseforTCPwhichutilizesxedstepsizes. (b) Thechaoticdynamicsofthetransportprotocolsdohaveasignicantimpactonpracticaltransfers,andtheprotocoldesignmustbecognizantofitseects.Inparticular,itmightbeworthwhiletoinvestigateprotocolsthatdonotcontaindominantchaoticregimes,particularlyforremotecontrolandsteeringapplications. Thetraditionaltransportprotocolsarenotdesignedtoexplicitlyaddresstheabovetwoissues,butjustiablysosincetheiroriginalpurposeisdatatransportratherthannercontrolofdynamics. Intherestofthissection,weshallrstbrieydescribethecwnddatastudiedhere.Thenweshallbrieyexplaintheanalysisprocedureanddescribetypicalresultsoftheanalysis. WehavecollectedanumberofcwndtracesusingsingleandtwocompetingTCPstreams.ThisdatawascollectedontwodierentconnectionsfromOakRidgeNationalLaboratory(ORNL)toGeorgiaInstituteofTechnology(GaTech)andtoLouisianaState 70


4-5 (a-d)segmentsoffourdatasetsforORNL-GaTechconnection.Powerspectralanalysisofthesedatadoesnotshowanydominantpeaks,andhence,thedynamicsarenotsimplyoscillatory.SinceourdatawasmeasuredontheInternetwith\live"backgroundtrac,itisapparentlymorecomplicatedandrealisticthanthetracesobtainedbynetworksimulation. Nextletusanalyzecwnddata,x(i);i=1;;n,usingtheconceptsoftimedependentexponent(k)curves[ 21 22 37 ].Ashasbeendiscussedearlier,werstneedtoemploytheembeddingtheorem[ 34 { 36 98 ]toconstructvectorsoftheform:Vi=[x(i);x(i+L);:::;x(i+(m1)L)],wheremistheembeddingdimensionandLthedelaytime.Theembeddingtheorem[ 34 { 36 ]statesthatwhentheembeddingdimensionmislargerthantwicetheboxcountingdimensionoftheattractor,thenthedynamicsoftheoriginalsystemcanbestudiedfromasinglescalartimeseries.Wenotethattheembeddingdimensionusedin[ 84 ]isonlytwo,whichhastobeconsiderednotlargeenough(thismaycallforacloserexaminationoftheirconclusions). Beforewemoveon,wepointoutafewinterestingfeaturesofthe(k)curves: (i) Fornoise,onlyforkuptotheembeddingwindowsize(m1)Lwill(k)increase.Thus,whenever(k)increasesforamuchlargerrangeofk,itisanindicationofnon-trivialdeterministicstructureinthedata. (ii) Forperiodicsignals,(k)isessentiallyzeroforanyk. 71


TimeseriesforthecongestionwindowsizecwndforORNL-GaTechconnection. (iii) Forquasi-periodicsignals,(k)isperiodicwithanamplitudetypicallysmallerthan0.1,hence,forpracticalpurposes,(k)canbeconsideredverycloseto0. (iv) Whenachaoticsignaliscorruptedbynoise,thenthe(k)curvesbreakthemselvesawayfromthecommonenvelope.Thestrongerthenoiseis,themorethe(k)curvesbreakawaytilltheenvelopeisnotdenedatall.Thisfeaturehasactuallybeenusedtoestimatetheamountofbothmeasurementanddynamicnoise[ 99 ]. Nowwearereadytocomputeandunderstandthe(k)curvesforcwndtraces.Thesetof(k)curves,correspondingtoFig. 4-5 areplottedinFigs. 4-6 .Inthecomputations,3104pointsareused,andm=10;L=1.Thesevencurves,fromthebottomto 72


Time-dependentexponent(k)vs.kcurvesforcwnddatacorrespondingtoFig. 4-5 .Inthecomputations,3104pointsareused,andm=10;L=1.Curvesfromthebottomtotopcorrespondtoshellswithsizes(2(i+1)=2;2i=2),i=2;3;;8. top,correspondtoshellsofsizes(2(i+1)=2;2i=2),i=2;3;;8.Wemakethefollowingobservations: (i) Thedynamicsarecomplicatedandcannotbedescribedaseitherperiodicorquasi-periodicmotions,since(k)ismuchlargerthan0. (ii) Thedynamicscannotbecharacterizedaspuredeterministicchaos,sinceinnocasecanweobserveawell-denedlinearenvelope.Thustherandomcomponentofthedynamicsduetocompetingnetworktracisevidentandcannotsimplybeignored. 73


Thedataisnotsimplynoisy,sinceotherwiseweshouldhaveobservedthat(k)isalmostatwhenk>(m1)L.Thus,thedeterministiccomponentofdynamicswhichisduetothetransportprotocolplaysanintegralroleandmustbecarefullystudied.Thefeatures(ii)and(iii)indicatethattheInternettransportdynamicscontainsbothchaoticandstochasticcomponents. (iv) Thereareconsiderabledierencesbetweenthedatawithonly1TCPsourceandwith2competingTCPsources.Inthelattercase,the(k)curvessharplyrisewhenkjustexceedstheembeddingwindowsize,(m1)L.Ontheotherhand,(k)forFig. 4-6 (a)withonlyoneTCPsourceincreasesmuchslowerwhenkjustexceeds(m1)L.Alsoimportantisthatthe(k)curvesinFigs. 4-6 (c,d)aremuchsmootherthanthoseinFigs. 4-6 (a,b).Hence,wecansaythatthedeterministiccomponentofthedynamicsismorevisiblewhentherearemorethan1competingTCPsources(alongthelinesof[ 84 ]). Sincetheincreasingpartofthe(k)curvesarenotverylinear,letusnextexamineifthecwndtracescanbebettercharacterizedbythemoregeneralizedconceptofpower-lawsensitivitytoinitialconditions(PSIC).Figure 4-7 showsthe(k)vs.lnkcurvesforORNL-GaTechconnection.Veryinterestingly,wehavenowindeedobservedabunchofbetterdenedlinearlines,especiallyforsmallscales.Inparticular,thetransportdynamicswithmorethan1competingTCPsourcesarebetterdescribedbytheconceptofPSIC. Tosummarize,byanalyzinganumberofhighqualitycongestionwindow-sizedatameasuredontheInternet,wehavefoundthatthetransportdynamicsarebestdescribedbytheconceptofPSIC,especiallyforsmallscales.Itisinterestingtofurtherexaminehowonemightsuppressthestochasticityofthenetworkbyexecutingmorecontrolsonthenetworkwhenmakingmeasurements,suchasusingalargenumberofcompetingTCPsourcestogetherwithaconstantUDPow.OuranalysismotivatesthatbothchaoticandstochasticaspectsbepaidacloseattentiontoindesigningInternetprotocolsthatarerequiredtoprovidethedesiredandtractabledynamics. 74


Time-dependentexponent(k)vs.lnkcurvesforcwnddatacorrespondingtoFig. 4-5 .Inthecomputations,3104pointsareused,andm=10;L=1.Curvesfromthebottomtotopcorrespondtoshellswithsizes(2(i+1)=2;2i=2),i=2;3;;8. 75


Understandingthenatureofseaclutteriscrucialtothesuccessfulmodelingofseaclutteraswellastofacilitatetargetdetectionwithinseaclutter.Tothisend,animportantquestiontoaskiswhetherseaclutterisstochasticordeterministic.Sincethecomplicatedseacluttersignalsarefunctionsofcomplex(sometimesturbulent)wavemotionsontheseasurface,whilewavemotionsontheseasurfaceclearlyhavetheirowndynamicalfeaturesthatarenotreadilydescribedbysimplestatisticalfeatures,itisveryappealingtounderstandseaclutterbyconsideringsomeoftheirdynamicalfeatures.Inthepastdecade,Haykinetal.havecarriedoutanalysisofsomeseaclutterdatausingchaostheory[ 100 101 ],andconcludedthatseaclutterwasgeneratedbyanunderlyingchaoticprocess.Recently,theirconclusionhasbeenquestionedbyanumberofresearchers[ 102 { 108 ].Inparticular,Unsworthetal.[ 105 106 ]havedemonstratedthatthetwomaininvariantsusedbyHaykinetal.[ 100 101 ],namelythe\maximumlikelihoodofthecorrelationdimensionestimate"andthe\falsenearestneighbors"areproblematicintheanalysisofmeasuredseaclutterdata,sincebothinvariantsmayinterpretstochasticprocessesaschaos.Theyhavealsotriedanimprovedmethod,whichisbasedonthecorrelationintegralofGrassbergerandProcaccia[ 109 ]andhasbeenfoundeectiveindistinguishingstochasticprocessesfromchaos.Still,noevidenceofdeterminismorchaoshasbeenfoundinseaclutterdata. Toreconcileevergrowingevidenceofstochasticityinseaclutterwiththeirchaoshypothesis,recently,Haykinetal.[ 110 ]havesuggestedthatthenon-chaoticfeatureofseacluttercouldbeduetomanytypesofnoisesourcesinthedata.Totestthispossibility,McDonaldandDamini[ 111 ]havetriedaseriesoflow-passlterstoremovenoise;butagaintheyhavefailedtondanychaoticfeatures.Furthermore,theyhavefoundthatthecommonlyusedchaoticinvariantmeasuresofcorrelationdimensionandLyapunovexponent,computedbyconventionalways,producesimilarresultsformeasuredsea 76


Whiletheserecentstudieshighlysuggestthatseaclutterisunlikelytobetrulychaotic,anumberoffundamentalquestionsarestillunknown.Forexample,mostofthestudiesin[ 102 { 106 ]areconductedbycomparingmeasuredseaclutterdatawithsimulatedstochasticprocesses.Wecanask:canthenon-chaoticnatureofseaclutterbedirectlydemonstratedwithoutresortingtosimulatedstochasticprocesses?Recognizingthatsimplelow-passlteringdoesnotcorrespondtoanydenitescalesinphasespace,canwedesignamoreeectivemethodtoseparatescalesinphasespaceandtotestwhetherseacluttercanbedecomposedassignalsplusnoise?Finally,willstudiesalongthislinebeofanyhelpfortargetdetectionwithinseaclutter? Inthischapter,weemploythedirectdynamicaltestfordeterministicchaosdiscussedinSec. 2.1 toanalyze280seaclutterdatameasuredundervariousseaandweatherconditions.Themethodoersamorestringentcriterionfordetectinglow-dimensionalchaos,andcansimultaneouslymonitormotionsinphasespaceatdierentscales.However,nochaoticfeatureisobservedfromanyofthesedataatallthedierentscalesexamined.Butinterestingly,weshowthatscale-dependentexponentcorrespondingtolargescaleappearstobeusefulfordistinguishingseaclutterdatawithandwithouttargets.Thissuggeststhatseacluttermaycontaininterestingdynamicfeaturesandthatthescale-dependentexponentmaybeanimportantparameterfortargetdetectionwithinseaclutter.Furthermore,wendthatseacluttercanbeconvenientlycharacterizedbythenewconceptofpower-lawsensitivitytoinitialconditions(PSIC).Weshowthatforthepurposeofdetectingtargetswithinseaclutter,PSICoersamoreeectiveframework. Below,weshallrstbrieydescribetheseaclutterdata.Thenwestudytheseaclutterdatabyemployingthedirectdynamicaltestforlow-dimensionalchaosdevelopedbyGaoandZheng[ 21 22 ].Andweshowthatthescale-dependentexponentcorrespondingtolargescalecanbeusedtoeectivelydetectlowobservabletargetswithin 77


5-1 .Thedistancebetweentwoadjacentrangebinswas15m.Oneorafewrangebins(say,Bi1,BiandBi+1)hitatarget,whichwasasphericalblockofstyrofoamofdiameter1m,wrappedwithwiremesh.Thelocationsofthethreetargetswerespeciedbytheirazimuthalangleanddistancetotheradar.Theywere(128degree,2660m),(130degree,5525m),and(170degree,2655m),respectively.Therangebinwherethetargetisstrongestislabeledastheprimarytargetbin.Duetodriftofthetarget,binsadjacenttotheprimarytargetbinmayalsohitthetarget.Theyarecalledsecondarytargetbins.Foreachrangebin,therewere217complexnumbers,sampledwithafrequencyof1000Hz.Amplitudedataoftwopolarizations,HH(horizontaltransmission,horizontalreception)andVV(verticaltransmission,verticalreception)wereanalyzed. Fig. 5-2 showstwoexamplesofthetypicalseaclutteramplitudedatawithoutandwithtarget.However,carefulexaminationoftheamplitudedataindicatesthat4datasetsareseverelyaectedbyclipping.ThiscanbereadilyobservedfromFigs. 5-3 (a,b).We 78


Collectionofseaclutterdata Figure5-2. Typicalseaclutteramplitudedata(a)withoutand(b)withtarget. discardthose4datasets,andanalyzetheremaining10measurements,whichcontain280seacluttertimeseries. 21 22 ],toanalyzeseaclutterdata.Themethodoersamorestringentcriterionforlow-dimensionalchaos,andcansimultaneouslymonitormotionsinphasespaceatdierentscales.Themethodhasfoundnumerousapplicationsinthe 79


Twoshortsegmentsoftheamplitudeseaclutterdataseverelyaectedbyclipping. studyoftheeectsofnoiseondynamicalsystems[ 23 24 26 27 ],estimationofthestrengthofmeasurementnoiseinexperimentaldata[ 28 29 ],pathologicaltremors[ 30 ],shear-thickeningsurfactantsolutions[ 31 ],diluteshearedaqueoussolutions[ 32 ],andserratedplasticows[ 33 ].Inparticular,thismethodwasusedbyGaoetal.[ 107 108 ]toanalyzeonesinglesetofseaclutterdata.Whilechaoswasnotobservedfromthatdataset,nodenitegeneralconclusioncouldbereached,duetolackoflargeamountofdataatthattime.Herewesystematicallystudy280seaclutterdatameasuredundervariousseaandweatherconditions,andexaminewhetheranychaoticfeaturescanbefoundfromtheseseaclutterdata. TheexplicitincorporationofscalesintheGaoandZheng'stest[ 21 22 ]enablesustosimultaneouslymonitormotionsinphasespaceatdierentscales.Wehavesystematicallyanalyzed280amplitudeseacluttertimeseriesmeasuredundervariousseaandweatherconditions.However,nochaoticfeaturehasbeenobservedfromanyofthesedataatallthescalesexamined.Typicalexamplesofthe(k)vs.kcurvesfortheseaclutteramplitudedatawithouttargetsareshowninFigs. 5-4 (a,c,e)andthecurvesforthedatawithtargetsshowninFigs. 5-4 (b,d,f),respectively.Wehavesimplychosenm=6and 80


Examplesofthetime-dependentexponent(k)vs.kcurvesfortheseaclutterdata(a,c,e)withoutand(b,d,f)withthetarget.Sixcurves,frombottomtotop,correspondtoshells(2(i+1)=2;2i=2)withi=13,14,15,16,17,and18.Thesamplingtimefortheseaclutterdatais1msec,andembeddingparametersarem=6,L=1.217datapointsareusedinthecomputation. 5-4 aregenericamongallthe280seaclutterdataanalyzedhere.Hence,wehavetoconcludethatnoneoftheseaclutterdataischaotic. 81


112 ],log-normal[ 113 ],K[ 114 115 ],andcompound-Gaussiandistributions[ 116 ],aswellasusingchaostheory[ 100 { 104 107 108 110 111 ],wavelets[ 117 ],neuralnetworks[ 118 119 ],wavelet-neuralnetcombinedapproaches[ 120 121 ],andtheconceptoffractaldimension[ 122 ]andfractalerror[ 123 124 ].However,nosimplemethodhasbeenfoundtorobustlydetectlowobservableobjectswithinseaclutter[ 125 ]. Inthissection,weexaminewhethertheLyapunovexponentestimatedbyconventionalmethodsandthescale-dependentexponentmaybehelpfulfortargetdetectionwithinseaclutter. WerstcheckwhethertheLyapunovexponentestimatedbyconventionalmethodscanbeusedfordistinguishingseaclutterdatawithandwithouttargets.Aspointedoutearlier,theconventionalwayofestimatingtheLyapunovexponentistocompute(k)=kt,where(k)isdenedasinEq. 2{4 ,subjecttotheconditionsthatkXiXjk

ThemodiedapproachforestimatingtheLyapunovexponentusestheleast-squaresttotherstfewsamplesofthe(k)curvewherethe(k)curveincreaseslinearlyorquasi-linearly.Thisisdoneforallthe280seacluttertimeseries.TobetterappreciatethevariationoftheLyapunovexponentamongthe14rangebinsoftheseaclutterdata,wehaverstsubtractedtheparameterofeachbinbytheminimumofvaluesforthatmeasurement,andthennormalizedtheobtainedvaluesbyitsmaximum.Thevariationsoftheparameterswiththe14rangebinsforthe10HHmeasurementsareshowninFigs. 5-5 (a)-(j),respectively,whereopencirclesdenotetherangebinswiththetarget,andasteroidsdenotethebinswithoutthetarget.Theprimarytargetbinisexplicitlyindicatedbyanarrow.Similarresultshavebeenobtainedforthe10VVmeasurements.WhileFigs. 5-5 (a)and(h)indicatethattheprimarytargetbincanbeseparatedfrombinswithoutthetarget,ingeneral,wehavetoconcludethattheLyapunovexponentestimatedbymethodsequivalentorsimilartoconventionalmeanscannotbeusedfordistinguishingseaclutterdatawithandwithouttargets. Nextletusexaminewhetherthescale-dependentexponentcorrespondingtolargescalemaybeusefulfordetectingtargetswithinseaclutterdata.Byscale-dependentexponent,wemeanthatifweusetheleast-squaresttothelinearlyorquasi-linearlyincreasingpartofthe(k)curvesofdierentshells,theslopesofthoselinesdependonwhichshellsareusedforcomputingthe(k)curves.Inotherwords,ifweplottheexponentestimatedthiswayagainstthesizeoftheshell,thentheexponentvarieswiththeshellsizeorthescale.Thescale-dependentexponenthasbeenusedinthestudyofnoise-inducedchaosinanopticallyinjectedsemiconductorlasermodeltoexaminehownoiseaectsdierentscalesofdynamicsystems[ 26 ].Wenowfocusontheexponentcorrespondingtolargeshellorscale.Figures 5-6 (a)-(j)showthevariationsofthescale-dependentexponentcorrespondingtolargescalewiththe14rangebinsforthe 83


VariationsoftheLyapunovexponentestimatedbyconventionalmethodsvs.the14rangebinsforthe10HHmeasurements.Opencirclesdenotetherangebinswithtarget,while*denotethebinswithouttarget.Theprimarytargetbinisexplicitlyindicatedbyanarrow. same10HHmeasurementsasstudiedinFig. 5-5 .Weobservethattheprimarytargetbincanbeeasilyseparatedfromtherangebinswithoutthetarget,sincethescale-dependentexponentfortheprimarytargetbinismuchlargerthanthoseforthebinswithoutthetarget. Itisthusclearthatthescale-dependentexponentcorrespondingtolargescaleisveryusefulfordistinguishingseaclutterdatawithandwithouttargets.Thissuggeststhatseacluttermaycontaininterestingdynamicfeaturesandthatthescale-dependentexponentmaybeanimportantparameterfortargetdetectionwithinseaclutter.This 84


Variationsofthescale-dependentexponentcorrespondingtolargescalevs.therangebinsforthe10HHmeasurements.Opencirclesdenotetherangebinswithtarget,while*denotethebinswithouttarget.Theprimarytargetbinisexplicitlyindicatedbyanarrow. exampleclearlysigniestheimportanceofincorporatingtheconceptofscaleinameasure.Infact,theconceptofscaleisonlyincorporatedinthetimedependentexponent(k)curves[ 21 22 ]inastaticmanner.InChapter 6 ,weshallseethatwhenameasuredynamicallyincorporatestheconceptofscale,itwillbecomemuchmorepowerful. Notethatonedicultyofusingthescale-dependentexponentfortargetdetectionwithinseaclutteristhatwemayneedtochooseasuitablescaleforestimatingtheexponent.Thismakesthemethodnoteasytouse,anditishardtomakethemethodautomatic. 85


5-7 (a)and(b),respectively,wherethecurvesdenotedbyasteroidsarefordatawithoutthetarget,whilethecurvesdenotedbyopencirclesarefordatawiththetarget.Weobservethatthecurvesarefairlylinearfortherstafewsamples.Alsonoticethattheslopesofthecurvesforthedatawiththetargetaremuchlargerthanthoseforthedatawithoutthetarget.Forconvenience,wedenotetheslopeofthecurvebytheparameter.Tobetterappreciatethevariationoftheparameterwiththerangebins,wenormalizeofeachbinbythemaximalvalueofthe14rangebinswithinthesinglemeasurement.Fig. 5-8 showsthevariationoftheparameterwiththe14rangebins.Itisclearthattheparametercanbeusedtodistinguishseaclutterdatawithandwithouttarget.Interestingly,thefeatureshowninFig. 5-8 isgenericallytrueforallthemeasurements.Itisworthpointingoutthatusuallythe(k)vs.lnkcurvesfordierentscalesarealmostparallel,thustheestimatedparameterisrelativelylessdependentonthescale.Thisisaverynicefeatureofthismethod. Letusexamineifarobustdetectorfordetectingtargetswithinseacluttercanbedevelopedbasedontheparameter.Wehavesystematicallystudied280timeseriesoftheseaclutterdatameasuredundervariousseaandweatherconditions.Tobetterappreciatethedetectionperformance,wehaverstonlyfocusedonbinswithprimarytargets,butomittedthosewithsecondarytargets,sincesometimesitishardtodeterminewhetherabinwithsecondarytargetreallyhitsatargetornot.Afteromittingtherangebindatawithsecondarytargets,thefrequenciesfortheparameterunderthetwohypotheses(thebinswithouttargetsandthosewithprimarytargets)forHHdatasets 86


Examplesofthe(k)vs.lnkcurvesforrangebins(a)withoutand(b)withtarget.Opencirclesdenotetherangebinswithtarget,while*denotethebinswithouttarget. Figure5-8. Variationoftheparameterwiththe14rangebins.Opencirclesdenotetherangebinswithtarget,while*denotethebinswithouttarget. areshowninFig. 5-9 .WeobservethatthefrequenciescompletelyseparatefortheHHdatasets.Thismeansthedetectionaccuracycanbe100%. 87


FrequenciesofthebinswithouttargetsandthebinswithprimarytargetsfortheHHdatasets. 88


InChapter 1 ,wehaveemphasizedtheimportanceofdevelopingscale-dependentmeasurestosimultaneouslycharacterizebehaviorsofcomplexmultiscaledsignalsonawiderangeofscales.Inthischapter,weshalldevelopaneectivealgorithmtocomputeanexcellentmultiscalemeasure,thescale-dependentLyapunovexponent(SDLE),andshowthattheSDLEcanreadilyclassifyvarioustypesofcomplexmotions,includingtrulylow-dimensionalchaos,noisychaos,noise-inducedchaos,processesdenedbypower-lawsensitivitytoinitialconditions(PSIC),andcomplexmotionswithchaoticbehavioronsmallscalesbutdiusivebehavioronlargescales.Finally,weshalldiscusshowtheSDLEcanhelpdetecthiddenfrequenciesinlargescaleorderlymotions. 126 ].ThealgorithmforcalculatingtheFSLEisverysimilartotheWolfetal'salgorithm[ 127 ].Itcomputestheaverager-foldtimebymonitoringthedivergencebetweenareferencetrajectoryandaperturbedtrajectory.Todoso,itneedstodene\nearestneighbors",aswellasneedstoperform,fromtimetotime,arenormalizationwhenthedistancebetweenthereferenceandtheperturbedtrajectorybecomestoolarge.Suchaprocedurerequiresverylongtimeseries,andtherefore,isnotpractical.Tofacilitatederivationofafastalgorithmthatworksonshortdata,aswellastoeasediscussionofcontinuousbutnon-dierentiablestochasticprocesses,wecastthedenitionoftheSDLEasfollows. Consideranensembleoftrajectories.Denotetheinitialseparationbetweentwonearbytrajectoriesby0,andtheiraverageseparationattimetandt+tbytandt+t,respectively.Beingdenedinanaveragesense,tandt+tcanbereadilycomputedfromanyprocesses,eveniftheyarenon-dierentiable.Nextweexaminetherelationbetweent


where(t)istheSDLE.Itisgivenby Givenatimeseriesdata,thesmallesttpossibleisthesamplingtime. ThedenitionoftheSDLEsuggestsasimpleensembleaveragebasedschemetocomputeit.Astraightforwardwaywouldbetondallthepairsofvectorsinthephasespacewiththeirdistanceapproximately,andthencalculatetheiraveragedistanceafteratimet.Thersthalfofthisdescriptionamountstointroducingashell(indexedask), whereVi;Vjarereconstructedvectors,k(theradiusoftheshell)andk(thewidthoftheshell)arearbitrarilychosensmalldistances.Suchashellmaybeconsideredasadierentialelementthatwouldfacilitatecomputationofconditionalprobability.Toexpeditecomputation,itisadvantageoustointroduceasequenceofshells,k=1;2;3;.Notethatthiscomputationalprocedureissimilartothatforcomputingtheso-calledtime-dependentexponent(TDE)curves[ 21 22 37 ]. Withalltheseshells,wecanthenmonitortheevolutionofallofthepairsofvectors(Vi;Vj)withinashellandtakeaverage.Wheneachshellisverythin,byassumingthattheorderofaveragingandtakinglogarithminEq. 6{2 canbeinterchanged,wehave wheretandtareintegersinunitofthesamplingtime,andtheanglebracketsdenoteaveragewithinashell.NotethatcontributionstotheSDLEataspecicscalefromdierentshellscanbecombined,withtheweightforeachshellbeingdeterminedbythe 90


Intheaboveformulation,itisimplicitlyassumedthattheinitialseparation,kViVjk,alignswiththemostunstabledirectioninstantly.Forhigh-dimensionalsystems,thisisnottrue,especiallywhenthegrowthrateisnon-uniformand/ortheeigenvectorsoftheJacobianarenon-normal.Fortunately,theproblemisnotasseriousasonemightbeconcerned,sinceourshellsarenotinnitesimal.WhencomputingtheTDEcurves[ 21 22 37 ],wehavefoundthatwhendicultiesarise,itisoftensucienttointroduceanadditionalcondition, whenndingpairsofvectorswithineachshell.SuchaschemealsoworkswellwhencomputingtheSDLE.Thismeansthat,aftertakingatimecomparabletotheembeddingwindow(m1)L,itwouldbesafetoassumethattheinitialseparationhasevolvedtothemostunstabledirectionofthemotion. Beforeproceedingon,wewishtoemphasizethemajordierencebetweenouralgorithmandthestandardmethodforcalculatingtheFSLE.Aswehavepointedout,tocomputetheFSLE,twotrajectories,oneasreference,anotherasperturbed,havetobedened.Thisrequireshugeamountsofdata.Ouralgorithmavoidsthisbyemployingtwocriticaloperationstofullyutilizeinformationaboutthetimeevolutionofthedata:(i)Thereferenceandtheperturbedtrajectoriesarereplacedbytimeevolutionofallpairsofvectorssatisfyingtheinequality( 6{5 )andfallingwithinashell,and(ii)introductionofasequenceofshellsensuresthatthenumberofpairsofvectorswithintheshellsislargewhiletheensembleaveragewithineachshelliswelldened.LetthenumberofpointsneededtocomputetheFSLEbystandardmethodsbeN.Thesetwooperationsimplythatthemethoddescribedhereonlyneedsaboutp 91


2{5 .Figure 6-1 (a)showsvecurves,forthecasesofD=0;1;2;3;4.Thecomputationsaredonewith10000pointsandm=4;L=2.Weobserveafewinterestingfeatures: Nextweconsidernoise-inducedchaos.Toillustratetheidea,wefollow[ 23 ]andstudythenoisylogisticmap 92


Scale-dependentLyapunovexponent()curvesfor(a)thecleanandthenoisyLorenzsystem,and(b)thenoise-inducedchaosinthelogisticmap.Curvesfromdierentshellsaredesignatedbydierentsymbols. whereisthebifurcationparameter,andPnisaGaussianrandomvariablewithzeromeanandstandarddeviation.In[ 23 ],wereportedthatat=3:74and=0:002,noise-inducedchaosoccurs,andthoughtthatitmaybediculttodistinguishnoise-inducedchaosfromcleanchaos.InFig. 6-1 (b),wehaveplottedthe(t)forthisparticularnoise-inducedchaos.Thecomputationwasdonewithm=4;L=1and10000points.WeobservethatFig. 6-1 (b)isverysimilartothecurvesofnoisychaosplottedinFig. 6-1 (a).Hence,noise-inducedchaosissimilartonoisychaos,butdierentfromcleanchaos. Atthispoint,itisworthmakingtwocomments:(i)Onverysmallscales,theeectofmeasurementnoiseissimilartothatofdynamicnoise.(ii)The()curvesshowninFig. 6-1 arebasedonafairlysmallshell.ThecurvescomputedbasedonlargershellscollapseontherightpartofthecurvesshowninFig. 6-1 .Becauseofthis,forchaoticsystems,oneorafewsmallshellswouldbesucient.Ifonewishestoknowthebehaviorofoneversmallerscales,onehastouselongerandlongertimeseries. 93


Scale-dependentLyapunovexponent()curvefortheMackey-Glasssystem.Thecomputationwasdonewithm=5;L=1,and5000pointssampledwithatimeintervalof6. Finally,weconsidertheMackey-Glassdelaydierentialsystem[ 39 ],denedbyEq. 2{6 .ThesystemhastwopositiveLyapunovexponents,withthelargestLyapunovexponentcloseto0.007[ 21 ].HavingtwopositiveLyapunovexponentswhilethevalueofthelargestLyapunovexponentofthesystemisnotmuchgreaterthan0,onemightbeconcernedthatitmaybediculttocomputetheSDLEofthesystem.Thisisnotthecase.Infact,thissystemcanbeanalyzedasstraightforwardlyasotherdynamicalsystemsincludingtheHenonmapandtheRosslersystem.Anexampleofthe()curveisshowninFig. 6-2 ,wherewehavefollowed[ 21 ]andusedm=5;L=1,and5000pointssampledwithatimeintervalof6.Clearly,weobserveawelldenedplateau,withitsvaluecloseto0.007.ThisexampleillustratesthatwhencomputingtheSDLE,onedoesnotneedtobeveryconcernedaboutnon-uniformgrowthrateinhigh-dimensionalsystems. Uptillnow,wehavefocusedonthepositiveportionof(t).Itturnsoutthatwhentislarge,()becomesoscillatory,withmeanabout0.Denotethecorrespondingscalesby 94


3 ,wehaveseenthatpower-lawsensitivitytoinitialconditions(PSIC)providesacommonfoundationforchaostheoryandrandomfractaltheory.Here,weexaminewhethertheSDLEisabletocharacterizeprocessesdenedbyPSIC.Pleasingly,thisisindeedso.Below,wederiveasimpleequationrelatingtheqandqofPSICtotheSDLE. First,werecallthatPSICisdenedbyt=limjjV(0)jj!0jjV(t)jj jjV(0)jj=1+(1q)(1)qt1=(1q) 6{2 ,wendthat Whent!0,1+(1q)qt(1q)qt.SimplifyingEq. 6{7 ,weobtain Wenowconsiderthreecases: 95


3{28 )and(1)q=H(Eq. 3{29 ).Therefore,(t)=H1 3{38 )and(1)q=1=(Eq. 3{39 ).Therefore,(t)=1 128 ]andstudythefollowingmap, where[xn]denotestheintegerpartofxn,tisanoiseuniformlydistributedintheinterval[1;1],isaparameterquantifyingthestrengthofnoise,andF(y)isgivenby ThemapF(y)isshowninFig. 6-3 asthedashedlines.ItgivesachaoticdynamicswithapositiveLyapunovexponentln(2+).Ontheotherhand,theterm[xn]introducesarandomwalkonintegergrids. Itturnsoutthissystemisveryeasytoanalyze.When=0:4,withonly5000pointsandm=2;L=1,wecanresolveboththechaoticbehavioronverysmallscales,andthenormaldiusivebehavior(withslope2)onlargescales.SeeFig. 6-4 (a). Wenowaskaquestion:Givenasmalldataset,whichtypeofbehavior,thechaoticorthediusive,isresolvedrst?Toanswerthis,wehavetriedtocomputetheSDLEwithonly500points.TheresultisshowninFig. 6-4 (b).Itisinterestingtoobservethatthechaoticbehaviorcanbewellresolvedbyonlyafewhundredpoints.However,the 96


ThefunctionF(x)(Eq. 6{10 )for=0:4isshownasthedashedlines.ThefunctionG(x)(Eq. 6{11 )isanapproximationofF(x),obtainedusing40intervalsofslope0.Inthecaseofnoise-inducedchaosdiscussedinthepaper,G(x)isobtainedfromF(x)using104intervalsofslope0.9. Scale-dependentLyapunovexponent()forthemodeldescribedbyEq. 6{9 .(a)5000pointswereused;forthenoisycase,=0:001.(b)500pointswereused. 97


Wehavealsostudiedthenoisymap.TheresultingSDLEfor=0:001isshowninFig. 6-4 (a),assquares.Wehaveagainused5000points.WhilethebehavioroftheSDLEsuggestsnoisydynamics,with5000points,wearenotabletowellresolvetherelation()ln.Thisindicatesthatforthenoisymap,onverysmallscales,thedimensionisveryhigh. Map( 6{9 )canbemodiedtogiverisetoaninterestingsystemwithnoise-inducedchaos.ThiscanbedonebyreplacingthefunctionF(y)inmap( 6{9 )byG(y)toobtainthefollowingmap[ 128 ], wheretisanoiseuniformlydistributedintheinterval[1;1],isaparameterquantifyingthestrengthofnoise,andG(y)isapiecewiselinearfunctionwhichapproximatesF(y)ofEq. 6{10 .AnexampleofG(y)isshowninFig. 6-3 .Inournumericalsimulations,wehavefollowedCencinietal.[ 128 ]andused104intervalsofslope0.9toobtainG(y).WithsuchachoiceofG(y),intheabsenceofnoise,thetimeevolutiondescribedbythemap( 6{11 )isnonchaotic,sincethelargestLyapunovexponentln(0:9)isnegative.Withappropriatenoiselevel(e.g.,=104or103),theSDLEforthesystembecomesindistinguishabletothenoisySDLEshowninFig. 6-4 forthemap( 6{9 ).Havingadiusiveregimeonlargescales,thisisamorecomplicatednoise-inducedchaosthantheonewehavefoundfromthelogisticmap. Beforeproceedingon,wemakeacommentonthecomputationoftheSDLEfromdeterministicsystemswithnegativelargestLyapunovexponents,suchasthemap( 6{11 )withoutnoise.Atransient-freetimeseriesfromsuchsystemsisaconstanttimeseries.Therefore,thereisnoneedtocomputetheSDLEorothermetrics.Whenthetimeseries 98


2{5 .InFigs. 6-5 (a-c),wehaveshownthepower-spectraldensity(PSD)ofthex;y;zcomponentsofthesystem.WeobservethatthePSDofx(t)andy(t)aresimplybroad.However,thePSDofz(t)showsupaverysharpspectralpeak.Recallthatgeometrically,theLorenzattractorconsistsoftwoscrolls(seeFig. 6-6 ).ThesharpspectralpeakinthePSDofz(t)oftheLorenzsystemisduetothecircularmotionsalongeitherofthescrolls. TheaboveexampleillustratesthatifthedynamicsofasystemcontainsahiddenfrequencythatcannotberevealedbytheFouriertransformofameasuredvariable(say,x(t)),theninordertorevealthehiddenfrequency,onehastoembedx(t)toasuitablephasespace.Thisideahasledtothedevelopmentoftwointerestingmethodsforidentifyinghiddenfrequencies.OnemethodisproposedbyOrtega[ 129 130 ],bycomputingthetemporalevolutionofdensitymeasuresinthereconstructedphasespace.AnotherisproposedbyChernetal.[ 131 ],bytakingsingularvaluedecompositionoflocalneighbors.TheOrtega'smethodhasbeenappliedtoanexperimentaltimeseriesrecorded 99

PAGE 100

Power-spectraldensity(PSD)for(a)x(t),(b)y(t),(c)z(t),(d)(x)1(t),(e)(y)1(t),and(f)(z)1(t). fromafar-infraredlaserinachaoticstate.Thelaserdatasetcanbedownloadedfromthelink 6-7 (a)showsthelaserdataset.ThePSDofthedataisshowninFig. 6-7 (b),whereoneobservesasharppeakaround1.7MHz.Fig. 6-7 (c)showsthePSDofthedensitytimeseries,whereonenotesanadditionalspectralpeakaround37kHz.Thispeakisduetotheenvelopeofchaoticpulsations,whichisdiscernablefromFig. 6-7 (a). 100

PAGE 101

Lorenzattractor. ToappreciatethestrengthoftheadditionalspectralpeakinFig. 6-7 (c),wehavetonotethattheunitsofthePSDinbothFigs. 6-7 (b,c)arearbitrary|thelargestpeakistreatedas1unit,aswasdonebyOrtega.Infact,theactualenergyoftheoriginallaserintensitytimeseriesismuchgreaterthanthatofthedensitytimeseries.Therefore,theadditionalspectralpeakof37kHzin(c)isreallyquiteweak.Inotherwords,althoughtheOrtega'smethodisabletorevealhiddenfrequencies,itisnotveryeective.Infact,Chernetal.havepointedoutthattheeectivenessoftheirmethodaswellastheOrtega'smethodtoexperimentaldataanalysisremainsuncertain. Oneofthemajordicultiesoftheabovetwomethodsisconceptual|bothmethodsarebasedonverysmallscaleneighbors,whilethephenomenonitselfislargescale.Recognizingthis,onecanreadilyunderstandthattheconceptoflimitingscale,developedinSec. 6.2.1 ,providesanexcellentsolutiontotheproblem.Inotherwords,onecansimplytaketheFouriertransformofthetemporalevolutionofthelimitingscale,andexpectadditionalwell-denedspectralpeaks,ifthedynamicscontainhiddenfrequencies.To 101

PAGE 102

Hiddenfrequencyphenomenonoflaserintensitydata.(a)laserintensitydata,(b)thepowerspectraldensity(PSD)fortheoriginallaserintensitytimeseries,(c)thePSDforthedensitytimeseriesofthereconstructedphasespacefromthelaserdata,and(d)thePSDforthelimitingscaleofthelaserdata.In(c)and(d),theembeddingparametersarem=4;L=1. 102

PAGE 103

6-5 (d-f)thePSDofthelimitingscalesofx(t);y(t)andz(t)oftheLorenzsystem.Weobserveverywell-denedspectralpeaksinallthreecases.Infact,nowthezcomponentnolongerplaysamorespecialrolethanthexandycomponents.Themethodissimilarlyamazinglyeectiveinidentifyingthehiddenfrequencyfromthelaserdata,asshowninFig. 6-7 (d):wenownotonlyobserveanadditionalspectralpeakaround37KHzin(d),butalsothatthispeakisevenmoredominantthanthepeakaround1.7MHz.Thereasonfortheexchangeofthisdominanceisduetofairlylargesamplingtime|although80nsissmall,itisonlyabletosampleabout8pointsineachoscillation.Whentheembeddingwindow,(m1)L,becomesgreaterthan4,theNyquistsamplingtheoremisviolatedinthereconstructedphasespace,thenthispeakmayevendiminish. Itisimportanttonotethatthelimitingscaleisnotaectedmuchbyeithermeasurementnoiseordynamicnoise.Therefore,onecanexpectthatthehiddenfrequenciesrevealedbylimitingscaleswillnotbeaectedmuchbynoiseeither. 103

PAGE 104

PAGE 105

[1] P.Fairley,\Theunrulypowergrid,"IEEESpectrum,vol.41,pp.22-27,2004. [2] C.Tsallis,\PossiblegeneralizationofBoltzmann-Gibbsstatistics,"JStatPhys,vol.52,pp.479-487,1988. [3] A.R.PlastinoandA.Plastino,\Stellarpolytropesandtsallisentropy,"Phys.Lett.A,vol.174,pp.384-386,1993. [4] A.Lavagno,G.Kaniadakis,M.Rego-Monteiro,P.Quarati,andC.Tsallis,\Non-extensivethermostatisticalapproachofthepeculiarvelocityfunctionofgalaxyclusters,"Astrophys.Lett.Commun.,vol.35,pp.449-455,1998. [5] M.AntoniandS.Ruo,\ClusteringandrelaxationinHamiltonianlong-rangedynamics,"Phys.Rev.E,vol.52,pp.2361-2374,1995. [6] V.Latora,A.Rapisarda,andC.Tsallis,\Non-Gaussianequilibriuminalong-rangeHamiltoniansystem,"Phys.Rev.E,vol.64,pp.056134,2001. [7] M.L.LyraandC.Tsallis,\Nonextensivityandmultifractalityinlow-dimensionaldissipativesystems,"Phys.Rev.Lett.,vol.80,pp.53-56,1998. [8] T.ArimitsuandN.Arimitsu,\Tsallisstatisticsandfullydevelopedturbulence,"J.Phys.A,vol.33,pp.L235-L241,2000. [9] C.Beck,\Applicationofgeneralizedthermostatisticstofullydevelopedturbulence,"PhysicaA,vol.277,pp.115-123,2000. [10] C.Beck,G.S.Lewis,andH.L.Swinney,\MeasuringnonextensitivityparametersinaturbulentCouette-Taylorow,"Phys.Rev.E,vol.63,pp.035303,2001. [11] C.Tsallis,A.R.Plastino,andW.M.Zheng,\Power-lawsensitivitytoinitialconditions-Newentropicrepresentation,"Chaos,Solitons,&Fractals,vol.8,pp.885-891,1997. [12] U.M.S.Costa,M.L.Lyra,A.R.Plastino,andC.Tsallis,\Power-lawsensitivitytoinitialconditionswithinalogisticlikefamilyofmaps:Fractalityandnonextensivity,"Phys.Rev.E,vol.56,245-250,1997. [13] V.Latora,M.Baranger,A.Rapisarda,andC.Tsallis,\Therateofentropyincreaseattheedgeofchaos,"Phys.Lett.A,vol.273,pp.97-103,2000. [14] E.P.Borges,C.Tsallis,G.F.J.Ananos,andP.M.C.deOliveira,Phys.Rev.Lett.,vol.89,pp.254103,2002. 105

PAGE 106

U.Tirnakli,C.Tsallis,andM.L.Lyra,\Circular-likemaps:sensitivitytotheinitialconditions,multifractalityandnonextensivity,"Eur.Phys.J.B,vol.11,pp.309-315,1999. [16] U.Tirnakli,\Asymmetricunimodalmaps:Someresultsfromq-generalizedbitcumulants,"Phys.Rev.E,vol.62,pp.7857-7860,2000. [17] U.Tirnakli,G.F.J.Ananos,andC.Tsallis,\GeneralizationoftheKolmogorov-Sinaientropy:logistic-likeandgeneralizedcosinemapsatthechaosthreshold,"Phys.Lett.A,vol.289,pp.51-58,2001. [18] U.Tirnakli,\Dissipativemapsatthechaosthreshold:numericalresultsforthesingle-sitemap,"PhysicaA,vol.305,pp.119-123,2001. [19] U.Tirnakli,C.Tsallis,andM.L.Lyra,\Asymmetricunimodalmapsattheedgeofchaos,"Phys.Rev.E,vol.65,pp.036207,2002. [20] U.Tirnakli,\Two-dimensionalmapsattheedgeofchaos:NumericalresultsfortheHenonmap,"Phys.Rev.E,vol.66,pp.066212,2002. [21] J.B.GaoandZ.M.Zheng,\Directdynamicaltestfordeterministicchaosandoptimalembeddingofachaotictimeseries,"Phys.Rev.E,vol.49,pp.3807-3814,1994. [22] J.B.GaoandZ.M.Zheng,\Directdynamicaltestfordeterministicchaos," [23] J.B.Gao,C.C.Chen,S.K.Hwang,andJ.M.Liu,\Noise-inducedchaos,"Int.J.Mod.Phys.B,vol.13,pp.3283-3305,1999. [24] J.B.Gao,S.K.Hwang,andJ.M.Liu,\Whencannoiseinducechaos?"Phys.Rev.Lett.,vol.82,pp.1132-1135,1999. [25] J.B.Gao,S.K.Hwang,andJ.M.Liu,\Eectsofintrinsicspontaneous-emissionnoiseonthenonlineardynamicsofanopticallyinjectedsemiconductorlaser,"Phys.Rev.A,vol.59,pp.1582-1585,1999. [26] S.K.Hwang,J.B.Gao,andJ.M.Liu,\Noise-inducedchaosinanopticallyinjectedsemiconductorlaser,"Phys.Rev.E,vol.61,pp.5162-5170,2000. [27] J.B.Gao,W.W.Tung,andN.Rao,\NoiseInducedHopfBifurcation-typeSequenceandTransitiontoChaosintheLorenzEquations,"Phys.Rev.Lett.,vol.89,pp.254101,2002. [28] J.Hu,J.B.Gao,andK.D.White,\Estimatingmeasurementnoiseinatimeseriesbyexploitingnonstationarity,"Chaos,Solitons,&Fractals,vol.22,pp.807-819,2004. 106

PAGE 107

C.J.Cellucci,A.M.Albano,P.E.Rapp,R.A.Pittenger,R.C.Josiassen,\Detectingnoiseinatimeseries,"Chaos,vol.7,pp.414-422,1997. [30] J.B.GaoandW.W.Tung,\Pathologicaltremorsasdiusionalprocesses,"BiologicalCybernetics,vol.86,pp.263-270,2002. [31] R.BandyopadhyayandA.K.Sood,\Chaoticdynamicsinshear-thickeningsurfactantsolutions,"Europhys.Lett.,vol.56,pp.447-453,2001. [32] R.Bandyopadhyay,G.Basappa,andA.K.Sood,\ObservationofchaoticdynamicsindiluteshearedaqueoussolutionsofCTAT,"Phys.Rev.Lett.,vol.84,pp.2022-2025,2000. [33] S.Venkadesan,M.C.Valsakumar,K.P.N.Murthy,andS.Rajasekar,\Evidenceforchaosinanexperimentaltimeseriesfromserratedplasticow,"Phys.Rev.E,vol.54,pp.611-616,1996. [34] N.H.Packard,J.P.Crutcheld,J.D.Farmer,andR.S.Shaw,\Geometryfromatimeseries,"Phys.Rev.Lett.,vol.45,pp.712-716,1980. [35] F.Takens,inDynamicalsystemsandturbulence,LectureNotesinMathematics,vol.898,pp.366,1981,editedbyD.A.RandandL.S.Young,Springer-Verlag,Berlin. [36] T.Sauer,J.A.Yorke,andM.Casdagli,\Embedology,"J.Stat.Phys.,vol.65,pp.579-616,1991. [37] J.B.GaoandZ.M.Zheng,\Localexponentialdivergenceplotandoptimalembeddingofachaotictimeseries,"Phys.Lett.A,vol.181,pp.153-158,1993. [38] S.Sato,M.Sano,andY.Sawada,\PracticalmethodsofmeasuringthegeneralizeddimensionandthelargestLyapunovexponentinhighdimensionalchaoticsystems,"Prog.Theor.Phys.,vol.77,pp.1-5,1987. [39] M.C.MackeyandL.Glass,\Oscillationandchaosinphysiologicalcontrolsystems,"Science,vol.197,pp.287-288,1977. [40] J.B.Gao,W.W.Tung,Y.H.Cao,J.Hu,andY.Qi,\Power-lawsensitivitytoinitialconditionsinatimeserieswithapplicationstoepilepticseizuredetection,"PhysicaA,vol.353,pp.613-624,2005. [41] J.P.EckmannandD.Ruelle,\Ergodictheoryofchaosandstrangeattractors,"Rev.Mod.Phys.,vol.57,pp.617-656,1985. [42] C.Beck,\DynamicalFoundationsofNonextensiveStatisticalMechanics,"Phys.Rev.Lett.,vol.87,pp.180601,2001. [43] B.B.Mandelbrot,TheFractalGeometryofNature,SanFrancisco,Freeman,1982. 107

PAGE 108

D.Saupe,\Algorithmsforrandomfractals",in:H.Peitgen,D.Saupe(Eds.),TheScienceofFractalImages,Springer-Verlag,Berlin,1988,pp.71-113. [45] W.H.Press,\Flickernoisesinastronomyandelsewhere",Commentson Astrophysics,vol.7.pp.103-119,1978. [46] P.Bak,HowNatureWorks:theScienceofSelf-OrganizedCriticality,Copernicus,1996. [47] G.M.Wornell,Signalprocessingwithfractals:awavelet-basedapproach,PrenticeHall,1996. [48] W.E.Leland,M.S.Taqqu,W.Willinger,andD.V.Wilson,\Ontheself-similarnatureofEthernettrac(extendedversion)",IEEE/ACMTrans.onNetworking,vol.2,1-15,1994. [49] J.Beran,R.Sherman,M.S.Taqqu,andW.Willinger,\Long-range-dependenceinvariable-bit-ratevideotrac",IEEETrans.onCommun.,vol.43,pp.1566-1579,1995. [50] V.PaxsonandS.Floyd,\WideAreaTrac|ThefailureofPoissonmodeling,"IEEE/ACMTrans.onNetworking,vol.3,pp.226-244,1995. [51] W.LiandK.Kaneko,\Long-RangeCorrelationandPartial1/fSpectruminaNon-CodingDNASequence",Europhys.Lett.,vol.17,pp.655-660,1992. [52] R.Voss,\Evolutionoflong-rangefractalcorrelationsand1/fnoiseinDNAbasesequences",Phys.Rev.Lett.,vol.68,pp.3805-3808,1992. [53] C.K.Peng,S.V.Buldyrev,A.L.Goldberger,S.Havlin,F.Sciortino,M.SimonsandH.E.Stanley,\Long-RangeCorrelationsinNucleotideSequences,"Nature,vol.356,pp.168-171,1992. [54] D.L.Gilden,T.Thornton,andM.W.Mallon,\1/fnoiseinhumancognition",Science,vol.267,pp.1837-1839,1995. [55] Y.H.Zhou,J.B.Gao,K.D.White,I.Merk,andK.Yao,\PerceptualDominanceTimeDistributionsinMultistableVisualPerception,"BiologicalCybernetics,vol.90,pp.256-263,2004. [56] J.B.Gao,V.A.Billock,I.Merk,W.W.Tung,K.D.White,J.G.Harris,andV.P.Roychowdhury,\Inertiaandmemoryinambiguousvisualperception,"CognitiveProcessing,vol.7,pp.105-112,2006. [57] Y.Chen,M.Ding,andJ.A.S.Kelso,\LongMemoryProcesses(1/falphaType)inHumanCoordination",Phys.Rev.Lett.,vol.79,pp.4501-4504,1997. 108

PAGE 109

J.J.CollinsandC.J.DeLuca,\RandomWalkingduringQuietStanding,"Phys.Rev.Lett.,vol.73,pp.764-767,1994. [59] V.A.Billock,\Neuralacclimationto1/fspatialfrequencyspectrainnaturalimagestransducedbythehumanvisualsystem,"PhysicaD,vol.137,pp.379-391,2000. [60] V.A.Billock,G.C.deGuzman,andJ.A.S.Kelso,\Fractaltimeand1/fspectraindynamicimagesandhumanvision",PhysicaD,vol.148,pp.136-146,2001. [61] M.Wolf,\1/fnoiseinthedistributionofprimes,"PhysicaA,vol.241,pp.493-499,1997. [62] J.B.Gao,YinheCao,andJae-MinLee,\PrincipalComponentAnalysisof1=fNoise,"Phys.Lett.A,vol.314,pp.392-400,2003. [63] B.B.MandelbrotandV.Ness,\FractionalBrownianmotions,fractionalnoisesandapplications,"SIAMRev.,vol.10,pp.422-437,1968. [64] M.S.Taqqu,W.Willinger,andR.Sherman,Proofofafundamentalresultinself-similartracmodeling,ACMSIGCOMMComputerCommunicationReview,vol.27,pp.5-23,1997. [65] B.B.Mandelbrot,FractalsandScalinginFinance,NewYork:Springer,1997. [66] T.Solomon,E.Weeks,andH.Swinney,\ObservationofAnomalousDiusionandLevyFlightsinaTwoDimensionalRotatingFlow,"Phys.Rev.Lett.,vol.71,pp.3975-3979,1993. [67] T.Geisel,J.Nierwetberg,andA.Zacherl,\AcceleratedDiusioninJosephsonJunctionsandRelatedChaoticSystems,"Phys.Rev.Lett.,vol.54,pp.616-620,1985. [68] S.Stapf,R.Kimmich,andR.Seitter,\ProtonandDeuteronField-cyclingNMRrelaxometryofLiquidsinPorousGlasses:EvidenceofLevy-walkStatistics,"Phys.Rev.Lett.,vol.75,pp.2855-2859,1993. [69] G.M.Viswanathan,V.Afanasyev,S.V.Buldyrev,E.J.Murphy,P.A.Prince,andH.E.Stanley,\Levyightsearchpatternsofwanderingalbatrosses,"Nature,vol.381,pp.413-415,1996. [70] G.M.Viswanathan,S.V.Buldyrev,S.Havlin,M.G.E.daLuz,E.P.Raposo,andH.E.Stanley,\Optimizingthesuccessofrandomsearches,"Nature,vol.401,pp.911-914,1999. [71] J.B.Gao,J.Hu,W.W.Tung,Y.H.Cao,N.Sarshar,andV.P.Roychowdhury,\Assessmentoflongrangecorrelationintimeseries:Howtoavoidpitfalls,"Phys.Rev.E,vol.73,pp.016117,2006. 109

PAGE 110

C.K.Peng,S.V.Buldyrev,S.Havlin,M.Simons,H.E.Stanley,andA.L.Goldberger,\Mosaicorganizationofdnanucleotides,"Phys.Rev.E,vol.49,pp.1685-1689,1994. [73] A.JanickiandA.Weron,\Canonesee-stablevariablesandprocesses,"StatisticalScience,vol.9,pp.109-126,1994. [74] M.Kanter,\Stabledensitiesunderchangeofscaleandtotalvariationinequalities,"AnnalsofProbability,vol.3,pp.697-707,1975. [75] J.M.Chambers,C.L.Mallows,B.W.Stuck,\Methodforsimulatingstablerandom-variables,"JournaloftheAmericanStatisticaLAssociation,vol.71,pp.340-344,1976. [76] S.ResnickandG.Samorodnitsky:Fluidqueues,on/oprocesses,andteletracmodelingwithhighlyvariableandcorrelatedinputs,inSelf-SimilarTracandPerformanceEvaluation,eds,K.ParkandW.Willinger,John-Wiley&Sons(2000),pp.171-192. [77] J.B.GaoandI.Rubin,\MultifractalmodelingofcountingprocessesofLong-Range-DependentnetworkTrac,"ComputerCommunications,vol.24,pp.1400-1410,2001;\MultiplicativeMultifractalModelingofLong-Range-Dependent(LRD)TracinComputerCommunicationsNetworks,"J.NonlinearAnalysis,vol.47,pp.5765-5774,2001;\MultiplicativemultifractalmodelingofLong-Range-Dependentnetworktrac,"Int.J.Comm.Systems,vol.14,pp.783-801,2001. [78] R.AlbertandA.L.Barabasi,\Statisticalmechanicsofcomplexnetworks,"Rev.ModernPhys.,vol.74,pp.47-97,2002. [79] C.Labovitz,G.R.Malan,andF.Jahanian,\Internetroutinginstability,"IEEE/ACMTransactionsonNetworking,vol.6,pp.515-528,1998. [80] J.C.R.Bennett,C.Partridge,andN.Shectman,\Packetreorderingisnotpathologicalnetworkbehavior,"IEEE/ACMTransactionsonNetworking,vol.7,pp.789-798,1999. [81] A.ErramilliandL.J.Forys,\TracsynchronizationEectsinTeletracSystems,"Proc.ITC-13,Copenhagen,1991. [82] R.Srikant,TheMathematicsofInternetCongestionControl,Birkhauser,2004. [83] C.Li,G.Chen,X.Liao,andJ.Yu,\Experimentalqueueinganalysiswithlong-rangedependentpackettrac,"Chaos,Solitons,&Fractals,vol.19,pp.853-868,2004. [84] A.VeresandM.Boda,\ThechaoticnatureofTCPcongestioncontrol,"Proc.ofInfocom,Piscataway,NJ,pp.1715-1723,2000. 110

PAGE 111

N.S.V.RaoandL.O.Chua,\Ondynamicsofnetworktransportprotocols,"Proc.WorkshoponSignalProcessing,Communications,ChaosandSystems,2002. [86] P.Ranjan,E.H.Abed,andR.J.La,\NonlinearinstabilitiesinTCP-RED,"IEEE/ACMTransactionsonNetworking,vol.12,pp.1079-1092,2004. [87] N.S.V.Rao,J.GaoandL.O.Chua,\Ondynamicsofnetworktransportprotocols,"inComplexDynamicsinCommunicationNetworks,editedbyL.KocarevandG.Vattay,2004. [88] J.GaoandN.S.V.Rao,\ComplicatedDynamicsofInternetTransportProtocols,"IEEECommun.Lett.,vol.9,pp.4-6,2005. [89] D.R.Figueiredo,B.Liu,V.Misra,andD.Towsley,\OntheAutocorrelationStructureofTCPTrac,"ComputerNetworks,vol.40,pp.339-361,2002;SpecialIssueon\AdvancesinModelingandEngineeringofLong-RangeDependentTrac",ComputerSci.Tech.Report00-55,Univ.ofMassachusetts,Nov.2002. [90] J.B.GaoandN.S.V.Rao,\SynchronizedOscillations,Quasi-Periodicity,andChaosinTCP,"Dept.ofElectricalandComputerEngineering,UniversityofFlorida,Technicalreport. [91] S.FloydandE.Kohler,inFirstWorkshoponHotTopicsinNetworks,28-29October2002,Princeton,NewJersey,USA. [92] K.ParkandW.Willinger,Self-SimilarNetworkTracandPerformance Evaluation,Wiley,NewYork,2000. [93] V.FiroiuandM.Borden,\Astudyofactivequeuemanagementforcongestioncontrol,"Proc.ofInfocom,vol.3,pp.1435-1444,2000. [94] V.Misra,W.Gong,andD.Towsley,\Astudyofactivequeueforcongestioncontrol,"Proc.ofSigcomm,2000. [95] C.Hollot,V.Misra,D.Towsley,andW.Gong,\AcontroltheoreticanalysisofRED,"Proc.ofInfocom,vol.3,pp.1510-1519,2001. [96] P.Kuusela,P.Lassilaand,andJ.Virtamo,\StabilityofTCPREDcongestioncontrol,"inProc.ofITC-17,pp.655-666,Elsevier. [97] N.S.V.Rao,Q.WuandS.S.Iyengar,\Onthroughputstabilizatinofnetworktransport,"IEEECommun.Lett.,vol.8,pp.66-68,2004. [98] T.Sauer,\ReconstructionofDynamical-SystemsfromInterspikeIntervals,"Phys.Rev.Lett.,vol.72,pp.3811-3814,1994. 111

PAGE 112

J.B.Gao,\Recognizingrandomnessinatimeseries,"PhysicaD,vol.106,pp.49-56,1997. [100] S.HaykinandS.Puthusserypady,\Chaoticdynamicsofseaclutter,"Chaos,vol.7,pp.777-802,1997. [101] S.Haykin,Chaoticdynamicsofseaclutter,JohnWiley,1999. [102] J.L.Noga,\BayesianState-SpaceModellingofSpatio-TemporalNon-GaussianRadarReturns,"Ph.Dthesis,CambridgeUniversity,1998. [103] M.Davies,\LookingforNon-LinearitiesinSeaClutter,"IEERadarandSonarSignalProcessing,Peebles,Scotland,UK,July1998. [104] M.R.CowperandB.Mulgrew,\NonlinearProcessingofHighResolutionRadarSeaClutter,"Proc.IJCNN,vol.4,pp.2633-2638,1999. [105] C.P.Unsworth,M.R.Cowper,S.McLaughlin,andB.Mulgrew,\Falsedetectionofchaoticbehaviorinthestochasticcompoundk-distributionmodelofradarseaclutter,"Proc.10thIEEEWorkshoponSSAP,August2000,pp.296-300. [106] C.P.Unsworth,M.R.Cowper,S.McLaughlin,andB.Mulgrew,\Re-examiningthenatureofseaclutter",IEEProc.RadarSonarNavig.,vol.149,pp.105-114,2002. [107] J.B.GaoandK.Yao,\Multifractalfeaturesofseaclutter,"IEEERadarConference,LongBeach,CA,April,2002. [108] J.B.Gao,S.K.Hwang,H.F.Chen,Z.Kang,K.Yao,andJ.M.Liu,\Canseaclutterandindoorradiopropagationbemodeledasstrangeattractors?"The7thExperimentalChaosConference,SanDiego,USA,Augustpp.25-29,2002. [109] P.GrassbergerandI.Procaccia,\Measuringthestrangenessofthestrangeattractor,"PhysicaD,vol.9D,pp.189-208,1983. [110] S.Haykin,R.Bakker,andB.W.Currie,\Uncoveringnonlineardynamics-thecasestudyofseaclutter,"Proc.IEEE,vol.90,pp.860-881,2002. [111] M.McDonaldandA.Damini,\Limitationsofnonlinearchaoticdynamicsinpredictingseaclutterreturns,"IEEProc.RadarSonarNavig.,vol.151,pp.105-113,2004. [112] F.A.Fay,J.Clarke,andR.S.Peters,\Weibulldistributionappliedtosea-clutter,"Proc.IEEConf.Radar,London,U.K.,pp.101-103,1977. [113] F.E.Nathanson,Radardesignprinciples,McGrawHill,pp.254-256,1969. 112

PAGE 113

E.JakemanandP.N.Pusey,\AmodelfornonRayleighseaecho,"IEEETransAntennas&Propagation,vol.24,pp.806-814,1976. [115] S.SayamaandM.Sekine,\Log-normal,log-WeibullandK-distributedseaclutter,"IEICETrans.Commun.,vol.E85-B,pp.1375-1381,2002. [116] F.Gini,A.Farina,andM.Montanari,\Vectorsubspacedetectionincompound-Gaussianclutter,PartII:performanceanalysis,"IEEETransactionsonAerospaceandElectronicSystems,vol.38,pp.1312-1323,2002. [117] G.DavidsonandH.D.Griths,\Waveletdetectionschemeforsmalltargetsinseaclutter,"ElectronicsLett.,vol.38,pp.1128-1130,2002. [118] T.BhattacharyaandS.Haykin,\NeuralNetwork-basedRadarDetectionforanOceanEnvironment,"IEEETrans.AerospaceandElectronicSystems,vol.33,pp.408-420,1997. [119] N.Xie,H.Leung,andH.Chan,\Amultiple-modelpredictionapproachforseacluttermodeling,"IEEETran.GeosciRemote,vol.41,pp.1491-1502,2003. [120] C.P.Lin,M.Sano,S.Sayama,andM.Sekine,\Detectionofradartargetsembeddedinseaiceandseaclutterusingfractals,wavelets,andneuralnetworks,"IEICETrans.Commun.,vol.E83B,pp.1916-1929,2000. [121] C.P.Lin,M.Sano,S.Obi,S.Sayama,andM.Sekine,\DetectionofoilleakageinSARimagesusingwaveletfeatureextractorsandunsupervisedneuralclassiers,"IEICETrans.Commun.,vol.E83B,pp.1955-1962,2000. [122] T.Lo,H.Leung,J.LitvaandS.Haykin,\Fractalcharacterisationofsea-scatteredsignalsanddetectionofsea-surfacetargets,"IEEProc.,vol.F140,pp.243-250,1993. [123] C.-P.Lin,M.Sano,andM.Sekine,\Detectionofradartargetsbymeansoffractalerror,"IEICETrans.Commun.,vol.E80-B,pp.1741-1748,1997. [124] B.E.Cooper,D.L.Chenoweth,andJ.E.Selvage,\Fractalerrorfordetectingman-madefeaturesinaerialimages,"ElectronicsLetters,vol.30,pp.554-555,1994. [125] J.Hu,W.W.Tung,andJ.B.Gao,\Detectionoflowobservabletargetswithinseaclutterbystructurefunctionbasedmultifractalanalysis,"IEEETransactionsonAntennas&Propagation,vol.54,pp.135-143,2006. [126] E.Aurell,G.Boetta,A.Crisanti,G.Paladin,andA.Vulpiani,\Predictabilityinthelarge:AnextensionoftheconceptofLyapunovexponent,"J.PhysicsA,vol.30,pp.1-26,1997. [127] A.Wolf,J.B.Swift,H.L.Swinney,andJ.A.Vastano,\DeterminingLyapunovexponentsfromatimeseries,"PhysicaD,vol.16,pp.285-317,1985. 113

PAGE 114

M.Cencini,M.Falcioni,E.Olbrich,H.Kantz,andA.Vulpiani,\Chaosornoise:Dicultiesofadistinction,"Phys.Rev.E,vol.62,pp.427-437,2000. [129] G.J.Ortega,\Anewmethodtodetecthiddenfrequenciesinchaotictimeseries,"Phys.Lett.A,vol.209,pp.351-355,1995. [130] G.J.Ortega,\InvariantmeasuresasLagrangianvariables:Theirapplicationtotimeseriesanalysis,"Phys.Rev.Lett.,vol.77,pp.259-262,1996. [131] J.L.Chern,J.Y.Ko,J.S.Lih,H.T.Su,andR.R.Hsu,\Recognizinghiddenfrequenciesinachaotictimeseries,"Phys.Lett.A,vol.238,pp.134-140,1998. 114

PAGE 115

JingHureceivedtheB.S.andM.E.degreesfromtheDepartmentofElectronicsandInformationEngineeringfromHuazhongUniversityofScienceandTechnology(HUST),Wuhan,China,in2000and2002,respectively.InMay2007,ShewasawardedthePh.D.degreefromtheDepartmentofElectricalandComputerEngineeringattheUniversityofFlorida,Gainesville,Florida.Herresearchinterestsincludesignalandimageprocessing,informationscience,nonlineartimeseriesanalysis,biologicaldataanalysis,bioinformatics,andbiomedicalengineering. 115