Cross Recurrence Plot Based Synchronization of Time Series

8 Pages·0.2 MB·English

by  N. Marwan

2
0
0
2

n
a
J

9
Cross Recurrence Plot Based Synchronization of Time Series
2

]
h
p
- N.Marwan1, M.Thiel1,and N.R. Nowaczyk2
o
e
1InstituteofPhysics,UniversityofPotsdam,Germany
g
. 2GeoForschungsZentrumPotsdam,Germany
s
c
i
s
y
h
p
[

Camera-ready Copy for
1
v
NonlinearProcesses inGeophysics
2
6
Manuscript-No. 1040
0
1
0
2 Offsetrequests to:
0
/ N.Marwan
s
c InstituteofPhysics
si UniversityofPotsdam
y
14415Potsdam
h
p Germany
:
v
i
X
r
a
NonlinearProcessesinGeophysics(20**)*:101–107
Nonlinear Processes
in Geophysics
©EuropeanGeophysicalSociety20**
Cross Recurrence Plot Based Synchronization of Time Series
N.Marwan1,M.Thiel1,andN.R.Nowaczyk2
1InstituteofPhysics,UniversityofPotsdam,Germany
2GeoForschungsZentrumPotsdam,Germany
Received:???–Accepted:???
Abstract. Themethodofrecurrenceplotsisextendedtothe cross recurrence plots (CRP), we have found an interesting
cross recurrence plots (CRP), which among others enables featureofit. Besidesthepossibilityoftheapplicationofthe
thestudyofsynchronizationortimedifferencesintwotime recurrencequantificationanalysis(RQA)ofWebberandZbi-
series. This is emphasized in a distorted main diagonal in lutonCRPs(1998),thereisamorefundamentalrelationbe-
thecrossrecurrenceplot,thelineofsynchronization(LOS). tweenthestructuresintheCRPandtheconsideredsystems.
A non-parametrical fit of this LOS can be used to rescale Finally,thisfeaturecanbeusedforthetaskofthesynchro-
the time axis of the two data series (whereby one of it is nizationofdatasets. Althoughthefirststepsofthismethod
e.g.compressedorstretched)sothattheyaresynchronized. are similar to the sequence slotting method, their roots are
Anapplicationofthismethodtogeophysicalsedimentcore different.
dataillustratesitssuitabilityforrealdata.Therockmagnetic First we give an introductionin CRPs. Then we explain
dataoftwodifferentsedimentcoresfromtheMakarovBasin the relationship between the structures in the CRP and the
can be adjusted to each other by using this method, so that systems and illustrate this with a simple model. Finally we
theyarecomparable. apply the CRP on geophysicaldata in order to synchronize
variousprofilesandtoshowtheirpracticalavailability.Since
wefocusonthesynchronizationfeatureoftheCRP,wewill
notgiveacomparisonbetweenthedifferentalignmentmeth-
1 Introduction
ods.
The adjustmentofdata sets with varioustime scales occurs
inmanyoccasions,e.g.datapreparationoftreeringsorgeo-
2 TheRecurrencePlot
physicalprofiles. In geology,often a large set of geophysi-
cal data series is gained at various locations (e.g. sediment
Recurrence plots (RP) were firstly introduced by Eckmann
cores). Thatiswhythesedataserieshaveadifferentlength
etal. (1987)in orderto visualizetime dependentbehaviour
andtimescale.Beforeanytimeseriesanalysiscanbestarted,
of orbits~x in the phase space. A RP represents the recur-
the data series have to be synchronized to the same time i
renceofthephasespacetrajectorytoastate. Therecurrence
scale. Usually, this is done visually by comparingand cor-
ofstatesisafundamentalpropertyofdeterministicdynami-
relating each maximum and minimum in both data sets by
calsystems(Argyrisetal.,1994;Casdagli,1997;Kantzand
hand(wigglematching),whichincludesthehumanfactorof
Schreiber, 1997). The main step of the visualization is the
subjectivenessandisalengthyprocess. Anautomaticaland
calculationoftheN ×N-matrix
objectivemethodforverificationshouldbeverywelcome.
Inthelastdecadessometechniquesforthiskindofcorre-
R =Θ ε−k~x −~x k , i,j =1...N, (1)
lation andadjustmentwere suggested. Theyspan graphical i,j i j
(cid:0) (cid:1)
methods (Prell et al., 1986), inverse algorithms, e.g. using
whereεisapredefinedcut-offdistance,k·kisthenorm
Fourierseries(Martinsonetal.,1982)andalgorithmsbased
(e.g. the Euclidean norm) and Θ(x) is the Heaviside func-
onsimilarityofdata,e.g.sequenceslotting(Thompsonand
tion. The values one and zero in this matrix can be simply
Clark,1989).
visualizedbythecoloursblackandwhite. Dependingonthe
However,we focuson a methodbased onnonlineartime
kind of the application, ε can be a fixed value or it can be
series analysis. During our investigationsof the method of
changedforeachiinsuchawaythatintheballwiththera-
Correspondenceto: N.Marwan diusεapredefinedamountofneighboursoccurs. Thelatter
101
102
will providea constantdensity of recurrencepointsin each
columnoftheRP.
6
Therecurrenceplotexhibitscharacteristicpatternsfortyp-
icaldynamicalbehaviour(Eckmannetal.,1987;WebberJr.
and Zbilut, 1994): A collectionof single recurrencepoints, 5
homogeneously and irregularly distributed over the whole
plot, reveals a mainly stochastic process. Longer, parallel
4
diagonals formed by recurrence points and with the same
g
distance between the diagonals are caused by periodic pro- n
cesses. A palingoftheRP awayfromthemaindiagonalto e i 3
m
the corners reveals a drift in the amplitude of the system. Ti
Vertical and horizontal white bands in the RP result from
2
states which occur rarely or represent extreme states. Ex-
tendedhorizontalandverticalblacklinesorareasoccurifa
statedoesnotchangeforsometime,e.g.laminarstates. All 1
thesestructureswereformedbyusingthepropertyofrecur-
rence of states. It should be pointed out that the states are
onlythe“same”andrecurinthesenseofthevicinity,which 1 2 3 4 5 6
is determined by the distance ε. RPs and their quantitative Time in f
analysis(RQA)becamemorewellknowninthelastdecade
(e.g. Casdagli, 1997). Their applicationsto a wide field of Fig.1.Crossrecurrenceplotsofsinefunctionsf(t)=sin(ϕt)andg(t)=
miscellaneousresearchshowtheirsuitabilityintheanalysis sin(ϕt+asin(ψt)),whereata = 0fortheblackCRP,a = 0.5forthe
greenCRPanda = 1fortheredCRP.Thevariationinthetimedomain
ofshortandnon-stationarydata.
leadstoadeformingofthesynchronizationline.
3 TheCrossRecurrencePlot
dendrochronology. Sedimentcoresmighthaveundergonea
numberofcoringdisturbancessuchascompressionorstretch-
Analogous to Zbilut et al. (1998), we have expanded the
ing. Moreover,coresfromdifferentsiteswithdifferingsed-
method of recurrenceplots (RP) to the method of cross re-
imentation rates would have different temporal resolutions.
currenceplots. IncontrasttotheconventionalRP,twotime
Allthesefactorsrequireamethodofsynchronizing.
seriesaresimultaneouslyembeddedinthesamephasespace.
ACRPofthetwocorrespondingtimeserieswillnotcon-
Thetestforclosenessofeachpointofthefirsttrajectoryx
i
tain a main diagonal. But if the sets of data are similar
(i = 1...N) with each point of the second trajectory y
j
e.g.onlyrescaled,amoreorlesscontinuouslineintheCRP
(j =1...M)resultsinaN ×M array
that is like a distorted main diagonal can occur. This line
CR =Θ ε−k~x −~y k . (2) containsinformationontherescaling.Wegiveanillustrative
i,j i j
example.ACRPofasinefunctionwithitself(i.e.thisisthe
(cid:0) (cid:1)
Itsvisualizationiscalledthecrossrecurrenceplot(CRP). RP)containsamaindiagonal(blackCRPinFig.1). Hence,
Thedefinitionoftheclosenessbetweenbothtrajectoriescan theCRPsintheFig.1arecomputedwithembeddingsofdi-
be varied as described above. Varying ε may be useful to mension one, further diagonal lines from the upper left to
handlesystemswithdifferentamplitudes. thelowerrightoccur. Theselinestypifythesimilarityofthe
TheCRPcomparestheconsideredsystemsandallowsus phase space trajectoriesin positiveand negativetime direc-
to benchmark the similarity of states. In this paper we fo- tion.
cusonthebowed“maindiagonal”intheCRP,becauseitis Nowwerescalethetimeaxisofthesecondsinefunction
relatedtothefrequenciesandphasesofthesystemsconsid- inthefollowingway
ered.
sin(ϕt)−→ sin ϕt+asin(ψt) (3)
(cid:0) (cid:1)
4 TheLineofSynchronizationintheCRP We will henceforth use the notion rescaling only in the
mentionofthe rescalingofthe timescale. Therescalingof
Regarding the conventional RP, Eq. 1, one always finds a the secondsine functionwith differentparametersϕresults
maindiagonalintheplot,becauseoftheidentityofthe(i,i)- in a deformation of the main diagonal (green and red CRP
states. The RP can be considered as a special case of the inFig.1). Thedistortedlinecontainstheinformationonthe
CRP,Eq.2,whichusuallydoesnothaveamaindiagonalas rescaling,whichwewillneedinordertore-synchronizethe
the(i,i)-statesarenotidentical. two time series. Therefore, we call this distorted diagonal
In data analysis one is often faced with time series that lineofsynchronization(LOS).
are measured on varying time scales. These could be sets Inthefollowing,wepresentatoyfunctionto explainthe
from borehole or core data in geophysics or tree rings in procedure. If we consider a one dimensional case without
103
embedding,theCRPiscomputedwith In the Appendix we describe a simple algorithm for esti-
mating the LOS. Its determination will be better for higher
CR(t1,t2)=Θ ε−kf(t1)−g(t2)k . (4) embeddings because the vertical and cross-diagonal struc-
(cid:0) (cid:1)
Ifwesetε = 0tosimplifythecondition,Eq.(4)delivers tureswillvanish. We donotconcealthatthe embeddingof
arecurrencepointif the time series is connected with difficulties. The Takens
EmbeddingTheoremholdsforclosed,deterministicsystems
f(t1)=g(t2). (5) without noise, only. If noise is present one needs its real-
ization to finda reasonableembedding. For stochastictime
Ingeneral,thisisanimplicitconditionthatlinksthevariable
series itdoesnotmake sense to considera phasespace and
t1 tot2. Consideringthephysicalexamplesofabove,itcan
soembeddingisingeneralnotjustifiedhereeither(Romano,
beassumedthatthetimeseriesareessentiallythesame–that
tobepublished;Takens,1981).
meansthatf =g–uptoarescalingfunctionoftime. Sowe
The choice of a special embedding lag could be correct
canstate
foronesection ofthe data, butincorrectforanother(foran
f(t1)=f φ(t1) . (6) exampleseebelow). Thiscanbethecaseifthedataisnon-
(cid:0) (cid:1) stationary. Furthermore,thechoiceofmethodofcomputing
Ifthefunctionsf(·)andg(·)arenotidenticalourmethodis
theCRPandthethresholdεwillinfluencethequalityofthe
ingeneralnotcapableofdecidingifthedifferenceinthetime
estimatedLOS.
seriesisduetodifferentdynamics(f(·)6=g(·))orifitisdue
Thenextsectionswillbededicatedtoapplication.
tosimplerescaling. Sothe assumptionthatthedynamicsis
alike upto a rescalingin time is essential. Eventhoughfor
somecaseswheref 6=gitcanbeappliedinthesameway.If 5 ApplicationtoaSimpleExample
weconsiderthefunctionsf(·)=a·f¯(·)+bandg(·)=g¯(·),
whereby f(·) 6= g(·) are the observations and f¯(·) = g¯(·) Atfirst, weconsidertwosinefunctionsf(t) = sin(ϕt)and
arethestates,normalizationwithrespecttothemeanandthe g(t) = sin(ψt2), where the time scale of the second sine
standarddeviationallowsustouseourmethod. differs from the first by a quadraticterm and the frequency
ψ = 0.01ϕ. Sediment parameters are related to such kind
f(·)−hf(·)i
f(·) = a·f¯(·)+b−→f˜(·)= (7) of functions because gravity pressure increases nonlinearly
σ(f(·))
withthedepth. Itcanbeassumedthatbothdataseriescome
g(·)−hg(·)i
from the same process and were subjected to different de-
g˜(·) = (8)
σ(g(·)) positalcompressions(e.g.asquaredorexponentialincreas-
ingofthecompression). TheirCRP containsa bowedLOS
Withg¯(·)=f¯(·)thefunctionsf˜(·)andg˜(·)arethesameaf-
(Fig.2). Wehaveusedtheembeddingparametersdimension
terthenormalization.Thenourmethodcanbeappliedwith-
m=2,delayτ =π/2andavaryingthresholdε,sothatthe
outanyfurthermodification.
CRPcontainsaconstantrecurrencedensityof20%. Assum-
InsomespecialcasesEq.(6)canberesolvedwithrespect
ingthatthetimescaleofgisnotthecorrectscale,wedenote
to t1. Such a case is a system of two sine functions with that scale by t′. In order to determinethe non-parametrical
differentfrequencies
LOS, we have implemented the algorithm described in the
f(t) = sin(ϕ·t+α) (9) Appendix.Althoughthisalgorithmisstillnotmature,weob-
tainedreliableresults(Fig.3). Theresultingrescalingfunc-
g(t) = sin(ψ·t+β) (10)
tion has the expected squared shape t = φ(t′) = 0.01t′2
UsingEq.(5)andEq.(6)wefind (redcurveinFig.3).Substitutingthetimescalet′inthesec-
onddataseriesg(t′)bythisrescalingfunctiont=φ(t′),we
sin(ϕt1+α)−sin(ψt2+β)=0 (11) get a set of synchronized data f(t) and g(t) with the non-
andoneexplicitsolutionofthisequationis parametric rescaling function t = φ(t′) (Fig.4). The syn-
chronizeddataseriesareapproximatelythesame. Thecause
ϕ
⇒ t2 =φ(t1)= t1+γ (12) ofsomedifferencesisthemeanderingoftheLOS,whichit-
(cid:18)ψ (cid:19) selfiscausedbypartialweakembedding. Nevertheless,this
with γ = α−β . In this special case the slope of the main canbeavoidedbymorecomplexalgorithmforestimatingthe
ψ LOS.
lineinacrossrecurrenceplotrepresentsthefrequencyratio
and the distance between the axes origin and the intersec-
tion of the line of synchronizationwith the ordinatereveals 6 ApplicationtoRealData
the phase difference. The function t2 = φ(t1) (Eq.6) is a
transferorrescalingfunctionwhichallowsustorescalethe In order to continue the illustration of the working of our
secondsystemtothefirstsystem. Iftherescalingfunctionis methodwehaveappliedittorealdatafromgeology.
notlineartheLOSwillalsobecurved. In the following we compare the method of cross recur-
FortheapplicationonehastodeterminetheLOS–usually renceplotmatchingwith theconventionalmethodofvisual
non-parametrically–andthenrescaleoneofthetimeseries. wigglematching(interactiveadjustment). Geophysicaldata
104
Reference Data
1
f ( t ) 0
80
−1
0 10 20 30 40 50 60 70 80 90
Time t
70
Rescaled Data
1
60 φg ( ( t‘ ) ) 0
g
m −1
50 0 10 20 30 40 50 60 70 80 90
e
st Time t
y
s
e in 40 F(riegd.)4a.nRdeaffetreernc(beladcakta)sreesriceasli(nugppbeyrupsainnegl)thaendrersecsaclainlegdfduantcatisoenrieosfbFeifgo.r3e
m
Ti (lowerpanel).
30
oftwosedimentcoresfromtheMakarovBasin,centralArc-
20
tic Ocean, PS 2178-3 and PS 2180-2, were analysed. The
taskshouldbetoadjustthedataofthePS2178-3data(data
10 lengthN = 436)tothescaleofthePS2180-2(datalength
N = 251)in orderto get a depth-depth-functionwhich al-
lowstosynchronizebothdatasets(Fig.5).
10 20 30 40 50 60 Wehaveconstructedthephasespacewiththenormalized
Time in system f six parameters low field magnetic susceptibility (κLF), an-
hysteretic remanent magnetization (ARM), ratio of anhys-
Fig.2. Crossrecurrenceplotsoftwosinefunctionsf(t) = sin(ϕt)and tereticsusceptibilitytoκLF (κARM/κLF),relativepalaeoin-
g(t) = sin(ψt2))whichisthebaseforthedeterminationoftherescaling tensity(PJA),mediandestructivefieldofARM(MDF )
ARM
functionbetweenbothdataseries.Theembeddingparametersweredimen-
and inclination (INC). A comprehensive discussion of the
sionm=2,delayτ =π/2andavaryingthresholdε,insuchawaythat
dataisgiveninNowaczyketal.(2001).Theembeddingwas
theCRPcontainsaconstantrecurrencedensityof20%.
combinedwith the time-delayedmethodaccordingto (Tak-
ens, 1981)in orderto increase furtherthe dimensionof the
phase-spacewiththefollowingrule:Ifwehavenparameters
a ,theembeddingwithdimensionmanddelayτ willresult
i
ina(m·n)-dimensionalphasespace:
Rescaling Function
70 ~x(t) = a1(t),...,an(t),
60 (cid:0)a1(t+τ),...,an(t+τ),
a1(t+2τ),...,an(t+2τ),...
50 a1(t+(m−1)τ),...,an(t+(m−1)τ (13)
φ ( t‘ ) 40 For our investigation we have used a dimension(cid:1)m = 3
=
e t 30 and a delay τ = 1, which finally led to a phase space of
m
Ti
20
250
ment
Fig.31.00T0herescaling2f0unction(blac4Tk0)imdeet e tr‘mined60fromtheCR8P0inFig.2. ARM of Core PS 2180−2 11225050500000 05112 0050 000 ARM of Core PS 2178−3 before adjust
0
It has the expected parabolic shape of the squared coherence in the time 00 220000 440000 660000 880000 11000000 11220000 11440000
Depth in Core PS 2180−2 [cm]
domain.Inredthesquarefunction.
Fig.5. ARMdataoftheboreholesPS2178-3GPCandPS2180-2GPCin
theCentralArcticOceanbeforeadjustment.
105
Depth in Core PS 2178−3 [cm] Depth in Core PS 2178−3 GPC [cm]
0 200 400 600 800 1000 1200 1400 00 200 400 600 800 1000 1200 1400
0
200 200
−2 [cm] 400 GPC [cm] 400
Depth in Core PS 2180 680000 Depth in Core PS 2180−2 680000
1000 1000 50 m]
1200 1200 −−0255205 Deviation [c
Fig.6. Crossrecurrenceplotbasedonsixnormalizedsedimentparameters Fig.7. Depth-depth-curves. InblackthecurvegainedwiththeCRP,inred
andanadditionalembeddingdimensionofm=3(τ =1,ε=0.05). themanuallymatchingresult.Thegreencurveshowsthedeviationbetween
bothresults.
7 Discussion
dimension18(3×6). Therecurrencecriterionwasε=5%
nearestneighbours.
Crossrecurrenceplots(CRP)revealsimilaritiesinthestates
TheresultingCRPshowsaclearLOSandsomeclustering ofthetwosystems. Asimilartrajectoryevolutiongivesadi-
ofblackpatches(Fig.6).Thelatteroccursduetotheplateaus agonalstructureintheCRP.Anadditionaltimedilatationor
inthedata. Thenextstepistofitanon-parametricfunction compressionofoneofthesesimilartrajectoriescausesadis-
(thedepth-depth-curve)totheLOSintheCRP(redcurvein tortionofthisdiagonalstructure(Fig.1). Thiseffectisused
Fig. 6). With this functionwe are able toadjustthe data of tolookintothesynchronizationbetweenbothsystems. Syn-
thePS2178-3coretothescaleofPS2180-2(Fig.8). chronizedsystemshavediagonalstructuresalongandinthe
direction of the main diagonal in the CRP. Interruptions of
The determination of the depth-depth-function with the
thesestructureswithgapsarepossiblebecauseofvariations
conventionalmethodofvisualwiggle matchingis based on
in the amplitudesof both systems. However, a loss of syn-
theinteractiveandparallelsearchingforthesamestructures
chronization is viewable by the distortion of this structures
in the different parameters of both data sets. If the adjust-
alongthe maindiagonal(LOS).Byfittinga non-parametric
ment does not work in a section of the one parameter, one
functiontotheLOSoneallowstore-synchronizationorad-
canuseanotherparameterforthissection,whichallowsthe
justment to both systems at the same time scale. Although
multivariateadjustmentofthe data sets. Therecognitionof
thismethodisbasedonprinciplesfromdeterministicdynam-
thesamestructuresinthedatasetsrequiresadegreeofexpe-
ics, no assumptionsaboutthe underlyingsystems has to be
rience. However,humaneyesareusuallybetterinthevisual
assessmentofcomplexstructuresthanacomputationalalgo-
rithm.
Table1. Correlationcoefficients̺1,2betweenadjusteddataandreference
Ourdepth-depth-curvediffersslightlyfromthecurvewhich dataandtheirχ2deviation. Thecorrelationoftheinteractiveadjusteddata
was gained by the visual wiggle matching (Fig. 7). How- variesmorethantheautomaticadjusteddata. ThedatalengthisN =170
ever, despite our (still) weak algorithm used to fit the non- (wigglematching)andN =250(CRPmatching). Thedifferencebetween
thebothcorrelationcoefficients ̺1 and̺2 issignificantata99%signifi-
parametric adjustment function to the LOS, we obtained a
cancelevel,whenthetestmeasurezˆisgreaterthanz0.01=2.576.
goodresultofadjusteddataseries. Iftheyarewelladjusted,
thecorrelationcoefficientbetweentheparametersofthead- Parameter ̺1,wigglematching ̺2,CRPmatching zˆ
justeddata andthereferencedata shouldnotvaryso much. ARM 0.8667 0.7846 6.032
The correlation coefficients between the reference and ad- MDFARM 0.8566 0.7902 4.791
justeddataseriesisabout0.70–0.80,wherethecorrelation κLF 0.7335 0.7826 2.661
κARM/κLF 0.8141 0.8049 0.614
coefficientsoftheinteractiverescaleddatavariesfrom0.71–
PJA 0.7142 0.6995 0.675
0.87(Tab.1). Theχ2measureofthecorrelationcoefficients INC 0.7627 0.7966 1.990
emphasizesmorevariationforthewigglematchingthanfor χ2 141.4 49.1
theCRPrescaling.
106
madeinorderforthemethodtowork.
250 Thefirstexampleshowstheobviousrelationshipbetween
ARM of Core PS 2178−3after interactive adjustment11205050000 trh(hiFaeeirsgm.L.o2TOn)h.SiecFasfnuiqndnuacalttrlhiyeoe,dnwticinmiactuhreesetadhssoiiasmngpLaaioOnrfasSbthoowefliecfthraLeerqOeucSoeanbnsclsheiyadptoeoerfeirtdnehsettchimsaeeleecCoRtsnhePde-
0 second function to the scale of the first harmonic function
250
51120050000 ARM of Core PS 2178−3after automatic adjustment (ctmhFaeauigisnC.e4Rdd)iPb.a.ygSHotohnomeawleaealdgvnioedfrrf,ietthorheeumnirrcurceesoselnaidcntiieontrhnnoesrhidsaiemptrowptolfiiottehucxdutthresaeoocttfnitmththheeeeLdrdoOeismsSutaloftirrnoatserm.de
0
250 Thesecondexampledealswithrealgeologicaldataandal-
ARM of Core PS 2180−2 11205050000 lojtouhfwsevteswisdauigdcagaoltlmwaeipsmghagoraliwtescoshmneaadwgtciodthhoaitdntahgce.(oFTrneichgsoe.url8vdtiaosanfunctadehlew9coc)i.otmhnTpvthahereneistridoeoefnnepaortelhfnm-tdcheeeetphaatnodhdd--
0 function differs up to 20 centimeters from the depth-depth-
0 200 400 600 800 1000 1200
functionofthewigglematching.Thecorrelationcoefficients
Depth in Core PS 2180−2 [cm]
betweentheCRPadjusteddataandthereferencedatavaries
Fig.8. ARMdataafteradjustmentbywigglematching(top)andbyauto- lessthanthecorrelationcoefficientsofthewigglematching.
maticadjustmentusingtheLOSfromFig.6. Thebottomfigureshowsthe However, the correlation coefficients for the CRP adjusted
referencedata. dataaresmallerthantheseforthewigglematcheddata. Al-
thoughtheircorrelationisbetter,itseemsthattheinteractive
methoddoesnotproducea balancedadjusting,whereasthe
automaticmatchinglooksforabetterbalancedadjusting.
40 Thesebothexamplesexhibitstheabilitytoworkwithsmooth
MDFARM234000 2300MDF_{{ARM} asbmnedhoaonntohdnld-esadmtabo.yoStthhmedaLlalOtaflS,uswcethaueracrtehiobinnygstiahnlegtohrereistunhlomtn.w-Tsimlhleobroeetfhobrdeeat,ttesarmcfoaoonrth-
ingstrategies,likesmoothingorparametricalfitoftheLOS,
400 are notnecessary. Thelatter would dampone advantageof
300
400 this method, that the LOS is yielded as a non-parametrical
κLF300 120000κLF function. A futuretask willbetheoptimizationoftheLOS
200
0 searchingalgorithm,in orderto geta clear LOSevenif the
100
data are non-smooth. Further, the influence of dynamical
30
noise to the result will be studied. Probably, this problem
κ / MLF2300 1200κ / MLF maybebypassedbyasuitableLOSsearchingalgorithmtoo.
κAR10 0κAR
0 Ourmethodhasconspicuoussimilaritieswiththemethod
ofsequenceslottingdescribedbyThompsonandClark(1989).
4 The first step in thismethodis the calculationofa distance
PJA4 02PJA matrix similar to our Eq. 2, which allows the use of multi-
2 variatedatasets. ThompsonandClark(1989)referredtothe
0 distancemeasureasdissimilarity. Itisusedtodeterminethe
100 alignmentfunctioninsuchawaythatthesumofthedissimi-
100 50 laritiesalongapathinthedistancematrixisminimized.This
C 50 0 C
IN 0 −50 IN approachisbasedondynamicprogrammingmethodswhich
−50 −100 weremainlydevelopedforspeechpatternrecognitioninthe
−100
00 220000 440000 660000 880000 11000000 11220000 70’s (e.g. Sakoe and Chiba, 1978). In contrast, RPs were
Depth in Core PS 2180−2 [cm] developedtovisualizethephasespacebehaviourofdynam-
icalsystems. Therefore,athresholdwasintroducedtomake
Fig.9. Theadjustedmarinesedimentparameters. Theconstructionofthe
recurrentstates visible. The involvingof a fixed amountof
CRPwasdonewiththenormalizedparameters. Inthisplotsweshowthe
parameters,whicharenotnormalized. nearestneighboursin the phasespace andthe possibility to
increasetheembeddingdimensionsdistinguishthisapproach
fromthesequenceslottingmethod.
107
8 Conclusion this algorithm.) The nextstep is to set the recurrencepoint
r toanewstartpointandtobeginwiththestepone
iα+1,jβ+1
Thecrossrecurrenceplot(CRP)cancontaininformationabout in order to find the next recurrence point. These steps are
the synchronization of data series. This is revealed by the repeateduntiltheendoftheRPisreached.
distortedmaindiagonal,whichiscalledlineofsynchroniza- Weknowthatthisalgorithmismerelyoneofmanypossi-
tion (LOS).After isolatingthis LOSfromthe CRP, oneob- bilities. Thefollowingcriteriashouldbemetinordertoob-
tainsanon-parametricrescalingfunction.Withthisfunction, tain a goodLOS. The amountof targeted recurrencepoints
one can synchronizethe time series. The underlyingmore- by the LOS N1 should converge to the maximum and the
dimensionalphasespaceallowstoincludemorethanonepa- amountofgapsintheLOSN0 shouldconvergetothemin-
rameterinthissynchronizationmethod,asitusuallyappears imum.AnanalysiswithvariousestimatedLOSconfirmsthis
in geological applications, e.g. core synchronization. The requirement.ThecorrelationbetweentwoLOS-synchronized
comparisonofCRPadjustedgeophysicalcoredatawiththe dataseriesariseswithN1andwith1/N0(thecorrelationco-
conventionallyvisualmatchingshowsanacceptablereliabil- efficientcorrelatesmoststronglywiththeratioN1/N0).
itylevelofthenewmethod,whichcanbefurtherimproved The algorithm for computation of the CRP and recogni-
by a better method for estimating the LOS. The advantage tionof theLOSareavailableasMatlab programmesonthe
is the automatic, objectiveand multivariateadjustment. Fi- WWW:http://www.agnld.uni-potsdam.de/~marwan.
nally,thismethodofCRPscanopenawiderangeofapplica-
tionsasscaleadjustment,phasesynchronizationandpattern
recognition for instance in geology, molecular biology and References
ecology.
Argyris,J.H.,Faust,G.,andHaase,M.,AnExplorationofChaos,North
Holland,Amsterdam,1994.
Acknowledgements. The authors thank Prof. Ju¨rgen Kurths and Dr. Udo
Casdagli,M.C.,Recurrenceplotsrevisited,PhysicaD,108,12–44,1997.
Schwarzforcontinuingsupportanddiscussion. Thisworkwassupported
Eckmann,J.-P.,Kamphorst,S.O.,andRuelle,D.,RecurrencePlotsofDy-
bythespecialresearchprogramme1097oftheGermanScienceFoundation
namicalSystems,EurophysicsLetters,5,973–977,1987.
(DFG).
Kantz, H. and Schreiber, T., Nonlinear Time Series Analysis, University
Press,Cambridge,1997.
Martinson,D.G.,Menke,W.,andStoffa,P.,AnInverseApproachtoSignal
Appendix: AnAlgorithmtoFittheLOS Correlation,JournalofGeophysicalResearchB6,87,4807–4818,1982.
Nowaczyk, N.R.,Frederichs, T.W.,Kassens,H.,Nørgaard-Pedersen, N.,
InordertoimplementarecognitionoftheLOSwehaveused Spielhagen, R.F.,Stein, R.,andWeiel, D.,Sedimentation rates inthe
MakarovBasin,CentralArcticOcean–Apaleo-androckmagneticap-
the following simple two-step algorithm. Denote all recur-
proach,Paleoceanography,2001.
rence points by r (α˜,β˜ = 1,2,...) and the recurrence
iα˜,jβ˜ Prell,W.L.,Imbrie,J.,Martinson,D.G.,Morley,J.J.,Pisias,N.G.,Shack-
points lying on the LOS by riα,jβ (α,β = 1,2,...). Be- leton,N.J.,andStreeter,H.F.,GraphicCorrelationofOxygenIsotope
fore the commonalgorithmstarts, find the recurrencepoint Stratigraphy Application totheLateQuaternary, Paleoceanography, 1,
r next to the axes origin. In the first step, the next re- 137–162,1986.
i1,j1 Romano,M.C.,TheDarkSideofEmbedding,tobepublished.
currencepointr afterapreviousdeterminedrecurrence
iα˜,jβ˜ Sakoe, H. and Chiba, S., Dynamic programming algorithm optimization
pointr istobedetermined.Thisiscarriedoutbyastep-
iα,jβ forspokenwordrecognition, IEEETrans.Acoustics, Speech, andSig-
wiseincreasingofasquared(w×w)sub-matrix,whereinthe nalProc.,26,43–49,1978.
previousrecurrencepointis at the (1,1)-location. The size Takens, F., Detecting Strange Attractors in Turbulence, pp. 366–381,
wofthissub-matrixincreasesstep-wiseuntilitmeetsanew Springer,Berlin,1981.
Thompson,R.andClark,R.M.,Sequenceslottingforstratigraphiccorre-
recurrencepointorthemarginofthe CRP. Whena nextre-
lationbetweencores:theoryandpractice,JournalofPaleolimnology,2,
currencepointriα˜,jβ˜ = riα+δi,jβ+δj (δi = worδj = w)in 173–184,1989.
thex-direction(y-direction)isfound,thesecondstep looks WebberJr.,C.L.andZbilut,J.P.,Dynamicalassessmentofphysiological
if there are following recurrence points in y-direction (x- systemsandstates usingrecurrence plotstrategies, JournalofApplied
direction). If this is true (e.g. there are a cluster of recur- Physiology,76,965–973,1994.
Zbilut,J.P.,Giuliani,A.,andWebberJr.,C.L.,Detectingdeterministicsig-
rence points) increase further the sub-matrix in y-direction
nalsinexceptionallynoisyenvironmentsusingcross-recurrencequantifi-
(x-direction) until a predefined size (w +dx˜)×(w +dy˜)
cation,PhysicsLettersA,246,122–128,1998.
(dx˜ < dx,dy˜ < dy) or until no new recurrence points are
met. This further increasing of the sub-matrix is done for
the both x- and y-direction. Using dx˜,dy˜we compute the
nextrecurrencepointr bydeterminationofthecen-
iα+1,jβ+1
ter of mass of the cluster of recurrencepoints with iα+1 =
iα+(dx˜+δi)/2andjβ+1 = jβ +(dy˜+δj)/2. Thelatter
avoidsthefactthatthealgorithmisdrivenaroundwidespread
areasofrecurrencepoints. Insteadofthis,thealgorithmlo-
cates the LOS within these areas. (However, the introduc-
ing of two additional parameter dx and dy is a disadvan-
tage which should be avoided in further improvements of