Blind Frequency Estimation and Symbol Recovery for the Analytically Solvable Chaotic System

The analytically solvable chaotic system (ASCS) is a promising chaotic system in chaos communication and radar fields. In this paper, we propose a maximum likelihood estimator (MLE) to estimate the frequency of ASCS, then a difference-integral (DI) detector is designed with the estimated frequency, and the symbols encoded in the signal are recovered. In the proposed method, the frequency parameter is estimated by an MLE based on the square power of the received signal. The Cramer-Rao lower bound in blind frequency estimation and the bit error performance in symbol detection are analyzed to assess the performance of the proposed method. Numerical results validate the analysis and demonstrate that the proposed symbol detector achieves the error performance with a little cost of 1 dB compared to the coherent detector. The robustness of the proposed method towards parameters is also verified through simulations.

to provide superior reliability and security compared with the conventional chaos direct spread spectrum system.
Although ASCS shows superior performance in communication and radar applications, the ASCS signal can be considered as simply a BPSK communication waveform with an unconventional basis function so that it has no inherent security. The dynamical characteristics and the exact solution of such a class of ASCS signal make it possible to incoherently recover the symbols conveyed in the chaotic signal even if it is severely contaminated by noise. In this paper, we propose an efficient method for estimating the frequency of ASCS, and for detecting the symbols encoded in the waveform. The system parameters of ASCS and the binary symbols encoded in the continuous chaotic waveform can be estimated directly from the received signal via a difference-integral (DI) detector. For a potential attacker, this method process requires only the knowledge of the form of ASCS, which is practical in implementation.
Based on the square power spectrum of the received signal, we first proposed a maximum likelihood estimator (MLE) to obtain the frequency parameter of the ASCS. For the symbol detection, we proposed a DI detector in which the binary symbols are retrieved by cumulatively summing the difference between the received signal and its time-delayed version. This detector requires no parameters of the transmitter other than the basis frequency, which can be estimated. The advantages of the proposed method are two-fold: first, it does not require prior knowledge at the transmitter other than the form of ASCS; second, although the error performance of the proposed symbol detector is around 1 dB worse than that of the coherent detector for ASCS, it provides a new tool for a potential adversary hoping to incoherently detect an ASCS signal.
The rest of this paper is organized as follows: In Section 2, we provide a brief review of the system description and exact solution of ASCS. In Section 3, the MLE based on the square power spectrum is proposed to estimate the basis frequency of ASCS. The performance of the estimator is also evaluated. In Section 4, a DI detector is proposed with detailed theoretical derivations. In Section 5, a closed-form expression for error probability in symbol detection is derived and verified through simulations. The robustness of the detector is also discussed. Finally, some conclusions are made in Section 6.

The Analytically Solvable Chaotic System
The standard form of ASCS is described with a differential equation as [1] u − 2βu + (ω 2 + β 2 )(u − s) = 0 (1) where ω determines the symbol period as T = 2π/ω and β ∈ (0, T −1 ln 2) is the damping coefficient of the linear part of ASCS system. The waveform u(t) ∈ R is continuous and state s ∈ {±1} is discrete. The state can be controlled by a binary symbol sequence s(t) which is set as the sign of current u(t) once the switching condition meets and is held constant until the next time when the switching condition triggers. The switching condition and state change process are described aṡ The chaotic signal u(t) can be obtained by the numerical solution of the Equation (1) with an adjustable step size Runge-Kutta integrator. Suppose the initial conditions are |u(0)| < 1 andu(0) = 0, the exact analytic local solution of Equation (1) for the segment of nT ≤ t < (n + 1)T is derived as with the iteration relation u n+1 = e 2βπ/ω u n − (e 2βπ/ω − 1)s n , where n is an integer. Figure 1 shows the corresponding phase-space projection and the return map of the iterated shift relation at regular return times nT where n is an integer. The system parameters are T = 1 and β = ln 2. As Figure 1a suggests, ASCS has an attractor topologically similar to the well-known Lorenz attractor. As Figure 1b shows, the sampled values at nT are divided into two groups by the binary symbol sequence s(t) and follow the standard iteration relation lines with the fixed slope of e 2βπ/ω . Hence, the piecewise linear return map implies the dynamical system generated from Equation (1) is chaotic with Lyapunov exponent λ = β. Based on the iteration relation of ASCS, the Nth return point u N and its former symbols determine the initial value u 0 , which is given by From the analytic solution Equation (3), u(t) can be derived as a linear convolution of a basis function and a binary symbol sequence when N approaches to infinity: where m is an integer and the basis function P(t) is defined as Figure 2 shows the shape of the basis function P(t) for T = 1 and β = ln 2. As Equation (5) implies, the continuous waveform of u(t) can be generated by the sum of binary symbols and the time delayed basis function P(t). Thus, binary symbols can be encoded into the continuous chaotic waveform u(t).

Frequency Estimation for Chaotic Signal
In ASCS, the frequency parameter ω determines the symbol rate, which is significantly important in symbol detection. In the following, we derive the amplitude spectrum of the square power of u(t) and then propose a frequency estimator based on it.

Amplitude Spectrum of the Square Power of the Chaotic Signal
Assume the basis radian frequency of ASCS is ω 0 . Based on the analytical solution Equation (3), the local waveform of u(t) for mT ≤ t ≤ (m + 1)T is where D m = u m − s m ; A = 1 + β 2 /ω 0 ; ϕ = arctan β/ω 0 . The feature of the amplitude spectrum of the ASCS has been calculated and analyzed in References [4,15]. The spectrum of the chaotic waveform u(t) shows zeros at the fundamental frequency and its harmonics. Although this conclusion is determined, the frequency component with zero amplitude is difficult to locate accurately as the zero zones are often too vague to detect, especially when the signal-to-noise-ratio (SNR) is low. In the following, we derive the amplitude spectrum of the square power of u(t) and then propose a frequency estimator based on it. The square power of u(t) is Since s 2 m = 1 is constant and contributes to the spike in the zero frequency component, we subtract it from v(t) and have w(t) = v(t) − 1 as where w 1 (t) = 2s m D m e β(t−mT) A cos(ω 0 t + ϕ); w 2 (t) = 1 2 D 2 m A 2 e 2β(t−mT) (cos(2ω 0 t + 2ϕ) + 1). Then we derive the Fourier transform of each part. The Fourier transform for w 1 (t) is Since the integral is difficult to derive, we calculate the Fourier coefficients as and By computing the integral, we have a(ω) as As m changes, the variables θ m , λ m , δ m , and η m change gradually. Take the cos θ m as an example. Since θ m is ergodic When ω = nω 0 , however, the variables θ m , λ m , δ m , and η m keep constant since mTω = 2mnπ. The analysis of b(ω) is similar to that of a(ω) (see Equation (A1) in Appendix A). After obtaining the Fourier transform coefficients, the amplitude spectrum is then calculated by |W 1 (ω)| = a 2 (ω) + b 2 (ω) and the approximate result for ω = nω 0 is where Since s m D m is uniformly distributed in [−1, 0], the average of |W 1 (ω)| for ω = nω 0 is fixed at AI 1 (n). The expression of I 1 (n) implies the existence of spikes at the integer multiplies of the basis frequency ω 0 and the amplitude of the spikes shrinks as n increases. Hence, the single-side amplitude spectrum |W 1 (ω)| shows the maximum spike at n = 1.
The analysis of the amplitude spectrum of w 2 (t) is similar to the above and the exact analytic expression of |W 2 (ω)| for ω = nω 0 is given by where The detailed derivation is provided in Appendix A. Note that the random variable D 2 m follows a uniform distribution in [0, 1], hence the average of |W 2 (ω)| is fixed at 2 ) β 2 +ω 2 0 I 2 (n). As the expression of I 2 (n) implies, the spikes exist at the integer multiplies of the basis frequency. When n = 0, the partial spectrum sees a maximum spike which is caused by the direct current component in w 2 (t). When n = 1, the spectrum sees a second maximum spike.
By summing the |W 1 (ω)| and |W 2 (ω)| together, we have the final amplitude spectrum of w(t). Figure 3a,b show a typical amplitude spectrum and the theoretical amplitude spectrum of the chaotic signal respectively, with the basis frequency f 0 as 100 Hz. From the comparison, the amplitude spectrum obtained by fast Fourier transform (FFT) matches the theoretical one well, proving the derivations above are reasonable and valid. The spikes occur at n f 0 (n = 1, 2, . . . ), that is, the inter multiplies the basis frequency and shows a damping amplitude as n increases. Despite the zero frequency part, the maximum spike occurs exactly at 100 Hz and the second highest shows at 200 Hz. Obviously, the basis frequency can be estimated from the maximum spike in the amplitude spectrum of the square power of the chaotic signal.

Maximum Likelihood Estimator and Performance Analysis
Based on the analysis of the square power of the chaotic signal, we proposed a maximum likelihood estimator (MLE). Suppose the observed signal with a time length of MT is r = (r 0 , r 1 , . . . , r N ), with for k = 0, 1, 2, . . . , N. The sample is u k = u(k f s t) and f s is the sampling frequency; N = MT f s is the number of samples; n k is the additive Gaussian noise with zero mean and variance of σ 2 . We calculate the square power of the observed signal and subtract its direct current component as where the first term contains signal part while the second and the third term are noise-related. The last term E(r 2 k ) is the average value of the square power of the received signal. Since the term 2u k n k has zero mean, we have E( The basis frequency of the chaotic signal is estimated by an MLE as According to the Fourier transform of x(t), we have E( ω 0_MLE ) = ω 0 . Hence, the estimator is unbiased. To assess whether the proposed unbiased estimator is asymptotically efficient or not, we continue to derive the Cramer-Rao lower bound (CRLB) in the following. Note that the segment of the where D m = u m − s m . For simplicity in derivation, we assume the noise term distributes normally with zero mean and variance of σ 2 x = 2σ 4 + 2σ 2 un , where σ 2 un is the variance of the term 2u k n k . Like the similar process in parameter estimation with white Gaussian noise, we have the Fisher information of the estimated frequency as where the cumulative sum part is derived as Then the CRLB is expressed as

Symbol Detection with a Difference-Integral Detector
At the friendly receiver of ASCS, the symbols are recovered with a coherent matched-filter-based detector which requires precise knowledge of the frequency parameter f and the damping coefficient β. However, the parameter β, which is related to the linear part of the ASCS system, is difficult to estimate since the received signal is severely contaminated by noise. In this section, we propose a DI detector to incoherently retrieve the binary symbols that the chaotic signal is conveying, which requires only the frequency parameter estimation.
The whole process of the frequency estimation and symbol detection is described in Figure 5. The chaotic signal is generated by the ASCS oscillator with the initial value u 0 and the basis frequency f 0 . After transmitting in the channel, the chaotic signal is contaminated by the noise. Then we conduct a blind frequency estimation from the received signal and obtain the estimated basis frequencyf 0 . With the estimated frequencyf 0 , the binary symbols that the chaotic signal conveying can be retrieved by a DI detector. As mentioned above, the main characteristic of ASCS is that this class of chaotic signal can be described by a weighted sum of basis functions. Note that continuous waveform u(t) contains the relations among binary symbols. Hence, the symbols conveyed in the ASCS waveform can be retrieved directly from the received signal by utilizing this feature.
The block diagram of the DI detector scheme is shown in Figure 6. Assume a segment of continuous chaotic signal u(t) carries a sequence of binary symbols {s k }. The intermediate state d(t) is obtained by differencing the received signal v(t) and the delayed one v(t − T). Then d(t) is sent to the integrator. The integral output η(t − T) is computed numerically by t −∞ d(t + τ)dτ. Note that the practical signal generated from the chaotic oscillator is not infinite in length. Suppose the received signal transmit N symbols and thus the time duration is NT. Then the whole integral region is from 0 to NT. By extracting the integral output at sampling time t k = kT(k = 0, 1, . . . , N − 1), we compare the samples η k with a flexible threshold θ k , calculated with the previously recovered symbols to detect the polarity, which is given by θ k = −λAB ∑ k−1 m=0 s m e −mβT . If the η k is larger than θ k , then the estimated symbol is 1. Otherwise, the estimated symbol is −1. Hence, the estimated symbol sequence is expressed asŝ k = sgn{η k − θ k }. In the following, we provide detailed derivations of the proposed detector. The received signal in the AWGN channel is modeled as where the white Gaussian noise term n(t) has zero mean and power density N 0 /2. The difference between the received signal and its time-shift version can be expressed in a convolution form as Define an integral from 0 to t as where χ(t) is noise term; (t) is computed as where λ = e βT + e −βT − 2; A = β 2 ω 2 +β 2 ; B = 2 β ; ϕ = arctan 2βω ω 2 −β 2 . To guarantee the accuracy of the integral, the integral step should be small and thus the sampling rate of the received signal v(t) is supposed to be high. The sample of ε(t) at t = kT is derived as where the upper term is the cumulative interference caused by the previous symbols; the bottom term is the interference caused by the future symbols. Hence, we can rewrite the samples as where ξ is the sampled noise term χ(t) in Equation (28) and it follows a Gaussian distribution with zero mean and σ 2 ξ variance; Σ is the inter-symbol interference defined as where I p = −λAB ∑ k−1 m=0 s m e −mβT represents the intersymbol interference of the past symbols; I f = λAB(e kβT − 1) ∑ ∞ m=k+1 s m e −mβT represents the intersymbol interference of the future symbols. Since I p can be calculated but I f is unavailable, we use I p as the detection threshold and the detection process is described asŝ with θ k = I p as the threshold.
To illustrate the detection process above, we generate a segment of chaotic signal with m truncated at 50 and the period T is set as 1 s, where the chaotic signal is not contaminated by noise. The continuous waveform and symbol sequence are depicted in Figure 7a. Both the numerical η(t) and theoretical η th (t) integrals are calculated and compared in Figure 7b, where the sampling rate is 0.01 s. From the comparison of the waveforms, we find that the theoretical integral computed by Equation (29) is quite close to the numerical integral result, proving that the previous calculation and approximation are reasonable. The samples η k of the numerical integral result η(t) stay around ±1 which are depicted with two horizontal dash lines. Hence, the symbol sequence can be easily retrieved by Equation (33). Compared with the original symbol sequence provided in Figure 7a, the extracted one is exactly the same.

Performance Evaluation and Simulation
In this section, we first evaluated the symbol error performance theoretically and derived an exact closed-form expression of error probability in the symbol detection. Then the robustness of the DI detector was verified by numerical simulations with various errors in parameter estimation.
Under the noise background, errors in symbol detection occur with rising probability as the variance of noise increases. Hence, it is necessary to evaluate the error performance for symbol detection with different SNR levels. In this paper, we define the BER as the ratio of the number of errors to the total number of detected symbols. As implied in Equation (33), the detection is a classical hypothesis test. Then the probability of detecting an incorrect bit is P e = P(η k > θ k |s k = −1)P(s k = −1) + P(η k < θ k |s k = 1)P(s k = 1). (34) Note that the source symbols {s k } are sent with equal probability, that is, P(s k = −1) = P(s k = 1) = 1/2. Assume that "+1" and "−1" in the future symbols distribute with the equal probability of 1/2, then where erfc(·) denotes the complementary error function; Φ = −λABe −kβT + (e −βT − 1)AB + T.
Consider that e −kβT decreases rapidly as k increases, we have the approximate Φ = (e −βT − 1)AB + T. The integral is computed as where . We implement our symbol detecting method in different frequency settings and illustrate the BER curves in Figure 8. The theoretical BER curves of the symbol detecting method with optimal and suboptimal matched filter introduced in References [4,16] are also depicted as a reference. The SNR is defined as the ratio of the power of the ASCS signal to the noise power, where the sample rate is 100 times faster than the symbol rate. In fact, the relative SNR calculated with in-band noise should be the actual SNR added to 10 log 10 50 dB. As can be noticed, the simulated curve for the proposed DI detector matches the theoretical one well. Compared with the theoretical error performance of the suboptimal matched-filter-based detector, our method performs a degradation around 1 dB. Nevertheless, it is acceptable for an attacker or an eavesdropper to obtain enough information. As implied in the theoretical BER expression Equation (36), the BER is related to the damping coefficient β. To evaluate the BER performance under different damping coefficients, we set the damping coefficient as β = a f ln 2 where the coefficient a ∈ {0.6, 0.8, 1}. The simulated BER curves of the proposed DI-HCC over different damping coefficients are shown in Figure 9. The theoretical BER for a = 1 is depicted as a reference. From the comparisons of BER curves over different settings, the BER performance of the DI detector slightly degrades as the damping coefficient β decreases. However, the deterioration is not significant since β has a limited impact on the constant values in Equation (36). Hence, the proposed DI detector is robust to the damping coefficient of the ASCS. Since the samples are extracted at the integer multiplies of T, the imprecise sampling time is supposed to have effects on the BER performance. Assume that the integral η(t) are sampled at t = kT + ∆T, where ∆T is the sampling lag. Figure 10 shows the BER curves over different sampling lags with ∆T ∈ {0.10, 0.15, 0.20, 0.25} ms where the basis frequency of the ASCS is set as 1 kHz. From the comparison of the curves, the performance deteriorates gradually as the lag increases. However, the deterioration is slight when the lag is less than 0.15 ms, meaning the proposed detector keeps valid within a limited range of sampling lag.  The key parameter in symbol detection is the symbol period T, which is defined as T = 1/ f 0 and the basis frequency f 0 is estimated by the FFT-based MLE. Unavoidably, the frequency estimation error has a negative impact on the symbol detection. Figure 11 illustrates BERs for various estimation error ∆ f with different numbers of symbols in the received signal with a fixed SNR as −5 dB, where the basis frequency is set as f 0 = 10 Hz and the estimated frequency is f 0 + ∆ f . As the absolute value of the frequency estimation error increases, the symbol detection error probability is larger when N is 200 while almost no error is seen for N = 50 and N = 100 in the tested estimation error range. It can be explained that the DI detection is a process with error accumulation over time. As time passes, the sampling time will gradually deviate from the original exact sampling positions, which brings deterioration to the symbol detection. Assume the difference between the estimated symbol period T = 1/( f 0 + ∆ f ) and the exact period T is ∆T = T − T. Suppose the error of sampling time increases to the half of T after N symbols, then we have N∆T = T/2 and the estimation error of the frequency is correspondingly ∆ f = f /(2N − 1). Hence, the maximum of the absolute value of the frequency estimation error should be less than ∆ f . In the simulation in Figure 11, as derived theoretically, the absolute value of ∆ f should be controlled within 0.025, 0.05, and 0.1 for N is set as 200, 100, 50 respectively, to avoid the accumulated error. The obtained BERs show that the theoretical bounds of the frequency estimation error for different N match the simulated results well. In the practical implementation of the proposed method, the frequency estimation error can be evaluated by CRLB with the known SNR level. With the theoretical frequency estimation error, we can compute the upper bound of N. Then the received signal is divided into segments with the limited time length of NT to avoid error caused by the accumulation effect.

Conclusions
In this paper, an efficient method of blind frequency estimation and symbol recovery for ASCS is proposed. To obtain the frequency parameter of the chaotic signal, we first proposed a blind frequency estimator based on the square power of the received signal. Then we designed a DI detector by integrating the difference between the received signal and its time-delayed version. The binary symbols conveyed in the chaotic signal can be retrieved directly from the sampled integral results. The theoretical analysis of the proposed scheme is provided with detailed derivations and performance evaluation. For the frequency estimation, the CRLB is derived to assess the lower bound on the error variance. For the symbol detection, the close-form BER expression is derived to evaluate the error in symbol detection. Moreover, the robustness of our proposed method in various parameter settings is verified through numerical experiments. Both the theoretical derivation and simulations prove that the proposed method performs well against strong noise.

Conflicts of Interest:
The authors declare no conflict of interest.