#### Abstract

A simple and effective algorithm for the identification of optimal time delays based on the geometrical properties of the embedded attractor is presented in this paper. A time series synchronization measure based on optimal time delays is derived. The approach is based on the comparison of optimal time delay sequences that are computed for segments of the considered time series. The proposed technique is validated using coupled chaotic Rössler systems.

#### 1. Introduction

The reconstruction of geometry from a time series is one of the paradigmatic algorithms used for the computational analysis of nonlinear dynamical systems. There is a number of algorithms for the analysis of time series that consider phase spaces, that is, the set of all possible system states, of dynamical systems. For a deterministic dynamical system, the phase space is known from the mathematical model (equations of motion). Thus, with the knowledge of the system state at a particular time, it is always possible to determine the future states. However, real-world or experimental dynamical systems are usually too complex, and phase spaces of such systems are unknown due to the chaotic nature of the analyzed data. Therefore, phase space reconstruction methods had to be developed. A branch of applied dynamical system theory, analyzing such algorithms, that is, extracting the information about geometrical and topological properties of the phase space of time series obtained by performing measurements on an evolving system, is called embedology [1]. One of the first methods for reconstructing the phase space system from a time series by using a time delay variable was proposed several decades ago [2]. A complex time series has an inherent geometry, and Packard et al. [2] were first to show that a representative geometry of a dynamical system can be obtained by using a time series of one of its observables.

The much celebrated Takens embedding theorem proved that it is possible to reconstruct the attractor of a dynamical system using only a single time sequence of scalar measurements [3]. The vectors accomplishing this reconstruction have the form , where are the observations on the system, the integer is the embedding dimension, and is the time difference between consecutive components (also called the time lag or the delay time) and is a multiple of the sampling time step. Takens assumed an infinite sequence of noise-free measurements and proved the existence of a diffeomorphism between the original and the reconstructed attractors for almost any choice of positive delay times and a sufficiently big dimension .

Although the application of the Takens time delay technique is straightforward for almost any time series, the selection of optimal reconstruction parameters and is a nontrivial problem due to their dependence on the nature of the analyzed data. Thus, numerous studies have been made in this area. The selection of optimal time delay value is usually based on the optimization of a particular target function that corresponds to some measures of the quality of phase space reconstruction. Classical statistical approaches for the construction of the target function include the autocorrelation function method and the mutual information method [4–6]. Such techniques are based on maximizing the independence between the coordinates of the delay vector. This approach is heavily reliant on the geometry of the phase space reconstruction, which is appropriate for the study of parameters such as fractal dimension and Lyapunov exponents. However, the selection of such target function does not always yield best results for nonlinear time series prediction problems [7].

Another class of methods for the selection of an optimal is based on the concept of phase space expansion. The “fill factor” as a spatial measure of the phase space was first introduced by Buzug and Pfister in [8, 9]. The fill factor quantifies an attractor’s utilization of the embedding space as a function of . While this method does not possess the drawbacks of statistical techniques with respect to nonlinear data, it is computationally intensive (for ) and is prone to overfolding of the attractor [10]. A more computationally effective algorithm based on the expansion of the area for all pairwise planar projections of the embedded attractors is presented in [11]. Similar approaches based on the spatial properties of phase spaces include the average displacement method [10], the SVF method [12], and the wavering product [13] method.

A different viewpoint considers both phase space reconstruction parameters simultaneously since some experiments indicate that an irrelevant relationship between and may impact the congruence between the original system and the reconstructed phase space. This class of methods includes the small-window solution [14] and the C-C method [15].

Novel approaches to the selection of embedding parameters have been developed in recent years. It is shown in [16] that an improved phase space reconstruction method based on manifold embedding and Laplacian eigenmaps can reveal the structure of the chaotic attractor. A new algorithm for the estimation of the dimension of chaotic dynamical systems using neural networks and robust location estimate is proposed in [17]. Statistical analysis of the nearest neighbors is exploited for the reliable estimation of minimum embedding dimension for noisy time series in [18].

It has been demonstrated in [11, 19–21] that nonuniform embedding (when all time delays are not equal) performs better than uniform embedding in a variety of applications—typical examples are causality and coupling detection and time series prediction. However, the procedure for the selection of time delays is usually related to a particular application and is implemented by introducing a specific target function which determines the utility generated by a concrete nonuniform embedding.

This paper has two main objectives. The first objective is to introduce a simple and effective algorithm for the identification of optimal time delays which is based on the geometrical properties of the embedded attractor and is applicable for short time series.

The second objective is to introduce a time series synchronization measure based on optimal time delays. A short review of commonly applied methods for the detection of synchronization is given below.

The first linear synchronization measures based on correlation analysis, namely, cross correlation and coherence functions, are widely applied due to their computational effectiveness [22]. However, such measures can only detect the most straightforward regimes of synchronization.

Nonlinear synchronization measures were introduced in order to quantify more complex synchronization effects. The mutual information measure is based on Shannon entropy and takes into account not only linear but also nonlinear dependencies [23]. Phase synchronization measure is used to quantify similarity between cyclic signals and time series. Two approaches—based on the Hilbert or wavelet transform—can be used to implement the phase synchronization detection algorithm [22].

A new form of synchronization between coupled chaotic oscillators called amplitude envelope synchronization has been discovered in [24]. Generalized synchronization of dynamical systems occurs if dynamical variables from one subsystem are a function of the variable of another subsystem [25]. A nonlinear interdependence of dynamical systems based on state space reconstruction is a similar approach to generalized synchronization but does not require a functional relationship between the dynamics of the underlying systems [26]. The function projective synchronization in relay coupled systems is studied in [27]. A novel sort of synchronization called complex antilag synchronization is introduced in [28]. A new approach for the investigation of hybrid chaos synchronization in discrete-time hyperchaotic dynamical systems based on stability theory of linear discrete-time systems and Lyapunov stability theory was proposed in [29].

This paper presents a novel time series synchronization measure based on the geometrical approach towards optimal time delays. This approach is based on the determination of time lags that maximize the volume of the state space occupied by the embedded segments of the time series. The sequences of obtained time lags reveal the level of synchronization between 2 time series. The presented algorithm for the identification of synchronization between two time series is validated using coupled chaotic Rössler systems.

#### 2. Preliminaries

##### 2.1. Two-Dimensional Delay Coordinate Space

Let us consider a harmonic function , where is the angular frequency and is the phase of harmonic oscillations, and the amplitude is set to 1. A pair of function values and , where is the time lag, is mapped into a point in the phase plane , where and are the coordinates of the delay coordinate plane. Two-dimensional time delay embedding maps the harmonic function into an ellipse; the equation of that ellipse reads

The geometrical shape of the ellipse can be exploited for the quantification of the quality of the embedding. Such geometric approach (for a two-dimensional phase plane) was firstly proposed by Buzug and Pfister in 1992 [9]. The radiuses of the embedded ellipse and can be directly expressed from (1):

Thus, the area of the ellipse reads

The maximum area of the ellipse is (when the ellipse becomes a circle). Now, the function representing the quality of embedding into a two-dimensional delay coordinate space can be defined as a ratio between the area of the ellipse and the area of the circle:

##### 2.2. -Dimensional Delay Coordinate Space

###### 2.2.1. The Embedding Quality Function

The geometrical approach for embedding a harmonic function into a two-dimensional phase plane is generalized for the -dimensional delay coordinate space in [11]. The coordinates of the reconstructed point in the -dimensional delay coordinate space read where is the time delay vector; is the total embedding window. It can be observed that there exist planar projections of the embedded ellipse. The function representing the quality of embedding into the -dimensional delay coordinate space is constructed as an arithmetic average of all quality functions is all possible planar projections [30]:

In the case of the uniform embedding (when ), (6) reads

The equality implies that the harmonic function is compressed into line segments in all possible planar projections. The maximization of in respect of time lags results in a maximum average area of ellipses in all possible planar projections.

###### 2.2.2. A Nonharmonic Function

Let us consider a function defined in the time interval , where and are the limits of the observation window. Every harmonic component of is affected by the quality function when the appropriate harmonic signal is embedded into the -dimensional delay coordinate space. Harmonic components with frequencies where is small will be suppressed (in average) in all possible planar projections more than those harmonic components where is large. Such motivation did help to construct the optimization problem for the identification of an optimal set of time delays which do result in the richest representation of the attractor in the delay coordinate space: where

is the optimal set of time delays and is the Fourier amplitude spectrum of the signal in the observation window [11]. The integral in the denominator is used in order to normalize the target function in respect of the signal (this integral can be computed once at the beginning of the optimization process). The multiplier is used to normalize the target function in respect of the white noise (the value of the target function for the white noise signal now becomes equal to 1 for any time delays not equal to 0) [11].

#### 3. The Proposed Algorithm

The technique proposed in [11] enables a fast and effective determination of the optimal set of time delays by assessing the geometrical properties of the embedded attractor. However, the optimization problem in [11] does not asses the phase spectrum of the Fourier transform. In other words, [11] is an approximate algorithm for the determination of optimal time delays when every discrete harmonic component of the Fourier amplitude spectrum is treated as a separate individual harmonic component. Nevertheless, such an approximating approach appears to be useful in practical applications—especially in time series forecasting algorithms based on fuzzy neural networks [31–33]. However, it is clear that the Fourier phase spectrum does have an impact on the geometrical representation of the reconstructed attractor in the phase space. In other words, a more accurate assessment of the geometrical properties of the reconstructed attractor is required. Let us consider a continuous function in the observation window and the -dimensional delay coordinate space with the set of time lags . Then, instead of computing any projections of the embedded attractor into a 2-dimensional phase planes, we compute the average distance of the attractor points to the origin of the embedding frame:

##### 3.1. Properties of : Embedding a Harmonic Function

Let us consider that the embedding is uniform and the embedded function is a harmonic function , where , , and are the amplitude, the circular frequency, and the phase, accordingly. Then,

The change of variables yields

Now,

Let us denote . Then,

Note that when . Thus, is the maximal value for the harmonic function. On the opposite, the pole is a removable singularity at and is the minimal value for the harmonic function.

##### 3.2. The Comparison between and

We continue with the harmonic function and uniform embedding. Then, and (9) yields , where is the Dirac delta function. Without the loss of generality, we assume that and . Then, ; . Graphs of and are illustrated in Figure 1.

**(a)**

**(b)**

**(c)**

**(d)**

The first observation is that graphs of and are similar at (Figure 1(a)). The best embedding of the harmonic signal into the 2-dimensional delay coordinate space is produced at , (then the harmonic function is embedded into a circle). The harmonic function is embedded into a line interval on one of the diagonals of the frame at , . The difference is that now , but (note that the scale of the *y*-axis for is shown on the left, but the scale for is shown on the right side of Figure 1(a)).

The differences between and are more clear at (Figure 1(b)). The harmonic function cannot be mapped into circles in all 3 projections at a fixed ; therefore, the maximum value of at is lower than that at (note that the scales of the axes are the same for different in Figure 1). But remains the same and does not depend on the dimension of the delay coordinate space. Both and yield the best embedding at , where , .

The maximum value of at becomes smaller than the one at (Figure 1(c)). But now a major difference between and can be observed in terms of the best embedding. yields the best embedding at , but yields the best embedding at , where , . Clearly, misses the time delay . That can be explained by the fact that measures the average area of the mapped attractor in all possible plane projections of the -dimensional delay coordinate space. On the contrary, directly measures the geometrical volume occupied by the attractor without considering any projections of that attractor.

The maximum value of tends to 1 when (Figure 1(d)). Note that when the embedded time series is the white noise [11]. However, as mentioned previously, is an invariant value in respect of . What is more important, yields the best embedding at , where , , which is in stark contrast to (Figure 1(d)).

##### 3.3. Embedding a Discrete Time Series

Let us consider a finite scalar discrete time series , where is the number of observation points uniformly spread in the time axis. Let the sampling rate of the time series (the time interval between adjacent observations) be . Then, all time delays are multiple values of . Let us denote . Then, (10) reads

Note that does not play any role on the right-hand side of (17)—because the definite integral in (10) is replaced by the limit sum in (17). Now, the problem of the identification of the optimal set of time delays is formulated as follows: where is the upper limit for every individual time delay. The total embedding window (the time delay between the first and the last time delay coordinates) is , and it could be also used for the estimation of . It is clear that must be greater than . However, must be at least 10 times greater than in order to produce a meaningful optimization result (18).

##### 3.4. Properties of : Embedding the Random Noise

It is well known that embedding of the random noise is invariant to time delays. In other words, the geometric shape of the reconstructed attractor does not depend on time delays (if only the embedded time series is an ergodic random process and the realization of that process is long enough). This geometric feature is exploited in one of the tests used to measure the quality of random number generators. For example, if a random number generator yields a sequence of numbers distributed uniformly in the interval , then the density of the embedded points should be the same in any location of a two-dimensional square for a sufficiently long sequence and any time delay greater than 1 [34]. It is clear that must be invariant to time delays if only the embedded time series is random. We demonstrate the properties of for the random time series distributed uniformly in the interval and for the Gaussian time series. Of course, the invariance in respect to time delays does hold, but some features of are interesting and are presented below.

###### 3.4.1. Embedding the Uniform Noise

Let us consider a discrete random variable distributed uniformly in the interval . Initially, let us consider and . In fact, this situation cannot be considered an embedding because all values of the time series are mapped on a single one-dimensional axis without any time delays. However, the formal result reads

Now, let us consider the situation when (for any ):

Thus, finally, if the embedded random time series is distributed uniformly in the interval .

###### 3.4.2. Embedding the Gaussian Noise

Now, let us assume that a discrete random variable is the Gaussian noise, with the mean value and standard deviation . Firstly, let us consider and :

Now, let us consider the case when (for any ):

Thus, finally, if the embedded random time series is the Gaussian noise, with the mean and the standard deviation .

##### 3.5. Computational Experiments: Optimal Embedding of Mackey-Glass Time Series

Equation (18) is a nonlinear integer optimization problem. The full sort method could be used for the identification of the optimal set of time delays—only if , , and are not very large and computational resources are not strictly limited. Otherwise, soft computing techniques could be exploited for the determination of the near-optimal set time delays [30]. Let us consider Mackey-Glass time series [35] the discrete numerical solution to chaotic Mackey-Glass delay differential equation that reads where , , and ; the observation window is and the time step (Figure 2). The standard FNN algorithm [36] yields (5) as the minimum embedding dimension for the investigated chaotic Mackey-Glass series.

Now, we use (18) to determine the optimal uniform time lag in the interval (Figure 3). The maximum value of is reached at (Figure 3). Note that this value of can be used as a lower bound for the integer programming problem when the near-optimal set of nonuniform time delays is sought using branch-and-bound [37] or other soft optimization algorithms.

The full sort method is used to determine the optimal time delay vector for nonuniform embedding. 100^{4} different combinations are assessed in 8.14 hours on a two-kernel 2.20 GHz processor with 8 GB RAM. As mentioned previously, soft computing techniques can be effectively exploited for the determination of the near-optimal set of time delays. Finally, the optimal set of time delays reads ; . All possible planar projections of the embedded Mackey-Glass attractor are illustrated in Figure 4.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

##### 3.6. Synchronization Measure Based on a Geometric Approach to Attractor Embedding

The proposed technique for the determination of optimal time delays is based on a straightforward geometric approach. Such an approach can be effectively exploited for comparing geometric properties of attractors embedded from different time series. Moreover, we will demonstrate that such an approach is working well with short time series.

The basic idea of the proposed synchronization feature is based neither on the direct comparison between geometric shapes of the reconstructed attractors nor on global characteristics of those attractors (like averaged maximum Lyapunov exponent). We will use the optimal time delay as a single numerical characteristic describing the geometry of the reconstructed attractor in a 2-dimensional delay coordinate space. The reasoning for such a measure is straightforward. It follows from the properties of the reconstructed attractor (Figure 4). Such synchronization measure is much more general compared to synchronization measures used to detect identical, phase, amplitude envelope, lag, or even generalized synchronization. Two attractors are considered to be similar if their optimal time lags are the same. The higher the difference between optimal time lags, the more different these attractors are considered (from the geometrical point of view).

We will demonstrate that such an approach for detecting synchronization between two different time series is highly feasible—especially if those time series are far from being similar.

Let us consider two paradigmatic chaotic Rössler systems with diffusive coupling: where , , , and ; is the coupling parameter. The initial condition is set to , , and .

Attractors of both Rössler systems are illustrated in Figure 5. The difference between the uncoupled attractors can be seen by a naked eye (Figures 5(a) and 5(b)). However, the attractors become similar as the coupling increases (Figures 5(e) and 5(f)).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Analogously, attractors of coupled Rössler systems (26) are visualized in the embedding space (Figure 6). Note that a short observation window (400 scalar data points) is used to reconstruct the attractor in the planar delay-coordinate space. One time delay is used for every embedding; optimal values of time delays are depicted in Figure 6. Note that both the geometric shape of the embedded attractors and the optimal time delays become closer as the coupling parameter grows (Figure 6).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

It is also interesting to compare the maximal Lyapunov exponents calculated for the system (26) and the one in the embedding space. Averaged maximal Lyapunov exponents for the uncoupled Rössler systems are computed separately: for the first system and for the second system. However, the coupled Rössler systems (26) are considered a single system. The 6-dimensional Jacobian matrix of the coupled system yields (after averaging) a single maximal Lyapunov exponent: at and at . It can be observed that the single maximal Lyapunov exponent of the coupled system at approaches the average of and . But the maximal Lyapunov exponent decreases (the coupled system becomes less chaotic) at .

The maximal averaged Lyapunov exponents computed directly from time series and in the embedded space are 0.0485 and 0.0541 for , 0.0367 and 0.0461 for , 0.0326 and 0.0344 for . Note that the maximal Lyapunov exponents computed directly from time series and are underestimated if compared to the maximal exponents for the coupled systems in (26). This effect is well known for synchronized differential equations [38]. Nevertheless, both methods demonstrate that the coupled systems become less chaotic when the coupling parameter is increased.

Figure 7(a) shows the evolution of and and the difference between and when two Rössler systems are uncoupled (the length of the observation window is 8000 points). Figure 7(b) shows , , and when , and Figure 7(c) corresponds to . Note that the scales of the vertical axes in the graphs in Figure 7(c) are different. Time series and become more similar as the coupling parameter increases; however, the difference remains chaotic even at (Figure 7). The whole observation window is now split into 20 equal intervals; each interval contains 400 points (both for and for ). The optimal time lag (2-dimensional delay coordinate space) is identified for time series in each individual interval according to (18); . Optimal time delays for and are represented by a solid and a dashed line correspondingly in Figures 8(a)–8(c); the correlation between two sequences of optimal time delays is , , at , and , respectively. It appears that the proposed algorithm based on the geometric approach to optimal attractor embedding is able to detect the synchronization between two chaotic series in an effective and efficient way. Note that the difference between and is not equal to zero at (Figure 7). However, it appears that (the correlation between two sequences of optimal time delays) is a good indicator of the geometrical synchronization between two sequences. It represents a geometric measure of synchronization between different time series. But instead of comparing geometric shapes of the embedded attractors, it compares sequences of optimal time lags in a predefined number of observation windows.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

Computational experiments are continued with the coupled Rössler systems with the additive Gaussian random noise. Note that different realizations of the Gaussian random noise are added to and . Moreover, the noise is added after the integration of (26) and thus is not smoothed by the numerical integrator. Computational experiments are repeated for different values of the coupling coefficient (, , and ) and with the different parameters of the additive Gaussian random noise. The Gaussian noise with the zero mean and the standard deviation equal to 2 is used in Figure 9, and that with the zero mean and the standard deviation equal to 18 is used in Figure 10.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

Again, the whole observation window is split into 20 equal intervals (each interval contains 400 points both for and for ). Optimal time lags for the coupled Rössler systems with the additive noise are represented by a solid and a dashed line correspondingly in Figures 8(d)–8(f) (for the standard deviation equal to 2) and in Figures 8(g)–8(i) (for the standard deviation equal to 18).

The correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 2) reads −0.39219, 0.36979, 0.79285 at , , and , respectively (Figure 8). The differences between and are large, random, and unpredictable (Figure 9). Nevertheless, the proposed algorithm is capable of identifying the synchronization between two sequences even in such a complicated computational experiment.

Analogously, the correlation between sequences of optimal time delays for the coupled Rössler systems with the additive noise (with standard deviation equal to 18) reads 0.094789, 0.23079, and 0.25094 at , , and , respectively (Figure 8). Now, the proposed algorithm for the identification of the geometrical synchronization between and reaches the limits of applicability. Simply, the signal-noise ratio is so low that it is impossible to trace the presence of any synchronization between two time series.

##### 3.7. Computational Experiments: Comparisons and Extensions

This paper presents a new technique for the quantification of geometric synchronization between different time series. This technique is based on the identification of optimal time lags which maximize the state space occupied by the embedded segments of scalar time series. A natural question is whether other existing methods for the determination of optimal time lags could be used instead.

Computational experiments are repeated with the coupled Rössler systems, but the autocorrelation method for the determination of optimal time lags [5, 6] is used instead of the proposed method. All parameters of the computational setup are kept the same—including the coupling parameter, system parameters, and the length of the observation window (Figure 11).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

The first observation is that the autocorrelation method yields incorrect correlation for the coupled Rössler systems at ; it appears that the correlation between uncoupled systems (Figures 11(a) and 11(b)) is higher than that between coupled systems (Figures 11(d) and 11(e)). Secondly, the autocorrelation method yields wrong results for the coupled Rössler systems with the additive noise—the correlation between the coupled systems with the additive Gaussian noise is higher than that between the coupled systems without the additive noise (Figures 11(b)–11(h)).

Computational experiments are now repeated with the mutual information method for the determination of optimal time lags [4] (Figure 12). Again, correlations between the coupled systems with the additive noise (Figures 12(b)–12(e)) are higher those between the coupled systems without noise (Figures 12(c)–12(f)). Also, correlations between strongly coupled systems (Figures 12(e) and 12(f)) are lower than correlations between weakly coupled systems (Figures 12(h) and 12(i)).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

These computational experiments demonstrate that the identification of optimal time lags based on the maximization of the state space occupied by the embedded segments of scalar time series is a perfect approach for the presented synchronization measure.

Computational experiments are continued with a completely different system. The proposed synchronization measure is tested using a coupled system of periodically forced nonlinear pendulums: where and are angular coordinates of the coupled systems, is the damping coefficient, and are the external forcing amplitudes of both systems, and is the frequency. A driven damped nonlinear pendulum is a paradigmatic model in nonlinear dynamics; the model yields chaotic behavior at , , and [39]. The parameters of the coupled model are set as follows: , , , and .

The proposed synchronization measure for (27) is presented in Figure 13. At first, the coupled system of nonlinear pendulums is considered without the additive noise (Figures 13(a)–13(c)). The synchronization measure for the uncoupled system () is . The coupling yields at and (complete synchronization) at .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

Computational experiments are continued by adding the additive Gaussian noise. The proposed method is still able to detect synchronization between the systems when the standard deviation of the additive Gaussian noise is equal to 0.6 (Figures 13(d)–13(f)). Finally, the method is not able to detect synchronization when the standard deviation of the additive noise is set to 5.6 (Figures 13(h)–13(i)).

The ability of the proposed technique to identify the coupling between two chaotic systems (even in the presence of a moderate additive noise) could be explained by the nature of the algorithm which is used to identify optimal time lags. We directly assess the geometric properties of the embedded attractor in the delay coordinate space. The proposed technique is indeed applicable to short time series; the length of the observation window must be sufficient to reconstruct at least one full orbit in the reconstructed delay coordinate space (Figure 6). This is in contrast with the autocorrelation method or the mutual information method requiring much longer observation windows [4–6].

#### 4. Concluding Remarks

This paper presents a novel approach towards geometric synchronization of chaotic time series. This approach is based on the correlation between sets of time lags used for the reconstruction of attractors from different time series in delay coordinate spaces. An efficient algorithm is proposed for the identification of optimal time delays in finite observation windows. Computational experiments with the Mackey-Glass time series, coupled Rössler systems, and coupled chaotic pendulums are used to demonstrate the efficacy of the presented techniques.

It is natural to expect that every method has its strengths and weaknesses. The No-Free-Lunch theorems [40] yield a conclusion that it is always possible to find a problem which will defeat any algorithm. The proposed method does not compare the shapes of the embedded attractors; instead, it compares optimal time delays resulting in a maximal expansion of the embedded attractors in the delay coordinate space. Therefore, it should be possible to design such two artificial time series which are not synchronized at all, but the presented technique would yield a positive synchronization measure (if only optimal time lags would coincide with appropriate observation windows). It is completely natural that every algorithm does have its bounds of applicability. In the present case, one should be cautious in testing the synchronization between time series of a completely different origin.

The proposed synchronization measure is based on a geometric approach to attractor embedding using finite observation windows. However, a short observation window is not the main strength of the presented method. The proposed technique uses a direct geometrical estimate of the embedded attractor in the delay coordinate space. In other words, the presented method compares the geometric features of the embedded attractors (in terms of the optimal delays). Such generalized geometric approach opens new possibilities for the design of new geometric synchronization measures between different time series.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by Jiangsu Provincial Recruitment Program of Foreign Experts (Type B, Grant no. JSB2017007).