A New 4 D Piecewise Linear Multiscroll Chaotic System with Multistability and Its FPGA-Based Implementation

Due to the complex behavior of a multiscroll chaotic system, it is a good candidate for the secure communications. In this paper, by adding an additional variable to the modiﬁed Lorenz-type system, a new chaotic system that includes only linear and piecewise items but can generate 4 n +4 scroll chaotic attractors via choosing the various values of natural number n is proposed. Its dynamics including bifurcation, multistability, and symmetric coexisting attractors, as well as various chaotic and periodic behaviors, are analyzed by means of attraction basin, bifurcation diagram, dynamic map, phase portrait, Lyapunov exponent spectrum, and C 0 complexity in detail. The mechanism of the occurrence for generating multiscroll chaotic attractors is presented. Finally, this multiscroll chaotic system is implemented by using the Altera Cyclone IV EP4CE10F17C8 FPGA. It is found that this FPGA-based design has an advantage of requiring less resources for 0% of the embedded multipliers and 0% of the PLLs of this FPGA are occupied.


Introduction
As indicated in many open literatures, the chaotic system that can generate multishape chaotic attractors has complex dynamical behaviors so that it is difficult to decode its information when it is used in the field of secure communications. Hence, exploring the multishape chaotic system and analyzing its complex dynamical behaviors through theoretical analysis, numerical simulations, implementations, and applications is a hot topic all the time, and many researchers exert much effort to investigate this research area. To date, many kinds of multiscroll chaotic systems have been presented [1][2][3][4][5][6][7]. For example, according to Chua's circuit and using the sine function, a multiscroll chaotic system was introduced and implemented by an electronic circuit, which consists of the commercial trigonometric function chip AD639 and the corresponding auxiliary chips and basic circuit elements [1]. Based on the cellular neural networks and using the trigonometric function, a multiscroll chaotic system was given and analyzed in [2]. By using the hyperbolic tangent function series as the unique nonlinear function, a multiscroll chaotic system was presented and confirmed by an electronic circuit which is constructed by a unity gain voltage buffer, a single current-feedback operational amplifier, and a transconductor in [3]. In [4], the theories, methods, and applications of generating multiscroll chaotic attractors were summarized. Based on the saturated function series, a multiscroll chaotic system that can generate 1D n-scroll, 2D n×m-grid scroll, and 3D n×m×l-grid scroll chaotic attractors was presented and implemented by using an electronic circuit in [8]. By using the Chen chaotic system and the sine function, a multiscroll chaotic and hyperchaotic system was given in [9]. An improved highorder Chua's circuit that can generate multiscroll chaotic attractor by introducing the signum function series was presented and analyzed in [10]. Also, multiwing chaotic systems that can generate multiwing chaotic attractors have also been concerned intensively and many good topologies were presented, such as multiwing Lorenz chaotic system [11], grid multiwing chaotic system [12], 3D grid multiwing chaotic system [13], 2N-butterfly wing chaotic system [14], and fractional-order multiwing chaotic system [15]. All these achievements enrich the chaos theory and the nonlinear circuit and precede their potential application in practical engineering.
In addition, because of the applications of chaotic system in digital secure communication, the realization of chaotic system via digital signal processing circuit attracts people's attention. For example, the multiscroll chaotic systems were implemented in DSP [16], FPGA [17,18], ARM [19,20], and Arduino [21,22]. Generally speaking, as indicated in [23], the implementation based on FPGA has more merits including more flexibility and lower cost. erefore, in this paper, based on the modified Lorenz-type chaotic system that includes only linear and piecewise items, a new 4D piecewise linear chaotic system that can generate 4n + 4 scroll chaotic attractors is proposed and analyzed. By employing attraction basin, bifurcation diagram, dynamic map, phase portrait, Lyapunov exponent spectrum, and C 0 complexity, its dynamical behaviors including multistability, bifurcation, and chaotic and period behaviors are investigated. Also, the mechanism of generating multiscroll chaotic attractors is described theoretically. It is found that symmetrical coexisting attractors will appear in the system with respect to different initial conditions. In addition, the multiscroll chaotic signal generator is implemented by using the digital chip FPGA. Notably, since 0% of the embedded multipliers and 0% of the PLLs of the Altera Cyclone IV EP4CE10F17C8 FPGA are occupied in the realization of the proposed chaotic generator, the new multiscroll chaotic system has the advantage of low FPGA resource required. e paper is organized as follows. In Section 2, the mathematical model of 4D linear piecewise chaotic system is presented. In Section 3, the dynamical behaviors of the proposed chaotic system are analyzed. In Section 4, the mechanism of generating multiscroll chaotic attractors is presented and some numerical simulations are provided. e FPGA implementation of the proposed multiscroll chaotic system is shown in Section 5. Finally, some concluding remarks are presented in Section 6.

A New 4D Piecewise Linear Chaotic System
In [24], based on the traditional Lorenz system, the modified 3D Lorenz-type chaotic system that includes only linear and piecewise items was introduced. Its differential equations are where a and b are positive constant parameters and x, y, and z are variables. Here, by adding an additional variable w into equation (1), the following new 4D system that includes only linear and piecewise items can be obtained.
where sgn () is the signum function, a, b, and c are positive constant parameters, x, y, and z and w are variables, and n is a natural number. Taking n � 0, obviously, there are two modes in system (2). One mode is ẏ � ẏandż �ż + which is denoted by M (−, +), and the other mode is ẏ � ẏ + andż �żwhich is denoted by M(+, −). For f (z) � sgn (4z), its value depends on the sign of z. If z > 0, f (z)�1. If z < 0, f (z)�−1. However, if z � 0, f (z) � 0. It is the same for sgn (4x) and sgn (4y). us, from system (2), one can obtain that this chaotic system has the six following equilibrium points: where X, Y, Z, and W are denoted by the equilibrium points of system (2). e Jacobian matrix of system (2) at the equilibrium points is where δ() is the delta function and it is the differential of sgn (). For δ (4Z), its value equals zero for all Z, except Z � 0 where it is infinite. When Z � 0, the related equilibrium points are S + 0,-0 . Because the value of δ (0) equals infinity, the corresponding eigenvalues of equation ( where "inf" is positive infinite. en, the Routh array can be derived to know the roots' characteristics of equation (6) erefore, based on the Routh-Hurwitz criterion, for the characteristic equation of the Jacobian matrix (5) at S + 0,-0 , it has three roots with negative real part and one root with positive real number, since the number of changes in sign of the first column of the Routh array is one. Hence, these two equilibrium points (S +0,-0 ) are unstable point so that the system's trajectory will be departed rapidly at S + 0,-0 .
However, for S + 1, + 2 and S -1,-2 , even if the Jacobian matrices of system (2) under two modes are different from each other, the characteristic equation of the Jacobian matrix (5) in both two modes is identical and given by us, under these four equilibrium points (S + 1, + 2 and S − 1, −2 ), the sets of eigenvalues are the same. Based on the above parameters, the eigenvalues of the Jacobian matrix (5) of the system under two modes is (−37.7489, 2.8744 + 14.4482i, 2.8744-14.4482i, −16) which means that all these four equilibrium points are index-2 saddle foci. us, the system's trajectory will be in a spiral movement at these four equilibrium points which result in four-scroll attractors. For example, under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01), the time-domain waveforms, the phase portraits, and the Poincaré maps from the numerical simulations by using the Matlab software for system (2) under two modes which are shown in Figures 1 and 2, respectively, can be obtained.

Case 1.
e red color and the blue color in Figure 4(a) are mixed in both positive and negative of y. at is to say, the four-scroll chaotic attractors under two types of initial conditions will appear. For example, let a � 19.2; the four chaotic attractors are shown in Figure 4(b). Figure 4(a), the red color is only filled in the positive of y, while the blue color is only filled in the negative of y. Hence, under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, −0.01, −0.01, 0.01), it has two-scroll chaotic attractors in the positive of y. However, under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01), only twoscroll chaotic attractors will appear in the negative of y. For example, choose a � 40; the numerical simulations are shown in Figure 4(c). Figure 4(a) are filled in the negative of y. Hence, under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, -0.01, −0.01, 0.01) and (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01), only two-scroll chaotic attractors in the negative of y will appear. For example, take a � 48; the two-scroll chaotic attractors are shown in Figure 4(d).

Complexity
Note that the influence of parameters b and c on the dynamical behaviors can also be analyzed by using the above technique. Here, under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01) (blue case) and (x 0 , y 0 , z 0 , w 0 )�(0.01, −0.01, −0.01, 0.01) (red case), only the bifurcation diagram about the variation of parameters b and c is presented and shown in Figure 5.

C 0 Complexity Analysis on System Parameters.
As indicated in [27], C 0 complexity is defined as the proportion of irregular components in the sequence, and C 0 algorithm based on fast Fourier transform removes the regular part of the signal transformation domain and leaves the irregular part. e larger the proportion of the irregular part in the sequence is, the closer the corresponding time series is to the random, and the greater the complexity is. Hence, C 0 algorithm can be used to calculate the randomness of chaotic time series; that is to say, it can be adopted to describe the complexity of chaotic time series [28].
When the length of chaotic time series of the new system is selected as N � 8000 and the complex parameter is selected as r � 10, also C 0 algorithm is used to calculate the complexity of x sequence with respect to varying parameter a, as shown in Figure 6(a). From Figures 6(a) and 3(a), one can obtain that a large C 0 complexity value corresponds to a positive Lyapunov exponent, indicating that the complexity of x sequence is high. Comparing Figure 6(a) with Figure 4(a), the C 0 complexity curve has good consistency with its bifurcation diagram. Similarly, the C 0 complexity of z sequence with respect to varying parameter b is shown in Figure 6(b), and the C 0 complexity of w sequence with respect to varying parameter c is shown in Figure 6(c). By combining Figures 6(b) and 5(a), as well as Figures 6(c) and 5(b), one can see that the C 0 complexity of the proposed system in the chaotic state is greater than that of the proposed system in its other states.
Hence, from C 0 complexity, it can effectively distinguish the dynamic characteristics of continuous chaotic time series

Dynamic Map of the New Chaotic
System. Under a � 112 and initial condition being (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01), dynamic map that is used to describe the dynamics of the proposed chaotic system with respect to the common influence of the varying parameters b and c is described in Figure 8. Obviously, there are three different color regions suggesting three different dynamic behaviors in Figure 8. e cyan area marked with L indicates that when the values of parameters b and c are within this region, the proposed chaotic system will be in a chaotic state, and the maximum

Complexity
Lyapunov exponent corresponding to this area is positive. e blue area marked with H indicates that when the values of parameters b and c are within this region, the proposed chaotic system will be in a quasi-periodic state, and there are two maximum Lyapunov values, both of which are zero corresponding to the blue area. e green area marked with G indicates that the proposed chaotic system will be in a periodic state under these parameters values, while there is only one maximum Lyapunov exponent whose value is zero in the green region. Hence, the dynamics of the proposed chaotic system for various parameters b and c can be vividly described in the dynamics map. en different parameters b and c are selected in three regions of the dynamic map of the proposed chaotic system, and the maximum Lyapunov exponents for different parameters are calculated. e calculated results are summarized in Table 1, where the first point b � 2.4 and c � 16 is chosen in the L region, the second point b � 1.5 and c � 33 is chosen in the H region, and the third point b � 3 and c � 55 is chosen in the G region. According to the calculated results, one can judge that the proposed chaotic system moves in chaotic state, quasi-periodic state, or periodic state when the different parameters are chosen. In order to verify the maximum Lyapunov exponents of the proposed chaotic system under different parameters, we have given the timedomain diagrams and phase diagrams of this system in three different dynamical regions, as shown in Figure 7. Hence, based on the relationship [29] between the dynamical characteristics of the nonlinear system and the values of the Lyapunov exponents, one can see that the dynamic behaviors of the proposed chaotic system are consistent with the results judged by the calculated Lyapunov exponents. Besides, the analysis method described above can also be employed to investigate the dynamics of the proposed chaotic system with respect to other different parameters.

Attraction Basin and Multistability.
e attraction domain is an important tool to analyze coexisting attractors with respect to different initial states of the dynamical system. Here, under a � 112, b � 2.4, c � 16, x 0 � 0.01, and w 0 � 0.01, the attraction basin of the proposed chaotic system is described in Figure 9. One can see that the basin of attraction contains two different color regions, indicating the coexisting characteristic of different state attractors; that is, the proposed chaotic system is of multistability. Among them, the pink region marked with E suggests that if the initial state starts from this region, the proposed chaotic system will be in one chaotic state. However, the yellow region marked with F indicates that if the initial state starts from this region, the proposed chaotic system will be in another chaotic state, which is symmetric to chaotic state in region E. at is to say, the initial trajectories of the proposed chaotic system from the pink region or the yellow region will enter two the symmetric chaotic motion states eventually and respectively.
As shown in Case 5 of Section 3.1, where the initial condition y 0 �0.01 and z 0 �0.01 is chosen in region F, the initial condition y 0 �−0.01 and z 0 �−0.01 is chosen in region E, respectively; they finally display two symmetric coexisting chaotic attractors in Figure 4(f ). Hence, there is multistability in the proposed chaotic system whose motion trajectories starting from different initial states in chaotic region (pink and yellow regions) always enter one of the symmetric chaotic states.

Multiscroll Chaotic Attractors Generation
In this section, the system parameters are designed to generate multiscroll chaotic attractors. From the theoretical analysis of four-scroll chaotic attractors generation in system (2) under n � 0, one can see that system (2) has 6 equilibrium points where two of them are unstable saddle point and the remaining four are index-2 saddle foci. For the unstable point, it cannot form the scroll attractors, since the trajectory of the system will be away from it quickly. So, the four index-2 saddle foci are essential to generate the four-scroll chaotic attractors. us, in order to obtain the multiscroll chaotic attractors, it is necessary e Jacobian matrix of system (2) at the equilibrium points is From equation (10) and the definition of δ (), one can see that it has δ (4Z-0.5) whose value will be infinite under Z�0.125, δ (4Z) whose value will be infinite under Z � 0, and δ (4Z + 0.5) whose value will be infinite under Z�−0.125.
us, under these three cases, the corresponding eigenvalues of equation (10) can also not be calculated.
Of course, the same technique in Section 2 can be used and the same results can be obtained; that is, S + 01, + 02 , S + 0,-0 , and S -01,-02 are unstable point. Hence, the system's trajectory will be departed rapidly at these six equilibrium points. However, for S + 3, + 4 , S + 1, + 2 , S -1,-2 , and S -3,-4 , the corresponding eigenvalues of equation (10) can be obtained and the results are the same as those in Section 2. erefore, all these eight equilibrium points are index-2 saddle foci, which are necessary to form eight-scroll chaotic attractors. For example, based on a � 32, b � 2.4, c � 16, and n � 1 and under the initial condition (x 0 , y 0 , z 0 , w 0 )�(0.01, 0.01, 0.01, 0.01), the numerical simulations by using the ode45 algorithm with variable step size and the simulation time T �1000s with time step being 0.002s in Matlab software are shown in Figure 10. One can see that the proposed multiscroll chaotic system has eight-scroll chaotic attractors. In addition, for the same conditions of simulation algorithm and simulation time, the calculated Lyapunov exponents of the proposed multiscroll chaotic system under this case are LE 1 � 2.76, LE 2 � 0, LE 3 � −11.37, and LE 4 � −39.35.

FPGA Implementation
In this section, the proposed multiscroll chaotic system is implemented by using FPGA. Here, the second-order Runge-Kutta method and the fixed-point data format 32Q26, where 1 bit is the sign part, 5 bits is the integer part, and 26 bits is the decimal part, are applied to calculate the multiscroll chaotic system.

Two-Stage Iteration.
Based on the second-order Runge-Kutta method, the discretization about system (2) under M (−, +) can be derived for implementing it in FPGA and its results are where Δt is the discrete time step. Defining V k �(x k , y k , z k , w k ) T , the proposed multiscroll chaotic system can be expressed as follows: Based on the second-order Runge-Kutta method [30], there are two stage iteration routines. In stage 1, the state vector is updated as follows: where h is the discrete time gap and h � 1/128 is taken for the implementation in experimental results. V k + 1/2 indicates the intermediate results in a half step, that is, Δt � h/2.
In stage 2, V k + 1 is updated based on V k + 1/2 from stage 1, and its formula is as follows:

Verilog HDL Implementation.
Here, Verilog HDL language [31,32] is used to realize the digital chaotic system in FPGA. Notably, in programming, after initial state is defined, state machine method is applied to split the whole digital chaotic generator into several subprocesses, which includes nine states (S 0 , S 1 , S 2 , S 3 , S 4 , S 5 , S 6 , S 7 , and S 8 ) and the corresponding state machine flowchart is shown in Figure 11. We have the following: (i) S 0 : calculation about 4x, 4y, and 4z by using shift operation. (ii) S 1 : implementation of sgn (4x), sgn (4y), |x|, and sgn (4z + 0.5k) operations where k � -n:n. (iii) S 2 : add up all sgn (4z + 0.5k), parallel implementation of ax and ay operations, and parallel implementations to get the results of 16sgn (4x)z and 16sgn (4y)w. (iv) S 3 : parallel implementations to get the results of -cw and f (z). (v) S 4 : parallel implementations to get the results of ẋ, ẏ,ż, and ẇ. (vi) S 5 : select 32 bits to present ẋ, ẏ,ż, and ẇ. (vii) S 6 : according to the second-order Runge-Kutta method expressed in equations (13) and (14), the variables ẋ, ẏ,ż, and ẇ are iterated. Obviously, there are two iterative stages. When the iteration is in stage 1, implement the operation in equation (13) and then direct to state S 0 . When the iteration is in stage 2, implement the operation in equation (14) and then direct to state S 7 .   e top-layer RTL viewer of the realization is shown in Figure 12, where we have the following: "clock" (the clock for FPGA, the input clock is set as 50 MHz in our experiments), "rst1" (the reset for FPGA), "clock1" and "clock2" (the clock for external DAC chip), "wreset1" and "wreset2" (the reset for external DAC chip), "ouput_en" (the output is valid under the high level signal), "x[13..0]" (the digital output of x in each step), "y[13..0]" (the digital output of y in each step), "z[13..0]" (the digital output of z in each step), "w [13..0]" (the digital output of w in each step), "dac_y [13..0]" (the output of y [13..0] to the external DAC chip), and "dac_x [13..0]" (the output of z [13..0] to the external DAC chip).

Experimental Results.
After writing the Verilog HDL program for the multiscroll chaotic system successfully, this program can be downloaded into FPGA via USB Blaster so that it is converted to a hardware list and then be configured into FPGA to form the corresponding digital circuits. e Altera Cyclone IV EP4CE10F17C8 FPGA, which has available internal resources including 10320 logic elements, 423936 RAMs implement bits, 2 PLLs, and 46 hardware multipliers, is adopted. Also the external DAC chip AD9767 (14 bits) with two output channels is applied in the hardware platform. e architecture of experimental hardware platform for the proposed chaotic system is shown in Figure 13. Notably, the proposed chaotic system can generate chaotic attractors containing negative signals; however, the digital chip AD9767 whose conversion rate up to 125 Msps can only process positive signals, so that a positive number should be added for each output of FPGA (+ in Figure 13). In order to keep consistency with the simulation results, the DAC output signal needs to be subtracted by the same positive number; that is, a subtraction unit before being captured by the digital oscilloscope GDS 3254 should be needed. e experimental results based on FPGA technology are shown In our experiments, the input clock is set as 50 MHz and the total thermal power dissipation is obtained as 77.01 mW by employing the PowerPlay Power Analyzer Tool in the software Quartus II. e other hardware resources required to implement the new multiscroll chaotic system based on FPGA chip are listed in Table 2. In this FPGA implementation, 22% of the logic elements, 18% of the logic registers, 24% of the RAMs, 0% of the embedded multipliers, and 0% of the PLLs of the Altera Cyclone IV EP4CE10F17C8 FPGA are occupied, where the resource utilization rates are determined by the ratio of the required number of resources to the total number of resources in the Altera Cyclone IV EP4CE10F17C8 FPGA board. Besides, the throughput of this digital implementation is given as 100Mbits/s, the last two lines provide that there are 16 clock cycles required to evaluate a single iteration, and the delay time is given as 320 nanoseconds when using a 50 MHz clock signal. In addition, the comparisons on FPGA implementations between some existent chaotic systems given in [33,34] and the proposed chaotic system are shown in Table 3. In summary, one can see that the calculation accuracy of the second-order Runge-Kutta method is sufficient to meet the requirements of the proposed chaotic system, and it requires fewer FPGA's resources. Hence, in this paper, the FPGA realization is feasible, the resource consumption is reasonable, the operation is stable and the output speed is fast.

Conclusion
Based on the modified Lorenz-type system and adding an additional variable, a novel chaotic system that can generate multiscroll chaotic attractors is proposed. Attraction basin, bifurcation diagram, dynamic map, phase portrait, Lyapunov exponent spectrum, and C 0 complexity are carried out for this proposed chaotic system to analyze its dynamical behaviors. Research results show that there are some typical dynamics including multistability, bifurcation, chaotic attractors, quasi-periodic attractor, and periodic attractor in this proposed chaotic system. For the different initial conditions, symmetric coexisting attractors will occur in this system. Moreover, theoretical analysis, numerical simulations, and FPGA implementations together demonstrate that the 4n + 4 scroll chaotic attractors can be obtained by selecting the value of the natural number n. Compared with other existent FPGA implementations of the chaotic generators, the proposed multiscroll chaotic system has an advantage of requiring less resources, especially only requiring 0% of the embedded multipliers and 0% of the PLLs of the FPGA chip. Hence, this work will benefit chaotic secure communications, which will be our future work.

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

Conflicts of Interest
e authors declare that they have no conflicts of interest.