# A New Fractional-Order Map with Infinite Number of Equilibria and Its Encryption Application

Authors: Ahlem Gasri; Amina-Aicha Khennaoui; Adel Ouannas; Giuseppe Grassi; Apostolos Iatropoulos; Lazaros Moysis; Christos Volos
Journal: Complexity
 (2022)

---

## Abstract

The study of the chaotic dynamics in fractional-order discrete-time systems has received great attention over the last years. Some efforts have been also devoted to analyze fractional maps with special features. This paper makes a contribution to the topic by introducing a new fractional map that is characterized by both particular dynamic behaviors and specific properties related to the system equilibria. In particular, the conceived one dimensional map is algebraically simpler than all the proposed fractional maps in the literature. Using numerical simulation, we investigate the dynamic and complexity of the fractional map. The results indicate that the new one-dimensional fractional map displays various types of coexisting attractors. The approximate entropy is used to observe the changes in the sequence sequence complexity when the fractional order and system parameter. Finally, the fractional map is applied to the problem of encrypting electrophysiological signals. For the encryption process, random numbers were generated using the values of the fractional map. Some statistical tests are given to show the performance of the encryption.


---

## Body

## 1. Introduction

Fractional calculus is a topic which is developed more than 300 years. However, it is only the last decades that it has been extensively and intensively investigated, due to its wide application in signal mechanical controls and other fields [1]. Compared with integer order derivatives, fractional-order derivatives are more accurate as they provide excellent tool for the description of the memory effect in all kinds of materials and processing. Based on this consideration, the application of fractional-order systems have attracted more and more researchers attention [2]. At the same time, during the last decade, attention has been focused on discrete fractional calculus and fractional difference operators [3, 4]. Several papers regarding the presence of chaotic phenomena in fractional discrete-time systems (maps) have been published to date [5–9]. For example, in [10] the hyperchaotic dynamic of the fractional generalized Hénon map has been investigated, whereas in [11] the presence of chaos in the fractional discrete memristor system has been illustrated. In [6] the presence of chaos in the fractional sine map and in the fractional standard map has been analyzed in details. In [7], control laws for stabilizing the chaotic dynamics of the fractional Grassi-Miller map have been developed, whereas in [8] the fractional Hénon map and its chaotic attractors have been studied. In [12], the chaotic dynamics of three maps (i.e., the fractional flow map, the fractional Lozi map and the fractional Lorenz map) have been investigated, whereas in [13] the chaotic behavior of the fractional Tinkerbell map has been illustrated. To our knowledge, all of the above reported systems have a finite number of equilibria. In general, systems with more equilibria may bring unexpected stabilities to some extent. Therefore, several efforts have been devoted to the study of fractional chaotic maps with some special features related to the system equilibria [14]. Among these studies, Zambrano-Serrano et al. [15] analyzed the dynamic properties and projective synchronization of the fractional difference map with no equilibrium, whereas in [16] Almatroud et al. found rich chaotic behaviours of a novel two-dimensional (2D) hyperchaotic fractional map with infinite line of equilibrium As a result, the analysis of chaotic dynamical behaviours of the fractional-order discrete-time systems without equilibrium points is an interesting topic.In recent years, there has been a growing interest in discrete time systems with special complex dynamical behaviours, such as hidden attractors [17], coexisting multiple attractors [18] and hyperchaotic behaviours. For instance, in [19] the dynamic properties of a novel memristive hyperchaotic map and its application in secure communication have been illustrated. A 2D sine map was presented in [20] and several interesting behaviour like coexisting attractors and initial offset boosting were explored. General speaking, chaotic maps with coexisting attractors has drawn the attention of many researcher. Coexistence of attractors is a special phenomenon in nonlinear dynamical systems, which denotes that with fixed values of system parameters, a tiny disturbance in the initial condition can lead to the coexistence of different attractors. This property makes the chaotic maps very useful in the fields of secure communication and encryption. Since such phenomenon has not received enough attention with fractional discrete-time systems [21], this paper aims to make a contribution by introducing a new fractional map that is characterized by both particular dynamic behaviors and specific properties related to the system equilibria. Namely, the proposed map possesses infinite number of equilibria in a bounded domain, being this a new feature for fractional map, not published in literature to date. Dynamics and complexity of the conceived map are investigated in details. In particular, bifurcation diagrams, maximum Lyapunov exponents and 0-1 test are reported to highlight the coexistence of different periodic and chaotic attractors. Moreover, the map is applied to the problem of data encryption, which is a well established application of chaotic systems. As recent examples of fractional systems being applied to encryption, in the work [22] a technique to improve chaotic behavior in fractional maps is proposed, and applied to image encryption. An image encryption using the fractional logistic map is proposed in [23], and a different one in [24]. The 2D fractional Hénon map was also applied to image encryption in [25]. In [26] a graphical user interface is designed for random number generators based on integer and fractional order chaotic systems. In [27] another pseudo random number generator was designed based on the coupling of multiple fractional chaotic systems. The above works are just a small sample of the expanding use of continuous and discrete fractional systems to encryption. The use of fractional systems is drawing increasing attention by researchers, since they have more complex dynamics due to their memory effect, and a higher key space than their integer order analogues, since the fractional order and the finite memory order constitute additional key values.In this work, the encryption of electrophysiological signals [28–33] is considered. To do so, first, a pseudo-random bit generator (PRBG) is designed using the values of the chaotic map. This chaotic PRBG is the basis for the encryption design, since it is used as the source of deterministic randomness [34–38]. Here, to take advantage of the fractional nature of the map that is used as the basis of the PRBG, a technique is proposed that takes into account its memory effect. So first, the map is implemented using finite memory, in order to reduce its computational cost. Then, in each iteration, all previous values of the map are taken into account in generating the bits, which leads in producing 459 bits per iteration. So the proposed technique reduces the computational load of the PRBG, since fewer iterations of the map are required to reach a desired bistream length, and is specifically designed for fractional order systems.After the PRBG is designed, to encrypt a given electrophysiological signal, two rounds of masking are performed. First, the signal is modulated to mask its structure by combining it with the values of the fractional chaotic map. Then, the modulated signal is transformed into its binary representation, and combined with a bitstream generated from the PRBG, to yield the encrypted signal. This binary signal can then be safely transmitted through a communication channel. The original signal can be retrieved back at the receiver end, by following the reverse encryption procedure. The performance of the encryption is finally evaluated by a series of statistical tests performed on the original, modulated, and encrypted signals.Finally, the encryption process is realised in a microcontroller board. This implementation helps verify the feasibility of simulating fractional maps in low cost hardware devices, a task that is of high interest, due to the potential implementability of chaotic systems in IoT related devices [39, 40]. Such realizations have already been explored in the literature, with positive results. For example, in [41] a fractional macro-economic model was established, and a microcontroller implementation was designed on an Arduino Due. In [42] the problem of impulsive synchronization for fractional order discrete chaotic maps is considered, and the design was again implemented in two Arduino Mega boards that simulated the master and observer systems. In [43] a hyperchaotic fractional order system was proposed and discretized. The system was then simulated in an Arduino Uno board and applied to the problem of text encryption. In [44], a fractional order modified Chua’s circuit was designed and implemented in an Arduino Uno microcontroller. In [45] a generalized fractional logistic map was constructed and applied to random number generation and image encryption, implemented on a Virtex-5 field-programmable gate array FPGA. In the current work, the realization is done on an STM32F103 nucleo development board, and the encrypted signal from the microcontroller has the same statistical properties to the signal generated from Matlab, as indicated from all the statistical tests performed.


### 1.1. Fractional Calculus and Preliminaries

In this section, some preliminaries and basic concepts associated with discrete fractional calculus are presented here for completeness. In the following we consider that our domain is the time scaleℕb=b,b+1,b+2,… with b∈ℝ. Let X denote any function defined from ℕb, thus the μ-th fractional sum for μ>0 is defined by [46]:(1)Δb−μXs=1Γμ∑l=bs−μs−l−1μ−1Xl,where s∈ℕb+μ. Observe that the term sμ indicates to the so-called a falling function which may be defined via Gamma function, Γ, as follows:(2)sμ=Γs+1Γs+1−μ.Based on the above definition of theμ-th fractional sum, it is possible to define the μ-Caputo like difference operator. Let X denote any function defined from ℕb. The Caputo difference operator with order μ∉ℕ is defined by:Definition 1.
Forn=⌈μ⌉+1, theμ-th order Caputo-like operator can be defined as [47]:(3)ΔCbμXs=Δb−n−μ,ΔnXs=1Γn−μ∑l=bs−n−μs−l−1n−μ−1ΔlnXl.
Fors∈ℕb+n−μ.
Now a theorem is briefly summarized, in order to derive in the following the discrete formula of the new fractional map.Theorem 1.
[48] For the fractional difference equation(4)ΔCμbμXs=fs+μ−1,Xs+μ−1,ΔrXb=Xr,n=⌈μ⌉+1,r=0,1,…,n−1.Then, the discrete integral equation which is equivalent to equations in (4) is:(5)Xs=Xb+1Γμ∑l=b+n−μs−μs−σlμ−1fl+μ−1,Xl+μ−1,t∈ℕb+n.Note that, for the purpose of numerical calculation, (5) can be written as [48]:(6)Xn=X0+1Γμ∑i=0n−1Γn−i+μΓn−i+1fi,Xi.


## 1.1. Fractional Calculus and Preliminaries

In this section, some preliminaries and basic concepts associated with discrete fractional calculus are presented here for completeness. In the following we consider that our domain is the time scaleℕb=b,b+1,b+2,… with b∈ℝ. Let X denote any function defined from ℕb, thus the μ-th fractional sum for μ>0 is defined by [46]:(1)Δb−μXs=1Γμ∑l=bs−μs−l−1μ−1Xl,where s∈ℕb+μ. Observe that the term sμ indicates to the so-called a falling function which may be defined via Gamma function, Γ, as follows:(2)sμ=Γs+1Γs+1−μ.Based on the above definition of theμ-th fractional sum, it is possible to define the μ-Caputo like difference operator. Let X denote any function defined from ℕb. The Caputo difference operator with order μ∉ℕ is defined by:Definition 1.
Forn=⌈μ⌉+1, theμ-th order Caputo-like operator can be defined as [47]:(3)ΔCbμXs=Δb−n−μ,ΔnXs=1Γn−μ∑l=bs−n−μs−l−1n−μ−1ΔlnXl.
Fors∈ℕb+n−μ.
Now a theorem is briefly summarized, in order to derive in the following the discrete formula of the new fractional map.Theorem 1.
[48] For the fractional difference equation(4)ΔCμbμXs=fs+μ−1,Xs+μ−1,ΔrXb=Xr,n=⌈μ⌉+1,r=0,1,…,n−1.Then, the discrete integral equation which is equivalent to equations in (4) is:(5)Xs=Xb+1Γμ∑l=b+n−μs−μs−σlμ−1fl+μ−1,Xl+μ−1,t∈ℕb+n.Note that, for the purpose of numerical calculation, (5) can be written as [48]:(6)Xn=X0+1Γμ∑i=0n−1Γn−i+μΓn−i+1fi,Xi.


## 2. The New Fractional Map

Let us consider the following one dimensional map, which has been recently proposed in [49] as an example of an elegant map that can display chaotic behavior:(7)xn+1=AsinCxn+B,where A,B and C are three positive parameters. In this work, we extend the original integer–order system to the fractional-order case. Specifically, we consider the effect of the fractional–order in the system dynamics by introducing the following new fractional difference equation:(8)Δbμxs=AsinCxs+1−μ−xs+1−μ+B,where s∈ℕb+1−μ and μ is the fractional order with 0<μ≤1. This fractional map is invariant under transformation x⟶−x for all values of parameters A,B,C and order μ. Hence this map could display coexisting attractors for appropriate choice of initial conditions and fractional order values as well. In order to investigate this property, the numerical formula is designed as:(9)xn=x0+1Γμ∑i=0n−1Γn−1−i+μΓn−iAsinCxi−xi+B.According to the discrete (9), the proposed fractional map (8) has memory effect, which means that the iterated solution xn is determined by all the previous states. In order to find the equilibrium points xf of the fractional map (8), we solve the following equation(10)ApEn=ϕmr−ϕm+1r,Obviously, (10) is a trigonometrical equation that is very difficult to solve analytically, therefore, to analyze the equilibrium points we adopt the graphic analytic method. Taking the system parameters A=0.5,C=3,B=0 as an example and the integer order value μ=1, Figure 1 shows the phase portrait obtained by simulating (9) in the xn−xn+1 plane along with the y=x line. The black line in Figure 1, bisects the first and third quadrants at the interval −0.5,0.5, and its crosses with the map are equilibrium points. As one can see, Figure 1 visually demonstrates that there are infinite equilibrium points in this case. On the other hand, Figure 2 shows different phase portraits obtained by simulating (9) for different values of order μ. In comparison to the integer-order case, the shape of the bounded attractors does not change much with the decrease of the value of μ, while the interval where equilibrium points exist is changed from −0.5,0.5 into −0.6,0.6.Figure 1
Phase diagram of fractional map (8) for μ=1,A=0.5,B=0,C=1.Figure 2
Different chaotic attractors of the fractional map (8) for system parameters A=0.5,B=0,C=3 and different values of μ: (a)μ=0.9 and (b)μ=0.8. (c)μ=0.8.
(a)(b)(c)


## 3. Dynamics and Complexity of the New Fractional Map

In this Section the dynamic behaviors of the fractional-order map (8) are numerically investigated using bifurcation diagrams and computation of maximum Lyapunov exponents, 0-1 test and entropy. Namely, the influence of both fractional order and initial conditions on the dynamics of the novel map are investigated in detail, with the aim to highlight the coexistence of different chaotic attractors. Namely, the influence of both fractional order and initial conditions on the dynamical behaviour of the novel map are investigated in detail, with the aim to highlight the coexistence of different chaotic attractors. To give a finer analysis of our fractional map we model its dynamics for two parameter sets, i.e for the above values A=0.5, B=0,C=3 and for A=0.05,B=0, C=1.


### 3.1. Bifurcation Diagrams and MLE

By changing the system parameters, the fractional map (8) can undertake different dynamic scenarios. Figure 3 gives the three-dimensional view of the bifurcation diagrams of the fractional map (8) for different values of μ (i.e., μ=1,μ=0.6,μ=0.2). We consider the parameter A as the bifurcation parameter and we take the other parameters as B=0 and C=1. We consider two values of symmetric initial conditions: x0=0.1 and the negative one x0=−0.1. Clearly, these diagrams are different (see Figure 4). In particular, in Figure 4 the diagram in magenta color represents the dynamic behavior of the fractional map (8) for order μ=1 (i.e., for the integer-order case). As it can be seen, the fractional map is chaotic over most of the range and is changed from chaotic to periodic with the increase of A. It is worth noting that the phenomenon of coexisting attractor is not observed in the integer-order case, where the states of our system come in the form of individual trajectories. By changing the order μ from 1 to 0.6, the bifurcation diagrams become the ones depicted in brown and green colors. In this case, coexisting periodic orbits appear within the interval 0.1323,0.1723∪0.489,0.517, along with a symmetry breaking bifurcation.Figure 3
Different bifurcation diagrams of the fractional map (8) in three-dimensional space with the variation of system parameter B: blue diagram forμ=1; red diagram for μ=0.8; magneta diagram for μ=0.6; green diagram for μ=0.4.Figure 4
Different bifurcation diagrams of the fractional map (8) in three-dimensional space with the variation of system parameter A, and for two initial conditions: blue diagram for μ=0.2 and for initial condition x0=0.1; red diagram μ=0.2 and for initial condition x0=−0.1; brown and green diagram for order μ=0.6; magenta diagram for order μ=1.Meanwhile, the region of chaotic motion increases with the value ofA. Finally, the bifurcation diagrams with μ=0.2 are considered and plotted in red and blue colors for different initial conditions. Specifically, the states colored in blue are those starting from the positive initial condition, whereas the states colored in red are those starting from the negative initial condition. In this case, the fractional map (8) displays symmetric separate coexisting diagrams. Figure 5 shows some coexisting chaotic and periodic attractors for different values of the system parameter A. Specifically, when we select the parameter A=0.04, the fractional map (8) highlights the coexistence of two symmetric chaotic attractors corresponding to the positive and negative initial conditions, respectively. Similarly, when we choose the parameter A=0.2, the phase plot in Figure 5 shows the coexistence of two periodic symmetric attractors corresponding to the positive and negative initial conditions, respectively. The analysis conducted so far clearly highlights that the dynamics of the fractional map (8) become very complex when the value of the order μ is significantly decreased.Figure 5
Two coexisting chaotic attractors forA=0.04 and μ=0.2 (a); two coexisting periodic attractors for A=0.2 and μ=0.2 (b).
(a)(b)(b) In order to deeply investigate the properties of coexisting attractors in the fractional map (8), bifurcation diagrams and computations of maximum Lyapunov exponents (LE) with respect to the fractional order are carried out. Various dynamic behaviors can be observed by changing the value of μ at the interval 0,1. By considering the values B=0,A=0.5 and C=3, Figure 5(a) displays the bifurcation diagram of the state variable xn corresponding to the positive initial condition (blue diagram) and the negative initial condition (red diagram). For C=3, the trajectories of the fractional map change from chaos to coexisting chaotic attractors with some periodic orbits through period doubling route to chaos. When μ≤0.31 the chaotic behavior disappears and separate coexisting periodic orbits appear. This indicates that the dynamics of system (8) for C=3 become simpler as the value of μ decreases. This result is also confirmed by the maximum Lyapunov exponents reported in Figure 5(b). Now the bifurcation diagram and the MLE of the fractional map (8) are illustrated in Figure 7 for A=0.05,B=0 and C=1. Similarly, the dynamic behavior of the map (8) depends on the initial condition and the value of μ. The computation of the MLE shows that the states of the fractional map (8) go from chaotic to periodic with the decrease of order μ. From both the Figures 6 and 7 it can be deduced that the property of having coexisting attractors is observed when the fractional order assumes small values. In order to further investigate this property, we choose to plot the phase portraits of the fractional map (8) for multiple initial conditions and for the same parameters used in Figure 6. The obtained results are plotted in Figure 8. When μ=0.9, a chaotic attractor is observed. When μ=0.6, there are two symmetric separate chaotic attractors and a chaotic attractor in magenta color. On the other hand, four coexisting attractors are obtained when μ=0.3, as depicted in Figure 8(c). Suppose now that A=0.05,B=0,C=1 and μ=0.3, then the corresponding phase diagram for the map (8) is shown in Figure 9 using different initial conditions. In particular, Figure 9 reveals the presence of four coexisting chaotic attractors, i.e., the blue attractor for x0=0.1, the red attractor for x0=−0.1, the green attractor for x0=2 and the magenta attractor for x0=−2. Figure 9 confirms the rich dynamics of the conceived fractional map, indicating that a number of different coexisting chaotic attractors could be found by taking other suitable values of the fractional order as well as other proper initial conditions.Figure 6
Bifurcation and maximum LE of system (8) versus μ depicted in (a) and (b), respectively, for parameters A=0.5,B=0,C=3, and two symmetric initial conditions (x0=0.1 for the blue diagram and x0=−0.1 for the red diagram).
(a)(b)Figure 7
Bifurcation and maximum LE of system (8) versus μ depicted in (a) and (b), respectively, for parameters A=0.05,B=0,C=1, and two symmetric initial conditions (x0=0.1 for the blue diagram and x0=−0.1 for the red diagram).
(a)(b)Figure 8
The coexisting attractors of the fractional-order map (8) with parametersA=0.5,B=0,C=3 and with the initial condition x0=−0.1 for the red attractor and x0=0.1 for the blue attractor; (a) chaotic attractor for μ=0.9, (b) three coexisting chaotic attractors for μ=0.6 with initial condition x0=10 for the magenta attractor; (c) four coexisting chaotic attractors for μ=0.3 with initial condition x0=10 for the magenta attractor and x0=−10 for the black attractor.
(a)(b)(c)Figure 9
Four coexisting chaotic attractors for different initial conditions and system parametersA=0.05,B=0,C=1 and μ=0.3: blue attractor for x0=0.1, red attractor for x0=−0.1, green attractor for x0=2 and magenta attractor for x0=−2.


### 3.2. 0-1 Test

Another tool that can be used to study the influence of the fractional order on the dynamic of a fractional map is “0-1 test.” This test, proposed in [50] for fractional-order systems, is able to check the presence of chaos in a series of data that originates from a deterministic system. For the fractional-order model, consider a set of data xn where n=1,…N. Using the approach in [50], we transform the trajectories of the fractional-order map into p−q plots. Generally, unbounded p−q trajectories imply chaotic behavior, whereas bounded trajectories implies regular behavior. Herein, we apply the 0−1 test method directly to the solution xn that is obtained from the discrete formula (9). Herein, we simulate the translation components of the system (8) in the p−q plane. By taking A=0.5,B=0,C=3 and by varying the value of μ, the results of the application of the 0-1 test to the fractional map (8) are reported in Figure 10. In particular, Figure 10(a) depicts the Brownian–like trajectories for all initial conditions, indicating that the suggested map is chaotic for μ=0.9. On the other hand, Figure 10(b) depicts the Brownian–like trajectories for two different initial conditions, confirming the coexistence of chaotic attractors for μ=0.6. Finally, when μ=0.3 the coexistence of chaotic attractors and periodic orbits is confirmed by the plot in Figure 10(c), which depicts bounded–like trajectories for the initial state x0=0.1 and Brownian–like trajectories for the initial state x0=2. It is concluded that the 0-1 test proves to be a valuable tool for checking the coexistence of different chaotic attractors as well as the coexistence of periodic one in the proposed fractional map (8).Figure 10
0-1 test of the fractional map (8) for system parametersA=0.5,B=0 and C=3: (a) Brownian-like trajectories for μ=0.9; (b) Brownian-like trajectories for μ=0.6 and for two different initial condition (red color for x0=10 and blue color for x0=0.1); (c) Brownian-like trajectories for μ=0.3 and initial condition x0=10 (in red), along with bounded like trajectories for x0=0.1 (in blue).
(a)(b)(c)


### 3.3. Approximate Entropy

The approximate entropy (ApEn) [51] is the measurement of the degree of complexity of a series of data from multi-dimensional perspective. This method estimates the regularity by assigning a non-negative number, where higher values indicate higher complexity. By applying the technique in [52], we consider a set of points x1,…,xN that are obtained from the discrete formula (9). The value of the approximate entropy depends on two important parameters, i.e., m and r, where the input r is the similar tolerance whereas m is the embedding dimension. Here, to calculate ApEn values we set m=2 and r=0.2stdx where stdx presents the standard deviation of the data x. Those values are preferred values according to similar previous studies. We reconstruct a subsequence of x such that Xi=xi,…,xi+m−1, where m presents the points from xi to xi+m−1. Let K be the number of xi such that the maximum absolute difference of two vectors xi and Xj is lower or equal to the tolerance r. The relative frequency of xi being similar to Xj is given by(11)ϕmr=1n−m−1∑i=1n−m+1logCimr.TheApEn is calculated as following(12)b8n:8n+7=de2bimod⌊1015xn⌋,256,where ϕmr is denoted to be(13)xn=x0+1Γμ∑i=1mΓm−i+μΓm−i+1AsinCxn−m+i−1−xn−m+i−1+B,We apply theApEn method to simulated complexity of the fractional map (8) by varying the fractional order and control parameter from 0 to 10. For each value of order μ, we analyze a series of points of length 3500 with different values of C. Figure 11 presents the 3D plots of the approximate entropy for fractional map (8). It is observed that, the fractional map (8) can have a higher complexity with relatively larger parameters values C and fractional order μ; which consist with the previous results. The analysis results in Figure 11 indicate that the change of system parameter C has little effect on the complexity. Thus, in real applications, we should be careful for the range of parameters and fractional order μ.Figure 11
The approximate entropy ApEn of the fractional map (8) in three-dimensional space with the variation of system parameterC and fractional order μ, for A=0.5 and B=0.


## 3.1. Bifurcation Diagrams and MLE

By changing the system parameters, the fractional map (8) can undertake different dynamic scenarios. Figure 3 gives the three-dimensional view of the bifurcation diagrams of the fractional map (8) for different values of μ (i.e., μ=1,μ=0.6,μ=0.2). We consider the parameter A as the bifurcation parameter and we take the other parameters as B=0 and C=1. We consider two values of symmetric initial conditions: x0=0.1 and the negative one x0=−0.1. Clearly, these diagrams are different (see Figure 4). In particular, in Figure 4 the diagram in magenta color represents the dynamic behavior of the fractional map (8) for order μ=1 (i.e., for the integer-order case). As it can be seen, the fractional map is chaotic over most of the range and is changed from chaotic to periodic with the increase of A. It is worth noting that the phenomenon of coexisting attractor is not observed in the integer-order case, where the states of our system come in the form of individual trajectories. By changing the order μ from 1 to 0.6, the bifurcation diagrams become the ones depicted in brown and green colors. In this case, coexisting periodic orbits appear within the interval 0.1323,0.1723∪0.489,0.517, along with a symmetry breaking bifurcation.Figure 3
Different bifurcation diagrams of the fractional map (8) in three-dimensional space with the variation of system parameter B: blue diagram forμ=1; red diagram for μ=0.8; magneta diagram for μ=0.6; green diagram for μ=0.4.Figure 4
Different bifurcation diagrams of the fractional map (8) in three-dimensional space with the variation of system parameter A, and for two initial conditions: blue diagram for μ=0.2 and for initial condition x0=0.1; red diagram μ=0.2 and for initial condition x0=−0.1; brown and green diagram for order μ=0.6; magenta diagram for order μ=1.Meanwhile, the region of chaotic motion increases with the value ofA. Finally, the bifurcation diagrams with μ=0.2 are considered and plotted in red and blue colors for different initial conditions. Specifically, the states colored in blue are those starting from the positive initial condition, whereas the states colored in red are those starting from the negative initial condition. In this case, the fractional map (8) displays symmetric separate coexisting diagrams. Figure 5 shows some coexisting chaotic and periodic attractors for different values of the system parameter A. Specifically, when we select the parameter A=0.04, the fractional map (8) highlights the coexistence of two symmetric chaotic attractors corresponding to the positive and negative initial conditions, respectively. Similarly, when we choose the parameter A=0.2, the phase plot in Figure 5 shows the coexistence of two periodic symmetric attractors corresponding to the positive and negative initial conditions, respectively. The analysis conducted so far clearly highlights that the dynamics of the fractional map (8) become very complex when the value of the order μ is significantly decreased.Figure 5
Two coexisting chaotic attractors forA=0.04 and μ=0.2 (a); two coexisting periodic attractors for A=0.2 and μ=0.2 (b).
(a)(b)(b) In order to deeply investigate the properties of coexisting attractors in the fractional map (8), bifurcation diagrams and computations of maximum Lyapunov exponents (LE) with respect to the fractional order are carried out. Various dynamic behaviors can be observed by changing the value of μ at the interval 0,1. By considering the values B=0,A=0.5 and C=3, Figure 5(a) displays the bifurcation diagram of the state variable xn corresponding to the positive initial condition (blue diagram) and the negative initial condition (red diagram). For C=3, the trajectories of the fractional map change from chaos to coexisting chaotic attractors with some periodic orbits through period doubling route to chaos. When μ≤0.31 the chaotic behavior disappears and separate coexisting periodic orbits appear. This indicates that the dynamics of system (8) for C=3 become simpler as the value of μ decreases. This result is also confirmed by the maximum Lyapunov exponents reported in Figure 5(b). Now the bifurcation diagram and the MLE of the fractional map (8) are illustrated in Figure 7 for A=0.05,B=0 and C=1. Similarly, the dynamic behavior of the map (8) depends on the initial condition and the value of μ. The computation of the MLE shows that the states of the fractional map (8) go from chaotic to periodic with the decrease of order μ. From both the Figures 6 and 7 it can be deduced that the property of having coexisting attractors is observed when the fractional order assumes small values. In order to further investigate this property, we choose to plot the phase portraits of the fractional map (8) for multiple initial conditions and for the same parameters used in Figure 6. The obtained results are plotted in Figure 8. When μ=0.9, a chaotic attractor is observed. When μ=0.6, there are two symmetric separate chaotic attractors and a chaotic attractor in magenta color. On the other hand, four coexisting attractors are obtained when μ=0.3, as depicted in Figure 8(c). Suppose now that A=0.05,B=0,C=1 and μ=0.3, then the corresponding phase diagram for the map (8) is shown in Figure 9 using different initial conditions. In particular, Figure 9 reveals the presence of four coexisting chaotic attractors, i.e., the blue attractor for x0=0.1, the red attractor for x0=−0.1, the green attractor for x0=2 and the magenta attractor for x0=−2. Figure 9 confirms the rich dynamics of the conceived fractional map, indicating that a number of different coexisting chaotic attractors could be found by taking other suitable values of the fractional order as well as other proper initial conditions.Figure 6
Bifurcation and maximum LE of system (8) versus μ depicted in (a) and (b), respectively, for parameters A=0.5,B=0,C=3, and two symmetric initial conditions (x0=0.1 for the blue diagram and x0=−0.1 for the red diagram).
(a)(b)Figure 7
Bifurcation and maximum LE of system (8) versus μ depicted in (a) and (b), respectively, for parameters A=0.05,B=0,C=1, and two symmetric initial conditions (x0=0.1 for the blue diagram and x0=−0.1 for the red diagram).
(a)(b)Figure 8
The coexisting attractors of the fractional-order map (8) with parametersA=0.5,B=0,C=3 and with the initial condition x0=−0.1 for the red attractor and x0=0.1 for the blue attractor; (a) chaotic attractor for μ=0.9, (b) three coexisting chaotic attractors for μ=0.6 with initial condition x0=10 for the magenta attractor; (c) four coexisting chaotic attractors for μ=0.3 with initial condition x0=10 for the magenta attractor and x0=−10 for the black attractor.
(a)(b)(c)Figure 9
Four coexisting chaotic attractors for different initial conditions and system parametersA=0.05,B=0,C=1 and μ=0.3: blue attractor for x0=0.1, red attractor for x0=−0.1, green attractor for x0=2 and magenta attractor for x0=−2.


## 3.2. 0-1 Test

Another tool that can be used to study the influence of the fractional order on the dynamic of a fractional map is “0-1 test.” This test, proposed in [50] for fractional-order systems, is able to check the presence of chaos in a series of data that originates from a deterministic system. For the fractional-order model, consider a set of data xn where n=1,…N. Using the approach in [50], we transform the trajectories of the fractional-order map into p−q plots. Generally, unbounded p−q trajectories imply chaotic behavior, whereas bounded trajectories implies regular behavior. Herein, we apply the 0−1 test method directly to the solution xn that is obtained from the discrete formula (9). Herein, we simulate the translation components of the system (8) in the p−q plane. By taking A=0.5,B=0,C=3 and by varying the value of μ, the results of the application of the 0-1 test to the fractional map (8) are reported in Figure 10. In particular, Figure 10(a) depicts the Brownian–like trajectories for all initial conditions, indicating that the suggested map is chaotic for μ=0.9. On the other hand, Figure 10(b) depicts the Brownian–like trajectories for two different initial conditions, confirming the coexistence of chaotic attractors for μ=0.6. Finally, when μ=0.3 the coexistence of chaotic attractors and periodic orbits is confirmed by the plot in Figure 10(c), which depicts bounded–like trajectories for the initial state x0=0.1 and Brownian–like trajectories for the initial state x0=2. It is concluded that the 0-1 test proves to be a valuable tool for checking the coexistence of different chaotic attractors as well as the coexistence of periodic one in the proposed fractional map (8).Figure 10
0-1 test of the fractional map (8) for system parametersA=0.5,B=0 and C=3: (a) Brownian-like trajectories for μ=0.9; (b) Brownian-like trajectories for μ=0.6 and for two different initial condition (red color for x0=10 and blue color for x0=0.1); (c) Brownian-like trajectories for μ=0.3 and initial condition x0=10 (in red), along with bounded like trajectories for x0=0.1 (in blue).
(a)(b)(c)


## 3.3. Approximate Entropy

The approximate entropy (ApEn) [51] is the measurement of the degree of complexity of a series of data from multi-dimensional perspective. This method estimates the regularity by assigning a non-negative number, where higher values indicate higher complexity. By applying the technique in [52], we consider a set of points x1,…,xN that are obtained from the discrete formula (9). The value of the approximate entropy depends on two important parameters, i.e., m and r, where the input r is the similar tolerance whereas m is the embedding dimension. Here, to calculate ApEn values we set m=2 and r=0.2stdx where stdx presents the standard deviation of the data x. Those values are preferred values according to similar previous studies. We reconstruct a subsequence of x such that Xi=xi,…,xi+m−1, where m presents the points from xi to xi+m−1. Let K be the number of xi such that the maximum absolute difference of two vectors xi and Xj is lower or equal to the tolerance r. The relative frequency of xi being similar to Xj is given by(11)ϕmr=1n−m−1∑i=1n−m+1logCimr.TheApEn is calculated as following(12)b8n:8n+7=de2bimod⌊1015xn⌋,256,where ϕmr is denoted to be(13)xn=x0+1Γμ∑i=1mΓm−i+μΓm−i+1AsinCxn−m+i−1−xn−m+i−1+B,We apply theApEn method to simulated complexity of the fractional map (8) by varying the fractional order and control parameter from 0 to 10. For each value of order μ, we analyze a series of points of length 3500 with different values of C. Figure 11 presents the 3D plots of the approximate entropy for fractional map (8). It is observed that, the fractional map (8) can have a higher complexity with relatively larger parameters values C and fractional order μ; which consist with the previous results. The analysis results in Figure 11 indicate that the change of system parameter C has little effect on the complexity. Thus, in real applications, we should be careful for the range of parameters and fractional order μ.Figure 11
The approximate entropy ApEn of the fractional map (8) in three-dimensional space with the variation of system parameterC and fractional order μ, for A=0.5 and B=0.


## 4. Encryption of Electrophysiological Signals

In this section, an electrophysiological signal will be encrypted using the proposed chaotic map. Here, the signal that will be considered as a sample is a 10 second EEG (electroencephalogram) recordingfp1 of a subject during arithmetic tasks. More details on the other types of EEG signals can be found in the PhysioBank database [53, 54].


### 4.1. Design a PRBG

The first step in the encryption scheme is to design a chaotic pseudo-random bit generator. Since the main computational drawback of fractional order systems is their memory effect, the system will be implemented using a finite memory. Thus, its finite form is given as(14)ℳi=εi+⌊105xi⌋103,i=1,…,ℓ,where m denotes the memory used. As one can see from (14), when i=1,…,m the fractional map states depends on the past m variable, i.e, x11,…,x1m. After some trial and error, the memory is chosen as m=50, which satisfyingly models the dynamical behavior of the system, as can be seen in Figure 12. Although the diagram is not identical to the one of the fractional map shown in Figure 4, it accurately showcases chaotic behavior in the range of around 0.2 to 0.8. Implementing fractional systems using finite memory is a standard approach in reducing computational cost and avoiding overflow, and it is used in most microcontroller implementations of fractional chaotic systems [41, 42, 55].Figure 12
Bifurcation diagram of (14) with respect to parameter A, for B=0, C=1, μ=0.2, m=50, and x0=0.1 (blue), x0=−0.1 (red).In order to take advantage of the fractional nature of the proposed map, each individual term inside the finite sum (14) is used in the bit generation. So in each iteration of the map, the PRBG generates bits using the following rules:(15)εenc−bin=ℳbin⊕ℬ,where ⌊⋅⌋ denotes the floor operation. So in each iteration, each individual term of the sum in (14), as well as the value xn of the map are multiplied by 1015 and the integer part of this product is taken modulo 512. The result of this operation is then transformed to binary representation. Since the outcome of the modulo is an integer from 0 to 511, transforming it to binary generates nine bits per iteration. Thus, each of the terms bn,1,…,bn,51 corresponds to nine bits. The resulting bitstream is taken by concatenating the bits as(16)ℬ=bn−1,1,….,bn−1,51,bn,1,….,bn,51,….With the algorithm, there are 459 bits generated in each iteration. So this approach brings an advantage over classic techniques that generate only one bit per iteration, takes advantage of the fractional nature of the map, and also reduces bit generation speed, since fewer iterations of the map are required to reach the required bitstream length.As noted above, since the sum in (14) is limited from i=1 to i=n when n<50, for n=0 there are 9 bits generated, for n=1 there are 18 bits, and so on, until n=49 where 450 bits are generated. Thus, to reach a bitstream of length N, there are ⌈N/459+9∑i=049i/459⌉ iterations required.To test the PRBG, a set of100⋅106 bitstreams is generated and tested through the National Institute of Standards and Technology (NIST) statistical test package SP 800-22 [56]. The test suite consists of 15 statistical tests that are used to verify if a bit sequence is indistinguishable from a random sequence. Each test returns a P value and if the value is higher than a significance level, chosen as the default value 0.01 here, the test is successful. For a PRBG to be classified as random, it should pass all 15 tests. This is verified from Table 1.Table 1
NIST test results (x0=0.1, A=0.7, B=0, C=1, μ=0.2, m=50).
No.TestChi-squareP-valueRate1Frequency0.971 699100/1002BlockFrequency0.437 274100/1003CumulativeSums0.213 309100/1004Runs0.494 39296/1005LongestRun0.319 08499/1006Rank0.275 70999/1007FFT0.366 91899/1008NonOverlappingTemplate0.554 420100/1009OverlappingTemplate0.897 76399/10010Universal0.897 763100/10011ApproximateEntropy0.455 93798/10012RandomExcursions0.031 49763/6413RandomExcursionsVariant0.834 30861/6414Serial0.051 94299/10015LinearComplexity0.657 93399/100


### 4.2. Encryption of Electrophysiological Signal

After the appropriate design of the PRBG, the encryption process is presented. The process consists of two rounds, and a different chaotic map of the form (14) is used in each round, each with parameters x0,A1,B1,C1 and y0,A2,B2,C2. In the first round, the source signal is modulated by adding to it the values of the first chaotic map, in order to mask its structure. Then, the modulated signal is transformed into its binary representation, and encrypted using the PRBG of the previous section, generated using the second chaotic map. The encryption is performed using the bitwise XOR operator between the information bitstream and the chaotic bitstream. The complete process is outlined in Algorithm 1.Algorithm 1: Chaotic Encryption of EEG signal.
Input: An EEC signal E of length ℓ. The key values of two chaotic maps of the form (14), x0,A1,B1,C1,μ1 and y0,A2,B2,C2,μ2.Output: An encrypted signal Eenc of the same length.Step 1. Generate ℓ values of the first chaotic map xi, i=1,…,ℓ. Modulate the EEG signal to mask its structure, asℳi=Ei+⌊105xi⌋/103,i=1,…,ℓ.Step 2. Transform the modulated signal ℳ into its binary representation ℳbin of length 18×ℓ.Step 3. Generate a bitstream ℬ of length 16×ℓ using the map yi as a basis for the PRBG of the previous section (15). Since 8 bits are generated per iteration, the map yi needs to be iterated ⌈2.25×ℓ⌉ times.Step 4. Encrypt the information bitstream ℳbin by combining it with the chaotic bistream asEenc−bin=ℳbin⊕ℬ.Step 5. Transform the binary encrypted signal Eenc−bin to decimal form, to obtain the encrypted signal Eenc.The transformation of the modulated signal into its binary representation in Step 2 of Algorithm 1 is performed as follows. The modulated signal takes values in the range−127.999,127.999. In order to sustain a low run time execution we chose to convert the sampled values by mixing a fixed-point-format with a two’s complement representation. So each sample is split into its sign, integer and a decimal parts. The first bit of the binary representation is used to denote the sign. Then, seven bits are used to represent the integer part and ten bits to represent the decimal part. Hence, there are 18 bits overall for each sample, resulting in a bitstream of length 18×ℓ. This transformation is outlined in Figure 13.Figure 13
Binary representation of the modulated samples using 18 bits.For the decryption process, the original EEG signal can be reconstructed at the receiver end by following the reverse procedure, that is, performing an XOR between the encrypted binary signal and the same bitstream from the PRBG, generated using the same parameter values, transforming the result back to its decimal format, and then performing the demodulation process, as described in Algorithm 2.Algorithm 2: Chaotic Decryption of EEG signal.
Input: An encrypted signal Eenc of length ℓ. The key values of two chaotic maps of the form (14), x0,A1,B1,C1 and y0,A2,B2,C2.Output: A decrypted signal E of the same length.Step 1. Transform the encrypted signal Eenc to its binary form, to obtain the encrypted bitstream Eenc−bin of length 18×ℓ.Step 2. Generate a bitstream ℬ of length 18×ℓ using the map yi as a basis for the PRBG of the previous section (15). Since 8 bits are generated per iteration, the map yi needs to be iterated ⌈2.25×ℓ⌉ times.Step 3. Decrypt the bitstream Eenc−bin by combining it with the chaotic bistream as.ℳbin=ℬ⊕Eenc−binStep 4. Transform the modulated binary signal ℳbin into its decimal representation ℳ of length ℓ.Step 5. Generate ℓ values of the chaotic map xi, i=1,…,ℓ. Demodulate the signal ℳ as to obtain the original signal E.Ei=ℳi−⌊105xi⌋/103,i=1,…,ℓNote that the only information that the receiver needs to know in order to decrypt the signal is the parameter values for the two chaotic maps used to modulate the signal and generate the PRBG. These parameters constitute the key values of the encryption design. Since each chaotic map has four parameters and one initial condition, there are overall ten key values,x0,A1,B1,C1,μ1,y0,A2,B2,C2,μ2. Assuming a 16-digit accuracy, an upper bound for the key space is 108⋅16=10128≈10342.6≈2426. This is higher than the bound of 2100 required to resist brute force attacks [57]. Note though, that since the system is not chaotic for all parameter values, the real key space is lower than 2426.


### 4.3. Microcontroller Implementation

The encryption and decryption process was simulated in a microcontroller. For this, an STM32F103 nucleo development board is used, which incorporates a high-performance ARM M3 32-bit RISC core operating at 72 MHz, depicted at Figure14.Figure 14
STM32F103 nucleo board with a 32bit RISC processor.Since an embedded system usually has low ram storage, the modulated signal was saved at the flash memory of the STM32F103RB which embeds an 128 K Byte of flash memory.The simulation results are shown in Figure15, where the original (plaintext), modulated, and encrypted signals are shown. Clearly, there is visually no relation between these signals, something that is verified in the next section, through a series of measures. The signals are plotted in Matlab, where the data are loaded from the txt files resulting from the microcontroller simulation. The numerical simulations for the encryption performed in Matlab yield similar results and are thus omitted.Figure 15
Original EEG signal, modulated signal, and encrypted signal.Note that the computation time for the complete procedure of modulating the signal, transforming it to binary format and encrypting it, takes around 8.5 seconds. To reduce computational load and increase speed, the termsΓm−i+μ/Γm−i+1 where precomputed and saved in the memory before execution begins. Additionally, further experimentation could lead to reduction in the execution time. For example, in the modulation step the fractional map could be implemented using a shorter memory.


### 4.4. Encryption Performance

To test the performance of the encryption algorithm, a series of tests are performed on the original, modulated and encrypted signals. The results are gathered on Table2, for the simulation performed in Matlab R2018b, as well as the microcontroller implementation. The key values of the maps are x0=0.111, A1=0.7, B1=0, C1=1, μ1=0.2 and y0=0.1, A2=0.7, B2=0, C2=1, μ2=0.2.Table 2
Performance measures of the encryption design.
SignalsEECMatlabMicrocontrollerModulatedEncryptedModulatedEncryptedSSIM—−0.11600.001 6−0.08940.006 8LLR—1.995 51.993 71.973 32.045 6SNR—−9.7873−14.6874−9.8402−14.6757rxy—0.307 6−0.01260.303 70.021 0SD—39.494 944.053 639.652 143.864 9ApEn0.497 01.577 42.270 21.579 82.228 2Note that apart from the Approximate Entropy measure that is computed for each individual signal, the rest of the measures are computed in each case between the original and modulated signal, and between the original and encrypted signal, and shown in their respective columns.In all cases, it can be seen that the measures for the Matlab and microcontroller implementations are very close to each other. This means that the microcontroller encryption performs equally well compared to Matlab.


#### 4.4.1. Histogram

Initially, the histogram of the original, modulated and encrypted signals are plotted in Figure16 for the microcontroller simulation. An encrypted signal should have a uniform histogram, so that no information on the distribution of the signal’s values are revealed. Indeed it can be seen that the histogram of the encrypted signal is much more uniform compared to the original signal and modulated ones, that have a shape closer to a normal distribution. The simulation performed in Matlab yields similar histograms, so both simulations perform the same.Figure 16
Histograms of the encrypted, modulated, and encrypted signals.


#### 4.4.2. Structural Similarity Index (SSIM)

The SSIM is a measure of the structural similarity between two signals [29], initially considered for images [58]. It is given by(17)SSIM=2μxμy+S12δxy+S2μx2+μy2+S1δx2+δy2+S2,where μx,μy the mean values of the original and encrypted (or modulated) signals respectively, δx2,δy2 their variances, and δxy their cross-covariance. The parameters S1,S2 take small values, to avoid unstable results when the denominator is close to zero.The SSIM value is between−1,1, with 1 for identical signals, and 0 for signals with no similarity. So, the closest the measure is to zero, the better the encryption. Indeed, the SSIM between the original and modulated signal performed in Matlab is equal to −0.1160, and between the original and encrypted signal is 0.001 6, which is indeed very close to zero.


#### 4.4.3. Log-Likelihood Ratio (LLR)

The LLR gives an estimation on the quality of encryption [29, 59–61]. It is based on the assumption that the segment can be represented by a p-th order all-pole linear predictive coding model(18)xi=∑m=1pamxi−m+Gxui,where xi is the ith signal sample, am, m=1,…,p the coefficients of the all-pole filter, Gx the gain of the filter and ui an appropriate excitation input to the signal. The LLR is then defined as:(19)LLR=logaxRyaxTayRzayT,where ax the vector of Linear Prediction Coefficients (LPCs) 1,a1,a2,…,am of the original signal, az the LPCs of the encrypted (or modulated) signal, and Ry the autocorrelation matrix of the encrypted (or modulated) signal. A higher LLR value indicates a good encryption. The LLR between the original and modulated signal is 1.995 5 and between the original and encrypted signals is 1.993 7. The values are both high, which indicates a good encryption.


#### 4.4.4. Signal to Noise Ratio (SNR)

The SNR is defined as [29, 60]:(20)SNR=10log10∑i=1Nxi2∑i=1Nxi−yi2,where x,y the encrypted and decrypted (or modulated) signals respectively, and N the number of samples. A low SNR indicates a good encryption. The SNR between the original and modulated signal is −9.7873, while the SNR between the original and encrypted signal is −14.6874. Both values are very low and indicate a good encryption.


#### 4.4.5. Correlation Coefficient

The correlation coefficient [29] between the original and encrypted (or modulated) signals is computed as:(21)rxy=cvx,yδx2δy2,where cvx,y the covariance of the two signals and δ2 their variances. For uncorrelated signals, the correlation coefficient should be close to zero. The correlation coefficient between original and modulated signal is 0.307 6, while the coefficient between original and encrypted is −0.0126.


#### 4.4.6. Spectral Distortion (SD)

The SD measures the difference between the spectrum of the original and encrypted (or modulated) signals [29, 60–62]. It is computed as:(22)SD=1M∑i=0M−1Vx,i−Vy,i,where Vx,i,Vy,i the spectrum of the original and encrypted (or modulated) signal in decibel at instance i. A higher value of SD indicates a better encryption. The SD between the original and moculated signals is 39.494 9, and the SD between the original and encrypted signals is 44.053 6, so both values are high.


#### 4.4.7. Approximate Entropy (ApEn)

The ApEn is applied to measure the complexity of the serie of data [51]. A series with a higher ApEn is considered as being more complex, so for the proposed scheme, the encrypted signal should have a higher value of ApEn compared to the original signal. Indeed the ApEn of the original signal is 0.497 0, for the modulated signal is 1.577 4, and for the encrypted signal is 2.270 2, so the encrypted signal has the highest ApEn.Overall, each of the measures computed indicates a good encryption performance. Moreover, the design can be accurately replicated on a microcontroller board.


## 4.1. Design a PRBG

The first step in the encryption scheme is to design a chaotic pseudo-random bit generator. Since the main computational drawback of fractional order systems is their memory effect, the system will be implemented using a finite memory. Thus, its finite form is given as(14)ℳi=εi+⌊105xi⌋103,i=1,…,ℓ,where m denotes the memory used. As one can see from (14), when i=1,…,m the fractional map states depends on the past m variable, i.e, x11,…,x1m. After some trial and error, the memory is chosen as m=50, which satisfyingly models the dynamical behavior of the system, as can be seen in Figure 12. Although the diagram is not identical to the one of the fractional map shown in Figure 4, it accurately showcases chaotic behavior in the range of around 0.2 to 0.8. Implementing fractional systems using finite memory is a standard approach in reducing computational cost and avoiding overflow, and it is used in most microcontroller implementations of fractional chaotic systems [41, 42, 55].Figure 12
Bifurcation diagram of (14) with respect to parameter A, for B=0, C=1, μ=0.2, m=50, and x0=0.1 (blue), x0=−0.1 (red).In order to take advantage of the fractional nature of the proposed map, each individual term inside the finite sum (14) is used in the bit generation. So in each iteration of the map, the PRBG generates bits using the following rules:(15)εenc−bin=ℳbin⊕ℬ,where ⌊⋅⌋ denotes the floor operation. So in each iteration, each individual term of the sum in (14), as well as the value xn of the map are multiplied by 1015 and the integer part of this product is taken modulo 512. The result of this operation is then transformed to binary representation. Since the outcome of the modulo is an integer from 0 to 511, transforming it to binary generates nine bits per iteration. Thus, each of the terms bn,1,…,bn,51 corresponds to nine bits. The resulting bitstream is taken by concatenating the bits as(16)ℬ=bn−1,1,….,bn−1,51,bn,1,….,bn,51,….With the algorithm, there are 459 bits generated in each iteration. So this approach brings an advantage over classic techniques that generate only one bit per iteration, takes advantage of the fractional nature of the map, and also reduces bit generation speed, since fewer iterations of the map are required to reach the required bitstream length.As noted above, since the sum in (14) is limited from i=1 to i=n when n<50, for n=0 there are 9 bits generated, for n=1 there are 18 bits, and so on, until n=49 where 450 bits are generated. Thus, to reach a bitstream of length N, there are ⌈N/459+9∑i=049i/459⌉ iterations required.To test the PRBG, a set of100⋅106 bitstreams is generated and tested through the National Institute of Standards and Technology (NIST) statistical test package SP 800-22 [56]. The test suite consists of 15 statistical tests that are used to verify if a bit sequence is indistinguishable from a random sequence. Each test returns a P value and if the value is higher than a significance level, chosen as the default value 0.01 here, the test is successful. For a PRBG to be classified as random, it should pass all 15 tests. This is verified from Table 1.Table 1
NIST test results (x0=0.1, A=0.7, B=0, C=1, μ=0.2, m=50).
No.TestChi-squareP-valueRate1Frequency0.971 699100/1002BlockFrequency0.437 274100/1003CumulativeSums0.213 309100/1004Runs0.494 39296/1005LongestRun0.319 08499/1006Rank0.275 70999/1007FFT0.366 91899/1008NonOverlappingTemplate0.554 420100/1009OverlappingTemplate0.897 76399/10010Universal0.897 763100/10011ApproximateEntropy0.455 93798/10012RandomExcursions0.031 49763/6413RandomExcursionsVariant0.834 30861/6414Serial0.051 94299/10015LinearComplexity0.657 93399/100


## 4.2. Encryption of Electrophysiological Signal

After the appropriate design of the PRBG, the encryption process is presented. The process consists of two rounds, and a different chaotic map of the form (14) is used in each round, each with parameters x0,A1,B1,C1 and y0,A2,B2,C2. In the first round, the source signal is modulated by adding to it the values of the first chaotic map, in order to mask its structure. Then, the modulated signal is transformed into its binary representation, and encrypted using the PRBG of the previous section, generated using the second chaotic map. The encryption is performed using the bitwise XOR operator between the information bitstream and the chaotic bitstream. The complete process is outlined in Algorithm 1.Algorithm 1: Chaotic Encryption of EEG signal.
Input: An EEC signal E of length ℓ. The key values of two chaotic maps of the form (14), x0,A1,B1,C1,μ1 and y0,A2,B2,C2,μ2.Output: An encrypted signal Eenc of the same length.Step 1. Generate ℓ values of the first chaotic map xi, i=1,…,ℓ. Modulate the EEG signal to mask its structure, asℳi=Ei+⌊105xi⌋/103,i=1,…,ℓ.Step 2. Transform the modulated signal ℳ into its binary representation ℳbin of length 18×ℓ.Step 3. Generate a bitstream ℬ of length 16×ℓ using the map yi as a basis for the PRBG of the previous section (15). Since 8 bits are generated per iteration, the map yi needs to be iterated ⌈2.25×ℓ⌉ times.Step 4. Encrypt the information bitstream ℳbin by combining it with the chaotic bistream asEenc−bin=ℳbin⊕ℬ.Step 5. Transform the binary encrypted signal Eenc−bin to decimal form, to obtain the encrypted signal Eenc.The transformation of the modulated signal into its binary representation in Step 2 of Algorithm 1 is performed as follows. The modulated signal takes values in the range−127.999,127.999. In order to sustain a low run time execution we chose to convert the sampled values by mixing a fixed-point-format with a two’s complement representation. So each sample is split into its sign, integer and a decimal parts. The first bit of the binary representation is used to denote the sign. Then, seven bits are used to represent the integer part and ten bits to represent the decimal part. Hence, there are 18 bits overall for each sample, resulting in a bitstream of length 18×ℓ. This transformation is outlined in Figure 13.Figure 13
Binary representation of the modulated samples using 18 bits.For the decryption process, the original EEG signal can be reconstructed at the receiver end by following the reverse procedure, that is, performing an XOR between the encrypted binary signal and the same bitstream from the PRBG, generated using the same parameter values, transforming the result back to its decimal format, and then performing the demodulation process, as described in Algorithm 2.Algorithm 2: Chaotic Decryption of EEG signal.
Input: An encrypted signal Eenc of length ℓ. The key values of two chaotic maps of the form (14), x0,A1,B1,C1 and y0,A2,B2,C2.Output: A decrypted signal E of the same length.Step 1. Transform the encrypted signal Eenc to its binary form, to obtain the encrypted bitstream Eenc−bin of length 18×ℓ.Step 2. Generate a bitstream ℬ of length 18×ℓ using the map yi as a basis for the PRBG of the previous section (15). Since 8 bits are generated per iteration, the map yi needs to be iterated ⌈2.25×ℓ⌉ times.Step 3. Decrypt the bitstream Eenc−bin by combining it with the chaotic bistream as.ℳbin=ℬ⊕Eenc−binStep 4. Transform the modulated binary signal ℳbin into its decimal representation ℳ of length ℓ.Step 5. Generate ℓ values of the chaotic map xi, i=1,…,ℓ. Demodulate the signal ℳ as to obtain the original signal E.Ei=ℳi−⌊105xi⌋/103,i=1,…,ℓNote that the only information that the receiver needs to know in order to decrypt the signal is the parameter values for the two chaotic maps used to modulate the signal and generate the PRBG. These parameters constitute the key values of the encryption design. Since each chaotic map has four parameters and one initial condition, there are overall ten key values,x0,A1,B1,C1,μ1,y0,A2,B2,C2,μ2. Assuming a 16-digit accuracy, an upper bound for the key space is 108⋅16=10128≈10342.6≈2426. This is higher than the bound of 2100 required to resist brute force attacks [57]. Note though, that since the system is not chaotic for all parameter values, the real key space is lower than 2426.


## 4.3. Microcontroller Implementation

The encryption and decryption process was simulated in a microcontroller. For this, an STM32F103 nucleo development board is used, which incorporates a high-performance ARM M3 32-bit RISC core operating at 72 MHz, depicted at Figure14.Figure 14
STM32F103 nucleo board with a 32bit RISC processor.Since an embedded system usually has low ram storage, the modulated signal was saved at the flash memory of the STM32F103RB which embeds an 128 K Byte of flash memory.The simulation results are shown in Figure15, where the original (plaintext), modulated, and encrypted signals are shown. Clearly, there is visually no relation between these signals, something that is verified in the next section, through a series of measures. The signals are plotted in Matlab, where the data are loaded from the txt files resulting from the microcontroller simulation. The numerical simulations for the encryption performed in Matlab yield similar results and are thus omitted.Figure 15
Original EEG signal, modulated signal, and encrypted signal.Note that the computation time for the complete procedure of modulating the signal, transforming it to binary format and encrypting it, takes around 8.5 seconds. To reduce computational load and increase speed, the termsΓm−i+μ/Γm−i+1 where precomputed and saved in the memory before execution begins. Additionally, further experimentation could lead to reduction in the execution time. For example, in the modulation step the fractional map could be implemented using a shorter memory.


## 4.4. Encryption Performance

To test the performance of the encryption algorithm, a series of tests are performed on the original, modulated and encrypted signals. The results are gathered on Table2, for the simulation performed in Matlab R2018b, as well as the microcontroller implementation. The key values of the maps are x0=0.111, A1=0.7, B1=0, C1=1, μ1=0.2 and y0=0.1, A2=0.7, B2=0, C2=1, μ2=0.2.Table 2
Performance measures of the encryption design.
SignalsEECMatlabMicrocontrollerModulatedEncryptedModulatedEncryptedSSIM—−0.11600.001 6−0.08940.006 8LLR—1.995 51.993 71.973 32.045 6SNR—−9.7873−14.6874−9.8402−14.6757rxy—0.307 6−0.01260.303 70.021 0SD—39.494 944.053 639.652 143.864 9ApEn0.497 01.577 42.270 21.579 82.228 2Note that apart from the Approximate Entropy measure that is computed for each individual signal, the rest of the measures are computed in each case between the original and modulated signal, and between the original and encrypted signal, and shown in their respective columns.In all cases, it can be seen that the measures for the Matlab and microcontroller implementations are very close to each other. This means that the microcontroller encryption performs equally well compared to Matlab.


### 4.4.1. Histogram

Initially, the histogram of the original, modulated and encrypted signals are plotted in Figure16 for the microcontroller simulation. An encrypted signal should have a uniform histogram, so that no information on the distribution of the signal’s values are revealed. Indeed it can be seen that the histogram of the encrypted signal is much more uniform compared to the original signal and modulated ones, that have a shape closer to a normal distribution. The simulation performed in Matlab yields similar histograms, so both simulations perform the same.Figure 16
Histograms of the encrypted, modulated, and encrypted signals.


### 4.4.2. Structural Similarity Index (SSIM)

The SSIM is a measure of the structural similarity between two signals [29], initially considered for images [58]. It is given by(17)SSIM=2μxμy+S12δxy+S2μx2+μy2+S1δx2+δy2+S2,where μx,μy the mean values of the original and encrypted (or modulated) signals respectively, δx2,δy2 their variances, and δxy their cross-covariance. The parameters S1,S2 take small values, to avoid unstable results when the denominator is close to zero.The SSIM value is between−1,1, with 1 for identical signals, and 0 for signals with no similarity. So, the closest the measure is to zero, the better the encryption. Indeed, the SSIM between the original and modulated signal performed in Matlab is equal to −0.1160, and between the original and encrypted signal is 0.001 6, which is indeed very close to zero.


### 4.4.3. Log-Likelihood Ratio (LLR)

The LLR gives an estimation on the quality of encryption [29, 59–61]. It is based on the assumption that the segment can be represented by a p-th order all-pole linear predictive coding model(18)xi=∑m=1pamxi−m+Gxui,where xi is the ith signal sample, am, m=1,…,p the coefficients of the all-pole filter, Gx the gain of the filter and ui an appropriate excitation input to the signal. The LLR is then defined as:(19)LLR=logaxRyaxTayRzayT,where ax the vector of Linear Prediction Coefficients (LPCs) 1,a1,a2,…,am of the original signal, az the LPCs of the encrypted (or modulated) signal, and Ry the autocorrelation matrix of the encrypted (or modulated) signal. A higher LLR value indicates a good encryption. The LLR between the original and modulated signal is 1.995 5 and between the original and encrypted signals is 1.993 7. The values are both high, which indicates a good encryption.


### 4.4.4. Signal to Noise Ratio (SNR)

The SNR is defined as [29, 60]:(20)SNR=10log10∑i=1Nxi2∑i=1Nxi−yi2,where x,y the encrypted and decrypted (or modulated) signals respectively, and N the number of samples. A low SNR indicates a good encryption. The SNR between the original and modulated signal is −9.7873, while the SNR between the original and encrypted signal is −14.6874. Both values are very low and indicate a good encryption.


### 4.4.5. Correlation Coefficient

The correlation coefficient [29] between the original and encrypted (or modulated) signals is computed as:(21)rxy=cvx,yδx2δy2,where cvx,y the covariance of the two signals and δ2 their variances. For uncorrelated signals, the correlation coefficient should be close to zero. The correlation coefficient between original and modulated signal is 0.307 6, while the coefficient between original and encrypted is −0.0126.


### 4.4.6. Spectral Distortion (SD)

The SD measures the difference between the spectrum of the original and encrypted (or modulated) signals [29, 60–62]. It is computed as:(22)SD=1M∑i=0M−1Vx,i−Vy,i,where Vx,i,Vy,i the spectrum of the original and encrypted (or modulated) signal in decibel at instance i. A higher value of SD indicates a better encryption. The SD between the original and moculated signals is 39.494 9, and the SD between the original and encrypted signals is 44.053 6, so both values are high.


### 4.4.7. Approximate Entropy (ApEn)

The ApEn is applied to measure the complexity of the serie of data [51]. A series with a higher ApEn is considered as being more complex, so for the proposed scheme, the encrypted signal should have a higher value of ApEn compared to the original signal. Indeed the ApEn of the original signal is 0.497 0, for the modulated signal is 1.577 4, and for the encrypted signal is 2.270 2, so the encrypted signal has the highest ApEn.Overall, each of the measures computed indicates a good encryption performance. Moreover, the design can be accurately replicated on a microcontroller board.


## 4.4.1. Histogram

Initially, the histogram of the original, modulated and encrypted signals are plotted in Figure16 for the microcontroller simulation. An encrypted signal should have a uniform histogram, so that no information on the distribution of the signal’s values are revealed. Indeed it can be seen that the histogram of the encrypted signal is much more uniform compared to the original signal and modulated ones, that have a shape closer to a normal distribution. The simulation performed in Matlab yields similar histograms, so both simulations perform the same.Figure 16
Histograms of the encrypted, modulated, and encrypted signals.


## 4.4.2. Structural Similarity Index (SSIM)

The SSIM is a measure of the structural similarity between two signals [29], initially considered for images [58]. It is given by(17)SSIM=2μxμy+S12δxy+S2μx2+μy2+S1δx2+δy2+S2,where μx,μy the mean values of the original and encrypted (or modulated) signals respectively, δx2,δy2 their variances, and δxy their cross-covariance. The parameters S1,S2 take small values, to avoid unstable results when the denominator is close to zero.The SSIM value is between−1,1, with 1 for identical signals, and 0 for signals with no similarity. So, the closest the measure is to zero, the better the encryption. Indeed, the SSIM between the original and modulated signal performed in Matlab is equal to −0.1160, and between the original and encrypted signal is 0.001 6, which is indeed very close to zero.


## 4.4.3. Log-Likelihood Ratio (LLR)

The LLR gives an estimation on the quality of encryption [29, 59–61]. It is based on the assumption that the segment can be represented by a p-th order all-pole linear predictive coding model(18)xi=∑m=1pamxi−m+Gxui,where xi is the ith signal sample, am, m=1,…,p the coefficients of the all-pole filter, Gx the gain of the filter and ui an appropriate excitation input to the signal. The LLR is then defined as:(19)LLR=logaxRyaxTayRzayT,where ax the vector of Linear Prediction Coefficients (LPCs) 1,a1,a2,…,am of the original signal, az the LPCs of the encrypted (or modulated) signal, and Ry the autocorrelation matrix of the encrypted (or modulated) signal. A higher LLR value indicates a good encryption. The LLR between the original and modulated signal is 1.995 5 and between the original and encrypted signals is 1.993 7. The values are both high, which indicates a good encryption.


## 4.4.4. Signal to Noise Ratio (SNR)

The SNR is defined as [29, 60]:(20)SNR=10log10∑i=1Nxi2∑i=1Nxi−yi2,where x,y the encrypted and decrypted (or modulated) signals respectively, and N the number of samples. A low SNR indicates a good encryption. The SNR between the original and modulated signal is −9.7873, while the SNR between the original and encrypted signal is −14.6874. Both values are very low and indicate a good encryption.


## 4.4.5. Correlation Coefficient

The correlation coefficient [29] between the original and encrypted (or modulated) signals is computed as:(21)rxy=cvx,yδx2δy2,where cvx,y the covariance of the two signals and δ2 their variances. For uncorrelated signals, the correlation coefficient should be close to zero. The correlation coefficient between original and modulated signal is 0.307 6, while the coefficient between original and encrypted is −0.0126.


## 4.4.6. Spectral Distortion (SD)

The SD measures the difference between the spectrum of the original and encrypted (or modulated) signals [29, 60–62]. It is computed as:(22)SD=1M∑i=0M−1Vx,i−Vy,i,where Vx,i,Vy,i the spectrum of the original and encrypted (or modulated) signal in decibel at instance i. A higher value of SD indicates a better encryption. The SD between the original and moculated signals is 39.494 9, and the SD between the original and encrypted signals is 44.053 6, so both values are high.


## 4.4.7. Approximate Entropy (ApEn)

The ApEn is applied to measure the complexity of the serie of data [51]. A series with a higher ApEn is considered as being more complex, so for the proposed scheme, the encrypted signal should have a higher value of ApEn compared to the original signal. Indeed the ApEn of the original signal is 0.497 0, for the modulated signal is 1.577 4, and for the encrypted signal is 2.270 2, so the encrypted signal has the highest ApEn.Overall, each of the measures computed indicates a good encryption performance. Moreover, the design can be accurately replicated on a microcontroller board.


## 5. Conclusions

Referring to fractional-order discrete-time systems with special features, this paper has introduced the first example of a fractional map with infinite number of equilibria in a bounded domain. The conceived map has also shown coexistence of different types of periodic and chaotic attractors. Dynamics and complexity of the conceived map have been analyzed in details via bifurcation diagrams, maximum Lyapunov exponents, 0-1 test and approximate entropy. Compared with the integer order map, the fractional map has more complexity when the fractional order values is smaller. That is the integer-order map do not have the property of coexisting attractors, while fractional map have it. Afterwards, the proposed map was applied to the problem of encrypting an electrophysiological signal. A PRBG was designed based on the values of the chaotic map, and two rounds of modulation and encryption were performed on the signal. A set of measures were then computed for the original and encrypted signals to verify the performance of the encryption. The design was realised in a microcontroller board, so to increase speed, the fractional map was implemented with finite memory. Future extensions of this work can consider modifications of the encryption design, for example different methods to modulate the signal.Also, in order to improve execution time, the use of a higher frequency CPU combined with hardware accelerators, such as CORDIC (COordinate Rotation DIgital Computer), will allow faster sinusoidal calculations and floating point arithmetic operations.Finally, FPGA based realizations of the encryption design are among the future goals set by the authors, since they have been proven an efficient low cost option for implementing fractional [45] and integer order maps [63].



---

*Source: 3592422-2022-02-24.xml*