# Efficient and Accurate Frequency Estimator under Low SNR by Phase Unwrapping

Authors: Shen Zhou; Liu Rongfang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In the case of low signal-to-noise ratio, for the frequency estimation of single-frequency sinusoidal signals with additive white Gaussian noise, the phase unwrapping estimator usually performs poorly. In this paper, an efficient and accurate method is proposed to address this problem. Different from other methods, based on fast Fourier transform, the sampled signals are estimated with the variances approaching the Cramer-Rao bound, followed with the maximum likelihood estimation of the frequency. Experimental results reveal that our estimator has a better performance than other phase unwrapping estimators. Compared with the state-of-the-art method, our estimator has the same accuracy and lower computational complexity. Besides, our estimator does not have the estimation bias.


---

## Body

## 1. Introduction

Frequency estimation of a complex sinusoid is a fundamental problem in signal processing and has applications in many areas including communications, power spectrum estimation, array, and radar signal processing [1–8]. The general signal model is(1)yn=Aexp⁡j2πθ+fn+wnn=0,1,…,N-1=rnexp⁡j2πxnwhere w(n)=wI(n)+jwQ(n). wI(n) and wQ(n) are independent and normally distributed with zero mean and variance σ2. r(n) is the absolute value of y(n) and x(n) is the argument of y(n). The frequency f, the phase θ, and the amplitude A are deterministic but unknown constants, and N is the number of samples. f and θ are in [-1/2,1/2). The problem is to estimate f with a low computational complexity and statistically efficient estimator.For all the frequency estimators, there is a signal-to-noise ratio (SNR) threshold. When the SNR is lower than the threshold, the mean square error (MSE) of the estimated frequency no longer converges to the Cramer-Rao bound (CRB) [7]. The classical periodogram estimator [6] is widely considered to have the best performance and the lowest SNR threshold. However, the implementation of this estimator is complicated and may suffer from the resolution problem [7, 8]. A commonly used phase unwrapping estimator was first suggested by Kay [9]. Through calculating the first-order difference of the phase signal, the resulting signal resembles a moving average process and the parameters can be estimated by standard linear techniques. Kay’s estimator can attain the CRB in high SNR, while it performs poorly from the moderate SNR. To change this situation, researchers presented many improved phase unwrapping estimators [10–17]. The main drawback of these phase unwrapping estimators is that the SNR threshold still begins from a relatively high SNR and the performance does depend on the value of the frequency. Among the phase unwrapping estimators [9–18], the least squares phase unwrapping estimator (LSPUE) [18, 19] performs well under low SNR, but its computational complexity is too high. Further, an iterative method requiring Nlog2⁡N operations was suggested in [20], showing a similar performance to that of the periodogram estimator.In this paper, we propose a new phase unwrapping estimator which has the same performance as the periodogram estimator. The asymptotic variance is given and the choice of parameters is analyzed. The main contribution of this paper is that we improve the performance of the phase unwrapping estimator by using fast Fourier transform (FFT) and derive the asymptotic estimation variance. Compared with other phase unwrapping estimators, both the SNR threshold and the accuracy are improved. Compared with the state-of-the-art method, our estimator has the same accuracy and lower computational complexity. Moreover, unlike the state-of-the-art method, our estimator does not have the estimation bias.In Section2, based on FFT, the estimated signals whose variances approach the CRB are given. Then, the frequency is estimated by phase unwrapping and the asymptotic variance is derived. Two methods are suggested. In Section 3, the simulations that display the statistical performance of our estimators alongside the periodogram estimator and some other phase unwrapping estimators are provided.


## 2. Methods

The phase unwrapping estimator like Kay’s [9] performs poorly in low SNR and the main reason is that this kind of estimator is no longer accurate from the medium SNR [15]. Considering (1), we have(2)yn=Aexp⁡j2πθ+fn+wn=Aexp⁡j2πθ+fn1+wnexp⁡-j2πθ+fnAAssuming w′(n)=w(n)(exp⁡-j2πθ+fn/A)=wI′(n)+jwQ′(n), wI′(n) and wQ′(n) are independent and normally distributed with zero mean and variance σ2/A2. Then, y(n) can be expressed as(3)yn=Aexp⁡j2πθ+fn1+wI′n+jwQ′nThe argument of y(n), denoted x(n)=∠y(n)/2π, has the form(4)xn=θ+fn+unwhere u(n) is the phase nosie. Considering (3), when the SNR is high enough, we have the approximation(5)1+wI′n+jwQ′n≈1+jwQ′n≈exp⁡jwQ′nAccording to (3), (4), and (5), u(n) can be approximated as u(n)≈wQ′(n)/2π. Then, we have the approximated linear phase model(6)xn≈θ+fn+wQ′n2πKay’s estimator is derived based on the model (6). However, from the medium SNR, the approximation (5) is no longer accurate and the phase noise u(n) can no longer be approximated as white Gaussian noise wQ′(n)/2π. Then, estimators based on model (6) will not be accurate.To address this problem, a commonly used method is to improve the SNR before using the phase unwrapping estimator [10–12, 16, 17]. However, for these estimators, the SNR threshold still begins from a relatively high SNR. Besides, when f is close to ±1/2, the performance is very poor. Different from these estimators, our estimator is realized by three steps. First, we do a coarse search to narrow the range of the frequency to be estimated. Then, we improve the SNR by using the moving average filter [16]. Finally, we do a fine search by using the phase unwrapping estimator to obtain the estimated frequency. We will show that, in this way, the estimator can achieve the optimal SNR threshold and its performance is no longer influenced by the value of the frequency.


### 2.1. The Coarse Search

The sequenceYm=FFT(yn) is the N-point FFT of yn. First, we do a coarse search and find the parameter m^N∈0,1,…,N-1(7)m^N=argmaxm⁡Ym2According to [7, 20], the frequency can be written as(8)f=m^N+δNIn [20], it has been shown that δ is almost surely in [-1/2,1/2] when the SNR is larger than the SNR threshold. Therefore, in the analysis of the asymptotic variance, similar to [20] we assume δ∈[-1/2,1/2], which will not influence the result. A complex exponential signal with the frequency -m^N/N is given as(9)zn=exp⁡-j2πm^NnNMultiplying y(n) by z(n), we have the signal s(n)(10)sn=ynzn=Aexp⁡j2πΔfn+θ+vnwhere v(n)=wnexp⁡-j2πm^Nn/N=vI(n)+jvQ(n) and Δf=δ/N. v(n) is also a complex white Gaussian noise with zero-mean and variance 2σ2. As m^N has been obtained by (7), according to (8), to estimate f, all we need to do is to estimate Δf from the sequence sn. Then the estimated frequency can be obtained by(11)f^=m^NN+Δf^


### 2.2. Improving the SNR by the Moving Average Filter

Before estimatingΔf, we use the moving average filter to improve the SNR. {s(n)} is divided into L subsequences with the length M=N/L. For the kth subsequence, an estimated signal can be obtained(12)s^k=∑l=0M-1skM+lk=0,…,L-1Substituting (10) into (12) yields(13)s^k=ΩM,δAexp⁡j2k+1M-1πΔf+2πθ+∑l=0M-1vkM+lwhere(14)ΩM,δ=sin⁡MπΔfsin⁡πΔf=sin⁡Mπδ/Nsin⁡πδ/NThe argument of s^(k), denoted ∠s^(k)/2π, can be written as(15)xk=2kM+M-1Δf2+θ+ukwhere uk is the phase noise. According to the phase model in [21], uk can be approximated as(16)uk=1ΩM,δ∑l=0M-1vQkM+l2πAObviously, the variance of uk is(17)var⁡uk=MΩM,δ2σ22π2A2In the appendix, it is demonstrated that, on the condition of M≤N/8, var⁡(uk) is quite close to σ2/M2π2A2 which is the CRB for phase estimation. For y(n), the phase noise is approximately wQ(n)/2πA with the variance σ2/2π2A2 [21]. Therefore, for s^(k), the SNR is improved by 10log10⁡(M)dB approximately. Correspondingly, for the frequency estimation based on s^k, the SNR threshold can also be 10log10⁡(M)dB lower. That is why we can obtain a much better performance. In addition, var⁡(uk) increases with the decreasing δ and attains the upper bound when δ=±1/2.


### 2.3. The Fine Search Method

In the following, usings^k, Δf is estimated by Kay’s two phase unwrapping estimation methods [9]: weighted phase average (WPA) and weighted linear predictor (WLP). In this paper, we call our two methods FFT-based weighted phase average (FWPA) and FFT-based weighted linear predictor (FWLP), respectively.First, we introduce the FWPA. We realize the FWPA estimator through using the WPA estimator for{s^(k)}. According to (15), the first difference of the argument of s^(k) has the form(18)Δk=∠s^k+1s^∗k2π=MΔf+uk+1-ukk=0,…,L-2For (18), to estimate Δf, the minimum variance unbiased estimation is well known [9](19)Δf^=1M1TC-1Δ1TC-11=wTΔwhere 1=1,…,1T, w=1/M(1TC-1/1TC-11)=w0,…,wL-2T, Δ=[Δ0,…,ΔL-2], and C=[Ci,j] is the (L-1)×(L-1) covariance matrix of uk+1-uk(20)Ci,j=MΩM,δ2σ22π2A21,i=j-12,i-j=10,othersThe variance of Δf^ is [9](21)var⁡Δf^=6M2NN2-M21ΩM,δ22σ22π2A2and the CRB for the estimated frequency is [6](22)CRBΔf^=6NN2-12σ22π2A2According to (21), (22), (A.2), and (A.3), there is an upper bound for var⁡(Δf^)(23)var⁡Δf^≤N2-1N2-M2M2CRBΔf^ΩM,±1/22≈11-M/N2CRBΔf^1-M/2N2π2/3!2It can be seen that the upper bound only has the relation with 1/L=M/N. To keep the balance between a low upper bound and a high SNR, we usually set M=N/8 or M=N/16, for which the upper bounds are, respectively, ≤1.029CRB and ≤1.007CRB. By doing an iteration (when we obtain the estimated frequency f^, we can take f^ as the result of the coarse search, let m^N/N=f^ and perform the frequency estimation again), δ converges to 0, and the limit of var⁡(Δf^) is(24)var⁡Δf^=N2-1N2-M2M2CRBΔf^ΩM,02≈CRBΔf^1-M/N2For L=8 and L=16, the limits of var⁡(Δf^) are, respectively, 1.015CRB and 1.004CRB. Our estimator has an asymptotic variance which is only a little larger than the CRB.Then we introduce the FWLP. For the frequency estimation problem, there are two kinds of relatively complicated operations, namely, thesin/cos operation and the arctangent operation, so these operations should be reduced as much as possible. In FFT and (9), calculating a complex exponential needs two sin/cos operations. However, all of these complex exponentials take fixed values that can be stored in memory in advance. Hence, our estimator does not need sin/cos operations and, in order to decrease the number of arctangent operations, it is possible to find another estimator(25)Δf^=∠∑k=0L-2wks^k+1s^∗k2πThis is the FWLP and the only one complicated arithmetic operation is the arctangent, which can further reduce the computational complexity. In [9, 15], it has been shown that, for Kay’s method, the linear predictor has the same performance as the phase average only in very high SNR. However, for our method, the SNR is high enough to make the linear predictor nearly have the same performance as the phase average, which is shown in Section 3. The algorithm for the FWPA and the FWLP is summarized in Algorithm 1.Algorithm 1:The algorithm for FWPA and FWLP.
Step1: Let Ym=FFTyn Step2: Find m^N=argmaxm⁡Ym2 Step3: For each n from 0 to N-1 do s ( n ) = y n exp ⁡ - j 2 π m ^ N n N Step4: For each k from 0 to L-1 do s ^ k = ∑ l = 0 M - 1 s k M + l Step5: FWPA:Δf^=∑k=0L-2wk∠s^k+1s^∗k/2π FWLP:Δf^=∠∑k=0L-2wks^k+1s^∗k/2π Step6: f^=m^N/N+Δf^


### 2.4. Analysis of Computational Complexity

We assume that theN-samples FFT requires Nlog2⁡N complex valued (CV) multiplications and additions. Locating the DFT maximum requires additional 2N real-valued (RV) multiplications and N RV additions (calculation of the squared modulus of the DFT) and N comparisons (the worst case). The FWPA needs (4Nlog2⁡N+4N) RV additions, (4Nlog2⁡N+6N+L) RV multiplications, N comparisons, and Larctangent operations. The FWLP needs (4Nlog2⁡N+4N) RV additions, (4Nlog2⁡N+6N+2L) RV multiplications, N comparisons, and only one arctangent operation. Among the previous estimators achieving the optimal threshold, based on FFT, the iterative estimators [4, 20] and the direct estimators [3] have a lower computational complexity. According to the recent result, among these estimators, the estimator in [4] has the lowest computational complexity. This estimator needs (4Nlog2⁡N+13N) RV additions, (4Nlog2⁡N+17N) RV multiplications, 6Nsin/cos calculations, and N comparisons. Obviously, compared with the state-of-the-art method [4, 20] our estimator has the same accuracy and a lower computational complexity. Besides, for the iterative estimators and the direct estimators, there is an inevitable estimation bias [22], while our estimator is unbiased.


## 2.1. The Coarse Search

The sequenceYm=FFT(yn) is the N-point FFT of yn. First, we do a coarse search and find the parameter m^N∈0,1,…,N-1(7)m^N=argmaxm⁡Ym2According to [7, 20], the frequency can be written as(8)f=m^N+δNIn [20], it has been shown that δ is almost surely in [-1/2,1/2] when the SNR is larger than the SNR threshold. Therefore, in the analysis of the asymptotic variance, similar to [20] we assume δ∈[-1/2,1/2], which will not influence the result. A complex exponential signal with the frequency -m^N/N is given as(9)zn=exp⁡-j2πm^NnNMultiplying y(n) by z(n), we have the signal s(n)(10)sn=ynzn=Aexp⁡j2πΔfn+θ+vnwhere v(n)=wnexp⁡-j2πm^Nn/N=vI(n)+jvQ(n) and Δf=δ/N. v(n) is also a complex white Gaussian noise with zero-mean and variance 2σ2. As m^N has been obtained by (7), according to (8), to estimate f, all we need to do is to estimate Δf from the sequence sn. Then the estimated frequency can be obtained by(11)f^=m^NN+Δf^


## 2.2. Improving the SNR by the Moving Average Filter

Before estimatingΔf, we use the moving average filter to improve the SNR. {s(n)} is divided into L subsequences with the length M=N/L. For the kth subsequence, an estimated signal can be obtained(12)s^k=∑l=0M-1skM+lk=0,…,L-1Substituting (10) into (12) yields(13)s^k=ΩM,δAexp⁡j2k+1M-1πΔf+2πθ+∑l=0M-1vkM+lwhere(14)ΩM,δ=sin⁡MπΔfsin⁡πΔf=sin⁡Mπδ/Nsin⁡πδ/NThe argument of s^(k), denoted ∠s^(k)/2π, can be written as(15)xk=2kM+M-1Δf2+θ+ukwhere uk is the phase noise. According to the phase model in [21], uk can be approximated as(16)uk=1ΩM,δ∑l=0M-1vQkM+l2πAObviously, the variance of uk is(17)var⁡uk=MΩM,δ2σ22π2A2In the appendix, it is demonstrated that, on the condition of M≤N/8, var⁡(uk) is quite close to σ2/M2π2A2 which is the CRB for phase estimation. For y(n), the phase noise is approximately wQ(n)/2πA with the variance σ2/2π2A2 [21]. Therefore, for s^(k), the SNR is improved by 10log10⁡(M)dB approximately. Correspondingly, for the frequency estimation based on s^k, the SNR threshold can also be 10log10⁡(M)dB lower. That is why we can obtain a much better performance. In addition, var⁡(uk) increases with the decreasing δ and attains the upper bound when δ=±1/2.


## 2.3. The Fine Search Method

In the following, usings^k, Δf is estimated by Kay’s two phase unwrapping estimation methods [9]: weighted phase average (WPA) and weighted linear predictor (WLP). In this paper, we call our two methods FFT-based weighted phase average (FWPA) and FFT-based weighted linear predictor (FWLP), respectively.First, we introduce the FWPA. We realize the FWPA estimator through using the WPA estimator for{s^(k)}. According to (15), the first difference of the argument of s^(k) has the form(18)Δk=∠s^k+1s^∗k2π=MΔf+uk+1-ukk=0,…,L-2For (18), to estimate Δf, the minimum variance unbiased estimation is well known [9](19)Δf^=1M1TC-1Δ1TC-11=wTΔwhere 1=1,…,1T, w=1/M(1TC-1/1TC-11)=w0,…,wL-2T, Δ=[Δ0,…,ΔL-2], and C=[Ci,j] is the (L-1)×(L-1) covariance matrix of uk+1-uk(20)Ci,j=MΩM,δ2σ22π2A21,i=j-12,i-j=10,othersThe variance of Δf^ is [9](21)var⁡Δf^=6M2NN2-M21ΩM,δ22σ22π2A2and the CRB for the estimated frequency is [6](22)CRBΔf^=6NN2-12σ22π2A2According to (21), (22), (A.2), and (A.3), there is an upper bound for var⁡(Δf^)(23)var⁡Δf^≤N2-1N2-M2M2CRBΔf^ΩM,±1/22≈11-M/N2CRBΔf^1-M/2N2π2/3!2It can be seen that the upper bound only has the relation with 1/L=M/N. To keep the balance between a low upper bound and a high SNR, we usually set M=N/8 or M=N/16, for which the upper bounds are, respectively, ≤1.029CRB and ≤1.007CRB. By doing an iteration (when we obtain the estimated frequency f^, we can take f^ as the result of the coarse search, let m^N/N=f^ and perform the frequency estimation again), δ converges to 0, and the limit of var⁡(Δf^) is(24)var⁡Δf^=N2-1N2-M2M2CRBΔf^ΩM,02≈CRBΔf^1-M/N2For L=8 and L=16, the limits of var⁡(Δf^) are, respectively, 1.015CRB and 1.004CRB. Our estimator has an asymptotic variance which is only a little larger than the CRB.Then we introduce the FWLP. For the frequency estimation problem, there are two kinds of relatively complicated operations, namely, thesin/cos operation and the arctangent operation, so these operations should be reduced as much as possible. In FFT and (9), calculating a complex exponential needs two sin/cos operations. However, all of these complex exponentials take fixed values that can be stored in memory in advance. Hence, our estimator does not need sin/cos operations and, in order to decrease the number of arctangent operations, it is possible to find another estimator(25)Δf^=∠∑k=0L-2wks^k+1s^∗k2πThis is the FWLP and the only one complicated arithmetic operation is the arctangent, which can further reduce the computational complexity. In [9, 15], it has been shown that, for Kay’s method, the linear predictor has the same performance as the phase average only in very high SNR. However, for our method, the SNR is high enough to make the linear predictor nearly have the same performance as the phase average, which is shown in Section 3. The algorithm for the FWPA and the FWLP is summarized in Algorithm 1.Algorithm 1:The algorithm for FWPA and FWLP.
Step1: Let Ym=FFTyn Step2: Find m^N=argmaxm⁡Ym2 Step3: For each n from 0 to N-1 do s ( n ) = y n exp ⁡ - j 2 π m ^ N n N Step4: For each k from 0 to L-1 do s ^ k = ∑ l = 0 M - 1 s k M + l Step5: FWPA:Δf^=∑k=0L-2wk∠s^k+1s^∗k/2π FWLP:Δf^=∠∑k=0L-2wks^k+1s^∗k/2π Step6: f^=m^N/N+Δf^


## 2.4. Analysis of Computational Complexity

We assume that theN-samples FFT requires Nlog2⁡N complex valued (CV) multiplications and additions. Locating the DFT maximum requires additional 2N real-valued (RV) multiplications and N RV additions (calculation of the squared modulus of the DFT) and N comparisons (the worst case). The FWPA needs (4Nlog2⁡N+4N) RV additions, (4Nlog2⁡N+6N+L) RV multiplications, N comparisons, and Larctangent operations. The FWLP needs (4Nlog2⁡N+4N) RV additions, (4Nlog2⁡N+6N+2L) RV multiplications, N comparisons, and only one arctangent operation. Among the previous estimators achieving the optimal threshold, based on FFT, the iterative estimators [4, 20] and the direct estimators [3] have a lower computational complexity. According to the recent result, among these estimators, the estimator in [4] has the lowest computational complexity. This estimator needs (4Nlog2⁡N+13N) RV additions, (4Nlog2⁡N+17N) RV multiplications, 6Nsin/cos calculations, and N comparisons. Obviously, compared with the state-of-the-art method [4, 20] our estimator has the same accuracy and a lower computational complexity. Besides, for the iterative estimators and the direct estimators, there is an inevitable estimation bias [22], while our estimator is unbiased.


## 3. Results and Discussion

This section shows the simulation results to illustrate the performance of our estimators. First, we compare the performance of our FWPA estimator, our FWLP estimator, the periodogram estimator [6], the LSPUE [18], the WPA estimator [9], the WLP estimator [9], and the hybrid estimator [17]. Among the improved estimators [10–17] based on Kay’s estimator [9], the hybrid estimator has the best performance. Compared with other phase unwrapping estimators, the LSPUE has a much better performance. Therefore, we compare our estimator with the hybrid estimator and the LSPUE. The performance was evaluated by computer simulation in complex white Gaussian noise and the SNR is 10log10⁡A2/2σ2dB. The SNR was incremented from -20 dB to 30 dB in steps of 2 dB. To ensure the accuracy, 10000 trials were run for each SNR value.Figures1 and 2 show the MSE of different estimators. In Figure 1, the parameters are f=0.012, θ=0.35 and the number of samples is N=64. In Figure 2, the parameters are f=0.4, θ=0.15 and the number of samples is N=1024. It is clear that the periodogram estimator, the LSPUE, our FWPA estimator, and our FWLP estimator have the best performance, while our two estimators have a lower computation complexity. The WPA estimator, the WLP estimator, and the hybrid estimator perform comparatively poorly especially when N is large or f is close to ±1/2. From the medium SNR, the LSPUE has an asymptotic variance which is a little larger than the CRB, while the variances of the periodogram estimator and our FWPA estimator approach the CRB accurately. Besides, when N is small (as is shown in Figure 1, N is 64.), the SNR threshold of the LSPUE is larger than that of the periodogram estimator while our estimators still possess the optimal threshold. The WLP has a much worse performance than WPA, while, for our estimators, the performance of FWLP is only marginally worse than that of FWPA.Figure 1
The MSE in frequency versus SNR for different estimators whenf=0.012,θ=0.35,N=64.Figure 2
The MSE in frequency versus SNR for different estimators whenf=0.4,θ=0.15,N=1024.Then, we compare the performance of our FWPA estimator and the A&M estimator [4, 20]. Among the iterative estimators, the A&M estimator has the best performance and is widely considered to be the state-of-the-of-art method. Considering that the performance of the two kinds of estimators can be influenced by the value of δ, we compare the performance when δ varied from -0.5 to 0.5 in particular SNR. Figures 3 and 4 show the MSE of different estimators when δ varied from -0.5 to 0.5. In Figure 3, N=1024 and SNR=-5dB. In Figure 4, N=1024 and SNR=5dB. Q is the number of iterations. It can be seen that the performance of FWPA estimator is not influenced by frequency variation. For the A&M estimator, an iteration must be taken, while, for the FWPA estimator, we do not need to take an iteration. Besides, the FWPA estimator has a better performance than the A&M estimator when δ is close to ±0.5.Figure 3
The MSE in frequency for different estimators whenδ varied from -0.5 to 0.5, N=1024, and SNR=-5dB.Figure 4
The MSE in frequency for different estimators whenδ varied from -0.5 to 0.5, N=1024, and SNR=5dB.


## 4. Conclusion

Phase unwrapping frequency estimator usually has a bad performance in low SNR. To solve this problem, we propose a new estimator. By improving the SNR before using the phase unwrapping estimator, the new estimator performs well in low SNR and has the optimal threshold. Compared with the LSPUE, it has a better performance and the computational complexity is reduced greatly. Compared with other phase unwrapping estimators, it has a much better performance and can well solve the problem of bad performance under low SNR. Compared with the state-of-the-art method, it has the same accuracy and a lower computational complexity. Moreover, unlike the state-of-the-art method, our estimator does not have the estimation bias. Due to its simplicity, efficiency, and low computational complexity, the proposed estimator represents a viable solution for real-time practical applications.



---

*Source: 7396074-2019-04-14.xml*