Maximum Likelihood Timing and Carrier Synchronization in Burst-Mode Satellite Transmissions

This paper investigates the joint maximum likelihood (ML) estimation of the carrier frequency offset, timing error, and carrier phase in burst-mode satellite transmissions over an AWGN channel. The synchronization process is assisted by a training sequence appended in front of each burst and composed of alternating binary symbols. The use of this particular pilot pattern results into an estimation algorithm of affordable complexity that operates in a decoupled fashion. In particular, the frequency offset is measured first and independently of the other parameters. Timing and phase estimates are subsequently computed through simple closed-form expressions. The performance of the proposed scheme is investigated by computer simulation and compared with Cramer-Rao bounds. It turns out that the estimation accuracy is very close to the theoretical limits up to relatively low signal-to-noise ratios. This makes the algorithm well suited for turbo-coded transmissions operating near the Shannon limit.


INTRODUCTION
Burst transmission of digital data and voice is widely adopted in satellite time-division multiple-access (TDMA) networks. In these applications the propagation medium can be reasonably modeled as an additive white Gaussian noise (AWGN) channel and knowledge of carrier frequency, symbol timing, and carrier phase is necessary for coherent demodulation of the received waveform. In the presence of nonnegligible phase noise and/or oscillator instabilities, differential detection is often employed to overcome the inherent difficulty posed by the phase estimation process. Even with differential detection, however, the problem of timing and frequency offset recovery still remains.
Depending on their topology, synchronization circuits can be divided into two main categories: feedback and feedforward schemes [1,2]. The former have good tracking capabilities but exhibit comparatively long acquisitions due to hang-up phenomena [3][4][5]. The latter have shorter acquisitions and, accordingly, are better suited for burst-mode transmissions. In many cases, a preamble of known symbols is appended at the beginning of each burst to assist the synchronization process. Actually, the use of a preamble allows data-aided (DA) operation and provides better estima-tion accuracy as compared to a non-data-aided (NDA) approach. Even so, however, synchronization may prove difficult, especially with turbo-coded modulations operating at relatively low signal-to-noise ratios (SNRs). Clearly, very efficient synchronization algorithms are needed in these conditions [6].
The common approach to solve the synchronization problem in burst-mode transmissions is to estimate the timing error first, and then use the time-synchronized samples for frequency and phase recovery. Two prominent feedforward schemes for NDA timing estimation are investigated in [7,8]. In particular, timing estimates are derived in [7] by searching for the maximum of an approximate version of the likelihood function while in [8] the received signal is sampled at some multiple of the symbol rate and a squarelaw nonlinearity (SLN) is employed to wipe the modulation out. As shown in [2], the method in [8] is an efficient way of maximizing the likelihood function of [7] as long as the bandwidth of the complex envelope of the transmitted signal does not exceed the signaling rate. Since the use of a SLN exhibits poor performance in the presence of narrowband signaling, alternative methods employing absolute value or fourth order-based nonlinearities have been devised [9]. The main advantage of the timing estimators in [7][8][9] is that they can operate correctly even in the presence of carrier frequency offsets (CFOs) as large as 20% of the symbol rate.
Frequency estimation is usually performed by exploiting the received time-synchronized samples. A large number of schemes proposed in the past operate in either the frequency or time domain. The Rife and Boorstyn (R&B) algorithm [10] belongs to the former class and provides maximum likelihood (ML) estimates of frequency and phase errors by looking for the peak of a periodogram. Interpolation techniques may be employed to find an explicit expression of the peak location [11]. In the time-domain approach, suitable correlations of the received samples are exploited to compute the frequency estimates. A representative selection of schemes derived along this line of reasoning can be found in [12][13][14][15]. These methods attain the Cramer-Rao lower bound (CRB) at intermediate/high SNRs, but exhibit different performance in terms of estimation range and threshold, that is, the SNR below which large estimation errors are likely to occur.
A possible drawback of conventional frequency estimation schemes as those discussed in [10,[12][13][14][15] is that they all assume ideal timing synchronization. Their performance is thus limited by the accuracy of the timing estimator. A DA algorithm for the joint estimation of the carrier phase, frequency offset, and timing error has been proposed in [16] by resorting to ML arguments. In order to work properly, however, the demodulated signal must incur negligible phase rotations during the preamble duration. This poses a stringent limit to the maximum tolerable CFO, which may prevent the application of this method to many practical situations.
In the present paper we are concerned with the joint estimation of all synchronization parameters for a burst-mode satellite system operating over AWGN channels. Since one distinct feature of packetized transmissions is that synchronization must be achieved as fast as possible, in the following we only focus our attention to a feedforward structure. Also, we assume that a preamble of alternating binary symbols is transmitted at the beginning of each burst to facilitate the timing estimation task [17]. Our approach is based on ML methods and leads to a three-step procedure. In the first step frequency recovery is accomplished through a mono-dimensional grid search. The estimated CFO is then exploited in the second step to obtain a closed-form expression of the timing estimate. The final step is devoted to phase estimation and can be skipped in case of differential data detection. Surprisingly, no complicated multidimensional searches are needed to jointly estimate all the unknown synchronization parameters. Simulations indicate that the proposed estimator is well suited for turbo-coded transmissions since its accuracy approaches the relevant CRBs even at low SNR values. However, it should be observed that this advantage is achieved at the price of a higher computational complexity as compared to other existing alternatives.
The paper is organized as follows. In Section 2 we introduce the signal model and formulate the synchronization problem. Section 3 illustrates the joint ML estimation of the unknown parameters and discusses in some detail the practical implementation of the frequency estimator. In Section 4 we derive CRBs to characterize the ultimate accuracy of fre-quency, timing, and phase estimates. Simulation results are presented in Section 5 while some conclusions are drawn in Section 6.

Statement of the problem
We consider the reverse link of a satellite communication system and assume a time-division multiple-access (TDMA) scheme where each earth station transmits bursts of data. The structure of a burst is detailed in Figure 1. Essentially, it consists of two parts: a header section followed by a payload. The header is further divided in two portions, namely, a synchronization preamble and a unique word (UW). The preamble is made of a sequence of training symbols which are exploited by the receiver for carrier and symbol timing recovery. The UW is located just after the preamble and is used for burst identification as well as to establish the start of the payload.
The first task of the receiver is the start of burst (SoB) detection, that is, the recognition of the time-of-arrival (ToA) of a generic burst. This is normally performed through a simple noncoherent energy-detection scheme which provides a coarse estimate of the position of each burst. Once the SoB has been identified, the preamble is exploited for carrier and symbol timing synchronization. This is the second task of the receiver and represents the focus of our paper. In order to explain how synchronization can be achieved, we concentrate on a single burst and assume that the SoB detection algorithm has provided a ToA estimate with an error τ, as shown in Figure 2. The offset τ can be decomposed as follows: where T is the symbol period, η is an integer (integer delay), and ε ∈ [−0.5, 0.5) is a real-valued parameter (fractional delay). During the preamble we are interested in the estimation of the fractional delay, because the integer delay is recovered later by searching for the location of the UW within the burst. The estimation of the synchronization parameters (fractional delay, carrier phase, and frequency offsets) is performed by observing a portion of the preamble of length NT (N is a design parameter) at the right of the assumed SoB, as shown in Figure 2. Clearly, the total duration of the preamble has to be larger than τ + NT. Since τ is a random variable, this condition can be practically met by a proper design of the preamble length. Since we are not concerned with the estimation of the integer delay, in the following we set η = 0.

Signal model
We consider a linearly modulated digital signal transmitted over an AWGN channel. The complex envelope of the received waveform is modeled as where s(t) is the useful signal, ϕ and f d are the carrier phase and frequency offset, respectively, and w(t) is thermal noise with independent real and imaginary components, each with two-sided power spectral density N 0 . Signal s(t) is expressed as where {a n } are modulation symbols taken from a PSK or QAM constellation and g(t) has a root-raised-cosine Fourier transform with some roll-off α. To facilitate the timing estimation process, during the preamble we assume a pilot pattern composed by alternating BPSK symbols +1 and −1 [17]. Accordingly, s(t) is given by with E s denoting the signal energy per symbol interval, and r(t) may be rewritten in the form In order to produce a discrete-time signal, the received waveform is fed to an anti-aliasing filter (AAF) and sampled at some rate f c . The filter bandwidth B AAF and the sampling rate are chosen such that the signal component is passed undistorted (even for the maximum frequency offset) and no aliasing occurs. Assuming that the CFO is less in magnitude than 0.5/T, from (5) it follows that we can set B AAF = 1/T and f c = 2/T. For simplicity, in the ensuing discussion the AAF is assumed with a brick-wall transfer function, even though the rectangular shape is not strictly necessary and could be easily made more realistic [18].
For normalization purposes, the output of the AAF is scaled by a factor T/2E s and we call x(k) the correspond-ing samples taken at t = kT/2, with 0 ≤ k ≤ 2N − 1. As the signal is not distorted in passing through the filter, we have where ν = f d T is the CFO normalized to the symbol-rate 1/T and n(k) = n R (k) + jn I (k) is the noise contribution. Due to the previous hypotheses, {n R (k)} and {n I (k)} are independent and white random sequences with the same variance As the signal component in (6) depends on ν, ε, and ϕ, we may estimate all these parameters from the observation of {x(k)}. This problem is addressed in the next section by resorting to ML methods.

Maximization of the likelihood function
Bearing in mind (6), the log-likelihood function for the unknown parameters is given by where ν, ε, and ϕ are trial values of ν, ε, and ϕ, respectively. The joint ML estimate of (ν, ε, ϕ) is the location where Λ( ν, ε, ϕ) achieves its global maximum. Skipping irrelevant factors and additive terms independent of ( ν, ε, ϕ), it turns out that Λ( ν, ε, ϕ) may equivalently be replaced by where e{·} denotes the real part of the enclosed quantity. Function Ψ( ν, ε, ϕ) can also be rewritten as To ease the search for the maximum of Ψ( ν, ε, ϕ), we rewrite (9) in the form where Z( ν, ε) is a function of ν and ε defined as while ψ( ν, ε) = arg{Z( ν, ε)} is the argument of Z( ν, ε). Clearly, for fixed ν and ε, the maximum of Ψ( ν, ε, ϕ) is achieved when the cosine factor in (11) is equal to unity, which occurs for In this case the right-hand side of (11) reduces to |Z( ν, ε)| and the ML estimates of ν and ε are found by maximizing the following function: with respect to ν and ε, where the factor 2 in the right-hand side of (14) has only been inserted to avoid a factor 1/2 in the subsequent equations.
To proceed further, we substitute (12) into (14) and obtain where A( ν) is defined as Observing with θ( ν) = arg{A( ν)}. For a given ν, the maximum of Γ( ν, ε) is achieved by setting Substituting this result into the right-hand side of (17) yields from which it follows that the ML estimate of the frequency offset is given by while the timing estimate is obtained from (18) in the form In case of coherent detection, an estimate of the carrier phase ϕ is also necessary. This is computed as indicated in (13) after replacing ( ν, ε) by ( ν, ε) and reads In the sequel the algorithm based on (20)-(22) is called the ML estimator (MLE).

Remarks
(1) Contrarily to what one might fear, the maximization of the likelihood function Λ( ν, ε, ϕ) needs not be made on a three-dimensional domain. Actually, the location ( ν, ε, ϕ) of the maximum can be found through simple steps, each involving a single synchronization parameter. In particular, the first step requires maximizing the function P( ν) defined in (19) in order to get the CFO estimate ν. As discussed later, this can be done through a grid search over the interval where ν is expected to lie. Once ν has been obtained, timing and phase estimates are computed in closed form as indicated in (21) and (22), respectively. In summary, the difficult and time-consuming part in the estimation of (ν, ε, ϕ) is the one that locates the maximum of P( ν). Once this is done, the computation of ε and ϕ becomes a trivial task.
(2) Maximizing function P( ν) may pose some difficulty due to the presence of many local maxima. This is clearly evident from Figure 3, which illustrates a typical realization of P( ν) as obtained by simulation with N = 64, E s /N 0 = 10 dB, and ν = 0.1. As discussed in [10], the global maximum can be sought in two steps. The first one (coarse search) calculates P( ν) for a set of ν values, say { ν n }, covering the uncertainty range of ν and determines the location ν M of the maximum over this set. The second step (fine search) makes an interpolation between the samples P( ν n ) and computes the local maximum nearest to ν M . It should be noted that the shape of P( ν) is occasionally so badly distorted by noise that its highest peak may be far from the true ν. When this happens, the MLE makes large errors (outliers) and the system performance is highly degraded. The SNR below which the outliers start to occur is referred to as the threshold of the estimator.
(3) In practice the coarse search can be efficiently performed using fast Fourier transform (FFT) techniques, as it is now explained. Starting from the observed samples {x(k)}, we first compute the following zero-padded sequences of length KN: where K is a design parameter called pruning factor. Next, the FFTs of {y e (k)} and {y o (k)} are evaluated at the points This produces the quantities {Y e ( ν n )} and {Y o ( ν n )}, which are next exploited to get {P( ν n )} as indicated in (19). Finally, the largest P( ν n ) is sought and this provides the coarse frequency estimate.
(4) Collecting (10) and (19), it is seen that P( ν) is periodic of unit period. This means that MLE gives unambiguous frequency estimates as long as ν is confined within the interval [−1/2, 1/2). This is the frequency estimation range of MLE.
(5) Compared to the R&B algorithm [10], the MLE is more complex to implement as it requires the computation of two FFTs instead of a single FFT. In addition to carrier synchronization, however, the MLE also provides timing recovery. Actually, from (16) and (21) we see that computing the timing estimate only requires knowledge of Y e ( ν) and Y o ( ν). Since these quantities can easily be obtained by {Y e ( ν n )} and {Y o ( ν n )} through interpolation, the timing estimation task is accomplished with a relatively low computational effort once ν is available. However, it should be observed that the complexity associated to the synchronization process is negligible compared to that of iterative data decoding [19]. So, the requirement for an additional FFT has only a marginal impact on the overall receiver complexity.

CRB ANALYSIS
By invoking the asymptotic efficiency property of the MLE, we expect that the accuracy of the estimates (20)-(22) approaches the corresponding CRBs for relatively large values of N and E s /N 0 . For this reason, it is of interest to derive the CRB for the joint estimation of the set of parameters η = (ν, ε, ϕ) based on the model (6).
We begin by computing the entries of the Fisher information matrix F. They are defined as [20] [ where Λ(η) is the log-likelihood function in (7) and η denotes the th entry of η. Substituting (7) into (25), after some manipulations we obtain The CRB for the estimation of η is given by [F −1 ] , . Skipping the details, it is found that Interestingly, for large data records we can approximate (27) as which represents the CRB for the estimation of the frequency of a complex sinusoid embedded in AWGN [10].

SIMULATION RESULTS
In this section we report on simulation results illustrating the performance of MLE over an AWGN channel. Unless otherwise specified, the synchronization parameters vary at each new simulation run and are modeled as statistically independent random variables with a uniform distribution. In particular, ν and ε are confined within [−0.5, 0.5) while ϕ takes values in the interval [−π, π). A pruning factor K = 4 is used to compute the quantities {Y e ( ν n )} and {Y o ( ν n )}. Also, a parabolic interpolation is chosen in the implementation of the fine search. This yields a frequency estimate in the form where δ ν = 1/KN is the distance between two adjacent samples P( ν n ) while ν M is the output of the coarse search. The observation length is fixed to either N = 32 or 64. For comparison, in the ensuing discussion we also consider a synchronization scheme in which timing recovery is first accomplished by resorting to the Oerder and Meyr (O&M)  marks. We see that MLE has the best accuracy, especially at low SNRs where a significant gain is observed with respect to O&M. In particular, for N = 64 the MLE is close to the CRB down to E s /N 0 values of 0 dB, while O&M approaches the bound only for E s /N 0 > 10 dB. This feature of the MLE is of great importance as it makes this estimator suitable for turbo-coded modulations operating at very low SNRs. Figure 8 illustrates the accuracy of the frequency estimates provided by MLE and R&B with either N = 32 or 64. Again, the simulation results are compared with the relevant CRBs. We see that the estimation accuracy keeps close to the CRB down to a certain value of E s /N 0 that depends on the adopted scheme and observation length. If the SNR is decreased further, a rapid increase in the MSE is observed. The abscissa at which the slope of the curve starts to change indicates the estimator threshold and is a manifestation of the occurrence of outliers. Since large errors have disabling effects on the system performance, the frequency estimator must operate above threshold. The results in Figure 8 reveal that MLE has a lower threshold than R&B, especially for N = 64, which translates into an increased robustness against outliers. This fact reinforces the idea that for low SNR applications the MLE is more efficient than other conventional synchronization schemes. As expected, the threshold is a decreasing function of the observation length. Actually, we see that doubling N results into a threshold decrease of approximately 3 dB with MLE, while a gain of 2 dB is observed with R&B.
As mentioned previously, in case of coherent detection phase recovery is required in addition to frequency and timing synchronization. The MSE of the phase estimates provided by MLE and R&B is illustrated in Figure 9 as a function of E s /N 0 . These results are qualitatively similar to those in Figure 8. In particular, it turns out that both schemes approach the CRB at intermediate/high SNR values, but MLE exhibits a lower threshold than R&B.

CONCLUSIONS
We have addressed the joint ML estimation of the carrier frequency, timing error, and carrier phase in burst-mode satellite transmissions. Thanks to a suitably designed pilot pattern composed of alternating binary symbols (which produces two spectral lines at ±1/2T), the estimation process can be divided into three separate steps, each devoted to the recovery of a single synchronization parameter. In particular, timing and phase recovery is accomplished in closed form, whereas the measurement of the frequency offset involves a grid search which represents the time-consuming part of the overall synchronization procedure.
Comparisons have been made with a conventional scheme in which timing recovery is accomplished in an NDA fashion and carrier synchronization is next achieved by exploiting the time-synchronized samples. Computer simulations indicate that the proposed ML algorithm provides more accurate timing estimates at low SNR values. In addition, it exhibits increased robustness against the occurrence of outliers in the frequency estimates. It is fair to say that these advantages are achieved at the price of a certain increase of the processing load as compared to the conventional scheme.
The question of which method is better is not easily answered because it depends on the different weights that may be given to the various performance indicators, including estimation accuracy, computational complexity, and constraints on the pilot pattern. It is likely that the choice will depend on the specific application. For example, the proposed algorithm seems attractive for coded transmissions as it approaches the relevant CRBs down to very low SNRs. On