Parameter estimation for 1D PWL chaotic maps using noisy dynamics

Many physical situations involve chaotic systems implemented in hardware. Among them one-dimensional piecewise linear maps are popular candidates for such applications because of their property of generating robust chaos. In physical implementations, the control parameter of these maps may deviate from its ideal value due to hardware imprecision. Since the dynamics of a chaotic map is completely defined by its control parameter, one needs to know the value of the parameter in a hardware realisation. In this paper, we show that it is possible to determine the parameter, through the realisation of the unstable fixed point of the map, by utilising noise that is always present in the system. We present this in the form of an algorithm and demonstrate its efficacy through simulated results. We also determine the bounds on the signal-to-noise ratio required for successful parameter estimation. The proposed approach is expected to be beneficial to the existing noise reduction techniques and time series recovery algorithms that require a reasonably accurate knowledge of the map.


Introduction
Chaotic maps have found use in applications like data encryption, cipher key generation in communication, dynamical systems synchronisation and control, and signal processing, to name a few. One-dimensional (1D) piecewise linear (PWL) maps, a sub-class of chaotic functions, are especially popular in practical applications owing to their simplicity and ease of implementation in the physical domain. In particular, certain maps belonging to this class are extensively studied and applied for their uniformly dense chaotic behaviour over a wide range of parametric values, referred to as 'robust chaos' [1]. Such maps have therefore received considerable attention for use in applications involving physical implementations [2][3][4][5][6].
The control parameter is one of the central elements of the chaotic functions, whose estimation is necessary for the specific knowledge of the resultant dynamical behaviour. There are some robust techniques available that use either time series trajectories [7,8] or the equivalent symbolic sequences [9] to decipher the dynamic behaviour of the map, thereby determining the control parameter that result in such behaviour. However, such techniques often do not consider the inevitability of the noise introduced in certain implemented scenarios that results into complete digression of system trajectories from the expected one [10]. In fact, in some cases, the severity of the noise can lead to outcomes bearing little or no resemblance to the actual system behaviour [11,12].
Distinguishing meaningful information from the distorted dynamics is crucial and can be computationally intensive. The contributions in [13][14][15] offer significant reduction of the noise with efficient reconstruction of the dynamical phase space manifolds, consistently highlighting the importance of either determining the noise level in the system or approximating the source function as a priori information so that either of the two information can be utilised to distinguish the other. In this paper, we study the dynamics of noise-affected 1D PWL maps and propose an innovative technique utilising Cartesian coordinate geometry to identify the control parameter of the map, consequently this knowledge of the control parameter can significantly reduce the computational steps involved in approximation of the source function.
Many studies have been conducted to broadly understand the effect of noise on the chaotic dynamics [16][17][18]. Noise and chaotic trajectories are both wideband signals; therefore, linear filtering techniques (e.g. lowpass, bandpass filters) or classical Fourier approaches cannot be applied to reduce noise, as meaningful trajectories might also be barred by the filtration [19,20]. Thus, eliminating noise in order to discern the actual system trajectories is a non-trivial problem and is approached algorithmically.
Works by Kostelich [21], Schreiber [13] and Grassberger et al. [15] formed the basis of the subsequent developments in the noise treatment approaches and have been well summarised by Kostelich and Schreiber in [22] through a qualitative comparison assessing the suitability of the most prevalent techniques for various application scenarios. One such technique is source function approximation which uses least squares approach for local linear approximations (Eckmann-Ruelle linearisation) [23] of the dynamics at each point of the attractor. The accuracy of such approaches depends on how appropriately the dynamics can be described by the chosen linear function. Similar approximation techniques like global function fits also depend on the appropriate choice of the basis functions [22,24].
Adaptive thresholding is also a commonly used approach that uses threshold functions to determine wavelet coefficients of the noisy signal, depending on the energy distributions in each set of samples, since energies of chaotic segments are comparable to noise. Such a method has been described in [25] that applies different threshold coefficients adaptively according to the detailing of the decomposed signal-scale coefficients. Adaptive thresholds can also be estimated with optimisation techniques like genetic algorithm as shown by Han and Chang [26].
Another popular approach is to determine maximum likelihood, where the noisy trajectories are compared with a noise free reference orbit, as demonstrated by Marteau and Abarbanel in [27]. Similarly works by Schweizer and Schimming [28] shows how least square estimates between the expectation and the observation can be used as likelihood cost function. However, to determine the expectation, the knowledge of either the actual trajectory or the noise level in the system is necessary. Also, most of the noise reduction and phase space reconstruction problems are simultaneously addressed through multidimensional delay embedding techniques that are fundamentally based on Takens' embedding theorem [19], which is only applicable to the systems with higher (greater than two) Euclidean dimensions [11,29].
From the mentioned techniques, it is understood that a number of the approaches depend either on good approximation of the source function, or determination of the thresholds between noisy and deterministic data that involve several layers of pre-processing stages, thus may add to the cost of computation as a trade-off. On the other hand, the likelihood estimation approaches demand prior information of the noise level or complete knowledge of the source function where by "complete knowledge" it is expected of one to know the chaotic system in full detail including the control parameter. Hence, if the information of the control parameter can be extracted from the noisy trajectories, then a considerable amount of information regarding the actual dynamics and the initial conditions can be successfully retrieved and it can be of additional support to the existing noise reduction techniques to improve accuracy.
In this paper, we investigate the behaviour of the noisy dynamics of 1D PWL maps under different parametric conditions. From the collective observations of highly distributed trajectories over the state space, as is the case with the noise-affected dynamics, we notice a unique dynamical property demonstrated by the iterates originating within a wide neighbourhood of the fixed point of the map. Knowing the magnitude and time coordinates of the iterative samples, we propose a geometric view of the iterates on a two-dimensional Cartesian coordinate system. We apply linear constructions between consecutive iterates within the set of collected samples whose solutions correspond to a unique point that we mathematically establish as having direct correspondence with the fixed point and control parameter of the map. Hence, the property can be utilised to determine the system parameter from the collection of noisy dynamics of the system.
In Sect. 2, we present the contextual overview of 1D PWL maps and suitable sampling techniques that will be required for the observations regarding the noisy dynamics. In Sect. 3, we show in detail how the dynamical noise within the system affects the trajectories and makes them spread all over the state space through the dynamics of tent map considering it as a typical example of 1D PWL maps. We establish the observed properties in Sect. 4, which we further utilise in Sect. 5 to formulate a method to estimate the parameter of the maps considered. We verify the efficacy of the proposed technique for the candidates of 1D PWL maps through simulation results. In Sect. 6, we add some final comments and concluding remarks about the proposed technique.

Context
1D PWL maps given by f are composed of a set of linear functions defined over segments of a onedimensional state space. Each of the linear segments is limited by contiguous restrictions that together constitute the invariant interval or the state space I . The dynamics of f (x n ) = x n+1 ∈ I ⊂ R is an iterative self-map such that f : I → I , whose resulting discrete time series is given by X = {x n | n = 0, 1, . . . , N − 1} for N iterations and can be treated as the trajectory of the map over time.
Among the 1D PWL maps, the tent map (both symmetric and skew), the zigzag map and Bernoulli shift map are the are the common examples, as these maps are widely used in various areas of application. First, we consider symmetric tent map T (x n ) [3] to observe the behaviour of noise-affected chaotic dynamics and investigate the properties of the noisy iterates reflecting upon the maps fixed points and the parameter. Later we adapt the estimation technique for the other 1D PWL maps. The tent map is defined by (1) over I = [0, 1] and the dynamics of T depends on the height or the control parameter of the map, given by μ ∈ [0, 1], and critical point of T is x c = 0.5 ∈ I .
T (x n ) exhibits chaotic behaviour for μ ∈ (0.5, 1]. The discrete time trajectory of an input or initial condition x 0 ∈ I has the following properties. 1.
There are two fixed points in the state space I such that x n+1 = x n . One such fixed point is evidently The other is the nonzero fixed point is given by The folding nature of the tent map ensures that every point in the invariant interval I maps arbitrarily close to every other point in I [30]. When chaotic maps are implemented physically, the feedback process of the original dynamics of X trajectories becomes significantly affected by the noise induced through the system hardware. Noise-affected chaotic trajectories consist of a deterministic part and a random part, intrinsic to the physical system that affects the entire dynamics at every iterative state. Such types of noise can be described as dynamical noise [10], whose evolution is given by ï n+1 = T (ï n ) + ó n , where ó n is the uncorrelated random perturbations at every step of iteration. The trajectory of such noisy dynamics is defined as The distribution of the random variables ó n from an unknown source can be best emulated by Gaussian distribution, commonly referred to as additive white Gaussian noise (AWGN) due to its intrinsically additive nature. The AWGN with a zero mean is characterised by signal-to-noise ratio (SNR) of 10 log 10 (σ 2 x /σ 2 ó ) measured in dB, with σ x and σ ó as the standard deviations of the signal and the noise, respectively. To investigate the behaviour of the noisy dynamics, we collect M samples of η trajectories. Therefore, each of the sampled trajectories (as shown in Fig. 1) can be For any m, we record each nth iterate through all the N iterations and then for the next (m + 1)th sample, we start from ï m+1 0 again through to ï m+1 N −1 . As summarised in [22], a noise reduction method must estimate the dynamical quantities of the function under operation. Even when the chaotic function is known, the true dynamics of the system can only be ascertained when some knowledge about the control parameter is also available.

Behaviour of the chaotic trajectories when affected by noise
For any perturbed point x n + n iterative transformation x n+1 + n+1 = f (x n + n ) leads to diverging paths, that is further deviated from the actual x n+1 = f (x n ), at the rate given by the Lyapunov exponent λ = ln( n+1 / n ). For the tent map, as long as the control parameter belongs within its ergodic range, i.e. μ ∈ (0.5, 1] and n+1 > n i.e. λ is positive, the trajectories will be divergent in nature, similar to what is also experienced when random perturbations ó n is introduced in every iteration. Figure 2 shows this divergence within a short length (N = 5) trajectory of an arbitrary initial condition for M = 10 samples. The η m sample trajectories diverge away from each other and the ï m n points become eventually distributed over the state space I = [0, 1]. As can be seen from Fig. 3, when larger set of M trajectories are observed collectively for a higher N , the η m are found to be densely distributed over I . We notice that when consecutive iterates ï m n−1 and ï m n are connected by straight lines, dense clusters of intersections appear in a con- We find that such clusters of intersections (or crossovers) majorly appear at about the same level on the Y -axis, conveniently around the nonzero fixed point x f . We also verified that any variation in μ is reflected on the positions of the clusters. Figure 4 shows that the positions of the intersection clusters have appeared at different levels for the two cases of μ, simulated for the same x 0 perturbed with same level of noise (SNR = 25 dB). Based on these observations, we study the intersections of the straight lines that are formed between the iterates and investigate their appearance

Dynamics in the neighbourhood of the nonzero fixed point
In order to realise the properties of the intersection clusters (crossovers), we explore the dynamics of the map around x f . For any parameter μ ∈ (0.5, 1], there exists a preimage of will not exceed x f , thus resulting into the mapping On the contrary, for any x n ∈ [x p , 1] ⊂ I , the corresponding x n+1 will map past x f on either sides, as given by the mappings If we have access to any (x n , x n+1 ) pair, we can also represent the points through a two-dimensional cartesian coordinate system XY, where the X -axis represents n and n + 1, and the Y -axis represents x n and x n+1 , and one can construct a straight line through coordinates (n, x n ) and (n + 1, x n+1 ). We therefore observe, how the straight lines formed by each corresponding pair of x n and x n+1 interact, for all the x n points originating within I . We show the line plots for a few points x n ∈ [0, 1] and their corresponding x n+1 iterates in Fig. 5a. Evidently, the mapping T : x f ], and from the observation in Fig. 5a, it can be stated that the straight lines connecting x n ∈ [x c , 1] ∈ I and their corresponding x n+1 intersect at a single point (that is the x f of the map), on the Y -axis. For the remaining x n ∈ [x p , x c ) points, even though the corresponding x n+1 maps beyond x f , any intersections with them are spread over a wide range on  1] the XY-plane instead of being concentrated on a single one.
To validate our statement, we consider a further generic situation. Let there be two arbitrary points x n ∈ [x c , 1] ∈ I and x n ∈ [x c , 1] ∈ I , with their corresponding next iterates given by x n+1 = T (x n ) and x n+1 = T x n . If we construct two straight lines by connecting the coordinates (n, x n ) with (n + 1, x n+1 ) and (n,x n ) with (n + 1, x n+1 ) on the XY-plane, then we can propose the following: Theorem 1 For any two points x n ∈ [x c , 1] and x n ∈ [x c , 1] such that x n = x n , and their respective iterates x n+1 and x n+1 , suppose n and n + 1 are plotted on the X -axis and x n , x n+1 , x n , and x n+1 are plotted on the Yaxis. Then the lines joining the coordinates (n, x n ) with (n + 1, x n+1 ), and (n, x n ) with (n + 1, x n+1 ) always intersect each other at a point whose Y -coordinate will be equal to the value of x f , the nonzero fixed point.
Proof Let there be two points x n = (x c + d) and x n = (1−d ); arbitrary distances d and d away from x c and 1, respectively, such that x n ∈ [x c , 1], x n ∈ [x c , 1] and x n = x n . Therefore, the iterates of the chosen points will be, x n+1 = T (x c + d) = μ − 2μd, and x n+1 = T (1 − d ) = 2μd from the properties of the tent map described in Sect. 2.
Since the ordinate value of the point, where the straight lines intersect, gives us the view of the phenomenon, solving only for Y solutions will suffice for our purpose.
Equating (3) and (4), we can write Putting x c = 1/2 in (5) and solving for Y Hence, the proof is complete. It is thus confirmed that for x n ∈ [x c , 1] ∈ I , i.e. 50% of the possible points within the state space I , will have the straight lines formed with their corresponding x n+1 intersecting at the same point x f .

Parameter estimation from the noisy trajectories
From Theorem 1 and previous observations, it is now evident that location of the nonzero fixed point of the map can be identified from the intersections (Y solutions of (5)) exhibited by straight lines formed between ï m n and ï m n+1 samples of the noisy distribution of trajectories. We apply the proposed idea to determine the parameter from the noisy dynamics of the tent map and eventually adapt it for the other maps.
For the tent map, we intend to determine a point within the concentrated neighbourhood of crossovers, as the closest estimate of x f , which in turn, is indicative of the map parameter as given by (2). To determine the intersections of the straight lines formed between the nth and (n + 1)th iterates, we select the set of points, Thus, for all i and j such that i = j, there will be M (M − 1)/2 number of intersections whose solutions of the ordinate value Y k is given by rewriting (5) in terms of ï i n , ï i n+1 and ï j n , ï j n+1 where k = 1, 2, . . ., M (M − 1)/2. The Y k solutions form a cluster of the points of intersection between the lines joining the n and (n + 1) time steps. To determine a central point within such a unidimensional cluster, we calculate the arithmetic meanȲ n of all the Y k solutions between each n and (n + 1) time steps. Despite the selection criterion ï m n ∈ [x c , 1], in order to have at least one Y k solution for an intersection, there must be at least two elements in H n ; hence, we exclude any such |H n | < 2.
To illustrate the concentration of crossovers, we consider the case: x 0 = 0.383, iterated for N = 50 with μ = 0.715 and the iterates were perturbed by AWGN with a noise level of SNR = 20 dB, and M = 200 samples of η m trajectories were collected. We show the distribution of Y k through histograms for n = 16 and n = 20, respectively, in Fig. 6a, b, with calculated meanȲ n and standard deviation D n .
For our chosen case of μ = 0.715, the corresponding x f = 2μ/(1 + 2μ) = 0.5884770. It has been noticed that all theȲ n values are the close approximations of the said x f . Hence, we can determine a single point ξ by calculating the arithmetic mean of the collection ofȲ n values, to derive the closest estimate of x f . Using ξ, we therefore estimate the control parameter of the tent map. The estimated parameter μ is given by rewriting (2) in terms of ξ and μ For the above-mentioned case, we determined the ξ = 0.5888251 and standard deviation SD = 0.0086253 for the collection ofȲ n values, and we find the estimated μ = 0.7160277 which is significantly close to the actual μ = 0.715. Next, we explore the effectiveness of the technique further by applying it to the asymmetric (skew) tent map, zigzag map and Bernoulli shift map which are more practical and widely implemented cases of 1D PWL maps.

Skew tent map
Skew tent map S (x n ) can be defined [4] as where ν ∈ (0, 1) defines the map partition that controls the position of the peak-and thereby the slope on either side of the partition-causing the map to appear as skewed or asymmetric when ν = 0.5 as given in Fig. 7. In Fig. 8

Zigzag map
The zigzag map Z (x n ) is defined [5] as where ω is the parameter of the map ranging from ω ∈ [−3, 3], and the dynamical state space of x n being [−1, 1].
In the zigzag map, the map parameter ω can vary within the range of positive and negative values. In case of a positive parameter, along with the fixed point x f Z = 0, the map has a pair of nonzero fixed points: x f N = − 2 (1+ω) for the negative domain and x f P = 2 (1+ω) for the positive domain as shown in Fig. 10a. When ω is negative, the map intersects the diagonal only at x f Z (shown in Fig. 10b). As a result, for + ω, the crossovers between the iterates appear at both negative [−1, 0) and positive [0, 1] ranges of the map (Fig. 11a), while that in case of − ω appear only around zero, as can be seen in Fig. 11b. Such crossovers around zero, however, are rendered ineffective for our purpose because, when the parameter is changed, the resulting function simply readjusts the slope of the map about zero, and hence, no change in parameter values is reflected by the position of x f Z .
This problem can be resolved by addressing the fact that the negative parameter reverses the sign in every odd count of the map iteration [5]. As a result, − ω causes ∓x n+1 = Z(± x n ) and ∓x n+3 = Z(Z(Z(± x n ))) for odd counts, while generating ± x n+2 = Z(Z (± x n )) for the even count. Therefore, the trajectory generated with − ω can be transformed into its positive equivalent by force-reversing the signs of the odd count of iterates. This can be verified through Fig. 12, where the time series of an arbitrary initial condition iterated with ω = −2 has been shown; every alternating iterate was multiplied by −1 for the forcedreversal, effectively converting the recorded time series to be congruent to the dynamics by Z (x n ) for ω = 2. As a result, in case of noisy dynamics, the meaningful crossovers around x f P and x f N can re-emerge when the signs of the alternate iterate in the time series generated by negative parameter is reversed. Since, whether the actual parameter is positive or negative will not be known-to be able to suitably adapt the proposed parameter estimation approach for the zigzag mapwe record the original samples of the noisy dynamics and also produce an alternately inverted copy of the recorded data. In either of the two datasets, the desired crossovers will emerge about the nonzero fixed points depending on the sign of the parameter. If the nonzero crossover solutions are found in the inverted data set, then it can be assured that the parameter was negative, and the magnitude of ω can be determined.
The simultaneous presence of x f P and x f N may instigate further modifications to enhance the performance of the proposed method. From the map definition in Eq. (9), it is seen that x f P = x f N (Fig. 10a) which results into symmetrically mirrored crossovers around x f P and x f N in the noisy dynamics. Therefore, Since the Y k solutions are generated on both positive and negative domain, the negative solutions are inverted before calculating the respectiveȲ n average. Once theȲ n points for all n have been computed, the positive fixed point of the map was estimated as ξ = 0.5694459. Replacing x f P with ξ in x f P = 2 (1+ω) , the parameter is estimated to be ω = 2−ξ ξ = 2.5121858, which is in close agreement with the chosen ω = 2.5 with −1.2% estimation error.
When the parameter is negative, the nonzero Y k solutions are found in the force-reversed dataset. For a case chosen with ω = −2.5, the fixed point was estimated from the alternately reversed dataset as ξ = 0.5684757. As a result, the sign of the parameter that has been determined must be reversed back, and therefore, the estimated parameter is given by ω = − 2−ξ ξ = − 2.5181797.

Bernoulli shift map
The Bernoulli shift map B(x n ) is defined [6] by (10).
where β is the parameter of the map, which has a chaotic range of β ∈ [0.7, 1] while x n ∈ [− 0.5, 0.5].
In case of Bernoulli shift map, it does not intersect the diagonal anywhere except at 0 and 1; therefore, the desired cluster of crossovers within the noisy trajectories do not appear in the neighbourhood of any dynamical fixed point, as given by the map definition which can also be verified from Fig. 14a. However, inverting the sign of the alternate iterates results in the map behaving as −B(x n ), where it intersects the diagonal at the points x f P = 0.5 1+2β and x f N = −0.5 1+2β , shown in Fig. 14b.
As a result, prominent crossovers appear at x f P and x f N in the noisy samples when odd iterates are inverted, which can be treated in a similar fashion as the zigzag map.
To demonstrate, M = 50 samples of noisy (SNR = 20 dB) time series for the Bernoulli map is generated with initial condition x 0 = − 0.23684 and parameter β = 0.85, for N = 50 iterates, and the initial 20 iterations (after inverting the odd iterates) are shown in Fig. 15. The ± Y k solutions are found, from which the −Y k solutions are inverted to reinforce the +Y k values. The concentration of Y k between iterations 12 and 13 is shown in Fig. 16

Confidence interval
In many cases, the harsher conditions of dynamical noise can cause the system to depart from normal distribution [10], as it may introduce some systematic error in the statistical estimates. Therefore, to further ensure the efficacy of the approach, the estimation experiment has been repeated for the individual maps over a range of SNR values in order to determine the error bar for each case of SNR. For each case of the maps, we notice that the estimated parameter values at SNR values 10 dB or less, are quite apart from the chosen test parameters. However, better results are achieved for slightly improved SNR values as relatively lower  for each map μ = 0.85, ν = 0.25, ω = 3, and β = 0.9, respectively. It shows the mean value of the estimated μ for all 25 attempts with 95% confidence intervals for each SNR condition noise may still preserve the qualitative properties of the trajectories [31]. We show in Fig. 17a-d the error bar estimates for tent map (μ = 0.85), skew tent (ν = 0.25), zigzag (ω = 3) and Bernoulli maps (β = 0.9), respectively, with the corresponding mean of parameter estimates.
Let us say for tent map, μ mean = 1 Q Q q=1 μ q of Q = 25 independent estimation attempts; μ q for q = 1, 2, . . . , Q is computed for each case of noise over a range of SNR values 15-30 dB. To estimate the standard error bar, we calculate the upper and lower bound with 95% confidence interval given by μ mean ± 1.96 μ SD √ Q , respectively, where μ SD is the standard deviation of Q estimation attempts. Similarly, we determine the standard error bars ν mean ± 1.96 ν SD √ Q , ω mean ± 1.96 ω SD √ Q and β mean ± 1.96 β SD √ Q for the skew tent, zigzag and Bernoulli shift maps, respectively. It can be seen that the estimated outcomes in the chosen cases of the maps start to improve from SNR 15 dB onwards.

Conclusion
In this paper, we proposed a method to estimate the control parameter of the 1D PWL maps from noisy dynamics of the system. Due to the inherent noise in physically implemented chaotic systems the dynamical behaviour of the function becomes greatly affected. Even a small amount of noise distorts the map trajectories to a great extent. For the effective retrieval of the meaningful trajectories, existing noise reduction techniques require additional information about either the dynamical system definition or the noise level within the system. Given the computational complexity that it might add while approximating the system model or determining the noise level, we chose to estimate the control parameter of the system, whose value fully defines the behaviour of the map.
We have shown that the iterates of the dynamical time series can be treated as points on a twodimensional Cartesian coordinate system. As the trajectories become highly distributed over the state space due to the presence of dynamical noise in the system, we show that the straight lines connecting the (X, Y ) coordinates of the consecutive iterates of all the sampled time series form a cluster of intersections between the nth and (n + 1)th iterates. We have shown that such clusters appear in the close neighbourhood of the unstable nonzero fixed point of the maps (tent map, skew tent map, zigzag map and Bernoulli shift map) considered. We have established the fact that the straight lines connecting the points (for tent map originating from [0.5, 1]) with their consecutive iterates will always intersect each other at the nonzero fixed point. Utilising this property, we have shown how the value of the fixed point can be estimated from the clusters and the parameter of the 1D PWL maps can be determined.
The effectiveness of the proposed method has been studied through numerical simulations where the esti-mations have been tried with various cases of parameters and noise levels for each case of 1D PWL maps. We have conducted a statistical analysis for the behaviour of the algorithm by showing the consistency of the estimated outcomes through numerous attempts, for a range of SNR values for each map described. We have shown for a range of SNR values (e.g. [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], the parameter estimates are satisfactorily close and that can further improve the system approximations through the mentioned noise reduction techniques.