# Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT

April 13, 2017

## Introduction

This is an article to hopefully give a better understanding to the Discrete Fourier Transform (DFT) by deriving exact formulas for the frequency of a complex tone in a DFT. This time it is three bin versions. Although the problem is similar to the two bin version in my previous blog article "A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT"[1], a slightly different approach is taken using linear algebra concepts. Because of an extra degree of freedom gained by including an additional bin value, many formulas are possible. This article derives and evaluates a few of them.

## A Pure Complex Tone Sequence

The same variables and equation for describing the tone from my previous articles [1] [2] will be used. The variables are amplitude ($M$), frequency ($f$), phase ($\phi$), and number of samples ($N$). The equation is: $$S_n = M e^{i( f \cdot \frac{ 2\pi }{ N } n + \phi )} \tag {1}$$ Also the same angle variables: $$\alpha = f \cdot \frac{ 2\pi }{ N } \tag {2}$$ $$\beta_k = k \cdot \frac{ 2\pi }{ N } \tag {3}$$ Finally, the same simpler formula for the signal: $$S_n = M \cdot { e^{ i(\alpha n + \phi) } } \tag {4}$$

## Bin Value Formula for a Pure Complex Tone

This is a formula to calculate the bin values in a DFT for a pure complex tone: $$Z_k = \frac{ M }{ N } e^{ i \phi } \cdot \frac{ 1 - e^{ i\alpha N } }{ 1 - e^{ i(\alpha - \beta_k ) } } \tag {5}$$ The derivation of this formula can be found in my previous blog article titled "DFT Bin Value Formulas for Pure Complex Tones." [2]

## Solving for the Frequency Term

Rather than solve for an equation using three arbitrary bins, which can be done, it will be more efficient to solve for a formula using three adjacent bins. This is the formula that will be needed for practical use. Assume that $Z_k$ is the peak bin value with $Z_{k-1}$ and $Z_{k+1}$ as the other two bin values. The formula will work on any three bin set of a DFT for a non-integer frequency pure complex tone, but it will be assumed the three bins are centered on the peak.

Also, rather than solving for the frequency term as the two bin version did [1], this version will solve directly for the distance of the frequency from the peak bin. Two special variable substitutions are introduced to greatly simplify the equations in the derivation. $$a = e^{ i(\alpha - \beta_k ) } \tag {6}$$ $$r = e^{ -i\beta_1 } \tag {7}$$ The first one, which represents the frequency term to be solved for, comes from my original derivation of this formula. The second one, also from the original derivation, shows up in the end product found in Julien Arzi's comparison paper[6]. In the real tone version of this article [3], the variable $r$ is referred to as $R_1$, the first Root of Unity. The simpler $r$ will make the formulas in this article a little more comprehensible.

From these, the following equations are true: $$e^{ i(\alpha - \beta_{k-1} ) } = a / r \tag {8}$$ $$e^{ i(\alpha - \beta_{k+1} ) } = a \cdot r \tag {9}$$ Multiplying (5) by the denominator on the right side, and applying the previous substitutions to all three bin equations yields this set: $$Z_{k-1} - \frac{a}{r} Z_{k-1} = \frac{ M }{ N } e^{ i \phi } \cdot \left( 1 - e^{ i\alpha N } \right) \tag {10}$$ $$Z_k - a Z_k = \frac{ M }{ N } e^{ i \phi } \cdot \left( 1 - e^{ i\alpha N } \right) \tag {11}$$ $$Z_{k+1} - a r Z_{k+1} = \frac{ M }{ N } e^{ i \phi } \cdot \left( 1 - e^{ i\alpha N } \right) \tag {12}$$ The three bin set of the DFT can be thought of as a vector. $$\vec Z = \begin{bmatrix} Z_{k-1} \\ Z_k \\ Z_{k+1} \end{bmatrix} \tag {13}$$ This allows the three equation set to be combined into one matrix equation. $$\vec Z - a \begin{bmatrix} \frac{1}{r} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & r \end{bmatrix} \vec Z = \frac{ M }{ N } e^{ i \phi } \cdot \left( 1 - e^{ i\alpha N } \right) \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \tag {14}$$ The next step is to find a vector $\vec K$ that is perpendicular to $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$. There are an infinite number to choose from. All that matters now is that at least one such vector exists. The diagonal matrix on the left hand side can be called $D$.

Dotting (14) with $\vec K$ collapses it down to a scalar equation. $$\vec K \vec Z - a \cdot \vec K D \vec Z = 0 \tag {15}$$ The right hand side becomes zero because $\vec K \cdot \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} = 0$. Solving for $a$ is now trivial. $$a = \frac{ \vec K \vec Z }{ \vec K D \vec Z } \tag {16}$$

## Calculating the Frequency

Once a $\vec K$ is selected, $a$ can be easily calculated. From there it is just a few steps to calculate the frequency ($f$). $$i(\alpha - \beta_k) = \ln(a) \tag {17}$$ The defintion used for the natural log of a complex number is shown in my previous article[1]. $$\alpha - \beta_k = \operatorname{ atan2 }( Im[a], Re[a] ) \tag {18}$$ Converting back to a bin scale and bringing the base bin index back to the right hand side produces the final version. $$f = k + \operatorname{ atan2 }( Im[a], Re[a] ) \cdot \frac{ N }{ 2\pi } \tag {19}$$ This approach is the same as the one used in the two bin article.

## Possible Values of $\vec K$

There are a lot of candidates for $\vec K$ to choose from. The main constraint is that $\vec K \cdot \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} = 0$. Since $\vec K$ is a three dimensional vector, this constraint cuts down its degrees of freedom to two. $\vec K$ has to lie in the plane that is perpendicular to $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$. The vector appears in both the numerator and the denominator in (16) so its length doesn't matter (as long as it's non-zero). This reduces its degree of freedom by one more leaving one. The entire set of possible vectors is then a one complex valued parameter family. It can be parameterized in multiple ways, or none at all, which is how it will be done in this article.

Any choice of $\vec K$ within the family will work in the noiseless case as long as $\vec K \cdot \vec Z \ne 0$ and $\vec K \cdot D \cdot \vec Z \ne 0$. Therefore the criteria for the best value of $\vec K$ is how robust it is in the presence of noise. The effect of other tones is not as important because their influence can be removed by other means which I will demonstrate in future blog articles.

## $\vec K$s with Zeroes

A quick inspection will confirm that $\begin{bmatrix} 0 & 1 & -1 \end{bmatrix}$ for $\vec K$ gives the two bin formula from my previous blog article[1]. Switching it around to $\begin{bmatrix} -1 & 1 & 0 \end{bmatrix}$ gives the corresponding two bin formula for the $k-1$ and $k$ bins. The choice of $\begin{bmatrix} -1 & 0 & 1 \end{bmatrix}$ is clearly not going to be as good since it ignores the peak bin which has the most signal information. Any two of these three vectors can be chosen as a basis set to span the entire set of possible $\vec K$s.

## The Von Hann $\vec K$

Another obvious candidate is the Von Hann kernel of $\begin{bmatrix} -.25 & .5 & -.25 \end{bmatrix}$. Since length doesn't matter it can be conveniently rescaled to $\begin{bmatrix} -1 & 2 & -1 \end{bmatrix}$. Perhaps it should be called Pascal's instead since this is an alternating row of Pascal's Triangle. The vector is, of course, simply the sum of both of the two bin versions $\begin{bmatrix} -1 & 1 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 1 & -1 \end{bmatrix}$ which is how Pascal's Triangle is constructed.

The nice thing about the Von Hann vector is that not only does it take out constant trends $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$, which is the criteria for being a valid $\vec K$, but it also zeroes out linear trends $\begin{bmatrix} m-d & m & m + d \end{bmatrix}$. This means it should be effective in taking out the effect of nearby tones as well as reducing random noise by subtracting out the best fit line to the noise.

## The $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ Squasher

When $\vec Z$ gets multiplied by $D$ in the denominator of (16), the $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ component of $\vec Z$ no longer gets zeroed by $\vec K$. There is a slight "twisting" that happens to $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ due to the multiplication. $$D \cdot \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} = \begin{bmatrix} \frac{1}{r} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & r \end{bmatrix} \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} = \begin{bmatrix} \frac{1}{r} \\ 1 \\ r \end{bmatrix} \tag {20}$$ This "twisted" vector is going to be linearly independent of $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ so the cross product can be used to find a vector that is perpendicular to both of them. $$\begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \times \begin{bmatrix} \frac{1}{r} & 1 & r \end{bmatrix} = \begin{bmatrix} r-1 & -r + \frac{1}{r} & 1 - \frac{1}{r} \end{bmatrix} \tag {21}$$ If this vector is used for $\vec K$, all vestiges of the $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ from the $\vec Z$ will be zeroed out. Since $\vec K$ can be rescaled arbitrarily, it can be multiplied by $-r$ to get rid of the fractions and reverse the signs. $$\begin{bmatrix} r-1 & -r + \frac{1}{r} & 1 - \frac{1}{r} \end{bmatrix} \cdot (-r) = \begin{bmatrix} -r^2+r & r^2 - 1 & -r + 1 \end{bmatrix} \tag {22}$$ In this new form it is easy to see that a factor of $r-1$ can be taken out. $$\frac{\begin{bmatrix} -r^2+r & r^2 - 1 & -r + 1 \end{bmatrix} }{(r-1)} = \begin{bmatrix} -r & r + 1 & -1 \end{bmatrix} \tag {23}$$ This leaves the simplest form of this $\vec K$ selection. One thing to notice is that for larger sample sizes ($N$) the value of $r$ gets closer to 1 and the value of this version of $\vec K$ gets closer to the Von Hann one.

Since $\vec K$ and $D$ are known ahead of time, their product can be precomputed. $$\begin{bmatrix} -r & r + 1 & -1 \end{bmatrix} \begin{bmatrix} \frac{1}{r} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & r \end{bmatrix} = \begin{bmatrix} -1 & r + 1 & -r \end{bmatrix} \tag {24}$$ An interesting thing happens in this operation, the result is the reverse of the $\vec K$. A neat little fact if you are trying to memorize the formula. It is also interesting that this is the same vector that appears in my version of an exact frequency formula for a real tone[3].

Plugging this result into (16) gives the following equation: $$a = \frac{ \begin{bmatrix} -r & r + 1 & -1 \end{bmatrix} \vec Z } { \begin{bmatrix} -1 & r + 1 & -r \end{bmatrix} \vec Z } \tag {25}$$ The inverse tangent function can be used instead of the atan2 function in (19) because the argument is expected to be less than a bin width on the angle scale. $$f = k + \tan^{-1} \left( \frac{Im[a]}{Re[a]} \right) \cdot \frac{ N }{ 2\pi } \tag {26}$$ Compare these results with the following two equations: $$q = \frac{ \begin{bmatrix} 1 & 0 & -1 \end{bmatrix} \vec Z } { \begin{bmatrix} -1 & 2 & -1 \end{bmatrix} \vec Z } \tag {27}$$ $$f = k + \tan^{-1} \left( Re[q] \tan \left(\frac{ \pi }{ N } \right) \right) \cdot \frac{ N }{ \pi } \tag {28}$$ The first (27) is the delta term for Jacobsen's Estimator [4]. The second (28) is Candan's 2013 [5] frequency formula based on Jacobsen's formula. The reason it is included here is because it happens to be mathematically equivalent to (26). This is not obvious by any means just looking at them. It takes several pages of algebra to prove it.

A performance analysis of all these formulas can be found in Julien Arzi's comparison paper[6].

## The Projection of $\vec Z$ and $D \vec Z$

The previous candidates for $\vec K$ are all independent of $\vec Z$. That is to say, they are all predetermined. This method attempts to select the best $\vec K$ based on the value of $\vec Z$.

First, assume that the DFT bins you are working with are the true value of $\vec Z$ plus some sort of error. If the error is due to white noise in the signal then the expected value of $\vec E$ will be zero, but there is a variance to this expectation. In other words, it might be the expectation that $\vec E$ is zero, but don't expect it to be zero. $$a = \frac{ \vec K ( \vec Z + \vec E ) }{ \vec K D ( \vec Z + \vec E ) } \tag {29}$$ $$a = \frac{ \vec K \vec Z + \vec K \vec E }{ \vec K D \vec Z + \vec K D \vec E } \tag {30}$$ If the noise is truly white, or nothing is known about it, then $\vec E$ can go in any direction and cannot be used to help decide what the best $\vec K$ should be.

Even though it is not possible to minimize the contribution of the error to the quotient, it is possible to try to maximize the contribution of the signal information by maximizing the magnitudes of $\vec K \vec Z$ and $\vec K D \vec Z$. It would be nice to maximize both of them, but that can't be done. To maximize either of the values, simply find the projection of the conjugate of the vector onto the plane that is orthogonal to $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$. Maximizing one over the other is not likely to lead to optimal results so finding a good compromise is necessary. A good guess is to use the average of the two. Since the length ends up being irrelevant, a little bit of computation can be saved by using the sum instead. What is important is the direction of the vector. $$\vec G = \overline{\vec Z + D \vec Z} \tag {31}$$ The projection onto the plane is going to be the vector minus some multiple of the normal to the plane. The projection is what will be used for $\vec K$. $$\vec K = \vec G - s \cdot \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \tag {32}$$ To find the multiplier ($s$), the equation needs to be dotted with the normal to the plane. $$\vec K \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} = \vec G \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} - s \cdot \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} = 0 \tag {33}$$ The multiplier can now be solved for. $$s = \frac{1}{3} \cdot \vec G \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \tag {34}$$ Substituting this value back into (32) gives the final result. $$\vec K = \vec G - \frac{1}{3} \left( \vec G \cdot \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \right) \begin{bmatrix} 1 & 1 & 1 \end{bmatrix} \tag {35}$$ This equation translates into an easy to understand computation: Take the average of the three values in $\vec G$ and subtract it from each value.

In the case of an integer frequency tone, the $\vec Z$ will be in the form of $\begin{bmatrix} 0 & Z_k & 0 \end{bmatrix}$. When it is run through the projection calculation the result is $\frac{2}{3} \overline{Z_k}\begin{bmatrix} -1 & 2 & -1 \end{bmatrix}$ which is equivalent to the Von Hann approach.

## Some Numerical Examples

Unlike traditional estimators, such as Jacobsen's, that perform better with larger sample sizes, all the $\vec K$ selections will give exact answers in the absence of noise. What happens in the presence of noise is what differentiates them. The differences are dependent on the sample size as the following sets of numbers show.

These tests consist of taking a sweep of frequencies from the low end of a bin to the high end. For each frequency a complex tone is generated. The phase is fixed at zero. Unlike the real tone case, a variable phase wouldn't change the results and the same is true for which bin number is tested. The numbers that are given are the standard deviation of the distribution of frequency errors. Large runs make the results a lot more even and close to what theoretical values would be.

The Target Noise Level is the standard deviation of the added noise. A level of 0.355 corresponds to signal to noise ratio (SNR) of 9 dB. The error levels are shown multiplied by 100 for better presentation.

The first data set is for a very small sample size.

The sample count is 10
and the run size is 100000

Errors are shown at 100x actual value

Target Noise Level = 0.355

Low      High      Von      Squash   Project  Candan
Freq  Two Bin  Two Bin   Hann                       2013
----  -------  --------  -------  -------  -------  -------
3.50   4.700    33.440    9.163    9.362    5.195    9.348
3.55   4.722    27.709    8.515    8.736    5.152    8.725
3.60   4.833    23.089    7.887    8.131    5.168    8.122
3.65   5.013    19.839    7.448    7.710    5.212    7.703
3.70   5.262    17.104    7.028    7.311    5.294    7.306
3.75   5.561    14.849    6.659    6.945    5.384    6.941
3.80   5.918    13.070    6.348    6.644    5.466    6.641
3.85   6.411    11.637    6.179    6.473    5.603    6.470
3.90   6.926    10.307    5.981    6.282    5.651    6.280
3.95   7.648     9.267    5.894    6.196    5.713    6.193
4.00   8.345     8.364    5.849    6.141    5.714    6.139
4.05   9.285     7.614    5.899    6.203    5.710    6.201
4.10  10.291     6.951    5.977    6.279    5.654    6.276
4.15  11.615     6.419    6.171    6.468    5.601    6.465
4.20  13.055     5.928    6.371    6.666    5.489    6.662
4.25  14.843     5.572    6.660    6.945    5.400    6.941
4.30  17.097     5.236    7.032    7.295    5.288    7.290
4.35  19.606     5.010    7.406    7.671    5.211    7.664
4.40  23.115     4.826    7.916    8.182    5.149    8.174
4.45  27.515     4.734    8.503    8.740    5.156    8.730
4.50  33.488     4.699    9.148    9.373    5.200    9.358


The earlier claim that the Squash $\vec K$ and Candan's formula are mathematically identical seems suspect from this data. The discrepency can be explained by the fact that the two different forms have different sensitivity to the small sample size.

As noted in the comments of my previous post [1], the two bin solution fares way better in the frequency region between two bins. The high end of the Low Two Bin case, and the low end of the High Two Bin case, show that the two bin solution needs to be on the two peak bins to be robust. What is suprising about this data is how well the Von Hann $\vec K$ does compared to the Squash and Candan values. This was unexpected.

Not surprisingly, the Projection $\vec K$ outperformed the Von Hann, Squash, and Candan formulas. However, the Projection formula does require considerably more calculation, so the extra robustness may not be worth it.

Here is a larger sample size.

The sample count is 100
and the run size is 100000

Errors are shown at 100x actual value

Target Noise Level = 0.355

Low      High      Von      Squash   Project  Candan
Freq  Two Bin  Two Bin   Hann                       2013
----  -------  --------  -------  -------  -------  -------
3.50   1.404     9.460    2.770    2.771    1.514    2.771
3.55   1.416     8.010    2.583    2.583    1.495    2.583
3.60   1.448     6.845    2.414    2.415    1.499    2.415
3.65   1.503     5.880    2.260    2.262    1.522    2.262
3.70   1.563     5.116    2.132    2.133    1.546    2.133
3.75   1.663     4.482    2.023    2.023    1.596    2.023
3.80   1.783     3.955    1.942    1.942    1.652    1.942
3.85   1.921     3.504    1.871    1.872    1.697    1.872
3.90   2.099     3.126    1.827    1.827    1.746    1.827
3.95   2.294     2.810    1.792    1.794    1.767    1.794
4.00   2.516     2.528    1.788    1.789    1.782    1.789
4.05   2.806     2.289    1.798    1.799    1.773    1.799
4.10   3.132     2.097    1.826    1.826    1.744    1.826
4.15   3.497     1.931    1.871    1.872    1.705    1.872
4.20   3.929     1.779    1.935    1.936    1.648    1.936
4.25   4.494     1.661    2.030    2.030    1.595    2.030
4.30   5.117     1.566    2.131    2.133    1.547    2.133
4.35   5.912     1.494    2.263    2.264    1.513    2.264
4.40   6.876     1.445    2.424    2.425    1.496    2.425
4.45   8.030     1.413    2.586    2.587    1.492    2.587
4.50   9.519     1.397    2.781    2.782    1.507    2.782


At this size, the Squash and Candan formulas look identical. The $N$ is large enough so that the Von Hann solution is almost the same as well. The Projection formula is still the best of the three bin formulas.

## A Conundrum

Which formula is best to use depends on where the frequency is relative to the nearest bin, but that isn't known until the formula is calculated. There are several ways to go about it. One would be to choose one of the formulas, say the Von Hann or the Projection, and use it universally. Another would be to start with the Von Hann, and if the frequency was sufficiently away from the nearest bin to apply the two bin formula. The reverse could also be done by starting with the two bin and if the frequency is close to an integer use a three bin at that bin.

It is also possible to get a rough estimate of the frequency by the ratio of the magnitude of the peak bin to the next largest neighbor. This is somewhat similar to the two bin formula, but since cutoffs are all that are needed, doing a quick quotient of the squares of the magnitudes will allow for a more computationally efficient calculation. No square roots need to be taken and a real division can be used instead of a complex one.

Here is a chart of the ratio of the squares of the magnitudes of the peak bin to the next bin for a pure complex tone by the fraction of the frequency from the peak bin for various sample sizes. The values are independent of the bin number and the phase of the tone. For a rough estimate, the values across the various sample counts are close enough that no special consideration is needed.

Frac     8       16      32      64     128
----   -----   -----   -----   -----   -----
0.10   77.72   80.17   80.79   80.95   80.99
0.15   30.97   31.82   32.04   32.09   32.11
0.20   15.51   15.88   15.97   15.99   16.00
0.25    8.77    8.94    8.99    9.00    9.00
0.30    5.33    5.42    5.44    5.44    5.44
0.35    3.40    3.44    3.45    3.45    3.45
0.40    2.23    2.24    2.25    2.25    2.25
0.45    1.49    1.49    1.49    1.49    1.49
0.50    1.00    1.00    1.00    1.00    1.00


If the prior bin to the peak is larger than the following one, use this chart in the negative direction.

A rule of thumb for which frequency formulas should be used can be developed by using this chart. If your strategy is to use the Von Hann close to bin values and the two bin formula between bins then the break even point is around 0.17. This means you should use a ratio cutoff around 29 or 30. If the ratio is lower than this, use the two bin formula and if the ratio is higher use the three bin. The break even point for the Projection method is around .3 which translates into a cutoff value of about 5.

## Conclusion

The three bin solution for frequency calculation comes with choices that have tradeoffs. There is no compelling choice other than, ironically, to use a two bin formula when the frequency in the middle of two bins. At the center frequency, the standard deviation is about half that of the three bin choices which makes the variance about one fourth. This is significantly more robust. This article has not proven that the sum used in (31) is the optimal choice. Nor has it tested the use of precomputed projection vectors based on a rough frequency estimate. In the end though, it should give a better idea of how the bin values of a DFT can be used to find the frequency of complex tones in the signal. Two new formulas (Von Hann and Projection) that beat the current standard (Candan 2013) are introduced.

## References

[1] Dawg, Cedron, A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT

[2] Dawg, Cedron, DFT Bin Value Formulas for Pure Complex Tones

[3] Dawg, Cedron, Exact Frequency Formula for a Pure Real Tone in a DFT

[4] Jacobsen, Eric and Kootsookos, Peter, "Fast, Accurate Frequency Estimators", IEEE SIGNAL PROCESSING MAGAZINE, pp. 123-125, MAY 2007

[5] Candan, Cagatay, "A Method For Fine Resolution Frequency Estimation From Three DFT Samples", IEEE SIGNAL PROCESSING LETTERS, VOL. 18, NO. 6, pp. 351-354, JUNE 2011

[6] Arzi, Julien, Comparison of different frequency estimation algorithms

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.