Harmonic Retrieval of CFO and Frame Misalignment for OFDM-based Inter-Satellite Links

The standardization activities in the integrated 5G networks and low earth orbits (LEO) mobile satellite communication (SatCom) systems and the emerging plans about LEO constellations, introduce a renewed research interest in the performance of inter-satellite links (ISLs) [1, 2, 3]. With the technological advances in small satellites, such as cubesats, ISLs will serve as the essential link to enable the satellites to work together to accomplish the same task and satellite formation flying missions [4, 5, 6].

Orthogonal frequency division multiplexing (OFDM) emerges as the leading waveform candidate for ISLs due to the inherent advantages to provide high data rates at relatively low complexity. OFDM is a multicarrier modulation technique widely adopted by modern mobile communication systems such as WiFi networks (IEEE 802.11a-802.11be), Long Term Evolution (LTE), and LTE-Advanced. It is also considered as the basic waveform for the 5G New Radio (NR). Although OFDM has very desirable traits including high spectral efficiency, achievable high data rates, and robustness in the presence of multipath fading, it is also very sensitive to synchronization errors such as the carrier frequency offset (CFO) and frame misalignments [7]. Due to the high speed movement of the LEO satellites, synchronization is expected to be a major challenge in ISLs with the high Doppler shifts [1, 8]. In the presence of a high Doppler shift, the subcarriers of an OFDM system with CFO may no longer be orthogonal and the resulting inter-carrier interference (ICI) may degrade the performance of the system. Another synchronization issue arises when the starting position of the discrete Fourier transform (DFT) window at the receiver is misaligned. Depending on the position of the misalignment, the effects of the frame misalignment may range from a simple phase offset to ICI.

Although the research on carrier and frame timing synchronization for OFDM systems is very mature [9, 10, 11, 12, 13, 14], in this paper, we aim to provide a fresh perspective to OFDM’s synchronization problem by also aiming to combat the residual Doppler shifts that may be encountered in ISLs due to their high mobility. Estimating the Doppler shifts is not only important for dealing with CFO effects but the Doppler shifts can also be used for securing the ISLs between spacecrafts, by using them as a mutually shared secret [15].

The proposed method transforms the estimation of the CFO and the frame misalignment in an OFDM based inter-satellite link into a 2-D harmonic retrieval problem. The representation is in the frequency domain and this allows the estimation of the frame misalignment unlike [11, 12] where the frame misalignment is assumed to be corrected prior. When compared to the methods that utilize the cyclic prefix (CP) portion of the OFDM symbol [9, 14], the proposed method relies on the use of pilot symbols. Furthermore, the proposed method does not require the knowledge of the noise statistics unlike the methods that use the CP [9, 14]. We also derive the Cramér-Rao lower bound (CRLB) for the joint estimation of the CFO and frame misalignment. Numerical results show that for the error range of $[10^{-2},10^{-4}]$, the difference between the proposed 2-D ESPRIT based method and the PSS method [14] is less than 5dB at its worse.

The contributions of the paper can be summarized as follows:

• •

  The proposed method represents the estimation of the CFO and the frame misalignment in an OFDM-based ISL as a 2-D harmonic retrieval problem. Unlike [11], the representation is in the frequency domain and this allows the estimation of the frame misalignment.

• •

  When compared to the methods that utilize the cyclic prefix symbols [9, 14], the proposed method relies on the use of pilot symbols. So the length of the CP symbols is not a factor that affects the estimation performance of the proposed method.

• •

  The proposed method does not require knowledge of noise statistics, unlike the methods that use cyclic prefix symbols

• •

  A disadvantage of the proposed method is that due to relying on sending constant pilot symbols, the peak-to-average ratio is high. [9, 14].

The organization for the rest of this paper is in the following way; Section 2 introduces the signal model for the OFDM-based ISLs. Section 3 reformulates the signal model as a 2-D estimation problem. Section 4 shows the CRLB of the parameters in the reformulated model. We compare the simulation results of the proposed method against the state-of-art methods in Section 5. Finally, Section 6 presents the conclusions and directions for future work.

In an OFDM system, the symbols are transmitted with the sampling period $T_{s}$ in a series of frames denoted by $X_{m}[k]$ where $m$ indicates the $m$-th OFDM frame and $k$ is the subcarrier index. The time-domain symbols are modulated by applying the inverse DFT to the frequency-domain symbols for $k=0,\ldots,N-1$ where $N$ is the total number of subcarriers. Then $N_{g}$ number of CP samples are appended to the front of the time-domain frame as in

where $0\leq n\leq N_{t}-1$ and $N_{t}=N+N_{g}$. The discrete-time received frames in baseband can be written as

where $\circledast$ denotes the cyclic convolution operation, $h_{m}[n]$ is the overall channel that is the result of the convolution of the multipath fading channel and the pulse shaping filters of both the transmitter and the receiver, the starting position of the DFT window at the receiver is shown as $\tau_{p}x_{m}[n]=x_{m}[n+p]$, and the additive white Gaussian noise (AWGN), $z_{m}[n]$, is modeled as a circularly symmetric Gaussian random variable, i.e. $z_{m}[n]\sim\mathcal{CN}(0,N_{0})$. Due to the clock inaccuracies of the transmitter and the receiver oscillators, and the Doppler spread caused by the mobility of the satellites, the conversion from passband to baseband generates the unwanted multiplicative terms, $e^{i2\pi\epsilon mN_{t}/N}e^{i2\pi\epsilon n/N}$, in (2) where $\epsilon$ is the total CFO term normalized by the subcarrier spacing, $1/NT_{s}$, as

$\Delta_{f_{c}}$ and $f_{d}$ in (3) shows the clock difference and the Doppler spread, respectively. Since the free-space loss and thermal noise of the electronics are enough to characterize the ISLs, an AWGN channel model, i.e. $h_{m}[n]=\delta_{0}[n]$, can be used in (2) [5]. Thus the received signal for the discrete baseband equivalent model can be written as

where $\alpha=N_{g}/N$. The normalized CFO, $\epsilon=\varepsilon+\ell$, consists of a fractional part, i.e. $|\varepsilon|\leq 0.5$, and an integer part $\ell\in\mathbb{Z}^{+}$. The OFDM receiver first removes the cyclic prefix samples and then applies the DFT to the remaining samples resulting in

where $C_{m}[k]$ is given as

and $Z_{m}[k]$ is also circularly symmetric Gaussian that is $Z_{m}[k]\sim\mathcal{CN}(0,NN_{0})$ since the DFT is a linear transformation and the circular symmetry is invariant to linear transformations.

The DFT of the received samples, $R_{m}[k]$, in (5) can be rewritten as a 2-D signal model by using the OFDM frame index, $m$, as a second dimension taking values in the range $0\leq m\leq M-1$

where the frequencies and the phase terms are respectively

The transmitted OFDM symbols in (7), $X[m,k-\ell-r]$, depend on both the frame index, $m$, and the subcarrier index, $k$. This dependency can be removed by sending the same pilot symbol on each subcarrier that is $X[m,k]=X$ for all $m$ and $k$, and then the DFT of the received samples (7) can be formulated as a harmonic retrieval problem

where $\phi=e^{i2\pi f_{1}}$, $\theta=e^{i2\pi f_{2}}$ and $c$ is a complex coefficient

The signal in (11) has a single 2-D mode defined by the frequencies $\{f_{1},f_{2}\}$, and a complex coefficient $c=\lambda e^{i\varphi}$ where $\lambda=|c|$ and $\varphi=\angle c$. The ESPRIT-based estimation methods identify the pairs $\{f_{1},f_{2}\}$ from the observed data $R[m,k]$ by turning the 2-D estimation problem into two 1-D estimation problems and exploiting the shift-invariance structure of the signal subspace for each 1-D problem [16]. Since the unknown 2-D signal mode is undamped, the forward-backward prediction [16] can be applied to increase the estimation accuracy by using an extended data matrix $\mathbf{R}_{ee}$ as

where complex conjugation without transposition is shown by the $(\cdot)^{*}$ symbol, and $\mathbf{\Pi}$ is a permutation matrix with ones on its antidiagonal and zeroes elsewhere. $\mathbf{R}_{e}$ (13) is the enhanced Hankel block structured matrix that is constructed by applying an observation window of size $P\times Q$ through the rows of the noisy received samples of which one of the dimensions is fixed. $\mathbf{R}_{e}$ matrix for the first dimension is given as

where $\mathbf{R}_{(m)}$ is a Hankel matrix of size $Q\times(N-Q+1)$ as given below

The extended data matrix constructed according to (13) for the first dimension, $\mathbf{R}_{ee1}$, can be decomposed in terms of a signal and a noise subspace as

where $\mathbf{Z}_{1}$ is the Hankel block structured matrix constructed from the noise samples $Z[m,k]$ in the same way as $\mathbf{R}_{e1}$ from $R[m,k]$. $\mathbf{s}_{L1}$ is a vector of size $PQ\times 1$ given as $\mathbf{s}_{L1}=[\boldsymbol{\theta}_{1}\,\boldsymbol{\theta}_{1}\phi\,\ldots\boldsymbol{\theta}_{1}\phi^{P-1}]^{T}$ where $\boldsymbol{\theta}_{1}$ is $\boldsymbol{\theta}_{1}=[1\,\,\theta\,\ldots\,\theta^{Q-1}]^{T}$. The construction of the vector $\mathbf{s}_{R1}$ is similar to that of $\mathbf{s}_{L1}$. The rank of the signal subspace, $c\,\mathbf{s}_{L1}\mathbf{s}^{T}_{R1}$, must be equal to the number of signal components which is equal to one since there is only one user. In order for the rank of the signal part be at least equal to one, $P$ and $Q$ must satisfy the following inequalities [16]

The singular value decomposition (SVD) of $\mathbf{R}_{ee1}$ yields

While the singular vectors and the singular value of the signal mode is contained respectively in $\mathbf{U}_{S1}$, $\mathbf{V}^{H}_{S1}$, and $\mathbf{\Sigma}_{S1}$, the singular vectors and the singular values of the remaining noise components are respectively in $\mathbf{U}_{Z1}$, $\mathbf{V}^{H}_{Z1}$, and $\mathbf{\Sigma}_{Z1}$ (18).

The frequencies related to the first dimension, $f_{1}$, can be estimated the $\mathbf{F}_{1}$ matrix

where $\underline{\mathbf{U}}_{S1}$ (resp. $\overline{\mathbf{U}}_{S1}$) is constructed with the first (resp. last) $(P-1)Q$ rows of the matrix $\mathbf{U}_{S1}$. $\underline{\mathbf{U}}^{\dagger}_{S1}$ denotes the pseudo-inverse matrix of $\underline{\mathbf{U}}_{S1}$ and the total least squares can be applied to calculate $\mathbf{F}_{1}$ in (19). The frequency $f_{1}$ is the eigenvalue of the matrix $\mathbf{F}_{1}$ that is

where $\mathbf{T}_{1}$ is the eigenvector matrix that diagonalizes $\mathbf{F}_{1}$. For the frequency related to the second dimension, $f_{2}$, the $\mathbf{F}_{2}$ matrix can be obtained by using the extended data matrix for the second dimension, $\mathbf{R}_{ee2}$, and following the same steps. Finally the pairing method of [16] must also be employed to calculate the eigenvalue decomposition of a linear combination of $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$

where $\beta$ is a scalar. The diagonalizing transformation $\mathbf{T}$ is applied to both $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$

which yields the ordered frequencies.

The unknown parameters of (11) can be collected in a vector of size $4\times 1$ as

where $\omega_{1}=2\pi f_{1}$ and $\omega_{2}=2\pi f_{2}$ represent the angular frequencies. If the DFT of the received signal samples (11) are written as an $MN\times 1$ column vector

then the joint probability density function (PDF) of the multivariate circularly symmetric Gaussian random vector $\mathbf{r}\sim\mathcal{CN}(\boldsymbol{\mu},NN_{0}\mathbf{I})$ is given as

where $\mathbf{I}$ is an identity matrix of size $MN\times MN$. The $j$-th entry of the mean vector $\boldsymbol{\mu}$ (26) is

where

The $(j,j^{\prime})$ entry of the Fisher information matrix for the multivariate Gaussian PDF (26) is

where $\Re\{\cdot\}$ takes the real parts of the entries of the matrix [17]. The diagonal entries of the inverse of the Fisher matrix are the CRLB of the parameters. The CRLB for $\omega_{1}$ and $\omega_{2}$ are written respectively in (31) and (32).

The simulated inter-satellite communication system uses the OFDM modulation which consists $N=64$ subcarriers and $N_{g}=16$ CP samples. The CFO is generated independently for each OFDM frame as a uniform random variable $\epsilon\sim\mathcal{U}[0.2,0.25]$ and the frame misalignment is fixed at $p=2$. The signal-to-noise (SNR) ratio is defined as $\text{SNR}=E_{s}/N_{0}$ where the symbol energy is normalized to unity, $E_{s}=1$. The performance of the proposed method is measured by calculating the mean squared error (MSE) for both the CFO and frame misalignment estimations and a total of 2000 Monte Carlo simulations are run for each SNR value.

The selected methods are the maximum likelihood estimator, also known as Beek’s method [9] that allows the transmission of data symbols by using the CP samples for estimation, Minn’s improved Schmidl & Cox estimator [10], and the cross-correlation based method that utilizes primary synchronization symbols (PSS method) [14]. The parameters for the proposed method are as follows; the number of the pilot symbols is $M=2$, the size of the observation window applied to the received samples is chosen according to (17) as $P=Q=2$ and the pairing coefficient is $\beta=8$. The results of the frame misalignment estimation (Figure 1) show that while Beek’s method and Minn’s method perform worse than the proposed method at each SNR value, the performance of the PSS method surpasses the performance of the proposed method.

The MSE regarding the CFO estimation results is shown in Figure 2. The MSE of the proposed 2-D ESPRIT method is closest to the CRLB bound at 0 dB SNR and for the SNR values less than 10 dB, the proposed method performs superior compared to the rest of the methods. The PSS method is the second-best performer for SNR values less than 5 dB.

The computational complexities of all the methods and their numerical evaluations for the parameter values used in the simulations are in Table I. The complexities are in terms of real floating-point operations (flops). One complex multiplication is counted as 6 real flops while one complex addition is counted as 2 real flops. The SVD step determines the complexity of the proposed 2-D ESPRIT estimation method and the complexity is written in terms of the size of the extended data matrix, $\mathbf{R}_{ee}\in\mathbb{C}^{L\times K}$, where $L=PQ=4$ and $K=2(N-Q+1)(M-P+1)=126$. The symbol $G$ used in the complexity of the PSS method shows the size of the search grid that is used in the CFO estimation step. $D$, the number of repetition of the training symbol used in Minn’s method, is $D=4$. In the numerical results, $M$ and $G$ are chosen as $M=2$ and $G=500$ respectively. While the complexity of the proposed method is less than the PSS’s method, the complexities of both Beek’s and Minn’s methods are lower.

Synchronization in the ISLs of the LEO systems is a critical issue due to the Doppler spread caused by the high mobility of the satellites. We presented a novel method for the estimation of the CFO and frame misalignment in OFDM-based ISLs. We reformulate the synchronization problem as a 2-D harmonic retrieval problem in the frequency domain and apply the 2-D ESPRIT method to estimate the parameters. Since this new approach is in the frequency domain unlike the previous harmonic retrieval approach, this allows the joint estimation of the CFO and the frame misalignment.

When compared to other synchronization methods that also rely on the pilot symbols like the well-known Schmidl & Cox method and the PSS method, the proposed method requires transmitting frames of pilot symbols for the reformulation to work. The CRLB for the joint estimation of the CFO and the frame misalignment is derived and the performance of the proposed method is compared against the Beek’s, Schmidl & Cox, and the PSS methods. Numerical results show that respectable estimation performance can be achieved by using the proposed method with only two consecutive pilot frames.