A Two Bin Solution
Introduction
This is an article to hopefully give a better understanding of the Discrete Fourier Transform (DFT) by showing an implementation of how the parameters of a real pure tone can be calculated from just two DFT bin values. The equations from previous articles are used in tandem to first calculate the frequency, and then calculate the amplitude and phase of the tone. The approach works best when the tone is between the two DFT bins in terms of frequency.
The Coding Language
The included sample code is written in Gambas, a modern Linux based BASIC dialect. The purpose of the code is to demonstrate the math in an unambiguous manner, readable by programmer and non-programmer alike. The choice of Gambas is good because it is readable and supports complex values and the needed vector operations in nice clean syntax. Non-programmers should be able to follow the calculations the code does and programmers should be able to translate the code into the language/platform of their choice.
The Simple Testing
To prove the math works, a pure real tone is generated, the DFT bins calculated, and the generating parameters are recovered from the bin values in a two step process.
The test program is straightforward and doesn't require much explanation.
' Gambas module file '============================================================================= Public Sub Main() Dim N, m As Integer, theM, theAlpha, thePhi As Float '---- Set Test Parameters N = 16 theM = 1.234567 thePhi = 0.56789 theAlpha = 3.456789 * Pi(2) / N '---- Establish the DFT Calculator Dim theDFT As New DftClass(N) '---- Generate a Test Signal Dim theSignal As New Vector(N, False) For m = 0 To theSignal.Count - 1 theSignal[m] = theM * Cos(theAlpha * m + thePhi) Next '---- Show the DFT bin Values Dim b3, b4 As Complex b3 = theDFT.SingleNormalizedBin(3, theSignal) b4 = theDFT.SingleNormalizedBin(4, theSignal) Print "b3", b3, 2 * b3.Abs(), b3.Arg() Print "b4", b4, 2 * b4.Abs(), b4.Arg() Print '---- Calculate the Frequency Dim f As Float f = theDFT.TwoBinRealFreq(3, b3, b4) Print f '---- Calculate the Amplitude and Phase Dim theAmpAndPhase As Float[] theAmpAndPhase = theDFT.TwoBinRealAmpAndPhase(3, b3, b4, f) Print theAmpAndPhase[0] Print theAmpAndPhase[1] End '=============================================================================The "False" in the "New Vector" creation means don't make it complex, i.e. theSignal is real valued.
Gambas can't return multiple values, but it can return arrays so theAmpAndPhase is an array of the two values.
The Normalized DFT
A $1/N$ normalized DFT can be defined as: $$ Z_k = \frac{ 1 }{ N } \sum_{n=0}^{N-1} { S_n e^{ -i \frac{2\pi}{N} k n } } \tag {1} $$ Think of a pizza with N slices. Okay, now back to the subject at hand. $$ Z_k = \frac{1}{N} \sum_{n=0}^{N-1} { S_n \left[ e^{ -i \frac{2\pi}{N} k } \right] ^ n } = \frac{1}{N} \sum_{n=0}^{N-1} { S_n \left( R_k \right) ^ n } \tag {2} $$ You can think of the "R" as standing for "Rotate" because that's how it works. $$ R_k = e^{ -i \frac{2\pi}{N} k } = \left( e^{ -i \frac{2\pi}{N} } \right) ^ k = \left( R_1 \right) ^ k \tag {3} $$ You can also think of it as "Root" because it is one of the N roots of unity: $$ 1^{1/N} = \left( e^{ -i 2\pi k } \right)^{1/N} = e^{ -i \frac{2\pi}{N} k } = e^{ -i \frac{k}{N} 2\pi } \tag {4} $$ The latter two forms show the exponent formatted to be interpreted as "Rotate this many pizza slices" and "Fraction of the circle rotated".
The negative sign doesn't matter as the $k$ can be any integer. I just included it to make it more similar to the DFT and explain the negative sign in the code.
Establishing The DftClass Object
The code being called, the code actually doing the work, also speaks for itself, but a little extra explaining won't hurt.
The "new" routine is called when theDft object is created. The only parameter is the number of bins in the DFT which is also the number of sample points in the signal. A sample of N points from a larger set is usually called a DFT frame, and the frequencies of the tones in the DFT are calculated in units of cycles per frame.
From the frequency formula article[5] you have: $$ \vec C = \begin{bmatrix} [ \cos( \beta_k ) - \cos( \beta_j ) ] / \sqrt{2} \\ \sin( \beta_k ) \\ \sin( \beta_j ) \end{bmatrix} \tag {5} $$ $$ \vec P_{k,j} = \frac{ \vec C }{ \| \vec C \| } = \frac{ \vec C }{ \sqrt{ \vec C\cdot \vec C } } \tag {6} $$ The $\vec P$ (pre-calculated) values are then normalized values of $\vec C$
The $\beta$s are the angles of the bin locations when arranged in a circle. $$ \beta_k = k \frac{2\pi}{N} \tag {7} $$ They are inversely related to their corresponding root values. $$ R_k = e^{-i \beta_k } = \cos( \beta_k ) - i \sin( \beta_k ) \tag {8} $$ The code gets the Cosine and Sine values for the $\beta$s from the real and imaginary parts of the myRoots elements[1].
' Gambas class file Static Private ourSqr2Recip As Float Private myN As Integer Private myRoots As Complex[] Private myP As Vector[] Private myBinToAngle As Float Private myAngleToBin As Float '============================================================================= Static Public Sub _init() ourSqr2Recip = 1 / Sqr(2) End '============================================================================= Public Sub _new(ArgN As Integer) Dim k As Integer, theAngle As Float myN = argN '---- Set Frequency Conversion Factors (Angle = Radians per Sample) myBinToAngle = Pi(2) / ArgN myAngleToBin = 1.0 / myBinToAngle '---- Set the Roots myRoots = New Complex[ArgN] theAngle = 0 For k = 0 To myRoots.Max myRoots[k] = Complex.Polar(1.0, -theAngle) theAngle += myBinToAngle ' One Bin Slice Width Next '---- Precalculated Normalized C's Dim C As Vector myP = New Vector[ArgN / 2 + 1] For k = 0 To myP.Max - 1 C = New Vector(3, False) C[0] = (myRoots[k].Real - myRoots[k + 1].Real) * ourSqr2Recip C[1] = -myRoots[k].Imag C[2] = -myRoots[k + 1].Imag myP[k] = C / Sqr(C.Dot(C)) Next End '=============================================================================The conversion factors convert between pizza slices (bin spacing) and the distance travelled on the circumference of a unit pizza.
A Single DFT Bin Calculation
This is where the natures of the roots of unity as the rotation factor for a finite set of rotations comes into play. For the $k$th bin, the rotation size of the DFT $R_k$ is $k$ pizza slices. So for each additional sample point being added in, the multplier (rotator) advances $k$ slices, which means the subscript in the array advances $k$ spots. When it goes too far, it gets wrapped around to the start.
If you can understand the relationship between the DFT definition and this code implementation of it, they you will really understand how the DFT actually works. I try to explain it for a real valued signal in my second blog article[2].
'============================================================================= Public Sub SingleNormalizedBin(ArgK As Integer, ArgSignal As Vector) As Complex Dim theSum As Complex, r, n As Integer theSum = 0 r = 0 For n = 0 To myRoots.Max theSum += ArgSignal[n] * myRoots[r] r += ArgK If r >= myN Then r -= myN Next Return theSum / myN End '=============================================================================
The Frequency Calculation
Why or how this works is not going to be obvious by inspection. You really need to read the article[5] for that. How the code reflects the math is straightforward. $$ \vec A = \begin{bmatrix} ( x_k - x_j ) / \sqrt{2} \\ y_k \\ y_j \end{bmatrix} \tag {9} $$ $$ \vec B = \begin{bmatrix} [ \cos( \beta_k ) \cdot x_k - \cos( \beta_j ) \cdot x_j ] / \sqrt{2} \\ \cos( \beta_k ) \cdot y_k \\ \cos( \beta_j ) \cdot y_j \end{bmatrix} \tag {10} $$ $$ \vec D = \vec A + \vec B \tag {11} $$ $$ \vec K = \vec D - \left( \vec D \cdot \frac{ \vec C }{ \| \vec C \| } \right) \frac{ \vec C }{ \| \vec C \| } = \vec D - ( \vec D \cdot \vec P_{k,j} ) \vec P_{k,j} \tag {12} $$ The cosine of the frequency term is what is calculated. $$ \cos( \alpha ) = \frac{ \vec K \cdot \vec B }{ \vec K \cdot \vec A } \tag {13} $$ The frequency is then recovered in units of radians per sample, which is then converted to cycles per frame as the returned result. $$ f = \alpha \cdot \frac{ N }{ 2\pi } = \cos^{-1} \left( \frac{ \vec K \cdot \vec B }{ \vec K \cdot \vec A } \right) \cdot \frac{ N }{ 2\pi } \tag {14} $$
'============================================================================= Public Sub TwoBinRealFreq(ArgLowerIndex As Integer, ArgLowerValue As Complex, ArgUpperValue As Complex) As Float Dim A, B, D, K As New Vector(3, False) A[0] = (ArgLowerValue.Real - ArgUpperValue.Real) * ourSqr2Recip A[1] = ArgLowerValue.Imag A[2] = ArgUpperValue.Imag B[0] = (myRoots[ArgLowerIndex].Real * ArgLowerValue.Real - myRoots[ArgLowerIndex + 1].Real * ArgUpperValue.Real) * ourSqr2Recip B[1] = myRoots[ArgLowerIndex].Real * ArgLowerValue.Imag B[2] = myRoots[ArgLowerIndex + 1].Real * ArgUpperValue.Imag D = A + B K = D - D.Dot(myP[ArgLowerIndex]) * myP[ArgLowerIndex] Return ACos(K.Dot(B) / K.Dot(A)) * myAngleToBin End '=============================================================================
The Amplitude and Phase Calculation
This section is the complicated one. It would be possible to combine this routine into the previous one, but leaving them separate increases flexibility and the name clashes of the A and B vectors might be confusing. Separated, the routines can retain naming that closely matches the source articles.
The A and B in this routine don't exactly match the A and B from article[4] either. In the article, a patch of the DFT is treated as a complex vector. In the code, the patch is unfurled into a real valued vector, double the length. In this case, a pair of complex bins becomes a length four real valued vector. This change has been made so the calculated coefficients (a and b in the article, ca and cb in the code) which are expected to be real valued are forced to be real values. Since the two unused imaginary parts have been eliminated they can't "absorb the error" and the results are much improved in the presence of noise.
This routine suffers greatly if the frequency of the measured tone is very close to a bin. If the frequency is at a bin, then the values of U and V for both A and B will be zero. This is the degenerate case for this math. It is the ideal case for reading a DFT directly. In cases where it is close, the formulas here will work well for bins away from the frequency bin, but there is a variation of the formula that allows accurate numerically stable calculations for the value of the bin nearest the frequency.
The real pure tone, a sinusoidal, is defined as: $$ \begin{aligned} S_n &= M \cos( \alpha n + \phi ) \\ &= M \cos( \phi ) \cos( \alpha n ) - M \sin( \phi ) \sin( \alpha n ) \\ &= a \cdot \cos( \alpha n ) + b \cdot \sin( \alpha n ) \end{aligned} \tag {15} $$ Where: $$ \begin{aligned} a &= M \cos( \phi ) \\ b &= -M \sin( \phi ) \end{aligned} \tag {16} $$ The DFT is a linear transform. So the following is true: $$ \begin{aligned} Z &= F(S) \\ &= F[ a \cos( \alpha k ) ] + b \sin( \alpha k ) ] \\ &= a \cdot F[ \cos( \alpha k ) ] + b \cdot F[ \sin( \alpha k ) ] \\ &= a \cdot F[ \cos( \alpha k ) ] + b \cdot F[ \cos( \alpha k - \pi/2) ] \end{aligned} \tag {17} $$ The DFT values of the two basis signals can be calculated using the Bin Value formulas[3] for pure tones: $$ \begin{aligned} U &= \cos( \alpha N + \phi ) - \cos(\phi) \\ V &= \cos( \alpha N -\alpha + \phi ) - \cos( -\alpha + \phi ) \end{aligned} \tag {18} $$ $$ Z_k = \frac{M}{2N} \left[ \frac{ Ue^{ i\beta_k } - V } { \cos( \alpha ) - \cos( \beta_n ) } \right] \tag {19} $$ This formula gets used twice. Once for A and once for B. The $Z_k$ becomes $A$ and $B$ to form the two basis vectors. If the frequency is really close to an integer number of cycles per frame, then the value at the frequency bin should be calculated using the formulas in this article[6] instead. This is not done in the code. The omission is for overall clarity. For robustness, the code should be enhanced. $$ \begin{aligned} A &= F[ \cos( \alpha n ) ] \\ B &= F[ \sin( \alpha n ) ] = F[ \cos( \alpha n - \pi / 2 ) ] \end{aligned} \tag {20} $$ This will calculate A and B. Now the problem is to solve for $a$ and $b$. $$ Z = a A + b B \tag {21} $$ Here the Z is the DFT of source signal. This vector is then dotted with the basis vectors. $$ \begin{aligned} A \cdot Z &= a ( A \cdot A ) + b ( A \cdot B ) \\ B \cdot Z &= a ( B \cdot A ) + b ( B \cdot B ) \end{aligned} \tag {22} $$ After taking the dot products, this pair of equations is a simple case of two linear equations with two unknowns. The solution is: $$ a = \frac{ ( B \cdot B )( A \cdot Z ) - ( A \cdot B )( B \cdot Z ) } { ( B \cdot B )( A \cdot A ) - ( A \cdot B )( B \cdot A ) } \tag {23} $$ $$ b = \frac{ ( A \cdot A )( B \cdot Z ) - ( B \cdot A )( A \cdot Z ) } { ( A \cdot A )( B \cdot B ) - ( B \cdot A )( A \cdot B ) } \tag {24} $$ Thanks to unfurling the complex bin values into real valued vectors, the need to deal with the conjugates is gone.
Finally, once $a$ and $b$ are calculated, they can be converted into polar form. $$ M = \sqrt { a^2 + b^2 } \tag {25} $$ $$ \phi = atan2( -b, a ) \tag {26} $$ These are amplitude and phase of the tone.
'============================================================================= Public Sub TwoBinRealAmpAndPhase(ArgLowerIndex As Integer, ArgLowerValue As Complex, ArgUpperValue As Complex, ArgFrequency As Float) As Float[] Dim A, B, Z As New Vector(4, False) Dim theAlpha, theAlphaN, UA, VA, UB, VB, fj, fk As Float Dim AB, AZ, BZ, AA, BB, d, ca, cb As Float Dim theResults As New Float[2] Dim theCosineOfAlpha, theCosineOfBetaJ, theCosineOfBetaK As Float Dim theSineOfBetaJ, theSineOfBetaK As Float '---- Set local variables theAlpha = ArgFrequency * myBinToAngle theAlphaN = theAlpha * myN theCosineOfAlpha = Cos(theAlpha) theCosineOfBetaJ = myRoots[ArgLowerIndex].Real theSineOfBetaJ = -myRoots[ArgLowerIndex].Imag theCosineOfBetaK = myRoots[ArgLowerIndex + 1].Real theSineOfBetaK = -myRoots[ArgLowerIndex + 1].Imag '---- Calculate the Bin Value Parameters UA = Cos(theAlphaN) - 1.0 VA = Cos(theAlphaN - theAlpha) - theCosineOfAlpha UB = Sin(theAlphaN) VB = Sin(theAlphaN - theAlpha) + Sin(theAlpha) '---- Calculate the Basis Bin Values fj = 1.0 / (2.0 * myN * (theCosineOfAlpha - theCosineOfBetaJ)) ' Factor fk = 1.0 / (2.0 * myN * (theCosineOfAlpha - theCosineOfBetaK)) A[0] = fj * (UA * theCosineOfBetaJ - VA) A[1] = fj * (UA * theSineOfBetaJ) A[2] = fk * (UA * theCosineOfBetaK - VA) A[3] = fk * (UA * theSineOfBetaK) B[0] = fj * (UB * theCosineOfBetaJ - VB) B[1] = fj * (UB * theSineOfBetaJ) B[2] = fk * (UB * theCosineOfBetaK - VB) B[3] = fk * (UB * theSineOfBetaK) '---- Unfurl the Passed Bin Values as a real valued Vector Z[0] = ArgLowerValue.Real Z[1] = ArgLowerValue.Imag Z[2] = ArgUpperValue.Real Z[3] = ArgUpperValue.Imag '---- Calculate the necessary dot products AB = A.Dot(B) AZ = A.Dot(Z) BZ = B.Dot(Z) AA = A.Dot(A) BB = B.Dot(B) '---- Calculate the Linear Coefficients d = AA * BB - AB * AB ca = (BB * AZ - AB * BZ) / d cb = (AA * BZ - AB * AZ) / d '---- Convert Cartesian Results to Polar Form theResults[0] = Sqr(ca * ca + cb * cb) ' Amplitude theResults[1] = ATan2(-cb, ca) ' Phase Return theResults End '=============================================================================
The Results
The thing to notice about the bin values is that they are nearly completely out of phase. This is what is expected when a pure tone is between two peaks. The results are accurate to the precision of the variables. This is reflected with a little bit of an error in the phase value. In any real life application, the level of noise will exceed the variable precision limit by orders of magnitude so this isn't a concern at all.
b3 -0.113598594199752+0.375122610206239i 0.783891863185482 1.86484794794942 b4 0.217236372698119-0.327922570624235i 0.786701605306211 -0.985710020979189 3.456789 1.234567 0.567890000000001The frequency, amplitude, and phase are correctly calculated, even though there were only sixteen signal points to work with and only two DFT bins were used. Even fewer signal points could have been used, but why push the point, the more samples the less noise matters.
Conclusion
The parameter values of an unknown real pure tone can be calculated from just two of its DFT bins. This solution does not require a large number of signal points to be sampled unlike most other approaches to the same problem. Also, this one is distinct in being exact, not an approximation, in the noiseless case. As long as a peak is distinct in a DFT (away from other peaks, above tne noise floor), this technique will deliver good results.
References
[1] Dawg, Cedron, The Exponential Nature of the Complex Unit Circle
[2] Dawg, Cedron, DFT Graphical Interpretation: Centroids of Weighted Roots of Unity
[3] Dawg, Cedron, DFT Bin Value Formulas for Pure Real Tones
[4] Dawg, Cedron, Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1
[5] Dawg, Cedron, Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
[6] Dawg, Cedron, An Alternative Form of the Pure Real Tone DFT Bin Value Formula
The source code can be found here: https://forum.gambas.one/viewtopic.php?f=4&t=728
- 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: