# Two-Step Frequency Estimation of GNSS Signal in High-Dynamic Environment

Authors: Chao Wu; Yafeng Li; Liyan Luo; Jian Xie; Ling Wang
Journal: International Journal of Antennas and Propagation
 (2022)

---

## Abstract

To improve frequency accuracy which is affected by two parameters in high-dynamic acquisition, we propose a two-step frequency estimation method based on the mean frequency (MF) model for high-dynamic parameters estimation. The first step is based on the discrete chirp-Fourier transform (DCFT) for coarse MF estimation, where the MF accuracy and frequency search step are derived. In the second step, the maximum likelihood estimation process (MLEP) is adopted for fine MF estimation. Compared with state-of-art methods, it is verified that the two-step method can improve the detection probability in coarse MF estimation and improve the MF accuracy with low computational burdens under conditions with a moderate signal-to-noise ratio (SNR).


---

## Body

## 1. Introduction

In global navigation satellite system (GNSS) receiver techniques, the acquisition is the most important process for estimating code phase and carrier frequency [1]. For fast acquisition with lower computational complexity, the method [2–4] based on fast Fourier transform (FFT) was proposed. Since the pull-in range of the tracking loop is only a few hertz, the number of FFT points should be increased [5]. Generally, the coarse-to-fine acquisition methods are used to reduce the computation costs [6]. The code phase and coarse carrier frequency parameters can be obtained from the coarse acquisition, and the carrier frequency can be refined in a specific fine acquisition process with the code stripped off.For low-dynamic acquisition, the carrier Doppler can be estimated in the fine acquisition process. Tang et al. [7] proposed an accurate estimation method for residual Doppler. However, this method has a restriction on the initial Doppler search step, and more computation is required to obtain an accurate Doppler. To reduce computational load, a method [8] was proposed based on the coarse Doppler and sampling frequency in moderate SNR. However, to improve the Doppler accuracy in low SNR, a long-time correlation process is typically needed, which costs a lot of computations. To reduce the computations with long integration, Mohamed and Aboelmagd [5] proposed the Schmidt method which utilizes orthogonal searching. To further reduce the computations, article [9] proposed the zero-forcing and a double FFT-based method to improve Doppler frequency accuracy without increasing the computational load. However, because of the trade-off between the Doppler frequency resolution and the computational complexity, the maximum error of carrier frequency estimation depends on the number of FFT points. To improve the Doppler frequency accuracy, Nguyen et al. [10] proposed a residual frequency estimation method with differential processing. Due to the differential processing, it performs not well in the low SNR.Above all, the articles listed only focus on Doppler frequency accuracy. However, both initial frequency and chirping rate [11] affect the correlation peak in high-dynamic applications. Moreover, with a long integration time, the influence of these two parameters cannot be ignored. Among methods for initial frequency and chirping rate estimation, the authors in [12, 13] proposed frequency estimation methods based on Fractional Fourier transform (FRFT) for high-dynamic applications. In addition, we proposed a frequency estimation method based on discrete chirp-Fourier transform (DCFT) [14]. However, in some high-dynamic applications, more accurate frequency is usually desired.To further improve frequency accuracy for high-dynamic applications, a two-step frequency parameter estimation method is proposed in this paper. An MF model has been derived to improve the frequency accuracy. In the first step for coarse MF estimation, the chirping rate and initial frequency for MF estimation have been estimated based on the DCFT. A maximum likelihood estimation process (MLEP) has been proposed for the fine MF estimation in the 2nd step. With the two-step processing, the computational burdens can be reduced when the peak value is smaller than the configured threshold, and high-frequency accuracy can be obtained when the signal is present. Simulation results show that for coarse MF estimation, the proposed method has a higher detection probability compared with conventional methods, and for fine MF estimation, the two-step method has a higher frequency accuracy and lower computational burdens than the compared methods.


## 2. Signal Model

After correlating with a one-period local code and a coarse Doppler bin [11], the postcorrelation signal can be obtained and depicted as follows:(1)Sn=Abnexpj2πf0nTs+μn2Ts2+Wn,where f0 represents residual Doppler frequency or initial frequency, μ represents the chirping rate or Doppler rate, Ts represents the sampling frequency, bn represents bit sign, A represents signal amplitude, and Wn denotes a zero-mean additive white Gaussian noise (AWGN) process. When the received signal is not aligned with local code or a wrong Doppler bin is detected, signal amplitude A≈0, which is called the signal-absent situation in the following analysis. Or, A≠0 is assumed as a constant 11. For GPS L1 CA signal with 1-ms code period, it is typical that f0  = (−250, 250) Hz, and μ  = (−500, 500) Hz/s. It is assumed that bn can be obtained by some auxiliary means 13, and bn equals 1 in the following analysis.


## 3. Proposed Method

In this section, the process based on DCFT has been proposed for coarse MF estimation. Then, MLEP has been adopted for fine MF estimation. Finally, the two-step method which combines the two processes has been proposed. In this two-step method, signal detection and MF accuracy improvement are realized through the first and second steps respectively.


### 3.1. Coarse Search of MF Based on DCFT

T transform of the postcorrelation signal can be written as follows:(2)STkT,αT=AT∑n=0N−1bnexpj2πΔTnTs+δTn2Ts2+w,where kT=0,±1,±2,... represents a searching range, and αT represents transform factor. w=wi+jwr. wi and wr both obey normal distribution N0,σw2. When T = f, it represents FRFT [15]. Af=Aexp−jπsgnsinαf/4+jαT/2/sinαf0.5expj1/2cotαfkF2. F=2π/NTscscαf. Δf=f0−kf/NTs. δf=μ+1/4πcotαf. When T = d, it represents DCFT. Ad=A. Δd=f0−kd/NTs. δd=μ−αd/N2Ts2, where αd represents the chirping rate factor. Based on (2), the correlation peak of FRFT is an unilinear function with αf, which may degrade the detection peak. Therefor DCFT is chosen in the following discussion.When the bit signs can be obtained by assisted means, the formula above can be simplified into(3)Sdkd,αd≈Ad∑n=0N−1expjω¯nTs+w=AdNsincω¯NTs2expjω¯N−12Ts+w,where N represents the integration time and ω¯=2πΔT+2πδTN−1/2Ts represents the MF from 0Ts to N−1Ts. It is assumed that A0=AdN. Based on the derivations above and Taylor expansion, the peak Ad∑n=0N−1expjω¯nTs can be approximated as(4)Ad∑n=0N−1expjω¯nTs=A0sincω¯NTs2≈A01−16ω¯NTs22,where Ad∑n=0N−1expjω¯nTs represents the amplitude of the signal Ad∑n=0N−1expjω¯nTs. It is assumed that the unit kd0,αd0 is corresponding to the peak, and Ad∑n=0N−1expjω¯nTs≥γA0, where γ is set based on the criterion that one search bin contains mostly useful energy [16]. Then, we can obtain(5)ω¯≤261−γ/NTs.Based on the equations above, the search step of MF in the first step is set to261−γ/NTs. Due to the influence between the initial frequency and chirping rate, the search step of the initial frequency needs to be configured first. Based on the influence of residual Doppler frequency on the correlation peak value of the postcorrelation signal in low-dynamic applications [4], the search step of initial frequency is set to 1/2NTs. Moreover, based on the relationship between the initial frequency search step and the MF search step in the first step, the search step of the chirping rate can be obtained.Above all, in the coarse search, the coarse estimations of both initial frequency and chirping rate can be obtained. Based on the coarse estimation, the coarse MF estimation can be obtained.


### 3.2. MLEP for Fine MF Estimation

In the MLEP, firstly, the signal amplitude is estimated based on ML; then, due to the fine MF range, the restricted search criteria for the fine MF estimation are adopted. Finally, combining the criteria and ML function, the fine MF is estimated based on the estimated signal amplitude.Based on (3), the observed peak can be written as follows:(6)SdA0,ω¯0=A0sincω¯0NTs2expjω¯0N−12Ts+w,where ω¯0 represents residual MF, which ranges from −61−γ/NTs to 61−γ/NTs based on Section 3.1 analysis. The joint probability density function of wi,wr of wi,wr can be written as(7)fwi,wr=12πσw2exp−wi2+wr22σw2.Then, based on ML estimation, the optimized objective function can be obtained as(8)minJA0,ω¯0=min−lnfwi,wr≈minwi2+wr22σw2,where minJA0,ω¯0 represents the minimum of the objective function JA0,ω¯0. It is assumed that ∂JA0,ω¯0/∂A0=0. Then, we can obtain(9)A0c=ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tssincω¯0NTs/2cosω¯0N−1Ts,where A0c represents the optimized value of A0 based on ML. Then, substituting (9) into (6), we can obtain(10)SdA0c,ω¯0=ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Ts·expjω¯0N−12Ts+w,where SdA0c,ω¯0 represents the correlation peak in presence of noise. However, when cosω¯0N−1Ts=0, ω¯0=π2k+1/2N−1Ts represents singular points of (10), k=0,±1,±2,... represents segmentation variable, Lk represents the length of k. Consequently, the segmentation optimization is taken based on the singular points. When π2k+1/2N−1Ts<ω¯0<π2k+3/2N−1Ts, the segmentation optimization is based on the objective function as follows:(11)JA0c,ω¯0+Δω=wiω¯0+Δω2+wrω¯0+Δω22σw2=wiω¯0+JiΔω2+wrω¯0+JrΔω22σw2,where J represents the objective function. It is assumed that ∂JA0c,ω¯0+Δω/∂Δω=0, Δω can be obtained as follows:(12)Δω=−JiTwiω¯0−JrTwrω¯0JiTJi+JrTJr,where Ji=∂wiω¯0/∂ω¯0. Jr=∂wrω¯0/∂ω¯0. The several simulations show that the number of iterations It can be chosen to be 10. The local optimal solution ω¯0,k segmentation variable k can be obtained based on the range of ω¯0. The restricted search criterion for choosing ω¯0,k is given as follows:(13)k0=minkJA0c,ω¯0,k,sincω¯0,kNTs2>0.1,A0c|ω¯0,k>0,where ω¯0,k0 is the estimated frequency parameter in the 2nd step. Above all, MF can be obtained from the fine MF estimation of the 2nd step.


### 3.3. Two-Step Frequency Estimation Method

The two-step frequency parameters estimation method is shown in Figure1. The method can be depicted in more detail as follows:(a)
Calculating the search step of the initial frequency and chirping rate based on (5).(b)
In coarse search, the estimated valueα˜0,f˜0 can be obtained based on the threshold γ.(c)
In fine search, an iterative approach based on ML is adopted:With frequency errorΔω¯0 initialized, the amplitude can be obtained after calculating the amplitude function (9). Moreover, based on ΔT,δT, the range of ω¯0 can be obtained.Based on the amplitude and mean frequency error, the value of the differential functionsJi and Jr can be obtained as follows:(14)Ji=Ts∗TT1∗ImSd∗Tc+TT1∗ReSd∗TsT2c−TT1∗Tc∗ReSd∗Tc−ImSd∗TsT2c−2∗TT1∗Ts∗T2s∗ReSd∗Tc−ImSd∗Ts,T2c2Jr=Tc∗TT1∗ImSd∗Tc+TT1∗ReSd∗TsT2c+TT1∗Ts∗ReSd∗Tc−ImSd∗TsT2c−2∗TT1∗Tc∗T2s∗ReSd∗Tc−ImSd∗TsT2c2,whereTc  =  cosTT1, Ts  =  sinTT1, T2c  =  cos2TT1, T2s  =  sin2TT1, and TT1  =  N−1/2Tsω¯0. The peak error functions wiω¯0 and wrω¯0 can be written as follows:(15)wiω¯0=ImSdk0,α0−ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Tssinω¯0N−12Ts,wrω¯0=ReSdk0,α0−ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Tscosω¯0N−12Ts,where the unitk0,α0 is corresponding to α˜0,f˜0.Based on differential functions, peak error function and (11), feedback error Δω can be obtained. The feedback function in Figure 1 can be written as follows:(16)ω¯0=ω¯0−Δω.(d)
Based on step (c),ω¯0,k can be obtained after 10 iterations. Then, based on (13), the final MF can be estimated. Above all, the frequency accuracy can be improved.Figure 1
The two-step method’s flow diagram combining DCFT for coarse frequency estimation and MLEP for fine frequency estimation.


## 3.1. Coarse Search of MF Based on DCFT

T transform of the postcorrelation signal can be written as follows:(2)STkT,αT=AT∑n=0N−1bnexpj2πΔTnTs+δTn2Ts2+w,where kT=0,±1,±2,... represents a searching range, and αT represents transform factor. w=wi+jwr. wi and wr both obey normal distribution N0,σw2. When T = f, it represents FRFT [15]. Af=Aexp−jπsgnsinαf/4+jαT/2/sinαf0.5expj1/2cotαfkF2. F=2π/NTscscαf. Δf=f0−kf/NTs. δf=μ+1/4πcotαf. When T = d, it represents DCFT. Ad=A. Δd=f0−kd/NTs. δd=μ−αd/N2Ts2, where αd represents the chirping rate factor. Based on (2), the correlation peak of FRFT is an unilinear function with αf, which may degrade the detection peak. Therefor DCFT is chosen in the following discussion.When the bit signs can be obtained by assisted means, the formula above can be simplified into(3)Sdkd,αd≈Ad∑n=0N−1expjω¯nTs+w=AdNsincω¯NTs2expjω¯N−12Ts+w,where N represents the integration time and ω¯=2πΔT+2πδTN−1/2Ts represents the MF from 0Ts to N−1Ts. It is assumed that A0=AdN. Based on the derivations above and Taylor expansion, the peak Ad∑n=0N−1expjω¯nTs can be approximated as(4)Ad∑n=0N−1expjω¯nTs=A0sincω¯NTs2≈A01−16ω¯NTs22,where Ad∑n=0N−1expjω¯nTs represents the amplitude of the signal Ad∑n=0N−1expjω¯nTs. It is assumed that the unit kd0,αd0 is corresponding to the peak, and Ad∑n=0N−1expjω¯nTs≥γA0, where γ is set based on the criterion that one search bin contains mostly useful energy [16]. Then, we can obtain(5)ω¯≤261−γ/NTs.Based on the equations above, the search step of MF in the first step is set to261−γ/NTs. Due to the influence between the initial frequency and chirping rate, the search step of the initial frequency needs to be configured first. Based on the influence of residual Doppler frequency on the correlation peak value of the postcorrelation signal in low-dynamic applications [4], the search step of initial frequency is set to 1/2NTs. Moreover, based on the relationship between the initial frequency search step and the MF search step in the first step, the search step of the chirping rate can be obtained.Above all, in the coarse search, the coarse estimations of both initial frequency and chirping rate can be obtained. Based on the coarse estimation, the coarse MF estimation can be obtained.


## 3.2. MLEP for Fine MF Estimation

In the MLEP, firstly, the signal amplitude is estimated based on ML; then, due to the fine MF range, the restricted search criteria for the fine MF estimation are adopted. Finally, combining the criteria and ML function, the fine MF is estimated based on the estimated signal amplitude.Based on (3), the observed peak can be written as follows:(6)SdA0,ω¯0=A0sincω¯0NTs2expjω¯0N−12Ts+w,where ω¯0 represents residual MF, which ranges from −61−γ/NTs to 61−γ/NTs based on Section 3.1 analysis. The joint probability density function of wi,wr of wi,wr can be written as(7)fwi,wr=12πσw2exp−wi2+wr22σw2.Then, based on ML estimation, the optimized objective function can be obtained as(8)minJA0,ω¯0=min−lnfwi,wr≈minwi2+wr22σw2,where minJA0,ω¯0 represents the minimum of the objective function JA0,ω¯0. It is assumed that ∂JA0,ω¯0/∂A0=0. Then, we can obtain(9)A0c=ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tssincω¯0NTs/2cosω¯0N−1Ts,where A0c represents the optimized value of A0 based on ML. Then, substituting (9) into (6), we can obtain(10)SdA0c,ω¯0=ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Ts·expjω¯0N−12Ts+w,where SdA0c,ω¯0 represents the correlation peak in presence of noise. However, when cosω¯0N−1Ts=0, ω¯0=π2k+1/2N−1Ts represents singular points of (10), k=0,±1,±2,... represents segmentation variable, Lk represents the length of k. Consequently, the segmentation optimization is taken based on the singular points. When π2k+1/2N−1Ts<ω¯0<π2k+3/2N−1Ts, the segmentation optimization is based on the objective function as follows:(11)JA0c,ω¯0+Δω=wiω¯0+Δω2+wrω¯0+Δω22σw2=wiω¯0+JiΔω2+wrω¯0+JrΔω22σw2,where J represents the objective function. It is assumed that ∂JA0c,ω¯0+Δω/∂Δω=0, Δω can be obtained as follows:(12)Δω=−JiTwiω¯0−JrTwrω¯0JiTJi+JrTJr,where Ji=∂wiω¯0/∂ω¯0. Jr=∂wrω¯0/∂ω¯0. The several simulations show that the number of iterations It can be chosen to be 10. The local optimal solution ω¯0,k segmentation variable k can be obtained based on the range of ω¯0. The restricted search criterion for choosing ω¯0,k is given as follows:(13)k0=minkJA0c,ω¯0,k,sincω¯0,kNTs2>0.1,A0c|ω¯0,k>0,where ω¯0,k0 is the estimated frequency parameter in the 2nd step. Above all, MF can be obtained from the fine MF estimation of the 2nd step.


## 3.3. Two-Step Frequency Estimation Method

The two-step frequency parameters estimation method is shown in Figure1. The method can be depicted in more detail as follows:(a)
Calculating the search step of the initial frequency and chirping rate based on (5).(b)
In coarse search, the estimated valueα˜0,f˜0 can be obtained based on the threshold γ.(c)
In fine search, an iterative approach based on ML is adopted:With frequency errorΔω¯0 initialized, the amplitude can be obtained after calculating the amplitude function (9). Moreover, based on ΔT,δT, the range of ω¯0 can be obtained.Based on the amplitude and mean frequency error, the value of the differential functionsJi and Jr can be obtained as follows:(14)Ji=Ts∗TT1∗ImSd∗Tc+TT1∗ReSd∗TsT2c−TT1∗Tc∗ReSd∗Tc−ImSd∗TsT2c−2∗TT1∗Ts∗T2s∗ReSd∗Tc−ImSd∗Ts,T2c2Jr=Tc∗TT1∗ImSd∗Tc+TT1∗ReSd∗TsT2c+TT1∗Ts∗ReSd∗Tc−ImSd∗TsT2c−2∗TT1∗Tc∗T2s∗ReSd∗Tc−ImSd∗TsT2c2,whereTc  =  cosTT1, Ts  =  sinTT1, T2c  =  cos2TT1, T2s  =  sin2TT1, and TT1  =  N−1/2Tsω¯0. The peak error functions wiω¯0 and wrω¯0 can be written as follows:(15)wiω¯0=ImSdk0,α0−ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Tssinω¯0N−12Ts,wrω¯0=ReSdk0,α0−ReSdcosω¯0N−1/2Ts−ImSdsinω¯0N−1/2Tscosω¯0N−1Tscosω¯0N−12Ts,where the unitk0,α0 is corresponding to α˜0,f˜0.Based on differential functions, peak error function and (11), feedback error Δω can be obtained. The feedback function in Figure 1 can be written as follows:(16)ω¯0=ω¯0−Δω.(d)
Based on step (c),ω¯0,k can be obtained after 10 iterations. Then, based on (13), the final MF can be estimated. Above all, the frequency accuracy can be improved.Figure 1
The two-step method’s flow diagram combining DCFT for coarse frequency estimation and MLEP for fine frequency estimation.


## 4. Algorithm Performance

In this section, Cramér–Rao bound of the estimated MF of the proposed method is derived. Then, computational burdens and detection probability of the proposed method are analyzed for performance evaluation.


### 4.1. Cramér–Rao Bound (CRB) of the Proposed Method

Based on the theory [17], the CRB of ω¯0 can be written as follows:(17)CRω¯0=1E∂2wi2+wr2/2σw2/∂2ω¯0=2σw2E∂2wi2+wr2∂2ω¯0,where wi=ImSd−A0sincω¯0NTs/2sinω¯0N−1/2Ts, and wr=ReSd−A0sincω¯0NTs/2cosω¯0N−1/2Ts. When sincω¯0NTs/2≈1,(18)CRω¯0≈2σw2N−12Ts2A02,where we can obtain the final CRB of ω¯0. When A0≈N, Ts=0.001s, N  = 200 ms, 400 ms and 600 ms, the CRB of MF is shown in Figure 2.Figure 2
CRB ofω¯0 under different integration times.In Figure2, CRB curves are shown with the integration time being 100, 200, and 300 ms, respectively. Under the same SNR, the longer the integration time is, the higher the frequency estimation accuracy is. This is because based on (19) the integration time is long, and the integration peak value is large.


### 4.2. Computations Analysis

Based on (2) and (10), the computations of the proposed method can be obtained. Here, the BASIC [11] is chosen as the benchmark for the coarse MF estimation. The method estimates the frequency parameters based on the differential signal as follows:(19)Sdn=S∗nSn+M0=A2bnbn+M0expj2π2nμM0Ts2+f0M0Ts+μM02Ts2=A2bnbn+M0expj2π2nμM0Ts2+θ0,where θ0=f0M0Ts+μM02Ts2. Based on the Fourier transform, the frequency parameters can be estimated.In Table1, methods are chosen for coarse or fine MF estimation. Tα represents the number of chirping rate search bins and Tf represents the number of initial frequency search. CP,M and CP,A can be calculated as follows.Table 1
Computational burdens comparison.
Coarse MF estimationComplex multiplicationsComplex additionsBASICTα+Tα/2log2Tα+Tf+Tf/2log2TfTαlog2Tα+Tflog2TfProposed method (1st step)TαTf/2log2TfTαTflog2TfFine MF estimationComplex multiplicationsComplex additionsProposed method(2nd step)N+CP,MCP,ASchmidt method (appendix table2) [5]2N+1M+1+1+2M!+M−1!+N+1∑m=1Mm!/20.5N−1M+2+MN+1+M−1!+M−2!+N∑m=1Mm!Table 2
Computational burdens for the Schmidt method algorithm [5].
MultiplicationsAdditionsR1,1=P12NN−1Rk,1=PkP1¯k=1,...,M1+MNMN−1C1=yP1¯N+1N−1fork=1,…,Mfori=1,…,kk>iRk,i=PkPi¯−∑j=1i−1αjiRk.jk=iRk,k=PkPk¯−∑j=1k−1αjk2Rj.j∑m=1Mm!N+1∑m=1Mm!Nfork=1,…,MCk=yPk¯−∑j=1k−1αjkCjMN+M−1!MN+M−2!gk=Ck/Rk,kM0θi=gi−∑k=i+1mRk,i/Ri,iθki=1,…,M2M!M−1!+MM represents the number of vectors, and N represents the length of the vector. gk represents coefficients of Schmidt orthogonalization.Firstly, A is calculated from (9). Calculating (9) involves 12 multiplications and 1 addition. Moreover, the value of cos () or sin () function can be realized by a look-up table, and their computations can be ignored. Then, calculating differential function costs 18∗2 multiplications and 2∗2 additions. Calculating the peak error function costs 5∗2 multiplications and 1∗2 additions. Hereafter, (12) costs 5 multiplications and 2 additions. Calculating the feedback function costs 1 addition. In addition, calculating (13) costs 19 multiplications and 3 additions.Above all, one complex multiplication equals two multiplications [4]. So, CP,M=Lk51Ii+19/2 and CP,A=Lk9Ii+3/2, where It represents the number of iterations and Lk represents the number of segmentations. Besides, Table 1 shows that the frequency accuracy of the Schmidt method is also dependent on the number of vectors.Since the computations simulation needs to set lots of simulation parameters, the simulation will be conducted in Section5.


### 4.3. Detection Performance

Since the signal is detected in the first step of the proposed method, the section is to discuss the detection probability of the first step in the proposed method. The theoretical simulation will be conducted in Section5.1.The detection variableSdk0,α02 obeys the chi-square distribution. When the signal is absent or a wrong frequency bin is searched, Jd=Sdk0,α02 obeys the central chi-square distribution with the variance σw2, and the probability density function can be written as follows:(20)p0Jd=12σw2exp−Jd2σw2,where p0 represents probability density function when a wrong bin is searched. When the right frequency unit is detected, Jd=Sdk0,α02 obeys the noncentral chi-square distribution with the variance σw2, and the probability density function can be written as follows:(21)p1Jd=12σw2exp−Jd+a22σw2I0a2Jdσw4,where p1 represents the probability density function when the right frequency bin is detected. Under an incorrect frequency hypothesis H0, the false alarm probability Pfa can be written as:(22)Pfa=Py≥γD|H0=1−∫0γDp0ydyTfTα−1≈1−1−exp−γD2σw2TfTα−1,where γD represents the detection threshold. Under the correct frequency hypothesis H1, the detection probability PD can be written as follows:(23)PD=Px≥γD|H1=∫γD+∞p1x∫0xp0ydyTfTα−1dx,where Tf represents the number of initial frequency search bins in the first step of the proposed method and Ta represents the number of chirping rate search bins. When the configured Pfa is small, ∫0Jdp0JddJdTfTα−1≈1, and the detection probability PD can be simplified into(24)PD=PJd≥γD/H1≈Qaσw,γDσw,where a equals to Sdk0,α0 in the absence of noise. Based on the definition of miss detection probability [18], the miss probability of the proposed method can be written as:(25)PM=∫0γDp1x∫0γDp0ydyTfTα−1dx.When the setPfa is small, ∫0γDp0JddJdTfTα−1≈1, PD+PM=1. Above all, the detection probability of coarse MF estimation of the proposed method can be obtained. When the signal is detected, the second step of the proposed method for fine MF estimation can be performed.


## 4.1. Cramér–Rao Bound (CRB) of the Proposed Method

Based on the theory [17], the CRB of ω¯0 can be written as follows:(17)CRω¯0=1E∂2wi2+wr2/2σw2/∂2ω¯0=2σw2E∂2wi2+wr2∂2ω¯0,where wi=ImSd−A0sincω¯0NTs/2sinω¯0N−1/2Ts, and wr=ReSd−A0sincω¯0NTs/2cosω¯0N−1/2Ts. When sincω¯0NTs/2≈1,(18)CRω¯0≈2σw2N−12Ts2A02,where we can obtain the final CRB of ω¯0. When A0≈N, Ts=0.001s, N  = 200 ms, 400 ms and 600 ms, the CRB of MF is shown in Figure 2.Figure 2
CRB ofω¯0 under different integration times.In Figure2, CRB curves are shown with the integration time being 100, 200, and 300 ms, respectively. Under the same SNR, the longer the integration time is, the higher the frequency estimation accuracy is. This is because based on (19) the integration time is long, and the integration peak value is large.


## 4.2. Computations Analysis

Based on (2) and (10), the computations of the proposed method can be obtained. Here, the BASIC [11] is chosen as the benchmark for the coarse MF estimation. The method estimates the frequency parameters based on the differential signal as follows:(19)Sdn=S∗nSn+M0=A2bnbn+M0expj2π2nμM0Ts2+f0M0Ts+μM02Ts2=A2bnbn+M0expj2π2nμM0Ts2+θ0,where θ0=f0M0Ts+μM02Ts2. Based on the Fourier transform, the frequency parameters can be estimated.In Table1, methods are chosen for coarse or fine MF estimation. Tα represents the number of chirping rate search bins and Tf represents the number of initial frequency search. CP,M and CP,A can be calculated as follows.Table 1
Computational burdens comparison.
Coarse MF estimationComplex multiplicationsComplex additionsBASICTα+Tα/2log2Tα+Tf+Tf/2log2TfTαlog2Tα+Tflog2TfProposed method (1st step)TαTf/2log2TfTαTflog2TfFine MF estimationComplex multiplicationsComplex additionsProposed method(2nd step)N+CP,MCP,ASchmidt method (appendix table2) [5]2N+1M+1+1+2M!+M−1!+N+1∑m=1Mm!/20.5N−1M+2+MN+1+M−1!+M−2!+N∑m=1Mm!Table 2
Computational burdens for the Schmidt method algorithm [5].
MultiplicationsAdditionsR1,1=P12NN−1Rk,1=PkP1¯k=1,...,M1+MNMN−1C1=yP1¯N+1N−1fork=1,…,Mfori=1,…,kk>iRk,i=PkPi¯−∑j=1i−1αjiRk.jk=iRk,k=PkPk¯−∑j=1k−1αjk2Rj.j∑m=1Mm!N+1∑m=1Mm!Nfork=1,…,MCk=yPk¯−∑j=1k−1αjkCjMN+M−1!MN+M−2!gk=Ck/Rk,kM0θi=gi−∑k=i+1mRk,i/Ri,iθki=1,…,M2M!M−1!+MM represents the number of vectors, and N represents the length of the vector. gk represents coefficients of Schmidt orthogonalization.Firstly, A is calculated from (9). Calculating (9) involves 12 multiplications and 1 addition. Moreover, the value of cos () or sin () function can be realized by a look-up table, and their computations can be ignored. Then, calculating differential function costs 18∗2 multiplications and 2∗2 additions. Calculating the peak error function costs 5∗2 multiplications and 1∗2 additions. Hereafter, (12) costs 5 multiplications and 2 additions. Calculating the feedback function costs 1 addition. In addition, calculating (13) costs 19 multiplications and 3 additions.Above all, one complex multiplication equals two multiplications [4]. So, CP,M=Lk51Ii+19/2 and CP,A=Lk9Ii+3/2, where It represents the number of iterations and Lk represents the number of segmentations. Besides, Table 1 shows that the frequency accuracy of the Schmidt method is also dependent on the number of vectors.Since the computations simulation needs to set lots of simulation parameters, the simulation will be conducted in Section5.


## 4.3. Detection Performance

Since the signal is detected in the first step of the proposed method, the section is to discuss the detection probability of the first step in the proposed method. The theoretical simulation will be conducted in Section5.1.The detection variableSdk0,α02 obeys the chi-square distribution. When the signal is absent or a wrong frequency bin is searched, Jd=Sdk0,α02 obeys the central chi-square distribution with the variance σw2, and the probability density function can be written as follows:(20)p0Jd=12σw2exp−Jd2σw2,where p0 represents probability density function when a wrong bin is searched. When the right frequency unit is detected, Jd=Sdk0,α02 obeys the noncentral chi-square distribution with the variance σw2, and the probability density function can be written as follows:(21)p1Jd=12σw2exp−Jd+a22σw2I0a2Jdσw4,where p1 represents the probability density function when the right frequency bin is detected. Under an incorrect frequency hypothesis H0, the false alarm probability Pfa can be written as:(22)Pfa=Py≥γD|H0=1−∫0γDp0ydyTfTα−1≈1−1−exp−γD2σw2TfTα−1,where γD represents the detection threshold. Under the correct frequency hypothesis H1, the detection probability PD can be written as follows:(23)PD=Px≥γD|H1=∫γD+∞p1x∫0xp0ydyTfTα−1dx,where Tf represents the number of initial frequency search bins in the first step of the proposed method and Ta represents the number of chirping rate search bins. When the configured Pfa is small, ∫0Jdp0JddJdTfTα−1≈1, and the detection probability PD can be simplified into(24)PD=PJd≥γD/H1≈Qaσw,γDσw,where a equals to Sdk0,α0 in the absence of noise. Based on the definition of miss detection probability [18], the miss probability of the proposed method can be written as:(25)PM=∫0γDp1x∫0γDp0ydyTfTα−1dx.When the setPfa is small, ∫0γDp0JddJdTfTα−1≈1, PD+PM=1. Above all, the detection probability of coarse MF estimation of the proposed method can be obtained. When the signal is detected, the second step of the proposed method for fine MF estimation can be performed.


## 5. Simulation Results

In this section, BASIC [11] and FRFT [13] are chosen as the benchmark for the 1st step of the proposed method, and Schmidt [5] is chosen as the benchmark for the 2nd step. The simulation parameters are listed in Table 3 where ⌈η⌉ represents the smallest integer that is larger than η.Table 3
Simulation parameters.
ParametersValueInitial frequency range(−250, 250) HzChirping rate range(−500, 500) Hz/sInitial frequency search stepΔd1/2T0=1/2NTsSampling timeTs0.001 sFactorγ0.5Number of iterationsIt10Number of segmentationLk⌈4N−1Tsω¯max/π⌉Monte Carlo simulation5000Max mean frequencyω¯max2πΔd+2πδdN−1/2TsFalse alarm probabilityPfa2×10−10


### 5.1. Coarse Frequency Detection Performance Comparison

Although complex multiplications of the proposed method for coarse MF estimation are larger than that of BASIC in Figure3, when the SNR of the postcorrelation signal is larger than −10 dB, the detection probability of the proposed method is almost 100%, which is larger than other FRFT and BAISIC in Figure 4. This is because FRFT has a search bin α and BAISIC adopts a differential process, which may degrade the correlation peak and lead to lower detection probability.Figure 3
Complex multiplications comparison for coarse MF estimation methods.Figure 4
Detection probabilities comparison when is 200 ms, chirping rate is 200 Hz/s, and initial frequency is 100 Hz.


### 5.2. Fine MF Accuracy and Complexity Comparison

Based on the simulation above, the DCFT method, two-step method, and Schmidt method are adopted for the fine MF search. In Figure5, the complex multiplications of Schmidt vary greatly with the change of MF search step and postcorrelation signal length.Figure 5
Complex multiplications comparison for fine MF estimation methods.In Figure6, even though 1 rad/s of MF search step is adopted for MF estimation, the two-step method based on MLEP gains higher precision than MF search based on Schmidt.Figure 6
MF square error comparison. The SNR of the postcorrelation signal is 5 dB. The MF search step is 1 rad/s. The integration time is 200 ms in coarse MF estimation.


## 5.1. Coarse Frequency Detection Performance Comparison

Although complex multiplications of the proposed method for coarse MF estimation are larger than that of BASIC in Figure3, when the SNR of the postcorrelation signal is larger than −10 dB, the detection probability of the proposed method is almost 100%, which is larger than other FRFT and BAISIC in Figure 4. This is because FRFT has a search bin α and BAISIC adopts a differential process, which may degrade the correlation peak and lead to lower detection probability.Figure 3
Complex multiplications comparison for coarse MF estimation methods.Figure 4
Detection probabilities comparison when is 200 ms, chirping rate is 200 Hz/s, and initial frequency is 100 Hz.


## 5.2. Fine MF Accuracy and Complexity Comparison

Based on the simulation above, the DCFT method, two-step method, and Schmidt method are adopted for the fine MF search. In Figure5, the complex multiplications of Schmidt vary greatly with the change of MF search step and postcorrelation signal length.Figure 5
Complex multiplications comparison for fine MF estimation methods.In Figure6, even though 1 rad/s of MF search step is adopted for MF estimation, the two-step method based on MLEP gains higher precision than MF search based on Schmidt.Figure 6
MF square error comparison. The SNR of the postcorrelation signal is 5 dB. The MF search step is 1 rad/s. The integration time is 200 ms in coarse MF estimation.


## 6. Conclusion

To improve frequency accuracy in high-dynamic acquisition, we propose a two-step frequency estimation method. The proposed method combines a coarse frequency estimation method based on DCFT and MLEP for fine frequency estimation. In the 1st step, the search step of initial frequency and chirping rate is configured based on Taylor expansion, and coarse MF is obtained. In the 2nd step, due to low-frequency error, fine MF is estimated by MLEP. Although DCFT costs much more computation in Figure3 compared with BASIC, it improves the detection probability in Figure 4. Moreover, the proposed MLEP obtains higher mean frequency accuracy and lower complex multiplications compared with the conventional method Schmidt. Furthermore, in practice, the proposed two-step method can provide a theoretical basis for open-loop frequency tracking.



---

*Source: 1017206-2022-07-07.xml*