A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT

• Basics
• DFT_FFT
• Tutorials

This is an article to hopefully give a better understanding to the Discrete Fourier Transform (DFT) by deriving an exact formula for the frequency of a complex tone in a DFT. It is basically a parallel treatment to the real case given in Exact Frequency Formula for a Pure Real Tone in a DFT. Since a real signal is the sum of two complex signals, the frequency formula for a single complex tone signal is a lot less complicated than for the real case.

Similar to the real case, the bin value formula is used on two DFT bins to create a system of two equations (implicitly four, if you count the real and imaginary parts) with three unknowns. A single combination of the two equations eliminates two of the unknowns leaving only the frequency term as unknown. This term is easily isolated and expressed in terms of the known values. From this, the frequency can be calculated. Since the most practical application of this equation will be at the two consecutive bins with the largest magnitude, a special version of the formula is derived for that usage.

A pure complex tone is an exponential function that can be defined by three parameters. The same three as for the real case, they are amplitude ($ M $), frequency ($ f $), and phase ($ \phi $). A discrete pure complex tone is a set of complex signal values within a frame. A frame consists of a defined number of points, $ N $, and the signal values are indexed from $ 0 $ to $ N-1 $. The signal can be extended beyond the frame by using out of range subscripts. Here is an equation defining a pure complex tone: $$ S_n = M e^{i( f \cdot { 2\pi \over N } n + \phi )} \tag 1 $$ The frequency is expressed in units of cycles per frame. For practical purposes the frequency should be constrained to be greater than zero and less than $ N $. It is also possible to use a range of $ -N/2 $ to $ N/2 $. Negative frequencies make a lot more sense for complex tones than for real tones due to the sign determining the chirality of the tone in the time domain. Frequencies outside of the ranges will appear to be frequencies within the range in the DFT. This phenomenon is known as aliasing.

The frequency and bin indexes can be rescaled to a distance along the circumference which is also an angular scale. Defining angle variables makes working with the equations simpler. $$ \alpha = f \cdot { 2\pi \over N } \tag 2 $$ $$ \beta_k = k \cdot { 2\pi \over N } \tag 3 $$ Substituting (2) into (1) gives a simpler formula for the signal. $$ S_n = M \cdot { e^{ i(\alpha n + \phi) } } \tag 4 $$

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

For non-integer frequencies a second bin is needed to solve for the frequency. The first bin has the index $ k $ and the second bin has the index $ j $. The formula for the $ j $ bin is the same as that for the $ k $ bin, except $ j $ is substituted for $ k $. $$ Z_j = { M \over N } e^{ i \phi } \cdot { 1 - e^{ i\alpha N } \over 1 - e^{ i(\alpha - \beta_j ) } } \tag {6} $$ (5) can be divided by (6) to give this equation: $$ { Z_k \over Z_j } = { 1 - e^{ i(\alpha - \beta_j ) } \over 1 - e^{ i(\alpha - \beta_k ) } } \tag {7} $$ This operation eliminates $ M $ and $ \phi $, leaving $ \alpha $ as the only unknown. The next step is to cross multiply and separate the exponent terms. $$ Z_k - Z_k e^{ i\alpha } e^{ -i\beta_k } = Z_j - Z_j e^{ i\alpha } e^{ -i\beta_j } \tag {8} $$ The frequency related term $( e^{ i\alpha } )$ occurs twice in the equation. It can be solved for by using ordinary algebra. $$ e^{ i\alpha } = { { Z_k - Z_j } \over { e^{ -i\beta_k } Z_k - e^{ -i\beta_j } Z_j } } = Q_1 \tag {9} $$ The result is that the only unknown is on the left side of the equation and the right side contains known values and can be evaluated.

The integer frequency case is a special case where all the DFT bin values will be zero except for the bin value corresponding to the frequency. The formula still works as long as $ Z_k $ or $ Z_j $ is the non-zero bin. If they are both zero, of course the equation won't work. Assuming, without loss of generality, that the $ k $ bin has the non-zero value. It can be seen by inspection that the $ Z_k $ value will cancel out of the numerator and the denominator leaving a negative exponent in the denominator. Flipping the reciprocal means flipping the sign in the exponent leaving $ e^{ i\alpha } = e^{ i\beta_k } $. Thus $ \alpha $ has to be equal to $ \beta $ plus an integer multiple of $ 2 \pi $. Assuming the signal is not an alias, this means $ f = k $. To translate to a negative frequency range, if $ k > N/2 $ then $ f = k - N $.

In order to solve for $ \alpha $ in (9) it is necessary to take the natural log of a complex number. To understand what this means it is useful to consider a complex number in both its Cartesian form and its polar form and how the two expressions are related. $$ a + ib = r e^{ i\theta } \tag {10} $$ A point in the complex plane $(a,b)$ corresponds to the complex value $ a + ib $. This representation is unique. It can also be located by the angle from the real axis $( \theta )$ and the distance from the origin $(r)$. This representation is not unique. Any integer multiple of $2\pi$ added to the angle value yields the same angle. The conversion from polar to Cartesian is straightforward trigonometry. $$ a = r \cdot \cos( \theta ) \tag {11} $$ $$ b = r \cdot \sin( \theta ) \tag {12} $$ It follows from Euler's equation: $$ e^{ i\theta } = \cos( \theta ) + i \cdot \sin( \theta ) \tag {13} $$ $$ r e^{ i\theta } = r \cos( \theta ) + i \cdot r \sin( \theta ) = a + ib \tag {14} $$ An explanation, and derivation, of Euler's equation can be found in my first blog article titled "The Exponential Nature of the Complex Unit Circle".

The conversion from Cartesian to polar uses the Pythagorean theorem for the distance to the origin. $$ r = ( a^2 + b^2 )^{ 1/2 } \tag {15} $$ The angle is usually solved for using the Arctangent function. However, this method may require an adjustment of +/- $\pi$ depending upon in what quadrant $(a,b)$ lies. $$ \theta = \tan^{-1} \left( { b \over a } \right) \tag {16} $$ In most math library implementations, there is an alternative function called "atan2". Instead of taking the ratio of $ b/a $ as an argument, it takes $a$ and $b$ as separate arguments. The function properly handles the quadrant adjustments and handles the $a=0$ cases. $$ \theta = \operatorname{ atan2 }( b, a ) \tag {17} $$ The natural log of (10) can now be taken. $$ \ln( a + ib ) = \ln( r ) + i\theta \tag {18} $$ Substituting in (15) and (17) leaves a working definition which uses only $a$ and $b$. $$ \ln( a + ib ) = \ln( a^2 + b^2 ) / 2 + i \cdot \operatorname{ atan2 }( b, a ) \tag {19} $$

Taking the natural log of (9) leaves this equation: $$ i\alpha = \ln(Q_1) \tag {20} $$ If the bin values correspond to a pure tone the magnitude of $ Q_1 $ will be 1 so the real part of the logarithm will be zero. If the tone's amplitude is growing or decaying exponentially, the rate will be reflected in $ \ln(\|Q_1\|) $, which is the real part of $\ln(Q_1)$. This will be covered in a future blog article.

Since a zero value is expected, the real part of the natural log is discarded and we are left equating the imaginary parts of the logarithm which correspond with the frequency. $$ \alpha = \operatorname{ atan2 }( Im[Q_1], Re[Q_1] ) \tag {21} $$ $Im[Q_1]$ and $Re[Q_1]$ are the imaginary and real parts of $Q_1$. There is an implicit integer multiple of $2\pi$ which can be added to the angle. This accounts for all the alias frequencies that are possible. Converting back from an angle scale to a bin index scale results in the frequency expressed as cycles per frame. The atan2 function will return a value within the range of $-\pi$ to $\pi$, so the frequency will be from $-N/2$ to $N/2$. $$ f = \alpha \cdot { N \over 2\pi } = \operatorname{ atan2 }( Im[Q_1], Re[Q_1] ) \cdot { N \over 2\pi } \tag {22} $$ If a frequency range of 0 to $N$ is desired, negative frequencies can have N added to them. This corresponds to a single $2\pi$ on the angle scale.

These numerical examples demonstrate the validity of the formula. The source code can be found below in Appendix A. First is a listing of the DFT values and the Roots of Unity used in the DFT calculation and the frequency formula. Following this are three calculations using three different bin pairs. The first pair are the DC bin and the Nyquist bin. The second pair are an arbitrary selection of far apart spaced bins. The third pair are the adjacent bins to the peak value. The calculations are done to a greater precision than the display values.

               DFT                    e^(-i*beta_k)
 k       Real        Imag            Real       Imag
--  --------------------------  --------------------------
 0   0.005657918  0.070641787    1.000000000 -0.000000000
 1  -0.009392193  0.080305352    0.923879533 -0.382683432
 2  -0.030314566  0.093739453    0.707106781 -0.707106781
 3  -0.064623699  0.115769094    0.382683432 -0.923879533
 4  -0.139603486  0.163913062    0.000000000 -1.000000000
 5  -0.500205999  0.395453272   -0.382683432 -0.923879533
 6   0.567743787 -0.290269399   -0.707106781 -0.707106781
 7   0.207141274 -0.058729189   -0.923879533 -0.382683432
 8   0.132161487 -0.010585221   -1.000000000 -0.000000000
 9   0.097852354  0.011444420   -0.923879533  0.382683432
10   0.076929981  0.024878521   -0.707106781  0.707106781
11   0.061879870  0.034542087   -0.382683432  0.923879533
12   0.049722484  0.042348255   -0.000000000  1.000000000
13   0.038948748  0.049265991    0.382683432  0.923879533
14   0.028589040  0.055917882    0.707106781  0.707106781
15   0.017815305  0.062835619    0.923879533  0.382683432

***************************************************
Example Case:   k = 0   j = 8

Num = (0.0057+0.0706i) - (0.1322-0.0106i)
    = (-0.1265+0.0812i)

Den = (1.0000-0.0000i)(0.0057+0.0706i) - (-1.0000-0.0000i)(0.1322-0.0106i)
    = (0.0057+0.0706i) - (-0.1322+0.0106i)
    = (0.1378+0.0601i)

Q1 = Num/Den = (-0.5556+0.8315i)

||Q1|| = 1.0000   Alpha = atan2(Im[Q1],Re[Q1]) = 2.1598

f = Alpha*N/(2Pi) = 5.5000

***************************************************
Example Case:   k = 2   j = 14

Num = (-0.0303+0.0937i) - (0.0286+0.0559i)
    = (-0.0589+0.0378i)

Den = (0.7071-0.7071i)(-0.0303+0.0937i) - (0.7071+0.7071i)(0.0286+0.0559i)
    = (0.0448+0.0877i) - (-0.0193+0.0598i)
    = (0.0642+0.0280i)

Q1 = Num/Den = (-0.5556+0.8315i)

||Q1|| = 1.0000   Alpha = atan2(Im[Q1],Re[Q1]) = 2.1598

f = Alpha*N/(2Pi) = 5.5000

***************************************************
Example Case:   k = 5   j = 6

Num = (-0.5002+0.3955i) - (0.5677-0.2903i)
    = (-1.0679+0.6857i)

Den = (-0.3827-0.9239i)(-0.5002+0.3955i) - (-0.7071-0.7071i)(0.5677-0.2903i)
    = (0.5568+0.3108i) - (-0.6067-0.1962i)
    = (1.1635+0.5070i)

Q1 = Num/Den = (-0.5556+0.8315i)

||Q1|| = 1.0000   Alpha = atan2(Im[Q1],Re[Q1]) = 2.1598

f = Alpha*N/(2Pi) = 5.5000

Since the best use of the formula is on adjacent bins, it makes sense to derive a variation for that case. Let $k$ be the lower bin and $j$ be the upper one. $$ j = k + 1 \tag {21} $$ Substituting this into (9) and factoring out a $e^{ -i\beta_k }$ from the denominator leaves this: $$ e^{ i\alpha } = e^{ i\beta_k } { { Z_k - Z_{k+1} } \over { Z_k - e^{ -i\beta_1 } Z_{k+1} } } \tag {22} $$ Both sides can be divided by $e^{ i\beta_k }$ leaving another quotient expression on the right hand side of the equation. $$ e^{ i( \alpha - \beta_k ) } = { { Z_k - Z_{k+1} } \over { Z_k - e^{ -i\beta_1 } Z_{k+1} } } = Q_2 \tag {23} $$ Once again, take the natural log of both sides and discard the real part which will be zero anyway in the pure tone case. $$ \alpha - \beta_k = \operatorname{ atan2 }( Im[Q_2], Re[Q_2] ) \tag {24} $$ 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[Q_2], Re[Q_2] ) \cdot { N \over 2\pi } \tag {25} $$ The same issues about frequency range and aliases apply to this equation as well. Calculation wise, a complex multiply has been swapped for a real add, so it is marginally more efficient computationally. Because the expected angle from the atan2 function is going to be within a bin width, the Arctangent function can be used instead. This means for certain applications, it may be more efficient to evaluate the first few terms of the Arctangent Taylor series when a good approximation is all that is required.

The derivation of these formulas is not that difficult. The result is a two bin formula for complex tones that is relatively simple. In applications where there is foreknowledge of the approximate frequency of a signal, only the DFT bins in the neighborhood of the frequency need to be calculated.

#include <math.h>
#include <stdio.h>

struct complex
{
    double Real;
    double Imag;
};

//============================================================================
void complexProduct( complex* z1, complex* z2, complex* product )
{
        product->Real = z1->Real * z2->Real - z1->Imag * z2->Imag;
        product->Imag = z1->Real * z2->Imag + z1->Imag * z2->Real;
}
//============================================================================
void complexQuotient( complex* z1, complex* z2, complex* quotient )
{
        double den2 = z2->Real * z2->Real + z2->Imag * z2->Imag;

        quotient->Real = ( z1->Real * z2->Real + z1->Imag * z2->Imag ) / den2;
        quotient->Imag = (-z1->Real * z2->Imag + z1->Imag * z2->Real ) / den2;
}
//============================================================================
void showCalcTwoBins( int N, int k, int j, complex* bins, complex* roots )
{
        complex Num, Den, Product_j, Product_k, Q1;
        
        Num.Real = bins[k].Real - bins[j].Real;
        Num.Imag = bins[k].Imag - bins[j].Imag;

        complexProduct( &roots[k], &bins[k], &Product_k );
        complexProduct( &roots[j], &bins[j], &Product_j );

        Den.Real = Product_k.Real - Product_j.Real;
        Den.Imag = Product_k.Imag - Product_j.Imag;

        complexQuotient( &Num, &Den, &Q1 ); 
        
        double Q1_Mag = sqrt(  Q1.Real* Q1.Real + Q1.Imag * Q1.Imag );
        double Alpha = atan2( Q1.Imag, Q1.Real );
        
        double f = Alpha * N / ( 2.0 * M_PI );
                
        printf( "\n" );
        printf( "***************************************************\n" );
        printf( "Example Case:   k = %d   j = %d\n", k, j );
        printf( "\n" );

        printf( "Num = (%6.4f%+6.4fi) - (%6.4f%+6.4fi)\n", 
                bins[k].Real, bins[k].Imag, 
                bins[j].Real, bins[j].Imag );

        printf( "    = (%6.4f%+6.4fi)\n", Num.Real, Num.Imag );

        printf( "\n" );

        printf( "Den = (%6.4f%+6.4fi)(%6.4f%+6.4fi) - (%6.4f%+6.4fi)(%6.4f%+6.4fi)\n", 
                roots[k].Real, roots[k].Imag, 
                bins[k].Real,  bins[k].Imag, 
                roots[j].Real, roots[j].Imag, 
                bins[j].Real,  bins[j].Imag );

        printf( "    = (%6.4f%+6.4fi) - (%6.4f%+6.4fi)\n", 
                Product_k.Real, Product_k.Imag, 
                Product_j.Real, Product_j.Imag );

        printf( "    = (%6.4f%+6.4fi)\n", Den.Real, Den.Imag );

        printf( "\n" );

        printf( "Q1 = Num/Den = (%6.4f%+6.4fi)\n", Q1.Real, Q1.Imag );

        printf( "\n" );

        printf( "||Q1|| = %6.4f   Alpha = atan2(Im[Q1],Re[Q1]) = %6.4f\n", 
                Q1_Mag, Alpha );

        printf( "\n" );

        printf( "f = Alpha*N/(2Pi) = %6.4f\n", f );

}
//============================================================================
void DFTbin( int N, int k, complex* signal, complex* roots, complex* result )
{
        complex Product;
        
        result->Real = 0.0;
        result->Imag = 0.0;
        
        int nk = 0;

        for( int n = 0; n < N; n++ )
        {
            complexProduct( &signal[n], &roots[nk], &Product );

            result->Real += Product.Real;
            result->Imag += Product.Imag;
            
            nk += k;
            if( nk >= N ) nk -= N;
        }

        result->Real /= (double) N;
        result->Imag /= (double) N;
}
//============================================================================
int main( int argc, char *argv[] )
{
        int    N    = 16;

        double freq = 5.5;
        double phi  = 1.0;
        double M    = 1.0;

        complex signal[N];

        double alpha = freq * 2.0 * M_PI / (double) N;

        for( int n = 0; n < N; n++ )
        {
            double angle = alpha * (double) n + phi;

            signal[n].Real = M * cos( angle );
            signal[n].Imag = M * sin( angle );
        }

        complex roots[N];

        for( int k = 0; k < N; k++ )
        {
            double beta_k = (double) k * 2.0 * M_PI / (double) N;
            
            roots[k].Real =  cos( beta_k );
            roots[k].Imag = -sin( beta_k );
        }

        printf( "               DFT                    e^(-i*beta_k)\n" );
        printf( " k       Real        Imag            Real       Imag\n" );
        printf( "--  --------------------------  --------------------------\n" );

        complex bins[N];

        for( int k = 0; k < N; k++ )
        {
            DFTbin( N, k, signal, roots, &bins[k] );
           
            printf( "%2d  %12.9f %12.9f   %12.9f %12.9f\n", 
                    k, 
                    bins[k].Real, 
                    bins[k].Imag,
                    roots[k].Real, 
                    roots[k].Imag );
        }

        showCalcTwoBins( N, 0,  8, bins, roots );
        showCalcTwoBins( N, 2, 14, bins, roots );
        showCalcTwoBins( N, 5,  6, bins, roots );

        return 0;
}
//============================================================================

Sign in