Algorithm for frequency capture and rectification in a low-orbit satellite IoT communication network

1.   Introduction

Terrestrial cellular communication systems have developed rapidly in recent years, and the growing demand for mobile communications by users. In order to comply with this demand ^[1,2], the improvement of terrestrial cellular communication systems focuses on the quality of service and 5G research of existing standards. However, in some special periods or locations, it is difficult for terrestrial base stations to guarantee their communication services. Satellite Internet of Things (IoT) can be used as a supplement to ground communication and the means of emergency communication. ^[3,4].

Because the relative moving speed between the ground mobile terminal and the LEO satellite is very large, it generates very large Doppler frequency shift, and seriously affects the performance of the LEO satellite communication system ^[5]. Therefore, an important problem to be solved in low-orbit satellite mobile communication systems is the estimation and elimination of Doppler frequency offset. The frequency offset estimation algorithm mainly has a data assist algorithm and a data-free assistant algorithm. The difference between the two is whether to add known auxiliary data in the burst. Most data-assisted frequency offset algorithms are based on maximum likelihood estimation, among which are classical algorithms such as Fitz algorithm ^[6], Kay algorithm ^[7], Luise M and Reggianini R (L & R) algorithm ^[8] and Mengali U and Morelli M (M & M) algorithm ^[9]. The non-data auxiliary frequency offset algorithm is mainly based on the improvement of the data-assisted algorithm ^[10]. The article ^[11] improved the M & M algorithm, and proposed the A & C (A. B. Awoseyila and C. Kasparis) algorithm. After using the A & C algorithm for rough estimation, the paper uses the L & R algorithm for fine estimation to achieve a large range of frequency offset estimation at low signal-to-noise ratio (SNR). But the time complexity is high for the new proposed two-step estimation method. Figure 1 shows the flow chart of the receive signal processing for the satellite communication. In the literature ^[12], some detection and practical algorithms are given, which can be used to adaptively deal with different types of wireless fading channels.

The different frequency offset correction algorithms mentioned above have faced two problems in practical applications. One is how to adaptively process different satellite wireless fading channels, and the other is how to use parallel algorithms with low complexity to process multi-users at the same time. Aiming at solve above problems and the inability of current data-free offset estimation algorithm, when balancing the estimation range and the accuracy at low SNR, we propose a new frequency offset algorithm in this paper.

The main innovations and contributions of this article are: 1) The delay difference is used to obtain the phase increment information. After smoothing and linear weighting, the frequency offset is obtained. 2) The algorithm realizes the fast capture of the code phase by the FFT + IFFT method, and uses the serial search method to realize the Doppler frequency search. 3) In the frequency offset tracking and compensation phase, the forward-based open-loop algorithm is used for a rough estimation of frequency, and the forward-carrier phase-offset estimation algorithm is used to correct the remaining frequency offset and phase offset. The simulation results show that the proposed algorithm maintains high estimation accuracy under low SNR and has an extensive estimation range. It has been applied in the actual satellite IoT constellation.

The paper is organized as follows: first section is introduction The first section is an introduction; the second section is the mathematical theory of the Doppler frequency offset algorithm; the specific implementation methods of the application theory on the satellite communication system are described in the third and fourth sections respectively. They are frequency capture implementations and tracking scheme algorithm. In section five, we compare the proposed algorithm and other algorithms under different wireless scenes. At the last section, the conclusion is drawn from simulation results.

2.   Mathematical theoretical model

2.1.   Theoretical framework

Assuming that the received signal is accurately synchronized, for PSK modulated signals, the model can be expressed as:

where $i$ is the symbol position, $A(i)$ represents the magnitude of the complex number, $T_{c}$ is the symbol period, $\Delta f$ is the frequency offset, $\varphi_{m}$ is the modulation phase, $\varphi_{p}$ is the initial phase offset, and $N$ is the number of symbols. If $n(i)$ is a complex Gaussian white noise sequence, where both the real and imaginary parts have a mean value of 0 and the variance is $\sigma^{2}$. First, take the differential phase of $r(i)$:

where $r^{*}(i)$ represents the conjugate of $r(i), \arg []$ represents the phase operation, the amplitude range $(-\pi, \pi), \varphi_{n}(i)$ represents the phase noise caused by $n(i), \Delta \varphi_{m}(i)$ represents the phase difference of the modulation of the two symbols. For the data-less auxiliary frequency offset estimation algorithm of PSK modulated signals, the influence of modulation information can be eliminated by using the $M$-th power non-linear transformation ^[13].

Define the target function:

Vector $\Delta \varphi = [\Delta \varphi(1), \Delta \varphi(2), \ldots, \Delta \varphi(n)]^{T}, e = [1, 1, \ldots, 1]^{T}, [\cdot]^{T}$ represents transpose operation, and $\boldsymbol{C}$ is $\Delta \varphi(i)$ covariance matrix. According to the Kay algorithm principle, when the objective function $\boldsymbol{J}$ takes the minimum value, using Least squares principle, an estimated value $\Delta f$ can be obtained, and the solution given by it is:

The Kay algorithm can estimate a wide range of frequency offsets when the signal-to-noise ratio is high. But when estimating the low-signal-to-noise signal frequency offset, the noise threshold is high, result the poor estimation accuracy. When estimating the frequency offset of an actual signal, the frequency offset value is often estimated in a signal with a low signal-to-noise ratio. In order to improve the estimation accuracy at low signal-to-noise ratios, existing estimation algorithms usually reduce the estimation range. Taking Fitz as an example, the equation for calculating the frequency offset of the Fitz algorithm without data assistance is given directly below:

Among them, the signal to be estimated is modulated by the power of $M. R(m)$ represents an autocorrelation function with a delay of m symbols. The expression is:

Fitz can obtain high estimation accuracy by increasing the delay value $L$ when the signal-to-noise ratio is low. At this time, the estimation range will decrease; decreasing the $L$ value will increase the estimation range, but it will reduce the estimation accuracy. When using the Fitz algorithm for frequency offset estimation, it is necessary to evaluate the estimation accuracy and the estimation range, and select an appropriate delay value $L$. While use a two-step method to roughly estimate the approximate frequency offset first, and then to make an accurate estimate after compensating the original signal. The following comparison test between the Fitz algorithm and the new algorithm will choose the appropriate delay value $L$ according to the needs of the comparison, to make full use of the performance of a single aspect of the Fitz algorithm.

For non-AWGN channels, the level pass rate or zero crossing rate (LCR) estimation algorithm is used to estimate the envelope level pass rate of the received signal to estimate the Doppler shift.

Set the envelope level crossing rate of the received signal $r(t)$ as $L_{R}$, which is the rate at which the envelope passes through the prescribed level $R$ with a positive or negative slope. The literature ^[14] gives the expression of the level pass rate, in the Rice channel, the expression is:

The expression under Rayleigh channel is

Among them, $K$ Rice factor, $I_{0}$ () is the zero-order modified Bessel function, $\rho = \frac{R}{R_{r m s}}$ and $R_{r m s}$ is the root mean square envelope level.

In fading channels, the autocorrelation function of Eq (6) can be replaced by the following equation.

2.2.   The novel algorithm

Considering that the system uses DPSK modulation, such as π/4-CQPSK, its signal characteristics are essentially two QPSK constellation diagrams with phase difference π/4 alternating modulation. Then, when the conventional M-th power modulation method is directly used, there requires $M = 8$. Compared with the QPSK modulation frequency offset estimation range will be half smaller, and the non-linear changes will amplify the effect of noise ^[15], and improve the noise threshold of the algorithm. In addition, in the auto-correlation operation, the delayed odd number of symbols and those of even number will have a phase difference of π/4. Hence, the frequency offset algorithm that directly uses the high-delay autocorrelation operation is used for π/4-CQPSK is more difficult when modulating a signal.

The novel algorithm uses the signal model of Eq (1). After the symbol conjugate complex multiplication, we get:

where $\Delta \varphi_{m}$ is the phase difference between the modulation of adjacent symbols. In π/4-CQPSK modulation mode, its value is $\pm \frac{\pi }{4}, \pm \frac{{3\pi }}{4}$.

We take the phase of Eq (11) and have:

Regardless of the influence of the noise $\varphi_{n}(i)$, and assuming a phase deviation $2 \pi \Delta f T_{c} < \pm \frac{\pi}{4}$ caused by the frequency deviation $\Delta f$, the actual modulation phase difference $\Delta \varphi_{m}^{\prime}(i)$ can be obtained by directly judging the phase difference $\Delta \varphi(i)$, the judgment is as follows:

Considering only the modulation phase of the symbol $r(i)$, it can be considered as the sum of the modulation phase of the symbol $r(i-1)$ and the modulation phase difference of the two symbols. It can be known by recursion that the modulation phase of the current symbol can be regard ed as the sum of the modulation phase difference of the first symbol and the modulation phase differences of all preceding and following symbols before the current symbol. Then, the modulation phase of each symbol can be removed by removing the accumulated modulation phase:

here, $\Delta \varphi_{m}^{\prime}(0)$ is obtained by taking the phase of symbol $r(1)$ and judging by Eq (13). With the above method, the modulation phase can be removed without increasing the influence of noise. Here we rewrite the phase of each symbol after removing the modulation phase:

From Eq (15), it can be seen that the phase $\varphi^{\prime}(i)$ includes the frequency offset phase, and at the same time, the phase offset and noise phase remain. The frequency offset phase has a linear relationship with the symbol position $i$. In order to remove the effect of phase deviation while smoothing the noise, the phase increment $R(m)$ is defined as:

where

The phase in $R(m)$ is composed of frequency offset phase and noise. Among them, the frequency offset phase has a linear relationship with the delay value $m$; since LEO satellite is a carrier system, its fading is flat. Hence the noise term is an additive Gaussian white noise with a mean value of 0, and $R(m)$ is uniformly distributed for all $m$ values. It can be seen that at low $m$ values, the phase increasement $\varphi^{\prime}(i)-\varphi^{\prime}(i-m)$ the noise term in Eq (12) has a greater impact, but since $R(m)$ has more accumulated data, it can be smoothed after averaging the effect of noise terms. As the value of $m$ increases, the accumulated data decreases, but the multiplicative relationship between $\Delta f$ and $m$ in the high-delay $R(m)$ enhances the frequency offset phase term, giving it better noise immunity. At this time, the frequency offset estimation actually estimates $\Delta f$ from the phase increment $R(m)$ after mixed noise. Regardless the influence of noise, $R(m)$ and $m$ satisfy $R(m) = k m+b$, where $k = 2 \pi \Delta f T_{c}$. Therefore, the linear offset can be used to obtain the estimated frequency offset.

Before linear fitting, different weights of the phase increment $R(m)$ of different $m$ values need to be added to further reduce the phase error caused by the low $m$ value noise term and the accumulated data of high $m$ value are small and random the effect of largeness on estimation accuracy.

Let:

According to the linear fitting equation, the frequency offset value $\Delta f$ is:

3.   Frequency capture implementation

In order to achieve fast capture of the burst signal ^[16], the signal processing flow shown in Figure 1 is used to achieve fast demodulation of the spread burst signal.

● The digital down conversion module completes down-converting the intermediate frequency AD sampling data into a baseband signal;

● The downsampling module completes the data of the primary clock rate of the AD sampling data to a half chip rate to reduce the storage capacity requirement of the subsequent data buffering module;

● The data buffer module completes the buffering of the demodulated data during the working period of the capture module. The burst capture module needs to occupy a processing time for the spread spectrum signal capture ^[17], hence it is necessary to buffer a certain amount of half chip data to ensure that the subsequent modulated data will not be lost. The size of the data buffer RAM is related to the capture time required by the burst capture module.

The demodulation is responsible for implementing the subsequent DBPSK (Differentially coherent Binary PSK) differential demodulation.

3.1.   Method overview

Estimating the theoretical basis of the second section. For the satellite communication environment, based on the level pass rate estimation algorithm and the autocorrelation function estimation algorithm, an adaptive Doppler frequency shift estimation algorithm is proposed, which can be used in when the speed of the mobile station changes, the estimation length is adjusted adaptively to improve the estimation accuracy. We use MATLAB/Simulink as tools to simulate the efficiency of the proposed method and comparison it with other algorithms. The main parameters to estimate the effectiveness are relative speed, and the normalized mean square error.

The FFT+IFFT method is used to realize the fast capture of the code phase, and the serial search method is used to realize the Doppler frequency search ^[18]. We give mathematic derivation of this fast capture algorithm in the appendix.

The analysis and implementation of the capture method are as follows:

① The data is decelerated to 2 times the chip rate before the capture. Since the length of the spread code is 511, a total of 1024 chip intervals need to be searched. The Doppler search is implemented by means of FFT data serial shift. Since the FFT processing data rate is 2 times the chip rate, the data is analyzed by 2048 points FFT.

The System parameters for are listed in Table 1. The calculation according to the following equations:

1) Frequency resolution = chip rate/symbol interval * 2;

2) Doppler search times = maximum Doppler range / frequency resolution;

3) Each capture time = code phase search times × frequency search times/working clock;

4) Each fast capture requires 1 FFT calculation and multiple IFFT calculations, frequency search times = FFT times + IFFT times.

In this paper, Table 2 lists the sensitivity and Eb/N0 values of different symbol rates:

1) Working clock: 80 MHz (FPGA master clock);

2) Doppler search range: ± 4 KHz (after compensation);

② Non-coherent accumulation.

The probability of burst capture false alarm is $10^{-6}$ ^[19]. According to the theoretical curve of capture probability and false alarm probability (Figure 2), if it is necessary to achieve a detection probability of 0.95 or more, the detection signal to noise ratio needs to be higher than 14 dB.

Since the symbol signal-to-noise ratio is 9 dB, the detection signal-to-noise ratio requirement cannot be met, and the improvement of signal-to-noise ratio detection requires an inter-symbol cumulative averaging. There are two options for accumulating averages between symbols, one is coherent accumulation between symbols, and the other is non-coherent cumulative averaging between symbols. In order to increase the frequency scanning speed ^[20], a method of non-coherent cumulative averaging between symbols is employed. This design intends to improve the detection signal-to-noise ratio by using 16 inter-symbol accumulation averaging methods. The detection signal-to-noise ratio is 17.4 dB, which meet the detection probability requirement of 0.96.

3.2.   Capture process implementation block diagram and timing

Figures 3 and 4 show the capture process implementation block diagram and implementation timing.

The implementation steps of the capture process are as follows:

(ⅰ) Downsampling module

The down-converted data rate is 80 MHz and is downsampled to $2^{*} R_{c}$ at half chip rate.

(ⅱ) Cache module

The downsampled data needs to be buffered. The buffer size is 16 symbols long. According to the $2^{*} R_{c}$ rate, a total of 8 blocks of RAM are needed to store the down sampled data.

(iii) FFT-IFFT module

The buffered output data, such as the FFT-IFFT capture module, implements a 16-symbol non-coherent accumulation. The proposed working clock is 80 MHz, and 16 FFT calculations and 50/15 IFFT calculations are completed, occupying 21 ms/6.6 ms, corresponding to 7 symbols and 9 symbols, respectively ^[21]. The capture time is less than 16 symbols, meeting the time requirement for burst capture.

(iv) Capture judgment module

After the 16-symbol non-coherent accumulation of the FFT output value of each phase, the maximum value of the accumulated value is found to obtain the correct code phase and Doppler. Then the code phase and Doppler frequency difference acquired by the capture module are obtained. The information is forwarded to the subsequent demodulation module.

3.3.   Case analysis

3.3.1.   Code phase resolution

During the capture process, space where the signal exists is defined as the space of the frequency space and the code phase ^[22,23]. At each frequency point, the signal is searched according to a certain code phase interval. For the code phase interval, it is generally divided by $1 / 2 T_{c}$. Thus, if the code rate deviation is not considered, the code deviation is evenly distributed. The resulting signal loss is a maximum of 2.4 dB with an average of 1.13 dB.

3.3.2.   Doppler resolution

a) Carrier Doppler

Figure 5 shows the structure of the integration circuit. Let the receiving intermediate frequency signal be $c(t) e^{j\left(2 \pi f_{IF} t+\theta\right)}$, where $c(t)$ is the address code of the signal, and $f_{I F}$ is the signal containing the Doppler frequency offset. The IF frequency multiplies the local orthogonal reproduction carrier $e^{j(2 \pi _{IFL} ~~t)}$ and filters out the high-frequency components, then $r(t) = c(t) c_{L}(t) e^{-j 2 \pi f_{d} t+\theta}$.

Here $c_{L}(t)$ is the local recurring address code. If the above signals are integrated and accumulated according to a certain frequency (hereinafter referred to as the accumulated frequency) $f_{A}$, the integration time $\Delta t = \frac{1}{f_{A}}$, the output of the accumulator at the time $n \Delta t$ is

Set the chip duration $T_{C}$, integration time $\Delta T^{*} T_{C}$, Doppler frequency offset period $\frac{1}{f_{d}} * T_{C}$, then

Since $c_{L}(t) = c(t+\tau)$, then $A=\frac{R(\tau)}{T_{\text {epoch }}} T_{A}=\frac{T_{A}}{T_{\text {epoch }}} \quad \frac{T_{c}-\tau}{T_{c}}, \tau$ is the length of time between the reproduction code and the received signal address code. Further we can derive

Let the sampling rate be $f_{s}$, the corresponding sampling interval is $T_{s}, \Delta t = N T s$, then $f_{A} = \frac{1}{\Delta t} = \frac{1}{N} f_{s}, N = \frac{f_{s}}{f_{A}}$.

Since $T_{s} * \Delta t$, the integral realized by the digital accumulation can be expressed by the following equation:

Then the signal power is $C(n) = |r(n)|^{2} = A^{2} \frac{f_{s}^{2}}{f_{A}^{2}}, \operatorname{sinc}^{2}\left(\frac{f_{d}}{f_{A}}\right)$.

For narrowband Gaussian noise with a bandwidth of $B_{n}$ and a standard deviation of $\sigma$, the accumulator adds $\frac{f_{s}}{f_{A}}$ noise samples, so the power is $\sigma^{2}$ after accumulator $\frac{f_{s}}{f_{A}}$, so the SNR after the accumulator is

where $f_{s}$ is the down sampled sampling frequency, the value should be $B_{n}$ (i.e., 2 times the bandwidth of the baseband signal, at which time the noise in the sample is irrelevant and thus independent) ^[24], so the above equation is

It is observed that when $f_{d} = 0$ the integration time increases (with decreasing $f_{A}$), the SNR increases. While $f_{A}$ is constant, the SNR of $f_{d} \neq 0$ is less than the SNR of $f_{d} = 0$, the decreased part of SNR is the signal loss introduced by the frequency offset. Figure 6 shows the relationship between normalized SNR (i.e., set $\left.\frac{A^{2}}{\sigma^{2}} B_{n} T = 1\right)$ and ratio. According to Figure 6, the SNR loss introduced by this process increases as the $\frac{f_{d}}{f_{A}}$ increases.

This scheme is designed for $f_{d}$ is 167.6/665 Hz, and the SNR loss introduced is less than 1 dB.

b) Symbol Doppler

Figure 7 shows the case with large frequency difference. Here, only the case where the signal chip slides forward relative to the local reproduction code is considered ^[25], and the backward sliding situation is similar. When the capture is started, the corresponding chip offset of the signal chip and the reproduced code is $x_{0}$, and by the actual situation, $x_{0} \in\left[0, \frac{1}{4}\right]$.

During the entire capture process $\Delta t$ time, the signal chip moves forward by 1/2 chip, that is, the signal chip advance rate is $\frac{1 / 2}{\Delta t}$, so that the deviation between the signal and the reproduction code at time $t$ is $x(t) = x_{0}-\frac{1}{2 \Delta t} t.$ The correlation operation (that is, the coincidence of the partial product) has an amplitude of $1-|x(t)|$ and its energy is $(1-|x(t)|)^{2}$. Thus, for the case where the starting position is $x_{0}$, the average power during the whole process is

According to Figure 8, in the case of such a large frequency difference, the maximum loss is about 1.25 dB in the case where the maximum slip of the chip does not exceed 1/2 chip. Then for the capture process, this effect is equivalent to the average loss in the case of small frequency differences.

The maximum Doppler of this scheme is 167.6/665 Hz, and the integration time is 16 symbols long. The maximum code Doppler is 0.5 Hz, and the maximum slip of the chip is 0.008, which can neglect the loss.

c) Capture judgment

The energy of 1024 code branches is compared to find the maximum value of the correlation result. When the useful signal exists, the correct code phase and FFT spectrum can be judged. When there is no signal, the noise plus interference spectrum is sent.

After 16 symbols of non-coherent accumulation, each code branch outputs a 1024-point FFT result. Take out the largest of the line positions and analyze the Doppler.

Simulation conditions: ${E_b}/{N_0} = 7.2dB$, code phase offset 1/4 chip, Doppler take (−4 kHz, 4 kHz) point, 1000 Monte Carlo simulation. When the signal is present, the detection probability ${P_d}$ reaches 0.96.

4.   Tracking scheme algorithm

In actual communication, the instability of the oscillators of both transceivers, the Doppler shift during signal transmission, the mismatch of the transmission and reception filters, etc., cause the carrier frequency of the signal to be inconsistent with the carrier frequency of the receiver ^[26], so that the baseband is received. The signal produces a frequency offset that seriously affects the demodulation of the signal and must be eliminated. Carrier synchronization technology mainly includes carrier frequency synchronization and phase synchronization. Whether it is frequency synchronization or phase synchronization, its main algorithm can be classified into two types: closed-loop algorithm based on phase-locked loop and forward-opening algorithm. Both types of algorithms have their own advantages and disadvantages, as shown in Table 3. The closed-loop algorithm has high precision and simple calculation, but its acquisition time is long, which is not suitable for short-term burst signals. The forward algorithm has high complexity and relatively low precision, but its acquisition time is fast, suitable for short-term burst signals, carrier synchronization.

Because of the sampling feedback mechanism, the closed-loop algorithm has a convergence process. The convergence speed is determined by factors such as initial deviation and loop bandwidth. The open-loop algorithm directly estimates the frequency offset or phase offset from the received positive electricity sequence ^[27], and there is no convergence process. In theory, dozens or even a few symbols can get an estimate, but the more the number of samples, the higher the estimation accuracy.

For the blind reception of short-term burst signals, the carrier synchronization should first solve the fast capture problem. If the frequency offset and phase offset are not captured in a short period of time, other subsequent processing cannot be discussed at all ^[28]. Secondly, the case of a large frequency offset should be considered. Because it is blind reception, the carrier is usually estimated, allowing for large deviations. Third, on the basis of solving the first two problems, we should try to improve the synchronization capture accuracy, which have a significant impact on the final bit error rate. Fourth, the synchronization capture problem with low SNR should also be considered. Finally, the limits of payload on satellite trigger from the engineering implementation to minimize the computational complexity of the algorithm.

Based on the above analysis, the design uses the forward-based open-loop algorithm for a rough estimation of frequency, and the forward-carrier phase-offset estimation algorithm is used to correct the residual frequency offset and phase offset ^[29]. The forward open loop algorithm is used to a roughly estimate the frequency, which has quickly captured speed, and low tracking accuracy. At low SNR, the residual frequency offset is related to the symbol rate, hence the frequency offset is small, and the error is relatively large. If the phase change caused by the residual frequency offset is slow in a burst time of the burst signal, and is small enough. Then the influence of the residual frequency offset on the system performance can be reduced according to the carrier phase offset estimation algorithm ^[30].

4.1.   Improved L & R algorithm

The baseband signal received by the receiver that has obtained accurate symbol synchronization and has no crosstalk between codes can be expressed as

In the equation $a_{k}$, it is an independent and identically distributed isotropic data, $T$ is a symbol period, $f_{e}$ is an unknown carrier frequency offset, and $\theta_{0}$ is an unknown carrier phase, equals to the addition of ${\varphi _m}$ and ${\varphi _p}$ in Eq (1). For noise $n_{k}$.

$\left|A_{k}\right|, \eta_{k}$ are the magnitude and phase that related to $a_{k}.$ Here $\eta_{k}$ corresponds to $\Delta \varphi_{m}$ and $\varphi_{p}$ in Eq (7), and $\theta_{k}$ should correspond to $\Delta n$.

The phase $R(k)$ is set to:

wherein, the phase noise $v(k)$ corresponding to $\varphi_{n}$ in the autocorrelation function $R(k)$ of the demodulated signal is expressed, $M$ is MPSK, and "*" indicates that the conjugate operation is taken.

According to the above equation, the estimated carrier frequency offset is:

According to the above equation, the estimated carrier frequency offset is:

In the equation, the modulation mode $M$ is represented. For BPSK, $M = 2$, for QPSK, $M = 4$, $T$ the symbol period, $N$ is the symbol length, and when $L = N / 2$, the variance of the algorithm is the smallest. The larger the value $L$, the higher the estimation accuracy; the smaller $L$, the smaller the estimation range. In this design, $L = 1$, the frequency is estimated and the estimated range is $\left|\hat{f}_{e}\right| < \frac{1}{4 T}$.

The estimated range of this algorithm is:

a) Using the autocorrelation function of a short delay, we can get:

For the MPSK signal, $\eta_{k} = 0$

According to the above equation, the estimated carrier frequency offset is:

Since the range of $\arg [R(k)]$ is $[-\pi, +\pi]$, the noise term is ignored, and the frequency offset capture range for short delay is

b) Using the long-delay $L = l$ autocorrelation function for frequency offset estimation, Eqs (14) and (15) indicates:

For MPSK signal $[\eta(L)]_{-\pi}^{\pi} = 0$, then

If we ignore the influence of noise, the estimated value of the frequency offset under long delay can be expressed as:

When $L$ is $N / 2$, the estimation error is the smallest. Since the range of $\arg [R(L)]$ is $[-\pi, \pi]$, the noise influence is ignored, and the long-term delay frequency offset capture range is

c) From a) and b), It shows that the long-delay autocorrelation function has good anti-noise performance. But it is easy to produce phase folding, hence decrease the capture range. Estimation of the frequency offset by using the short-delay autocorrelation function can expand the capture range, but it is susceptible to noise. The estimation accuracy is not precise.

Therefore, firstly, the short-delay autocorrelation function $R(1)$ is used to estimate the approximate range of the frequency offset to obtain $\left|\hat{f}_{e 1}\right|$, and then the long-delay autocorrelation function $R(L)$ is used to estimate the exact value of the frequency offset, $\left|\hat{f}_{e 2}\right|$.

where $\lfloor x\rfloor$ represents the integer closest to $x$.

d) The advantage of the algorithm is that compared with the traditional algorithm, the algorithm takes into account the accuracy and estimation range of the frequency difference estimation, and the calculation amount is relatively small.

4.2.   Maximum likelihood estimation of carrier phase bias

After the carrier frequency offset is obtained by the frequency offset estimation algorithm, the frequency offset can be corrected as shown in Figure 9.

Due to the influence of noise and self-noise, the result obtained by the forward carrier frequency offset estimation algorithm described above is not completely equal to the actual frequency offset. Therefore, the received signal $x_{k}$ still undergoes the frequency offset correction and the signal y still contains the residual frequency offset. Because the residual frequency offset at this time is sufficiently small relative to the symbol rate, the phase change caused by the residual frequency offset in $y_{k}$ is slow during one burst of the burst signal. In one burst time, the carrier phase can be expressed as:

where $\hat{\varphi}_{i}$ denotes the average phase of the carrier in one burst time and $\eta(k)$ is noise. That is, there is no frequency deviation in a sudden period only the phase deviation indicated by $\hat{\varphi}_{i}$. In order to reduce the impact of residual frequency offset on system performance, it is necessary to estimate the $\hat{\varphi}_{i}$ information in the $y_{k}$, which is the carrier phase offset estimation problem.

If the burst time is short, it is enough to perform a phase estimation in one burst time. If the burst time is long, it is not reasonable to assume that there is no frequency offset within one burst time, for the burst time. In the longer case, the entire burst can be divided into several time periods, and each segment is assumed to have no frequency offset, and phase offset estimation is performed for each piece of data. The disadvantage of this processing method is that there is a phase jump between the segments, that is, the phase blur between the segments, and each phase estimation value $\hat{\varphi}_{i}$ needs to be de-jumped. Phase estimate of the de-jump process $\hat{\varphi}_{p}(i) = \hat{\varphi}_{i}+\xi_{i}$.

Introducing a counting period, assuming that the phase change caused by the residual frequency offset does not cause a phase jump during the counting period, hence

$L$ is the length of each phase estimate in a burst of data, $T$ is the symbol period, and $R$ is the number of segments in which the phase does not jump. When $phasChag = 1$, it indicates that a phase jump has occurred in the counting period $R$, and $phasChag = 0$ indicates that no phase jump has occurred in the counting period $R$. According to the definition of $R$, there can be at most one phase jump in each $R$ segment counting period.

5.   Simulation

In order to verify the accuracy of the design, we simulated the algorithm under the condition of $E_{b} / N_{0}=10 d B$, the symbol rate ${R_b} = 4K$, the carrier frequency offset $carr\_err=-0.23 R_{b}=-920 {~Hz}$, and the initial phase difference $=\pi / 9$. The burst signal length is 512 symbols.

For low-orbit satellites with a flying height of 590 km, with 45 degrees orbital inclination, the relative Doppler frequency offset and frequency offset rate at different times are calculated in Figures 12 and 13 under different orbital inclinations ^[31,32].

Figure 14 shows the burst signal with real and imaginary parts.

Figure 15 shows the data after frequency compensation, the blue line is the raw modulation data, and the red line is the frequency compensated data. The frequency estimation result is −916.2223 Hz, the frequency rough estimate is −919.3417 Hz, the fine estimate is 1.7465 Hz, and the adjusted frequency offset is-916.2223 Hz.

Figure 16 shows the demodulated output after carrier phase compensation, the blue line is the raw modulation data, and the red line is the phase compensated data. Phase estimation is performed every 16 symbols, and $R$ is 4, which means that there is no phase jump in $16 \times 4$ symbols. In theory, the maximum residual frequency difference $R_{b} /(64 \times 8)=8$ Hz can be tracked, and the actual simulation results show that 16 Hz residual frequency difference can be tracked. The number of symbols per phase estimate is reduced, or the value of $R$ can be increased to estimate the range of residual frequency differences ^[33,34].

It is evident that the demodulated output result is in good agreement with the original modulated data.

The simulation data symbol rate ${R_b} = 4K$, the filter roll-off coefficient is 0.25, and the optimal sampling one-way data is taken after synchronization. The simulation data takes the length $N{\rm{ = }}512$, and the SNR is 0 dB. Figure 17 compares the estimated range of different $L$-value Fitz algorithms with the new algorithm. Since the different $L$ values of the Fitz algorithm affect the estimation range of the algorithm, in order to represent the estimation range of the Fitz algorithm as much as possible, the Fitz algorithm with $L = 1$ and $L = 2$ are used for simulation after $M$-th demodulation. According to the simulation results, the estimation range of the new algorithm is larger than the Fitz algorithm of the $M$-th power demodulation.

Figure 18 is a graph under situation of normalized frequency offset $fT$ = 0.4, observed data length $N$ = 512, by comparison of different algorithms such as: classical L & R algorithm, Key algorithm, estimated variance of Fitz algorithm with different $L$ as a function of ${E_b}/{N_0}$ or SNR. The change in accuracy is estimated in the case of different SNR. Figure 18 shows that the estimation accuracy of the classical L & R algorithm is the highest in the whole SNR range, but it is at the cost of higher computational complexity and narrower estimation range, and poor performance at large frequency offset ^[35].

Compared with the Fitz algorithm and the Key algorithm, the new algorithm has noticeable improvements in estimation accuracy and without the signal-to-noise ratio threshold. The new algorithm uses FFT+IFFT to reduce the complexity of the algorithm. In the proposed method, the traditional L & R algorithm is divided into rough estimation and fine estimation, avoiding the disadvantage of the traditional method in the accuracy of large frequency offset ^[31], and can be adapted under different frequency offsets. The result shows that its performance can reach the Cramér–Rao bound.

6.   Conclusions

In this paper, the frequency offset estimation problem of satellite IoT communication system is studied. Aiming at the characteristics of DBPSK modulated signal, an algorithm based on FFT and IFFT for fast phase acquisition and Doppler frequency search is proposed, which preserves the frequency offset caused by each symbol. The relationship between phase and noise phase. Unlike the traditional $M$-th power demodulation method, this method does not introduce nonlinear changes. Thus, it does not amplify the effects of noise. The phase increment information is obtained by different delays. According to the relationship between the frequency offset phase and the noise phase, smoothing processing and function weighting are adopted to suppress the influence of noise. Finally, linear fitting is performed to obtain the frequency offset estimation value. The new algorithm has high estimation accuracy and an extensive estimation range. After comparison with simulation experiments, the new algorithm is superior to the $M$-th demodulation Fitz algorithm in estimating accuracy and estimation range, and can complete the frequency offset estimation of the satellite IoT system ^[36,37]. This algorithm is used for burst signal demodulation and has been implemented in hardware. The input data is bit-synchronized data, and the bit synchronization signal is stored in the RAM, and the output is demodulated data. First, the data in the RAM is roughly estimated by frequency. Subsequently, the frequency difference is used to estimate the frequency difference. Finally, the phase estimation and phase hopping are performed to obtain the demodulated data. The actual operation shows that the hardware implementation is consistent with the simulation results and can support the frequency offset correction of the satellite network.

Acknowledgments

This article was funded by the state grid's project "Key technology of scale engineering application of power battery for echelon utilization", Project No. 52010119002F.

Conflict of interest

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.

Authors contributions

Hanming Liu and Jie Shen contributed equally to this work, and should be considered as co-first authors.

Appendix

A.1. Remarks about the capture algorithm

The capture algorithm in this paper uses parallel code phases, and the search for all code phases can be completed at one time by Fast Fourier Transform (FFT), thereby improving the search efficiency. Figure A1 shows the circuit block diagram of fast capture algorithm.

After the FFT of the received signal is performed, the result of the transform is multiplied by the Fourier transform conjugate of the locally pseudo code. Inverse Fast Fourier Transform (IFFT) is performed on this value, then we obtain correlation values at each code phase on the current search band at once.

The parallel code phase search and capture algorithm uses FFT digital signal processing to replace digital correlators for digital operations. We set the correlation $z(n)$ between two periodic sequences $x(n)$ and $y(n)$ with length $N$ points as:

The Discrete Fourier Transform (DFT) is performed on the above correlation values $z(n)$, and the resulting DFT $Z(k)$ is as follows:

In the above equation, $X(k)$ and $Y(k)$ are DFT of $x(n)$ and $y(n)$ respectively, $\overline{Y(k)}$ representing the conjugate of the complex number $Y(k)$. The above equation shows that two sequences $x(n)$ and $y(n)$ perform correlation operations in the time domain, which is equivalent to the DFT of $X(k)$ and $\overline{Y(k)}$ product operations in the frequency domain. Conversely, the IDFT of the product $X(k) \overline{Y(k)}$ is exactly the correlation value $z(n)$ at each code phase that the receiver needs to detect. Once the receiver calculates the correlation value through the IFFT, then next signal detection is the same as the linear search capture method, that is, find the $|z(n)|$ of the autocorrelation amplitude in all search units, and compare the peak with the capture threshold. If the peak value exceeds the threshold, the receiver captures the signal and also obtains the frequency and code phase parameters from it.