Design and Implementation of a Digital Dual Orthogonal Outputs Chaotic Oscillator

Abstract

Discrete time dynamical chaotic systems obey a set of recurrence equations involving one or more variables. Many chaotic maps have been proposed. None that simultaneously provides two sine–cosine outputs has stationary mean and standard deviation, or is quite robust with respect to the data format used in the hardware implementation. Here, we propose a chaotic oscillator based on a complex phasor whose angular argument evolves according to a geometric progression that is independent of the instantaneous amplitude. In order to maintain the oscillations, the phasor magnitude is normalized at each iteration using an approximation factor. The statistical characteristics of this oscillator are stationary in the short term, and do not depend on the initial conditions. The mean and standard deviation of the two orthogonal sequences quickly approach 0 and $1/\sqrt{2}$, respectively. The resulting distribution is similar to that of a digital sine with a constant angular step. We also present an FPGA architecture and its implementation results. This oscillator can be used in modulation schemes, such as the chaotic shift keying one or for data and image encryption. Finally, we show an original application that exploits the orthogonality of the two chaotic signals for the simultaneous encryption of two digital images.

1. Introduction

Chaotic maps offer a mathematical approach that describes the dynamic behavior of systems that are extremely sensitive to the initial conditions. Their study began in the 19th century with the works of Cauchy, Poincaré, Van der Pol and Liapunov. In 1963, Lorenz made a practical use of it for the forecasting of certain meteorological phenomena [1]. In continuous time, they can be represented mathematically by a system of nonlinear differential equations of dimension ≥3. The dynamics of Lorenz, Chen, Rössler and Lu are well-known examples [2]. Discrete time dynamical systems, such as Bernoulli, Hénon and the logistic maps, obey a set of recurrence equations involving one or more variables. More than 100 different chaotic maps have been proposed in the literature.

These systems generate aperiodic signals whose evolution, seen from the outside, lies on the edges between random and deterministic worlds. Therefore, the chaotic oscillators are of major interest for the generation of pseudo-random signals [3,4], or for the encryption of data [5,6,7], speech [8] or images [9,10,11,12]. Various modulation schemes, like the chaos shift keying, use a chaotic carrier in order to secure the transmission and to reinforce its confidentiality [13,14]. Chaotic sequences are also used in compressive sensing to build a measurement matrix that satisfies the restricted isometry property [15,16,17]. In the literature [18], the authors show that a chaotic generator can be used to create a reliable and efficient wireless communication system. Finally, chaotic oscillators and, more generally, systems with nonlinear dynamics, are used to perform complex operations and design stochastic logic systems [19].

Despite their abundance, there is no known chaotic system that simultaneously meets the following criteria:

■

jointly generates two orthogonal signals, i.e., the outputs are point by point in phase quadrature and constitute a unit vector.

■

presents statistical characteristics, such as stationary mean and standard deviation, that are also independent of the initial conditions.

■

is robust when one reduces the sample quantification format.

In an attempt to solve this problem, we looked for systems that implement circular functions that are implicitly normalized. In 1979, Ikeda was one of the first researchers to propose a discrete-time attractor using sine and cosine functions [20]. Then came the chaotic sine map (CSM) and chaotic cosine map (CCM) [21,22]. In order to increase the chaotic complexity, recurrences of the type $x_k=G\left[F(x_{k-1})\right]$, where G is a CCM and F a usual chaotic function, were studied more recently [23].

In all of these systems, the instantaneous phase is a function of the previous sample. Here, we propose a chaotic oscillator based on a complex phasor whose angular argument evolves according to a geometric progression that is independent of the instantaneous amplitude.

In the first part of this article, we review the most commonly cited discrete chaotic sequences and highlight their drawbacks, with respect to the three previously mentioned criteria.

The synthesis process and the main properties of the proposed chaotic oscillator are presented in the second part.

The third part analyzes the influence of the quantification on the performance of the oscillator.

Finally, a practical application for fast image cryptography is presented.

2. Usual Digital Chaotic Sequences

2.1. Recurrence Equations

In order to make digital oscillators, one solution is to discretize the differential equations. Numerical differentiation according to the Euler transformation $\frac{du(t)}{dt}\to \frac{u_k-u_{k-1}}{h}$ is often used. This method requires a sufficiently fine step h, otherwise the sequence generated fades out or diverges very quickly.

One- or multi-dimensional numerical sequences with a more or less pronounced chaotic character according to the internal parameters have also been proposed in the literature.

Table 1. Usual chaotic numerical time series.

Dim.	Model	Recursive Equations	Chaos Conditions
1D	Logistic	$x_k=\mu x_{k-1}\left(1-x_{k-1}\right)$	$x_0\in \left]0, 1\right[$ $3.57\le \mu <4$
1D	Bernoulli	$x_k=\frac{2}{1-\alpha }\left(x_{k-1}+1\right)-1\hspace{1em} -1<x_{k-1}<\alpha$ $x_k=\frac{2}{1-\alpha }\left(x_{k-1}-1\right)+1\hspace{1em}\hspace{1em}\alpha <x_{k-1}<1$	$-1<\alpha <1$
2D	Hénon	$x_k=y_{k-1}+1-a{x}_{k-1}^2$ $y_k=b{x}_{k-1}$	a = 1.4 b = 0.3
3D	Lorenz discrete (Euler approx. dt = h)	$x_k=x_{k-1}-\sigma h(x_{k-1}-y_{k-1})$ $y_k=y_{k-1}+h(\rho x_{k-1}-y_{k-1}-x_{k-1}{z}_{k-1})$ $z_k=z_{k-1}+h(x_{k-1}{z}_{k-1}-\beta z_{k-1})$	$\sigma =10$$\rho =28$ $\beta =8/3$h < 0.025

provides a synthesis of usual discrete chaotic sequences.

2.2. Influence of the Data Format

Chaotic sequences are very sensitive to the data quantification format. In fact, without special precautions, certain sequences fade very quickly. The Bernoulli map (also named the doubling map in the case of α = 0) vanishes to zero after about 30 or 60 iterates with single precision or double precision floating point computations, respectively. This can be easily demonstrated, as the recurrence relation can be viewed as a left shift operation with the loss of the most significant bit [24].

The Logistic map does not suffer from the vanishing problem with floating point numbers of both single or double precision; but with fixed point representation, there are always cycles in the sequence. The cycles have a duration ranging from a few iterates to a maximum of about 300 iterates for μ = 4, as well as a 2.21 fixed point format.

The sequences obtained by discretized differential equations are also very sensitive to the precision of computations, as well as to rounding rules, which may differ from one processor to another. For example, we show in Figure 1 the evolution over 60,000 points of the signal $x_k$ of the Lorenz sequence discretized according to the Euler approximation. The parameters used are $\sigma =10$, $\rho =28$, $\beta =8/3$, h = 0.001, $x_0=1$, $y_0=1$ and $z_0=0.3$. The first computation is carried out by a MATLAB script with the double precision floating point format, using the three recurrence equations provided in Table 1 (Figure 1, left). A second computation is carried out according to the same mathematical approach and with the same computer, but this time we used a Simulink diagram in discrete mode.

According to Figure 2, very strong differences appear from the 45,000th iterate. This dependency creates a problem when such sequences are used, for example, for cryptographic applications.

2.3. Statistical Variability

The statistical characteristics of chaotic sequences such as the mean (m) and standard deviation (std) may depend on the intrinsic parameters and initial conditions. This is particularly true for discretized differential equations. The variations of the latter being slow with respect to h, and the estimates of m and std require a large number of samples.

Table 2. Mean (m) and standard deviation (std) of the discrete Lorenz sequence for various initial conditions.

x0; y0; z0	0.1; 0.1; 0.1	0.5; 0.5; 0.5	1; 1; 1	2; 2; 2
m	−2.19	−0.64	−1.00	−1.44
std	7.72	8.00	7,94	7.86

, for example, shows the high variability of the mean of the discrete Lorenz sequence on 100,000 $x_k$ samples.

The Logistic and Hénon sequences have m and std values, independent of the initial conditions. Their values are 0.64 and 0.06, respectively, for the logistic sequence, with μ = 3.8 and 0.26 and 0.72 for the Hénon sequence. Nevertheless, the standard deviation of the logistic suite is strongly influenced by the parameter μ.

3. Proposed Orthogonal Chaotic Sequence

3.1. Theoretical Background

We consider a digital oscillator that simultaneously generates at time k the outputs x_k and y_k, given by the following:

$$x_k=cos(k\Delta \theta ) \mathrm{and} y_k=sin(k\Delta \theta )$$

(1)

where $\Delta \theta$ is the angular step and ($x_0=1$ and $y_0=1$) are the initial conditions. Equation (1) corresponds to the complex phasor, as follows:

$$W_k=x_k+j{y}_k=e^{j{\theta }_k}=e^{jk\Delta \theta }$$

(2)

The recurrence relation $W_k=e^{j\Delta \theta }{W}_{k-1}$ leads to the following:

$$\begin{array}{c}{x}_k=cos(\Delta \theta )x_{k-1}-sin(\Delta \theta )y_{k-1}\\ y_k=sin(\Delta \theta )x_{k-1}+cos(\Delta \theta )y_{k-1}\end{array}$$

(3)

This oscillator uses the principle of the angular integration depicted in Figure 2.

We consider a phasor whose argument ${\theta }_k$ follows a geometric progression with a common ratio of two. It follows the recurrence relations of Equation (4), as follows.

$$W_k=W_{k-1}^2 ➔ \begin{array}{c}{x}_k=x_{k-1}{x}_{k-1}-y_{k-1}{y}_{k-1}\\ y_k=y_{k-1}{x}_{k-1}-y_{k-1}{y}_{k-1}\end{array}$$

(4)

This is equivalent to substituting functions $y_{k-1}$ and $y_{k-1}$ for coefficients $cos\left(\Delta \theta \right)$ and $sin\left(\Delta \theta \right)$, respectively, in the previous diagram. If the initial phase ${\theta }_0$ is exactly equal to 0 or $2\pi /2^N$, the system stabilizes itself in the states ($x_k=1$ and $y_k=0$) and ($x_k=0$ and $y_k=1$).

On the other hand, the oscillator is very sensitive to the initial phase conditions. For example, Figure 3 shows the first 60 samples of $x_k$ and $y_k$ obtained for the two initial phases separated by a 0.001 rad (computations performed with MATLAB with double precision numbers).

3.2. Normalized Recurrence Equations

In practice, the phase errors due to the quantization accumulate and lead to the fade-out of signals. In order to overcome this problem, it is necessary to ensure that the relation $r_k^2=x_k^2+y_k^2=1$ is satisfied at each iteration. This can be achieved by using the normalized Equations (5).

$$x_k=\frac{x_{k-1}^2-y_{k-1}^2}{x_{k-1}^2-y_{k-1}^2} \mathrm{and} y_k=\frac{{2x}_{k-1}{y}_{k-1}}{x_{k-1}^2+y_{k-1}^2}$$

(5)

Equations (5) allow for the use of initial conditions other than the pair ($cos{\theta }_0$ and $sin{\theta }_0$). The values of $x_0\ne 0$ and $y_0\ne 0$, such as $x_0\ne y_0$, can also ensure the start of the oscillator. The orthogonalization and the normalization of the signals appear by the second iteration.

3.3. Approximation of Normalized Recurrence Equations

To eliminate the dividers in Equations (5), the following approximation is used:

$$x_{k-1}^2+y_{k-1}^2\approx 1+\epsilon ➔ \frac{1}{x_{k-1}^2+y_{k-1}^2}\approx 1-\epsilon \approx 2-(x_{k-1}^2+y_{k-1}^2)$$

(6)

Using this relation, we get the following:

$$r_k^2=[{(x_{k-1}^2-y_{k-1}^2)}^2+{4x}_{k-1}^2{y}_{k-1}^2]{[2-(x_{k-1}^2+y_{k-1}^2\left)\right]}^2=\left[{(x_{k-1}^2+y_{k-1}^2)}^2\right]{[2-(x_{k-1}^2+y_{k-1}^2\left)\right]}^2$$

(7)

Thus,

$$r_k=r_{k-1}^2\left(2-r_{k-1}^2\right)$$

(8)

Finally,

$$\{\begin{array}{l}{x}_k=\left(x_{k-1}^2-y_{k-1}^2\right)\left(2-(x_{k-1}^2+y_{k-1}^2)\right)\\ y_k=2x_{k-1}{y}_{k-1}\left(2-(x_{k-1}^2+y_{k-1}^2)\right)\end{array}$$

(9)

The asymptotic algorithmic complexity of our map is therefore O(n).

The magnitude of the generated points obeys Equation (8). Its evolution strongly depends on the initial value $r_0$. Only a convergence of $r_k$ towards 1 makes it possible to maintain a steady chaotic state.

The continuous function $f\left(w\right)=w^2(2-w^2)$ corresponding to Equation (8) has three positive roots, namely, 0, 1 and $\frac{1}{\phi }=\frac{\sqrt{5}-1}{2}=0.618$, where $\phi =1.618$ is the golden ratio.

The dynamic behavior of Equation (6) is analyzed by the mean of the function f(w). Values 0 and 1 are attractive points, and $1/\phi$ is a repulsive point. The main convergence intervals are summarized in

Table 3. Convergence properties of the sequence of Equation (6).

Case	Initial Magnitude	Convergence	Behavior	Examples
1	$r_0<1/\phi$	$r_k\stackrel{k\to \infty }{\to }0$	Chaotic Transient	$x_0=0.13$$y_0=0.6$$r_0=0.614$
2	$0.618<r_0<1.272$	$r_k\stackrel{k\to \infty }{\to }1$	Chaotic steady state	$x_0=0.6$$y_0=0.6$$r_0=0.848$
3	$1.272<r_0<\phi$	$r_k\stackrel{k\to \infty }{\to }0$ or $r_k\stackrel{k\to \infty }{\to }1$	Chaotic Transient or Steady state	$r_0=1.44$ ➔ transient $r_0=1.58$ ➔ steady state
4	$r_0>\phi$	$r_k\stackrel{k\to \infty }{\to }\infty$	Divergence

.

In conclusion, the values $x_0=cos\left({\theta }_0\right)$ and $y_0=sin\left({\theta }_0\right)$, and more broadly, all the initial values $(x_0,y_0)$, whose magnitude corresponds to case number two in Table 3, make it possible to generate a maintained chaotic oscillator.

3.4. Statistic Properties of the Chaotic Oscillator

As the proposed oscillator generates sine and cosine functions, its statistical characteristics are stationary in the short term. The mean and standard deviation of the two sequences quickly approach 0 and $1/\sqrt{2}$, respectively. The trajectories of the oscillator signals in different planes are presented in Figure 4.

where $dx=x_k-x_{k-1}$ and $dy=y_k-y_{k-1}$ are the differentials of x and y, respectively.

The first two diagrams show butterfly-shaped trajectories. The circular envelope of the (x, y) points is characteristic of quadrature signals.

The observed dynamics are similar to the Ikeda attractor in the (x, y) plane and to the Hénon attractor in the (dx and x) plane.

The histogram on 1000 points of x and y shows a non-uniform distribution of output values (Figure 5). This distribution is similar to that of a digital sine with a constant angular step. On the other hand, the instantaneous angular phase ${\theta }_k$ is uniformly distributed between −π and π, as shown in the histogram of Figure 6, made from a million values.

3.5. Chaotic Behavior Analysis

The recognition and quantification of chaos in a digital signal has been the subject of different approaches. Among them, the estimate of Lyapunov’s largest exponent λ of the progression is one of the most reliable. This parameter measures the average exponential rate of divergence or convergence between two trajectories of the sequence resulting from very close initial conditions. Three situations can occur, namely:

■

λ < 0 corresponds to an orbital trajectory attracted by a stable fixed point. This situation is representative of strongly deterministic, harmonic or random sequences.

■

λ very close to 0 represents a steady state close to a chaotic transition.

■

λ > 0 characterizes an unstable and chaotic orbit. The larger λ, the more chaotic character is observed.

In order to study the dynamics of financial markets, A Bensaida developed a neuron network-based chaos test algorithm that mimics the chronological behavior of the sequence to be analyzed [25]. We used this algorithm because it requires no prior knowledge of the signal to be analyzed, and can also be used on noisy data.

Table 4. Comparative of the Lyapunov coefficients.

Logistic	Hénon	Proposed Chaotic Oscillator
μ = 3.8 $x_0$ = 0.2	a =1.4; b = 0.3 $x_0$ = 1; $y_0=0$	${\theta }_0=0.01$	$x_0=0.5$ $y_0=0.8$
λx = 1.29	λx = 0.89 λy = 1.21	λx = 2.14 λy = 0.30	λx = 1.63 λy = 0.68

compares the Lyapunov coefficients recorded on the Logistic progression, Hénon sequence and the signals delivered by the proposed oscillator. The estimates are made on 1000 points, with the same configuration of the neural network.

This comparison confirms the chaotic character of the signals $x_k$ and $y_k$, with a higher grade of $x_k$.

4. Hardware Implementation

The oscillator was also implemented on a Xilinx FPGA target (Zynq UltraScale+ MPSoC ZCU106). The implementation diagram is given in Figure 7.

As the dynamics of $x_k$ and $y_k$ are limited to the interval [−1; 1], a quantization format with 23 bits, of which 21 are reserved for the fractional part, is sufficient to sustain the oscillations and to provide a reasonably large enough period (more than 2820 iterations), compared to a single precision floating point implementation (more than 1920 iterations).

Table 5. Implementation results summary on the FPGA Xilinx board.

Data Format	LUT	FF	DSP	Dyn. Power	fs max
Float 64	2999	128	40	35 mW	33 MHz
Float 32	1453	64	17	20 mW	50 MHz
Fix 23_21	33	32	6	6 mW	100 MHz

summarizes the necessary hardware resources, power consumption and the maximum sampling frequency f_s obtained for three different calculation accuracies, namely, double precision floating point (IEEE754 64-bit), single precision floating point (IEEE754 32-bit) and the custom defined fixed-point format, with two bits for the integer part and 21 bits for the fractional part.

In this case, one can see that an algorithm architecture matching the methodology provides a resulting architecture that lowers the necessary hardware resources, consumes less power, has a greater throughput and has a similar or even better accuracy than its direct “translation” from a software model.

5. Application to the Fast Image Encryption

5.1. Image Encryption/Decryption With Bitwise XOR Operation

Many authors use chaotic systems to encrypt digital images. Basically, the encryption processes of an image plane of N*M pixels remain similar, and take place in three successive stages, namely:

■

Elaboration of a sequence of N*M numerical coefficients of the chosen chaotic system in which the precise initial conditions have been introduced.

■

Quantification of the coefficients in a format identical to that of the pixels of the original image (eight bits in most cases). Given the variability of the statistical characteristics of the sequence with the initial conditions, this operation can only be carried out after the computation of the entire sequence.

■

Application of a bitwise XOR operation between the pixel values and the quantized coefficients.

The decryption of the image is carried out according to the same approach.

The use of this coding strategy with the proposed chaotic system is also possible. In this case, the quantification is easy because of the statistical stability of the sequence. Indeed, the range of our sequence is [−1; 1], therefore it becomes very easy to convert it into, for example, a sequence of eight-bit words. The hardware resources needed are very low, as it only requires a fixed-point increment whose output is truncated to the eight most significant bits. The scheme of this method with the proposed chaotic oscillator is shown in Figure 8.

5.2. Simultaneous Images Mixing Encryption/Decryption

We propose a second process, which can be used if the confidentiality of the information is not the main criterion.

Let U and V be two different grayscale images of the same size N*M.

Let k be the position of a pixel in the image, according to the conventional scanning from left to right and from top to bottom.

The $u_k$ and $v_k$ pixels of each image can be considered as the Cartesian coordinates of a point $P(u_k,v_k)$, respectively.

Let us apply a rotation of an angle ${\theta }_k$ to each point P.

This operation leads to point $P^{\prime }(u_k^{\prime },v_k^{\prime })$, such that:

$$\begin{array}{l}{u}_k^{\prime }=u_{k}cos\left({\theta }_k\right)+v_{k}sin\left({\theta }_k\right)\\ v_k^{\prime }=-u_{k}sin\left({\theta }_k\right)+v_{k}cos\left({\theta }_k\right)\end{array} ➔ \begin{array}{l}{u}_k^{\prime }=u_k{x}_k+v_k{y}_k\\ v_k^{\prime }=-u_k{y}_k+v_k{x}_k\end{array}$$

(10)

where the ordered pairs ($x_k$, $y_k$) are generated by the orthogonal outputs of the proposed chaotic oscillator.

The values $u_k^{\prime }$ and $v_k^{\prime }$ provide two encrypted images, named U’ and V’, respectively.

Two uniform images represented by a single point P are located after this transformation in a different point P’, distributed randomly on a circle of the same radius.

The transformation is reversible. It is performed by a reverse rotation, as follows:

$$\begin{array}{l}{u}_k=u_k^{\prime }{x}_k-v_k^{\prime }{y}_k\\ v_k=u_k^{\prime }{y}_k+v_k^{\prime }{x}_k\end{array}$$

(11)

This principle is similar to that of a synchronous modulation/demodulation with a random frequency.

The algorithm for the encryption is as follows:

Step 1:

Read the two original images with size N*M, then convert them into one dimensional arrays U and V, and assign the initial conditions x₀ and y₀.

Step 2:

For each pair (u_k and v_k), compute the following:

$u{\prime }_k=u_k{x}_k+v_k{y}_k$

$v{\prime }_k=-u_k{y}_k+v_k{x}_k$

$x_{k+1}=\left(x_k^2-y_k^2\right)\left(2-(x_k^2+y_k^2)\right)$

$y_{k+1}=2x_k{y}_k\left(2-(x_k^2+y_k^2)\right)$

Step 3:

Convert the two arrays U’ and V’ into two images with the size N*M.

The algorithm for the decryption is as follows:

Step 1:

Read the two encrypted images with size N*M, then convert them into one-dimensional arrays U’ and V’, and assign initial conditions x₀ and y₀.

Step 2:

For each pair (u’_k and v’_k), compute the following:

$u_k=u{\prime }_k{x}_k-v{\prime }_k{y}_k$

$v_k=u{\prime }_k{y}_k+v{\prime }_k{x}_k$

$x_{k+1}=\left(x_k^2-y_k^2\right)\left(2-(x_k^2+y_k^2)\right)$

$y_{k+1}=2x_k{y}_k\left(2-(x_k^2-y_k^2)\right)$

Step 3:

Convert the two arrays of U and V into two images with the size N*M.

The asymptotic algorithm complexity of both the encryption and decryption is O(N*M).

We applied this encryption method using the digital architecture described in the previous paragraph. Figure 9 and Figure 10 illustrate the encryption and decryption, respectively, of the original images (Lena and the Pirate, respectively).

Upon reception, a visualization of the addition of the two images can be carried out very simply by performing the computation $\sqrt{{u^{\prime }}_k^2+{v^{\prime }}_k^2}$, which is equal to $\sqrt{u_k^2+v_k^2}$.

However, the separation of the two images cannot be achieved without knowledge of the chaotic sequence.

A practical application of this process is merging pairs of pages from a book in order to provide a quick overview to the reader, who, if interested, can request the decryption key from the publisher.

6. Conclusions

We propose a chaotic oscillator using a complex phasor with phase doubling at each sampling time. The recursive structure performs the angular integration and simultaneously generates two orthogonal signals representative of the pair (cosine and sine). The nonlinear evolution of the phase is responsible for the chaotic behavior that we confirmed by the positive Lyapunov coefficients evaluated on the generated signals.

Thanks to the trigonometric properties related to the double angle, no direct computation of the cosine and sine functions is necessary.

The permanence of the oscillations is obtained using a normalization of the magnitude of the ordered pairs ($x_k$ and $y_k$) carried out at each iteration.

In order to simplify the implementation of the system, the normalization ratio is replaced by an approximation factor. The recurrence equations thus take a purely polynomial form.

The analysis of the dynamics of the system shows that any initial pair ($x_0$ and $y_0$) whose magnitude is between 0.618 and 1.272 ensures a chaotic behavior.

The hardware implementation of this oscillator is presented and the influence of the quantification format on the chaotic complexity is studied. A fixed point 23-bit format, including 21 bits for the fractional part, allows for the generation of deterministic chaos with a periodicity greater than 2800.

This oscillator can be used in modulation schemes, such as the chaotic shift keying one and for data and image encryption. We finally show an original application that exploits the orthogonality of the two chaotic signals for the simultaneous encryption of two digital images.