Dynamic analysis of a fractional-order hyperchaotic system and its application in image encryption

By adding an additional state $w$ and modifying the structure and parameters of 3D chaotic system [5], we introduce the following 4D quadratic smooth autonomous system:

$$\begin{eqnarray}&&\left\{\begin{array}{l}\dot{x}=-{a}_{1}x+{a}_{2}y\\ \dot{y}={a}_{3}x-{a}_{4}xz+{a}_{5}w\\ \dot{z}={a}_{6}{x}^{2}-{a}_{7}z\\ \dot{w}=-{a}_{8}x\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where ${a}_{i}(i=1,\,2,\cdots ,8)$ is positive real parameters, and ${(x,y,z,w)}^{{\rm{T}}}\in {R}^{4}$ is the state vector. It is worth noting that the introduced state vector $w$ can be considered as a simple external feedback controller which is connected with the system state $x.$ Consequently, from the view of anticontrol of chaos, $w$ can enhance the chaotic dynamical behaviors.

Based on the definition of the Caputo fractional-order differentiation, system (1) can be expressed by

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}_{t}^{q}{x}_{1}=-{a}_{1}{x}_{1}+{a}_{2}{x}_{2}\\ {D}_{t}^{q}{x}_{2}={a}_{3}{x}_{1}-{a}_{4}{x}_{1}{x}_{3}+{a}_{5}{x}_{4}\\ {D}_{t}^{q}{x}_{3}={a}_{6}{x}_{1}^{2}-{a}_{7}{x}_{3}\\ {D}_{t}^{q}{x}_{4}=-{a}_{8}{x}_{1}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where $q$ is the fractional-order, ${x}_{1},{x}_{2},{x}_{3}$ and ${x}_{4}$ represent the state variables, ${a}_{i}(i=1,2,\cdots ,8)$ is the system parameters. It is worthy of note that the system is exactly an integer-order system when $q=1.$ And when $q\gt 1,$ the system is an improper fractional order system. In this paper, we mainly discuss the case of true fractional order, that is, $q\,\lt \,1.$

2.2.1. ADM

ADM is the latest proposed time-domain approximation algorithm suitable for solving fractional-order systems and is capable of handling linear and nonlinear problems excellently [54]. This algorithm can obtain high precision and fast convergent numerical approximate solution. Most importantly, it requires neither large amounts of computer memory nor discretization.

For the system

$$\begin{eqnarray}&&{D}_{{t}_{0}}^{q}x(t)=f(x(t))+g(t)\end{eqnarray} \tag{ 3 }$$

where ${D}_{{t}_{0}}^{q}$ is the Caputo differential operator of order $q,$ $f$ includes the linear and nonlinear parts, $x(t)\,={[{x}_{1}(t),{x}_{2}(t),\mathrm{...},{x}_{n}(t)]}^{{\rm{T}}}$ is the given function variable. And when system is autonomous, $g(t)={[{g}_{1}(t),\mathrm{...},{g}_{n}(t)]}^{T}$ is constant.

The system (3) is divided into three parts as follows

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}_{{t}_{0}}^{q}x(t)=Lx(t)+Nx(t)+g(t)\\ x(k)({t}_{0}^{+})={b}_{k},k=0,1,\mathrm{...},n-1\\ n\in N,n-1\lt q\leqslant n\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where $L$ and $N$ are the linear and nonlinear parts respectively. ${b}_{k}$ represents the initial values.

By integrating both sides of ${D}_{{t}_{0}}^{q}x(t)=Lx(t)+Nx(t)+g(t),$ one can get

$$\begin{eqnarray}&&x={J}_{{t}_{0}}^{q}(Lx+Nx+g)+{\rm{\Phi }}\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Phi }}=\displaystyle {\sum }_{k=0}^{n-1}{b}_{k}\tfrac{{(t-{t}_{0})}^{k}}{k!}$ is the initial condition.

According to the principle of ADM, the solution of system (3) can be shown as

$$\begin{eqnarray}&&x=\displaystyle \sum _{i=0}^{\infty }{x}^{i}\end{eqnarray} \tag{ 6 }$$

Decompose the nonlinear term according to the following formula

$$\begin{eqnarray}&&\left\{\begin{array}{l}{A}_{j}^{i}=\displaystyle \frac{1}{i!}{\left[\displaystyle \frac{{d}^{i}}{d{\lambda }^{i}}N({\nu }_{j}^{i}(\lambda ))\right]}_{\lambda =0}\\ {\nu }_{j}^{i}(\lambda )=\displaystyle \sum _{k=0}^{i}{(\lambda )}^{k}{x}_{j}^{k}\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where $i=0,1,\mathrm{...},$ $j=1,2,\mathrm{...},n.$ Then the nonlinear part is obtained by

$$\begin{eqnarray}&&Nx=\displaystyle \sum _{i=0}^{\infty }{A}^{i}\left({x}^{0},{x}^{1}\mathrm{...},{x}^{i}\right)\end{eqnarray} \tag{ 8 }$$

As a result, the solution of equation (3) is

$$\begin{eqnarray}&&x=\displaystyle \sum _{i=0}^{\infty }{x}^{i}={J}_{{t}_{0}}^{q}\left(L\displaystyle \sum _{i=0}^{\infty }{x}^{i}+\displaystyle \sum _{i=0}^{\infty }{A}^{i}+g\right)+\displaystyle \sum _{k=0}^{n-1}{b}_{k}\displaystyle \frac{{\left(t-{t}_{0}\right)}^{k}}{k!}\end{eqnarray} \tag{ 9 }$$

and the corresponding derivation relation is

$$\begin{eqnarray}&&\left\{\begin{array}{l}{x}^{0}={J}_{{t}_{0}}^{q}g+\displaystyle \sum _{k=0}^{n-1}{b}_{k}\displaystyle \frac{{(t-{t}_{0})}^{k}}{k!}\\ {x}^{1}={J}_{{t}_{0}}^{q}(L{x}^{0}+{A}^{0}({x}^{0}))\\ {x}^{2}={J}_{{t}_{0}}^{q}(L{x}^{1}+{A}^{1}({x}^{0},{x}^{1}))\\ \mathrm{....}\\ {x}^{i}={J}_{{t}_{0}}^{q}(L{x}^{i-1}+{A}^{i-1}({x}^{0},{x}^{1},\mathrm{...},{x}^{i-1}))\\ \mathrm{....}\end{array}\right.\end{eqnarray} \tag{ 10 }$$

2.2.2. Numerical solution

According to ADM, the system (2) can be decomposed as

$$\begin{eqnarray}&&\left[\begin{array}{l}L{x}_{1}\\ L{x}_{2}\\ L{x}_{3}\\ L{x}_{4}\end{array}\right]=\left[\begin{array}{c}-{a}_{1}{x}_{1}+{a}_{2}{x}_{2}\\ {a}_{3}{x}_{1}+{a}_{5}{x}_{4}\\ -{a}_{7}{x}_{3}\\ -{a}_{8}{x}_{1}\end{array}\right],\left[\begin{array}{l}N{x}_{1}\\ N{x}_{2}\\ N{x}_{3}\\ N{x}_{4}\end{array}\right]=\left[\begin{array}{c}0\\ -{a}_{4}{x}_{1}{x}_{3}\\ {a}_{6}{x}_{1}^{2}\\ 0\end{array}\right],\left[\begin{array}{l}{g}_{1}\\ {g}_{2}\\ {g}_{3}\\ {g}_{4}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right].\end{eqnarray} \tag{ 11 }$$

According to equation (7), decomposing the nonlinear terms $-{a}_{4}{x}_{1}{x}_{3}$ and ${a}_{6}{x}_{1}^{2},$ and intercepting the first 6 terms, one can get

$$\begin{eqnarray}&&\left\{\begin{array}{l}{A}_{-{a}_{4}{x}_{1}{x}_{3}}^{0}=-{a}_{4}{x}_{1}^{0}{x}_{3}^{0}\\ {A}_{-{a}_{4}{x}_{1}{x}_{3}}^{1}=-{a}_{4}({x}_{1}^{1}{x}_{3}^{0}+{x}_{1}^{0}{x}_{3}^{1})\\ {A}_{-{a}_{4}{x}_{1}{x}_{3}}^{2}=-{a}_{4}({x}_{1}^{2}{x}_{3}^{0}+{x}_{1}^{1}{x}_{3}^{1}+{x}_{1}^{0}{x}_{3}^{2})\\ {A}_{-{a}_{4}{x}_{1}{x}_{3}}^{3}=-{a}_{4}({x}_{1}^{3}{x}_{3}^{0}+{x}_{1}^{2}{x}_{3}^{1}+{x}_{1}^{1}{x}_{3}^{2}+{x}_{1}^{0}{x}_{3}^{3})\\ {A}_{-{a}_{4}{x}_{1}{x}_{3}}^{4}=-{a}_{4}({x}_{1}^{4}{x}_{3}^{0}+{x}_{1}^{3}{x}_{3}^{1}+{x}_{1}^{2}{x}_{3}^{2}+{x}_{1}^{1}{x}_{3}^{3}+{x}_{1}^{0}{x}_{3}^{4})\\ {A}_{-{a}_{4}{x}_{1}{x}_{3}}^{5}=-{a}_{4}({x}_{1}^{5}{x}_{3}^{0}+{x}_{1}^{4}{x}_{3}^{1}+{x}_{1}^{3}{x}_{3}^{2}+{x}_{1}^{2}{x}_{3}^{3}+{x}_{1}^{1}{x}_{3}^{4}+{x}_{1}^{0}{x}_{3}^{5})\end{array}\right..\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}{A}_{{a}_{6}{x}_{1}^{2}}^{0}={a}_{6}{\left({x}_{1}^{0}\right)}^{2}\\ {A}_{{a}_{6}{x}_{1}^{2}}^{1}={a}_{6}\left(2{x}_{1}^{0}{x}_{1}^{1}\right)\\ {A}_{{a}_{6}{x}_{1}^{2}}^{2}={a}_{6}{\left({\left({x}_{1}^{1}\right)}^{2}+2{x}_{1}^{0}{x}_{1}^{2}\right)}^{2}\\ {A}_{{a}_{6}{x}_{1}^{2}}^{3}={a}_{6}\left(2\right.\left({x}_{1}^{0}{x}_{1}^{3}+{x}_{1}^{1}{x}_{1}^{2}\right)\\ {A}_{{a}_{6}{x}_{1}^{2}}^{4}={a}_{6}\left({\left({x}_{1}^{2}\right)}^{2}+2\left({x}_{1}^{0}{x}_{1}^{4}+{x}_{1}^{1}{x}_{1}^{3}\right)\right)\\ {A}_{{a}_{6}{x}_{1}^{2}}^{5}={a}_{6}\left(2\left({x}_{1}^{0}{x}_{1}^{5}+{x}_{1}^{1}{x}_{1}^{4}+{x}_{1}^{2}{x}_{1}^{3}\right)\right.\end{array}\right..\end{eqnarray} \tag{ 13 }$$

Assume that the initial condition of the system (2) is

$$\begin{eqnarray}&&\left\{\begin{array}{l}{x}_{1}^{0}={x}_{1}({t}_{0})\\ {x}_{2}^{0}={x}_{2}({t}_{0})\\ {x}_{3}^{0}={x}_{3}({t}_{0})\\ {x}_{4}^{0}={x}_{4}({t}_{0})\end{array}\right..\end{eqnarray} \tag{ 14 }$$

Let ${c}_{1}^{0}={x}_{1}^{0},{c}_{2}^{0}={x}_{2}^{0},{c}_{3}^{0}={x}_{3}^{0},{c}_{4}^{0}={x}_{4}^{0}.$ According to ${x}^{1}={J}_{{t}_{0}}^{q}(L{x}^{0}+{A}^{0}({x}^{0}))$ in equation (10) and the properties of fractional calculus, it follows that

$$\begin{eqnarray}&&\left\{\begin{array}{l}{x}_{1}^{1}=(-{a}_{1}{c}_{1}^{0}+{a}_{2}{c}_{2}^{0})\displaystyle \frac{{(t-{t}_{0})}^{q}}{{\rm{\Gamma }}(q+1)}\\ {x}_{2}^{1}=({a}_{3}{c}_{1}^{0}-{a}_{4}{c}_{1}^{0}{c}_{3}^{0}+{a}_{5}{c}_{4}^{0})\displaystyle \frac{{(t-{t}_{0})}^{q}}{{\rm{\Gamma }}(q+1)}\\ {x}_{3}^{1}=({a}_{6}{({c}_{1}^{0})}^{2}-{a}_{7}{c}_{3}^{0})\displaystyle \frac{{(t-{t}_{0})}^{q}}{{\rm{\Gamma }}(q+1)}\\ {x}_{4}^{1}=-{a}_{8}{c}_{1}^{0}\displaystyle \frac{{(t-{t}_{0})}^{q}}{{\rm{\Gamma }}(q+1)}\end{array}\right..\end{eqnarray} \tag{ 15 }$$

The coefficients in equation (15) are denoted as the corresponding variables, as follows

$$\begin{eqnarray}&&\left\{\begin{array}{l}{c}_{1}^{1}=-{a}_{1}{c}_{1}^{0}+{a}_{2}{c}_{2}^{0}\\ {c}_{2}^{1}={a}_{3}{c}_{1}^{0}-{a}_{4}{c}_{1}^{0}{c}_{3}^{0}+{a}_{5}{c}_{4}^{0}\\ {c}_{3}^{1}={a}_{6}{({c}_{1}^{0})}^{2}-{a}_{7}{c}_{3}^{0}\\ {c}_{4}^{1}=-{a}_{8}{c}_{1}^{0}\end{array}\right..\end{eqnarray} \tag{ 16 }$$

Then, ${x}^{1}$ can be expressed as ${x}^{1}={c}^{1}\displaystyle \frac{{(t-{t}_{0})}^{q}}{{\rm{\Gamma }}(q+1)},$ where ${x}^{1}=[{x}_{1}^{1},{x}_{2}^{1},{x}_{3}^{1},{x}_{4}^{1}],$ ${c}^{1}=[{c}_{1}^{1},{c}_{2}^{1},{c}_{3}^{1},{c}_{4}^{1}].$

Similarly, the coefficients of the other four terms of $x$ are

$$\begin{eqnarray}&&\left\{\begin{array}{l}{c}_{1}^{2}=-{a}_{1}{c}_{1}^{1}+{a}_{2}{c}_{2}^{1}\\ {c}_{2}^{2}={a}_{3}{c}_{1}^{1}-{a}_{4}({c}_{1}^{1}{c}_{3}^{0}+{c}_{1}^{0}{c}_{3}^{1})+{a}_{5}{c}_{4}^{1}\\ {c}_{3}^{2}={a}_{6}(2{c}_{1}^{0}{c}_{1}^{1})-{a}_{7}{c}_{3}^{1}\\ {c}_{4}^{2}=-{a}_{8}{c}_{1}^{1}\end{array}\right..\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}{c}_{1}^{3}=-{a}_{1}{c}_{1}^{2}+{a}_{2}{c}_{2}^{2}\\ {c}_{2}^{3}={a}_{3}{c}_{1}^{2}-{a}_{4}\left({c}_{1}^{2}{c}_{3}^{0}+{c}_{1}^{0}{c}_{3}^{2}+{c}_{1}^{1}{c}_{3}^{1}\displaystyle \frac{{\rm{\Gamma }}(2q+1)}{{{\rm{\Gamma }}}^{2}(q+1)}\right)+{a}_{5}{c}_{4}^{1}\\ {c}_{3}^{3}={a}_{6}\left(2{c}_{1}^{0}{c}_{1}^{2}+{({c}_{1}^{1})}^{2}\displaystyle \frac{{\rm{\Gamma }}(2q+1)}{{{\rm{\Gamma }}}^{2}(q+1)}\right)-{a}_{7}{c}_{3}^{2}\\ {c}_{4}^{3}=-{a}_{8}{c}_{1}^{2}\end{array}\right..\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}{c}_{1}^{4}=-{a}_{1}{c}_{1}^{3}+{a}_{2}{c}_{2}^{3}\\ {c}_{2}^{4}={a}_{3}{c}_{1}^{3}-{a}_{4}\left({c}_{1}^{3}{c}_{3}^{0}+{c}_{1}^{0}{c}_{3}^{3}+({c}_{1}^{1}{c}_{3}^{2}+{c}_{1}^{2}{c}_{3}^{1})\displaystyle \frac{{\rm{\Gamma }}(3q+1)}{{\rm{\Gamma }}(q+1){\rm{\Gamma }}(2q+1)}\right)+{a}_{5}{c}_{4}^{3}\\ {c}_{3}^{4}={a}_{6}\left(2{c}_{1}^{0}{c}_{1}^{3}+2{c}_{1}^{1}{c}_{1}^{2}\displaystyle \frac{{\rm{\Gamma }}(3q+1)}{{\rm{\Gamma }}(q+1){\rm{\Gamma }}(2q+1)}\right)-{a}_{7}{c}_{3}^{3}\\ {c}_{4}^{4}=-{a}_{8}{c}_{1}^{3}\end{array}\right..\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}{c}_{1}^{5}=-{a}_{1}{c}_{1}^{4}+{a}_{2}{c}_{2}^{4}\\ {c}_{2}^{5}={a}_{3}{c}_{1}^{4}-{a}_{4}\left({c}_{1}^{4}{c}_{3}^{0}+{c}_{1}^{0}{c}_{3}^{4}+({c}_{1}^{1}{c}_{3}^{3}+{c}_{1}^{3}{c}_{3}^{1})\displaystyle \frac{{\rm{\Gamma }}(4q+1)}{{\rm{\Gamma }}(q+1){\rm{\Gamma }}(3q+1)}+{c}_{1}^{2}{c}_{3}^{2}\displaystyle \frac{{\rm{\Gamma }}(4q+1)}{{{\rm{\Gamma }}}^{2}(2q+1)}\right)+{a}_{5}{c}_{4}^{4}\\ {c}_{3}^{5}={a}_{6}\left(2{c}_{1}^{0}{c}_{1}^{4}+2{c}_{1}^{1}{c}_{1}^{3}\displaystyle \frac{{\rm{\Gamma }}(4q+1)}{{\rm{\Gamma }}(q+1){\rm{\Gamma }}(3q+1)}+{({c}_{1}^{2})}^{2}\displaystyle \frac{{\rm{\Gamma }}(4q+1)}{{{\rm{\Gamma }}}^{2}(2q+1)}\right)-{a}_{7}{c}_{3}^{4}\\ {c}_{4}^{5}=-{a}_{8}{c}_{1}^{4}\end{array}\right..\end{eqnarray} \tag{ 20 }$$

To sum up, the solution of system (2) can be expressed as

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{x}}_{j}(t)={c}_{j}^{0}+{c}_{j}^{1}\displaystyle \frac{{h}^{q}}{{\rm{\Gamma }}(q+1)}+{c}_{j}^{2}\displaystyle \frac{{h}^{2q}}{{\rm{\Gamma }}(2q+1)}+{c}_{j}^{3}\displaystyle \frac{{h}^{3q}}{{\rm{\Gamma }}(3q+1)}\\ \,\,+{c}_{j}^{4}\displaystyle \frac{{h}^{4q}}{{\rm{\Gamma }}(4q+1)}+{c}_{j}^{5}\displaystyle \frac{{h}^{5q}}{{\rm{\Gamma }}(5q+1)}.\end{array}\end{eqnarray} \tag{ 21 }$$

where $j=1,2,3,4,$ $h=t-{t}_{0}$ is the time step.

Dynamic behaviors of fractional-order chaotic systems are usually not only sensitive to system parameters but also closely related to the order. Consequently, the dynamic behaviors of system (2) are explored by using phase diagrams, bifurcation diagrams, Lyapunov exponential spectrum (LEs) and SE complexity varying with different parameters and order.

3.1.1. Phase diagram

In the finite phase space, the larger the space occupied by the system trajectory, the higher the randomness of the system and the better the ergodicity. Take the base parameters value as ${a}_{1}=10,$ ${a}_{2}=15,$ ${a}_{3}=40,$ ${a}_{4}=1,$ ${a}_{5}=1,$ ${a}_{6}=4,$ ${a}_{7}=2.5,$ ${a}_{8}=5,$initial values ${x}_{1}=2.2,$ ${x}_{2}=2.4,$ ${x}_{3}=2.6,$ ${x}_{4}=2.8,$ order $q=0.74$ and time step is $h=0.001.$ Then iterate the system (2) 50000 times and remove the first 60 percent. Some typical attractors of the constructed fractional-order hyperchaotic system with different planes are shown in figure 1. As is shown in the figure, the attractors of the system (2) are distributed over a considerable area, which indicates that the system has superior ergodicity and then sequences with better random characteristics can be obtained.

Zoom In Zoom Out Reset image size

3.1.2. Bifurcation diagram

The bifurcation diagram shows that the motion state of the system will change substantially with different control parameters, so the variation of system state with parameters can be observed intuitively. Let system parameters ${a}_{1}\in [0,25],$ ${a}_{2}\in [0,30],$ ${a}_{3}\in [0,40],$ ${a}_{4}\in [0,30],$ ${a}_{5}\in [0,30],$ ${a}_{6}\in [0,30],$ ${a}_{7}\in [0,12],$ ${a}_{8}\in [0,70],$order $q\in [0.7,1],$ then keep initial condition and time step the same as above. Figure 2 shows the bifurcation diagram with different parameters intervals. From the figure, it can be seen that the system (2) has a rich dynamical behavior and a large chaotic cross section in different parameter ranges.

Zoom In Zoom Out Reset image size

3.1.3. LEs

LEs of a system can effectively characterize the sensitivity of the system to initial values as it evolves in time. A positive Lyapunov exponent indicates that the phase volume of the system is expanding and collapsing in some direction such that the otherwise neighboring tracks in the attractor become increasingly uncorrelated, thus forming a chaotic attractor. In addition, hyperchaotic systems have at least two Lyapunov exponents greater than zero.

Take the base parameter value, then the Lyapunov exponents are obtained through numerical calculation as ${L}_{1}=11.5031,{L}_{2}=0.4258,{L}_{3}$ = $-0.0755,{L}_{4}=-94.2615$ and the Lyapunov dimension is ${d}_{L}=3.1258.$ It is obvious that system (2) is a hyperchaotic system. In addition, figure 3 shows the LEs with different parameters intervals. It is interesting to note that the fourth exponential spectra with different parameters are always less than zero, so they are not shown in the figure. It can be observed from the figure that the system (2) has two positive Lyapunov exponents ${L}_{1},{L}_{2},$ which illustrates that the extension and folding characteristics of its evolution orbit are more complex than conventional chaotic systems. Consequently, it can be inferred that the system is sensitive to initial conditions and the generated hyperchaotic sequences are quite random to be applied in secure communication. Furthermore, it can be observed that the LEs and the bifurcation diagr

Zoom In Zoom Out Reset image size

3.1.4. SE complexity

The SE complexity algorithm adopts Fourier transform, and the corresponding spectral entropy value is obtained through the energy distribution in the Fourier transform domain and Shannon entropy. The specific calculation details are as follows.

Step 1. To make the spectrum reflect the signal energy information more effectively, remove the direct current portion of ${\{x(n)\}}_{0}^{N-1}$ as

$$\begin{eqnarray}&&x(n)=x(n)-\displaystyle \frac{1}{N}\displaystyle \sum _{n=0}^{N-1}x(n)\end{eqnarray} \tag{ 22 }$$

Step 2 . Discrete Fourier transformation of $x(n)$ as

$$\begin{eqnarray}&&X(k)=\displaystyle \sum _{n=0}^{N-1}x(n){e}^{-j\displaystyle \frac{2\pi }{N}nk}=\displaystyle \sum _{n=0}^{N-1}x(n){W}_{N}^{nk}\end{eqnarray} \tag{ 23 }$$

where $k=0,1,\mathrm{...},N-1.$

Step 3. Take the first half of the transformed sequence $X(k)$ for calculation. According to Parseval theorem, the power spectrum value of a certain frequency point is

$$\begin{eqnarray}&&p(k)=\displaystyle \frac{1}{N}{\left|X(k)\right|}^{2}\end{eqnarray} \tag{ 24 }$$

where $k=0,1,\mathrm{...},N/2-1$ and the total power can be defined as

$$\begin{eqnarray}&&{p}_{tot}=\displaystyle \frac{1}{N}\displaystyle \sum _{k=0}^{N/2-1}{\left|X(k)\right|}^{2}\end{eqnarray} \tag{ 25 }$$

The relative power spectrum probability is

$$\begin{eqnarray}&&{P}_{k}=\displaystyle \frac{p(k)}{{p}_{tot}}=\displaystyle \frac{\displaystyle \frac{1}{N}{\left|X(k)\right|}^{2}}{\displaystyle \frac{1}{N}\displaystyle \sum _{k=0}^{N/2-1}{\left|X(k)\right|}^{2}}=\displaystyle \frac{{\left|X(k)\right|}^{2}}{\displaystyle \sum _{k=0}^{N/2-1}{\left|X(k)\right|}^{2}}\end{eqnarray} \tag{ 26 }$$

where $\displaystyle \sum _{k=0}^{N/2-1}{P}_{k}=1.$

Step 4 . Combing ${P}_{k}$ and the concept of Shannon entropy, the spectral entropy $se$ of the signal is obtained as

$$\begin{eqnarray}&&se=-\displaystyle \sum _{k=0}^{N/2-1}{P}_{k}\,\mathrm{ln}\,{P}_{k}\end{eqnarray} \tag{ 27 }$$

If $SE,$ then $SE(N)=\displaystyle \frac{se}{\mathrm{ln}(N/2)}.$ It can be shown that the magnitude of spectral entropy converges to $\mathrm{ln}(N/2).$ To facilitate comparative analysis, the normalized $SE$ is shown below

$$\begin{eqnarray}&&SE(N)=\displaystyle \frac{se}{\mathrm{ln}(N/2)}\end{eqnarray} \tag{ 28 }$$

Take the base parameters and initial values, and then the complexity of chaotic sequence is obtained through numerical simulation, as shown in figure 4. For two-dimensional complexity diagram, the lighter the image color, the simpler the system. On the contrary, the darker the image color, the more complex the system. It is obvious from the figure that the constructed fraction-order hyperchaotic system (2) has rich dynamical behaviors and the chaotic state occupies a large interval window, which indicates that generated sequences have excellent randomness to be applied to chaotic cryptography. In addition, it is important to note that the values of system parameters and order should avoid light color region to improve security.

Zoom In Zoom Out Reset image size

The hardware implementation of chaotic systems is of great significance in practical application. Compared with analog devices, DSP has been widely used in practical applications because of its unique stability, repeatability, large-scale integration, especially programmability and easy implementation.

The DSP chip can be either fixed-point or floating-point according to the system requirements. Its main frequency affects the speed of the operation, which in turn affects the speed of chaotic sequence generation. The width of the data bus has an impact on the numerical accuracy of the generated chaotic sequence, and the peripheral interface of the DSP chip needs to be considered, which is best to facilitate interfacing with other parts of the system. This section uses the typical TMS320F28335 DSP chip, the 16-bit dual channel DAC8552 chip as the D/A converter and the MAX3232 chip as the interface. Among them, TMS320F28335 is a floating-point 32-bit chip with a main frequency up to 150MHz and rich peripherals. The serial interface UART transmits the generated chaotic sequence to the computer accurately and for communication with DSP in real time. The D/A replacement part converts the chaotic sequence generated by DSP calculation into analog signal, and then the attractor phase diagram can be observed on the oscilloscope. Therefore, the D/A conversion part is required to have a dual-channel conversion function with synchronous output, and also has excellent conversion speed and conversion accuracy. In addition, figure 5 is a software flow chart of the system (2) based on DSP.

Zoom In Zoom Out Reset image size

Furthermore, data processing is the key to DSP implementation. In order to realize the fractional hyperchaotic system better, some data processing should be carried out and figure 6 shows the flow chart. Some specific steps are as follows. Because some of the generated hyperchaotic sequences are negative, the binary sequence calculated by DSP will be signed. In order to correct D/A conversion, a sufficiently large positive number needs to be added to all output chaotic sequence values, which is equivalent to translating the system output. The DAC8552 acts as a 16-bit D/A converter, the converted digital quantity must be a positive integer between 0 and 65535. Most of the values generated by DSP are decimal, which need to be multiplied by a positive integer to enlarge the sequence value in equal proportion, so as to make the digital quantity distributed in the whole installation and replacement range.

Zoom In Zoom Out Reset image size

Take the base parameters and initial values, then the DSP hardware implementation platform and some typical hyperchaotic attractors of the system (2) implemented on DSP platform are shown in figures 7 and 8. It is clear from figure 8 that the hyperchaotic attractor realized based on DSP platform and the numerical simulation results are in a consistent manner, which provides certain reference values for the system to be applied in the field of secure communication.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1.1. DNA encoding and decoding rules

4.1.2. DNA addition and subtraction rules

4.1.3. DNA complementary rule

For each nucleotide ${x}_{i},$ the DNA complementary rule [56] is subject to the following equation

$$\begin{eqnarray}&&\left\{\begin{array}{l}{x}_{i}\ne L({x}_{i})\ne L(L({x}_{i}))\ne L(L(L({x}_{i})))\\ {x}_{i}=L(L(L(L({x}_{i}))))\end{array}\right.\end{eqnarray} \tag{ 29 }$$

where ${x}_{i}$ and $L({x}_{i})$ are a pair of complementary nucleic acid bases, following single mapping. Six complementary base pairs can be obtained according to equation (29), as follows:

• (1)  
  ${L}_{1}(A)=T,{L}_{1}(T)=C,{L}_{1}(C)=G,{L}_{1}(G)=A$
• (2)  
  ${L}_{2}(A)=T,{L}_{2}(T)=G,{L}_{2}(G)=C,{L}_{2}(C)=A$
• (3)  
  ${L}_{3}(A)=C,{L}_{3}(C)=T,{L}_{3}(T)=G,{L}_{3}(G)=A$
• (4)  
  ${L}_{4}(A)=C,{L}_{4}(C)=G,{L}_{4}(G)=T,{L}_{4}(T)=A$
• (5)  
  ${L}_{5}(A)=G,{L}_{5}(G)=T,{L}_{5}(T)=C,{L}_{5}(C)=A$
• (6)  
  ${L}_{6}(A)=G,{L}_{6}(G)=C,{L}_{6}(C)=T,{L}_{6}(T)=A$

where ${L}_{i},$ $i=1,2,\mathrm{...},6$ is the ith complement rule. It is worth mentioning that after basic DNA coding, a complementary principle is randomly selected to carry out diffusion of pixel values deeply in this paper.

4.1.4. The improved Arnold algorithm

Expand the plaintext image $P$ into one-dimensional row vector $A,$with size of $MN.$Then the Arnold transformation is performed on any coordinate position $(1,j)$ of vector $A$ to obtain the new coordinate position $(p,q),$ as shown in the following equation:

$$\begin{eqnarray}&&\left[\begin{array}{l}p\\ q\end{array}\right]=\left[\begin{array}{cc}1 & a\\ b & ab+1\end{array}\right]\left[\begin{array}{l}1\\ j\end{array}\right]\end{eqnarray} \tag{ 30 }$$

Then,

$$\begin{eqnarray}&&p=1+aj\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&q=b+(ab+a)j\end{eqnarray} \tag{ 32 }$$

The improved Arnold algorithm [19] only considers the above equation (32) and can realize the transformation between pixel points $(1,j)$ and $(1,q)$ through $a$ and $b$ without considering the role of $p.$ Furthermore, $ab+1$ is regarded as a new random number and is still denoted as $a.$ Finally, equation (32) becomes as shown in equation (33):

$$\begin{eqnarray}&&q=b+aj\end{eqnarray} \tag{ 33 }$$

In this section, a color image encryption algorithm based on the novel fractional-order hyperchaotic system and DNA sequences operations is designed. Figure 9 shows the main process of the designed encryption scheme which includes two scrambling processes and DNA diffusion parts. Some specific details are described as follows.

Step 1. Input an original color image $P$ with the size of $3\times M\times N,$ and decompose it into red, green and blue parts to obtain R, G and B matrices with the size of $M\times N.$

Step 2. Set the parameters and initial values of the constructed system, and iterate it $t1+2MN$ times. Then the former $t1$ values are removed to avoid the transient effect and the hyperchaotic sequence ${\{{x}_{i}\}}_{i=1}^{2MN}$ is obtained.

Step 3. Transform all elements of the sequence ${\{{x}_{i}\}}_{i=1}^{2MN}$ through equation (34) to obtain the pseudo-random sequence $X.$

$$\begin{eqnarray}&&X(i)=mod(floor((abs({x}_{i}))\times {10}^{16}),256)+1\end{eqnarray} \tag{ 34 }$$

where $X(i)$ represents the ith element of chaotic sequence $X,$ $i\in \{1,2,\mathrm{...},2MN\}.$ $abs(x)$ means to calculate the absolute value of $x,$ and $floor(x)$ gives the maximum integer less than or equal to $x.$

Step 4. The pseudo-random number vectors $a,$ $b$ and $q$ are expressed by

$$\begin{eqnarray}&&a=X\left(1:{\rm{M}}N\right),b=X\left(MN+1:2MN\right)\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&q({\rm{i}})=\,\mathrm{mod}\left(b\left(j\right)+a\left(i\right)\times i,MN\right)+1\end{eqnarray} \tag{ 36 }$$

where $a(i),b(j)$ and $q(i)$ are the ith and jth element of vectors $a,b$ and $q$ respectively, $i\in \{1,2,\mathrm{...},MN\},$ $j\in \{MN+1,\mathrm{...},2MN\}.$

Step 5. R, G and B are transformed into one-dimensional row vectors, and then the element at position $(1,j)$ is exchanged with that at position $(1,q(j))$ according to the improved Arnold algorithm, where $j\in \{1,2,\mathrm{...},MN\}.$

Step 6. Restore the vectors R, G and B to the matrix of size $M\times N$ to obtain the result matrices R1, G1 and B1 of the first scrambling.

Step 7. Input the parameters and initial values again, iterate the system (2) $t2+MN/4$ times, and remove the former $t2$ values. Then the hyperchaotic sequences ${\{{x}_{i}\}}_{i=1}^{MN/4},$ ${\{{y}_{i}\}}_{i=1}^{MN/4},$ ${\{{z}_{i}\}}_{i=1}^{MN/4},$ ${\{{w}_{i}\}}_{i=1}^{MN/4}$ are obtained and they are recombined into a new chaotic sequence $S$ with length $MN$ by the following formula.

$$\begin{eqnarray}&&{\rm{S}}\left(1:4:{\rm{end}}\right)=\left\{{{\rm{x}}}_{i}| i=1,\mathrm{...},{\rm{MN}}/4\right\}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\rm{S}}\left(2:4:\mathrm{end}\right)=\left\{{{\rm{y}}}_{i}| i=1,\mathrm{...},{\rm{MN}}/4\right\}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\rm{S}}\left(3:4:{\rm{end}}\right)=\left\{{{\rm{z}}}_{i}| i=1,\mathrm{...},{\rm{MN}}/4\right\}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{\rm{S}}\left(4:4:{\rm{end}}\right)=\left\{{{\rm{w}}}_{i}| i=1,\mathrm{...},{\rm{MN}}/4\right\}\end{eqnarray} \tag{ 40 }$$

Step 8. The chaotic sequence $S$ is quantized by

$$\begin{eqnarray}&&S(i)=mod(floor((abs(S(i)))\times {10}^{16}),256)\end{eqnarray} \tag{ 41 }$$

where $S(i),$ $i\in \{1,2,\mathrm{...},MN\}$ is the ith element of chaotic sequence $S.$

Step 9. Restore the vector $S$ to the matrix with size of $M\times N,$ then the matrix $I$ is obtained by arranging each column of matrix $S$ in descending order. The scrambling matrices R1, G1 and B1 obtained in step 5 are scrambled for the second time by the following formula to obtain matrices R2, G2 and B2.

$$\begin{eqnarray}&&\left\{\begin{array}{l}R2(I(i,j-i),j-i)=R1(I(i,j),j),j\gt i\\ R2(I(i,N+j-i),N+j-i)=R1(I(i,j),j),j\leqslant i\end{array}\right.\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}G2(I(i,j-i),j-i)=G1(I(i,j),j),j\gt i\\ G2(I(i,N+j-i),N+j-i)=G1(I(i,j),j),j\leqslant i\end{array}\right.\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l}B2(I(i,j-i),j-i)=B1(I(i,j),j),j\gt i\\ B2(I(i,N+j-i),N+j-i)=B1(I(i,j),j),j\leqslant i\end{array}\right.\end{eqnarray} \tag{ 44 }$$

Step 10. The scrambling matrices R2, G2 and B2 with size of $M\times N$ are transformed into binary matrices DR2, DG2 and DB2 with size of $8\times M\times N.$ Then, the binary matrices are transformed into DNA matrices DNA_R2, DNA_G2 and DNA_B2 through DNA coding rule r1 and the size is $4\times M\times N.$

Step 11. According to the DNA complementary rule, the sequence $C\_R2=\left\{{c}_{i},i=1,2,\mathrm{...},4MN\right\}$ is obtained by

Zoom In Zoom Out Reset image size

If ${c}_{2i-2}=A,$ then ${c}_{2i-1}={L}_{l1}({d}_{2i-1}).$

If ${c}_{2i-2}=C,$ then ${c}_{2i-1}={L}_{l2}({d}_{2i-1}).$

If ${c}_{2i-2}=G,$ then ${c}_{2i-1}={L}_{l3}({d}_{2i-1}).$

If ${c}_{2i-2}=T,$ then ${c}_{2i-1}={L}_{l4}({d}_{2i-1}).$

If ${c}_{2i-1}=A,$ then ${c}_{2i}={L}_{l5}({d}_{2i}).$

If ${c}_{2i-1}=C,$ then ${c}_{2i}={L}_{l6}({d}_{2i}).$

If ${c}_{2i-1}=G,$ then ${c}_{2i}={L}_{l7}({d}_{2i}).$

If ${c}_{2i-1}=T,$ then ${c}_{2i}={L}_{l8}({d}_{2i}).$

where ${c}_{i}$ is the ith element of the sequence $C\_R2,$ $i\in \{1,\mathrm{...},4MN\}.$ ${d}_{i}$ is the ith element after the DNA matrix DNA_R2 is expanded into a one-dimensional row vector and ${l}_{i},$ $i\in \{1,2,\mathrm{...},8\}$ is any DNA complementarity principle.

Step 12. As in step 10, the other two sequences $C\_B2$ and $C\_G2$ are obtained.

Step 13. Input the parameters and initial values again, iterate the system (2) $t3+MN$ times, and remove the former $t3$ values. Then the hyperchaotic sequences ${\{{x}_{i}\}}_{i=1}^{MN},$ ${\{{y}_{i}\}}_{i=1}^{MN}$ and ${\{{z}_{i}\}}_{i=1}^{MN}$ are obtained.

Step 14. Transform all elements of the sequences in step 12 through the following equations to obtain the pseudo-random sequences $k1,$ $k2$ and $k3.$

$$\begin{eqnarray}&&k1(i)=mod(floor((abs({x}_{i}))\times {10}^{16}),256)\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&k2(i)=mod(floor((abs({y}_{i}))\times {10}^{16}),256)\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&k3(i)=mod(floor((abs({z}_{i}))\times {10}^{16}),256)\end{eqnarray} \tag{ 47 }$$

where $k1(i),$ $k2(i)$ and $k3(i),$ $i\in \{1,2,\mathrm{...},MN\}$ are the ith element of chaotic sequences $k1,$ $k2$ and $k3.$

Step 15. Transform the pseudo-random sequences $k1,$ $k2$ and $k3$ into binary sequences $K1,$ $K2$ and $K3$ with size of $8\times M\times N$ and then transform them into DNA sequences $DNA\_K1,$ $DNA\_K2$ and $DNA\_K3$ with size of $4\times M\times N$ by DNA coding rule r2.

Step 16. The DNA sequences DNA_R3, DNA_G3 and DNA_B3 are obtained by the following formula.

$$\begin{eqnarray}&&DNA\_R3(i)=C\_R2(i)+DNA\_K1(i)+DNA\_R3(i-1)\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&DNA\_G3(i)=C\_G2(i)+DNA\_K2(i)+DNA\_G3(i-1)\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&DNA\_B3(i)=C\_B2(i)+DNA\_K3(i)+DNA\_B3(i-1)\end{eqnarray} \tag{ 50 }$$

where $DNA\_R3(0)=C\_R2(4MN),$ $DNA\_G3(0)=C\_G2(4MN),$ $DNA\_B3(0)=C\_B2(4MN)$ and '+' represents DNA addition.

Step 17. Finally, the DNA sequences DNA_R3, DNA_G3 and DNA_B3 are restored to matrices, which are converted to binary matrices through the DNA decoding rule r3. And then restored to decimal matrices $C1,$ $C2$ and $C3$ which are combined to obtain the final encryption image matrix $C$ with size of $3\times M\times N.$

It is worth noting that, a symmetric encryption scheme is used in this paper, so the decryption algorithm is the inverse process of the encryption. That is to say, inverse diffusion and inverse dislocation are performed on the ciphertext image in turn.

In this paper, MATLAB R2018b is used to verify the designed encryption algorithm in a personal computer with memory 8.00GB, processor AMD A9-9410, CPU 2.90 GHz and the operating system is Windows 7. In addition, color images Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07 with size of $256\times 256$ are used for simulation experiments. Set the key of the proposed algorithm, in which the parameters are ${a}_{1}=10,$ ${a}_{2}=15,$ ${a}_{3}=40,$ ${a}_{4}=1,$ ${a}_{5}=1,$ ${a}_{6}=4,$ ${a}_{7}=2.5,$ ${a}_{8}=5,$ initial values ${x}_{1}=2.2,$ ${x}_{2}=2.4,$ ${x}_{3}=2.6,$ ${x}_{4}=2.8,$ order $q=0.74$ and iteration numbers are $t1=800,$ $t2=500$ and $t3=700.$ Furthermore, the initial nucleic acid base is selected as ${c}_{0}=A$ and DNA encoding rules are $r1=1,$ $r2=1$ $r3=3.$

After numerical simulation, figure 10(a) shows the plaintext images of Sailboat, Fruits, and Pepper, figures 10(b) and (c) are the corresponding encrypted and decrypted results respectively. It is clear from the figure that the encrypted images are all noise-like images, hiding the effective information of the original plaintext images perfectly, which indicates the effectiveness of encryption. In addition, the decrypted images are all no different from the original plaintext images, which illustrates the effectiveness and feasibility of the decryption algorithm.

Zoom In Zoom Out Reset image size

Key space refers to the total number of various keys that can be used in the encryption process. Generally, the larger the key space, the better the security. For the designed encryption algorithm, the key mainly includes: system parameters (${a}_{i},$ $i=1,2,\mathrm{...},8$), initial values(${x}_{i},$ $i=1,2,\mathrm{...},4$), order $q,$ the initial nucleic acid base(${c}_{0}\in A,T,C,G$), DNA encoding rule($r1,r2,r3\in [1,8]$), and DNA complementary rule(${l}_{i}\in [1,6],i=1,2,\mathrm{...},8$). The maximum error range of each parameter and initial value in the decryption process is obtained by numerical simulation, where $q$ is ${10}^{-16},$ ${a}_{i}$ ($i=1,2,4,\mathrm{...},8$), ${x}_{1},{x}_{2}$ are ${10}^{-15}$ and ${a}_{3},$ ${x}_{3},$ ${x}_{4}$ are ${10}^{-14}.$ Therefore, the key space is $4\times {8}^{3}\times {6}^{8}\times {10}^{16+15\times 9+14\times 3}\approx {2}^{672},$ which is much larger than the theoretical key space value ${2}^{100}$ [58] of the cryptosystem. In addition, the key space size of the designed scheme is compared with that in the existing literature [22, 23, 40], and the results are shown in table 3. It is clear that the encryption algorithm designed in this paper has a larger key space and is resistant to exhaustive attacks effectively.

Key sensitivity is a prominent index to detect whether encryption algorithms can resist violent attacks. For an excellent encryption scheme, it should be sensitive enough to the key. That is to say, when the key is changed slightly, the decrypted image will change greatly. In this section, Sailboat in figure 10(a) is selected for testing. With other keys unchanged, a small difference is added to the initial values ${x}_{1},{x}_{2},{x}_{3},{x}_{4}$ and parameters $q,{a}_{1},{a}_{8}.$ The corresponding decryption images are obtained as shown in figure 11. One can see that the decrypted image obtained after slightly changing the key is completely different from that obtained by the correct key, which indicates that the designed algorithm is extremely sensitive to the key so that it can resist violent attacks effectively.

Zoom In Zoom Out Reset image size

The histogram reflects distributions of image pixel values, which is mainly used to evaluate the ability of an algorithm to resist statistical attacks. In the encryption process, the closer the digital transmission distance of the encrypted data, the better the encryption performance. Through numerical simulations, histograms of plaintext Sailboat, Fruits and Pepper images in R, G and B channels are obtained, as shown in figure 12(a). In addition, figure 12(b) shows the corresponding histograms of encrypted images in each channel. It is clear that those histograms of the encrypted images are unified and do not expose any valuable data. Consequently, it is quite difficult for an attacker to obtain any statistical information.

Zoom In Zoom Out Reset image size

In addition, chi-square statistics (unilateral hypothesis test) are commonly used to measure the difference between plaintext and ciphertext image histograms numerically. For an image with gray level of 256 and size of $M\times N,$ the chi-square statistics are calculated by

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=0}^{255}\displaystyle \frac{{({f}_{i}-g)}^{2}}{g}\end{eqnarray} \tag{ 51 }$$

where ${f}_{i}$ represents the frequency of each gray value and $g=MN/256$ is the corresponding theoretical frequency distribution. When the significance level is 0.05, ${\chi }_{0.05}^{2}\left(255\right)=293.24783.$ After calculation, the chi-square statistics of color Sailboat, Fruits, Pepper plaintext images and the corresponding ciphertext images in each channel are shown in table 4. One can see that the calculated values of the chi-square statistics of plaintext images are significantly greater than ${\chi }_{0.05}^{2}\left(255\right),$ and that of the corresponding ciphertext images in each channel are all less than ${\chi }_{0.05}^{2}\left(255\right).$ Therefore, it can be considered that ciphertext image histograms are approximately evenly distributed, indicating that the designed algorithm can withstand statistical attacks effectively.

In general, the correlation between neighboring pixels in the original image is extremely high, so an excellent encryption scheme should minimize it.

In this section, $n$ pairs of adjacent pixels are randomly selected from the R, G and B channels of color images Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07. The corresponding gray values are denoted as $({u}_{i},{v}_{i}),i=1,2,\mathrm{...},MN,$ then the correlation coefficient between $u=\{{u}_{i}\}$ and $v=\{{v}_{i}\}$ is

$$\begin{eqnarray}&&{r}_{uv}=\displaystyle \frac{E((u-E(u))(v-E(v)))}{\sqrt{D(u)}\sqrt{D(v)}}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&D(u)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{({u}_{i}-E(u))}^{2}\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&E(u)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{u}_{i}\end{eqnarray} \tag{ 54 }$$

Let the coordinate of ${u}_{i}$ be $({x}_{i},{y}_{i}).$ If the coordinate of ${v}_{i}$ is $({x}_{i}+1,{y}_{i}),$ then the correlation coefficient in the horizontal direction is calculated. Similarly, if the coordinate is $({x}_{i},{y}_{i}+1)$ or $({x}_{i}+1,{y}_{i}+1),$ then the correlation coefficient in the horizontal or diagonal direction is calculated respectively.

Based on the above equations, take $n=10000$ to obtain the correlation coefficient among R, G and B channels of Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07 plaintext and ciphertext images in each direction. Figure 13 shows the pixel distributions of color Lena image and one can see clearly that the pixels of the original image are concentrated near $y=x,$ while that of the ciphertext image are evenly distributed in the whole pixel value range. That is to say, the proposed algorithm can destroy the correlation of original images effectively and the statistical information cannot be obtained through the pixel distribution. Table 5 lists the correlation coefficients of plaintext and the corresponding ciphertext images in different directions. In addition, the comparison results between the correlation coefficients of Lena ciphertext image in this paper and the existing literature [33, 38, 39] are shown in table 6. Apparently, the correlation coefficient between adjacent pixels of ciphertext image calculated by our scheme is closer to 0, which indicates that the designed scheme can weaken the correlation better and resist statistical attacks effectively.

Zoom In Zoom Out Reset image size

Tables 7 and 8 lists the correlation coefficients of Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07 at the identical and adjacent positions respectively. Moreover, the comparison results under different circumstances of Lena ciphertext image in this paper and the existing literature [26, 39, 44] are shown in table 9. It is clear that the correlation coefficients between each channel of the ciphertext image at the same and adjacent positions are all closer to 0, indicating that the designed scheme has a strong ability to destroy the correlation between different channels.

Information entropy is a prominent feature to evaluate the security performance of encryption schemes, which can reflect the uncertainty of image information. The higher the value of entropy, the less visual information obtained and the better the encryption effect. And it is calculated as shown below.

$$\begin{eqnarray}&&H=-\displaystyle \sum _{i=0}^{L}p(i){\mathrm{log}}_{2}p(i)\end{eqnarray} \tag{ 55 }$$

where $L$ represents the number of gray levels, and $p(i)$ is the probability of occurrence of $i.$ For an image with a gray level of 256, the theoretical value of $H$ is 8.

Table 10 lists the information entropy of the plaintext and corresponding ciphertext Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07 images in each channel. One can see that the information entropy of ciphertext images in each channel is all very close to 8, while that of each plaintext image is far from 8. It can be inferred that the designed encryption algorithm seriously interferes with the pixel value of the original image, and the encryption effect is quite prominent. In addition, the comparison results between the information entropy of Lena image in this paper and the existing literature [33, 40, 51] are shown in table 11. It is clear that the information entropy of ciphertext images obtained by our scheme is closer to 8, which indicates that our algorithm has better security and can withstand information entropy attack effectively.

Differential attack analysis refers to encrypting two slightly different plaintext images with the same key and comparing the differences between two ciphertext images. The number of pixels change rate (NPCR) and the unified average changing intensity (UACI) are two important indicators to test whether a cryptosystem can withstand the differential attack. The corresponding calculation formula is as follows:

$$\begin{eqnarray}&&\left\{\begin{array}{l}NPCR=\displaystyle {\sum }_{i,j}\displaystyle \frac{D(i,j)}{M\times N}\times 100 \% \\ UACI=\displaystyle {\sum }_{i,j}\displaystyle \frac{\left|{C}_{1}(i,j)-{C}_{2}(i,j)\right|}{255\times M\times N}\times 100 \% \end{array}\right.\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&D(i,j)=\left\{\begin{array}{l}0,{C}_{1}(i,j)={C}_{2}(i,j)\\ 1,{C}_{1}(i,j)\ne {C}_{2}(i,j)\end{array}\right.\end{eqnarray} \tag{ 57 }$$

where ${C}_{1}(i,j)$ and ${C}_{2}(i,j)$ are the pixel values of the original ciphertext image and the encrypted image obtained by slightly changing the plaintext image, respectively. For an image with a gray level of 256, the theoretical values of NPCR and UACI are 99.6094% and 33.4635% respectively [59].

In this section, one pixel of the plaintext Sailboat, Fruits, Pepper, Lena and 4.1.04–4.1.07 in each channel is selected randomly and changed by 1. Twenty times tests were conducted, and the average values of NPCRs and UACIs are shown in table 12. It is clear that the average NPCRs and UACIs calculated by our scheme are all close to the theoretical value 99.6094 and 33.4635, which illustrates that the designed scheme can withstand differential attacks effectively. In addition, the comparison between the average NPCRs and UACIs values of Lena image calculated in this paper and the existing literature [27, 28, 39] are shown in table 13. Obviously, the average NPCRs and UACIs values of the designed scheme are closer to 99.6094 and 33.4635, the plaintext sensitivity is stronger, and the ability to withstand differential attacks is superior.

5.7.1. Anti-noise performance

Since ciphertext images are easily disturbed by various noises during the communication transmission, an excellent cryptosystem should have enough ability to resist noise attack. To test the performance of the designed encryption scheme against noise attacks, Gaussian noise with intensities of 0.001, 0.003 and 0.005 are added to sailboat, fruit and pepper cipher images, respectively. Figure 14 shows the corresponding decrypted images. Moreover, salt and pepper noise with intensity of 0.1, 0.2 and 0.3 are added respectively, and figure 15 shows the decrypted images. The noise test results show that the quality of the recovered image gets worse as the noise intensity increases. However, the main information of original images can still be obtained. In conclusion, the designed algorithm can withstand Gaussian noise attack and salt and pepper noise attack perfectly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.7.2. Anti-shear performance

In addition to noise interference, ciphertext images are also easy to lose data in the process of transmission or storage. Therefore, a good encryption algorithm should also have sufficient resistance to clipping attack. To measure the ability of the designed algorithm to withstand shear attack, different degrees of data of Sailboat, Fruits and Pepper ciphertext images are cut out respectively. The ciphertext images with different degrees of data loss and the corresponding recovered images are given in figures 16 and 17, respectively. One can see that for different degrees of data loss, decrypted images can display main information of original images. The results show that the algorithm designed in this paper can withstand the clipping attack effectively and then can avoid the decryption difficulty caused by data loss perfectly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size