An Alternative Form of the Pure Real Tone DFT Bin Value Formula
Introduction
This is an article to hopefully give a better understanding of the Discrete Fourier Transform (DFT) by deriving alternative exact formulas for the bin values of a real tone in a DFT. The derivation of the source equations can be found in my earlier blog article titled "DFT Bin Value Formulas for Pure Real Tones"[1]. The new form is slighty more complicated and calculation intensive, but it is more computationally accurate in the vicinity of near integer frequencies. This may be important in high precision implementations of my earlier method for finding the amplitude and phase of a pure tone in a DFT described in my blog article "Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1"[2].
The Starting Equations
This is the same real pure tone formula I use in all the blog articles. The frequency term ($\alpha$) is expressed in units of radians per sample which makes the formula as simple as possible. $$ S_n = M \cos( \alpha n + \phi ) \tag {1} $$ The frequency in the formula can readily be converted from units of cycles per frame using this formula: $$ \alpha = f \cdot \frac{ 2\pi }{ N } \tag {2} $$ This is a formula for a bin value of the Discrete Fourier Transform (DFT) for any signal where $N$ is the number of samples in the frame. $$ Z_k = \frac{ 1 }{ N } \sum_{n=0}^{N-1} { S_n e^{ -i\beta_k n } } \tag {3} $$ Similar to the signal definition, the "frequency term" is the bin index in terms of cycles per frame which can be translated into a unit of radians per sample to simplify the DFT definition. $$ \beta_k = k \cdot \frac{ 2\pi }{ N } \tag {4} $$
Pure Real Tone Bin Value Formula
The bin value of a pure real tone can be found without the need to do a summation by plugging (1) into (3) and applying the formula for a summation of a geometric series. The result is the the following formula: $$ Z_k = \frac{ M }{ 2N } \left[ \frac{ U e^{ i\beta_k } - V } { \cos( \alpha ) - \cos( \beta_k ) } \right] \tag {5} $$ Where: $$ U = \cos( \alpha N + \phi ) - \cos(\phi) \tag {6} $$ $$ V = \cos( \alpha N -\alpha + \phi ) - \cos( -\alpha + \phi ) \tag {7} $$ The derivation can be found in my blog article appropriately titled "DFT Bin Value Formulas for Pure Real Tones" [1]. Please note that the formula has a pluggable discontinuity at integer frequencies in terms of cycles per frame.
The Difference of Two Cosines
The formula has three occurrences of the difference of two cosines. There is a handy trigonometric formula that covers this: $$ \cos( \theta_1 ) - \cos( \theta_2 ) = -2 \sin \left( \frac{ \theta_1 + \theta_2 }{2} \right) \sin \left( \frac{ \theta_1 - \theta_2 }{2} \right) \tag {8} $$
Applying the Differences Equations
The first application of the formula is done to the denominator of (5). $$ \cos( \alpha ) - \cos( \beta_k ) = -2 \sin \left( \frac{ \alpha + \beta_k }{2} \right) \sin \left( \frac{ \alpha - \beta_k }{2} \right) \tag {9} $$ The next two applications are done to (6), and (7). $$ U = -2 \sin \left( \frac{ \alpha N }{2} + \phi \right) \sin \left( \frac{ \alpha N }{2} \right) \tag {10} $$ $$ V = -2 \sin \left( \frac{ \alpha N }{2} - \alpha + \phi \right) \sin \left( \frac{ \alpha N }{2} \right) \tag {11} $$ Notice that there is a common factor. These last three equations can be now combined with (5) into a single equation for the bin value. $$ Z_k = \frac{ M }{ 2N } \cdot \frac{ \sin \left( \frac{ \alpha N }{2} \right) }{ \sin \left( \frac{ \alpha - \beta_k }{2} \right) } \cdot \left[ \frac{ \sin \left( \frac{ \alpha N }{2} + \phi \right) e^{ i\beta_k } - \sin \left( \frac{ \alpha N }{2} - \alpha + \phi \right) }{ \sin \left( \frac{ \alpha + \beta_k }{2} \right) } \right] \tag {12} $$ When the frequency ($f$) of the tone is very close to an integer value, the denominators in (5) and (12) become very small at the bin corresponding to the integer frequency. When the denominator is calculated as the difference of two cosines, that is two values that are not near zero, due to the nature of finite precision in computer calculations, the low order bits of the difference are lost. This is not true in (12) since the small difference has been rendered into a product of a very small number and another not near zero value. The numerator of the sine quotient outside the brackets will also be a very small number. Dividing two very small numbers represented in floating point values will retain much more precision.
Changing the Frame of Reference
As a matter of mathematical analysis and not practical use, the equation can be reframed with the frequency represented as the difference of the frequency parameter and the corresponding bin's frequency parameter. $$ \delta = \alpha - \beta_k = ( f - k ) \cdot \frac{ 2\pi }{ N } \tag {13} $$ $$ \alpha = \delta + \beta_k \tag {14} $$ The substitution is first made into the numerator of the sine quotient outside the brackets in (12). $$ \begin{aligned} \sin \left( \frac{ \alpha N }{2} \right) &= \sin \left( \frac{ ( \delta + \beta_k ) N }{2} \right) = \sin \left( \frac{ ( \delta + k \cdot \frac{ 2\pi }{ N } ) N }{2} \right) \\ &= \sin \left( \frac{ \delta N }{2} + k \cdot \pi \right) \\ &= \sin \left( \frac{ \delta N }{2} \right) \cos \left( k \cdot \pi \right) + \cos \left( \frac{ \delta N }{2} \right) \sin \left( k \cdot \pi \right) \\ &= \sin \left( \frac{ \delta N }{2} \right) \cdot (-1)^k + \cos \left( \frac{ \delta N }{2} \right) \cdot 0 \\ &= (-1)^k \sin \left( \frac{ \delta N }{2} \right) \end{aligned} \tag {15} $$ Similar substitutions can be made into both sine terms in the numerator inside the brackets. The $(-1)^k$ factored out from inside the brackets cancels the one from outside the brackets. The end result is that (12) is now expressed in terms of $\delta$. $$ Z_k = \frac{ M }{ 2 } \cdot \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) } \cdot \left[ \frac{ \sin \left( \frac{ \delta N }{2} + \phi \right) e^{ i\beta_k } - \sin \left( \frac{ \delta N }{2} - \delta - \beta_k + \phi \right) }{ \sin \left( \frac{ \delta }{2} + \beta_k \right) } \right] \tag {16} $$ The numerator can be simplified by applying another substitution variable. $$ \rho = \frac{ \delta N }{2} + \phi \tag {17} $$ The final result is this: $$ Z_k = \frac{ M }{ 2 } \cdot \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) } \cdot \left[ \frac{ \sin \left( \rho \right) e^{ i\beta_k } - \sin \left( \rho - \delta - \beta_k \right) }{ \sin \left( \frac{ \delta }{2} + \beta_k \right) } \right] \tag {18} $$
An Approximation for Near Integer Frequencies
When $\delta$ is extremely close to zero, the sine quotient outside the brackets can be replaced by an approximation. It's hard to imagine a real life case where this would be practical, but it is important for analyzing the discontinuity at integer frequencies. The first two terms of the Taylor series expansion for the sine function is: $$ \sin( \theta ) \approx \theta - \frac{ \theta ^ 3 }{ 3! } \tag {19} $$ This approximation can be plugged into the sine quotient to get an approximation for it. $$ Q = \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) } \approx \frac{ \frac{ \delta N }{2} - \frac{ \left( \frac{ \delta N }{2} \right) ^ 3 }{ 3! } }{ N \left[ \frac{ \delta }{2} - \frac{ \left( \frac{ \delta }{2} \right) ^ 3 }{ 3! } \right] } =\frac{ 24 - ( \delta N ) ^ 2 }{ 24 - \delta ^2 } \approx 1 \tag {20} $$
Plugging the Discontinuity
From the approximation, it is clear that: $$ \lim_{\delta \to 0} Q = 1 \tag {21} $$ Also, from (17): $$ \lim_{\delta \to 0} \rho = \phi \tag {22} $$ Using these, the limit for (18) can be found. $$ \lim_{\delta \to 0} Z_k = \frac{ M }{ 2 } \cdot \left[ \frac{ \sin \left( \phi \right) e^{ i\beta_k } - \sin \left( \phi - \beta_k \right) }{ \sin \left( \beta_k \right) } \right] \tag {23} $$ This limit is a complex value and can be expressed using its constituent values. $$ \lim_{\delta \to 0} Z_k = x_k + i y_k \tag {24} $$ The exponent in the numerator can be also be split into a real and imaginary part by using Euler's equation: $$ e^{ i\beta_k } = \cos( \beta_k ) + i \sin( \beta_k ) \tag {25} $$ Therefore, the following equations for the real and imaginary parts can be derived: $$ \begin{aligned} x_k &= \frac{ M }{ 2 } \cdot \left[ \frac{ \sin \left( \phi \right) cos( \beta_k ) - \sin \left( \phi - \beta_k \right) }{ \sin \left( \beta_k \right) } \right] \\ &= \frac{ M }{ 2 } \cdot \left[ \frac{ \sin ( \phi ) \cos( \beta_k ) - \sin ( \phi ) \cos( \beta_k ) + \cos ( \phi ) \sin( \beta_k ) }{ \sin \left( \beta_k \right) } \right] \\ &= \frac{ M }{ 2 } \cos ( \phi ) \end{aligned} \tag {26} $$ $$ y_k = \frac{ M }{ 2 } \cdot \left[ \frac{ \sin \left( \phi \right) \sin( \beta_k ) }{ \sin \left( \beta_k \right) } \right] = \frac{ M }{ 2 } \sin ( \phi ) \tag {27} $$ These results can be recombined into a complex value by once again applying Euler's equation. $$ \lim_{\delta \to 0} Z_k = x_k + i y_k = \frac{ M }{ 2 } \cos ( \phi ) + i \frac{ M }{ 2 } \sin ( \phi ) = \frac{ M }{ 2 } e^{ i \phi } \tag {28} $$ Not surprisingly, this is exactly the well known value for a bin value of a pure real tone with an integer frequency. Another proof that this is the correct value is equation (19) from my original bin value formula blog article[1].
Conclusion
An alternative equation for bin values for a pure tone (12) is derived that should be used for a bin value if the frequency is very close to that bin index. Since both forms of the equation are exact, the cutoff value of which to use can be rather arbitrary. The larger concept to be understood from all the algebraic gymnastics around the equations' discontinuity is that rather than being considered as an ideal case, tones with integer frequencies are actually more properly considered as degenerate cases.
References
[1] Dawg, Cedron, DFT Bin Value Formulas for Pure Real Tones
[2] Dawg, Cedron, Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1
- Comments
- Write a Comment Select to add a comment
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.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: