UFDC Home  myUFDC Home  Help 
PAGE 105
[1] P.Fairley,\Theunrulypowergrid,"IEEESpectrum,vol.41,pp.2227,2004. [2] C.Tsallis,\PossiblegeneralizationofBoltzmannGibbsstatistics,"JStatPhys,vol.52,pp.479487,1988. [3] A.R.PlastinoandA.Plastino,\Stellarpolytropesandtsallisentropy,"Phys.Lett.A,vol.174,pp.384386,1993. [4] A.Lavagno,G.Kaniadakis,M.RegoMonteiro,P.Quarati,andC.Tsallis,\Nonextensivethermostatisticalapproachofthepeculiarvelocityfunctionofgalaxyclusters,"Astrophys.Lett.Commun.,vol.35,pp.449455,1998. [5] M.AntoniandS.Ruo,\ClusteringandrelaxationinHamiltonianlongrangedynamics,"Phys.Rev.E,vol.52,pp.23612374,1995. [6] V.Latora,A.Rapisarda,andC.Tsallis,\NonGaussianequilibriuminalongrangeHamiltoniansystem,"Phys.Rev.E,vol.64,pp.056134,2001. [7] M.L.LyraandC.Tsallis,\Nonextensivityandmultifractalityinlowdimensionaldissipativesystems,"Phys.Rev.Lett.,vol.80,pp.5356,1998. [8] T.ArimitsuandN.Arimitsu,\Tsallisstatisticsandfullydevelopedturbulence,"J.Phys.A,vol.33,pp.L235L241,2000. [9] C.Beck,\Applicationofgeneralizedthermostatisticstofullydevelopedturbulence,"PhysicaA,vol.277,pp.115123,2000. [10] C.Beck,G.S.Lewis,andH.L.Swinney,\MeasuringnonextensitivityparametersinaturbulentCouetteTaylorow,"Phys.Rev.E,vol.63,pp.035303,2001. [11] C.Tsallis,A.R.Plastino,andW.M.Zheng,\PowerlawsensitivitytoinitialconditionsNewentropicrepresentation,"Chaos,Solitons,&Fractals,vol.8,pp.885891,1997. [12] U.M.S.Costa,M.L.Lyra,A.R.Plastino,andC.Tsallis,\Powerlawsensitivitytoinitialconditionswithinalogisticlikefamilyofmaps:Fractalityandnonextensivity,"Phys.Rev.E,vol.56,245250,1997. [13] V.Latora,M.Baranger,A.Rapisarda,andC.Tsallis,\Therateofentropyincreaseattheedgeofchaos,"Phys.Lett.A,vol.273,pp.97103,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,\Circularlikemaps:sensitivitytotheinitialconditions,multifractalityandnonextensivity,"Eur.Phys.J.B,vol.11,pp.309315,1999. [16] U.Tirnakli,\Asymmetricunimodalmaps:Someresultsfromqgeneralizedbitcumulants,"Phys.Rev.E,vol.62,pp.78577860,2000. [17] U.Tirnakli,G.F.J.Ananos,andC.Tsallis,\GeneralizationoftheKolmogorovSinaientropy:logisticlikeandgeneralizedcosinemapsatthechaosthreshold,"Phys.Lett.A,vol.289,pp.5158,2001. [18] U.Tirnakli,\Dissipativemapsatthechaosthreshold:numericalresultsforthesinglesitemap,"PhysicaA,vol.305,pp.119123,2001. [19] U.Tirnakli,C.Tsallis,andM.L.Lyra,\Asymmetricunimodalmapsattheedgeofchaos,"Phys.Rev.E,vol.65,pp.036207,2002. [20] U.Tirnakli,\Twodimensionalmapsattheedgeofchaos:NumericalresultsfortheHenonmap,"Phys.Rev.E,vol.66,pp.066212,2002. [21] J.B.GaoandZ.M.Zheng,\Directdynamicaltestfordeterministicchaosandoptimalembeddingofachaotictimeseries,"Phys.Rev.E,vol.49,pp.38073814,1994. [22] J.B.GaoandZ.M.Zheng,\Directdynamicaltestfordeterministicchaos," [23] J.B.Gao,C.C.Chen,S.K.Hwang,andJ.M.Liu,\Noiseinducedchaos,"Int.J.Mod.Phys.B,vol.13,pp.32833305,1999. [24] J.B.Gao,S.K.Hwang,andJ.M.Liu,\Whencannoiseinducechaos?"Phys.Rev.Lett.,vol.82,pp.11321135,1999. [25] J.B.Gao,S.K.Hwang,andJ.M.Liu,\Eectsofintrinsicspontaneousemissionnoiseonthenonlineardynamicsofanopticallyinjectedsemiconductorlaser,"Phys.Rev.A,vol.59,pp.15821585,1999. [26] S.K.Hwang,J.B.Gao,andJ.M.Liu,\Noiseinducedchaosinanopticallyinjectedsemiconductorlaser,"Phys.Rev.E,vol.61,pp.51625170,2000. [27] J.B.Gao,W.W.Tung,andN.Rao,\NoiseInducedHopfBifurcationtypeSequenceandTransitiontoChaosintheLorenzEquations,"Phys.Rev.Lett.,vol.89,pp.254101,2002. [28] J.Hu,J.B.Gao,andK.D.White,\Estimatingmeasurementnoiseinatimeseriesbyexploitingnonstationarity,"Chaos,Solitons,&Fractals,vol.22,pp.807819,2004. 106
PAGE 107
C.J.Cellucci,A.M.Albano,P.E.Rapp,R.A.Pittenger,R.C.Josiassen,\Detectingnoiseinatimeseries,"Chaos,vol.7,pp.414422,1997. [30] J.B.GaoandW.W.Tung,\Pathologicaltremorsasdiusionalprocesses,"BiologicalCybernetics,vol.86,pp.263270,2002. [31] R.BandyopadhyayandA.K.Sood,\Chaoticdynamicsinshearthickeningsurfactantsolutions,"Europhys.Lett.,vol.56,pp.447453,2001. [32] R.Bandyopadhyay,G.Basappa,andA.K.Sood,\ObservationofchaoticdynamicsindiluteshearedaqueoussolutionsofCTAT,"Phys.Rev.Lett.,vol.84,pp.20222025,2000. [33] S.Venkadesan,M.C.Valsakumar,K.P.N.Murthy,andS.Rajasekar,\Evidenceforchaosinanexperimentaltimeseriesfromserratedplasticow,"Phys.Rev.E,vol.54,pp.611616,1996. [34] N.H.Packard,J.P.Crutcheld,J.D.Farmer,andR.S.Shaw,\Geometryfromatimeseries,"Phys.Rev.Lett.,vol.45,pp.712716,1980. [35] F.Takens,inDynamicalsystemsandturbulence,LectureNotesinMathematics,vol.898,pp.366,1981,editedbyD.A.RandandL.S.Young,SpringerVerlag,Berlin. [36] T.Sauer,J.A.Yorke,andM.Casdagli,\Embedology,"J.Stat.Phys.,vol.65,pp.579616,1991. [37] J.B.GaoandZ.M.Zheng,\Localexponentialdivergenceplotandoptimalembeddingofachaotictimeseries,"Phys.Lett.A,vol.181,pp.153158,1993. [38] S.Sato,M.Sano,andY.Sawada,\PracticalmethodsofmeasuringthegeneralizeddimensionandthelargestLyapunovexponentinhighdimensionalchaoticsystems,"Prog.Theor.Phys.,vol.77,pp.15,1987. [39] M.C.MackeyandL.Glass,\Oscillationandchaosinphysiologicalcontrolsystems,"Science,vol.197,pp.287288,1977. [40] J.B.Gao,W.W.Tung,Y.H.Cao,J.Hu,andY.Qi,\Powerlawsensitivitytoinitialconditionsinatimeserieswithapplicationstoepilepticseizuredetection,"PhysicaA,vol.353,pp.613624,2005. [41] J.P.EckmannandD.Ruelle,\Ergodictheoryofchaosandstrangeattractors,"Rev.Mod.Phys.,vol.57,pp.617656,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,SpringerVerlag,Berlin,1988,pp.71113. [45] W.H.Press,\Flickernoisesinastronomyandelsewhere",Commentson Astrophysics,vol.7.pp.103119,1978. [46] P.Bak,HowNatureWorks:theScienceofSelfOrganizedCriticality,Copernicus,1996. [47] G.M.Wornell,Signalprocessingwithfractals:awaveletbasedapproach,PrenticeHall,1996. [48] W.E.Leland,M.S.Taqqu,W.Willinger,andD.V.Wilson,\OntheselfsimilarnatureofEthernettrac(extendedversion)",IEEE/ACMTrans.onNetworking,vol.2,115,1994. [49] J.Beran,R.Sherman,M.S.Taqqu,andW.Willinger,\Longrangedependenceinvariablebitratevideotrac",IEEETrans.onCommun.,vol.43,pp.15661579,1995. [50] V.PaxsonandS.Floyd,\WideAreaTracThefailureofPoissonmodeling,"IEEE/ACMTrans.onNetworking,vol.3,pp.226244,1995. [51] W.LiandK.Kaneko,\LongRangeCorrelationandPartial1/fSpectruminaNonCodingDNASequence",Europhys.Lett.,vol.17,pp.655660,1992. [52] R.Voss,\Evolutionoflongrangefractalcorrelationsand1/fnoiseinDNAbasesequences",Phys.Rev.Lett.,vol.68,pp.38053808,1992. [53] C.K.Peng,S.V.Buldyrev,A.L.Goldberger,S.Havlin,F.Sciortino,M.SimonsandH.E.Stanley,\LongRangeCorrelationsinNucleotideSequences,"Nature,vol.356,pp.168171,1992. [54] D.L.Gilden,T.Thornton,andM.W.Mallon,\1/fnoiseinhumancognition",Science,vol.267,pp.18371839,1995. [55] Y.H.Zhou,J.B.Gao,K.D.White,I.Merk,andK.Yao,\PerceptualDominanceTimeDistributionsinMultistableVisualPerception,"BiologicalCybernetics,vol.90,pp.256263,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.105112,2006. [57] Y.Chen,M.Ding,andJ.A.S.Kelso,\LongMemoryProcesses(1/falphaType)inHumanCoordination",Phys.Rev.Lett.,vol.79,pp.45014504,1997. 108
PAGE 109
J.J.CollinsandC.J.DeLuca,\RandomWalkingduringQuietStanding,"Phys.Rev.Lett.,vol.73,pp.764767,1994. [59] V.A.Billock,\Neuralacclimationto1/fspatialfrequencyspectrainnaturalimagestransducedbythehumanvisualsystem,"PhysicaD,vol.137,pp.379391,2000. [60] V.A.Billock,G.C.deGuzman,andJ.A.S.Kelso,\Fractaltimeand1/fspectraindynamicimagesandhumanvision",PhysicaD,vol.148,pp.136146,2001. [61] M.Wolf,\1/fnoiseinthedistributionofprimes,"PhysicaA,vol.241,pp.493499,1997. [62] J.B.Gao,YinheCao,andJaeMinLee,\PrincipalComponentAnalysisof1=fNoise,"Phys.Lett.A,vol.314,pp.392400,2003. [63] B.B.MandelbrotandV.Ness,\FractionalBrownianmotions,fractionalnoisesandapplications,"SIAMRev.,vol.10,pp.422437,1968. [64] M.S.Taqqu,W.Willinger,andR.Sherman,Proofofafundamentalresultinselfsimilartracmodeling,ACMSIGCOMMComputerCommunicationReview,vol.27,pp.523,1997. [65] B.B.Mandelbrot,FractalsandScalinginFinance,NewYork:Springer,1997. [66] T.Solomon,E.Weeks,andH.Swinney,\ObservationofAnomalousDiusionandLevyFlightsinaTwoDimensionalRotatingFlow,"Phys.Rev.Lett.,vol.71,pp.39753979,1993. [67] T.Geisel,J.Nierwetberg,andA.Zacherl,\AcceleratedDiusioninJosephsonJunctionsandRelatedChaoticSystems,"Phys.Rev.Lett.,vol.54,pp.616620,1985. [68] S.Stapf,R.Kimmich,andR.Seitter,\ProtonandDeuteronFieldcyclingNMRrelaxometryofLiquidsinPorousGlasses:EvidenceofLevywalkStatistics,"Phys.Rev.Lett.,vol.75,pp.28552859,1993. [69] G.M.Viswanathan,V.Afanasyev,S.V.Buldyrev,E.J.Murphy,P.A.Prince,andH.E.Stanley,\Levyightsearchpatternsofwanderingalbatrosses,"Nature,vol.381,pp.413415,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.911914,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.16851689,1994. [73] A.JanickiandA.Weron,\Canoneseestablevariablesandprocesses,"StatisticalScience,vol.9,pp.109126,1994. [74] M.Kanter,\Stabledensitiesunderchangeofscaleandtotalvariationinequalities,"AnnalsofProbability,vol.3,pp.697707,1975. [75] J.M.Chambers,C.L.Mallows,B.W.Stuck,\Methodforsimulatingstablerandomvariables,"JournaloftheAmericanStatisticaLAssociation,vol.71,pp.340344,1976. [76] S.ResnickandG.Samorodnitsky:Fluidqueues,on/oprocesses,andteletracmodelingwithhighlyvariableandcorrelatedinputs,inSelfSimilarTracandPerformanceEvaluation,eds,K.ParkandW.Willinger,JohnWiley&Sons(2000),pp.171192. [77] J.B.GaoandI.Rubin,\MultifractalmodelingofcountingprocessesofLongRangeDependentnetworkTrac,"ComputerCommunications,vol.24,pp.14001410,2001;\MultiplicativeMultifractalModelingofLongRangeDependent(LRD)TracinComputerCommunicationsNetworks,"J.NonlinearAnalysis,vol.47,pp.57655774,2001;\MultiplicativemultifractalmodelingofLongRangeDependentnetworktrac,"Int.J.Comm.Systems,vol.14,pp.783801,2001. [78] R.AlbertandA.L.Barabasi,\Statisticalmechanicsofcomplexnetworks,"Rev.ModernPhys.,vol.74,pp.4797,2002. [79] C.Labovitz,G.R.Malan,andF.Jahanian,\Internetroutinginstability,"IEEE/ACMTransactionsonNetworking,vol.6,pp.515528,1998. [80] J.C.R.Bennett,C.Partridge,andN.Shectman,\Packetreorderingisnotpathologicalnetworkbehavior,"IEEE/ACMTransactionsonNetworking,vol.7,pp.789798,1999. [81] A.ErramilliandL.J.Forys,\TracsynchronizationEectsinTeletracSystems,"Proc.ITC13,Copenhagen,1991. [82] R.Srikant,TheMathematicsofInternetCongestionControl,Birkhauser,2004. [83] C.Li,G.Chen,X.Liao,andJ.Yu,\Experimentalqueueinganalysiswithlongrangedependentpackettrac,"Chaos,Solitons,&Fractals,vol.19,pp.853868,2004. [84] A.VeresandM.Boda,\ThechaoticnatureofTCPcongestioncontrol,"Proc.ofInfocom,Piscataway,NJ,pp.17151723,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,\NonlinearinstabilitiesinTCPRED,"IEEE/ACMTransactionsonNetworking,vol.12,pp.10791092,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.46,2005. [89] D.R.Figueiredo,B.Liu,V.Misra,andD.Towsley,\OntheAutocorrelationStructureofTCPTrac,"ComputerNetworks,vol.40,pp.339361,2002;SpecialIssueon\AdvancesinModelingandEngineeringofLongRangeDependentTrac",ComputerSci.Tech.Report0055,Univ.ofMassachusetts,Nov.2002. [90] J.B.GaoandN.S.V.Rao,\SynchronizedOscillations,QuasiPeriodicity,andChaosinTCP,"Dept.ofElectricalandComputerEngineering,UniversityofFlorida,Technicalreport. [91] S.FloydandE.Kohler,inFirstWorkshoponHotTopicsinNetworks,2829October2002,Princeton,NewJersey,USA. [92] K.ParkandW.Willinger,SelfSimilarNetworkTracandPerformance Evaluation,Wiley,NewYork,2000. [93] V.FiroiuandM.Borden,\Astudyofactivequeuemanagementforcongestioncontrol,"Proc.ofInfocom,vol.3,pp.14351444,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.15101519,2001. [96] P.Kuusela,P.Lassilaand,andJ.Virtamo,\StabilityofTCPREDcongestioncontrol,"inProc.ofITC17,pp.655666,Elsevier. [97] N.S.V.Rao,Q.WuandS.S.Iyengar,\Onthroughputstabilizatinofnetworktransport,"IEEECommun.Lett.,vol.8,pp.6668,2004. [98] T.Sauer,\ReconstructionofDynamicalSystemsfromInterspikeIntervals,"Phys.Rev.Lett.,vol.72,pp.38113814,1994. 111
PAGE 112
J.B.Gao,\Recognizingrandomnessinatimeseries,"PhysicaD,vol.106,pp.4956,1997. [100] S.HaykinandS.Puthusserypady,\Chaoticdynamicsofseaclutter,"Chaos,vol.7,pp.777802,1997. [101] S.Haykin,Chaoticdynamicsofseaclutter,JohnWiley,1999. [102] J.L.Noga,\BayesianStateSpaceModellingofSpatioTemporalNonGaussianRadarReturns,"Ph.Dthesis,CambridgeUniversity,1998. [103] M.Davies,\LookingforNonLinearitiesinSeaClutter,"IEERadarandSonarSignalProcessing,Peebles,Scotland,UK,July1998. [104] M.R.CowperandB.Mulgrew,\NonlinearProcessingofHighResolutionRadarSeaClutter,"Proc.IJCNN,vol.4,pp.26332638,1999. [105] C.P.Unsworth,M.R.Cowper,S.McLaughlin,andB.Mulgrew,\Falsedetectionofchaoticbehaviorinthestochasticcompoundkdistributionmodelofradarseaclutter,"Proc.10thIEEEWorkshoponSSAP,August2000,pp.296300. [106] C.P.Unsworth,M.R.Cowper,S.McLaughlin,andB.Mulgrew,\Reexaminingthenatureofseaclutter",IEEProc.RadarSonarNavig.,vol.149,pp.105114,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.2529,2002. [109] P.GrassbergerandI.Procaccia,\Measuringthestrangenessofthestrangeattractor,"PhysicaD,vol.9D,pp.189208,1983. [110] S.Haykin,R.Bakker,andB.W.Currie,\Uncoveringnonlineardynamicsthecasestudyofseaclutter,"Proc.IEEE,vol.90,pp.860881,2002. [111] M.McDonaldandA.Damini,\Limitationsofnonlinearchaoticdynamicsinpredictingseaclutterreturns,"IEEProc.RadarSonarNavig.,vol.151,pp.105113,2004. [112] F.A.Fay,J.Clarke,andR.S.Peters,\Weibulldistributionappliedtoseaclutter,"Proc.IEEConf.Radar,London,U.K.,pp.101103,1977. [113] F.E.Nathanson,Radardesignprinciples,McGrawHill,pp.254256,1969. 112
PAGE 113
E.JakemanandP.N.Pusey,\AmodelfornonRayleighseaecho,"IEEETransAntennas&Propagation,vol.24,pp.806814,1976. [115] S.SayamaandM.Sekine,\Lognormal,logWeibullandKdistributedseaclutter,"IEICETrans.Commun.,vol.E85B,pp.13751381,2002. [116] F.Gini,A.Farina,andM.Montanari,\VectorsubspacedetectionincompoundGaussianclutter,PartII:performanceanalysis,"IEEETransactionsonAerospaceandElectronicSystems,vol.38,pp.13121323,2002. [117] G.DavidsonandH.D.Griths,\Waveletdetectionschemeforsmalltargetsinseaclutter,"ElectronicsLett.,vol.38,pp.11281130,2002. [118] T.BhattacharyaandS.Haykin,\NeuralNetworkbasedRadarDetectionforanOceanEnvironment,"IEEETrans.AerospaceandElectronicSystems,vol.33,pp.408420,1997. [119] N.Xie,H.Leung,andH.Chan,\Amultiplemodelpredictionapproachforseacluttermodeling,"IEEETran.GeosciRemote,vol.41,pp.14911502,2003. [120] C.P.Lin,M.Sano,S.Sayama,andM.Sekine,\Detectionofradartargetsembeddedinseaiceandseaclutterusingfractals,wavelets,andneuralnetworks,"IEICETrans.Commun.,vol.E83B,pp.19161929,2000. [121] C.P.Lin,M.Sano,S.Obi,S.Sayama,andM.Sekine,\DetectionofoilleakageinSARimagesusingwaveletfeatureextractorsandunsupervisedneuralclassiers,"IEICETrans.Commun.,vol.E83B,pp.19551962,2000. [122] T.Lo,H.Leung,J.LitvaandS.Haykin,\Fractalcharacterisationofseascatteredsignalsanddetectionofseasurfacetargets,"IEEProc.,vol.F140,pp.243250,1993. [123] C.P.Lin,M.Sano,andM.Sekine,\Detectionofradartargetsbymeansoffractalerror,"IEICETrans.Commun.,vol.E80B,pp.17411748,1997. [124] B.E.Cooper,D.L.Chenoweth,andJ.E.Selvage,\Fractalerrorfordetectingmanmadefeaturesinaerialimages,"ElectronicsLetters,vol.30,pp.554555,1994. [125] J.Hu,W.W.Tung,andJ.B.Gao,\Detectionoflowobservabletargetswithinseaclutterbystructurefunctionbasedmultifractalanalysis,"IEEETransactionsonAntennas&Propagation,vol.54,pp.135143,2006. [126] E.Aurell,G.Boetta,A.Crisanti,G.Paladin,andA.Vulpiani,\Predictabilityinthelarge:AnextensionoftheconceptofLyapunovexponent,"J.PhysicsA,vol.30,pp.126,1997. [127] A.Wolf,J.B.Swift,H.L.Swinney,andJ.A.Vastano,\DeterminingLyapunovexponentsfromatimeseries,"PhysicaD,vol.16,pp.285317,1985. 113
PAGE 114
M.Cencini,M.Falcioni,E.Olbrich,H.Kantz,andA.Vulpiani,\Chaosornoise:Dicultiesofadistinction,"Phys.Rev.E,vol.62,pp.427437,2000. [129] G.J.Ortega,\Anewmethodtodetecthiddenfrequenciesinchaotictimeseries,"Phys.Lett.A,vol.209,pp.351355,1995. [130] G.J.Ortega,\InvariantmeasuresasLagrangianvariables:Theirapplicationtotimeseriesanalysis,"Phys.Rev.Lett.,vol.77,pp.259262,1996. [131] J.L.Chern,J.Y.Ko,J.S.Lih,H.T.Su,andR.R.Hsu,\Recognizinghiddenfrequenciesinachaotictimeseries,"Phys.Lett.A,vol.238,pp.134140,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



Full Text  
NEW APPROACHES TO MULTISCALE SIGNAL PROCESSING By JING HU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 S2007 Jing Hu To my family, teachers and friends ACKENOWLED GMENTS 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 SuShing C'I. in~ Yuguang Fang, .Jose C. Principe, and K~eith D. White), for their interest in my work and valuable aHulMy 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 pillarlike support and their unwavering relief and confidence in me without this I do not think I would have made it this far! TABLE OF CONTENTS page ACK(NOWLEDGMENTS ......... . . 4 LIST OF FIGURES ......... .. . 7 ABSTRACT ......... ..... . 9 CHAPTER 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 Fractal Theories .. .. 2:3 1.6 Importance of the Concept of Scale . ..... .. 24 1.7 Structure of this Dissertation ........ .. .. 24 2 GENERALIZATION OF CHAOS: POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) ......... ... 27 2.1 Dynamical Test for C'!s I... .... .. . .. .. 27 2.2 General Computational Framework for Powerlaw Sensitivity to Initial Conditions (PSIC) .. . .. .. .. .. 3:3 2.3 C'll I l .terizing Edge of C'I!s I. by Powerlaw Sensitivity of Initial Conditions (PSIC) .... ........ ............ :36 :3 CHARACTERIZING RANDOM FRACTALS BY POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) . ..... .. 41 :3.1 C'll I l .terizing 1/ fP Processes by Powerlaw Sensitivity to Initial Conditions (PSIC) .... ........ .. .. ...... ..... 41 :3.1.1 C'll I l .terizing Fractional Brownian Motion (fBm) Processes by Powerlaw 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 Powerlaw Sensitivity to Initial Conditions (PSIC) 50 :3.2 C'll I l .terizing Levy Processes by Powerlaw Sensitivity to Initial Conditions (PSIC) .... ........ ............ 54 4 STITDY OF INTERNET DYNAMICS BY POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) ........ ... .. 60 4.1 Study of Transport Dynamics by Network Simulation .. .. .. .. 61 4.2 Complicated Dynamics of Internet Transport Protocols .. .. .. .. 69 5 STUDY OF SEA CLUTTER RADAR RETURNS BY POWERLAH SENSITIVITY TO INITIAL CONDITIONS (PSIC) .. .. .. .. 76 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 Powerlaw Sensitivity to Initial Conditions (PSIC) ........ . .. 86 6 MITLTISCALE ANALYSIS BY SCALEDEPENDENT LYAPITNOV EXPONENT (SDLE) ........... ......... .. 89 6.1 Basic Theory 6.2 Classification of Complex Motions. 6.2.1 C'! I ..~ Noisy (I I ..~ and Noiseinduced ChI ... . 6.2.2 Processes Defined by Powerlaw Sensitivity to Initial (PSIC). 6.2.3 Complex Motions with Multiple Scaling Behaviors. 6.3 C'!I. .) ..terizing Hidden Frequencies Conditions 7 CONCLUSIONS ......... ... .. 104 REFERENCES .. .......... ........... 105 BIOGRAPHICAL SK(ETCH ...... .. 115 LIST OF FIGURES Figure page 11 Example of a trajectory in the phase space ..... .. .. 14 12 Giant ocean wave (tsunami). ......... .. 17 13 Example of sea clutter data. ......... .. 18 14 Example of Heart rate variability data for a normal subject. .. .. .. 19 21 Timedependent exponent A(k) vs. evolution time k curves for Lorenz system. .32 22 Timedependent exponent A(k) vs. evolution time k curves for the MackeyGlass system. ......... ..... 33 23 Timedependent exponent A(t) vs. t curves for time series generated from logistic map. ......... ............ .... 37 24 Timedependent exponent A(t) vs. t curves for time series generated from Henon map. ............ ............... 40 31 White noise and Brownian motion . ...... .. 46 32 Several fBm processes with different H. . ... .. 48 33 Timedependent exponent A(k) vs. In k curves for fractional Brownian motion. 51 34 Example of ON/OFF processes. .... ... .. 51 35 Timedependent exponent A(k) vs. In k curves for ON/OFF processes with Pareto distributed ON and OFF periods. ... ... .. 53 36 Examples for Brownian motions and Levy motions. .. .. .. 57 37 Timedependent exponent A(k) vs. In k curves for Levy processes. .. .. .. 58 41 Example of quasiperiodic congestion window size W(i) time series. .. .. .. 63 42 Examples of chaotic T(i) time series. . ..... .. .. 65 43 Complementary cumulative distribution functions (CCDFs) for the two chaotic T(i) time series of Figs. 42(a,c). . ... ... .. 66 44 The time series T(i) extracted from a W(i) time series collected using net100 instruments over ORNLLSU connection. ..... .. .. 68 45 Time series for the congestion window size cwnd for ORNLGaTech connection. 72 46 Timedependent exponent A(k) vs. k curves for cwnd data corresponding to Figf.45. ......... ..... 73 47 Timedependent exponent A(k) vs. In k curves for cwnd data corresponding to Figf.45. ......... ..... 75 51 Collection of sea clutter data ......... ... 79 52 Typical sea clutter amplitude data. ........ .. .. 79 53 Two short segments of the amplitude sea clutter data severely affected by clipping. 80 54 Examples of the timedependent exponent A(k) vs. k curves for the sea clutter data ........ ... . .... 81 55 Variations of the Lyapunov exponent A estimated by conventional methods vs. the 14 range bins for the 10 HH measurements. .... .. .. 84 56 Variations of the scaledependent exponent corresponding to large scale vs. the range bins for the 10 HH measurements. ..... .. .. 85 57 Examples of the A(k) vs. In k curves for the sea clutter data. .. .. .. .. 87 58 Variation of the P parameter with the 14 range bins. .. .. ... .. 87 59 Frequencies of the bins without targets and the bins with primary targets for the HH datasets. .. ... .. .. 88 61 Scaledependent Lyapunov exponent A(e) curves for the Lorenz system and logistic map ........ ... . .... 93 62 Scaledependent Lyapunov exponent A(e) curve for the MackeyGlass system. .. 94 63 The functions F(x) and G(x). ......... ... .. 97 64 Scaledependent Lyapunov exponent A(e) for the model described by Eq. 69. .97 65 Powerspectral density (PSD) for time series generated from Lorenz system. 100 66 Lorenz attractor. .. ... .. .. 101 67 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 NEW APPROACHES TO MITLTISCALE SIGNAL PROCESSING By 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, powerlaw 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, scaledependent 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. CHAPTER 1 INTRODUCTION 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 l1~ we decisions made by market makers, as well as influenced by interactions of heterogeneous traders, including intradwl~i traders, shortperiod traders, and longperiod 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, nanosciences, 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: * Powerlaw 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. * Scaledependent 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 lowdintensional chaos, noisy chaos, noiseinduced chaos, random 1/ fP and ccstable 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 scales. 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 (11) 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 equation 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 (12) 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. 11) can be rewritten in the form of (Eq. 12), thanks to the following trick: we introduce new variables xl = x and x2 x Then xl = Z2. Hence the equivalent system (Eq. 12) is xi = x2 b k xa2=x 2 1i m m This system is said to be linear, because all the xi on the righthand 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 11. 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. 11. 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 system! The phase space for the general system (Eq. 12) is the space with coordinates xl, x,. Because this space is ndimensional, we will refer to Eq. 12 as an ndimensional system or an nthorder 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!~!sI 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: Largescale 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 longrange dependent properties in traffic, * Space, essentially due to topology and geography and again manifested by scalefree properties, * 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. 12. To be quantitative, in Fig. 13, 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 12. 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 13(e,\ where Xt(m) is the nonoverlapping 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 firstorder 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)) ~ ml, 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 powerlaw~. However, neither behavior is observed in Fig. 13(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) 8 log2 m U 2U 4U bU0 tU 1UU 12U Time (sec) Figure 13. 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) (Xtmm+l + + Xtm)/m, t = 1, 2, } is the nonoverlapping 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")) ~ m1 Here Var(Xe(m)) decatys 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 (b) (d) a,6 4 E 0 1 2 3 4 5 6 7 Index n4 x 104 100 110 (b) (c) 100 95 90 I80 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 102 10 102 10 f f Figure 14. 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. 14, for a normal young subject. Evidently, the signal is highly nonstationary and niultiscaled, appearing oscillatory for some period of time (Figs. 14(h,d)), and then varying as a powerlaw for another period of time (Figs. 14(c,e)). The latter is an example of the socalled 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 data. 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 twentiethcentury, 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 longterm behavior in a deterministic system that exhibits sensitive dependence on initial conditions. *"Aperiodic longterm b. !! lit.! I means that there are trajectories which do not settle down to fixed points, periodic orbits, or quasiperiodic 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 selfsimilar features on many or all scales. Selfsimilarity 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 in 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 powerlaw relation, which translates to a linear relation in loglog 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 stepsize, and different hikers may have different stepsizes a person riding a horse has a huge stepsize, while a group of people with a little child must have a tiny stepsize. The length of our route is L = N(e) (13) 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 powerlaw manner, N~e) ED e 0(14) with D being a noninteger, 1 < D < 2. Such a nonintegral D is often called the fractal dimension to emphasize the fragmented and irregular characteristics of the object under study. Let us now try to understand the meaning of the nonintegral 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) ~ e1 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) ~ E2 to cover the area, or NV(e) ~ e3 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 lost! 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, nanosciences, 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 scaledependent measure. The most ideal scenario is that a scaledependent 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 powerlaw 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 highdimensional case. We shall present a general computational framework for assessing PSIC from real world data. We also apply the proposed computational procedure to study noisefree 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 longrangecorrelations 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 scaledependent 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 lowdintensional chaos, noisy chaos, noiseinduced 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 CHAPTER 2 GENERALIZATION OF CHAOS: POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) The fornialism of nonextensive statistical niechanics (NESM) [2] has found numerous applications to the study of systems with longrangeinteractions [36] and niultifractal behavior [7, 8], and fully developed turbulence [810], 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 powerlaw 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 1D logisticlike maps and 2D Henon niap at the edge of chaos [7, 1120], 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 lowdintensional 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 highdintensional 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: noisefree 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 I.ng.", 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 lowdintensional deterministic chaos, we have d(t) ~ d(0)ex' (21) 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. 14. A nonintegral fractal dimension of an attractor contrasts sharply with the integervalued 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 nonintegral fractal dimension is chaotic. However, it has been found that this assumption may not he a sufficient indication of deterministic chaos. One counterexample is the socalled 1/ fP processes. Such processes have spectral density S(f ) f (2H+1) (22) 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 ahrl : 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 deterministiclike. 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 lowdimensional 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 [2327], estimation of the strength of measurement noise in experimental data [28, 29], pathological tremors [30], shearthickening 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), (.III11.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 [3436]: 1K = [x(i), x(i + L), ..., x(i + (m 1)L)] (23) 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 +kV (24) 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. 24 with Eq. 21, one notices that %~+k Vy,  phI i4 the role of d(t), while %~ Vy1  phi4 the role of d(0). Figure 21(a) shows the A(k) vs. k curves for the chaotic Lorenz system driven by stochastic forcing: dx =16(x y) + Drll(t) = xz + 45.92x y + Drl2(t) (25) = xy 4z + Dr/3( with < rli(t) > 0, < i6i,(')> b(tt'), 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. 21(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. 21(b) for the noisy Lorenz system with D = 4. We observe from Fig. 21(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 sniallscale dynamics is stronger than that on the largescale 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 MackeyGlass 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) (26) 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 AlackeyGlass "tI 3I O 50 200 200 Timedependent 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, 2i/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 21. 100 150 0 50 100 150 Evolution timne k 3.5  2.5 1.5 0.5 0 100 200 300 400 500 600 Evolution time k Figure 22. Timedependent exponent A(k) vs. evolution time k curves for the MackeyGlass system. Nine curves, from bottom to top, correspond to shells (2(i+1)/2, 2i/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. 22, 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 nonuniform growth rate in highdimensional systems . 2.2 General Computational Framework for Powerlaw Sensitivity to Initial Conditions (PSIC) To understand the essence of powerlaw sensitivity to initial conditions (PSIC), let us first consider the 1D case and consider Eq. 27, where Ax(0) is the infinitesimal discrepancy in the initial condition, Ax(t) is the discrepancy at time t > 0. ((t) =lim (27) When the motion is chaotic, then ((t) = exlt (28) ((t) satisfies Tsallis and coworkers [7, 11] have generalized Eq. 29 to di = A,((t)V (210) where q is called the entropic index, and X, is interpreted to be equal to K,, the generalization of the K~olmogorovSinai entropy. Eq. 210 defines the PSIC in the 1D case. Obviously, PSIC reduces to ESIC when q 1 The solution to Eq. 210 is s'(t) =[1 +(1 q)Agt] 'il) (211) When t is large and q / 1, ((t) increases with t as a powerlaw, ((t) Ct 7 ')(212) where C = [(1 q) A,] 1 (V). For Eq. 212 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. 23. Eq. 211 can then be generalized to highdimensional 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 qLyapunov exponent, corresponding to the powerlaw 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~olmogorovSinai entropy is the sum of all the positive Lyapunov exponents. We conjecture a similar relation may hold between the qLyapunov 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. 213 again gives a powerlaw 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 < %~ g1 < r2, Where rl and T2 arT ClOSe. These arguments so~ 1 that the timedependent exponent curves defined by Eq. 24 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) (214) 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 noisefree and noisy logistic and Henon maps. 2.3 Characterizing Edge of Chaos by Powerlaw 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, (215) 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 transientfree time series. In Figs. 23(ac), we have plotted the A(t) vs. t curves for parameter values al, a,, and a2, TOSpectively. We observe from Fig. 23(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 coworkers [11] that at the edge of chaos for the logistic map, ((t) is given by Eq. 211 with q m 0.2445. Surprisingly, we do not observe such a divergence in Fig. 23(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. 24(a)). The more interesting pattern is the one that is shown in Fig. 23(c), where we observe a linearly increasingly A(t) vs. t curve. In fact, shown in Fig. 23(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 6 (d) 1 05 6 (e) O 4 1 O 0 1 23 456 In Timedependent exponent A(t) vs. t curves for time series generated from the noisefree 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 (df) 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 (cf) 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 23. 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. 23(c) are less smooth than those reported earlier [2123, 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 2.0 1.5  0.5 0 50 100 150 200 250 2.0 ' 1.5 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. 23(df), 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. 23(df) 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~ Vy1 < 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 104 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. 213 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(df) 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 noisefree and noisy nmap for parameter values listed above, and found very similar results to those presented in Fig. 2:3. In Figs. 24(a,h), we have plotted A(t) vs. In t curves for the noisefree and noisy nmap for a,. = 1.:386:37288 b = 0.01. As can he observed clearly, Fig. 24(b) is very similar to Figs. 2:3(df), while the slope for the curves in Fig. 24(a) is much smaller than that in Fig. 24(b). Again, we have observed (but not shown by figures here, since they are very similar to Figs. 2:3(df) and Fig. 24(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: noisefree and noisy logistic and Henon maps near the edge of chaos, we have found that when there is no 2.0 (a)''""'""''"" 6 (b 1.5 S0.5  0.0 1 2 0.5 1 1 1.0O ~......... ......... ......... ......... .......... ............ O . 0123456 0123456 Int Int Figure 24. Timedependent exponent A(t) vs. t curves for time series generated from (a) the noisefree 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. CHAPTER 3 CHARACTERIZING RANDOM FRACTALS BY POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) 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 longrangecorrelations 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 Powerlaw Sensitivity to Initial Conditions (PSIC) A continuoustime stochastic process, X = {X(t), t > 0}, is said to be selfsimilar if X(At) d AHX(t), t > 0 (31) for A > 0, O < H < 1, where u~~denotes equality in distribution. H is called the selfsimilarity parameter, or the Hurst parameter. Before proceeding on, we note that, more rigorously speaking, processes defined by Eq. 31 should be called selfaffine processes [43] instead of selfsimilar 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 selfsimilar, to follow the convention in some of the engineering disciplines. The following three properties can be easily derived from Eq. 31: E[X(At)] E [X(t)] Mean (32) Var [X(At)] Var[X(t)]= 2HVariance (33) R, (t, s) = 2HAutocorrelation (34) Note that Var[X(t)] = R,(t, t); hence, Eq. 33 can be simply derived from Eq. 34; we have listed it as a separate equation for convenience of future reference. If we consider only secondorder statistics, or if a process is Gaussian, Eq. 32 to Eq. 34 can be used instead of Eq. 31 to define a selfsimilar process. A very useful way of describing a selfsimilar 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 (35) otherwise and take the Fourier transform of X(t, T): F(f, T) = lX(t)e2xjftd (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) T>oo Noting that the PSD for AHX(At) is A2H limX(t2ft _2S/) T>oo AT$ I 'o X.X)ZiLt 2 X"sfX and that the PSD for AHX(At) and X(t) are the same [44], we have S(f) = S(f/A)A2H1 The solution to the above equation is S( f) ~ f P (37) with /7 = 2H + 1 (38) Processes with power spectral densities as described by Eq. 37 are called 1/ f processes. Typically, the powerlaw 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 [4850], DNA sequence [5153], 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 powerlaw 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 powerlaw 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 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 A1/2B(At) have the same distribution. Thus, we see from Eq. 31 that Bm is a selfsimilar 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 zeromean and unitvariance 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 22 2. D1 is simply the third term of Eq. 311, and the coefficient 22 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 23 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 = S1/2 aS s 0. Although B(t) is almost surely not differentiable in t, symbolically one still often writes B(t) = w~(T)dv (312) 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 RiemannStieltjes integral sense. This has profound consequence when numerically integrating Eq. 312. 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 severely underestimating the strength of the noise. (a) White noise (b) Brownian motion Figure 31. 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. 313 becomes B(t) = w~7(i) (3 14) i= 1 Eq. 314 provides perhaps the simplest method of generating a sample of Bm. An example is shown in Fig. 31. Eq. 314 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 unitstep function 1t>0 P(t) = O < (315) 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) (316) 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. 31 for a selfsimilar stochastic process. Fractional 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. 32 to Eq. 34, the above three properties completely determine its selfsimilar character. Fig. 32 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"te~d 7 It is easy to verify that fBm satisfies the following scaling property: P(BHt + ) BH~t I) = P 8H(X + ) BXH XHX} (320) In other words, this process is invariant for transformations conserving the similarity variable X/tH (c) H=0.75 (d) H=0.90 Figure 32. Several fBm processes with different H. The integral defined by Eq. 319 diverges. The more precise definition is BH t) BH(0) =r( K1 7dB0 (321) where the kernel K(t r) is given by [63] K(t r) =l l (322) 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. 318, one obtains the autocovariance function y(k) for the fGn process: q~~~~k)2` =. E(4ik/(X)= ( )2H 2k2H + k 12H, 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) = 222 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 wellknown 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 antipersistent and persistent correlations [43], respectively. 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~ 12H = k~2H[(1 + 1/k)2H + (1 1/k~)2H] = k~2H[2 + 2H(2H 1)k2] hence, y(k) ~ H(2H 1)k2H2, Sk 00(324) 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 (325) Therefore , dt)= HtH1 (326) Expressing t in terms of (, we have = H( $ H327) Comparing with the defining equation of PSIC (Eq. 210), we find that 1 2 q =1 1 (328) H P1 X( ) = H= (329) Noticing 0 < H < 1, from Eq. 328, we have oo < q < 0. Computationally, the key is to demonstrate behaviors described by Eq. 325. This can be readily done by calculating the time dependent exponent curves defined by Eq. 214 (see also Eq. 24). In Fig. 33 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 Powerlaw 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. 34 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 heavytailed distribution. Time Figure 34. ON/OFF sources with NV = 3, X (t), X2(t), X3(t), and their summation S3 ) SON OFF ON OFF ON OFF H = 0.25 H = 0.50 H = 0.75 ___  6.5 2 2.5 3 3.5 4 Ink Figure 33. Timedependent exponent A(k) vs. In k curves for fractional Brownian motion. I I X ,(t) I I I I X 2t) X3,t S3(t A heavytailed distribution is commonly expressed as f (x) ~ x"1, x 00 Equivalently, one may write: P X > x] ~ x". x 0 o 330) This expression is often called the complementary cumulative distribution function (CCDF) of a ]!. l.ir I.l:d distribution. In fact, it is more popular for expressing a heavytailed distribution, since it emphasizes the tail of the distribution. It is easy to prove that when p < 2, the variance and all higher than 2ndorder moments do not exist. Furthermore, when p < 1, the mean also diverges. When the powerlaw relation extends to the entire range of the allowable x, we have the Pareto distribution: P[X J > ~ I x]= >biU > 2> 2> (331) 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 (332) 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. 35 we have shown a few examples for p = 1.2, 1.6, and 2.0. Noticing Eq. 332, 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. (325) and (328) 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 2.8 3.4 .53 3.5 4 4.5 Ink Figure 35. Timedependent 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 powerlaw 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 processes, The latter equation thus provides another simple means of estimating the Hurst parameter. 3.2 Characterizing Levy Processes by Powerlaw Sensitivity to Initial Conditions (PSIC) We now consider Levy processes, another important type of random fractal models that have found numerous applications [6570]. 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 powerlaw 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 acstable 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 acstable distribution with characteristic function e(ts) 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. 333, 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 heavytailed, 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 acstable Levy process. The key is to simulate its increment process, which follows an acstable 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+(a1)8] 1a/ 34 (cos000cos8)1/a 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 (335) yZ+( S"qln n+6 a~= 1 where Z is given by Eq. 334. We have simulated a number of symmetric acstable Levy motions with different a~. Figure 36 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 flightslike 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 Levyflightslike patterns have been used as a type of screen saver for Linux operating systems. The symmetric acstable Levy flight is 1/a~ selfsimilar. That is, for c > 0, the processes {L,(ct), t > 0} and {cl/oL,(t), t > 0} have the same finitedimensional distributions. By this argument as well as Eq. 333, it is clear that the length of the (a) (b) (c) (d) (e) (f) Figure 36. Examples for Brownian motions and Levy motions. 1D and 2D Brownian motions (a,b) vs. Levy motions with a~ = 1.5 and 1 for (c,d) and (e,f), respectively. o a =1.0 4.5 0 a =1.5 5 5.5 6 06~.5 1 1.5 2 2.5 3 Ink Figure 37. Timedependent 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 (336) 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~ (338) A 0 = 1/al (339) Noticing 0 < a~ < 2, from Eq. 338, 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. 37. The slopes of the straight lines correctly estimate the values of a~ used in the simulations. CHAPTER 4 STUDY OF INTERNET DYNA1\ICS BY POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) 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 powerlaw sensitivity to initial conditions (PSIC) provides an interesting framework for time series analysis. We illustrate this by studying the complicated dynamics of Internet transport protocols. The Internet is one of the most complicated systems that man has ever made. In recent years, two types of fascinating multiscale behaviors have been found in the Internet. One is the temporaldomain 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 spatialdomain scalefree 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 widearea connections. Typically, the transport dynamics over the Internet connections are a result of the nonlinear dynamics of the widelydeploi II Transmission Control Protocol (TCP) interacting with the Internet traffic, which is stochastic and often selfsimilar. 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 ns2 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 windowsize 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 steadystate analysis and study of convergence to a steadystate 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" nonconvergent transport dynanxics. Although recently a lot of research has been carried out [8288] 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 1D 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 [3436 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 quasiperiodic motions (Figs. 41(ad)). When NV is large, such as NV = 17, chaos is observed (see Fig. 42(a,b)). We thus conclude that one route to chaos in TCP dynamics is the wellknown quasiperiodic route. In order to show more diversified dynamical behavior such as quasiperiodic 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 quasiperiodic W(i) is shown in Fig. 41(a), where C = 0.1 Mbps, delay d = 10 ms, B = 10, NV = 3. Its power spectral density (PSD) is shown in Fig. 41(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. 41(a). The start of an exponential backoff indicates triggering of a loss episode or heavy congestion. Hence, the T(i) is related to the wellknown 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. 41(a) is shown in Fig. 41(c), together with its 1015 1010 0 105 100 10 0 0.05 0.1 Frequency 10'" (d) 0 1000 2000 3000 4000 5000 index i 50 index i Frequency 105 100 105 0 0.1 0.2 0.3 0.4 0.5 Frequency Figure 41. Example of quasiperiodic 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. 41(d). We observe that this T(i) time series is a simple periodic sequence, and hence, the IT(i) time series is quasiperiodic. However, T(i) time series can still be quasiperiodic. An example is shown in Fig. 41(e), with C 0.5 Mbps, delay d ms, B = 10 packets, NV = :3, together with its PSD is Fig. 41(f). Note that there are still many discrete sharp peaks in Fig. 41(f). To appreciate how many independent frequencies 100 90 80 0 50 100 150 index i Iruilu~ may be contained in the T(i) time series of Fig. 41(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 quasiperiodic with four independent frequencies. Before we proceed to discuss the chaotic TCP dynamics, we note that simple periodic or quasiperiodic 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 quasiperiodic with two independent frequencies. Furthermore, the probability distributions for the constant, periodic, and quasiperiodic T(i) time series are only composed of a few discrete peaks. This II__ ; that the observed quasiperiodic T(i) time series constitutes a discrete torus. Let us now turn to the discussion of chaotic TCP dynamics. Shown in Figs. 42(a,c) are two irregular T(i) time series. The parameters C, d, and B used for Fig. 42(a) are the same as those for Figs. 41(ad). Noting from Fig. 42 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 illsuited for the study of the underlying irregular dynamics. To show that the T(i) time series in Figs. 42(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. 42(a,c) are shown in Figs. 42(b,d), where the four curves, from bottom to top, correspond to shells of sizes (2(i+1)/2, 2i/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. 42(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. (c) 5 1 5 x14 8 5 2 15 xr 101 (d) 5 4 3 2 0 0 0 2000 4000 6000 0 10 20 30 40 50 Figure 42. 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. 42(a,c)? They are powerlawlike, P(T < t) ~ t Y, for almost two orders of magnitude in t, as shown in Figs. 43(a,b). The exponent y is not a universal quantity, however. We now develop a simple 1D 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 103~ 103 141 14 102 103 104 105 102 103 104 105 t t Figure 43. Complementary cumulative distribution functions (CCDFs) for the two chaotic T(i) time series of Figs. 42(a,c). Flow control ensures that the sender sends no more than the size of the receiver's last advertised flowcontrol 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 flowcontrol window ( fwnd) advertised by the receiver to the sender and a congestioncontrol 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 wup~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. 41) 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. 42 by assuming Wmax to be periodic and quasiperiodic, 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. 42(a,c) an artifact of Tahoe version of TCP S10 400 200 0 CIVYUU yIYYUU"UU VRl V'IIIIYUII IIIVV 111 10 0 50 100 150 200 250 300 350 100 101 102 103 Figure 44. The time series T(i) extracted from a W(i) time series collected using net100 instruments over ORNLLSU connection. (a) the time series T(i), (b) the complementary cumulative distribution function (CCDF) for the T(i) time series. we used or more generally of the ns2 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. 44(a), whose profile is qualitatively quite similar to Fig. 42(a). In fact, it also follows a (truncated) powerlaw, as shown in Fig. 44(b). The exponent for the powerlaw part is even smaller than the time series of Fig. 42. The truncation in the powerlaw 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 powerlawlike 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 ns2 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 nonlinear behavior of TCP's Additive Increase and Multiplicative Decrease (AIMD) method emploix & for congestion control. The stochastic component is due to the often selfsimilar 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 windowsize, denoted by cwnd, in response to acknowledgments, and multiplicative decreases in response to inferred losses (ignoring the initial slowstart 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 chaoslike 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 [9396]. 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 Robbins1\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 analysis. 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 highbandwidth (OC192 at 10Gbps) with relatively low backbone traffic and a roundtrip 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 triptime 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 ORNLGaTech traces, and we only discuss the results for the latter. To appreciate better what these data look like, we have plotted in Figs. 45(ad) segments of four datasets for ORNLGaTech 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 [3436, 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 [3436] 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 nontrivial deterministic structure in the data. (ii) For periodic signals, A(k) is essentially zero for any k. x 104 4 rl qX104 (b) 1 TCP source  3 1 0 0 2000 4000 6000 2000 4000 6000 x 104 3~ i X 104 41(d) 2 TCP sources [Hdi n8 n, Ellrfli rfl TCP sources 1 2000 4000 sample n 6000 2000 4000 sample n 6000 Figure 45. Time series for the congestion window size cwnd for ORNLGaTech connection. (iii) For quasiperiodic 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. 45 are plotted in Figs. 46. 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 3 1 0 0 50 100 150 200 250 4 (c 12 TC urces 3 2 1 0 50 100 150 200 250 50 100 150 200 Evolution time k 50 100 150 200 Evolution time k 250 Figure 46. Timedependent exponent A(k) vs. k curves for cwnd data corresponding to Fig. 45. 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, 2i/2) i = 2, 3, ,8. top, correspond to shells of sizes (2(i+1)/2, 2i/2), i 2, 3, 8. We make the following observations: (i) The dynamics are complicated and cannot be described as either periodic or quasiperiodic 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 welldefined 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. 46(a) with only one TCP source increases much slower when k just exceeds (m1)L. Also important is that the A(k) curves in Figs. 46(c,d) are much smoother than those in Figs. 46(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 powerlaw sensitivity to initial conditions (PSIC). Figure 47 shows the A(k) vs. In k curves for ORNLGaTech 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 windowsize 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 Timedependent exponent A(k) vs. In k curves for cwnd data corresponding to Fig. 45. 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, 2i/2) i = 2, 3, ,8. Figure 47. CHAPTER 5 STITDY OF SEA CLITTTER RADAR RETITRNS BY POWERLAW SENSITIVITY TO INITIAL CONDITIONS (PSIC) 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 [102108]. 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_tr 1.. that the nonchaotic 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 lowpass 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 [102106] are conducted by comparing measured sea clutter data with simulated stochastic processes. We can ask: can the nonchaotic nature of sea clutter he directly demonstrated without resorting to simulated stochastic processes'? Recognizing that simple lowpass 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 lowdintensional 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 scaledependent 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 scaledependent 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 powerlaw 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 lowdintensional chaos developed by Gao and Zheng [21, 22]. And we show that the scaledependent 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 2530 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. 51. 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. 52 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. 53(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 Target BB..B. Bi B.+ ... B, 1 2 11 1 111 (Secondary)(Primary)(Secondary) I~RRi Figure 51. Collection of sea clutter data 25 20 o 15 E 10 5 (b) E 1 0 20 rrl Figure 52. 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 lowdimensional 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 clipping. 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], shearthickening 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. 54(a, c, e) and the curves for the data with targets shown in Figs. 54(h, d, f), respectively. We have simply chosen m = 6 and 4 3 4r 3 r2 1 0 4 3 1 (a) 5 10 15 20 25 30 (b) 5 10 15 20 25 30 (c) 5 10 15 20 25 30 (d) 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 54. Examples of the timedependent 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, 2i/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. 54 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, lowflying 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 nonGaussian, 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], lognormal [113], K( [114, 115], and compoundGaussian distributions [116], as well as using chaos theory [100104, 107, 108, 110, 111], wavelets [117], neural networks [118, 119], waveletneural 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 scaledependent exponent may be helpful for target detection within sea clutter. 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. 24, subject to the conditions that XiXj < 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 leastsquares fit to the first few samples of the A(k) curve where the A(k) curve increases linearly or quasilinearly. 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. 55(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. 55(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 scaledependent exponent corresponding to large scale may be useful for detecting targets within sea clutter data. By scaledependent exponent, we mean that if we use the leastsquares fit to the linearly or quasilinearly 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 scaledependent exponent has been used in the study of noiseinduced 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 56(a)(j) show the variations of the scaledependent exponent corresponding to large scale with the 14 range hins for the Primary JTarget Bin &z 05 0 5 10 15 (d) 05 0 5 10 15 1x (f)~* 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 55. 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. 55. We observe that the primary target bin can be easily separated from the range bins without the target, since the scaledependent exponent for the primary target bin is much larger than those for the bins without the target. It is thus clear that the scaledependent 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 scaledependent 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 56. Variations of the scaledependent 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 scaledependent 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 automatic. 5.4 Target Detection within Sea Clutter by Powerlaw 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. 57(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. 58 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. 58 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.5' 0 0.5 1 1.5 2 2.5 3 3.5 In k In k Figure 57. 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. Primary Target Bin 0.9 co. 0.85 0.75 x 0.7' 0 2 4 8 Bin Number 10 12 14 Figure 58. 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. 59. We observe that the frequencies completely separate for the HH datasets. This means the detection accuracy can be 1011' . no targets 12 primary targets 10 a, LL n mm t0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 59. Frequencies of the bins without the HH datasets. targets and the bins with primary targets for CHAPTER 6 1\ULTISCALE ANALYSIS BY SCALEDEPENDENT LYAPUNOV EXPONENT (SDLE) In C'!. Ilter 1, we have emphasized the importance of developing scaledependent 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 scaledependent Lyapunov exponent (SDLE), and show that the SDLE can readily classify various types of complex motions, including truly lowdimensional chaos, noisy chaos, noiseinduced chaos, processes defined by powerlaw 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 scaledependent 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 rfold time by monitoring the divergence between a reference trajectory and a perturbed trajectory. To do so, it needs to define 1.. I.rest 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 nondifferentiable 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 nondifferentiable. Next we examine the relation between et and et+nt, where at is small. When at 0 we have, et+nt = etex(tt)nt (61) where A(et) is the SDLE. It is given by In et+nt In et A(et) = (62) 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(63) 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 socalled timedependent 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. 62 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 highdimensional systems, this is not true, especially when the growth rate is nonuniform and/or the eigenvectors of the Jacobian are nonnormal. 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 L65 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 (65) 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 Noiseinduced Chaos Obviously, for truly lowdimensional 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. 25. Figure 61(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 eentropy. In fact, such a relation can he readily proven for the eentropy. 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 noiseinduced 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 (66) 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 ~D=4 0.5 v 0 1.5 1 0.5 2 1.5 1 0.5 0 log ,E log ,E Figure 61. Scaledependent Lyapunov exponent A(e) curves for (a) the clean and the noisy Lorenz system, and (b) the noiseinduced 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, noiseinduced chaos occurs, and thought that it may be difficult to distinguish noiseinduced chaos from clean chaos. In Fig. 61(b), we have plotted the A(et) for this particular noiseinduced chaos. The computation was done with m = 4, L = 1 and 10000 points. We observe that Fig. 61(b) is very similar to the curves of noisy chaos plotted in Fig. 61(a). Hence, noiseinduced 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. 61 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. 61. 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 35 5 10 100 Figure 62. Scaledependent Lyapunov exponent A(e) curve for the MackeyGlass system. The computation was done with ni = 5, L = 1, and 5000 points sampled with a time interval of 6. Finally, we consider the MackeyGlass delay differential system [39], defined by Eq. 26. 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. 62, 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 nonuniform growth rate in highdimensional 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 Powerlaw Sensitivity to Initial Conditions (PSIC) In C'!s Ilter 3, we have seen that powerlaw 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 SDLE. 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. 62, we find that (()=In $t+nt In (a(67 When at 0, 1+ (1 q)Agt > (1 q)XAft. Simplifying Eq. 67, we obtain A((s) = Aq tq1 (68) 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. :338) and A,"l = 1/0~ (Eq. :339). 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 (69) 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) (610) (2 + A) y (1 + a) if y E (1/2, 1] The map Fly) is shown in Fig. 63 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. 64(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. 64(b). It is interesting to observe that the chaotic behavior can he well resolved by only a few hundred points. However, the (30.6 S0.4 LL. 1 T 0.2' 0 Figure 63. The function F(x) (Eq. 610) for a = 0.4 is shown as the dashed lines. The function G(x) (Eq. 611) is an approximation of F(x), obtained using 40 intervals of slope 0. In the case of noiseinduced chaos discussed in the paper, G(x) is obtained from F(x) using 104 HinerValS Of Slope 0.9. 100 102 100 102 clean data noisy data 104 103 102 10 1 104 103 102 Scaledependent Lyapunov exponent A(e) for (a) 5000 points were used, for the noisy case, used. Figure 64. the model described by Eq. 69. 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. 64(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 (69) can be modified to give rise to an interesting system with noiseinduced chaos. This can be done by replacing the function Fly) in map (69) by G(y) to obtain the following map [128], xt+l = [Xt] + G(xt [xt]) + arlt (611) 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. 610. An example of G(y) is shown in Fig. 63. 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 (611) is nonchaotic, since the largest Lyapunov exponent In(0.9) is negative. With appropriate noise level (e.g., a = 104 or 103), the SDLE for the system becomes indistinguishable to the noisy SDLE shown in Fig. 64 for the map (69). Having a diffusive regime on large scales, this is a more complicated noiseinduced 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 (611) without noise. A transientfree 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 ext, 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 socalled 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 welldefined spectral peak, indicating oscillatory motions. An example is the chaotic Lorenz system described by Eq. 25. In Figs. 65 (ac), we have shown the powerspectral 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. 66). 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 " 104 102 104~ 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) 106 ,8100 4 104 1021 104 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 104 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 65. P~owerspectral density (PSD>) for (a) x(t), (b) y(t), (c) z(t), (d) ei)(t), (e) from a farinfrared laser in a chaotic state. The laser dataset can be downloaded from the link http://wwwpsy ch. stanf ord .edu/~ andreas/Time Ser ies /Sant aFe .htm1. It contains 10000 points, sampled with a time interval of 80 ns. Fig. 67(a) shows the laser dataset. The PSD of the data is shown in Fig. 67(b), where one observes a sharp peak around 1.7 MHz. Fig. 67(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. 67(a). 