Blind Estimation Methods for BPSK Signal Based on Duffing Oscillator

To realize the blind estimation of binary phase shift keying (BPSK) signal, this paper describe a new relational expression among the state of Duffing oscillator excited by BPSK signal, the pseudo-random code of BPSK signal, and the difference frequency between the to-be-detect signal and internal drive force signal of Duffing oscillator. Two output characteristics of Duffing oscillators excited by BPSK signals named implied periodicity and pilot frequency array synchronization are presented according to the different chaotic states of Duffing oscillator. Then two blind estimation methods for the carrier frequency and pseudo-random sequence of the BPSK signal are proposed based on these two characteristics, respectively. These methods are shown to have a significant effect on the parameter estimation of BPSK signals with no prior knowledge, even at very low signal-to-noise ratios (SNRs).


Introduction
In modern radar electronic warfare, reconnaissance and anti-reconnaissance efforts, as well as the evaluations of the reconnaissance effect, are of great importance. Obtaining useful information from enemy radar signals can provide guidance for radar jamming and fire destruction and enable effective countermeasures to be taken against enemy interference [1,2]. In the field of radar and communication, binary phase shift keying (BPSK) signals are commonly used because of their high frequency band utilization, strong anti-noise interference ability, wide signal bandwidth, and resistance to detection [3][4][5][6]. In actual battlefield situations, there is often a complex electromagnetic environment and electronic countermeasures under non-cooperative conditions, and so useful signals are often drowned in strong noise, complicating the tasks of efficient intercepted efficiently and accurate estimated. For these reasons, blind detection of the weak BPSK signals at low signal-to-noise ratios (SNRs) is an important subject.
Part of the literatures focuses on the estimation of carrier frequency. Methods based on the cyclic spectrum density [7], second-order cyclic statistics [8], and stochastic resonance (SR) [9] have been developed to estimate carrier frequency and symbol rate of BPSK signal, but these methods cannot calculate the pseudo-random code, and their estimation accuracy is poor in low SNR environments (the SNR should be greater than −6 dB, −2 dB, and 0 dB, respectively). The other part of the research estimates the pseudo-random sequence. Wang et al. [10] estimate the period and staring bit of pseudo-random sequence through the method of reprocessing the power spectrum density, segmentation processing, and an average cross-correlation calculation, which can be realized even the SNR = −15 dB, but the carrier frequency needs to be known in advance.
Synchronous demodulation method can identify pseudo-random codes effectively and provide easier parameter estimation algorithm [11], but this approach needs to exact the carrier frequency before estimating the pseudo-random code. Some estimators based on cyclic spectral density function were proposed in [12,13], which has less computationally intensive, but the chip time width and carrier frequency should be known. In the above methods, only carrier frequency estimation is not enough, but the estimation of pseudo-random code requires prior knowledge or additional calculation, which requires time accumulation or multiple processing. Moreover, most of these methods have poor performance at low SNR and cannot be applied to the complex electromagnetic environment.
Because of the sensitivity to regular signals and immunity to noise under certain conditions, Duffing oscillators have often been used to detect weak signals [14][15][16][17]. The intermittent chaotic states in chaotic systems are known to enhance the practicability of weak signal detection in chaotic systems [18]. Recently, many methods for estimating the carrier frequency of BPSK signals using a Duffing oscillator have been developed [19,20], which can have great performance at very low SNR due to the Duffing oscillator. However, the pseudo-random sequence cannot be detected. The pseudo-random sequence estimation using the output characteristics of the Duffing oscillator excited by a known carrier frequency signal are studied in [21,22]. These methods are more easily to be calculated than the usual, but the carrier frequency of the BPSK signal at low SNR is not generally known under non-cooperative conditions. In short, these methods based on Duffing oscillator are limited in some way, but there have been few attempts to estimate both the carrier frequency and pseudo-random code of BPSK signals without prior knowledge using the Duffing oscillator.
In this paper, the relational expression among the state of Duffing oscillator, the phase code and the cosine function of the difference frequency in the intermittent chaotic state of the Duffing chaotic oscillator excited by BPSK signal is derived. Two parameter estimation methods for the joint estimation of carrier frequency and pseudo-random code based on this output characteristics are presented. These methods enable exact parameter estimations to be generated, even when no prior knowledge is available. Simulation and experiment results show that these methods have low computation complexity and high estimation precision, remain feasible in low SNR environments, and eliminate the high dependence of the traditional BPSK signal estimation method on the signal carrier frequency.
The remainder of paper is organized as follows. The next section analyzes the output characteristics of the BPSK signal in the Duffing oscillator system and deduces their relationship. In Section 3, the two blind parameter estimation methods are described. The simulation and experiment results in Section 4 demonstrate the feasibility of the proposed methods. Finally, the conclusions to this study are presented in Section 5.

Relationship among Functions in Duffing Oscillator System under Intermittent Chaotic State Excited by BPSK Signal
Traditional weak signal detection systems based on Duffing oscillator are excited by a standard sinusoidal signal, which can be expressed as . x = ωy . y = ω(−ky + ax − bx 3 + γ r cos(ωt) + γ s cos((ω + ∆ω)t + ϕ 0 ) + n(t)) (1) where x is displacement, k is the damping coefficient, and −ax(t) + bx 3 (t) is the nonlinear restoring force. As shown in previous works, if we fix k = 0.5, a = b = 1, the Duffing oscillator system is more stable and represents better chaotic states [11][12][13][14]. γ r cos(ωt) is a periodic internal driving force with an amplitude of γ r ≈ γ c , where γ c is the critical threshold. γ s cos((ω+∆ω)t + ϕ 0 ) is the to-be-detected signal. γ s is the amplitude, ϕ 0 is the initial phase, and ∆ω is the frequency difference between the internal drive force signal and to-be-detected signal. n(t) is the stochastic disturbance, which is considered to be Gaussian white noise in this letter. According to [18], the value of ∆ω can influence the equivalent driving force of the system, and then the chaotic states of Duffing oscillator can also change. If ∆ω = 0, after adding the signal to be measured, the system transforms from a chaotic state to a large-scale periodic state. If ∆ω 0, the system will be intermittently chaotic. When ∆ω ≤ 0.03ω, the intermittent chaotic state can be maintained regularly and stably. When ∆ω exceeds this limit, the intermittent chaotic state may be broken due to the insufficient maintenance time of the equivalent policy force. In weak signal detection, there is usually only a rough estimate for the frequency of the signal to be measured, which is often different from the drive signal frequency in the system. Thus, weak signal detection based on the intermittent chaotic state is universal and significant. Figure 1 illustrates the intermittent chaotic state of the Duffing oscillator excited by the standard sinusoidal signal. Here, ω = 100 MHz, ∆ω = 3 MHz, γ r = 0.826, γ s = 0.1, and ϕ 0 = 0.
Sensors 2020, 20, x FOR PEER REVIEW 3 of 20 maintained regularly and stably. When △ω exceeds this limit, the intermittent chaotic state may be broken due to the insufficient maintenance time of the equivalent policy force. In weak signal detection, there is usually only a rough estimate for the frequency of the signal to be measured, which is often different from the drive signal frequency in the system. Thus, weak signal detection based on the intermittent chaotic state is universal and significant. Figure 1 illustrates the intermittent chaotic state of the Duffing oscillator excited by the standard sinusoidal signal. Here, ω = 100 MHz, △ω = 3 MHz, γr = 0.826, γ s= 0.1, and φ0 = 0. BPSK signals can realize phase modulation using pseudo-random codes containing some sequence of 1 s and −1 s. They can be expressed as where φi = (0, π) is the phase code of the BPSK signal. Replacing the sinusoidal signal to be detected in Equation (1) with the BPSK signal, the new state equation of the Duffing oscillator is Similar to the sinusoidal signal, BPSK signals can also cause the Duffing oscillator to transit intermittently between the chaotic state and large-scale periotic state if 0 < |△ω| ≤ 0.03ω. However, the duration of the different states is affected by both the difference frequency △ω and the phase code φi. Hence, the Duffing oscillator does not have a stable and regular intermittent chaotic period, unlike for the sinusoidal signal. The state in this condition is illustrated in Figure 2. Like the sinusoidal signal, we also fix ω = 100 MHz, △ω = 3 MHz, γr = 0.826, γs = 0.1, and φ0 = 0. To estimate the parameters of the BPSK signal, we need further analysis of the output characteristics of the Duffing oscillator. For Equation (3), the equivalent driving force of the system can be written as γecos(ωt + θ), where 2 2 e r r s s Here, φ(t) = φ0 + φi. Because the to-be-detected signal is weak, so γs ≈ 0, γr ≈ γc, and γs << γr, the value of θ(t) can be neglected. Therefore, the γe(t) can be written as BPSK signals can realize phase modulation using pseudo-random codes containing some sequence of 1 s and −1 s. They can be expressed as where ϕ i = (0, π) is the phase code of the BPSK signal. Replacing the sinusoidal signal to be detected in Equation (1) Similar to the sinusoidal signal, BPSK signals can also cause the Duffing oscillator to transit intermittently between the chaotic state and large-scale periotic state if 0 < |∆ω| ≤ 0.03ω. However, the duration of the different states is affected by both the difference frequency ∆ω and the phase code ϕ i . Hence, the Duffing oscillator does not have a stable and regular intermittent chaotic period, unlike for the sinusoidal signal. The state in this condition is illustrated in Figure 2. Like the sinusoidal signal, we also fix ω = 100 MHz, ∆ω = 3 MHz, γ r = 0.826, γ s = 0.1, and ϕ 0 = 0.
Sensors 2020, 20, x FOR PEER REVIEW 3 of 20 maintained regularly and stably. When △ω exceeds this limit, the intermittent chaotic state may be broken due to the insufficient maintenance time of the equivalent policy force. In weak signal detection, there is usually only a rough estimate for the frequency of the signal to be measured, which is often different from the drive signal frequency in the system. Thus, weak signal detection based on the intermittent chaotic state is universal and significant. Figure 1 illustrates the intermittent chaotic state of the Duffing oscillator excited by the standard sinusoidal signal. Here, ω = 100 MHz, △ω = 3 MHz, γr = 0.826, γ s= 0.1, and φ0 = 0. BPSK signals can realize phase modulation using pseudo-random codes containing some sequence of 1 s and −1 s. They can be expressed as where φi = (0, π) is the phase code of the BPSK signal. Replacing the sinusoidal signal to be detected in Equation (1) with the BPSK signal, the new state equation of the Duffing oscillator is Similar to the sinusoidal signal, BPSK signals can also cause the Duffing oscillator to transit intermittently between the chaotic state and large-scale periotic state if 0 < |△ω| ≤ 0.03ω. However, the duration of the different states is affected by both the difference frequency △ω and the phase code φi. Hence, the Duffing oscillator does not have a stable and regular intermittent chaotic period, unlike for the sinusoidal signal. The state in this condition is illustrated in Figure 2. Like the sinusoidal signal, we also fix ω = 100 MHz, △ω = 3 MHz, γr = 0.826, γs = 0.1, and φ0 = 0.
According to Equation (7) and Table 1, although the intermittent chaotic state is not regular, it has some implicit periodicity, and the relationship among the state of the system at some point, the phase code ϕ i and the cosine of the difference frequency cos(∆ωt + ϕ 0 ) can be expressed as If the system is in the large-scale periodic state, γ e (t) > γ c , sgn(γ e (t) > γ c ) = 1; if it is in the chaotic state or critical state, γ e (t) ≤ γ c , and the value of sgn(γ e (t) ≤ γ c ) is −1 or 0. Thus, we can define Sys(t) as the state of Duffing oscillator: We fix Df (t) as the output state of the difference frequency when the input is the sinusoidal signal, that is D f (t) = 1, sgn(cos(∆ωt + ϕ 0 )) = 1 −1, else The phase code can be represented as According to Equation (8), the relationship of these functions can be written as As Sys(t), Df (t) and Pc(t) are all bi-valued functions with outputs of 1 or -1, if two of these three functions are known, the third can be determined by a simple multiplication. The equivalent forms of Equation (12) are D f (t) = Sys(t) × Pc(t) (13) Sensors 2020, 20, 6412

of 21
Pc(t) = Sys(t) × D f (t) (14) We wish to simulate an X-band radar BPSK signal with ω + ∆ω = 10 GHz; in this case, the data volume is very large and requires significant computation. Hence, we set ω + ∆ω = 100 MHz (which can be regarded as a radar BPSK signal after down-conversion), increase the calculation speed by a factor of 100, and set ∆ω = 3 MHz, ϕ 0 = 0. The pseudo-random sequence is m sequence with the symbol width of 300 ns. It can be represented as 10000111110101001100 over a period of 6 × 10 −6 s. Then, we obtain (a) ϕ i , (b) cos(∆ωt + ϕ 0 ) and (c) the output of Duffing oscillator y(t) in Figure 3. After binarization of the values 1 and -1, Figure 4 presents the time-domain waveforms of (a) Pc(t), (b) Df (t), and (c) Sys(t). This figure provides a visual representation of the multiplicative relation in Equation (12). can be regarded as a radar BPSK signal after down-conversion), increase the calculation speed by a factor of 100, and set △ω = 3 MHz, φ0 = 0. The pseudo-random sequence is m sequence with the symbol width of 300 ns. It can be represented as 10000111110101001100 over a period of 6 × 10 −6 s. Then, we obtain (a) φi, (b) cos(△ωt+φ0) and (c) the output of Duffing oscillator y(t) in Figure 3. After binarization of the values 1 and -1, Figure 4 presents the time-domain waveforms of (a) Pc(t), (b) Df(t), and (c) Sys(t). This figure provides a visual representation of the multiplicative relation in Equation (12).  This relationship provides a theoretical basis for blind parameter estimation of BPSK signal. In actual applications, we can easily find Sys(t) from the Duffing oscillator. Therefore, the key to parameter estimation is to obtain Df(t) or Pc(t). If we can find either one of them, we can also determine the other based on Equations (12)- (14). On the basis of this theory prerequisite, we propose two estimation methods in the next section. can be regarded as a radar BPSK signal after down-conversion), increase the calculation speed by a factor of 100, and set △ω = 3 MHz, φ0 = 0. The pseudo-random sequence is m sequence with the symbol width of 300 ns. It can be represented as 10000111110101001100 over a period of 6 × 10 −6 s. Then, we obtain (a) φi, (b) cos(△ωt+φ0) and (c) the output of Duffing oscillator y(t) in Figure 3. After binarization of the values 1 and -1, Figure 4 presents the time-domain waveforms of (a) Pc(t), (b) Df(t), and (c) Sys(t). This figure provides a visual representation of the multiplicative relation in Equation (12).  This relationship provides a theoretical basis for blind parameter estimation of BPSK signal. In actual applications, we can easily find Sys(t) from the Duffing oscillator. Therefore, the key to parameter estimation is to obtain Df(t) or Pc(t). If we can find either one of them, we can also determine the other based on Equations (12)- (14). On the basis of this theory prerequisite, we propose two estimation methods in the next section.  This relationship provides a theoretical basis for blind parameter estimation of BPSK signal. In actual applications, we can easily find Sys(t) from the Duffing oscillator. Therefore, the key to parameter estimation is to obtain Df (t) or Pc(t). If we can find either one of them, we can also determine the other based on Equations (12)- (14). On the basis of this theory prerequisite, we propose two estimation methods in the next section.

Parameter Estimation Method for BPSK Signals Based on Output Characteristics Including Implied Periodicity and Array Synchronization of Duffing Oscillator
In previous articles, BPSK signal detection and parameter estimation methods based on the Duffing oscillator have used special signals to increase the known information. However, these approaches require preprocessing to obtain more a priori information, which affects the real-time capability of the reconnaissance system, and are not applicable in all cases. Therefore, completing the blind estimation of the BPSK signal with less prior information, sometimes when only Sys(t) is known, is the focus of this study.

Parameter Estimation Method for BPSK Signals Based on Implied Periodicity
The output of Duffing oscillator excited by BPSK signal is affected by difference frequency ∆ω and phase code ϕ i . The vector diagram of driving forces of the Duffing oscillator under the influence of BPSK signal is shown in Figure 5. The γ e changes alternately periodically based on ∆ω. When ϕ i is converted from 0 to π, γ s is transferred to γ s ' , which causes the changes of numerical value relation between the new equivalent drive force γ e ' (t) and threshold γ c ; that is the reason why the intermittent chaotic state changes in this system.

Parameter Estimation Method for BPSK Signals Based on Output Characteristics Including Implied Periodicity and Array Synchronization of Duffing Oscillator
In previous articles, BPSK signal detection and parameter estimation methods based on the Duffing oscillator have used special signals to increase the known information. However, these approaches require preprocessing to obtain more a priori information, which affects the real-time capability of the reconnaissance system, and are not applicable in all cases. Therefore, completing the blind estimation of the BPSK signal with less prior information, sometimes when only Sys(t) is known, is the focus of this study.

Parameter Estimation Method for BPSK Signals Based on Implied Periodicity
The output of Duffing oscillator excited by BPSK signal is affected by difference frequency △ω and phase code φi. The vector diagram of driving forces of the Duffing oscillator under the influence of BPSK signal is shown in Figure 5. The γe changes alternately periodically based on △ω. When φi is converted from 0 to π, γs is transferred to γs ' , which causes the changes of numerical value relation between the new equivalent drive force γe ' (t) and threshold γc; that is the reason why the intermittent chaotic state changes in this system. Taking the case of Figure 5a as an example, letting the value of γs to be γ0 when γs = γc, we can get the vector diagram in Figure 6 by studying the moment of system state transition. As shown in Figure 6a, before the change of the value of φi, γe > γc, the system is in large-scale motion; the value of γe will become smaller as it rotates clockwise, and after the time of △α/△ω, γe = γc, the system will move to a chaotic motion. In Figure 6b, after the change of φi, γe ' < γc, the system is inchaotic motion; while after the time of △α/△ω, and it will become large scale motion. In summary, the state of Duffing oscillator system will transform before and after the conversion of φi, but it does not break the periodic vector motion of the system, the system will still change its state at the same time as the Duffing oscillator driven by the sinusoidal signal after the time of △α/△ω. According to these, although the intermittent chaotic state is not regular, it also has implicit periodicity. Taking the case of Figure 5a as an example, letting the value of γ s to be γ 0 when γ s = γ c , we can get the vector diagram in Figure 6 by studying the moment of system state transition. As shown in Figure 6a, before the change of the value of ϕ i , γ e > γ c , the system is in large-scale motion; the value of γ e will become smaller as it rotates clockwise, and after the time of ∆α/∆ω, γ e = γ c , the system will move to a chaotic motion. In Figure 6b, after the change of ϕ i , γ e ' < γ c , the system is inchaotic motion; while after the time of ∆α/∆ω, and it will become large scale motion. In summary, the state of Duffing oscillator system will transform before and after the conversion of ϕ i , but it does not break the periodic vector motion of the system, the system will still change its state at the same time as the Duffing oscillator driven by the sinusoidal signal after the time of ∆α/∆ω. According to these, although the intermittent chaotic state is not regular, it also has implicit periodicity.
while after the time of △α/△ω, and it will become large scale motion. In summary, the state of Duffing oscillator system will transform before and after the conversion of φi, but it does not break the periodic vector motion of the system, the system will still change its state at the same time as the Duffing oscillator driven by the sinusoidal signal after the time of △α/△ω. According to these, although the intermittent chaotic state is not regular, it also has implicit periodicity.   Figure 7 compares the numerical variation of Sys(t) and Df (t). The length between the adjacent red fine line is a period T of Df (t). No direct periodicity can be observed in Sys(t), whereas the seemingly irregular numerical change implies the periodicity of Df (t). The length of different segments in Sys(t) can be marked as 1~16 (the same lengths are ignored), and the sum of 2 + 3 + 4, 2 + 5 + 6, 2 + 7 + 8, 2 + 9 + 10, 2 + 11 + 12, 2 + 13 + 14, 2 + 15 + 16 (where, for example, 2 refers to all segments with lengths equal to the segment marked as 2) in this figure is exactly the period of Df (t). Thus, we first calculate the length of time that Sys(t) spends in the 1 or -1 states, then calculate the sum of lengths of every two and three adjacent segments (T 2i and T 3i , respectively) and find their minimum T 3min . Some of these summed lengths of adjacent segments will be less than the period T of Df (t), and some of them will exceed T. Statistically, however, it is more likely that the summed lengths of these segments will approximate the period T (we set ∆ω = 2π × 3 MHz, so the points of T is 3.333 × 10 −7 s), see Figure 8. Therefore, the period T of Df (t) can be obtained by summing the time length and using statistical methods. In this paper, we discard values of T 2i that are less than 0.9 × T 3min and obtain the period T of Df (t) by averaging the remaining T 2i .
Sensors 2020, 20, x FOR PEER REVIEW 7 of 20 Figure 7 compares the numerical variation of Sys(t) and Df(t). The length between the adjacent red fine line is a period T of Df(t). No direct periodicity can be observed in Sys(t), whereas the seemingly irregular numerical change implies the periodicity of Df(t). The length of different segments in Sys(t) can be marked as 1~16 (the same lengths are ignored), and the sum of 2 + 3 + 4, 2 + 5 + 6, 2 + 7 + 8, 2 + 9 + 10, 2 + 11 + 12, 2 + 13 + 14, 2 + 15 + 16 (where, for example, 2 refers to all segments with lengths equal to the segment marked as 2) in this figure is exactly the period of Df(t). Thus, we first calculate the length of time that Sys(t) spends in the 1 or -1 states, then calculate the sum of lengths of every two and three adjacent segments (T2i and T3i, respectively) and find their minimum T3min. Some of these summed lengths of adjacent segments will be less than the period T of Df(t), and some of them will exceed T. Statistically, however, it is more likely that the summed lengths of these segments will approximate the period T (we set △ω = 2π × 3 MHz, so the points of T is 3.333 × 10 −7 s), see Figure 8. Therefore, the period T of Df(t) can be obtained by summing the time length and using statistical methods. In this paper, we discard values of T2i that are less than 0.9 × T3min and obtain the period T of Df(t) by averaging the remaining T2i.  Assuming △ω > 0, the frequency difference can be obtained as △ω = 2π/T. A binary signal Df0(t) with the same period T can be constructed: Sensors 2020, 20, x FOR PEER REVIEW 7 of 20 Figure 7 compares the numerical variation of Sys(t) and Df(t). The length between the adjacent red fine line is a period T of Df(t). No direct periodicity can be observed in Sys(t), whereas the seemingly irregular numerical change implies the periodicity of Df(t). The length of different segments in Sys(t) can be marked as 1~16 (the same lengths are ignored), and the sum of 2 + 3 + 4, 2 + 5 + 6, 2 + 7 + 8, 2 + 9 + 10, 2 + 11 + 12, 2 + 13 + 14, 2 + 15 + 16 (where, for example, 2 refers to all segments with lengths equal to the segment marked as 2) in this figure is exactly the period of Df(t). Thus, we first calculate the length of time that Sys(t) spends in the 1 or -1 states, then calculate the sum of lengths of every two and three adjacent segments (T2i and T3i, respectively) and find their minimum T3min. Some of these summed lengths of adjacent segments will be less than the period T of Df(t), and some of them will exceed T. Statistically, however, it is more likely that the summed lengths of these segments will approximate the period T (we set △ω = 2π × 3 MHz, so the points of T is 3.333 × 10 −7 s), see Figure 8. Therefore, the period T of Df(t) can be obtained by summing the time length and using statistical methods. In this paper, we discard values of T2i that are less than 0.9 × T3min and obtain the period T of Df(t) by averaging the remaining T2i.  Assuming △ω > 0, the frequency difference can be obtained as △ω = 2π/T. A binary signal Df0(t) with the same period T can be constructed: Assuming ∆ω > 0, the frequency difference can be obtained as ∆ω = 2π/T. A binary signal Df 0 (t) with the same period T can be constructed: sgn(cos(∆ωt)) = 1 else (15) From Equations (10) and (15) we can see Df (t) and Df 0 (t) differ by an initial phase ϕ 0 . To eliminate all the numerical jumps in Sys(t) caused by Df 0 (t), it is necessary to find and compensate this initial phase ϕ 0 . We find the position in Sys(t) for which the length of the two adjacent segments is T, select one of these segments, and set its starting point as t 0 . Taking this point as a reference, then ϕ 0 satisfies where k is an arbitrary integer. Taking an appropriate value of k to shift Df 0 (t) to the right by ϕ 0 , the reconstructed signal Df rc (t) can be obtained. We can then replace Df (t) with Df rc (t) and obtain Pc rc (t) (the estimation of Pc(t)).
In the above process, we assume ∆ω > 0. However, if we want to determine the accurate carrier frequency, we need to make a judgment about the sign of ∆ω. In this letter, a simultaneous pseudo-random sequence estimation method for multiplex Duffing oscillators with different internal driving force is applied. Taking two channels as an example, two difference frequencies |∆ω 1 | and |∆ω 2 | can be obtained. By comparing their value and the different frequencies of the two internal driving forces, the positive and negative values of ∆ω 1 and ∆ω 2 can be determined, and then two estimated values of carrier frequency (ω + ∆ω 1 and ω + ∆ω 2 ) can be obtained. Averaging these values gives the estimated carrier frequency of the BPSK signal. The estimation accuracy of the pseudo-random sequence can be improved by increasing the number of oscillators.
The estimation method for the pseudo-random sequence and carrier frequency of the BPSK signal based on implied periodicity of Duffing oscillator is illustrated in Figure 9. As shown in Figure 9, we can get Sys(t) after binarization of the state identification of Duffing oscillator excited by the unknown BPSK signal. According to Sys(t), ∆ω can be obtained by the implied periodicity of Duffing oscillator, and then the carrier frequency can be estimated after judging the plus-minus sign of ∆ω. Besides, based on ∆ω and the phase, duty cycle of Sys(t), we can reconstruct Df (t). Finally, the pseudo-random sequence can be estimated based on the multiplication formula in Equation (14). This approach is used for the pseudo-random sequence estimation of BPSK signals in Figure 3. The time-domain waveform of several key nodes in the process is shown in Figure 10. Figure 10a shows the output of the Duffing oscillator system under the action of BPSK signal; Figure 10b is the binarized output Sys(t); Figure 10c is the reconstructed output state of frequency difference named Df rc (t) based on this method; and Figure 10d is the estimated result, Pc rc (t), given by this method.
A deburring method for Pc rc (t) is now studied. The burrs in Pc rc (t) are mainly caused by the small errors of the reconstructed signal Df rc (t) and the binary system output Sys(t) for the same intermittent chaotic state. Compared with the time-domain waveform diagrams of Pc rc (t) and Df rc (t), almost all the burrs in Pc rc (t) correspond to a state transition of Df rc (t). Therefore, as shown in Figure 11, for each state transition moment of Df rc (t), if Pc rc (t) has a pair of state transitions in a relatively short time around this moment, the two-state transitions are eliminated at the same time, resulting in pseudo-random sequences without burrs. It is clear that the sequence is 10000111110101001100, which is consistent with the pseudo-random sequence given in Figure 3.
Equation (14). This approach is used for the pseudo-random sequence estimation of BPSK signals in Figure 3. The time-domain waveform of several key nodes in the process is shown in Figure 10. Figure 10a shows the output of the Duffing oscillator system under the action of BPSK signal; Figure 10b is the binarized output Sys(t); Figure 10c is the reconstructed output state of frequency difference named Dfrc(t) based on this method; and Figure 10d is the estimated result, Pcrc(t), given by this method.    A deburring method for Pcrc(t) is now studied. The burrs in Pcrc(t) are mainly caused by the small errors of the reconstructed signal Dfrc(t) and the binary system output Sys(t) for the same intermittent chaotic state. Compared with the time-domain waveform diagrams of Pcrc(t) and Dfrc(t), almost all the burrs in Pcrc(t) correspond to a state transition of Dfrc(t). Therefore, as shown in Figure 11, for each state transition moment of Dfrc(t), if Pcrc(t) has a pair of state transitions in a relatively short time around this moment, the two-state transitions are eliminated at the same time, resulting in pseudo-random sequences without burrs. It is clear that the sequence is Dfrc(t), almost all the burrs in Pcrc(t) correspond to a state transition of Dfrc(t). Therefore, as shown in Figure 11, for each state transition moment of Dfrc(t), if Pcrc(t) has a pair of state transitions in a relatively short time around this moment, the two-state transitions are eliminated at the same time, resulting in pseudo-random sequences without burrs. It is clear that the sequence is 10000111110101001100, which is consistent with the pseudo-random sequence given in Figure 3. For carrier frequency estimation, two-channel Duffing oscillators are adopted with internal driving force signal frequencies of 98 MHz and 103 MHz, respectively. The measured frequency difference of oscillator 1 is 1.99854 MHz, and that of oscillator 2 is 2.99814 MHz. According to the numerical relation, the signal to be tested should be between the frequencies of oscillator 1 and oscillator 2. The estimated value of the carrier frequency is the average of the estimated values of the two oscillators, which in this case is 100.0002 MHz (oscillator 1: 99.99854 MHz, oscillator 2: 100.00186 MHz), an error of only 0.0002%. According to the results, the parameter estimation method for BPSK signals based on implied periodicity can realize high accurate blind estimation for a carrier frequency and pseudo-random sequence. For carrier frequency estimation, two-channel Duffing oscillators are adopted with internal driving force signal frequencies of 98 MHz and 103 MHz, respectively. The measured frequency difference of oscillator 1 is 1.99854 MHz, and that of oscillator 2 is 2.99814 MHz. According to the numerical relation, the signal to be tested should be between the frequencies of oscillator 1 and oscillator 2. The estimated value of the carrier frequency is the average of the estimated values of the two oscillators, which in this case is 100.0002 MHz (oscillator 1: 99.99854 MHz, oscillator 2: 100.00186 MHz), an error of only 0.0002%. According to the results, the parameter estimation method for BPSK signals based on implied periodicity can realize high accurate blind estimation for a carrier frequency and pseudo-random sequence.

Parameter Estimation Method for BPSK Signals Based on Pilot Frequency Array Synchronization
In a weak signal detection system, the use of a Duffing oscillator array is of great importance in many ways. By setting different internal driving force signals for each Duffing oscillator, the array can achieve multiply frequency measuring ranges and increased detection and estimation precision while eliminating the phase and frequency of blind areas.
When BPSK signals pass the Duffing oscillator array with different internal driving force signals, each oscillator output corresponding to different intermittent chaos period. However, all the oscillators are controlled by the same phase code ϕ i , the states of them will change at the same time. Therefore, the intermittent chaos states of array have a certain regularity, which can be used for the blind estimation of pseudo random sequences.
For each Duffing oscillator array, the internal driving force signal is γ r cos(ω j t). At this time, the binarization cosine function of the frequency difference between the internal driving force signal and the BPSK signal Df (t) varies with the frequency of the internal driving force signal sgn(cos(∆ω j t + ϕ 0 )) = 1 else (17) After inputting the BPSK signal, the output characteristics are converted into binary Sys j (t). Taking a simple array of two Duffing oscillators as an example, the binarization of the output characteristics Sys 1 (t) and Sys 2 (t) involves numerical conversion between 1 and -1 simultaneously affected by Pc(t), Df 1 (t), and Df 2 (t). Due to the difference in the period of Df 1 (t) and Df 2 (t), the time-domain waveforms of Sys 1 (t) and Sys 2 (t) are quite different. However, as the time of the output state transformation caused by Pc(t) in each Duffing oscillator in the array is always the same, there will always be a time when the numerical transformation occurs synchronously in Sys 1 (t) and Sys 2 (t). This time corresponds to the numerical transformation in Pc(t), as shown in Figure 12. Figure 12a is the pseudo-random sequence Pc(t). Figure 12b,c show the binarization results of Sys 1 (t) and Sys 2 (t) obtained after the BPSK signal modulated by Pc(t) has passed through two different Duffing oscillators in the array. The synchronous numerical transformations are marked with red ellipses, and these are consistent with the jump law of Pc(t). This constitutes the array synchronization of the output characteristics of the Duffing oscillator for the BPSK signals. The array synchronization obtained through the Duffing oscillator array with different frequencies is called pilot frequency array synchronization. 0 sgn(cos( )) 1 1 D ( )= 1 else After inputting the BPSK signal, the output characteristics are converted into binary Sysj(t). Taking a simple array of two Duffing oscillators as an example, the binarization of the output characteristics Sys1(t) and Sys2(t) involves numerical conversion between 1 and -1 simultaneously affected by Pc(t), Df1(t), and Df2(t). Due to the difference in the period of Df1(t) and Df2(t), the time-domain waveforms of Sys1(t) and Sys2(t) are quite different. However, as the time of the output state transformation caused by Pc(t) in each Duffing oscillator in the array is always the same, there will always be a time when the numerical transformation occurs synchronously in Sys1(t) and Sys2(t). This time corresponds to the numerical transformation in Pc(t), as shown in Figure 12. Figure 12a is the pseudo-random sequence Pc(t). Figure 12b,c show the binarization results of Sys1(t) and Sys2(t) obtained after the BPSK signal modulated by Pc(t) has passed through two different Duffing oscillators in the array. The synchronous numerical transformations are marked with red ellipses, and these are consistent with the jump law of Pc(t). This constitutes the array synchronization of the output characteristics of the Duffing oscillator for the BPSK signals. The array synchronization obtained through the Duffing oscillator array with different frequencies is called pilot frequency array synchronization. According to this property, we can input the BPSK signal into a group of Duffing oscillator arrays that have different internal driving force signals with different frequencies. After obtaining According to this property, we can input the BPSK signal into a group of Duffing oscillator arrays that have different internal driving force signals with different frequencies. After obtaining the binarization output characteristics, we can determine the position at which the numerical conversion takes place synchronously according to the array synchronization and obtain Pc(t). Subsequently, we can calculate the Df i (t) of the array according to Equation (12). Averaging Df i (t), we then obtain the estimated carrier frequency Df (t). For a pilot frequency Duffing oscillator array, because the frequency range of a single oscillator is about 0.03 ω, we set the frequency difference between each oscillator and the center frequency ω to be less than 0.03 ω.
We take the pilot frequency Duffing array as an example and propose a parameter estimation method for BPSK signals based on the synchronization of the array, as shown in Figure 13. As in Figure 9, we can get Sys j (t) after binaryzation of the state identification of Duffing oscillator array excited by the unknown BPSK signal. Then the pseudo-random sequence can be estimated by the pilot frequency array synchronization of Duffing oscillators. After obtaining Df i (t) based on the multiplication formula in Equation (13), we can complete the estimation of the carrier frequency.
According to the parameter estimation method in Figure13, the carrier frequency and seudo-random sequence of the BPSK signal in Figure 3 can be estimated. The signal to be etected is processed through an array formed by four Duffing oscillators with internal driving force requencies of 97 MHz (oscillator 1), 98 MHz (oscillator 2), 101 MHz (oscillator 3), and 103 MHz oscillator 4). The output of this system is shown in Figure 14.  According to the parameter estimation method in Figure 13, the carrier frequency and pseudo-random sequence of the BPSK signal in Figure 3 can be estimated. The signal to be detected is processed through an array formed by four Duffing oscillators with internal driving force frequencies of 97 MHz (oscillator 1), 98 MHz (oscillator 2), 101 MHz (oscillator 3), and 103 MHz (oscillator 4). The output of this system is shown in Figure 14.
Sensors 2020, 20, x FOR PEER REVIEW 11 of 20 the binarization output characteristics, we can determine the position at which the numerical conversion takes place synchronously according to the array synchronization and obtain Pc(t). Subsequently, we can calculate the Dfi(t) of the array according to Equation (12). Averaging Dfi(t), we then obtain the estimated carrier frequency Df(t). For a pilot frequency Duffing oscillator array, because the frequency range of a single oscillator is about 0.03 ω, we set the frequency difference between each oscillator and the center frequency ω to be less than 0.03 ω. We take the pilot frequency Duffing array as an example and propose a parameter estimation method for BPSK signals based on the synchronization of the array, as shown in Figure 13. As in Figure 9, we can get Sysj(t) after binaryzation of the state identification of Duffing oscillator array excited by the unknown BPSK signal. Then the pseudo-random sequence can be estimated by the pilot frequency array synchronization of Duffing oscillators. After obtaining Dfi(t) based on the multiplication formula in Equation (13), we can complete the estimation of the carrier frequency.
According to the parameter estimation method in Figure13, the carrier frequency and pseudo-random sequence of the BPSK signal in Figure 3 can be estimated. The signal to be detected is processed through an array formed by four Duffing oscillators with internal driving force frequencies of 97 MHz (oscillator 1), 98 MHz (oscillator 2), 101 MHz (oscillator 3), and 103 MHz (oscillator 4). The output of this system is shown in Figure 14.  After binarizing this array, the pseudo-random sequence estimation can be obtained based on pilot frequency synchronization. The results are shown in Figure 15. After binarizing this array, the pseudo-random sequence estimation can be obtained based on pilot frequency synchronization. The results are shown in Figure 15. The cosine function of binarized frequency difference obtained from the Duffing oscillator array according to Equation (12) is shown in Figure 16. As for the analysis of Figure 11, the arrows in this figure represent the time at which the value of the pseudo-random sequence shown in Figure 15e changes. If two state changes occur in the vicinity of this time within a relatively short period, they will be eliminated at the same time, and the estimated results after deburring can be obtained, as shown in Figure 17.  The cosine function of binarized frequency difference obtained from the Duffing oscillator array according to Equation (12) is shown in Figure 16. As for the analysis of Figure 11, the arrows in this figure represent the time at which the value of the pseudo-random sequence shown in Figure 15e changes. If two state changes occur in the vicinity of this time within a relatively short period, they will be eliminated at the same time, and the estimated results after deburring can be obtained, as shown in Figure 17. After binarizing this array, the pseudo-random sequence estimation can be obtained based on pilot frequency synchronization. The results are shown in Figure 15. The cosine function of binarized frequency difference obtained from the Duffing oscillator array according to Equation (12) is shown in Figure 16. As for the analysis of Figure 11, the arrows in this figure represent the time at which the value of the pseudo-random sequence shown in Figure 15e changes. If two state changes occur in the vicinity of this time within a relatively short period, they will be eliminated at the same time, and the estimated results after deburring can be obtained, as shown in Figure 17.

Simulation Experiment
To verify the feasibility of the two parameter estimation methods for BPSK signal, their performance in a Duffing oscillator weak signal detection system needs to be investigated. Thus, experiments were conducted to examine the estimation accuracy of the carrier frequency and the similarity of pseudo-random sequences under various SNRs. As a performance reference, the pseudo-random sequence was also estimated directly based on the known carrier frequency.
As in the simulation described in Section 3, the carrier frequency of the to-be-detected BPSK signal is 100 MHz, the symbol width of pseudo-random sequence is 300 ns, and signal amplitude is 0.6. The two-way Duffing oscillators adopted for the Duffing oscillator array based on implied periodicity have internal driving force frequencies of 103 MHz and 98 MHz. The Duffing oscillator array based on pilot frequency array synchronism is composed of four different Duffing oscillators with frequencies of 97 MHz, 98 MHz, 101 MHz, and 103 MHz, respectively. For the Duffing oscillator with the known carrier frequency, the driving force frequency is 100 MHz. The dynamic amplitudes of all dynamic amplitudes of Duffing oscillators are set to be close to the critical value, that is, γr = 0.826. These parameters are shown in Table 2.  Figures 18 and 19, when the SNR is −10 dB or −20 dB, the three methods all obtain good pseudo-random sequence estimation results. As SNR decreases, the difference between the chaotic state and the large-scale periodic state becomes smaller, and the

Simulation Experiment
To verify the feasibility of the two parameter estimation methods for BPSK signal, their performance in a Duffing oscillator weak signal detection system needs to be investigated. Thus, experiments were conducted to examine the estimation accuracy of the carrier frequency and the similarity of pseudo-random sequences under various SNRs. As a performance reference, the pseudo-random sequence was also estimated directly based on the known carrier frequency.
As in the simulation described in Section 3, the carrier frequency of the to-be-detected BPSK signal is 100 MHz, the symbol width of pseudo-random sequence is 300 ns, and signal amplitude is 0.6. The two-way Duffing oscillators adopted for the Duffing oscillator array based on implied periodicity have internal driving force frequencies of 103 MHz and 98 MHz. The Duffing oscillator array based on pilot frequency array synchronism is composed of four different Duffing oscillators with frequencies of 97 MHz, 98 MHz, 101 MHz, and 103 MHz, respectively. For the Duffing oscillator with the known carrier frequency, the driving force frequency is 100 MHz. The dynamic amplitudes of all dynamic amplitudes of Duffing oscillators are set to be close to the critical value, that is, γ r = 0.826. These parameters are shown in Table 2.  Figures 18 and 19, when the SNR is −10 dB or −20 dB, the three methods all obtain good pseudo-random sequence estimation results. As SNR decreases, the difference between the chaotic state and the large-scale periodic state becomes smaller, and the identification becomes more difficult. Figures 20 and 21 show the gradual appearance of several burrs that are difficult to remove. However, although the carrier frequency of the BPSK signal is unknown, the accuracy of the two blind estimation methods proposed in this paper approaches that of the traditional parameter estimation method with known carrier frequency. The pseudo-random sequence estimation also achieves good accuracy under the −35 dB SNR.
In parameter estimation, the correlation similarity coefficient is widely used to characterize the estimation accuracy of the pseudo-random sequence. This is the ratio of the peak value of the cross-correlation function between the obtained pseudo-random sequence and the original sequence to the peak value of the autocorrelation function of the original sequence. Table 3 presents the correlation similarity coefficients between the obtained pseudo-random sequences and the original sequences estimated by the three methods at different SNR. The pseudo-random sequences estimated by the three methods under −35 dB all have good cross-correlation performance with respect to the original sequence (correlation similarity coefficient >0.9). Relatively speaking, when the SNR is relatively high, the BPSK signal parameter estimation method based on array synchronization is closer to (and sometimes even better than) the estimation result of the special case in which the carrier frequency is known. When SNR is less than −35 dB, the correlation similarity coefficients obtained by the three methods decrease and take similar values.
The carrier frequency estimation results of the two parameter estimation methods proposed in the paper at different SNR are shown in Table 4. Under an SNR of −35 dB, the parameter estimation methods based on implicit periodicity and array synchronization achieve high precision. Parameter estimation method based on the implicit periodicity produces more accurate carrier frequency estimations. This probably because the implicit periodicity method estimates carrier frequency directly, whereas array synchronization first estimates the pseudo-random sequence, and then performs additional multiplication operations to get the carrier frequency, which will produce more errors. identification becomes more difficult. Figures 20 and 21 show the gradual appearance of several burrs that are difficult to remove. However, although the carrier frequency of the BPSK signal is unknown, the accuracy of the two blind estimation methods proposed in this paper approaches that of the traditional parameter estimation method with known carrier frequency. The pseudo-random sequence estimation also achieves good accuracy under the −35 dB SNR.   In parameter estimation, the correlation similarity coefficient is widely used to characterize the estimation accuracy of the pseudo-random sequence. This is the ratio of the peak value of the cross-correlation function between the obtained pseudo-random sequence and the original sequence to the peak value of the autocorrelation function of the original sequence. Table 3 presents the correlation similarity coefficients between the obtained pseudo-random sequences and the original sequences estimated by the three methods at different SNR. The pseudo-random sequences estimated by the three methods under −35 dB all have good cross-correlation performance with respect to the original sequence (correlation similarity coefficient >0.9). Relatively speaking, when the SNR is relatively high, the BPSK signal parameter estimation method based on array synchronization is closer to (and sometimes even better than) the estimation result of the special case in which the carrier  The simulation results show that the two parameter estimation methods proposed in this paper for the detection of weak BPSK signals based on the Duffing oscillator system can achieve high precision under −35 dB SNR. Under the condition of an unknown carrier frequency, these methods achieve similar estimation precision to the traditional approach in which the carrier frequency is known.

Semi-physical Simulation Experiment
In order to realize the blind parameter estimation verification of BPSK signal based on Duffing oscillator, a semi-physical experiment system was established. A BPSK signal can be sent by the pseudo-random code phase modulation prototype and received by an antenna. After down-conversion, the signal with intermediate frequency is fed into the computer and estimated by Duffing oscillator system. Subsequently, the BPSK signal can be reconstruct according to the detected parameters. Finally, we test the effect of parameter estimation by correlative processing between the received intermediate frequency signal and the reconstructed one. Images of spectrum of emitted BPSK signal, time-domain of the signal after down-conversion, time-domain of the reconstructed signals based on Duffing oscillator methods and the experiment result of correlation processing are provided in       According to Figure 22, the frequency of emitted BPSK signal is approximately 10 GHz. We down-convert it to adjust the frequency to 100 MHz. The pseudo-random sequence of the emitted BPSK signal can be expressed as 100000111010110111110001010010 over a period of 6 × 10 -6 s. As shown in Figure 23, each place where the phase flips represents a change in the pseudo-random   According to Figure 22, the frequency of emitted BPSK signal is approximately 10 GHz. We down-convert it to adjust the frequency to 100 MHz. The pseudo-random sequence of the emitted BPSK signal can be expressed as 100000111010110111110001010010 over a period of 6 × 10 -6 s. As shown in Figure 23, each place where the phase flips represents a change in the pseudo-random sequence. After parameter estimation based on Duffing oscillator implied periodicity and pilot frequency array synchronization, the estimated carrier frequencies are 99.992 MHz and 99.9873 MHz, respectively. The estimated pseudo-random sequences and reconstructed BPSK signals obtained by these two Duffing oscillator methods are presented in Figures 24 and 25. As shown in Figure 26, the maximum correlation coefficient based on two Duffing oscillator methods between the received signal and the reconstructed signal is 0.9401 and 0.9351, respectively, which proves the validity of these two methods.

Conclusions
In this letter, we have analyzed the output characteristics of BPSK signals through the Duffing oscillator and derived the multiplication relationship among the output of the Duffing oscillator, pseudo-random sequence, and cosine of the difference frequency, providing a theoretical basis for the blind estimation of BPSK signal. Based on it, we propose two blind estimation methods. Simulation and experiment results show that these two methods can estimate pseudo-random sequences and carrier frequencies with high accuracy and offer strong performance when SNR = -35 dB.
However, in our methods, pseudo-random sequences and carrier frequencies are closely related, if one parameter estimation has an error, it affects the other. Small errors have great effect on the estimation results. In addition, the binarization precision of the Duffing oscillator output signal in this paper is not high. In the future, we need to find a method to distinguish chaotic states in intermittent chaotic state of Duffing oscillator with less error, which can improve the overall estimation accuracy.