# Sinc Interpolation for Further Improvement in Frequency Estimation Using Three DFT Samples

Authors: Hasan Bayazit; Erdoğan Dilaveroğlu
Journal: Journal of Sensors
 (2022)

---

## Abstract

This study presents three different sinc estimators and a method to estimate complex-valued exponential tone signal frequency using three DFT samples. The proposed method suggests using sinc interpolation together with the well-known Jacobsen estimator. According to simulation results, the root mean square error (RMSE) of the proposed algorithm is lower than those of the Jacobsen estimator and its improved version suggested by Candan. The price paid for improvements in the RMSE is a slight increase in computation time.


---

## Body

## 1. Introduction

The frequency estimation problem plays an essential role in communications, audio, medical, instrumentation, and other applications. Nonparametric and parametric techniques are the main tools used in spectrum estimation methods. Nonparametric techniques do not assume a particular functional form but allow the form of the estimator to be determined entirely by the data. These methods are based on the discrete Fourier transform of either the signal segment or its autocorrelation sequence. In contrast, parametric methods assume that the available signal segment has been generated by a specific parametric model [1] (p. 195). Many nonparametric and parametric methods, such as DFT-based frequency shifting and filtering, phase-locked loop, maximum likelihood, wavelet, Prony’s, Taylor, interpolation, and of multi-DFT-bins, artificial neural network design, have been published [2–5].In digital signal processing, the discrete Fourier transform (DFT) is a frequently used tool to detect and estimate the frequencies of a signal. Several estimation methods based on DFT have been proposed. The well-known method for frequency estimation is to use N-point Fourier transform after sampling a signal. After applying an N-point DFT, the frequency of interest rarely resides exactly on a DFT bin center. In this case, frequency estimation resolution depends on the bin spacing of the transform. However, we always need better frequency estimation resolution. There are lots of frequency estimation methods to improve that resolution. These methods have been described in [6–21].Equation (1) is the best-known Jacobsen peak location estimator and has been used in the DSP business for decades, which uses three DFT samples [22] (p. 733). This spectral peak location estimator δ which is given by
(1)δ=ReXkp−1−Xkp+12Xkp−Xkp−1−Xkp+1,considers the bin with the maximum DFT value, Xkp and its adjacent right and left bin DFT values, Xkp−1 and Xkp+1, respectively, to calculate a correction factor δ for the peak location (see Figure 1).Figure 1
DTFT and DFT of a complex tone in the absence of noise: (a)delta>0; (b) delta<0.
(a)(b)The frequency is estimated askp+δ with the units of DFT bins or f^=kp+δ/N in terms of the normalized frequency. A bias-corrected version of the Jacobsen estimator, proposed by Candan [10], is given in
(2)δ=tanπ/Nπ/NReXkp−1−Xkp+12Xkp−Xkp−1−Xkp+1.Clearly, for a large number of data samples (asN⟶∞), the two estimators coincide.These estimators perform well when the actual signal frequency happens to be in the center of the bin. When the actual frequency moves toward either end of the bin, RMS error of the estimators increases. In this paper, we developed estimators based on sinc function interpolations and measured their performances to find an answer to the question, how could these sinc function-based estimators be used to increase the performance of the Jacobsen estimator (1) and its bias-corrected version (2).


## 2. Preliminaries

A complex signal with an amplitudeVm>0, a phase θ∈0,2π, and a normalized frequency f can be expressed as
(3)xn=Vm.ej2πfn+θ,n=0,⋯,N−1,where N is the total number of samples available. The signal frequency f in bin values is defined as
(4)f=kpeak/N,(see Figure 1). The DFT of the sequence xn is
(5)Xk=∑n=0N−1xn.e−j2π/Nkn,and can be expressed as
(6)Xk=Vm.ejθ.ejπN−1/Nkpeak−k.sinπkpeak−ksinπ/Nkpeak−k.DFT of thexn is given by Equation (6). The spectrum would be continuous if N samples were increased infinitely. But the DFT works for finite-length signals. Selecting more samples means taking a longer signal and multiplying the signal by a rectangular pulse. Since the multiplication of a signal in the time domain results in the frequency domain being convolved with a sinc function, this reduces frequency resolution [23–27]. Therefore, some resolution improvements are needed for better estimation.The DFT of the signal forms equally spacedN/2 samples in the frequency spectrum. Each spectral value is at a different interval of k.fs/N. These frequency intervals fs/N define frequency resolution and k.fs/N are denoted frequency bins. Frequency bins are integer values from 0 up to N−1. Bin kp has the largest DFT magnitude value (Figure 1). Bins next to the kp bin, kp+1, and the previous bin, kp−1, are considered in estimating the tone frequency. The actual frequency could be anywhere in the bin kp, depending on the bin kp+1 and kp−1 magnitude values. And this exact frequency can be found by adding or subtracting the δ correction value to the bin kp, since δ changes in the range of ±1/2.The DFT magnitudes ofXkp−1, Xkp, and Xkp+1 for Vm=1V,θ=0, N=32, and frequency f sitting anywhere in the tenth bin are shown in Figure 2. The signal frequencies in bin 10±1/2 are all represented in bin 10. If the actual frequency f happens to be in the middle of bin 10, the magnitude of Xkp becomes maximum, and the previous and next bin’s magnitude values become zero. On the other hand, when the actual frequency goes toward either end of the bin, the magnitude of Xkp decreases while magnitudes of Xkp+1 and Xkp−1 increase. A closer look will notice how symmetry exists between Xkp−1 and Xkp+1. Since the magnitude margins increase toward the bin ends, any estimator which calculates the correction factor value just by using the magnitude values of Xkp−1 and Xkp+1 would work well, while there is no doubt that it will not work well where the tone frequencies are close to the bin center.Figure 2
DFT magnitudesXkp−1, Xkp, and Xkp+1 vs. actual bin frequency in the absence of noise.


## 3. Sinc-Based Estimators

Interpolation is a method of finding new data points based on a set of known data. These data points represent the values of a function for a limited value of the independent variable. It is often required to estimate the value of the function for an intermediate value of the independent variable. The DTFT of a single tone has the characteristics of a sinc function. After performing the DFT of the signal, all the DFT values of transformation will be on the continuous spectrum of Fourier transformation in the noise-free cases. We can find the correction factorδ by solving the sinc function equations defined by the maximum and the second maximum DFT magnitudes, which is the main idea behind sinc function-based estimators [28–30].


### 3.1. The First Sinc Estimator

The location of the signal frequency we are trying to estimate is between the largest amplitude and the second biggest one adjacent to it. The position of the second largest magnitude could be before or next to the largest amplitude. Both cases should be investigated to find a solution. Let us start to analyze the situation depicted in Figure1(a), where the second biggest amplitude comes after the maximum. In this case, DFT values at instances kp and kp+1 can be found from Equation (6) after substitutions k=kpeak−δ, and k=kpeak+1−δ, respectively. The following equations are obtained:
(7)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπδsinπδ/N,(8)Xkp+1=Vm.ejθ.ejπ/NN−1.δ−1.sinπδ−1sinπ/Nδ−1.The magnitude values of (7) and (8) are as follows:
(9)Xkp=Vm.sinπδsinπδ/N,(10)Xkp+1=Vm.sinπδ−1sinπ/Nδ−1.Equations (9) and (10) can be written as
(11)Xkp.sinπδN−Vm.sinπδ=0,(12)Xkp+1.sinπNδ−1−Vm.sinπδ−1=0.If we solve Equations (11) and (12) for δ, we get the following equation:
(13)tanπ.δN=Xkp+1.sinπ/NXkp+Xkp+1.cosπ/N.Let us take the Taylor series expansion of the functiontanπ.δ/N around δ=0 to find a closed-form equation for δ.
(14)tanπ.δN=π.δN+π3.δ33.N3+2.π5.δ515.N5+⋯.Assuming smallδ and/or large N, we obtain the following estimator
(15)δ^1≅Nπ.Xkp+1.sinπ/NXkp+Xkp+1.cosπ/N,(16)δ^1≅Xkp+1Xkp+Xkp+1asN⟶∞,after taking the first term of the expansion.On the other end, let us find out how the estimator behaves in the situation depicted in Figure1(b), where the second biggest magnitude comes before the largest one. In this case, DFT values at kp and kp−1 can be found from Equation (6) after substitutions k=kpeak−δ, and k=kpeak−δ+1, respectively. δ is defined as a negative variable in this region. The following equations are obtained:
(17)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπδsinπδ/N,(18)Xkp−1=Vm.ejθ.ejπ/NN−1.δ+1.sinπδ+1sinπ/Nδ+1.After taking the magnitudes of Equations (17) and (18) and solving for δ, we obtain the following equations:
(19)δ^1≅Nπ.−Xkp−1.sinπ/NXkp+Xkp−1.cosπ/N,δ^1≅−Xkp−1Xkp+Xkp−1,asN⟶∞.


### 3.2. The Second Sinc Estimator

Now, let us find an estimator that calculates the correction factor by looking at the DFT magnitude values adjacent to the maximum DFT value. We have to define the magnitude values ofXkp+1 and Xkp−1 for δ>0 to get a solution
(20)Xkp+1=Vm.sinπδ−1sinπ/Nδ−1,(21)Xkp−1=−Vm.sinπδ+1sinπ/Nδ+1.If we solve Equations (20) and (21) for δ, we get the following equation:
(22)δ^2≅Nπ.Xkp+1−Xkp−1Xkp+1+Xkp−1.tanπN,(23)δ^2≅Xkp+1−Xkp−1Xkp+1+Xkp−1,asN⟶∞.The magnitude values ofXkp+1 and Xkp−1 for δ<0 are defined as follows:
(24)Xkp+1=−Vm.sinπδ−1sinπ/Nδ−1,(25)Xkp−1=Vm.sinπδ+1sinπ/Nδ+1.The solution of Equations (24) and (25) for δ gives the same result as in Equation (22).


### 3.3. The Third Sinc Estimator

The previous two estimators compute the correction factorδ from the DFT magnitude values. Now, let us find an estimator that calculates the correction factor by looking at the maximum and second maximum complex DFT values. We have to go through similar steps as before. We obtain the following equations for the situation depicted in Figure 1(a). where δ is defined as positive
(26)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπ.δsinπδ/N=Um.sinπ.δsinπδ/N,(27)Xkp+1=Vm.ejθ.ejπ/NN−1.δ−1.sinπδ−1sinπ/Nδ−1≅−Um.sinπδ−1sinπ/Nδ−1,asN⟶∞.Equations (26) and (27) can be rewritten as follows:
(28)Xkp.sinπδN−Um.sinπ.δ=0,(29)Xkp+1.sinπNδ−1+Um.sinπ.δ−1=0.We obtain the third estimatorδ3 after solving Equations (28) and (29) as follows:
(30)δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/N,(31)δ^3≅ReXkp+1−Xkp+Xkp+1,asN⟶∞.Now, let us find the estimator for the situation depicted in Figure1(b). We should repeat the same steps to find DFT complex values at kp and kp−1, which can be found from substitutions k=kpeak−δ and k=kpeak−1+δ, where δ is negative. The following equations are obtained:
(32)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπ.δsinπδ/N=Um.sinπ.δsinπδ/N,(33)Xkp−1=Vm.ejθ.ejπ/NN−1.δ+1.sinπδ+1sinπ/Nδ+1≅−Um.sinπδ+1sinπ/Nδ+1,asN⟶∞.After solving Equations (32) and (33) for δ, we obtain the following equations:
(34)δ^3≅ReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/N,δ^3≅Re−Xkp−1−Xkp+Xkp−1,asN⟶∞,which are a similar estimator as in (30) except for a sign difference.


### 3.4. Decision Criteria

Table1 summarizes all the sinc-function-based estimators. The first and third estimators consist of two parts. Before computing the correction factor, there must be a decision criterion for identifying which part of the expression should be used. One of the criteria would be a comparison of magnitude values of Xkp−1 and Xkp+1 as
(35)α1=Xkp+1>Xkp−1.Table 1
Summary of estimators.
Estimatorsδ^1≅Nπ.Xkp+1.sinπ/NXkp+Xkp+1.cosπ/Nforα2>0;ifα1istrueNπ.−Xkp−1.sinπ/NXkp+Xkp−1.cosπ/Nforα2<0;ifα1isfalseδ^2≅Nπ.Xkp+1−Xkp−1Xkp+1+Xkp−1.tanπNforanyα2;foranyα1δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/Nforα2>0;ifα1istrueReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/Nforα2<0;ifα1isfalseDecision criteriaα1=Xkp+1>Xkp−1;α2=ReXkp−1−Xkp+1/XkpWhenXkp+1is equal to Xkp−1, δ equals to a bin frequency.If the signal frequency is near the bin center, the error of the estimators would be high due to the declining margin betweenXkp−1 and Xkp+1, as shown in Figure 2. However, the criterion becomes more stable when the signal frequency goes toward the bin ends.Figure3 depicts how the DFT complex values of Xkp−1, Xkp, and Xkp+1 varies vs. bin frequency for a complex signal having an amplitude of Vm=1V, θ=0, and a frequency changing from 9.5 to 10.5. Figure 3 shows that Xkp has the value of 32 when N=32, while Xkp−1 and Xkp+1 are zero when the signal frequency is in the bin center. Xkp starts decreasing as the actual signal frequency goes from bin 10 to bin 10+1/2, which means that the second biggest DFT magnitude value comes after the maximum value. Meanwhile, Xkp−1 and Xkp+1 start increasing. When bin frequency goes from bin 10 to bin 10.5, the contour line of Xkp begins at point A and goes to point B. At the same time, Xkp−1 and Xkp+1 contour lines start at the origin, as the first one follows the path toward D, and the second one goes toward point C in the opposite direction. When bin frequency goes from bin 10 to bin 9.5, the contour line of Xkp starts at point A, but goes toward point C. Similarly, Xkp−1 and Xkp+1 contour lines start at the origin, while the first one follows the path toward B, and the second one goes toward point E. In Figure 3, the complex DFT values of the maximum bin and its adjacent bins are depicted in circles for the actual frequency of bin 10.25. On the other hand, the complex DFT values are indicated in squares for the actual frequency of bin 9.75. The phase of the DFT values vs. actual bin frequency is depicted in Figure 4, which shows that the rate of phase change for all DFT values is equal to each other. As δ goes toward a positive end, the phase change of Xkp and Xkp−1 is identical but has opposite signs with Xkp+1. The same situation can be seen as the δ goes toward a negative end.Figure 3
The variation in the complex DFT values as bin frequency changes from bin frequency 9.5 to bin frequency 10.5 in the absence of noise.Figure 4
Phase of the complex DFT values vs. the actual bin frequency in the absence of noise.Let us analyze Figure3 to find another decision criterion for the estimators δ^1 and δ^3. When the actual signal frequency is in bin 10.25, the complex DFT values of bin ten and its adjacent bins are shown in Figure 3. They have the following complex values: Xkp=28.81∠0.7608,Xkp−1=5.77∠0.6626, and Xkp+1=9.61∠−2.2826. If we add the negative of Xkp+1 to Xkp−1 and divide the sum by Xkp, we get 0.531∠0.024, which has a positive real value. The result of this operation is always positive for all δ>0 because both the nominator and denominator have similar signs. If we make a similar calculation for the actual frequency of bin 9.75, we end up with negative real values for all negative values of δ because both the nominator and denominator have opposite signs. For this reason, we propose the following decision criterion:
(36)α2=ReXkp−1−Xkp+1Xkp.Figure5 shows a comparison of the decision criteria between α1 and α2 for N=8 at SNR=8 dB. This comparison was made by looking at the values of Expressions (35) and (36). If α1 and α2 are true and positive, a positive sign is assigned for this operation; meanwhile, a negative sign is assigned for false and minus values. Both decision criteria should give a positive sign for δ>0 but a negative sign for δ<0. Figure 5 shows that false decisions mainly concentrate around δ=0, which is a cause of RMSE deterioration. The decision criterion α2 outperforms α1 in numbers of correct decisions even in small sample sizes and low SNR levels.Figure 5
Comparison of decision criteriaα1 and α2 vs. δ.Figure6 shows a comparison of estimators δ^1 and δ^3 with the two different criteria in terms of RMSE vs. SNR. Two different criteria were applied to each estimator simultaneously. The RMS errors of the estimators which use the criterion α1 are shown with dashed lines. Solid lines represent the RMS errors of the estimators which use the criterion α2. This comparison indicates that estimators δ^1 and δ^3 have better performance with decision criterion 2.Figure 6
RMSE of sinc function estimators with the two different decision criteria vs. SNR forN=8, δ=0.10.


## 3.1. The First Sinc Estimator

The location of the signal frequency we are trying to estimate is between the largest amplitude and the second biggest one adjacent to it. The position of the second largest magnitude could be before or next to the largest amplitude. Both cases should be investigated to find a solution. Let us start to analyze the situation depicted in Figure1(a), where the second biggest amplitude comes after the maximum. In this case, DFT values at instances kp and kp+1 can be found from Equation (6) after substitutions k=kpeak−δ, and k=kpeak+1−δ, respectively. The following equations are obtained:
(7)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπδsinπδ/N,(8)Xkp+1=Vm.ejθ.ejπ/NN−1.δ−1.sinπδ−1sinπ/Nδ−1.The magnitude values of (7) and (8) are as follows:
(9)Xkp=Vm.sinπδsinπδ/N,(10)Xkp+1=Vm.sinπδ−1sinπ/Nδ−1.Equations (9) and (10) can be written as
(11)Xkp.sinπδN−Vm.sinπδ=0,(12)Xkp+1.sinπNδ−1−Vm.sinπδ−1=0.If we solve Equations (11) and (12) for δ, we get the following equation:
(13)tanπ.δN=Xkp+1.sinπ/NXkp+Xkp+1.cosπ/N.Let us take the Taylor series expansion of the functiontanπ.δ/N around δ=0 to find a closed-form equation for δ.
(14)tanπ.δN=π.δN+π3.δ33.N3+2.π5.δ515.N5+⋯.Assuming smallδ and/or large N, we obtain the following estimator
(15)δ^1≅Nπ.Xkp+1.sinπ/NXkp+Xkp+1.cosπ/N,(16)δ^1≅Xkp+1Xkp+Xkp+1asN⟶∞,after taking the first term of the expansion.On the other end, let us find out how the estimator behaves in the situation depicted in Figure1(b), where the second biggest magnitude comes before the largest one. In this case, DFT values at kp and kp−1 can be found from Equation (6) after substitutions k=kpeak−δ, and k=kpeak−δ+1, respectively. δ is defined as a negative variable in this region. The following equations are obtained:
(17)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπδsinπδ/N,(18)Xkp−1=Vm.ejθ.ejπ/NN−1.δ+1.sinπδ+1sinπ/Nδ+1.After taking the magnitudes of Equations (17) and (18) and solving for δ, we obtain the following equations:
(19)δ^1≅Nπ.−Xkp−1.sinπ/NXkp+Xkp−1.cosπ/N,δ^1≅−Xkp−1Xkp+Xkp−1,asN⟶∞.


## 3.2. The Second Sinc Estimator

Now, let us find an estimator that calculates the correction factor by looking at the DFT magnitude values adjacent to the maximum DFT value. We have to define the magnitude values ofXkp+1 and Xkp−1 for δ>0 to get a solution
(20)Xkp+1=Vm.sinπδ−1sinπ/Nδ−1,(21)Xkp−1=−Vm.sinπδ+1sinπ/Nδ+1.If we solve Equations (20) and (21) for δ, we get the following equation:
(22)δ^2≅Nπ.Xkp+1−Xkp−1Xkp+1+Xkp−1.tanπN,(23)δ^2≅Xkp+1−Xkp−1Xkp+1+Xkp−1,asN⟶∞.The magnitude values ofXkp+1 and Xkp−1 for δ<0 are defined as follows:
(24)Xkp+1=−Vm.sinπδ−1sinπ/Nδ−1,(25)Xkp−1=Vm.sinπδ+1sinπ/Nδ+1.The solution of Equations (24) and (25) for δ gives the same result as in Equation (22).


## 3.3. The Third Sinc Estimator

The previous two estimators compute the correction factorδ from the DFT magnitude values. Now, let us find an estimator that calculates the correction factor by looking at the maximum and second maximum complex DFT values. We have to go through similar steps as before. We obtain the following equations for the situation depicted in Figure 1(a). where δ is defined as positive
(26)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπ.δsinπδ/N=Um.sinπ.δsinπδ/N,(27)Xkp+1=Vm.ejθ.ejπ/NN−1.δ−1.sinπδ−1sinπ/Nδ−1≅−Um.sinπδ−1sinπ/Nδ−1,asN⟶∞.Equations (26) and (27) can be rewritten as follows:
(28)Xkp.sinπδN−Um.sinπ.δ=0,(29)Xkp+1.sinπNδ−1+Um.sinπ.δ−1=0.We obtain the third estimatorδ3 after solving Equations (28) and (29) as follows:
(30)δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/N,(31)δ^3≅ReXkp+1−Xkp+Xkp+1,asN⟶∞.Now, let us find the estimator for the situation depicted in Figure1(b). We should repeat the same steps to find DFT complex values at kp and kp−1, which can be found from substitutions k=kpeak−δ and k=kpeak−1+δ, where δ is negative. The following equations are obtained:
(32)Xkp=Vm.ejθ.ejπ/NN−1.δ.sinπ.δsinπδ/N=Um.sinπ.δsinπδ/N,(33)Xkp−1=Vm.ejθ.ejπ/NN−1.δ+1.sinπδ+1sinπ/Nδ+1≅−Um.sinπδ+1sinπ/Nδ+1,asN⟶∞.After solving Equations (32) and (33) for δ, we obtain the following equations:
(34)δ^3≅ReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/N,δ^3≅Re−Xkp−1−Xkp+Xkp−1,asN⟶∞,which are a similar estimator as in (30) except for a sign difference.


## 3.4. Decision Criteria

Table1 summarizes all the sinc-function-based estimators. The first and third estimators consist of two parts. Before computing the correction factor, there must be a decision criterion for identifying which part of the expression should be used. One of the criteria would be a comparison of magnitude values of Xkp−1 and Xkp+1 as
(35)α1=Xkp+1>Xkp−1.Table 1
Summary of estimators.
Estimatorsδ^1≅Nπ.Xkp+1.sinπ/NXkp+Xkp+1.cosπ/Nforα2>0;ifα1istrueNπ.−Xkp−1.sinπ/NXkp+Xkp−1.cosπ/Nforα2<0;ifα1isfalseδ^2≅Nπ.Xkp+1−Xkp−1Xkp+1+Xkp−1.tanπNforanyα2;foranyα1δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/Nforα2>0;ifα1istrueReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/Nforα2<0;ifα1isfalseDecision criteriaα1=Xkp+1>Xkp−1;α2=ReXkp−1−Xkp+1/XkpWhenXkp+1is equal to Xkp−1, δ equals to a bin frequency.If the signal frequency is near the bin center, the error of the estimators would be high due to the declining margin betweenXkp−1 and Xkp+1, as shown in Figure 2. However, the criterion becomes more stable when the signal frequency goes toward the bin ends.Figure3 depicts how the DFT complex values of Xkp−1, Xkp, and Xkp+1 varies vs. bin frequency for a complex signal having an amplitude of Vm=1V, θ=0, and a frequency changing from 9.5 to 10.5. Figure 3 shows that Xkp has the value of 32 when N=32, while Xkp−1 and Xkp+1 are zero when the signal frequency is in the bin center. Xkp starts decreasing as the actual signal frequency goes from bin 10 to bin 10+1/2, which means that the second biggest DFT magnitude value comes after the maximum value. Meanwhile, Xkp−1 and Xkp+1 start increasing. When bin frequency goes from bin 10 to bin 10.5, the contour line of Xkp begins at point A and goes to point B. At the same time, Xkp−1 and Xkp+1 contour lines start at the origin, as the first one follows the path toward D, and the second one goes toward point C in the opposite direction. When bin frequency goes from bin 10 to bin 9.5, the contour line of Xkp starts at point A, but goes toward point C. Similarly, Xkp−1 and Xkp+1 contour lines start at the origin, while the first one follows the path toward B, and the second one goes toward point E. In Figure 3, the complex DFT values of the maximum bin and its adjacent bins are depicted in circles for the actual frequency of bin 10.25. On the other hand, the complex DFT values are indicated in squares for the actual frequency of bin 9.75. The phase of the DFT values vs. actual bin frequency is depicted in Figure 4, which shows that the rate of phase change for all DFT values is equal to each other. As δ goes toward a positive end, the phase change of Xkp and Xkp−1 is identical but has opposite signs with Xkp+1. The same situation can be seen as the δ goes toward a negative end.Figure 3
The variation in the complex DFT values as bin frequency changes from bin frequency 9.5 to bin frequency 10.5 in the absence of noise.Figure 4
Phase of the complex DFT values vs. the actual bin frequency in the absence of noise.Let us analyze Figure3 to find another decision criterion for the estimators δ^1 and δ^3. When the actual signal frequency is in bin 10.25, the complex DFT values of bin ten and its adjacent bins are shown in Figure 3. They have the following complex values: Xkp=28.81∠0.7608,Xkp−1=5.77∠0.6626, and Xkp+1=9.61∠−2.2826. If we add the negative of Xkp+1 to Xkp−1 and divide the sum by Xkp, we get 0.531∠0.024, which has a positive real value. The result of this operation is always positive for all δ>0 because both the nominator and denominator have similar signs. If we make a similar calculation for the actual frequency of bin 9.75, we end up with negative real values for all negative values of δ because both the nominator and denominator have opposite signs. For this reason, we propose the following decision criterion:
(36)α2=ReXkp−1−Xkp+1Xkp.Figure5 shows a comparison of the decision criteria between α1 and α2 for N=8 at SNR=8 dB. This comparison was made by looking at the values of Expressions (35) and (36). If α1 and α2 are true and positive, a positive sign is assigned for this operation; meanwhile, a negative sign is assigned for false and minus values. Both decision criteria should give a positive sign for δ>0 but a negative sign for δ<0. Figure 5 shows that false decisions mainly concentrate around δ=0, which is a cause of RMSE deterioration. The decision criterion α2 outperforms α1 in numbers of correct decisions even in small sample sizes and low SNR levels.Figure 5
Comparison of decision criteriaα1 and α2 vs. δ.Figure6 shows a comparison of estimators δ^1 and δ^3 with the two different criteria in terms of RMSE vs. SNR. Two different criteria were applied to each estimator simultaneously. The RMS errors of the estimators which use the criterion α1 are shown with dashed lines. Solid lines represent the RMS errors of the estimators which use the criterion α2. This comparison indicates that estimators δ^1 and δ^3 have better performance with decision criterion 2.Figure 6
RMSE of sinc function estimators with the two different decision criteria vs. SNR forN=8, δ=0.10.


## 4. Simulation Results

This section compares sinc estimators, as stated in Table1, with the Jacobsen estimator (1) and its modified version suggested by Candan (2) in terms of complex and real signals.


### 4.1. Simulation Results for Complex Signals

During numerical simulations, a complex signal with an amplitude of 1 V, a randomly changing phase between0 and 2π, and a frequency f=kpeak/N have been used. The signal has also been contaminated with a white Gaussian noise
(37)xn=Vm.ej2πkpeak/Nn+θ+wn,n=0,⋯,N−1.Herewn is a circularly symmetric complex white Gaussian noise with mean zero and variance σw2. The signal-to-noise ratio is defined as SNR=Vm2/σw2.Figures7 and 8 examine how sample size N and SNR levels affect estimators’ behavior for the signal frequency at bin ten with an offset of δ=0.35. Both figures show that the RMSE decreases as the sample size increases. The first and third sinc estimators outperform the other estimators at all levels of SNR values.Figure 7
Comparison of RMS errors of the estimators forN=64 and δ=0.35.Figure 8
Comparison of RMS errors of the estimators forN=32 and δ=0.35.Figure9 compares the sinc estimators with the estimators of Jacobsen and Candan. The number of samples was reduced to 8, but the signal frequency was kept the same as in Figures 7 and 8. Figure 9 shows that the sinc estimators perform well even at a small sample size of N=8. Also, RMSE improvements are so distinct at low SNR values.Figure 9
RMS errors of the estimators vs. SNR forN=8 and δ=0.35.Table2 shows how much execution time is needed for each estimator at bin frequency 10 for a sample size of N=32. Execution times have been measured in a loop and then averaged to find the time for a single run. The price to be paid for improvements in RMSE is roughly a microsecond.Table 2
The proposed frequency estimation method for small sample values.
Step 1(i) Search the peak sample from the N-point DFT and its left and right neighborsXkp,Xkp−1, and Xkp+1Step 2(i) Ifδ calculated from (1) is δ≥0.2, recalculate the correction factor from the following equation:δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/NStep 3(i) Ifδ calculated from (1) is δ≤−0.2, recalculate the correction factor from the following equation:δ^3≅ReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/NStep 4(i) Calculatekpeak and f as follows:kpeak=kp+δ and f=kpeak/N, for normalized frequencyFigure10 expresses the RMS errors as a function of the correction factor δ when the sample size N=8 and SNR=8 dB. It shows that as the frequency, we are trying to estimate is at a location toward the bin edges; the RMS errors of the sinc estimators get smaller when δ is greater than +0.20 or less than -0.20. We examined RMSE variation with δ for sample sizes ranging from N=8 to N=64 at different SNR values to get a whole picture of the behavior of the estimators. It has been observed that when the sample size is large, the RMS errors of the sinc estimators become lower than those of the Jacobsen and Candan estimators when δ is greater than +0.15 and less than -0.15. When the sample size is reduced to as small as N=8, as in Figure 10, the RMS errors of the sinc estimators go below those of the others when δ is greater than +0.20 or less than -0.20.Figure 10
RMS errors of the estimators with respect toδ for N=8 and SNR=8 dB.Figure11 shows how the bias of the estimators changes with respect to SNR for N=8 and δ=0.3. The sinc estimators δ^1 and δ^3 have lower biases than those of the Jacobsen estimator, while their biases are almost the same as those of the Candan estimator.Figure 11
Biases of the estimators vs. SNR forN=8 and δ=0.3.Figure12 shows how the RMS errors of the estimators’ changes vs. number of periods for N=32, while the SNR level was kept at 10 dB. The sinc estimators δ^1 and δ^3 have lower RMS errors than the Jacobsen and its bias-corrected version estimator for all spectrum ranges. The number of periods does not significantly affect the RMS errors of the estimators.Figure 12
RMS error vs. number of periods.


### 4.2. Simulation Results for Real Signals

In this section, numerical simulations have been carried out in order to determine how the proposed estimators behave in real signals. A real signal with an amplitude of 1 V, a randomly changing phase between0 and 2π, and a frequency of f=kpeak/Nhave been used.The signalxn(38)xn=Vm.cos2πkpeak/Nn+wnA,n=0,⋯,N−1,has also been contaminated with a white Gaussian noise, wn which has zero mean and variance σw2. The signal-to-noise ratio is defined as SNR=Vm2/σw2.Figures13 and 14 examine how the estimators behave in real signals for different sample sizes N and SNR levels for the signal frequency at bin ten with an offset of δ=0.35. Both figures show that δ^1 and δ^3 estimators give estimation errors of 0.0196 and 0.011 bin widths for SNR=20dB, N=32, and N=64, respectively. These RMS errors indicate that δ^1 and δ^3 estimators developed for complex signals perform as well as for real signals.Figure 13
Comparison RMS errors of the estimators vs. SNR forN=16 and δ=0.35, in case of real signal.Figure 14
Comparison RMS errors of the estimators vs. SNR forN=64 and δ=0.35, in case of real signal.


## 4.1. Simulation Results for Complex Signals

During numerical simulations, a complex signal with an amplitude of 1 V, a randomly changing phase between0 and 2π, and a frequency f=kpeak/N have been used. The signal has also been contaminated with a white Gaussian noise
(37)xn=Vm.ej2πkpeak/Nn+θ+wn,n=0,⋯,N−1.Herewn is a circularly symmetric complex white Gaussian noise with mean zero and variance σw2. The signal-to-noise ratio is defined as SNR=Vm2/σw2.Figures7 and 8 examine how sample size N and SNR levels affect estimators’ behavior for the signal frequency at bin ten with an offset of δ=0.35. Both figures show that the RMSE decreases as the sample size increases. The first and third sinc estimators outperform the other estimators at all levels of SNR values.Figure 7
Comparison of RMS errors of the estimators forN=64 and δ=0.35.Figure 8
Comparison of RMS errors of the estimators forN=32 and δ=0.35.Figure9 compares the sinc estimators with the estimators of Jacobsen and Candan. The number of samples was reduced to 8, but the signal frequency was kept the same as in Figures 7 and 8. Figure 9 shows that the sinc estimators perform well even at a small sample size of N=8. Also, RMSE improvements are so distinct at low SNR values.Figure 9
RMS errors of the estimators vs. SNR forN=8 and δ=0.35.Table2 shows how much execution time is needed for each estimator at bin frequency 10 for a sample size of N=32. Execution times have been measured in a loop and then averaged to find the time for a single run. The price to be paid for improvements in RMSE is roughly a microsecond.Table 2
The proposed frequency estimation method for small sample values.
Step 1(i) Search the peak sample from the N-point DFT and its left and right neighborsXkp,Xkp−1, and Xkp+1Step 2(i) Ifδ calculated from (1) is δ≥0.2, recalculate the correction factor from the following equation:δ^3≅ReNπ.Xkp+1.sinπ/N−Xkp+Xkp+1.cosπ/NStep 3(i) Ifδ calculated from (1) is δ≤−0.2, recalculate the correction factor from the following equation:δ^3≅ReNπ.−Xkp−1.sinπ/N−Xkp+Xkp−1.cosπ/NStep 4(i) Calculatekpeak and f as follows:kpeak=kp+δ and f=kpeak/N, for normalized frequencyFigure10 expresses the RMS errors as a function of the correction factor δ when the sample size N=8 and SNR=8 dB. It shows that as the frequency, we are trying to estimate is at a location toward the bin edges; the RMS errors of the sinc estimators get smaller when δ is greater than +0.20 or less than -0.20. We examined RMSE variation with δ for sample sizes ranging from N=8 to N=64 at different SNR values to get a whole picture of the behavior of the estimators. It has been observed that when the sample size is large, the RMS errors of the sinc estimators become lower than those of the Jacobsen and Candan estimators when δ is greater than +0.15 and less than -0.15. When the sample size is reduced to as small as N=8, as in Figure 10, the RMS errors of the sinc estimators go below those of the others when δ is greater than +0.20 or less than -0.20.Figure 10
RMS errors of the estimators with respect toδ for N=8 and SNR=8 dB.Figure11 shows how the bias of the estimators changes with respect to SNR for N=8 and δ=0.3. The sinc estimators δ^1 and δ^3 have lower biases than those of the Jacobsen estimator, while their biases are almost the same as those of the Candan estimator.Figure 11
Biases of the estimators vs. SNR forN=8 and δ=0.3.Figure12 shows how the RMS errors of the estimators’ changes vs. number of periods for N=32, while the SNR level was kept at 10 dB. The sinc estimators δ^1 and δ^3 have lower RMS errors than the Jacobsen and its bias-corrected version estimator for all spectrum ranges. The number of periods does not significantly affect the RMS errors of the estimators.Figure 12
RMS error vs. number of periods.


## 4.2. Simulation Results for Real Signals

In this section, numerical simulations have been carried out in order to determine how the proposed estimators behave in real signals. A real signal with an amplitude of 1 V, a randomly changing phase between0 and 2π, and a frequency of f=kpeak/Nhave been used.The signalxn(38)xn=Vm.cos2πkpeak/Nn+wnA,n=0,⋯,N−1,has also been contaminated with a white Gaussian noise, wn which has zero mean and variance σw2. The signal-to-noise ratio is defined as SNR=Vm2/σw2.Figures13 and 14 examine how the estimators behave in real signals for different sample sizes N and SNR levels for the signal frequency at bin ten with an offset of δ=0.35. Both figures show that δ^1 and δ^3 estimators give estimation errors of 0.0196 and 0.011 bin widths for SNR=20dB, N=32, and N=64, respectively. These RMS errors indicate that δ^1 and δ^3 estimators developed for complex signals perform as well as for real signals.Figure 13
Comparison RMS errors of the estimators vs. SNR forN=16 and δ=0.35, in case of real signal.Figure 14
Comparison RMS errors of the estimators vs. SNR forN=64 and δ=0.35, in case of real signal.


## 5. Conclusions

The estimator proposed by Jacobsen has drawn wide attention because of its fast and accurate computation. There have been contributions regarding its further improvements. This study is another attempt to increase its accuracy.All the simulations performed showed that the improvements suggested by the sinc estimators,δ^1 and δ^3, are remarkable even for small sample sizes and in a wide range of SNR values. When the sinc estimators δ^1 and δ^3 are compared with Jacobsen and Candan estimators for SNR=20 dB, N=8, and an offset of δ=0.45, estimation errors are roughly 0.0150, 0.0327, and 0.0285 bin widths for the sinc, Jacobsen, and Candan estimators, respectively. So, we propose a frequency estimation method for the Jacobsen estimator (1), as stated in Table 3. The price paid for the reduced RMS errors is a slight increase in the computational workload.Table 3
Execution times of the estimators.
EstimatorsExecution timeJacobsen estimator (1)2.2 microsecondsJacobsen bias-corrected (2)3.2 microsecondsδ^1 (16)1.1 microsecondsδ^3 (31)1.8 microseconds



---

*Source: 5200230-2022-11-28.xml*