A Two Bin Solution

• Basics
• DFT_FFT
• Tutorials

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
'=============================================================================

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
'=============================================================================

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.567890000000001

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

Sign in