1. Introduction
Our society now hugely depends on space technologies and space-based capabilities. Space, however, is becoming increasingly congested due to numerous spacecraft and orbital debris, presenting severe challenges for monitoring significant space assets [1]. Space Situational Awareness (SSA) is thus becoming increasingly critical in protecting space assets and safeguarding the space environment.
At present, the most widely used methods for monitoring, surveillance, and observation in SSA is radar and optical telescope [2]. Nevertheless, limited by cost, atmospheric conditions, and operation locations [3], all of these methods are inadequate and could not provide all-weather, all-time measurement. In contrast, spacecraft monitoring technology based on passive interferometry, especially the Connected Element Interferometry (CEI) has the unique advantages, which is very suitable for enhancing the existing monitoring means in SSA. If CEI is integrated with conventional means such as radar and optical telescope, prominent improvement will be achieved on the surveillance and monitoring of high-value spacecraft.
Realization of high precision measurement through CEI requires the estimation of carrier differential phase, which could be modeled as the frequency estimation of sinusoidal signals, where many scholars have carried out research [4]-[8]. Such methods can be divided into two categories: non-iterative and iterative ones. Existing non-iterative methods generally rely on the peak DFT coefficient and its neighbors, such as Quinn [7], Candan [5] and Liang [8]. These algorithms, however, suffer from problems like lower SNR threshold or high computational load. Comparatively, the AM iterative methods proposed by Aboutanios and Mulgrew [6] have attracted many attention. The existing original and improved AM algorithms, however, either lack computational efficiency, or could not reach the CRLB. Therefore, a fast and accurate frequency estimation algorithm for CEI is urgently required to enhance the capability of high-quality monitoring and urgent measurement during SSA.
Considering the urgent need of CEI in monitoring and measurement during SSA and the imperfection of the above algorithms, a high precision frequency estimation algorithm is proposed in this letter. In Sect. 2, the nominal CEI in SSA problem is formulated as a frequency estimation problem for sinusoidal signals. On this basis, improved AM frequency estimation algorithm based on more interpolation information and its iteration form is proposed accordingly in Sects. 3 and 4. Finally, in Sect. 5, the simulation is conducted to demonstrate the performance of the proposed method.
Denotion \({\delta _0}\): the original estimate value; \({\delta _{\rm{1}}}\): estimate value after first iteration; \({\delta _p}\): estimate value after \(p\)th iteration.
2. System Model
The basic idea of CEI carrier differential phase estimation is to extract the same side tone \({s_1}(t)\) and \({s_2}(t)\) from the signals received by the two stations:
| \[\begin{equation*} \left\{ \begin{array}{@{}l@{}} {s_1}(t) = {e^{j2\pi {f_0}t}},{S_1}(f) = FFT[{s_1}(t)]\\ {s_2}(t) = {s_1}(t - \tau ) = {e^{j2\pi {f_0}(t - \tau )}},{S_2}(f) = FFT[{s_2}(t)] \end{array} \right. \tag{1} \end{equation*}\] |
The above process could be modeled as a frequency estimation problem for sinusoidal signals as:
| \[\begin{equation*} {x_i}(n) = S(n) + w(n) = A{e^{j{\phi _0}}}{e^{j2\pi {f_0}n\Delta t}} + w(n) \tag{2} \end{equation*}\] |
Where, \({\rm i} = 1,2,\ldots,k\) correspond to stations \(1, 2,\ldots,k\) respectively, \(n = 0,1,\ldots,N - 1\), \(N\) is the sampling points, \(A\) is the signal amplitude, \({\phi _0}\) is the initial phase of the signal, \({f_0}\) is the signal frequency, \(\Delta t\) is the sampling period, and \(w(n)\) is the complex gaussian white noise.
For signal \(X\), its fractional Fourier coefficient is calculated as:
| \[\begin{equation*} {X_{n + p}} = \frac{1}{N}\sum\limits_{k = 0}^{N - 1} x (k){e^{ - j2\pi k\frac{{n + p}}{N}}} \tag{3} \end{equation*}\] |
The Fourier coefficients thus obtained follow a \(N(0,\frac{{{\sigma ^2}}}{N})\) distribution. Let \(m\) be the index returned by the Maximum Bin Search (MBS). The signal frequency \(f\) is given by: \(f = (m + {\delta _0})\frac{{{f_s}}}{N}\). Let \(n = m + l\), and for small x, \({e^x} = 1 + x + O({x^2}) \approx 1 + x\). Ignoring the serial number \(m\), then the fractional Fourier coefficient is: \({X_{l + p}} = {b_0}\frac{{{\delta _0}}}{{{\delta _0} - l - p}} + {W_{n + p}}\), where \({b_0} = - \frac{{{e^{j2\pi {\delta _0}}} + 1}}{{j2\pi {\delta _0}}}\).
The DFT form of noiseless signal is:
| \[\begin{equation*} S(n) = \frac{{A{e^{j\phi }}}}{N}\frac{{{e^{j2\pi N\frac{f}{{{f_s}}}}} - 1}}{{{e^{j2\pi \left( {\frac{f}{{{f_s}}} - \frac{n}{N}} \right)}} - 1}} \tag{4} \end{equation*}\] |
Thus, the residual of MBS is: \(\left| {\Delta f} \right| = \left| {f - m\frac{{{f_s}}}{N}} \right| \le \frac{{{f_s}}}{{2N}}\). And the Mean Square Error (MSE) of frequency estimation \(\mathop f\limits^ \wedge\) is: \({\sigma_f}^2 = {\mathop{\rm var}} [\mathop f\limits^\wedge]=\frac{{{f_s}^2}}{{12{N^2}}}\). Consequently, it can be found from above modulation that the residual of MBS crude search method is \(O({N^{ - 1}})\), while that of asymptotic Cram\(\acute{\mathrm{e}}\)r-Rao Lower Bound (ACRB) is \(O({N^{ - \frac{3}{2}}})\) (which will be analyzed in detail below). Therefore, there is much room for improvement in the residual of frequency estimation, which lays a foundation for the high-precision frequency estimation method of CEI signals based on interpolation.
3. Improved AM Frequency Estimation Algorithm Based on More Interpolation Information
The ratio of the resulting asymptotic variance to the ACRB goes to \(1\) as the number of samples used in the interpolation goes to infinity [9]. Motivated by this, we increased the number of interpolations in the original AM estimator from \(2\) to \(4\) for better estimation performance.
In the absence of noise, consider the following estimator:
| \[\begin{equation*} \psi = \frac{1}{2}\frac{{({X_{0.5}} + {X_q}) + ({X_{ - 0.5}} + {X_{ - q}})}}{{({X_{0.5}} + {X_q}) - ({X_{ - 0.5}} + {X_{ - q}})}} \tag{5} \end{equation*}\] |
Where \({\delta _0} \in [ - 0.5,0.5],q \in [0,1]\). Search for such a \(q\) that satisfies:
| \[\begin{equation*} \frac{{({\delta _0}^2 - {q^2}) + ({\delta _0}^2 - {{0.5}^2})}}{{({\delta _0}^2 - {q^2}) + 2q({\delta _0}^2 - {{0.5}^2})}} \equiv C \tag{6} \end{equation*}\] |
Where \(C\) is a constant. So that we could meet the requirement of \(\psi = k{\delta _0}\) to satisfy the unbiased property of the estimation. Furthermore, the more closer \(C\) is to \(1\), the easier it is to iterate in the following Algorithm 1. To solve the above problems, we selected \({\delta _0} \in [ - 0.5,0.5]\) and \(q \in [ - 1,1]\) with \(0.01\) intervals, and conducted \(1000\) Monte Carlo simulations, and the value of \(C\) varies with \({\delta _0}\) and \(q\). To further study the selection of optimal \(q\), we solved the standard deviation of \(C\) with different \(q\) when \({\delta _0}\) changes, and the above results are shown in Fig. 1.
 |
Fig. 1 (a) Value of \(C\) with the variation of \(q\) and \({\delta _0}\); (b) standard deviation of \(C\) with different \(q\) (\(q \in [ - 1,1]\)) when \({\delta _0}\) changes. |
It is easy to find that when \(q = 0.51\), \(C\) is mostly approaches to \(1\), even with the variation of \({\delta _0}\). What is more, when \(q = 0.51\), considering the overall range of \({\delta_0}\) (\({\delta _0} \in [ - 0.5,0.5]\)), we can get the smallest std of \(C\), which further prove the correctness of our choice. Then the above estimator could be converted into: \(\psi = {\delta _0}\). Therefore, Aalgorithm 1 can be obtained.
 |
As \({\mathop{\rm var}} (\Re \{ A\} ) = {\mathop{\rm var}} (\Re \{ B\} ) = \frac{{{\sigma ^2}}}{N}\), \(|{b_0}{|^2} = \frac{{{{\cos }^2}(\pi {\delta _0})}}{{{{(\pi {\delta _0})}^2}}}\), \({\sigma _{MFFCI}}\) can be given by:
| \[\begin{equation*} \begin{array}{@{}l@{}} {\sigma _{MFFCI}}^2\\ = \frac{{{f_s}^2}}{{2{N^3}\rho }}\frac{{{\pi ^2}{{({\delta _0}^2 - {{0.51}^2})}^2}{{({\delta _0}^2 - {{0.5}^2})}^2}}}{{{{\cos }^2}(\pi {\delta _0}){{[({\delta _0}^2 - {{0.51}^2}) + ({\delta _0}^2 - {{0.5}^2})]}^2}}}(4{\delta _0}^2 + 1) \end{array} \tag{7} \end{equation*}\] |
Where \({f_s}\) is the sampling frequency, \(N\) is the sampling number, and the signal-to-noise ratio (SNR) \(\rho\) (in \(dB\)) of the signal to be estimated is: \(\rho = 10{\log _{10}}(\frac{{{A^2}}}{{{\sigma ^2}}})\). Where \(\rho\) and \(\sigma\) are the amplitudes of signal and noise respectively.
Signal amplitude \(A\) in this letter is set as \(1\). When phase and frequency are unknown, and \(N\) is large enough, the asympototic frequency estimation CRLB(ACRB) of the signal is: \(ACRB = \frac{{6{f_s}^2}}{{{{(2\pi )}^2}\rho {N^3}}}\). Thus, the ratio of the asymptotic variance of the estimator to ACRB, \(R_{ACRB}^{MFFCI}\), is shown as follows:
| \[\begin{equation*} \begin{array}{@{}l@{}} R_{ACRB}^{MFFCI}\\ = \frac{{{\pi ^4}}}{3}\frac{{{{({\delta _0}^2 - {{0.51}^2})}^2}{{({\delta _0}^2 - {{0.5}^2})}^2}}}{{{{\cos }^2}(\pi {\delta _0}){{[({\delta _0}^2 - {{0.51}^2}) + ({\delta _0}^2 - {{0.5}^2})]}^2}}}(4{\delta _0}^2 + 1) \end{array} \tag{8} \end{equation*}\] |
4. Improved Fractional Fourier Coefficient Interpolation Estimation Algorithm Based on Multiple Iterations
The iterative implementation of an estimator aims at improving the performance by essentially making it uniform over the entire bin. However, only those estimators that have a maximum variance at \({\delta _0} = 0\), rather than a minimum one could be improved by iterations. In this section, we consider the iterative implementation of the MFFCI Algorithm 2.
 |
Firstly, consider the estimated interpolation function \(h(\delta)\), which is obtained by Algorithm 1:
| \[\begin{equation*} h(\delta ) = \frac{1}{2}\Re \left\{ {\frac{{({X_{0.5}} + {X_{0.51}}) + ({X_{ - 0.5}} + {X_{ - 0.51}})}}{{({X_{0.5}} + {X_{0.51}}) - ({X_{ - 0.5}} + {X_{ - 0.51}})}}} \right\} \tag{9} \end{equation*}\] |
Substitute the expression for \({X_p},p = \pm 0.5, \pm 0.51\) into \(h(\delta )\) and simplify, we have: \(h({\delta _0}) = [1 + O({N^{ - \frac{1}{2}}}\sqrt {\ln N} )]O({10^{ - 11}})\). As it can be found that: when \(N \to \infty\), \(h({\delta _0}) \to \infty\), for the interpolation function \(h(\delta )\), applying Taylor expansion at \({\delta _0}\), the iterative function \(\psi (\delta ) = \delta + h(\delta )\) can be expanded as:
| \[\begin{equation*} \begin{array}{@{}c@{}} \psi (\delta ) = \delta + h({\delta _0}) + (\delta - {\delta _0})h'({\delta _0})\\ = \delta (1 + h'({\delta _0})) - {\delta _0}h'({\delta _0}) \end{array} \tag{10} \end{equation*}\] |
5. Simulation Results
During these simulation experiments, the hardware is set as follows: CPU: Intel Core i7-10750H 2.60 GHz; RAM 16.0 GB; System type: Windows 64-bit operating system; Matlab: 2019b. The other specific simulation parameter setting is shown in Table 1.
Table 1 Simulation parameter setting. |
 |
As has been discussed in Sects. 1 and 2, the targets of this letter are significant space assets in Earth orbit and Earth-Moon space orbit that could transmit downlink signals, including Communications satellites, manned spacecraft, and deep space probes, which calls for increasingly accurate tracking, measuring and orbiting in SSA missions.
5.1 Performance Analysis and Discussion
Figure 2(a) shows that after two iterations, the estimated asymptotic variance of the proposed algorithm (IMFFCI) can basically achieve uniform convergence, and the ratio of it to ACRB is \(1.00235\) in \({\delta _0} \in [ - 0.5,0.5]\), which is closer to 1 than the current best AM algorithms [6], [10] (FFCI, MSI, MOI, \(1.0147\)), significantly better than Quinn [7] (\(3.264\)), Maclcod [8] (\(1.678\)) and other algorithms whose performance deteriorates after iteration, with a performance improved by \(1.23\%\), \(225.63\%\) and \(67.41\%\) respectively. Although it is higher than that of QSE and HAQSE [10] when \({\delta _0}\) is quite approaching \(0\), it performs better from \({\delta _0} \in [ - 0.35,-0.5]\) and \({\delta _0} \in [ 0.35,0.5]\) (that of QSE and HAQSE will finally increase to \(1.0145\)), demonstrating the performance advantage of our algorithm after iteration. Through iteration, the phenomenon that the asymptotic variance of the proposed algorithm (MFFCI) slightly deteriorates before iteration (when \(\left| {{\delta _0}} \right| > 0.3\)) is improved, and the final frequency estimation achieve a better performance in the whole range of \({\delta _0} \in [ - 0.5,0.5]\).
 |
Fig. 2 (a) Ratio of asymptotic variances of different algorithms to ACRB before and after iteration (\(SNR=0,N=1024\)); (b) comparison of frequency estimation errors of different algorithms under different SNRs before and after iteration (\(N=1024\)). |
As Fig. 2(b) illustrates, before iteration, the estimation error of the proposed algorithm (MFFCI) is improved to a certain extent (only 0.7561 times of traditional AM series algorithms [6], [10] (FFCI, MOI, MSI)) within the range of \(SNR \in [ - 20\,dB,0\,dB]\). After iteration, in the interval of \(SNR \in [ - 14\,dB,0\,dB]\), the estimation error of the proposed algorithm (IMFFCI) decreases significantly, which is basically comparable to CRLB (maintains at 1.236 times) and also slightly better than traditional AM algorithms (only \(0.897\) times of them). Comparatively, when \(SNR \in [ - 20\,dB,- 6\,dB]\), that of RCTSL [9] algorithm is substantially higher than our algorithm, up to four orders of magnitude worse.
The above analysis proves that the proposed algorithms (MFFCI and IMFFCI) achieve small estimation error before and after iteration, and after iteration, the estimation error is closer to CRLB, showing the effectiveness and superiority of the proposed method. Compared with the existing spacecraft tracking and measurement technology such as radar and optic telescope, CEI has the following significant advantages: (i) It only needs downlink signals, which are more user-friendly and is easier for operation and maintenance; (ii) It carries out high-precision monitoring and early warning on space assets, especially when the radar is limited by geographical location or transmitting power, and the astronomical telescope is hindered by atmospheric factors; (iii) Its system can provide strong horizontal constraints and high-precision Angle information of the target. If combined with the distance and velocity information provided by optics and radar, the near-real-time and high-precision orbit determination of non-cooperative targets can be realized, effectively supporting future monitoring and cataloging missions.
5.2 Computational Complexity Analysis
Lower computational complexity is not only preferable, but also crucial for satellite communication and monitor systems. In Table 2, 2th I stands for 2th iterations, \(\times\) and\(+\) stand for complex multiplications and additions under \(N=1024\). Compared with previous algorithms, the proposed algorithms’ multiplication and addition complexity is relatively small. As the total performing time is determined by complex multiplication, it can be found that the proposed algorithms significantly improve the estimation performance without significant increase of the total operation time.
Table 2 Comparison of computational complexity. |
 |
6. Conclusion
Based on more interpolation information of fractional Fourier Coefficients, a MFFCI algorithm for accurate CEI during SSA is proposed in this letter. The proposed algorithms can estimate the frequency of sinusoidal signals with high precision and moderate complexity. Furthermore, the estimation accuracy of MFFCI is improved by iterations and its convergence characteristic is proved. Simulation results demonstrate the effectiveness of the proposed algorithms, indicating a huge potential in futuristic high-precision space assets monitoring of SSA missions.
References
[1] T.J. Muelhaupt, M.E. Sorge, J. Morin, and R.S. Wilson, “Space traffic Management in the new Space era,” Journal of Space Safety Engineering, vol.6, no.2, pp.80-87, 2019.
CrossRef
[2] G.H. McCall and J.H. Darrah, “Space situational awareness: Difficult, expensive and necessary,” Air and Space Power Journal, vol.28, no.6, p.6, 2014.
[3] B.J. Geldzahler, “Coherent uplink arraying techniques for next generation orbital debris, near earth object, and space situational awareness radar systems,” Active and Passive Signatures III, vol.8382, SPIE, 2012.
CrossRef
[4] L.Y. Fang, D.L. Duan, and L.Q. Yang, “A new DFT-based frequency estimator for single-tone complex sinusoidal signals,” Proc. IEEE Military Communications Conference, Orlando, America, pp.1-6, Oct. 2012.
CrossRef
[5] C. Candan, “A method for fine resolution frequency estimation from three DFT samples,” IEEE Signal Process. Lett., vol.18, no.6, pp.351-354, June 2011.
CrossRef
[6] E. Aboutanios and B. Mulgrew, “Iterative frequency estimation by interpolation on Fourier coefficients,” IEEE Trans. Signal Process., vol.53, no.4, pp.1237-1242, April 2005.
CrossRef
[7] B. Quinn, “Estimating frequency by interpolation using Fourier coefficients,” IEEE Trans. Signal Process., vol.42, no.5, pp.1264-1268, May 1994.
CrossRef
[8] X. Liang, A. Liu, X. Pan, Q. Zhang, and F. Chen, “A new and accurate estimator with analytical expression for frequency estimation,” IEEE Commun. Lett., vol.20, no.1, pp.105-108, Jan. 2016
CrossRef
[9] C. Yang and G. Wei, “A noniterative frequency estimator with rational combination of three spectrum lines,” IEEE Trans. Signal Process., vol.59, no.10, pp.5065-5070, 2011.
CrossRef
[10] A. Serbes, “Fast and efficient sinusoidal frequency estimation by using the DFT coefficients,” IEEE Trans. Commun., vol.67, no.3, pp.2333-2342, 2018.
CrossRef
