Exact Frequency Formula for a Pure Real Tone in a DFT

• Basics
• DFT_FFT
• Tutorials

This is an article to hopefully give a better understanding of the Discrete Fourier Transform (DFT) by deriving an exact formula for the frequency of a real tone in a DFT. According to current teaching, this is not possible, so this article should be considered a major theoretical advance in the discipline. The formula is presented in a few different formats. Some sample calculations are provided to give a numerical demonstration of the formula in use. This article is an extension of my previous blog article:

• DFT Bin Value Formulas for Pure Real Tones

You should really read it first before reading this article.

In my previous blog article "DFT Bin Value Formulas for Pure Real Tones", closed form analytic formulas were derived for DFT bin values for a pure tone. For any given bin and frame size, the value of the bin is strictly a function of the three parameters that define the signal. Given three different bin values, a system of three equations/three unknowns is formed. Fortunately, the equations are simple enough that this system can be solved analytically. This blog article shows the solution of the frequency value. The magnitude and phase can also be solved for. It turns out that the frequency solution is robust in the presence of noise and other signals whereas the other two are extremely sensitive. Therefore, the magnitude and phase solutions are not presented and an alternative approach, which is also robust, will be presented in a subsequent blog article.

A DFT is done on finite length signal consisting of $ N $ evenly space points in time. This is called a frame. A pure real tone consists of a signal that is sinusoidal. This can be modeled as the projection onto the real axis of a point travelling around the complex unit circle at a constant speed, then magnified by a real scalar factor. A mathematical equation for this is:

$$ S_n = M cos( f \cdot { 2\pi \over N } n + \phi ) \tag {1} $$

The number of times the point travels around the circle, called cycles, is called the frequency ($f$). The units of $f$ are therefore cycles per frame. The angle to the starting point is $\phi$ and this is measured in radians, which is also the distance along the circumference. The magnification factor ($M$) rescales the whole signal and is called the amplitude. The extremes of the cosine function are -1 and 1, so (M) defines the minimum and maximum possible values for $ S_n $. Because of discrete sampling, these extremes may or may not be reached in a particular frame. The value of $ f $ has to be rescaled to radians so the argument to the cosine function is the position on the unit circle of the travelling point at each $ n $th sampling.

These definitions come from the previous blog article.

The distance travelled by the point along the circumference of the complex unit circle from one sample to the next is $ \alpha $.

$$ \alpha = f \cdot { 2\pi \over N } \tag {2} $$

When the unit circle is divided into N equal pieces, the distance to each piece is given by the $ \beta_k $ values. They are indexed counter-clockwise on the circle.

$$ \beta_k = k \cdot { 2\pi \over N } \tag {3} $$

The points between the pieces are located at a set of Nth Roots of Unity. They are indexed clockwise on the circle.

$$ R_k = e^{ -i \beta_k } \tag {4} $$

The equation for the DFT bin values for a pure real tone can be stated like this:

$$ Z_k = { M \over { 2N } } \left[ { Ue^{ i\beta_k } - V \over cos \alpha - cos \beta_k} \right] \tag {5} $$

Where $ U $ and $ V $ are defined expressions to simplify the bin value equation.

$$ U = cos( \alpha N + \phi ) - cos(\phi) \tag {6} $$

$$ V = cos( \alpha N -\alpha + \phi ) - cos( -\alpha + \phi ) \tag {7} $$

There were no approximations taken during the derivation of these formulas. They are exact.

Solving for the frequency means solving for $ \alpha $ and then converting the value. Since the equation has $ \alpha $ as the argument of a cosine function, $ cos \alpha $ will actually be solved for and $ \alpha $ derived from there.

The bin value equation for bin $ k $ (5) can be cross multiplied to yield:

$$ 2N \left[ cos \alpha - cos \beta_k \right] Z_k = M Ue^{ i\beta_k } - M V \tag {8} $$

This same equation can also be expressed for bins $ k-1 $ and $ k+1 $. Actually any other two bins could be selected, but it makes sense to use three adjacent bins. This is because the DFT concentrates the information about a tone signal in the bins nearest to the frequency. Of course, for integer frequencies it all goes into one bin and a formula is unnecessary.

$$ 2N \left[ cos \alpha - cos \beta_{k-1} \right] Z_{k-1} = M Ue^{ i\beta_{k-1} } - M V \tag {9} $$

$$ 2N \left[ cos \alpha - cos \beta_{k+1} \right] Z_{k+1} = M Ue^{ i\beta_{k+1} } - M V \tag {10} $$

In order to do the next steps the exponential expressions on the right hand side need to be manipulated slightly.  A factor of $ e^{ -i\beta_k } $ is what distinguishes the U terms in the equations.  This also happens to be the first Root of Unity ($ R_1 $) which can be seen from (3) and (4).

$$ e^{ i\beta_{k-1} } = e^{ i\beta_k - i\beta_1 } = e^{ i\beta_k } \cdot e^{ -i\beta_1 } = e^{ i\beta_k } \cdot R_1 \tag {11} $$

$$ e^{ i\beta_{k+1} } = e^{ i\beta_k + i\beta_1 } = e^{ i\beta_k } / e^{ -i\beta_1 } = e^{ i\beta_k } / R_1 \tag {12} $$

With these changes applied, (9) and (10) can be subtracted from (8) to produce these two equations:

$$ 2N \left[ cos \alpha( Z_k - Z_{k-1} ) - ( cos \beta_k Z_k - cos \beta_{k-1} Z_{k-1} ) \right] \\ \quad = M Ue^{ i\beta_k }( 1 - R_1 ) \tag {13} $$

$$ 2N \left[ cos \alpha( Z_k - Z_{k+1} ) - ( cos \beta_k Z_k - cos \beta_{k+1} Z_{k+1} ) \right] \\ \quad = M Ue^{ i\beta_k }( 1 - 1/R_1 ) \tag {14} $$

By making these subtractions, the unknown value $ V $ as been eliminated.

Again, making a slight substitution will be the next step a little easier to see.

$$ 1 - 1/R_1 = ( R_1 - 1 ) / R_1 \tag {15} $$

The next step is to divide (13) by (14).

$$ { cos \alpha( Z_k - Z_{k-1} ) - ( cos \beta_k Z_k - cos \beta_{k-1} Z_{k-1} ) \over cos \alpha( Z_k - Z_{k+1} ) - ( cos \beta_k Z_k - cos \beta_{k+1} Z_{k+1} ) } = -R_1 \tag {16} $$

This one step removes two unknowns $ M $ and $ U $. It also almost eliminates $ N $. Only $ \alpha $ remains as an unknown. The value of $ cos \alpha $ appears twice in the equation and there are no other occurrences of $ \alpha $. By a little straightforward algebra, $ cos \alpha $ can be solved for.

$$ cos \alpha = \\ { -cos \beta_{k-1} Z_{k-1} + (1+R_1)cos \beta_k Z_k - R_1 cos \beta_{k+1} Z_{k+1} \over -Z_{k-1} + (1+R_1)Z_k - R_1 Z_{k+1} } \tag {17} $$

The goal has now been reached. All the terms on the right hand side are known and therefore $ cos \alpha $ is now known.

There is no requirement the $ k $ that is selected is near $ f $, unless $ f $ is an integer. In that case, one of the $ Z $s has to be non-zero. This is the degenerate case and will still work. This can be seen by inspection from (17). For non-integer frequencies, any three bin set can be chosen for a pure real tone DFT bin set. When noise or other tones are present, $ k $ should be chosen to the the closest to $ f $. This will be the bin with the greatest magnitude.

With a little rearrangement, (16) can be put in a bilinear matrix form. Indeed, by putting (8)-(10) in matrix form, the bilinear form can be solved for directly.

$$ cos \alpha = \\{ { \begin{bmatrix} -1 \\ 1+R_1 \\ -R_1 \\ \end{bmatrix} }^T \begin{bmatrix} cos \beta_{k-1} & 0 & 0 \\ 0 & cos \beta_k & 0 \\ 0 & 0 & cos \beta_{k+1} \\ \end{bmatrix} \begin{bmatrix} Z_{k-1} \\ Z_{k} \\ Z_{k+1} \\ \end{bmatrix} \over { \begin{bmatrix} -1 \\ 1+R_1 \\ -R_1 \\ \end{bmatrix} }^T \begin{bmatrix} Z_{k-1} \\ Z_{k} \\ Z_{k+1} \\ \end{bmatrix} } \tag {18} $$

In this format, how the formula actually works is a little clearer.

A different arrangement in matrix form has the advantage of showing the way to calculate the formula in the most efficient manner. The left vector in both the numerator and denominator is the same and needs only to be calculated once.

$$ cos \alpha = { { \begin{bmatrix} -Z_{k-1} \\ (1+R_1)Z_{k} \\ -R_1 Z_{k+1} \\ \end{bmatrix} }^T \begin{bmatrix} cos \beta_{k-1} \\ cos \beta_k \\ cos \beta_{k+1} \\ \end{bmatrix} \over { \begin{bmatrix} -Z_{k-1} \\ (1+R_1)Z_{k} \\ -R_1 Z_{k+1} \\ \end{bmatrix} }^T \begin{bmatrix} 1 \\ 1 \\ 1 \\ \end{bmatrix} } \tag {19} $$

The denominator is a scalar. It could actually be divided into every element of the left vector in the numerator to clearly show that $ cos \alpha $ is a weighted average of $ cos \beta_{k-1} $, $ cos \beta_k $, and $ cos \beta_{k+1} $.

Once the value of $ cos \alpha $ has been found, calculating the frequency ($ f $) is straightforward.

$$ f = \alpha \cdot { N \over 2\pi } = cos^{-1} \left( cos \alpha \right) \cdot { N \over 2\pi } \tag {20} $$

Please note that inverse cosine function can return multiple values. The assumption is that $ \alpha $ is between 0 and $ \pi $, making $ f $ between 0 and $ N/2 $. The possible solutions outside this range are called aliases. The fact that the solution can generate the possible aliases is a testament to its correctness.

In the presence of noise or other tones, the value of $ cos( \alpha ) $ may have an imaginary component.  This can be discarded by using just the real component for the $ cos^{-1} $ function.  In practice, it is also a good idea to range check the value to be between -1 and 1 inclusive.

The denominator in (18) deserves some special mention. The left vector appears both in the numerator and the denominator, so it can be rescaled arbitrarily. When $ N $ is large, the value of $ R_1 $ is approximately one.

$$ R_1 \approx 1 \tag {21} $$

As $ N $ gets larger, this becomes even more true.

$$ \lim \limits_{N \to \infty} R_1 = 1 \tag {22} $$

This same limit argument can be applied to the left vector:

$$ \lim \limits_{N \to \infty} \begin{bmatrix} -1 & 1+R_1 & -R_1  \end{bmatrix} \\ = \begin{bmatrix} -1 & 2 & -1  \end{bmatrix} = 4 \cdot \begin{bmatrix} -.25 & .50 & -.25  \end{bmatrix} \tag {23} $$

The rightmost vector in (23) is the same as the weighting associated with the Von Hann window function. The denominator in (18) is then approximately the value one would get in bin $ k $ by applying the Von Hann window function. Since the DFT bin values are needed in isolation in the numerator, there is no point in applying the window function. However, since the weighted multiplication does take place, it is as if the window function has been applied. The ramification of this is that the portion of the signal that is in the center of the frame is much more significant than the ends.

The formula in (19) is fairly easy to compute. Three examples are provided. The sample DFT values are from the previous blog article. The equation for the signal is: $$ S_n = 1.0 \cdot cos( 10.4 \cdot { 2\pi \over N } n + .6 ) \tag {24} $$

There are 32 ($ N $) points in the sequence and therefore that many DFT bins.

The value of $ { \begin{bmatrix} -1 & 1+R_1 & -R_1 \end{bmatrix} }^T $ is :

   -1.00000000000   0.00000000000 
    1.98078528040  -0.19509032202 
   -0.98078528040   0.19509032202

In the three calculations, the bin values and the cosine values are listed in the first block. The second block shows the evaluation of (19). The first two columns are the denominator's real and imaginary parts, the last two columns are the numerator's real and imaginary parts. The first three rows are the three elements of the vectors. The dashed line separates the sums, so the next row is the denominator and numerator as complex values. The last row has the quotient as a complex number followed by the frequency.

This first example is the DFT evaluated at the bin with the biggest magnitude. Notice that the frequency is exactly correct to the precision of the display. Since the bin is near the top of the unit circle ( $ N/4 = 8 $ ), the cosine values are nearly linear.

  9  -0.00032563186   0.10802118551      -0.19509032202 
 10  -0.07619790924   0.36944527683      -0.38268343237 
 11   0.10202082457  -0.23340312262      -0.55557023302
 
  0.00032563186  -0.10802118551    -0.00006352762   0.02107388787
 -0.07885649900   0.74665724091     0.03017707570  -0.28573335575
 -0.05452583268   0.24882162258     0.03029292957  -0.13823788684
---------------  --------------    --------------  --------------
 -0.13305669982   0.88745767798     0.06040647764  -0.40289735473

 -0.45399049974   0.00000000000            10.40000000000

This second example is an evaluation at the Nyquist bin ($ \beta_n = \pi $). This is six bins away from the peak bin. Since the values of the DFT drop off from the peak bin some error in the valuation of the frequency is introduced due to the truncation of the input values.

 15   0.04268851510  -0.01055994389      -0.98078528040 
 16   0.04218971842   0.00000000000      -1.00000000000 
 17   0.04268851510   0.01055994389      -0.98078528040 

 -0.04268851510   0.01055994389     0.04186826725  -0.01035703753
  0.08356877323  -0.00823080575    -0.08356877323   0.00823080575
 -0.04392841011  -0.00202892137     0.04308433802   0.00198993622
---------------  --------------    --------------  --------------
 -0.00304815198   0.00030021677     0.00138383205  -0.00013629556

 -0.45399050196  -0.00000000000            10.40000001267

The third and last example is evaluated at the DC bin. This is ten bins away from the peak bin. The DFT values fall off even more so there is a touch more inaccuracy in the results. Six zeroes still means a million.

 31   0.02331048640  -0.00387720744       0.98078528040 
  0   0.02337925966   0.00000000000       1.00000000000 
  1   0.02331048640   0.00387720744       0.98078528040 

 -0.02331048640   0.00387720744    -0.02286258194   0.00380270799
  0.04630929340  -0.00456106730     0.04630929340  -0.00456106730
 -0.02361898759   0.00074494231    -0.02316515536   0.00073062845
---------------  --------------    --------------  --------------
 -0.00062018059   0.00006108246     0.00028155610  -0.00002773086

 -0.45399050301   0.00000000000            10.40000001872

Bins that are further away from the peak are going to be more susceptible to noise. Therefore it is advisable to do the evaluation at the peak bin. If the noise level is high enough, the advantages of an exact formula may not be worth the extra calculations compared to a quick and dirty estimation formula.

An exact formula for the frequency of a pure real tone in a DFT exists.  Techniques for solving for the phase and magnitude, handling multiple tones, and noise reduction will be addressed in upcoming blog articles.

Afterword

http://www.tsdconseil.fr/log/scriptscilab/festim/index-en.html

Thanks to Julien Arzi for such an excellent analysis.

[Edit 2015-06-01: Swap n and k in subscripts, suggestions from Steve Pope (Thanks)]

Sign in