# A Sparse Temporal Synchronization Algorithm of Laser Communications for Feeder Links in 5G Nonterrestrial Networks

Authors: Lichen Zhu; Hangcheng Han; Xiangyuan Bu; Jichao Wang
Journal: Wireless Communications and Mobile Computing
 (2018)

---

## Abstract

To foster the rollout of 5G in unserved areas, 3GPP has kicked off a study item on new radio to support nonterrestrial networks (NTNs). Due to ultra-wideband of laser, laser communication is very promising for the feeder links of NTNs; however, imprecise temporal synchronization hinders its deployment, which results from a combination of propagation delay, velocity, acceleration, and jerk of NTN platform. The prior synchronization algorithms are inapplicable to the temporal synchronization in laser communications due to the extremely high data rate and Doppler shift. This paper is devoted to addressing the temporal synchronization problem in laser communications. In particular, we first observe the sparsity of laser signal in time-frequency domain. On top of this observation, we propose a new sparsity-aware algorithm for temporal synchronization without carrier aid through sparse discrete polynomial-phase transform and sparse discrete fractional Fourier transform. Subsequently, we implement the proposed algorithm via designing a hardware prototype. To further evaluate its performance, we conduct extensive simulations, and the results demonstrate the effectiveness of the proposed algorithm in terms of good accuracy, low power consumption, and low computational complexity, suggesting its attractiveness for the feeder links of 5G NTNs.


---

## Body

## 1. Introduction

With the exponential growth of wireless traffic volume, new radio (NR) becomes the foundation of 5G to provide universal coverage [1, 2]. To foster the rollout of 5G in unserved areas, nonterrestrial networks (NTNs) have been kicked off by 3GPP [3–5]. In this context, the major challenge is that the feeder links between the gateway and the NTN platforms must be broadband; i.e., the data rate should be on the order of gigabits per second (Gbps) [6–8]. Fortunately, laser is innately of the ultra-wideband, and the laser communication is thus very promising for the feeder links in satisfying the data rate requirement, while it is difficult for laser communications to achieve precise temporal synchronization because of the high data rate, e.g., the order of Gbps, and the high Doppler resulting from the movement of NTN platforms [9, 10].To address the temporal synchronization problem, a large body of works have been proposed, which can be classified into three categories: tracking loop-based methods, transformation-domain methods, and sparse transformation-domain methods.Specifically, tracking loop is a type of conventional temporal synchronization algorithms [11–14]. In [11], the authors developed a novel architecture for tracking loop to accurately measure the aircraft’s speed, which is helpful for temporal synchronization. In [12, 13], the performance of tracking loop on ultra-wideband impulse radio was analyzed. A common tracking loop called digital delay locked loop (DDLL) was discussed in [14] for tracking direct sequence spread spectrum signal with high Doppler. As analyzed in these works, tracking loop algorithms are suitable for systems with low Doppler and low data rate. While the convergence rate degrades rapidly with the increase of Doppler and Doppler changing rate. Furthermore, high data rate system requires parallel architecture of tracking loop, which needs more resources than the serial one because multiple samples have to be handled during one clock cycle [15, 16] and cannot meet the low power consumption requirement on NTN platforms [17].The transformation-domain algorithms [18–21] follow open loop principle and are able to maintain a constant processing delay in high Doppler environments. In [18], the discrete fractional Fourier transform (DFrFT) was adopted for acceleration and velocity estimation. In [19], the discrete polynomial-phase transform (DPT) was used to estimate the aircraft dynamic parameters. In [20], the authors studied the time-frequency characteristics from the perspective of signal time-frequency distribution, i.e., Wigner-Ville distribution, to facilitate the temporal synchronization. In [21], the keystone transform was employed to solve the target range migration problem in temporal synchronization. Although transformation-domain algorithms are suitable for the scenario with high Doppler, their computational complexity is usually very high, which violates the low power consumption requirement as well.To reduce the complexity, the sparse transformation-domain algorithm is introduced into high Doppler and high data rate systems [22–25]. The authors in [22] proposed a sparse algorithm based on the fast Fourier transform (FFT), namely, sparse FFT (SFFT). The applications of SFFT were introduced in [23], and the authors in [24] further evaluated the implementation performance of SFFT. The paper [25] designed a sparse using the discrete fractional Fourier transform (DFrFT) for the temporal synchronization in the scenario with high Doppler changing rate. With low complexity, the sparse algorithm is recommended to low power consumption platforms. However, these sparse transformation-domain algorithms require carrier recovery. It is pointed out that carrier recovery is difficult to implement on laser communications for NTN platforms because of the complex structure [26] and then restricts the application of the sparse transformation-domain algorithm on the feeder links of 5G NTNs.Motivated from the observations above, we employ an incoherent laser transmission system called intensity modulation direct detection (IMDD) for its simplicity and low cost [27, 28]. The IMDD system can be modeled as a Gaussian channel with the positive real input signal and input-dependent noise [26, 29]. This paper designs a sparse transformation-domain algorithm in IMDD-based laser communications for 5G NTNs, which is of good accuracy, low power consumption, and low computational complexity. The primary contributions of the paper are summarized as follows.(1)
We develop a temporal synchronization algorithm based on the sparse pilot. Armed with the sparse DPT (SDPT) and the sparse DFrFT (SDFrFT), the proposed algorithm is able to estimate the propagation delay, velocity, acceleration, and jerk without carrier aid in high dynamic environment.(2)
We analyze the accuracy and complexity of the proposed algorithm and compare its performance with the nonsparse algorithms and the DDLL algorithms. The analytical and experimental results show that the proposed algorithm performs similar to the nonsparse algorithm and superior to the DDLL algorithm in terms of accuracy. While the proposed algorithm can achieve lower complexity, making it the most suitable for scenarios with high Doppler and low power consumption among these three algorithms.(3)
We discuss the implementation issues including module reuse analysis and clock rate analysis. According to the analysis, the proposed algorithm can be implemented with less clock rate than the conventional DDLL algorithm.(4)
We implement the proposed algorithm via designing a hardware prototype, and the results agree well with the theoretical analysis, demonstrating that the algorithm can be implemented on resource constrained platforms, e.g., 5G NTN platforms.The rest of the paper is organized as follows. In Section2, the system model is described. The proposed temporal synchronization algorithm is then presented in Section 3 and analyzed in Section 4. The related implementation issues are shown in Section 5. Simulation and experiment results are given in Sections 6 and 7 to demonstrate the effectiveness of the proposed algorithm, respectively, which are followed by the conclusions drawn in Section 8.


## 2. System Model

The system model of NTNs is presented in Figure1, where the IMDD-based laser communications serve for the feeder link between spaceborne platform and gateway.Figure 1
The system model of NTNs.The laser signal is received by photodetection and is transformed to an electrical signal with positive and negative level, which is expressed as(1)st=A0gt-τt+wt=A0gt-a0+a1t+a2t2+a3t3+wt,where gt is the waveform of signal with amplitude A0, τt denotes the propagation delay, which contains four components related to distance (a0), velocity (a1), acceleration (a2), and jerk (a3), respectively, and wt is the zero-mean additive white Gaussian noise with variance σ2. The experimental results show that a1≫a2≫a3 [11]. The nature of temporal synchronization is to estimate τt, namely, a0, a1, a2, and a3, which can be obtained with the aid of pilot.Figure2 shows the pilot position in the transmitted frame, where Lfrm, Lh, and Lp denote the lengths of the entire frame, the frame header, and the pilot, respectively. The pilot waveform, denoted by gpt, is a sequence of periodic square waves and is expressed as(2)gpt=12,kTsyb<t≤k+12Tsyb,-12,k+12Tsyb<t≤k+1Tsyb,where Tsyb denotes the symbol duration and k=0,1,2,…,Lp-1. According to the Fourier series, we rewrite (2) as(3)gpt=∑k=1∞12k-1sin⁡22k-1πfct,where fc=1/2Tsyb. Through low-pass filtering, (3) is converted to(4)gpt=∑k=1∞Ak2k-1sin⁡22k-1πfct,Ak∈0,1,where Ak denotes the amplitude of the kth order harmonic, which decreases rapidly with the increase of k.Figure 2
The pilot in the transmitted frame.After analog-digital conversion, from (1) and (4), we have(5)sn=rn+wn=A0gnTs-τnTs+wn,gpnTs=∑k=1∞Ak2k-1sin⁡22k-1πfcnTs,where Ts denotes the sample interval.To simplify the computation, we assumeAk=0 when k>3; then gpnTs in (5) is rewritten as(6)gpnTs=∑k=13Ak2k-1sin⁡22k-1πfcnTs.


## 3. The Proposed Temporal Synchronization Algorithm

AsgpnTs is the sum of sinusoidal functions, the parameters estimation can be expressed as follows according to the maximum likelihood principle [30, 31],(7)a^1,a^2,a^3=argmaxa1,a2,a3⁡∑nscne-j2πfc-a1nTs-a2nTs2-a3nTs3,a^0=arg⁡∑nscne-j2πfc-a^1nTs-a^2nTs2-a^3nTs3,where scn=sn+j·sπ/4n, and sπ/4n represents sn with π/4 phase shift. Based on (7), we will elaborate how to perform the temporal synchronization in four steps, as shown in Figure 3. We conduct the preprocess in the first step to facilitate the following three steps. In the second step, as estimating a^3 is difficult within limited process duration in high Doppler environment, the parameter a^2 will be preferentially derived based on a two-order SDPT (SDPT2). In the third step, the parameters a^0 and a^1 will be obtained based on SDFrFT, and a^2 will also be improved accordingly. The parameter a^3 will be presented based on a three-order SDPT (SDPT3) after resampling in the fourth step. Finally, we output all the estimated results, namely, τ^nTs.Figure 3
The diagram of the sparsity-based scheme.


### 3.1. The First Step: Preprocess

To prepare data for estimation, preprocess contains three stages.First, the pilot is located with the aid of frame header. Based on the cross-correlation between the received signal and the template in header, the pilot is located by(8)n^p=arg⁡maxn⁡∑b=0fsLhTsyb-1sn·ghn-bTs+fsLhTsyb,where · denotes the floor function, fs=1/Ts denotes the sample rate, and ghnTs is the header waveform.Second, we derive the complex-form ofsn to meet the requirement of (7). Considering the π/4 phase shift of gpnTs, namely,(9)
g
p
-
π
/
4
n
T
s
=
g
p
n
T
s
-
T
s
y
b
2
=
∑
k
=
1
3
A
k
2
k
-
1
sin
⁡
2
2
k
-
1
π
f
c
n
T
s
-
T
s
y
b
2
=
∑
k
=
1
3
-
1
k
-
1
A
k
2
k
-
1
cos
⁡
2
2
k
-
1
π
f
c
n
T
s
,and the complex-form of gpnTs is expressed as(10)gcnTs=gp-π/4nTs+j·gpnTs=∑k=13-1k-1Ak2k-1exp⁡-1k-1j22k-1πfcnTs=ATP,where A=A1,-A2/3,A3/5T∈R3, and P=P1,P2,P3T∈C3 with Pk=exp⁡j-1k-122k-1πfcnTs-τnTs and k=1,2,3. From (10), the complex-form of sn is (11)scn=A0gcnTs+wn=A0TP+wn,where A0=A0A1,-A0A2/3,A0A3/5T∈R3.Third, since the pilot has been located, we obtainNseq entries by subsampling the pilot data start from n^p with subsample rate rsmp=βLfrmTsyb/Ts, where β is a positive integer. The challenge here lies in subsampling. The subsequent samples will drift away from the original position caused by Doppler, and the sample position may be out of the pilot range. Hence, a sufficient length of pilot, i.e., Lp, is required to prevent the sampling position out of range. In the following, we discuss the minimum value of Lp.LettingLm denote the maximum range of drift, we have(12)Lm=Nseqa1mβLfrmc,where a1m denotes the maximum velocity. The reason for ignoring a2 and a3 is that the platform cannot accelerate any more after it reaches the maximum velocity. Generally, as the drift direction is unknown, the lower limit of Lp is 2Lm, where · denotes the ceiling function. An improved method can reduce the lower limit of Lp to Lm by simultaneously obtaining 3 sample sets that start from n^p, n^p+Lm/2fsTsyb, and n^p+LmfsTsyb, respectively, among which at least a set will not be out of range. Then we select the best one among these sets with the following method. Since the frame header is exploited to locate the frame, it assists in determining the drift direction and magnitude, as shown in Figure 4, where L1st and Llast denote the intervals from the sample position to the header in the first frame and to the last frame, respectively, and Ld denotes the drift value with Ld=L1st-Llast. The first, or second, or third sample set is selected when Ld>Lp/2, Ld∈-Lp/2,Lp/2, or Ld<-Lp/2, respectively.Figure 4
The method of determining drift direction and magnitude.After subsampling,scn is converted to scm with m=0,1,…,Nseq-1, and the sample rate and sample interval turn into that fs′=fs/rsmp and that Ts′=rsmp·Ts, respectively. In the following section, we focus on the signal processing of the proposed temporal synchronization.


### 3.2. The Second Step: Derivea^2 Based on SDPT2

According to (7), it is possible to jointly estimate all the four unknown parameters. Unfortunately, the complexity of joint estimation is unaffordable. To reduce the complexity, we decompose the parameter estimation problem into several subproblems to estimate the parameters one by one, which is essentially the same as the joint estimation because the parameters to be estimated have a weak association with each other and the separate estimation results in no performance loss. Generally, the highest-order parameter, a3, is estimated firstly [32]. However, as mentioned above, high Doppler makes the sample position drift away and leads to a limited process duration. Moreover, as a3 represents the third derivative of distance, it cannot be accurately estimated in such short duration. In comparison, the acceleration estimation is much easier, and we thus preferentially estimate a2.As a reduced-order method to estimatea2, DPT2 is defined as [32](13)DPT2scm,f,ξ≜DFTDP2scm,ξ,f=∑m=0Nseq-1DP2scm,ξe-j2πfmTs′,where DFT denotes the discrete Fourier transform with rotation factor e-j2πfmTs′, ξ is a positive integer, and DP2 is given by(14)DP2scm,ξ≜scmsc∗m-ξ,where ·∗ denotes the conjugate operation. For simplicity, let xm=DP2scm,ξ, and we have (15)xm=scmsc∗m-ξ=rm+wmr∗m-ξ+w∗m-ξ=rmr∗m-ξ+rmw∗m-ξ+r∗m-ξwm+wmw∗m-ξ,where(16)rmr∗m-ξ=A0TPPξTA0=A0TPPξ∗TA0=A0TP1P1∗ξP1P2∗ξP1P3∗ξP2P1∗ξP2P2∗ξP2P3∗ξP3P1∗ξP3P2∗ξP3P3∗ξA0is the principal component of xm; Pξ and Pξ denote P and P with ξ-delay, respectively. The other components of xm are related to wm and can be regarded as the noise terms. Moreover, rmr∗m-ξ is dominated by P1P1∗ξ, P2P2∗ξ, and P3P3∗ξ. The amplitude of P1P1∗ξ is much larger than that of P2P2∗ξ and P3P3∗ξ. Consequently, the components of rmr∗m-ξ except for P1P1∗ξ are regarded as the disturbance terms, and (15) can be rewritten as(17)xm=A0A12P1P1∗ξ+lmTs′+ƛmTs′=A0A12ej2π-2a2ξTs′2+3a3ξ2Ts′3m-3a3ξTs′3m2-ϕ+lmTs′+ƛmTs′,where lmTs′ and ƛmTs′ denote the disturbance terms and noise terms, respectively, and ϕ is a constant. Letting X^k=DFTxm and k^0=argmax⁡X^k, since 3a3ξTs′3m2 and 3a3ξ2Ts′3m are small enough to be ignored, a^2 can be expressed as(18)a^2=-12ξTs′2k^0.Considering the sparsity of DFT result, SFFT can be employed to reduce the complexity, as shown in Algorithm 1.Algorithm 1:Estimate a^2 based on the SFFT algorithm.
Input: x m. Output: X ^ k. 1: Zero padding and the length of xm is changed from Nseq to M. 2: Let L=log2⁡M. 3: For  i=0; i<L;i++  do 4: Tear apart the spectrum randomly by x1m=xσi⋅mmodM and X1k=Xσi-1⋅mmodM, where Xk and X1k denote the expressions of xm and x1m in frequency-domain, respectively, and m,σi,k∈1,M. 5: Apply a flat window function to expand the spectrum range. (1) Let Gk∈1-δ,1+δ,k∈-ε′M,ε′M0,δ,k∉-εM,εM be the window function in frequency-domain, where ε′∈0,1 and ε∈0,1 denote the cutoff frequencies of the passband and stopband, respectively, and δ denotes the ripple amplitude, whose reference value is 1/Mc, where c is a positive integer. (2)  Compute x2m=gm⋅x1m, where gm denotes the expression of Gk in time-domain. The length of gm is ω=oBlog2⁡M/σi. 6: Subsample in frequency-domain. Letting B be the data length after subsampling, we have X3k=FFT∑j=0M/B-1x2m+j⋅B with k∈1,B. 7: Map with a hash function. (1) Define a hash function hσk=σi⋅kB/M. (2) Define an offset function oσik=σi⋅k-hσik⋅B/M. (3) Let Γi denote the support set of the largest l coefficients in X3k. The preimage set of Γi is Ιi, whose size is lM/B, where Ιi=k∈1,M∣hσik∈Γi and B=Ml/log2⁡M/δ. (4) Obtain the l largest spectrum coefficients as X ^ 4 , i k = X 3 h σ k e - j π o σ i k ω / M G o σ i k , k ∈ I i , 0 , k ∈ [ 1 , M ] ∩ I - i . (5)  Record the nonzero position of X^4,ik. 8: End for 9: Letting vk be the occurrence times of coordinate k in the sets and only retaining the coordinate whose occurrence times are larger than L/2, we have I′=k∈I1∪⋯∪IL∣vk>L/2, and the rest terms are zeroes. 10: For each coordinate in I′, we obtain the corresponding spectrum coefficients X^rk with k∈I′ and r=1,2,⋯L. 11: Return  X^k=medianX^rk,k∈I′0,k∉I′, where median⋅ denotes computing the median of a sequence.As the accuracy of SDPT2 is limited by the resolution of DFT, which is given by Δf=fs′/M, accurately estimating a2 cannot be achieved with small M. Thus, we adopt Rife and Newton methods to improve the estimation accuracy by updating k^0 and recalculating a^2 due to (18), as shown in Algorithm 2, whose effectiveness is demonstrated in Figure 5.Algorithm 2:Modify the DFT result based on Rife and Newton methods.
1: Modify the DFT result by Rife method. Since X^k denotes the output of SFFT, the modified result is k ^ 1 = 1 T s ′ k ^ 0 + r X ^ k ^ 0 + r X ^ k ^ 0 + X ^ k ^ 0 + r, where r = 1 , X ^ k ^ 0 + 1 ≥ X ^ k ^ 0 - 1 - 1 , X ^ k ^ 0 + 1 < X ^ k ^ 0 - 1 . 2: Further modify the result by the iteration of Newton method. For a positive integer i, the i+1-th iteration result is expressed as k ^ i + 1 = k ^ i - λ i X ^ ′ k ^ i / X ^ ′ ′ k ^ i, where X^′k^ and X^′′k^ denote the first-order and the second-order derivatives of X^k^, respectively, and λ denotes the tuning step. 3: Replace  k^0  by the latest iteration result,k^N.Figure 5
The effectiveness of Rife and Newton methods.


### 3.3. The Third Step: Obtaina^0, a^1, and the Improved a^2

In this part,a^0, a^1, and the improved a^2 are obtained by Pei sampling-type DFrFT [22], which is given in the following equation:(19)Xαk≜sin⁡α-jcos⁡α⁡Mejcot⁡α/2kUs2·∑m=0M-1scme-j2πkm/Mejcot⁡α/2mTs′2,α≥0jcos⁡α-sin⁡αMejcot⁡α/2kUs2·∑m=0M-1scme-j2πkm/Mejcot⁡α/2mTs′2,α<0≈1Mejcot⁡α/2kUs2·DFTscmejcot⁡α/2mTs′2.with m∈0,M-1 and k∈0,M-1.In (19), α and Us denote the rotation angle and the sample interval of the output, Xαk, respectively. Besides, Ts′×Us=2πsin⁡α/M shall be satisfied. Then we derive the optimal values of k and α by(20)k^m,α^m=arg⁡maxα,k⁡Xαk.From (19) and (20), we obtain the estimated results of a^0, a^1, and a^2 in the way (21)a^2=cot⁡α^m4π,a^1=k^mfs′M,a^0=argXα^mk^m-sgn⁡α^mej1/2cot⁡α^mk^mUs2+α^m+π/2.To reduce the computational complexity, we replace DFT in the implementation of Pei sampling-type DFrFT with SFFT due to the sparse spectrum. The implementation of DFrFT based on SFFT is called SDFrFT, which exhibits a similar performance to DFrFT, as shown in Figure6.Figure 6
The comparison between DFrFT and SDFrFT.In addition, Rife and Newton methods can be employed again to replacek^m with k^N, whose procedure is similar to that in Algorithm 2 and is omitted for brevity. Then we substitute k^N into (19), (20), and (21) to recalculate a^1 and a^0.


### 3.4. The Fourth Step: Resample and Derivea^3 Based on SDPT3

As mentioned above, high dynamic, especially high velocity brings the drift of sample position and results in an insufficient process duration, which is contrary to accurately estimatinga3. A method to mitigate the drift by updating the sample rate is introduced as follows.Note that the velocity,a^1, has been derived from (21); it can be substituted into the following formula to update the sample rate:(22)fs′′=2a^1+Rbfs′Rb.The sample interval turns into that Ts′′=1/fs′′ accordingly. As fs′′ matches with the change of data rate caused by Doppler, the drift of sample position is mitigated, and estimating a3 becomes easier with the increase of process duration.In the following, SDPT3 is introduced to estimate a3. Assuming Nseq′ entries are obtained after resampling, we have(23)SDPT3scm,f,ξ≜SFFTDP3scm,ξ,f=∑m=0Nseq′-1DP3scm,ξe-j2πfmTs′′,where DP3 is defined as(24)DP3scm,ξ≜DP2DP2scm,ξ,ξ=DP2scmsc∗m-ξ,ξ=scmsc∗m-ξ2scm-2ξ.Letting P2ξ and P2ξ denote P and P with 2ξ-delay, respectively, x2m=DP3scm,ξ, and ƛ′mTs′′ denote the noise terms, (24) can be further expanded as(25)x2m=scmsc∗m-ξ2scm-2ξ=A0TPPξ∗TA0·A0TPξ∗P2ξTA0+ƛ′mTs′′=A0T·PPξ∗T·A0A0T·Pξ∗P2ξT·A0+ƛ′mTs′′=A0T⊗A0T·PPξ∗T⊗Pξ∗P2ξT·A0⊗A0+ƛ′mTs′′,where ⊗ denotes the Kronecker product. Obviously, x2m is dominated by P1P1∗ξP1∗ξP12ξ, P2P2∗ξP2∗ξP22ξ, and P3P3∗ξP3∗ξP32ξ. Furthermore, as the amplitude of P1P1∗ξP1∗ξP12ξ is much larger than that of the other two components, (25) can be rewritten as(26)x2m=P1P1∗ξP1∗ξP12ξ+l′mTs′′+ƛ′mTs′′=A0A14e-j2π6a3ξ2Ts′′3m+2a2ξ2Ts′′2-6a3ξ3Ts′′3+l′mTs′′+ƛ′mTs′′,where l′mTs′′ denotes the disturbance terms. From (26), a^3 is obtained by(27)a^3=-16ξ2Ts′′3argmax⁡SFFTx2m.The estimation accuracy can be further improved by Rife and Newton methods, whose procedure is similar to that in Algorithm 2 and is omitted for brevity.In addition, as longer duration is obtained after resampling, we conduct SDFrFT again to further improve the accuracy ofa^2,a^1, and a^0, whose procedure is similar to that in the third step and is also omitted for brevity.


## 3.1. The First Step: Preprocess

To prepare data for estimation, preprocess contains three stages.First, the pilot is located with the aid of frame header. Based on the cross-correlation between the received signal and the template in header, the pilot is located by(8)n^p=arg⁡maxn⁡∑b=0fsLhTsyb-1sn·ghn-bTs+fsLhTsyb,where · denotes the floor function, fs=1/Ts denotes the sample rate, and ghnTs is the header waveform.Second, we derive the complex-form ofsn to meet the requirement of (7). Considering the π/4 phase shift of gpnTs, namely,(9)
g
p
-
π
/
4
n
T
s
=
g
p
n
T
s
-
T
s
y
b
2
=
∑
k
=
1
3
A
k
2
k
-
1
sin
⁡
2
2
k
-
1
π
f
c
n
T
s
-
T
s
y
b
2
=
∑
k
=
1
3
-
1
k
-
1
A
k
2
k
-
1
cos
⁡
2
2
k
-
1
π
f
c
n
T
s
,and the complex-form of gpnTs is expressed as(10)gcnTs=gp-π/4nTs+j·gpnTs=∑k=13-1k-1Ak2k-1exp⁡-1k-1j22k-1πfcnTs=ATP,where A=A1,-A2/3,A3/5T∈R3, and P=P1,P2,P3T∈C3 with Pk=exp⁡j-1k-122k-1πfcnTs-τnTs and k=1,2,3. From (10), the complex-form of sn is (11)scn=A0gcnTs+wn=A0TP+wn,where A0=A0A1,-A0A2/3,A0A3/5T∈R3.Third, since the pilot has been located, we obtainNseq entries by subsampling the pilot data start from n^p with subsample rate rsmp=βLfrmTsyb/Ts, where β is a positive integer. The challenge here lies in subsampling. The subsequent samples will drift away from the original position caused by Doppler, and the sample position may be out of the pilot range. Hence, a sufficient length of pilot, i.e., Lp, is required to prevent the sampling position out of range. In the following, we discuss the minimum value of Lp.LettingLm denote the maximum range of drift, we have(12)Lm=Nseqa1mβLfrmc,where a1m denotes the maximum velocity. The reason for ignoring a2 and a3 is that the platform cannot accelerate any more after it reaches the maximum velocity. Generally, as the drift direction is unknown, the lower limit of Lp is 2Lm, where · denotes the ceiling function. An improved method can reduce the lower limit of Lp to Lm by simultaneously obtaining 3 sample sets that start from n^p, n^p+Lm/2fsTsyb, and n^p+LmfsTsyb, respectively, among which at least a set will not be out of range. Then we select the best one among these sets with the following method. Since the frame header is exploited to locate the frame, it assists in determining the drift direction and magnitude, as shown in Figure 4, where L1st and Llast denote the intervals from the sample position to the header in the first frame and to the last frame, respectively, and Ld denotes the drift value with Ld=L1st-Llast. The first, or second, or third sample set is selected when Ld>Lp/2, Ld∈-Lp/2,Lp/2, or Ld<-Lp/2, respectively.Figure 4
The method of determining drift direction and magnitude.After subsampling,scn is converted to scm with m=0,1,…,Nseq-1, and the sample rate and sample interval turn into that fs′=fs/rsmp and that Ts′=rsmp·Ts, respectively. In the following section, we focus on the signal processing of the proposed temporal synchronization.


## 3.2. The Second Step: Derivea^2 Based on SDPT2

According to (7), it is possible to jointly estimate all the four unknown parameters. Unfortunately, the complexity of joint estimation is unaffordable. To reduce the complexity, we decompose the parameter estimation problem into several subproblems to estimate the parameters one by one, which is essentially the same as the joint estimation because the parameters to be estimated have a weak association with each other and the separate estimation results in no performance loss. Generally, the highest-order parameter, a3, is estimated firstly [32]. However, as mentioned above, high Doppler makes the sample position drift away and leads to a limited process duration. Moreover, as a3 represents the third derivative of distance, it cannot be accurately estimated in such short duration. In comparison, the acceleration estimation is much easier, and we thus preferentially estimate a2.As a reduced-order method to estimatea2, DPT2 is defined as [32](13)DPT2scm,f,ξ≜DFTDP2scm,ξ,f=∑m=0Nseq-1DP2scm,ξe-j2πfmTs′,where DFT denotes the discrete Fourier transform with rotation factor e-j2πfmTs′, ξ is a positive integer, and DP2 is given by(14)DP2scm,ξ≜scmsc∗m-ξ,where ·∗ denotes the conjugate operation. For simplicity, let xm=DP2scm,ξ, and we have (15)xm=scmsc∗m-ξ=rm+wmr∗m-ξ+w∗m-ξ=rmr∗m-ξ+rmw∗m-ξ+r∗m-ξwm+wmw∗m-ξ,where(16)rmr∗m-ξ=A0TPPξTA0=A0TPPξ∗TA0=A0TP1P1∗ξP1P2∗ξP1P3∗ξP2P1∗ξP2P2∗ξP2P3∗ξP3P1∗ξP3P2∗ξP3P3∗ξA0is the principal component of xm; Pξ and Pξ denote P and P with ξ-delay, respectively. The other components of xm are related to wm and can be regarded as the noise terms. Moreover, rmr∗m-ξ is dominated by P1P1∗ξ, P2P2∗ξ, and P3P3∗ξ. The amplitude of P1P1∗ξ is much larger than that of P2P2∗ξ and P3P3∗ξ. Consequently, the components of rmr∗m-ξ except for P1P1∗ξ are regarded as the disturbance terms, and (15) can be rewritten as(17)xm=A0A12P1P1∗ξ+lmTs′+ƛmTs′=A0A12ej2π-2a2ξTs′2+3a3ξ2Ts′3m-3a3ξTs′3m2-ϕ+lmTs′+ƛmTs′,where lmTs′ and ƛmTs′ denote the disturbance terms and noise terms, respectively, and ϕ is a constant. Letting X^k=DFTxm and k^0=argmax⁡X^k, since 3a3ξTs′3m2 and 3a3ξ2Ts′3m are small enough to be ignored, a^2 can be expressed as(18)a^2=-12ξTs′2k^0.Considering the sparsity of DFT result, SFFT can be employed to reduce the complexity, as shown in Algorithm 1.Algorithm 1:Estimate a^2 based on the SFFT algorithm.
Input: x m. Output: X ^ k. 1: Zero padding and the length of xm is changed from Nseq to M. 2: Let L=log2⁡M. 3: For  i=0; i<L;i++  do 4: Tear apart the spectrum randomly by x1m=xσi⋅mmodM and X1k=Xσi-1⋅mmodM, where Xk and X1k denote the expressions of xm and x1m in frequency-domain, respectively, and m,σi,k∈1,M. 5: Apply a flat window function to expand the spectrum range. (1) Let Gk∈1-δ,1+δ,k∈-ε′M,ε′M0,δ,k∉-εM,εM be the window function in frequency-domain, where ε′∈0,1 and ε∈0,1 denote the cutoff frequencies of the passband and stopband, respectively, and δ denotes the ripple amplitude, whose reference value is 1/Mc, where c is a positive integer. (2)  Compute x2m=gm⋅x1m, where gm denotes the expression of Gk in time-domain. The length of gm is ω=oBlog2⁡M/σi. 6: Subsample in frequency-domain. Letting B be the data length after subsampling, we have X3k=FFT∑j=0M/B-1x2m+j⋅B with k∈1,B. 7: Map with a hash function. (1) Define a hash function hσk=σi⋅kB/M. (2) Define an offset function oσik=σi⋅k-hσik⋅B/M. (3) Let Γi denote the support set of the largest l coefficients in X3k. The preimage set of Γi is Ιi, whose size is lM/B, where Ιi=k∈1,M∣hσik∈Γi and B=Ml/log2⁡M/δ. (4) Obtain the l largest spectrum coefficients as X ^ 4 , i k = X 3 h σ k e - j π o σ i k ω / M G o σ i k , k ∈ I i , 0 , k ∈ [ 1 , M ] ∩ I - i . (5)  Record the nonzero position of X^4,ik. 8: End for 9: Letting vk be the occurrence times of coordinate k in the sets and only retaining the coordinate whose occurrence times are larger than L/2, we have I′=k∈I1∪⋯∪IL∣vk>L/2, and the rest terms are zeroes. 10: For each coordinate in I′, we obtain the corresponding spectrum coefficients X^rk with k∈I′ and r=1,2,⋯L. 11: Return  X^k=medianX^rk,k∈I′0,k∉I′, where median⋅ denotes computing the median of a sequence.As the accuracy of SDPT2 is limited by the resolution of DFT, which is given by Δf=fs′/M, accurately estimating a2 cannot be achieved with small M. Thus, we adopt Rife and Newton methods to improve the estimation accuracy by updating k^0 and recalculating a^2 due to (18), as shown in Algorithm 2, whose effectiveness is demonstrated in Figure 5.Algorithm 2:Modify the DFT result based on Rife and Newton methods.
1: Modify the DFT result by Rife method. Since X^k denotes the output of SFFT, the modified result is k ^ 1 = 1 T s ′ k ^ 0 + r X ^ k ^ 0 + r X ^ k ^ 0 + X ^ k ^ 0 + r, where r = 1 , X ^ k ^ 0 + 1 ≥ X ^ k ^ 0 - 1 - 1 , X ^ k ^ 0 + 1 < X ^ k ^ 0 - 1 . 2: Further modify the result by the iteration of Newton method. For a positive integer i, the i+1-th iteration result is expressed as k ^ i + 1 = k ^ i - λ i X ^ ′ k ^ i / X ^ ′ ′ k ^ i, where X^′k^ and X^′′k^ denote the first-order and the second-order derivatives of X^k^, respectively, and λ denotes the tuning step. 3: Replace  k^0  by the latest iteration result,k^N.Figure 5
The effectiveness of Rife and Newton methods.


## 3.3. The Third Step: Obtaina^0, a^1, and the Improved a^2

In this part,a^0, a^1, and the improved a^2 are obtained by Pei sampling-type DFrFT [22], which is given in the following equation:(19)Xαk≜sin⁡α-jcos⁡α⁡Mejcot⁡α/2kUs2·∑m=0M-1scme-j2πkm/Mejcot⁡α/2mTs′2,α≥0jcos⁡α-sin⁡αMejcot⁡α/2kUs2·∑m=0M-1scme-j2πkm/Mejcot⁡α/2mTs′2,α<0≈1Mejcot⁡α/2kUs2·DFTscmejcot⁡α/2mTs′2.with m∈0,M-1 and k∈0,M-1.In (19), α and Us denote the rotation angle and the sample interval of the output, Xαk, respectively. Besides, Ts′×Us=2πsin⁡α/M shall be satisfied. Then we derive the optimal values of k and α by(20)k^m,α^m=arg⁡maxα,k⁡Xαk.From (19) and (20), we obtain the estimated results of a^0, a^1, and a^2 in the way (21)a^2=cot⁡α^m4π,a^1=k^mfs′M,a^0=argXα^mk^m-sgn⁡α^mej1/2cot⁡α^mk^mUs2+α^m+π/2.To reduce the computational complexity, we replace DFT in the implementation of Pei sampling-type DFrFT with SFFT due to the sparse spectrum. The implementation of DFrFT based on SFFT is called SDFrFT, which exhibits a similar performance to DFrFT, as shown in Figure6.Figure 6
The comparison between DFrFT and SDFrFT.In addition, Rife and Newton methods can be employed again to replacek^m with k^N, whose procedure is similar to that in Algorithm 2 and is omitted for brevity. Then we substitute k^N into (19), (20), and (21) to recalculate a^1 and a^0.


## 3.4. The Fourth Step: Resample and Derivea^3 Based on SDPT3

As mentioned above, high dynamic, especially high velocity brings the drift of sample position and results in an insufficient process duration, which is contrary to accurately estimatinga3. A method to mitigate the drift by updating the sample rate is introduced as follows.Note that the velocity,a^1, has been derived from (21); it can be substituted into the following formula to update the sample rate:(22)fs′′=2a^1+Rbfs′Rb.The sample interval turns into that Ts′′=1/fs′′ accordingly. As fs′′ matches with the change of data rate caused by Doppler, the drift of sample position is mitigated, and estimating a3 becomes easier with the increase of process duration.In the following, SDPT3 is introduced to estimate a3. Assuming Nseq′ entries are obtained after resampling, we have(23)SDPT3scm,f,ξ≜SFFTDP3scm,ξ,f=∑m=0Nseq′-1DP3scm,ξe-j2πfmTs′′,where DP3 is defined as(24)DP3scm,ξ≜DP2DP2scm,ξ,ξ=DP2scmsc∗m-ξ,ξ=scmsc∗m-ξ2scm-2ξ.Letting P2ξ and P2ξ denote P and P with 2ξ-delay, respectively, x2m=DP3scm,ξ, and ƛ′mTs′′ denote the noise terms, (24) can be further expanded as(25)x2m=scmsc∗m-ξ2scm-2ξ=A0TPPξ∗TA0·A0TPξ∗P2ξTA0+ƛ′mTs′′=A0T·PPξ∗T·A0A0T·Pξ∗P2ξT·A0+ƛ′mTs′′=A0T⊗A0T·PPξ∗T⊗Pξ∗P2ξT·A0⊗A0+ƛ′mTs′′,where ⊗ denotes the Kronecker product. Obviously, x2m is dominated by P1P1∗ξP1∗ξP12ξ, P2P2∗ξP2∗ξP22ξ, and P3P3∗ξP3∗ξP32ξ. Furthermore, as the amplitude of P1P1∗ξP1∗ξP12ξ is much larger than that of the other two components, (25) can be rewritten as(26)x2m=P1P1∗ξP1∗ξP12ξ+l′mTs′′+ƛ′mTs′′=A0A14e-j2π6a3ξ2Ts′′3m+2a2ξ2Ts′′2-6a3ξ3Ts′′3+l′mTs′′+ƛ′mTs′′,where l′mTs′′ denotes the disturbance terms. From (26), a^3 is obtained by(27)a^3=-16ξ2Ts′′3argmax⁡SFFTx2m.The estimation accuracy can be further improved by Rife and Newton methods, whose procedure is similar to that in Algorithm 2 and is omitted for brevity.In addition, as longer duration is obtained after resampling, we conduct SDFrFT again to further improve the accuracy ofa^2,a^1, and a^0, whose procedure is similar to that in the third step and is also omitted for brevity.


## 4. Performance Analysis

To further present the advantages of the proposed temporal synchronization algorithm, we will theoretically analyze its performance in terms of estimation accuracy and algorithm complexity in this section.


### 4.1. Estimation Accuracy Analysis

Asa^3, a^2, and a^1 are derived from DPT3, DPT2, and DFrFT, respectively, we discuss the performance of these three algorithms as follows.According to (17) and (26), the result of DPT contains the disturbance terms and the noise terms. The ratios of the principal components to the maximum disturbances in DPT2 and DPT3 are 9A12/A22 and 81A14/A24, respectively. As A1≫A2, the disturbance terms are negligible. The noise terms follow a zero-mean Gaussian distribution with variances σƛ,DPT22 and σƛ,DPT32, as shown in (28) and (29), respectively,(28)σƛ,DPT22=EA0TP+wmA0TP+w∗m-ξ-A0TP22=2A2σ2+σ4(29)σƛ,DPT32=EA0TP+wmA0TP+w∗m-ξ2·A0TP+wm-2ξ-A0TP42=6A6σ2+10A4σ4+4A2σ6+σ8with A=A0A1+A2/3+A3/5.Considering the impact of sample position of pilot signal on the valid signal length, we discuss the SNRs of DPT3, DPT2, and DFrFT when the pilot length is insufficient as follows. Ignoring the impact of a2 and a3, the valid signal length is (30)Nvalid=Lp·cLfrm·a1m.The rest of signal is approximate random and can be regarded as noise, whose length is Nseq-Nvalid. The SNRs of DPT3, DPT2, and DFrFT are(31)SNRDPT3=Nvalid2A0A18Nseq′6A6σ2+10A4σ4+4A2σ6+σ8+Nseq′-NvalidA0A18SNRDPT2=Nvalid2A0A14Nseq2A2σ2+σ4+Nseq-NvalidA0A14SNRDFrFT=Nvalid2A0A12Nseqσ2+Nseq-NvalidA0A12,respectively. The SNRs will be much lower than 0 dB when Nvalid≪Nseq, which leads to an inaccurate estimation result. Fortunately, Nvalid can be obtained by (30) with a given a1m, and Nvalid≈Nseq can be ensured by selecting the data length. Thus, we assume the pilot length is sufficient for estimation, and (31) can be simplified as follows in the high SNR regime:(32)SNRDPT3≈Nseq′A0A186A6σ2+10A4σ4+4A2σ6+σ8≈Nseq′A0A126σ2,SNRDPT2≈NseqA0A142A2σ2+σ4≈NseqA0A122σ2,SNRDFrFT≈NseqA0A12σ2.With the aid of Newton method, the accuracy of a^3, a^2, and a^1 based on DPT3, DPT2, and DFrFT is improved and given by(33)σa^32=62π2SNRDPT3Ts′′2Nseq′2-1≈36σ22π2A0A12Ts′′2Nseq′Nseq′2-1,σa^22=62π2SNRDPT2Ts′′2Nseq′2-1≈12σ22π2A0A12Ts′2NseqNseq2-1,σa^12=62π2SNRDFrFTTs′′2Nseq′2-1≈6σ22π2A0A12Ts′2NseqNseq2-1.With a given SNR, the accuracy of a^3, a^2, and a^1 can be improved by lengthening data or increasing sample interval. Assuming Nseq′=Nseq and Ts′′≈Ts′, the accuracy relationship among a^3, a^2, and a^1 is(34)σa^32≈3σa^22≈6σa^12.Note that the accuracy of a^2 can be further improved by α-search of DFrFT. Letting the estimation error before α-search and the search range of α-search be xa^2 and a^2-3σa^2,a^2+3σa^2, respectively, and assuming that the search times are sufficient, the estimation error after α-search is (35)xa^2′=xa^2+3σa^2,xa^2<-3σa^2,0,xa^2∈-3σa^2,3σa^2,xa^2-3σa^2,xa^2>3σa^2.The expectation and variance of xa^2′ are presented in the following equation:(36)Exa^2′=∫-∞-3σa^2xa^2+3σa^2·pxa^2dxa^2+∫-3σa^23σa^20·pxa^2dxa^2+∫3σa^2+∞xa^2-3σa^2·pxa^2dxa^2=0,varxa^2′=∫-∞-3σa^2xa^2+3σa^22·pxa^2dxa^2+∫-3σa^23σa^20·pxa^2dxa^2+∫3σa^2+∞xa^2-3σa^22·pxa^2dxa^2=2∫-∞-3σa^2xa^22·pxa^2dxa^2+6σa^2∫-∞-3σa^2xa^2·pxa^2dxa^2+9σa^22∫-∞-3σa^2pxa^2dxa^2=101-erf32σa^22≈0.027σa^22,with pxa^2=1/2πσa^22e-xa^22/2σa^22 and erfx=2/π∫0xe-t2dt. If the search range becomes to -2σa^2,2σa^2, the variance turns into 0.35σa^22 accordingly.Then we discuss the accuracy ofa^0, which is derived from DFrFT due to (21). In the high SNR regime, a^0 follows a zero-mean Gaussian distribution with variance(37)σa^02=σ22π2NseqA0A12sinc2⁡Δf,where Δf denotes the offset between the DFT result and the true value.Finally, the performance of SFFT is discussed as follows. LettingX^-Γk=X^k∣k∉Γ, the estimation error probability based on SFFT is [22](38)PX^k-Xk2≥εlX^-Γk22+3δ2X^k12<OlεB.SFFT performs similar to FFT with a sufficient small error probability if l≪B. Consequently, the accuracies of a^3, a^2, a^1, and a^0 based on SFFT are approximately equal to those based on FFT, which demonstrates the agreement between the proposed sparse algorithm and nonsparse algorithm.


### 4.2. Algorithm Complexity Analysis

The algorithm complexity is shown in Table1, where λitrt and λα denote the iteration times of Newton method and the α-search times, respectively. Referring to the values of ω, B, and L in Algorithm 1, the complexities of SFFT and FFT can be expressed as Olog2⁡MMllog2⁡M and OMlog2⁡M, respectively. In the large data size regime, the complexity of FFT is about M times more than that of SFFT.Table 1
The algorithm complexity.
Approach Multiplies Adds Correlation 2 f s L h T s y b 2 f s L h T s y b - 1 DP2 N s e q 0 DFT FFT M 2 log 2 ⁡ M M log 2 ⁡ M SFFT ω + B 2 log 2 ⁡ B + l M B L B log 2 ⁡ B + l M B L DFrFT DFrFT M + M 2 log 2 ⁡ M λ α M log 2 ⁡ M λ α SDFrFT M + ω + B 2 log 2 ⁡ B + l M B L λ α B log 2 ⁡ B + l M B L λ α Rife & Newton methods 3 N s e q + 1 λ i t r t + 2 2 N s e q - 1 λ i t r t + 2In the condition ofNseq=M, considering DP2 and Rife and Newton methods are employed three times, respectively, and SDFrFT and SFFT are employed twice, respectively, the total numbers of multiplications, Cmul, and additions, Cadd, in the proposed sparse algorithm are(39)Cmul≈9λitrt+2λαNseq+2λαlog2⁡NseqNseqllog2⁡Nseq,Cadd≈6λitrtNseq+2λαlog2⁡NseqNseqllog2⁡Nseq.Similarly, the total numbers of multiplications, C˙mul, and additions, C˙add, in the nonsparse algorithm are(40)C˙mul≈9λitrtNseq+λαNseqlog2⁡Nseq,C˙add≈6λitrtNseq+2λαNseqlog2⁡Nseq.Figure7 presents the complexity comparison of the proposed sparse algorithm with the nonsparse algorithm and the DDLL algorithm. The results indicate that the proposed algorithm has a much lower complexity than both the nonsparse algorithm and the DDLL algorithm.Figure 7
The comparison in terms of complexity.


## 4.1. Estimation Accuracy Analysis

Asa^3, a^2, and a^1 are derived from DPT3, DPT2, and DFrFT, respectively, we discuss the performance of these three algorithms as follows.According to (17) and (26), the result of DPT contains the disturbance terms and the noise terms. The ratios of the principal components to the maximum disturbances in DPT2 and DPT3 are 9A12/A22 and 81A14/A24, respectively. As A1≫A2, the disturbance terms are negligible. The noise terms follow a zero-mean Gaussian distribution with variances σƛ,DPT22 and σƛ,DPT32, as shown in (28) and (29), respectively,(28)σƛ,DPT22=EA0TP+wmA0TP+w∗m-ξ-A0TP22=2A2σ2+σ4(29)σƛ,DPT32=EA0TP+wmA0TP+w∗m-ξ2·A0TP+wm-2ξ-A0TP42=6A6σ2+10A4σ4+4A2σ6+σ8with A=A0A1+A2/3+A3/5.Considering the impact of sample position of pilot signal on the valid signal length, we discuss the SNRs of DPT3, DPT2, and DFrFT when the pilot length is insufficient as follows. Ignoring the impact of a2 and a3, the valid signal length is (30)Nvalid=Lp·cLfrm·a1m.The rest of signal is approximate random and can be regarded as noise, whose length is Nseq-Nvalid. The SNRs of DPT3, DPT2, and DFrFT are(31)SNRDPT3=Nvalid2A0A18Nseq′6A6σ2+10A4σ4+4A2σ6+σ8+Nseq′-NvalidA0A18SNRDPT2=Nvalid2A0A14Nseq2A2σ2+σ4+Nseq-NvalidA0A14SNRDFrFT=Nvalid2A0A12Nseqσ2+Nseq-NvalidA0A12,respectively. The SNRs will be much lower than 0 dB when Nvalid≪Nseq, which leads to an inaccurate estimation result. Fortunately, Nvalid can be obtained by (30) with a given a1m, and Nvalid≈Nseq can be ensured by selecting the data length. Thus, we assume the pilot length is sufficient for estimation, and (31) can be simplified as follows in the high SNR regime:(32)SNRDPT3≈Nseq′A0A186A6σ2+10A4σ4+4A2σ6+σ8≈Nseq′A0A126σ2,SNRDPT2≈NseqA0A142A2σ2+σ4≈NseqA0A122σ2,SNRDFrFT≈NseqA0A12σ2.With the aid of Newton method, the accuracy of a^3, a^2, and a^1 based on DPT3, DPT2, and DFrFT is improved and given by(33)σa^32=62π2SNRDPT3Ts′′2Nseq′2-1≈36σ22π2A0A12Ts′′2Nseq′Nseq′2-1,σa^22=62π2SNRDPT2Ts′′2Nseq′2-1≈12σ22π2A0A12Ts′2NseqNseq2-1,σa^12=62π2SNRDFrFTTs′′2Nseq′2-1≈6σ22π2A0A12Ts′2NseqNseq2-1.With a given SNR, the accuracy of a^3, a^2, and a^1 can be improved by lengthening data or increasing sample interval. Assuming Nseq′=Nseq and Ts′′≈Ts′, the accuracy relationship among a^3, a^2, and a^1 is(34)σa^32≈3σa^22≈6σa^12.Note that the accuracy of a^2 can be further improved by α-search of DFrFT. Letting the estimation error before α-search and the search range of α-search be xa^2 and a^2-3σa^2,a^2+3σa^2, respectively, and assuming that the search times are sufficient, the estimation error after α-search is (35)xa^2′=xa^2+3σa^2,xa^2<-3σa^2,0,xa^2∈-3σa^2,3σa^2,xa^2-3σa^2,xa^2>3σa^2.The expectation and variance of xa^2′ are presented in the following equation:(36)Exa^2′=∫-∞-3σa^2xa^2+3σa^2·pxa^2dxa^2+∫-3σa^23σa^20·pxa^2dxa^2+∫3σa^2+∞xa^2-3σa^2·pxa^2dxa^2=0,varxa^2′=∫-∞-3σa^2xa^2+3σa^22·pxa^2dxa^2+∫-3σa^23σa^20·pxa^2dxa^2+∫3σa^2+∞xa^2-3σa^22·pxa^2dxa^2=2∫-∞-3σa^2xa^22·pxa^2dxa^2+6σa^2∫-∞-3σa^2xa^2·pxa^2dxa^2+9σa^22∫-∞-3σa^2pxa^2dxa^2=101-erf32σa^22≈0.027σa^22,with pxa^2=1/2πσa^22e-xa^22/2σa^22 and erfx=2/π∫0xe-t2dt. If the search range becomes to -2σa^2,2σa^2, the variance turns into 0.35σa^22 accordingly.Then we discuss the accuracy ofa^0, which is derived from DFrFT due to (21). In the high SNR regime, a^0 follows a zero-mean Gaussian distribution with variance(37)σa^02=σ22π2NseqA0A12sinc2⁡Δf,where Δf denotes the offset between the DFT result and the true value.Finally, the performance of SFFT is discussed as follows. LettingX^-Γk=X^k∣k∉Γ, the estimation error probability based on SFFT is [22](38)PX^k-Xk2≥εlX^-Γk22+3δ2X^k12<OlεB.SFFT performs similar to FFT with a sufficient small error probability if l≪B. Consequently, the accuracies of a^3, a^2, a^1, and a^0 based on SFFT are approximately equal to those based on FFT, which demonstrates the agreement between the proposed sparse algorithm and nonsparse algorithm.


## 4.2. Algorithm Complexity Analysis

The algorithm complexity is shown in Table1, where λitrt and λα denote the iteration times of Newton method and the α-search times, respectively. Referring to the values of ω, B, and L in Algorithm 1, the complexities of SFFT and FFT can be expressed as Olog2⁡MMllog2⁡M and OMlog2⁡M, respectively. In the large data size regime, the complexity of FFT is about M times more than that of SFFT.Table 1
The algorithm complexity.
Approach Multiplies Adds Correlation 2 f s L h T s y b 2 f s L h T s y b - 1 DP2 N s e q 0 DFT FFT M 2 log 2 ⁡ M M log 2 ⁡ M SFFT ω + B 2 log 2 ⁡ B + l M B L B log 2 ⁡ B + l M B L DFrFT DFrFT M + M 2 log 2 ⁡ M λ α M log 2 ⁡ M λ α SDFrFT M + ω + B 2 log 2 ⁡ B + l M B L λ α B log 2 ⁡ B + l M B L λ α Rife & Newton methods 3 N s e q + 1 λ i t r t + 2 2 N s e q - 1 λ i t r t + 2In the condition ofNseq=M, considering DP2 and Rife and Newton methods are employed three times, respectively, and SDFrFT and SFFT are employed twice, respectively, the total numbers of multiplications, Cmul, and additions, Cadd, in the proposed sparse algorithm are(39)Cmul≈9λitrt+2λαNseq+2λαlog2⁡NseqNseqllog2⁡Nseq,Cadd≈6λitrtNseq+2λαlog2⁡NseqNseqllog2⁡Nseq.Similarly, the total numbers of multiplications, C˙mul, and additions, C˙add, in the nonsparse algorithm are(40)C˙mul≈9λitrtNseq+λαNseqlog2⁡Nseq,C˙add≈6λitrtNseq+2λαNseqlog2⁡Nseq.Figure7 presents the complexity comparison of the proposed sparse algorithm with the nonsparse algorithm and the DDLL algorithm. The results indicate that the proposed algorithm has a much lower complexity than both the nonsparse algorithm and the DDLL algorithm.Figure 7
The comparison in terms of complexity.


## 5. Implementation Issues

In this section, we discuss the implementation issues of the proposed algorithm including the module reuse analysis and the clock rate analysis.


### 5.1. Module Reuse Analysis

In NTN platforms, with the requirement of low power consumption and low cost, the chip size is strictly limited, which means that the reuse of system module is necessary.Figure8 presents the implementation procedure of the proposed algorithm. According to Section 3, it mainly contains four steps, among which some modules can be reused. First, DP3 can be conducted by twice DP2 due to (24). Thus, SDPT3 and SDPT2 share a DP2 module. Second, SDPT and SDFrFT share an SFFT module. Third, SDPT2, SDPT3, and SDFrFT share a “Rife and Newton methods” module. Consequently, a correlation module, a DP2 module, an SFFT module, and a “Rife and Newton methods” module are sufficient for the whole procedure. We just need fsLhTsyb real multipliers, about Blog2⁡B+lM/B+2Nseq+fsLhTsyb adders, and about ω+B/2log2⁡B+l/B+1M+4Nseq complex multipliers in total, which is easy to implement on chip, e.g., on the Field Programmable Gate Array (FPGA).Figure 8
The implementation of the proposed algorithm.


### 5.2. Clock Rate Analysis

The clock rate is a crucial parameter in high speed system. A 5 Gbps communication system requires a clock rate offs≥10 GHz, which is difficult to implement on FPGA as the maximum rate of FPGA is on the order of hundreds of MHz. To achieve this target clock rate, the parallel processing is required in the conventional DDLL algorithm, which would lead to high resource consumption. In comparison, the clock rate and the degree of parallelism can be decreased in the proposed sparse algorithm thanks to the subsampling, and the power consumption is accordingly decreased to a great extent.The comparison between the DDLL algorithm and the proposed algorithm in terms of clock rate and degree of parallelism are shown in Table2 with Lfrm = 8192 bit, β=1, fs=10 GHz, and rsmp=βLfrmTsyb/Ts=16384. The numerical results demonstrate the superiority of the proposed algorithm over the DDLL algorithm in terms of power consumption. From the above discussions, the proposed algorithm can be applied to the spaceborne platform of 5G NTNs.Table 2
The comparison between DDLL and the proposed algorithm.
Items The DDLL Algorithm The proposed Algorithm Ratio Process clock rate 156.25 MHz 610.35 kHz 256:1 Degree of parallelism 64 1 64:1


## 5.1. Module Reuse Analysis

In NTN platforms, with the requirement of low power consumption and low cost, the chip size is strictly limited, which means that the reuse of system module is necessary.Figure8 presents the implementation procedure of the proposed algorithm. According to Section 3, it mainly contains four steps, among which some modules can be reused. First, DP3 can be conducted by twice DP2 due to (24). Thus, SDPT3 and SDPT2 share a DP2 module. Second, SDPT and SDFrFT share an SFFT module. Third, SDPT2, SDPT3, and SDFrFT share a “Rife and Newton methods” module. Consequently, a correlation module, a DP2 module, an SFFT module, and a “Rife and Newton methods” module are sufficient for the whole procedure. We just need fsLhTsyb real multipliers, about Blog2⁡B+lM/B+2Nseq+fsLhTsyb adders, and about ω+B/2log2⁡B+l/B+1M+4Nseq complex multipliers in total, which is easy to implement on chip, e.g., on the Field Programmable Gate Array (FPGA).Figure 8
The implementation of the proposed algorithm.


## 5.2. Clock Rate Analysis

The clock rate is a crucial parameter in high speed system. A 5 Gbps communication system requires a clock rate offs≥10 GHz, which is difficult to implement on FPGA as the maximum rate of FPGA is on the order of hundreds of MHz. To achieve this target clock rate, the parallel processing is required in the conventional DDLL algorithm, which would lead to high resource consumption. In comparison, the clock rate and the degree of parallelism can be decreased in the proposed sparse algorithm thanks to the subsampling, and the power consumption is accordingly decreased to a great extent.The comparison between the DDLL algorithm and the proposed algorithm in terms of clock rate and degree of parallelism are shown in Table2 with Lfrm = 8192 bit, β=1, fs=10 GHz, and rsmp=βLfrmTsyb/Ts=16384. The numerical results demonstrate the superiority of the proposed algorithm over the DDLL algorithm in terms of power consumption. From the above discussions, the proposed algorithm can be applied to the spaceborne platform of 5G NTNs.Table 2
The comparison between DDLL and the proposed algorithm.
Items The DDLL Algorithm The proposed Algorithm Ratio Process clock rate 156.25 MHz 610.35 kHz 256:1 Degree of parallelism 64 1 64:1


## 6. Simulation Results

In this section, we evaluate the performance of the proposed sparse algorithm in comparison with the nonsparse algorithm and the DDLL algorithm. The main parameters are listed in Table3. In the simulation, we first verify the root-mean-square error (RMSE) of these three algorithms versus data size, Nseq. As shown in Figure 9, the theoretical results derived from (33), (36), and (37) are presented by the green lines in the figures. It is observed that the simulation results match well with theoretical results.Table 3
The critical parameters.
Parameter Value Modulation OOK Demodulation Direct detection Data rate 5 Gbps SNR 14 dB The maximum velocity 7 km/s The maximum acceleration 800 m/s2 The maximum jerk 60 m/s3 L f r m 8192 bit L s y n 32 bit L p 49 bitFigure 9
The comparison of the three algorithms versusNseq.
(a)
The RMSE of distance estimation (b)
The RMSE of velocity estimation (c)
The RMSE of acceleration estimation (d)
The RMSE of jerk estimationSpecifically, Figure9(a) plots the RMSE of distance estimation (a^0) versus Nseq. Although there exists a small gap between the results of DFrFT and SDFrFT, the proposed sparse algorithm can still accurately estimate a0. Moreover, it outperforms the DDLL algorithm, especially in the small data size regime. The reason is that DDLL needs long convergence time to achieve high-precision estimation in high dynamic environment. Although the convergence rate can be improved by increasing the loop bandwidth, it results in a high steady-state error and violates the accuracy requirement.Figure9(b) plots the RMSE of velocity estimation (a^1) versus Nseq. SDFrFT exhibits a similar performance to DFrFT on estimating a1, which is similar to that on estimating a0. Furthermore, with the aid of Newton method, SDFrFT outperforms the DDLL algorithm, especially in the small data size regime. From the result, a small Lp is sufficient for estimating a1, which indicates that the proposed algorithm can be implemented without enormous channel cost.Figure9(c) plots the RMSE of acceleration estimation (a^2) versus Nseq. As shown in the figure, the accuracy of DPT2-based result is improved by Newton method and further improved by the α-search of DFrFT, which can be explained by the analysis in Section 4. The improved estimated result has a better accuracy than the DDLL-based result. In addition, in the small data size regime, there exists a larger gap between the results of the proposed algorithm and the DDLL algorithm in Figure 9(c) than that in Figure 9(b). This is because the transformation-domain algorithm is designed for estimating the acceleration, and DDLL responds more slowly in tracking the acceleration compared with tracking the velocity.Figure9(d) plots the RMSE of jerk estimation (a^3) versus Nseq. The comparison among DPT3, SDPT3, and DDLL indicates that there are a huge gap and a minor gap between the results of the proposed algorithm and the DDLL algorithm and between the results of the proposed sparse algorithm and the nonsparse algorithm, respectively. Moreover, in the small data size regime, in comparison to the acceleration estimation results in Figure 9(c), there exists a larger gap between the results of the proposed algorithm and the DDLL algorithm in terms of the jerk estimation. The reason is that SDPT3 is designed for estimating the jerk, and DDLL responds more slowly in tracking the jerk than the acceleration, which presents the advantage of the proposed algorithm in high dynamic environment.Next, we focus on the comparison between the sparse algorithm and nonsparse algorithm versus parameterβ in the formula rsmp=βLfrmTsyb/Ts with Nseq=2048. The simulation results are plotted in Figure 10 and match well with the theoretical results.Figure 10
The comparison of the sparse algorithm and nonsparse algorithm versusβ with Nseq=2048.
(a)
The RMSE of distance estimation (b)
The RMSE of velocity estimation (c)
The RMSE of acceleration estimation (d)
The RMSE of jerk estimationSpecifically, Figure10(a) presents the RMSE of distance estimation (a^0) versus β. There exists no pronounced difference between the results of SDFrFT and of DFrFT, which agrees with the results in Figure 9(a). However, the variation tendency of the results in Figure 10(a) is different from that in Figure 9(a). The RMSE decreases with the increase of Nseq in Figure 9(a) and remains constant no matter the value of β in Figure 10(a). It can be explained by (37) that the accuracy of a^0 is related to data size and is independent of the sample interval.Figures10(b) and 10(c) plot the RMSEs of velocity estimation (a^1) and acceleration estimation (a^2) versus β, respectively. The results indicate that the performance is almost the same no matter the sparse method is used or not, which verifies the rationality of the sparse algorithm. Compared with lengthening data, increasing sample interval only has a slighter effect on estimating a1 and a2 in the condition of the same process duration. It can be explained by (33) that lengthening data contributes more to the accuracy than increasing sample interval. Fortunately, the difference between these two cases is so small that an accurate estimation still can be achieved by increasing the interval. As increasing sample interval improves the estimation accuracy without increasing computational burden, we prefer it rather than lengthening data for the sake of power consumption.From the results in Figures9(d) and 10(d), lengthening data performs better than increasing sample interval in the same process duration, which is similar to the comparison results of Figures 9(b), 9(c), 10(b), and 10(c). The gap can be filled up by slightly increasing the data size, e.g., as shown in Figure 11 with Nseq=4096, whose results indicate that the configuration of Nseq=4096 and β=16 exhibits a prior performance to the configuration of Nseq=32768 and β=1 in Figure 9(d), and a much better performance than that of Nseq=2048 and β=16 in Figure 10(d).Figure 11
The jerk estimation versusβ with Nseq=4096.From the above discussion, as increasing sample interval improves the estimation accuracy without increasing computation resource, it is a better choice than lengthening data in resource constrained platforms.


## 7. Experiment Results

Referring to the parameters listed in Table3, we have built a hardware platform for the experiment to verify the implementability of the proposed algorithm. The topology of hardware platform is given in Figure 12. The laser signal is received by telescope and is transformed to an electrical signal in the optical receiver. After that, the electrical signal is transmitted to channel simulator, which is employed to tune the transmission delay a0 and delay rate a1 with the control of PC, via the LAN extension for instrumentation (LXI) protocol. The output of channel simulator is sampled by an analog-digital converter (ADC). Then FPGA utilizes the samples to implement the proposed temporal synchronization algorithm according to Figure 8. The synchronization results are reported to PC, via the peripheral component interconnect express (PCI-E) protocol.Figure 12
The topology of hardware platform.We tune the transmission delay and delay rate by channel simulator and estimate them by the proposed algorithm withM=Nseq=2048 and β=2. The experiment is repeated 50 times for each configuration of a0 and a1, and the statistical results are presented in Table 4, where the RMSEs of transmission delay and delay rate are on the orders of 0.01 ns and 0.3 ns/s, respectively. Due to the performance loss caused by the nonideal hardware and the nonideal channel, the experimental result performs inferior to the simulation result. Fortunately, the accuracy of experimental result still meets the demand of 5G NTNs, and the performance can be further improved by lengthening the process duration.Table 4
The experiment results of temporal synchronization.
True value ofa0 (ns) True value ofa1 (ns/s) RMSE ofa^0 (ns) RMSE ofa^1 (ns/s) 1 100 0.0103 0.264 1 1000 0.0100 0.311 10 100 0.0168 0.192 10 1000 0.0270 0.220 100 100 0.0138 0.268 100 1000 0.0133 0.306 1000 100 0.0185 0.212 1000 1000 0.0265 0.109Next, we repeat the proposed temporal synchronization algorithm once per 50 milliseconds to continuously estimate the transmission delay, and the result is depicted in Figure13. The result in Figure 13 indicates that the proposed algorithm can be implemented to continuously estimate the temporal parameters in real time, and the accuracy achieves about 0.01 ns in experimental environment. As the employed channel simulator cannot simulate the acceleration and jerk currently, the performance on estimating a2 and a3 will be verified in follow-up experiments.Figure 13
The experiment result of continuous estimation.


## 8. Conclusions

In this paper, we have developed a temporal synchronization algorithm based on the sparse pilot and the sparse transform method in the scenario with high Doppler and high data rate. The performance of the proposed algorithm has been analyzed in comparison with the conventional DDLL algorithm and the nonsparse algorithm. The analytical and simulation results demonstrate that the proposed algorithm performs best in terms of a good accuracy and the lowest complexity. In addition, we have implemented the proposed algorithm, which confirms its implementability in resource constrained platforms, e.g., the IMDD-based laser communications for feeder links in 5G NTN platforms.



---

*Source: 8284617-2018-06-20.xml*