1. Introduction
Direct sequence spread spectrum (DSSS) systems are extensively used in the field of satellite telemetry, tracking and command (TT&C), navigation and communication [1], [2], due to its robust anti-interference ability, convenience of multiple access, and simultaneous communication and ranging capabilities. Band-limited DSSS systems use band-limited pulse waveforms such as root raised cosine (RRC), Gaussian, and other optimized designed pulse waveforms, and have advantages in spectrum efficiency, spectrum leakage, and ranging accuracy, compared with square pulse waveform systems.
Code tracking is the fine time synchronization process. In band-limited DSSS systems carrying critical applications, the stability is equally important in addition to the error variance performance in code tracking. In order to achieve higher stability, i.e., lower probability of loss of lock, when the error variance performance is the same, the code tracking method needs to have a larger maximum and linear code tracking ranges. However, the existing code tracking methods based on correlation delay-locked loop (DLL) have relatively small maximum and linear code tracking ranges, which are difficult to meet the high stability requirement for high demand band-limited DSSS systems. Therefore, the aim of this paper is to propose a high stability code tracking method with a larger maximum and linear code tracking ranges under the condition of constant error variance performance and appropriate computational complexity.
The code tracking problem has been extensively studied for square pulse waveform DSSS systems over the past few decades [3]-[5]. The first all-digital code tracking method applied to band-limited DSSS systems was proposed in [6], which used digital matched filter, interpolator and a \(2\) times decimator to generate the required data of forward and backward branches. To reduce the complexity of code tracking, [7], [8] use parallel matched filter to improve the matched filter part, [9], [10] use fast Fourier transform (FFT) to improve both the matched filter and interpolator parts and [11] realizes the matched filter, interpolator, and code despreading simultaneously in the frequency domain, achieving further complexity reduction. In addition, many methods based on Kalman filter [12] and adaptive loop bandwidth methods [13] were proposed to further improve the error variance performance of code tracking by improving the loop filter part. To cope with complex channel conditions, researchers have also improved the performance of code tracking in multipath environment [14] and narrowband interference [15], etc. Although these studies have improved the code tracking under different conditions, they are still based on the correlation DLL, which limits the maximum and linear code tracking ranges and the probability of loss of lock, as well as the stability. Simultaneously using the outputs of multiple correlators to obtain the error estimation can extend the maximum code tracking range, but additional noise accumulation inevitably degrades the error variance performance. Therefore, [16], [17] selects the optimal two out of multiple correlation outputs for time-delay variance estimation, i.e., the selection DLL method, but the presence of selection leads to having a severe loss of error variance performance under low signal-to-noise ratios (SNRs).
This paper presents a code tracking method for band-limited DSSS systems in order to achieve higher stability with constant error variance performance. The time-delay error adjustment and correlation are implemented in frequency domain to obtain the frequency domain vectors. In order to reduce the computational complexity of the subsequent processing, we realize the dimension reduction of the vectors by frequency domain integration and dump. Afterwards, in order to make full use of the frequency domain vectors to achieve a larger maximization and linear code tracking range, as well as to improve the stability, we use the modified root-MUSIC subspace algorithm to estimate the time-delay error. Finally the estimated time-delay error is passed through the loop to realize the code tracking. The frequency domain integration and dump is to ensure that the proposed method has appropriate computational complexity, which exploits the property that the autocorrelation of the band-limited DSSS signals aggregates at zero point during the code tracking process, and realizes the dimension reduction of the vectors in frequency domain in a way similar to low-pass filtering, thus the subsequent time-delay error estimation can be implemented with small length frequency domain vectors. We also analyze the code tracking range, steady-state error variance performance and computational complexity of the proposed method, and simulate the proposed method focusing on the first two.
The rest of the paper is organized as follows. We first model the received signal in Sect. 2. Then, in Sect. 3, we propose the high stability code tracking method. In Sect. 4, Sect. 5 and Sect. 6, we analyze the code tracking range, steady-state error variance performance and computational complexity of the proposed method, respectively. Simulation results are shown in Sect. 7. Finally, Sect. 8 gives the conclusion. The major variables adopted in this paper is listed in Table 1 for ease of reference.
Table 1 Variable list. |
 |
2. Signal Model
We adopt DSSS binary phase-shift keying (BPSK) modulation with RRC pulse shaping in this paper, and the transmitted signal can be modeled as [6]
| \[\begin{align} s(t) = \sum\limits_{m = - \infty }^\infty {{a_{\lfloor {m{T_c}/T} \rfloor }}{c_m}} {p_T}(t - m{T_c}), \tag{1} \end{align}\] |
where \(\lfloor \bullet \rfloor\) is the rounding toward negative infinity operation, \(a_i\) is the data symbol with a period of \(T\), \({c_m} = \pm 1\) is the spread spectrum code with chip period \(T_c\) and cycle length \(N\), and \({p_T}(t)\) is the pulse waveform of the transmit signal with its Fourier transform \({P_T}(f) = {T_c}\sqrt {{P_N}(f)}\), where \({P_N}(f)\) is the frequency response of Nyquist RC filter.
We assume that the spread spectrum code period is equal to the data symbol period, i.e., \(T = NT_c\) in this paper, and the method proposed here can be extended to the cases easily where \(T\) is an integer multiple of the code period.
In code tracking, we assume that the receiver has completed acquisition and the received signal is still affected by carrier phase, time delay, and noise. As a result, the complex baseband representation of the received signal can be expressed as
| \[\begin{align} r(t)&=e^{j\theta}s(t-\tau)+w(t) \notag \\ &=e^{j\theta}\sum\limits_{m = - \infty }^\infty {{a_{\lfloor {m{T_c}/T} \rfloor }}{c_m}} {p_T}(t - m{T_c}-\tau)+w(t), \tag{2} \end{align}\] |
where \(\theta\) is the carrier phase, uniformly distributed in \([ - \pi, \pi)\), \(\tau\) is the initial time delay, \(w(t)\) is the complex Gaussian white noise with power spectral density (PSD) \({N_0}/P\), where \(P\) is the power of the intermediate frequency signal.
Code tracking is the fine time synchronization process. Before this, the receiver has already completed the rough estimation of \(\tau\) through acquisition, and the initial time delay \(\tau\) is generally much smaller than the spread spectrum code period, i.e., \(\tau \ll NT_c\).
The received signal after sampling can be expressed as
| \[\begin{align} r_{\ell}=e^{j \theta} \sum_{m=-\infty}^{\infty} a_{\left\lfloor m T_c / T\right\rfloor} c_m p_T\left[\left(\ell-d \right)T_s-m T_c\right]+w_\ell, \tag{3} \end{align}\] |
where \(r_{\ell}=r(\ell T_s)\), \(w_\ell=w(\ell T_s)\), \(T_s\) is the sampling period, and \(d=\tau/T_s\) is the normalized initial time delay. We make \({N_s}\) be the number of sampling points of a single data period, so that \({T=NT_c=N_sT_s}\).
3. Proposed Method
The proposed code tracking method accomplishes time-delay adjustment and correlation in the frequency domain, dimension reduction of vectors by frequency domain integration and dump, and time-delay error estimation by a modified subspace algorithm, after which code tracking is achieved through the loop. It is able to fully utilize the generated frequency domain vectors to achieve time-delay error estimation instead of the two forward and backward correlation values in the correlation DLL. The structure of the proposed method is illustrated in Fig. 1.
 |
Fig. 1 Structure of the proposed code tracking method. |
The proposed method contains three main components. The first is the frequency domain vector generation, which completes time-delay adjustment of the received signal based on the previously feedback time delay, realizes correlation with the local waveform, and vector dimension reduction through frequency domain integration and dump. The second is the time-delay error estimation, which estimates the time-delay error using the dimension-reduced vector through the modified root-MUSIC subspace algorithm. The third is the time delay generation, which involves a multiplier and an integrator, used to obtain the estimated value of time delay for the next code tracking period.
In addition, following the frequency domain vector generation, an accumulator can be used to generate the signal for carrier recovery and data detection.
3.1 Frequency Domain Vector Generation
The frequency domain vector generation component primarily completes time-delay adjustment, correlation with the local waveform, and vector dimension reduction.
In order to make the data symbols almost constant within a single processing cycle, i.e., data period, the time-delay adjustment is divided into two parts: integer delay adjustment and fractional delay adjustment. Integer delay adjustment is completed by shifting the time domain signal after sampling directly, while fractional delay adjustment is realized by conjugate multiplication of signal and the fractional delay filter in the frequency domain.
To achieve correlation between the received signal and the local waveform, we perform a conjugate multiplication in the frequency domain. The local waveform includes information about the spread spectrum code and the pulse waveform and can be generated in advance and remain unchanged. By employing correlation, we complete the despreading and matched filtering of the received signal simultaneously, thereby enabling the accumulation of energy of the useful signal in the time domain.
The correlation output of the band-limited DSSS signal is concentrated near zero in the time domain. Similar to low-pass filtering, we can use integration and dump in the frequency domain to filter out the invalid data in time domain. The integration and dump greatly reduces the length of the frequency domain vectors, allowing the subsequent estimation of the time-delay error realized efficiently, making the computational complexity of the proposed code tracking method appropriate.
We assume that in the \({n}\)-th data period, the feedback normalized time delay is \({\hat d_n}\), with its initial value \({\hat d_{0}}=0\). We adjust the integer delay \({\lfloor {{{\hat d}_n}} \rfloor }\) by shifting the received signal in the time domain. The fractional delay \({{\hat d}_n-\lfloor {{{\hat d}_n}} \rfloor }\) is adjusted by frequency domain fractional delay filter. Since time domain shifting does not affect the statistical properties of the noise, we can assume the noise term unchanged. Since both the time delay and the pulse waveform period are much smaller than the data symbol period, \(\tau \ll T\) and \(T_c=T/N\), for simplicity we ignore the effect between different data symbols. The received signal of the \(n\)-th data period, after integer delay adjustment and serial-to-parallel conversion, is given by
| \[\begin{align} {r_{\ell,n }} &= {a_n}{e^{j\theta }}\sum\limits_{m = 0}^{N - 1} {{c_m}} {p_T}\left[ {\left( {\ell - d + \lfloor {{{\hat d}_n}} \rfloor } \right){T_s} - m{T_c}} \right] + w_{\ell,n}\notag \\ &\ {\rm{for}}\quad \ell = 0,1, \cdots ,{N_s} - 1. \tag{4} \end{align}\] |
The received signal can be transformed into frequency domain after \({N_s}\)-pt FFT. Subsequently, the frequency domain representation of (4) can be expressed as
| \[\begin{align} \label{R_{k,n}} {R_{k,n}} &= {a_n}{e^{j\theta }}H_k^{FD}\left( {d - \lfloor {{{\hat d}_n}} \rfloor } \right)G_k^{TX} + {W_{k,n}}\notag \\ &\ {\rm{for}} \quad k = 0,1, \cdots ,{N_s} - 1, \tag{5} \end{align}\] |
where \(H_k^{FD}\left( d \right)\) is the frequency domain fractional delay filter with length \(N_s\) [11], [18], [19], which can be used to model the time delay in frequency domain, and its expression is
| \[\begin{align} H_k^{FD}(d) = \left\{ {\begin{array}{ll@{}} {1,}&{k = 0}\\ {{e^{ - j\frac{{2\pi }}{{{N_s}}}kd}},}&{k = 1,2, \cdots ,{N_s}/2 - 1}\\ {cos(\pi d),}&{k = {N_s}/2}\\ {{e^{j\frac{{2\pi }}{{{N_s}}}({N_s} - kd)}},}&{k = {L_s}/2 + 1, \cdots ,{N_s} - 1}, \end{array}} \right. \tag{6} \end{align}\] |
and where \(G_k^{TX}\) is the frequency domain waveform of the transmitted signal, which contains the information of the spread spectrum code and the pulse waveform, and can be expressed as
| \[\begin{align} G_k^{TX} = {{\mathop{\rm FFT}\nolimits} _{{N_s}}}\left[ {\sum\limits_{m = 0}^{N - 1} {{c_m}} {p_T}(\ell {T_s} - m{T_c})} \right]. \tag{7} \end{align}\] |
After the integer delay adjustment has been completed, the fractional delay is given by \({d - \lfloor {\hat{d}_n} \rfloor}\). This is also reflected in (5). To adjust the fractional time delay, we perform a conjugate multiplication of \(R_{k,n}\) with a frequency domain fractional delay filter, denoted by \(H_k^{FD}({\hat{d}_n} - \lfloor{\hat{d}_n}\rfloor)\). The signal after fractional delay adjustment can be expressed as
| \[\begin{align} \label{R_{k,n}^{(1)}} R_{k,n}^{\alpha } &= {\left[ {H_k^{FD}({{\hat d}_n} - \lfloor {{{\hat d}_n}} \rfloor )} \right]^*}{R_{k,n}} \notag \\ &= {a_n}{e^{j\theta }}H_k^{FD}\left( {{\varepsilon _n}} \right)G_k^{TX} + {W_{k,n}}, \tag{8} \end{align}\] |
where \({\varepsilon_n} = d - {\hat{d}_n}\) represents the normalized time-delay error. The frequency domain fractional delay filter almost does not affect the power spectrum of the noise, but only causes a phase change. Therefore, we can assume that in (8), the noise term remains unchanged.
After the time delay adjustment, \(R_{k,n}^{\alpha }\) still contains information about the spread spectrum code and pulse waveform in \(G_k^{TX}\). To accumulate the energy of the useful signal, we perform the correlation by conjugate multiplication of \(R_{k,n}^{\alpha }\) with a local frequency domain waveform, denoted by \(G_k^{RX}\). This operation effectively performs the code despreading and matched filtering. The local frequency domain waveform is given by
| \[\begin{align} G_k^{RX} = {{\mathop{\rm FFT}\nolimits} _{{N_s}}}\left[ {\frac{{{T_s}}}{N}\sum\limits_{m = 0}^{N - 1} {{c_m}} {p_R}(\ell {T_s} - m{T_c})} \right], \tag{9} \end{align}\] |
where \({p_R}(t)\) is the local pulse waveform, and its Fourier transformation is \({P_R}(f) = \sqrt {{P_N}(f)}\). After completing the correlation, the result can be expressed as
| \[\begin{align} R_{k,n}^{\beta} &= {\left( {G_k^{RX}} \right)^*}R_{k,n}^{\alpha } \notag \\ &= {a_n}{e^{j\theta }}H_k^{FD}\left( {{\varepsilon _n}} \right){G_k} + {{\tilde W}_{k,n}}, \tag{10} \end{align}\] |
where \({G_k} = {\left( {G_k^{RX}} \right)^*}G_k^{TX}\), and \({\tilde W}_{k,n}=\left( {G_k^{RX}} \right)^*W_{k,n}\). Through the derivation presented in Appendix A, \({G_k}\) can be expressed as
| \[\begin{align} G_k=\frac{T_c}{T_s} P_N\left(\frac{k}{N_s T_s}\right). \tag{11} \end{align}\] |
We define \({{\mathbf{\tilde{W}}}_n} = \left[ \begin{array}{cccc} {\tilde{W}_{0,n}}&{\tilde{W}_{1,n}}& \cdots &{\tilde{W}_{N_s-1,n}} \end{array} \right]^T\), and the noise covariance matrix \(\overline {{{{\bf{\tilde W}}}_n}{\bf{\tilde W}}_n^H}\) is given in Appendix B.
When correlation is completed, the time domain signal associated with \(R_{k,n}^{\beta}\) becomes concentrated around the zero point, and its offset is related to the time-delay error \(\varepsilon_n\). If the frequency domain form of a time domain signal has non-zero values only in the low frequency part of the signal, we can realize a reduction in the data rate of the time domain signal by using time domain integration and dump, i.e., low-pass filtering and extraction operations. Because the time domain and frequency domain are two symmetrical representations, we can also use frequency domain integration and dump to reduce the length of the frequency domain signal, i.e., to remove the invalid time domain data. The result can be expressed as
| \[\begin{align} \label{X_{k,n}} {X_{k,n}} &= \frac{1}{{{N_I}}}\sum\limits_{m = k{N_I} - {N_I}/2}^{k{N_I} + {N_I}/2 - 1} {R_{m,n}^{\beta }} \notag \\ &\approx {a_n}{e^{j\theta }}H_k^{FD}\left( {{\varepsilon _n}} \right){{\tilde G}_k} + {V_{k}} \notag \\ &\ {\rm{for}}\quad k = 0,1, \cdots ,{N_b} - 1, \tag{12} \end{align}\] |
where \(N_I\) is the integration and dump number, and \(N_b\) is the length of output signal, which satisfies \(N_b=N_s/N_I\). \(H_k^{FD}\left( {{\varepsilon _n}} \right)\) is the frequency domain fractional delay filter with length \(N_b\) here. The RC pulse function’s attenuation characteristic and (11)–(12) allow us to express the signal waveform approximately as
| \[\begin{align} \tilde{G}_k&\approx\frac{T_c}{T_s} P_N\left(\frac{k}{N_b T_s}\right). \tag{13} \end{align}\] |
We define \({{\mathbf{V}}_n} = \left[ \begin{array}{cccc} {V_0}&{V_1}& \cdots &{V_{N_b-1}} \end{array} \right]^T\), and the noise covariance marix \(\overline{{\mathbf{V}}_n{\mathbf{V}}_n^H}\) is also given in Appendix B.
We rewrite (12) to simplify the expression as
| \[\begin{align} {X_{k,n}}&= {A_k}({\varepsilon _n}){F_n} + {V_{k}} \notag \\ &\ {\rm{for}}\quad k = 0,1, \cdots ,{N_b} - 1, \tag{14} \end{align}\] |
where
| \[\begin{align} & F_n=a_n e^{j\theta}, \notag \\ & {A_k}({\varepsilon _n})=H_k^{FD}\left( {{\varepsilon _n}} \right) {\tilde G}_k, \tag{15} \end{align}\] |
and the vectorial form is
| \[\begin{align} {{\mathbf{X}}_n} &= {\mathbf{{A}}}({\varepsilon _n}){F_n} + {{\mathbf{V}}_n}, \tag{16} \end{align}\] |
where
| \[\begin{align} & {{\mathbf{X}}_n} = \left[ \begin{array}{cccc} {{X_{0,n}}}&{{X_{1,n}}}& \cdots &{{X_{{N_b-1},n}}} \end{array} \right]^T, \notag \\ & \mathbf{{ A}}({\varepsilon _n}) = {\mathop{\rm diag}} \{ {\bf{\tilde G}}\} {{\bf{H}}^{FD}}(\varepsilon ), \notag \\ & \bf{\tilde G}= \left[ \begin{array}{cccc} {\tilde G}_0&{\tilde G}_1& \cdots &{\tilde G}_{N_b-1} \end{array} \right]^T, \notag \\ & {{\bf{H}}^{FD}}(\varepsilon )= \left[ \begin{array}{cccc} H_0^{FD}\left( {{\varepsilon}} \right)&H_1^{FD}\left( {{\varepsilon}} \right)& \cdots &H_{N_b-1}^{FD}\left( {{\varepsilon}} \right) \end{array} \right]^T. \tag{17} \end{align}\] |
The frequency domain vector \({\mathbf{X}}_n\) can be used for two purposes. First, it can generate the signal for carrier recovery and data detection through an \(N_b\) point accumulator. Second, it can be used to estimate the time-delay error \(\varepsilon _n\).
3.2 Time-Delay Error Estimation
In order to fully utilize the frequency domain vectors, we use a modified root-MUSIC subspace algorithm to estimate the time-delay error. The parameter estimation of the root-MUSIC algorithm can be divided into two parts, one is the covariance matrix estimation, and the other is the eigendecomposition and polynomial root finding.
In covariance matrix estimation, the root-MUSIC algorithm uses the method of averaging the sample covariance matrices, which is not suitable for code tracking with feedback. To address this limitation, we propose a new method for covariance matrix estimation by using RC (resistor-capacitor) integral filter, since it can be used in the presence of feedback and the sum of coefficients of the filter is 1. The system function of the RC (resistor-capacitor) integral filter can be expressed as
| \[\begin{align} H(s) = \frac{1}{{1 + s{\tau _1}}}, \tag{18} \end{align}\] |
where \(\tau _1\) is the time constant. By replacing the differential operator \(s\) with \((1-z^{-1})/T\), we can obtain the system function in the \(z\) domain as
| \[\begin{align} H(z) = \frac{T}{{T + {\tau _1} - {\tau _1}{z^{ - 1}}}}. \tag{19} \end{align}\] |
Then at the \(n\)-th data period, the estimated covariance matrix can be expressed as
| \[\begin{align} {{\mathbf{\hat R}}_n} = \frac{{{\tau _1}}}{{T + {\tau _1}}}{{\mathbf{\hat R}}_{n - 1}} + \frac{T}{{T + {\tau _1}}}{{\mathbf{X}}_n}{\mathbf{X}}_n^H, \tag{20} \end{align}\] |
where the initial value \({{\mathbf{\hat R}}_{-1}} = 0\). After eigendecomposition, (20) can be expressed as
| \[\begin{align} {{\bf{\hat R}}_n} = {{\bf{\hat { E}}}_S}{{{\hat \Lambda }}_S}{\bf{\hat E}}_S^H + {{\bf{\hat E}}_N}{{\bf{\hat \Lambda }}_N}{\bf{\hat E}}_N^H, \tag{21} \end{align}\] |
where \({{\bf{\hat E}}_S}\) is the estimated signal subspace, which is a vector in this paper, and \({{{\hat \Lambda }}_S}\) is the corresponding eigenvalue. \({{\bf{\hat E}}_N}\) is the estimated noise subspace, \({{\bf{\hat \Lambda }}_N}\) is the corresponding eigenvalue matrix. The signal subspace \({{\bf{\hat E}}_S}\) is eigenvector that corresponds to the largest eigenvalue of the covariance matrix. Conversely, the noise subspace \({{\bf{\hat E}}_N}\) is composed of the eigenvectors that correspond to the remaining eigenvalues and \({{\bf{\hat E}}_N}\) is a matrix of \(N_b\times (N_{b}-1)\).
We can use the estimated noise subspace matrix \({{\bf{\hat E}}_N}\) to construct the parameter spectrum of the time-delay error estimation, that is
| \[\begin{align} S(\varepsilon ) = \frac{1}{{{{\bf{{ A}}}^H}(\varepsilon ){{{\bf{\hat E}}}_N}{\bf{\hat E}}_N^H{\bf{{ A}}}(\varepsilon )}}. \tag{22} \end{align}\] |
Since the signal subspace and the noise subspace are nearly orthogonal, the parameter spectrum, denoted by \(S(\varepsilon)\), attains its maximum value when the time-delay error \(\varepsilon\) is optimally estimated.
In order to skip the search process and improve performance at the same time, we estimate the time-delay error by finding the root of \(S^{-1}(\varepsilon )\). By substituting (13), (15), \({S^{ - 1}}(\varepsilon )\) can be rewritten as
| \[\begin{align} {S^{ - 1}}(\varepsilon ) &= {{\bf{{ A}}}^H}(\varepsilon ){{{\bf{\hat E}}}_N}{\bf{\hat E}}_N^H{\bf{{ A}}}(\varepsilon ) \notag \\ &= {\left[ {{{\bf{H}}^{FD}}}(\varepsilon ) \right]^H}{{\mathop{ \rm diag}} ^H}\{ {\bf{\tilde G}}\} {{{\bf{\hat E}}}_N}{\bf{\hat E}}_N^H{\mathop{\rm diag}} \{ {\bf{\tilde G}}\} {{\bf{H}}^{FD}}(\varepsilon ), \tag{23} \end{align}\] |
According to (17) and (6), \(\mathbf{H}^{FD}(\varepsilon)\) can be decomposed as
| \[\begin{align} {{\mathbf{H}}^{FD}}(\varepsilon ) = {\mathbf{CE}}(\varepsilon ), \tag{24} \end{align}\] |
where the phase vector
| \[\begin{align} {\mathbf{E}}(\varepsilon ) = \left[ {\begin{array}{*{20}{c}} {{e^{j\frac{{2\pi }}{{{N_b}}}\frac{{{N_b}}}{2}\varepsilon }}}&{{e^{j\frac{{2\pi }}{{{N_b}}}\left( {\frac{{{N_b}}}{2} - 1} \right)\varepsilon }}}& \cdots &{{e^{ - j\frac{{2\pi }}{{{N_b}}}\frac{{{N_b}}}{2}\varepsilon }}} \end{array}} \right]^T, \tag{25} \end{align}\] |
and the transformation matrix
| \[\begin{align} {\bf{C}} = {\left[ {\begin{array}{ccccccccc} {}&{}&{}&{}&1&{}&{}&{}&{}\\ {}&{}&{}&{}&{}&1&{}&{}&{}\\ {}&{}&{}&{}&{}&{}& \ddots &{}&{}\\ {}&{}&{}&{}&{}&{}&{}&1&{}\\ {1/2}&{}&{}&{}&{}&{}&{}&{}&{1/2}\\ {}&1&{}&{}&{}&{}&{}&{}&{}\\ {}&{}& \ddots &{}&{}&{}&{}&{}&{}\\ {}&{}&{}&1&{}&{}&{}&{}&{} \end{array}} \right]_{{N_b} \times ({N_b} + 1).}} \tag{26} \end{align}\] |
Thus (23) can be further expressed as
| \[\begin{align} {S^{ - 1}}(\varepsilon ) = {\left[ {{\mathbf{E}}(\varepsilon )} \right]^H}{{\mathbf{C}}^H}{\text{dia}}{{\text{g}}^H}\{ {\mathbf{\tilde G}}\} {{\mathbf{\hat E}}_N}{\mathbf{\hat E}}_N^H{\text{diag}}\{ {\mathbf{\tilde G}}\} {\mathbf{CE}}(\varepsilon ). \tag{27} \end{align}\] |
We define
| \[\begin{align} {\bf{B}} = {{\bf{C}}^H}{{\mathop{\rm diag}} ^H}\{ {\bf{\tilde G}}\} {{\bf{\hat E}}_N}{\bf{\hat E}}_N^H{\mathop{\rm diag}} \{ {\bf{\tilde G}}\} {\bf{C}}, \tag{28} \end{align}\] |
and let \(B_{-N_b/2,-N_b/2}\) denote the first element in the upper-left corner of the matrix \(\mathbf{B}\). We can rewrite \({S^{ - 1}}(\varepsilon)\) as follows
| \[\begin{align} {S^{ - 1}}(\varepsilon ) &= \sum\limits_{{m_1} = - {N_b}/2}^{{N_b}/2} {\sum\limits_{{m_2} = - {N_b}/2}^{{N_b}/2} {{e^{j\frac{{2\pi }}{{{N_b}}}{m_1}\varepsilon }}} } {B_{m_1,m_2}}{e^{ - j\frac{{2\pi }}{{{N_b}}}{m_2}\varepsilon }}\notag \\ &= \sum\limits_{m = - {N_b}}^{{N_b}} {{b_m}} {e^{j\frac{{2\pi }}{{{N_b}}}m\varepsilon }}. \tag{29} \end{align}\] |
Additionally, \(b_m\) represents the sum of the \(m\)-th diagonal elements in the matrix \(\mathbf{B}\) and is given by the following equation
| \[\begin{align} {b_m} = \sum\limits_{{m_1} - {m_2} = m} B_{m_1,m_2}. \tag{30} \end{align}\] |
According to (29), we can define polynomial
| \[\begin{align} D(z) = \sum\limits_{m = - {N_b}}^{{N_b}} {{b_m}} {z^m}. \tag{31} \end{align}\] |
By finding the root of \(D(z)\) closest to the unit circle, the time-delay error estimation can be obtained. Assume that the root closest to the unit circle is
| \[\begin{align} z=z_1. \tag{32} \end{align}\] |
The time-delay error estimation, corresponding to the \(n\)-th estimated covariance matrix \({{\bf{\hat R}}_n}\), can be expressed as
| \[\begin{align} {\hat \varepsilon _n} = \frac{{{N_b}}}{{2\pi }}\arg ({z_1}). \tag{33} \end{align}\] |
3.3 Time Delay Estimation and Parameter Settings
In this subsection, we obtain the time delay estimation using a loop. The time-delay error estimation \({\hat \varepsilon _n}\) has been obtained, and we use it as input, passing it through a multiplier and an integrator to calculate the time delay estimation, which is then used as the feedback.
The estimated time delay for the \(n+1\) code tracking period, can be expressed as
| \[\begin{align} {\hat d_{n + 1}} = {\hat d_n} + {G_d}{\hat \varepsilon _n}, \tag{34} \end{align}\] |
where \({G_d}\) is the coefficient of the multiplier for two purposes. First, it consists the gain factor of the integrator. Second, it can controls the loop bandwidth.
The damping coefficient affects the oscillation characteristics and response speed of the tracking system, and in the RC (resistor-capacitor) integral filter tracking loop, the damping coefficient \(\zeta=\sqrt{T/G_d \tau_1}/2\) [20]. The loop bandwidth is the bandwidth of the loop’s equivalent ideal narrow-band filter, which affects the loop’s noise resistance. In the RC (resistor-capacitor) integral filter tracking loop, the loop bandwidth \(B_L=G_d/4T\) [20]. In order to balance the oscillation characteristics and response speed of the tracking system, we set the damping coefficient to a typical value, i.e., \({\zeta}=0.707\), and keep the tracking loop in a critical damping state. Based on the known damping coefficient \(\zeta\) and the required loop bandwidth \(B_L\), we can calculate the multiplication factor \(G_d\) and the time constant \(\tau_1\) of the RC (resistor-capacitor) integral filter required for loop operation as
| \[\begin{align} G_d&=4TB_L, \notag \\ \tau_1&=\frac{1}{8B_L}. \tag{35} \end{align}\] |
In the proposed band-limited DSSS code tracking method, the setting of the sampling period \(T_s\) needs to satisfy the Nyquist sampling theorem, and it is best to meet the requirements of efficient FFT at the same time, that is, the number of sampling points in a single data period is a power of \(2\). When the spread spectrum code adopts the M code or the Gold code, its length satisfies \(N=2^r-1\), where \(r\) is the number of the generation registers. The number of sampling points in a single data period and the sampling period can be set as
| \[\begin{align} &N_s=2(N+1), \notag \\ &T_s= \frac{T}{N_s }=\frac{NT_c}{2(N+1) }. \tag{36} \end{align}\] |
In contrast, the correlation DLL code tracking method proposed in [6] uses \(N_s=2N\) so that the forward and backward branches can be generated by a decimation of \(2\).
4. Code Tracking Range Analysis
The code tracking curve refers to the relationship between the time-delay error output of the code tracking method and the true time-delay error of the input signal with open loop and noiseless conditions, and its two mian observations are the maximum and linear code tracking ranges. The linear code tracking range is the error range in which the error output is linear or nearly linear with the actual error, and the maximum code tracking range is the error range in which the positive and negative signs of the error output are the same as the actual error. A larger maximum and linear code tracking ranges imply a lower probability of loss of lock i.e. higher stability when the steady-state error variance performance is the same.
In this section, we analyze the maximum and linear code tracking ranges of the proposed method and compare them with the correlation DLL and selection DLL methods [6], [17]. The comparison is done using a combination of theoretical and numerical analyses, where the selection DLL is only analyzed numerically due to the complexity of the theoretical analysis.
The proposed code tracking method uses subspace algorithm to estimate the time-delay error in frequency domain, and its maximum code tracking range is determined by the length of frequency domain vectors. At the same time, since the used subspace algorithm is asymptotically unbiased, the linear code tracking range can be considered equivalent to the maximum code tracking range. Referring to (16) and (17), the length of the frequency domain vector \({\mathbf{X}}_n\) is \(N_b\). As there are two time-delay error directions, positive and negative, the normalized maximum and linear code tracking ranges of the proposed method is
| \[\begin{align} \varepsilon_P&=\frac{N_b}{2}. \tag{37} \end{align}\] |
Therefore, the code tracking curve of the proposed method can be expressed as
| \[\begin{align} {\eta _P}(\varepsilon) = \left\{ {\begin{array}{*{20}{c}} {\varepsilon ,}&{ - {\varepsilon_P} \le \varepsilon \le {\varepsilon_P}}\\ {0,}&{{\rm{else}}.} \end{array}} \right. \tag{38} \end{align}\] |
The code tracking curve of the correlation DLL method is given in [6], which can be expressed as
| \[\begin{align} \eta_D(\varepsilon)=\frac{T_c}{A} \left[p^2\left(\varepsilon T_s -\frac{T_c}{2}\right)-p^2\left(\varepsilon T_s+\frac{T_c}{2}\right)\right], \tag{39} \end{align}\] |
where \(p(t)\) is the RC pulse waveform with frequency response \(P(f)=T_c P_N(f)\), and \(A\) is a known scale factor.
The code tracking curves for the proposed code tracing method, the correlation DLL method, and the selection DLL method are shown in Fig. 2, where the solid lines are the theoretical results, the dashed lines are the numerical results, and \(N_{corr}\) is the number of correlation outputs used in the selection DLL method. The signal pulse is a root raised cosine waveform with the roll-off factor \(\alpha=0.25\). The selection algorithm in the selection DLL method determines first the branch with maximal energy and then its neighbor with highest energy. We can discover that the maximum and linear code tracking ranges of the proposed method and the selection DLL method are larger than that of the correlation DLL method, which means that the proposed method and the selective DLL method have lower probability of loss of lock and higher stability under the same steady-state error variance performance. The figure only shows the code tracking curve of the proposed method when \(N_b=8\). When \(N_b\) takes a larger value, such as \(16\) or \(32\), the code tracking range will be increase. As \(N_{corr}\) becomes larger, the code tracking range of the selection DLL method can likewise be further increased.
 |
Fig. 2 Code tracking curves of the proposed method, the correlation DLL method and the selection DLL method. |
5. Steady-State Error Variance Analysis
In this section, we first analyze the output of time-delay error estimation and then perform loop analysis to obtain the expression of the steady-state time-delay error variance \(\overline{\vert\varepsilon_n\vert^2}\).
The frequency domain vector is (16), which can be rewritten as
| \[\begin{align} {{\mathbf{X}}_n} = {\mathbf{{A}}}({\varepsilon _n}){F_n} + {{\mathbf{V}}_n}. \tag{40} \end{align}\] |
Then the covariance matrix can be expressed as
| \[\begin{align} {{\mathbf{ R}}_n}&=\overline{{\mathbf{X}}_n{\mathbf{X}}_n^H} \notag \\ &=\overline{ \left( {\mathbf{{A}}}({\varepsilon _n}){F_n} + {{\mathbf{V}}_n} \right) \left( {\mathbf{{A}}}({\varepsilon _n}){F_n} + {{\mathbf{V}}_n} \right)^H} \notag \\ &={\mathbf{{A}}}({\varepsilon _n})\overline{F_nF_n^H}{\mathbf{{A}}}^H({\varepsilon _n})+\overline{{{\mathbf{V}}_n}{{\mathbf{V}}^H_n}} \notag \\ &={\mathbf{{A}}}({\varepsilon _n}){\mathbf{{A}}}^H({\varepsilon _n})+\overline{{{\mathbf{V}}_n}{{\mathbf{V}}^H_n}} \notag \\ &={\mathbf{{A}}}({\varepsilon _n}){\mathbf{{A}}}^H({\varepsilon _n})+\mathbf{E}\mathbf{\Lambda}\mathbf{E}^H, \tag{41} \end{align}\] |
where the orthonormal matrix \(\mathbf{E}\) and noise eigenvalue matrix \(\mathbf{\Lambda}\) are
| \[\begin{align} & \mathbf{E}=[\mathbf{S}_1,\mathbf{S}_2,\cdots,\mathbf{S}_{N_b}], \notag \\ & \mathbf{\Lambda}=\mathrm{diag}(\lambda_1,\lambda_2,\cdots,\lambda_{N_b}). \tag{42} \end{align}\] |
The first eigenvector in the orthonormal matrix \(\mathbf{E}\) can be defined as
| \[\begin{align} \mathbf{S}_1&=\frac{{\mathbf{{A}}}({\varepsilon _n})}{\Vert {\mathbf{{A}}}({\varepsilon _n}) \Vert_2}. \tag{43} \end{align}\] |
Thus, the first noise eigenvalue can be expressed as
| \[\begin{align} \lambda_1&=\overline{\mathbf{S}_1^H {\mathbf{V}}_n \left(\mathbf{S}_1^H {\mathbf{V}}_n\right)^H} \notag \\ &=\mathbf{S}_1^H \overline{\mathbf{V}_n \mathbf{V}_n^H} \mathbf{S}_1. \tag{44} \end{align}\] |
The sum of residual noise eigenvalues can be expressed based on the matrix analysis theory as
| \[\begin{align} \sum_{i=2}^{N_b} \lambda_i=\mathrm{tr}(\overline{\mathbf{V}_n \mathbf{V}_n^H})-\lambda_1, \tag{45} \end{align}\] |
where the specific expression of the noise covariance matrix \(\overline{\mathbf{V}_n \mathbf{V}_n^H}\) is given in Appendix B. The estimation of the covariance matrix is given by (20), which can be rewritten as
| \[\begin{align} {{\mathbf{\hat R}}_n} &= \frac{{{\tau _1}}}{{T + {\tau _1}}}{{\mathbf{\hat R}}_{n - 1}} + \frac{T}{{T + {\tau _1}}}{{\mathbf{X}}_n}{\mathbf{X}}_n^H \notag \\ &=\sum_{i=0}^{n}\frac{T\tau_1^i}{(T+\tau_1)^{i+1}}\mathbf{X}_{n-i}\mathbf{X}_{n-i}^H. \tag{46} \end{align}\] |
When the system is in a steady state, it can be assumed that \(\mathbf{X}_n\) has the same statistical characteristics at different times \(n\). Applying the central limit theorem, we can assume that \({{\mathbf{\hat R}}_n}\) obeys a complex Wishart distribution with its degree of freedom \(M\) expressed as
| \[\begin{align} & \frac{1}{M} = \mathop {\lim }\limits_{n \to + \infty } \sum\limits_{i = 0}^n {{{\left[ {\frac{{T\tau _1^i}}{{{{(T + {\tau _1})}^{i + 1}}}}} \right]}^2}} \notag \\ & M=1+\frac{2\tau_1}{T}. \tag{47} \end{align}\] |
Therefore, the error variance in the signal zeros of the modified root-MUSIC algorithm can be expressed as [21]
| \[\begin{align} \overline{\left| \Delta z_n \right|^2}=\frac{(\lambda_1^s+\lambda_1)\sum_{i=1}^{N_b}\lambda_i}{M (\lambda_1^s)^2 } \frac{\left| \mathbf{A}^H(\varepsilon_n)\mathbf{S}_1\right|^2}{\left( \mathbf{S}_1^{'}\right)^H \mathbf{P}_N \mathbf{S}_1^{'}}, \tag{48} \end{align}\] |
where \(\mathbf{S}_1^{'}\) is the derivative of \(\mathbf{S}_1\) with respect to \(2\pi \varepsilon_n/N_b\). Additionally, since the time-delay error \(\varepsilon_n\) is approximately equal to \(0\) in steady state, we can substitute \(\varepsilon_n=0\) into (48) when calculating the theoretical value. In this paper, the phase of root can be expressed with normalized time-delay error \(\varepsilon_n\) as \(2\pi \varepsilon_n/N_b\). Therefore, we can express the variance of the estimated normalized time-delay error using the error variance in signal zeros with \(\hat{\varepsilon}_n=\varepsilon_n+\tilde{N}_n\) as [21]
| \[\begin{align} \overline{\left| \tilde{N}_n\right|^2}=\frac{N_b}{8\pi^2}\overline{\left| \Delta z_n \right|^2}. \tag{49} \end{align}\] |
We can observe that \(\tilde{N}_n\) is not independent at different times \(n\) due to the use of RC (resistor-capacitor) integral filter, which is inconvenient for further processing, but we can model the noise term approximately as white noise passes through an RC (resistor-capacitor) integral filter. The time-delay error estimation can be expressed as
| \[\begin{align} \hat{\varepsilon}_n&=\varepsilon_n+\tilde{N}_n \notag \\ &=\varepsilon_n+N_n\ast h_n, \tag{50} \end{align}\] |
where \(h_n\) represents the impulse response of the RC (resistor-capacitor) integral filter. Then, the digital PSD of \({N}_n\) can be expressed as
| \[\begin{align} S_{N}(e^{j\omega})=\frac{2\pi\overline{\left| \tilde{N}_n\right|^2}}{\int_{-\pi}^{\pi}\left|H(e^{j\omega})\right|^2\mathrm{d}\omega}, \tag{51} \end{align}\] |
where \(H(e^{j\omega})\) is the discrete-time Fourier transform of \(h_n\).
Subsequently, we get the final expression of the steady-state time-delay error variance \(\overline{\vert\varepsilon_n\vert^2}\). Since the time-delay error is given by \(\varepsilon_n=d-\hat{d}_n\). Referring to (34) and (50), we can get the time-delay error as
| \[\begin{align} \varepsilon_{n+1}&=d-\hat{d}_{n+1} \notag \\ &=\varepsilon_{n}-G_d{\hat \varepsilon _n} \notag \\ &=\varepsilon_{n}-G_d \left( \varepsilon _n+N_n\ast h_n\right). \tag{52} \end{align}\] |
In order to get the variance of \(\varepsilon_{n}\), we can represent \(\varepsilon_{n}\) and \({N}_n\) in their respective frequency domain forms as [22]
| \[\begin{align} \varepsilon_n&=\frac{1}{2 \pi} \int_{-\pi}^\pi e^{j n \omega} \mathrm{d} \xi_{\varepsilon}\left(e^{j \omega}\right),\notag \\ \quad {N}_n&=\frac{1}{2 \pi} \int_{-\pi}^\pi e^{j n \omega} \mathrm{d} \xi_{{N}}\left(e^{j \omega}\right), \tag{53} \end{align}\] |
where \(\xi_{\varepsilon}\left(e^{j \omega}\right)\) and \(\xi_{N}\left(e^{j \omega}\right)\) are their respective orthogonal incremental processes. By substituting (53) into (52), when the code tracking is in a steady state, we can get
| \[\begin{align} \varepsilon_n&=\frac{1}{2 \pi} \int_{-\pi}^\pi e^{j n \omega} \mathrm{d} \xi_{\varepsilon}\left(e^{j \omega}\right) \notag \\ &=\frac{1}{2 \pi} \int_{-\pi}^\pi \frac{-G_d H(e^{j\omega})}{e^{j\omega}-1+G_d} e^{j n \omega} \mathrm{d} \xi_{{N}}\left(e^{j \omega}\right). \tag{54} \end{align}\] |
Since \(\overline{\mathrm{d} \xi_{{N}}\left(e^{j \omega_1}\right) \mathrm{d} \xi_{{N}}^*\left(e^{j \omega_2}\right)}=2 \pi \delta_{\omega_1, \omega_2} S_{{N}}\left(e^{j \omega_1}\right) \mathrm{d} \omega_1\), and the time-delay error variance \(\overline{\vert\varepsilon_n\vert^2}=\overline{\varepsilon_n\varepsilon_n^*}\), we can get the expression of the final steady-state time-delay error variance as
| \[\begin{align} \overline{\vert\varepsilon_n\vert^2}&=\frac{S_N(e^{j\omega})}{2\pi}\int_{-\pi}^{\pi}\left|\frac{G_d H(e^{j\omega})}{e^{j\omega}-1+G_d}\right|^2 \mathrm{d} \omega \notag \\ &=\frac{N_b \overline{\left| \Delta z_n \right|^2}}{8\pi^2\int_{-\pi}^{\pi}\left|H(e^{j\omega})\right|^2\mathrm{d}\omega}\int_{-\pi}^{\pi}\left|\frac{G_d H(e^{j\omega})}{e^{j\omega}-1+G_d}\right|^2 \mathrm{d} \omega. \tag{55} \end{align}\] |
The above equation provides a closed-form expression for the steady-state error variance of the proposed code tracking method. By substituting (48) into the above equation, we can calculate the steady-state time-delay error variance theoretically.
6. Computational Complexity Analysis
In this section, we analyze the computational complexity of the proposed code tracking method and compare it with the correlation DLL method and the selection DLL method [6], [17]. Specifically, we use the number of complex multiplications needed to process the received signal over a single data period as an indicator of their computational complexities. Due to the very small computational complexity, the time delay generation part is not accounted for all three methods.
The computational complexity of each part of the proposed method is shown in Table 2. \(N_s\) is the number of sampling points in a single data period, and \({N_b}\) is the length of frequency domain vector after integration and dump. The overall computational complexity of the proposed code tracking method is the sum of each part, which can be expressed as
| \[\begin{align} {C_P} = \frac{{{N_s}}}{2}{\log _2}({N_s}) + 2{N_s} + 10N_b^3 + 4N_b^2 + 2{N_b}. \tag{56} \end{align}\] |
Table 2 Computational complexity of each part of the proposed code tracking method. |
 |
For the correlation DLL method, we use FFT to accelerate the matched filter part by default, and the interpolator adopts the cubic interpolation. Then the computational complexity required by the correlation DLL method can be expressed as
| \[\begin{align} {C_D} = 2{N_s}{\log _2}(2M) + 6{N_s} + 2N, \tag{57} \end{align}\] |
where \(M\) is the length of the matched filter and \(N\) is the length of the spread spectrum code. Similarly, for the selection DLL method, the computational complexity can be expressed as
| \[\begin{align} {C_S} = 2{N_s}{\log _2}(2M) + 6{N_s} + N_{corr}N, \tag{58} \end{align}\] |
where \(N_{corr}\) is the number of correlation outputs used in the selection DLL method.
We can perform a specific comparative analysis of the computational complexity of the three methods by substituting typical values. We set the length of the spread spectrum code \(N=1023\). In the proposed method, we set the number of sampling points in a single data period \({N_s} = 2(N + 1)\). In the correlation DLL method and the selection method we set \(N_s=2N\), the length of matched filter \(M = 16\). At the same time, we set \({N_b}\) and \(N_{corr}\) so that the proposed method and the selected DLL method have approximate maximum and linear code tracking ranges. When the number of frequency domain vector after integration and dump \({N_b} = 8\) in the proposed method and the correlation number \(N_{corr}=4\) in the selection DLL method, the computational complexity ratio of the proposed method to the correlation DLL method is \({C_P}/{C_D} = 0.60\), the ratio of the proposed method to the selection DLL method is \({C_P}/{C_S} = 0.56\). When \({N_b} = 16\) and \(N_{corr}=8\), the computational complexity ratios are \({C_P}/{C_D} = 1.65\) and \({C_P}/{C_S} = 1.40\), respectively. When \({N_b} = 32\) and \(N_{corr}=16\), the ratios are \({C_P}/{C_D} = 9.98\) and \({C_P}/{C_S} = 7.07\), respectively.
We can find out that the proposed method’s computational complexity is comparable to that of the correlation DLL method and the selection DLL method when \({N_b}=8\) or \({N_b}=16\). Additionally, as the value of \({N_b}\) increases, the computational complexity of the proposed method increases significantly. In summary, we recommend using relatively small \(N_b\) when using the proposed method for code tracking in band-limited DSSS systems. If a larger maximum and linear code tracking ranges is required, a larger \(N_b\) can be chosen.
7. Simulations
In this section, we perform three simulations to analyze and verify the proposed code tracking method for band-limited DSSS systems. First, we perform simulation to analyze the advantages of the proposed method in terms of the maximum and linear code tracking ranges. Second, we verify the theoretical steady-state error variance closed-form expression. Finally, we compare the proposed method with the correlation DLL method and the selection DLL method [6], [17], in terms of steady-state error variance performance, using Monte Carlo method.
First, we perform simulation to analyze the advantages of the proposed method in terms of the maximum and linear code tracking ranges. We examine the differences between the proposed method, the correlation DLL method and the selection DLL method in their code tracking results under various initial time delays \(\tau\). We set the length of the spread spectrum code \(N=1023\), chip rate \(1/T_c=10.23\: \mathrm{MHz}\), roll-off factor of the RRC pulse waveform \(\alpha=0.25\), code period equal to data period, i.e., \(T=NT_c=N_sT_s\), loop bandwidth \({B_L} = 10\;{\rm{ Hz}}\), and \(\mathrm{SNR}=-10\;\mathrm{ dB}\). We consider five different initial time delays as \(\tau=T_s\), \(\tau=2T_s\), \(\tau=3T_s\), \(\tau=3.9T_s\), and \(\tau=7.9T_s\) where \(T_s\) is the sampling period of the proposed method. For the proposed method, we set sampling rate \(1/{T_s} = 20.48\: {\rm{ MHz}}\), \(N_b=16\) when \(\tau=7.9T_s\), and \(N_b=8\) for other values of \(\tau\). For the correlation DLL method, we set sampling rate \(1/{T_s} = 20.46\: {\rm{ MHz}}\), length of matched filter \(M=16\), and use cubic interpolater to adjust fractional time delay. For the selection DLL method, we set the correlation number \(N_{corr}=8\) when \(\tau=7.9T_s\), \(N_{corr}=4\) for other values of \(\tau\), other parameters are consistent with the correlation DLL method, and the selection algorithm determines first the branch with maximal energy and then its neighbor with highest energy. The code tracking results are presented in Fig. 3.
 |
Fig. 3 Code tracking results of the three methods under different initial time delays. |
In Fig. 3, the solid lines represent the code tracking results of the proposed method, the dashed lines represent the code tracking results using the correlation DLL method and the dashdotted lines represent the code tracking result of the selection DLL method. The x-axis is the time, united by second, and the y-axis is the estimation of the time delay, united by the sampling period of the proposed method.
For an initial time delay \(\tau=T_s\), all three methods establish a steady state quickly, and at this time, \(\tau\) falls within the approximate linear code tracking range of all three methods. For an initial time delay \(\tau=2T_s\), the proposed method and the selection DLL method can establish a steady state faster. At this time, \(\tau\) falls within the nonlinear code tracking range of the correlation DLL method, but falls within the linear code tracking range of the other two methods still. When we set the initial time delay to be \(\tau=3T_s\), \(\tau=3.9T_s\) and \(\tau=7.9T_s\), the correlation DLL method fails to track as the three initial time delays exceed its maximum tracking range, but the proposed method and the selection DLL method can track effectively still. Therefore, the proposed method and the selection DLL method both have larger maximum and linear code tracking ranges compared to the correlation DLL method, which means lower probability of loss of lock and higher stability for the same steady-state error variance performance.
Second, to verify the accuracy of the theoretical steady-state error variance closed-form expression of the proposed method in Sect. 5, we compare the theoretical and simulation results under different SNRs. We set \(N_b=8\), the loop bandwidths \({B_L} = 10\,{\rm{ Hz}}\) and \({B_L} = 100\,{\rm{ Hz}}\), and the other parameters remain unchanged. The steady-state error variance results are presented in Fig. 4.
 |
Fig. 4 Theoretical and simulated steady-state time-delay error variances. |
The x-axis of the Fig. 4 is SNR of the received signal, which is \(\mathrm{SNR}=10 \lg \left(P T_s / N_0\right)\). The y-axis is the normalized steady-state time-delay error variance. The theoretical results are calculated by (55), and the simulation result corresponding to each SNR under each loop bandwidth \(B_L\) is obtained by performing \(50\) Monte Carlo simulations, and the duration of the steady-state tracking used in each simulation is set to \(1\) second to obtain a stable error variance result.
As shown in Fig. 4, the theoretical and simulated steady-state time-delay error variances bear again a good agreement. The results demonstrate that the variance of the steady-state delay error decreases as the loop bandwidth decreases and the SNR of the received signal increases. Additionally, the simulation results are slightly greater than the theoretical values, which may be attributed to the approximations used in the intermediate steps during the theoretical derivation process.
Finally, we compare the steady-state time-delay error variance performance of the proposed code tracking method, the correlation DLL method and the selection DLL method [6], [17]. The parameters settings remain unchanged and the results are shown in Fig. 5. The error variance corresponding to each SNR under each method is obtained by performing \(50\) Monte Carlo simulations with the duration of the steady-state tracking used in each simulation is set to \(1\) second. The y-axis has been normalized using the sampling period of the proposed method.
 |
Fig. 5 Steady-state time-delay error variances of the proposed method, the correlation DLL method and the selection DLL method. |
As shown in Fig. 5, the correlation DLL method performs slightly worse than the proposed code tracking method under different SNRs and loop bandwidths. The reason may be that the correlation DLL method uses cubic interpolation for delay adjustment, which has a larger error than the frequency domain delay adjustment method. Besides that, under low SNRs, the steady-state time-delay error variance performance of the selection DLL method is worse than the other two methods. This is due to the fact that the probability of selecting the wrong correlation outputs is higher for the selective DLL method under low SNRs, leading to deteriorated error variance performance. Therefore, the steady-state time-delay error variance performance of the proposed method is similar to that of the correlation DLL method, and the selected DLL method has a severe performance deterioration under low SNRs.
In summary, compared with the correlation DLL method, the proposed code tracking method in this paper has advantages in maximum and linear code tracking ranges, which means lower probability of loss of lock and higher stability, with constant error variance performance and appropriate computational complexity. In contrast, the selection DLL method has a severe error variance performance deterioration under low SNRs.
8. Conclusion
This paper proposes a high stability code tracking method for band-limited DSSS systems, which has a larger maximum and linear code tracking ranges, constant error variance performance and appropriate computational complexity. Then we analyze the code tracking range, steady-state error variance and computational complexity of the proposed method, and give a closed-form expression for the steady-state error variance, which can be used to analyze the error variance performance theoretically of the proposed method and design proper band-limited DSSS systems. Beside, simulations are also made to analyze and verify the proposed method. The results of both the theoretical analysis and simulations show that the proposed code tracking method has a larger maximum and linear code tracking ranges which means lower probability of loss of lock, constant error variance performance and appropriate computational complexity, and it can be used in band-limited DSSS systems and improve its stability.
References
[1] A. Modenini and B. Ripani, “A tutorial on the tracking, telemetry, and command (TT&C) for space missions,” IEEE Commun. Surveys Tutss, vol.25, no.3, pp.1510-1542, 2023.
CrossRef
[2] K. Tomita, M. Okumura, and E. Okamoto, “Demonstration of chaos-based radio encryption modulation scheme through wired transmission experiments,” IEICE Trans. Commun., vol.E106-B, no.8, pp.686-695, Aug. 2023.
CrossRef
[3] K.J. Quirk and M. Srinivasan, “PN code tracking using noncommensurate sampling,” IEEE Trans. Commun., vol.54, no.10, pp.1845-1856, 2006.
CrossRef
[4] W. Liu, Y. Hu, T.H. Hsieh, J. Zhao, and S. Wang, “Quinary offset carrier modulations for global navigation satellite system,” IEICE Trans. Commun., vol.E104-B, no.5, pp.563-569, May 2021.
CrossRef
[5] J.W. Betz and K.R. Kolodziejski, “Generalized theory of code tracking with an early-late discriminator part I: Lower bound and coherent processing,” IEEE Trans. Aerosp. Electron. Syst., vol.45, no.4, pp.1538-1556, 2009.
CrossRef
[6] R. De Gaudenzi, M. Luise, and R. Viola, “A digital chip timing recovery loop for band-limited direct-sequence spread-spectrum signals,” IEEE Trans. Commun., vol.41, no.11, pp.1760-1769, 1993.
CrossRef
[7] F.J. Harris and M. Rice, “Multirate digital filters for symbol timing synchronization in software defined radios,” IEEE J. Sel Areas Commun., vol.19, no.12, pp.2346-2357, 2001.
CrossRef
[8] Z. Gao, M. Zhou, P. Reviriego, and J.A. Maestro, “Efficient fault-tolerant design for parallel matched filters,” IEEE Trans. Circuits Syst. II, Exp. Briefs, vol.65, no.3, pp.366-370, 2017.
CrossRef
[9] M. Srinivasan, C.C. Chen, G. Grebowsky, and A. Gray, “An all-digital, high data-rate parallel receiver,” JPL TDA Progress Report, vol.42, p.131, 1997.
[10] C. Wang, C. Lin, Q. Chen, B. Lu, X. Deng, and J. Zhang, “A 10-gbit/s wireless communication link using 16-QAM modulation in 140-GHz band,” IEEE Trans. Microw. Theory Techn., vol.61, no.7, pp.2737-2746, 2013.
CrossRef
[11] Z. Lu and Y. Jiao, “Efficiently all-digital code tracking for band-limited DSSS systems,” IEEE Commun. Lett., vol.27, no.2, pp.686-690, 2023.
CrossRef
[12] J.H. Won, “A novel adaptive digital phase-lock-loop for modern digital gnss receivers,” IEEE Commun. Lett., vol.18, no.1, pp.46-49, 2013.
CrossRef
[13] I. Cortés, J.R. Van Der Merwe, A. Rügamer, and W. Felber, “Adaptive loop-bandwidth control algorithm for scalar tracking loops,” 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), pp.1178-1188, IEEE, 2020.
CrossRef
[14] X. Chen, F. Dovis, S. Peng, and Y. Morton, “Comparative studies of GPS multipath mitigation methods performance,” IEEE Trans. Aerosp. Electron. Syst., vol.49, no.3, pp.1555-1568, 2013.
CrossRef
[15] H. Wang, Q. Chang, Y. Xu, and X. Li, “Adaptive narrow-band interference suppression and performance evaluation based on code-aided in GNSS inter-satellite links,” IEEE Syst. J., vol.14, no.1, pp.538-547, 2019.
CrossRef
[16] A. Wilde, “Extended tracking range delay-locked loop,” Proc. IEEE International Conference on Communications ICC’95, pp.1051-1054, IEEE, 1995.
CrossRef
[17] A. Wilde, “Correlation branch selection for extended tracking range delay-locked loops,” European Transactions on Telecommunications, vol.9, no.1, pp.57-64, 1998.
CrossRef
[18] S.C. Pei and Y.C. Lai, “Closed form variable fractional time delay using FFT,” IEEE Signal Process. Lett., vol.19, no.5, pp.299-302, 2012.
CrossRef
[19] M. Blok, “Comments on “closed form variable fractional time delay using FFT”,” IEEE Signal Process. Lett., vol.20, no.8, pp.747-750, 2013.
CrossRef
[20] F.M. Gardner, Phaselock Techniques, John Wiley & Sons, 2005.
CrossRef
[21] B.D. Rao and K.S. Hari, “Performance analysis of root-music,” IEEE Trans. Acoust., Speech, Signal Process., vol.37, no.12, pp.1939-1949, 1989.
CrossRef
[22] P.J. Brockwell and R.A. Davis, Time Series: Theory and Methods, Springer Science & Business Media, 2009.
CrossRef
Appendix A: Derivation of the Frequency Domain Waveform \({G_k}\)
We set \(g_{\ell} =\operatorname{IFFT}_{N_s}\left[G_k\right]\), and get the expression of \({G_k}\) by solving \(g_{\ell}\).
Finding the solution for \(g_{\ell}\) is generally a simplifying process of the correlation operation. We use the corresponding relationship between frequency domain conjugate multiplication and time domain circular convolution and \([\cdot]_{N_s}\) to represent the \(N_s\) points cycle extension of a signal. Then \(g_{\ell}\) can be represented as (A\(\cdot\)1).
| \[\begin{align} g_{\ell} & =\operatorname{IFFT}_{N_s}\left[G_{k}\right] \notag\\ & =\operatorname{IFFT}_{N_s}\left[{\left( {G_k^{RX}} \right)^*}G_k^{TX}\right] \notag\\ & =\frac{T_s}{N} \sum_{i=0}^{N_s-1}\left[\sum_{m_1=0}^{N-1} c_{m_1} p_R\left(i T_s-m_1 T_c\right)\right]_{N_s}\left[\sum_{m_2=0}^{N-1} c_{m_2} p_T\left[(i+\ell) T_s-m_2 T_c\right]\right]_{N_s} \notag\\ & =\frac{T_s}{N} \sum_{i=0}^{N_s-1} \sum_{m_1=0}^{N-1} c_{m_1}\left[p_R\left(i T_s-m_1 T_c\right)\right]_{N_s} \sum_{m_2=0}^{N-1} c_{m_2}\left[p_T\left[(i+\ell) T_s-m_2 T_c\right]\right]_{N_s} \notag\\ & =\frac{T_s}{N} \sum_{i=0}^{N_s-1} \sum_{m_1=0}^{N-1} c_{m_1} \sum_{x_1=-\infty}^{\infty} p_R\left[i T_s-\left(m_1+x_1 N\right) T_c\right] \sum_{m_2=0}^{N-1} c_{m_2} \sum_{x_2=-\infty}^{\infty} p_T\left[(i+\ell) T_s-\left(m_2+x_2 N\right) T_c\right] \notag\\ & =\frac{T_s}{N} \sum_{m_1=0}^{N-1} \sum_{m_2=0}^{N-1} c_{m_1} c_{m_2} \sum_{i=0}^{N_s-1} \sum_{x_1=-\infty}^{\infty} p_R\left[i T_s-\left(m_1+x_1 N\right) T_c\right] \sum_{x_2=-\infty}^{\infty} p_T\left[(i+\ell) T_s-\left(m_2+x_2 N\right)T_c\right]. \tag{A$\cdot $1} \end{align}\] |
Because the pulse waveforms \(p_R(t)\) and \(p_T(t)\) are approximately equal to \(0\) when \(\vert t\vert > T/2\), where \(T=NT_c=N_sT_s\). The finite accumulations about \(i\) and infinite accumulations about \(x_1\), \(x_2\) in (A\(\cdot\)1) can be rewritten as an infinite accumulation about \(i\). Meanwhile, the discrete convolution \(p_R\left(i T_s\right) * p_T\left(i T_s\right)=p\left(i T_s\right) / T_s\) corresponds to the continuous convolution \(p_R(t) * p_T(t)=p(t)\). \(p_T(t)\) is the RRC pulse waveform of the transmitted signal, \(p_R(t)\) is the local RRC pulse waveform, and \(p(t)\) is the RC pulse waveform, with the Fourier transform \(P(f) = {T_c}{P_N}(f)\). Thus, (A\(\cdot\)1) can be further simplified to (A\(\cdot\)2).
| \[\begin{align} g_{\ell}&=\frac{T_s}{N} \sum_{m_1=0}^{N-1} \sum_{m_2=0}^{N-1} c_{m_1} c_{m_2} \sum_{i=-\infty}^{\infty} p_R\left(i T_s-m_1 T_c\right)p_T\left[(i+\ell) T_s-m_2 T_c\right] \notag \\ &=\frac{1}{N} \sum_{m_1=0}^{N-1} \sum_{m_2=0}^{N-1} c_{m_1} c_{m_2} p\left[\ell T_s+\left(m_1-m_2\right) T_c\right] \notag \\ &= p\left(\ell T_s\right)+\frac{1}{N} \sum_{m_1=0}^{N-1} \sum_{\substack{m_2=0 \\ m_2 \neq m_1}}^{N-1} c_{m_1} c_{m_2} p\left[\ell T_s+\left(m_1-m_2\right) T_c\right]. \tag{A$\cdot $2} \end{align}\] |
Based on the orthogonality property of spread spectrum code, the latter term in (A\(\cdot\)2) is approximately equal to \(0\). Therefore, (A\(\cdot\)2) can be further simplified as
| \[\begin{align} g_{\ell} \approx p\left(\ell T_s\right). \tag{A$\cdot $3} \end{align}\] |
At the same time, the Fourier transform of \(p(t)\) is known to be \(P(f) = {T_c}{P_N}(f)\), then \(G_k\) can be expressed finally as
| \[\begin{align} G_k&=\operatorname{FFT}_{N_s}\left[ g_{\ell}\right] \notag \\ &=\operatorname{FFT}_{N_s}\left[ p\left(\ell T_s\right)\right] \notag \\ &=\frac{T_c}{T_s} P_N\left(\frac{k}{N_s T_s}\right). \tag{A$\cdot $4} \end{align}\] |
Appendix B: Derivation of the Frequency Domain Noise Covariance Matrices \(\overline{{\mathbf{\tilde{W}}}_n{\mathbf{\tilde{W}}}_n^H}\) and \(\overline{{{{\bf{V}}}_n}{\bf{ V}}_n^H}\)
We define \(\tilde w_{\ell,n}=\mathrm{IFFT}_{N_s}\left[ \tilde{W}_{k,n}\right]\), then \(\tilde w_{\ell,n}\) can be expressed as (A\(\cdot\)5), where \(w'(t) = {p_R}(t) * w(t)\).
| \[\begin{align} \tilde{w}_{\ell,n} &=\mathrm{IFFT}_{N_s}\left[ \tilde{W}_{k,n}\right] \notag \\ &=\mathrm{IFFT}_{N_s}\left[ \left(G_k^{RX}\right)^* W_{k,n}\right] \notag \\ & =\frac{{{T_s}}}{N} \sum_{i=0}^{N_s-1}\left[\sum_{m_1=0}^{N-1} c_{m_1} p_R\left(i T_s-m_1 T_c\right)\right]_{N_s} w\left[(i+\ell+nN_s) T_s\right] \notag\\ & = \frac{{{T_s}}}{N}\sum_{i=0}^{N_s-1} \sum_{m_1=0}^{N-1} c_{m_1} \sum_{x_1=-\infty}^{\infty} p_R\left[i T_s-\left(m_1+x_1 N\right) T_c\right] w\left[(i+\ell+nN_s) T_s\right] \notag \\ & =\frac{{{T_s}}}{N} \sum_{m_1=0}^{N-1} c_{m_1} \sum_{i=-\infty}^{\infty} p_R\left[i T_s-\left(m_1+x_1 N\right) T_c\right] w\left[(i+\ell+nN_s) T_s\right] \notag\\ & =\frac{1}{N} \sum_{m_1=0}^{N-1} c_{m_1} w^{\prime}\left[(\ell+nN_s) T_s\right], \tag{A$\cdot $5} \end{align}\] |
Furthermore, given that the PSD of \(w(t)\) is \({N_0}/P\), and the frequency response of \({p_R}(t)\) is \({P_R}(f) = \sqrt {{P_N}(f)}\), thus the PSD of \(w'(t)\) is \(N_0{P_N}(f)/P\).
As the spread spectrum code is uncorrelated with the noise term in (A\(\cdot\)5), the PSD of \(\tilde w(t)\) can be represented by \(N_0{P_N}(f)/\left(NP\right)\). Assuming the noise bandwidth of the baseband signal is \(1/T_s\), we have \({\tilde w_{\ell,n} } = \tilde w\left[(\ell+nN_s) {T_s}\right]\). Thus, the digital PSD of \({\tilde w_{\ell,n} }\) is
| \[\begin{align} {S_{\tilde w}}({e^{j\omega }}) = \frac{{{N_0}}}{{PNT_s}}{P_N}(\frac{\omega}{2\pi T_s}). \tag{A$\cdot $6} \end{align}\] |
Then we can get that
| \[\begin{align} \overline{\tilde{W}_{k,n} \tilde{W}_{k,n}^\ast}=\frac{{{N_s N_0}}}{{PNT_s}}{P_N}(\frac{k}{N_sT_s}). \tag{A$\cdot $7} \end{align}\] |
Since the noise at different frequency indices is uncorrelated, the frequency domain noise covariance can be expressed as
| \[\begin{align} \overline {{{{\bf{\tilde W}}}_n}{\bf{\tilde W}}_n^H} = \frac{{{N_s}{N_0}}}{{PN{T_s}}} \left[ \begin{array}{cccc} {P_N}(0)& & & \\ & {{P_N}(\frac{1}{{{N_s}{T_s}}})} & & \\ & & \ddots & \\ & & {{P_N}(\frac{{{N_s} - 1}}{{{N_s}{T_s}}})}& \end{array} \right]. \tag{A$\cdot $8} \end{align}\] |
The digital PSD of \(v_{\ell,n}\), where \(v_{\ell,n}=\mathrm{IFFT}_{N_b}\left[ {V}_{k,n}\right]\), is approximately equal to that of \({\tilde w_{\ell,n}}\), thus we can get that
| \[\begin{align} \overline{{V}_{k,n}{V}_{k,n}^*}=\frac{{{N_b N_0}}}{{PNT_s}}{P_N}(\frac{k}{N_bT_s}). \tag{A$\cdot $9} \end{align}\] |
Then the frequency domain noise covariance after integration and dump can be expressed as
| \[\begin{align} \overline {{{{\bf{V}}}_n}{\bf{ V}}_n^H} = \frac{{{N_b}{N_0}}}{{PN{T_s}}} \left[ \begin{array}{cccc} {P_N}(0)& & & \\ & {{P_N}(\frac{1}{{{N_b}{T_s}}})} & & \\ & & \ddots & \\ & & {{P_N}(\frac{{{N_b} - 1}}{{{N_b}{T_s}}})}& \end{array} \right]. \tag{A$\cdot $10} \end{align}\] |
Authors
Space Engineering University
received the B.S. degree from Space Engineering University, Beijing, China, in 2019, where he is currently pursuing the Ph.D. degree in information and communication engineering. His research interests include space TT&C technology, spread spectrum technology, and signal processing.
 |
Space Engineering University
received the B.S. degree in electronic and information technology from the Beijing University of Aeronautics and Astronautics, in 2008, and the M.S. and Ph.D. degrees from the Institution of Command and Technology of Equipment, Beijing, in 2010 and 2015, respectively. He is currently working as an Associate Professor at Space Engineering University, Beijing, China. His research interests include space TT&C technology, high-speed signal processing, and advanced signal processing.
 |
Space Engineering University
received the B.S. degree from Space Engineering University, Beijing, China, in 2018, where he is currently pursuing the Ph.D. degree. His research interests include TT&C of spacecraft, target tracking and signal processing.
 |
