Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

• Basics
• DFT_FFT
• Tutorials

This is an article to hopefully give a better understanding of the Discrete Fourier Transform (DFT) by deriving exact formulas for the phase and amplitude of a non-integer frequency real tone in a DFT. The linearity of the Fourier Transform is exploited to reframe the problem as the equivalent of finding a set of coordinates in a specific vector space. The found coordinates are then used to calculate the phase and amplitude of the pure real tone in the DFT. This article is an extension of, and depends on, my previous two blog articles:

• DFT Bin Value Formulas for Pure Real Tones
• Exact Frequency Formula for a Pure Real Tone in a DFT

It would be helpful to read them first before reading this article.

The phase and amplitude can be solved for from the bin value formulas of a pure real tone in a similar manner that the frequency was in the previous blog. However, unlike the frequency formula, these solutions turn out to be sensitive to noise and other tones in the DFT. A different, more robust, approach is shown in this blog article. A section of the DFT, preferably surrounding the peak value, is treated like a vector. Since the frequency can be calculated, the bin value formulas can be used to create two basis vectors for the vector subspace. Starting from the definition of a tone, the DFT vector is shown to be within the span of the two basis vectors having real valued coordinates. Since the basis vectors are neither normalized nor othogonal, the problem of finding the coordinate values is reduced to being an ordinary linear system of two equations, two unknowns. Once the coordinates are known, the phase and amplitude are readily calculated from their definitions.

A pure real tone of finite length is a sinusoidal signal consisting of $ N $ evenly space points in time. This mathematical definition of the signal is the same as in the previous blog articles:

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

This can be simplified as:

$$ S_n = M cos( \alpha n + \phi ) \tag {2} $$

Where:

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

Using the angle addition formula for cosines, the phase component ( $ \phi $ ) can be broken out from the cosine argument.

$$ S_n = M cos( \phi ) cos( \alpha n ) - M sin( \phi ) sin( \alpha n ) \tag {4} $$

Since $ M $ and $ \phi $ are constants, the tone can now be seen as a linear combination of a cosine and sine signal. Giving each coefficient its own label makes this a lot clearer and will simplify working with these expressions.

$$ S_n = a \cdot cos( \alpha n ) + b \cdot sin( \alpha n ) \tag {5} $$

Where:

$$ a = M cos( \phi ) \tag {6} $$

$$ b = -M sin( \phi ) \tag {7} $$

When the DFT is applied on the signal the result is a set of complex values known as bins. These bins can also be considered the elements of a complex valued vector.

$$ Z = F[ S_n ] \tag {8} $$

Since the DFT ( $ F[\,] $ ) is a linear transform, when it is applied to the linear sum form of the signal, the sum and the coefficients can be brought outside the argument.

$$ Z = a \cdot F[ cos( \alpha n ) ] + b \cdot F[ sin( \alpha n ) ] \tag {9} $$

Giving the transforms of the two component signals their own labels will make working with them a lot easier.

$$ Z = a A + b B \tag {10} $$

Where:

$$ A = F[ cos( \alpha n ) ] \tag {11} $$

$$ B = F[ sin( \alpha n ) ] = F[ cos( \alpha n - \pi / 2 ) ] \tag {12} $$

$ A $ and $ B $ are complex basis vectors and $ Z $ falls within their span. Furthermore, from their definitions, $ a $ and $ b $ are known to be real valued. Equation (10) holds true whether Z is considered to be the whole DFT or any subset of it. This is because the vector equation really represents a set of equivalent scalar equations in parallel.

The first thing that needs to be known in order to construct $ A $ and $ B $ is the value of $ \alpha $. How to find $ \alpha $ from a three bin set of $ Z $ is the primary subject of the Exact Frequency Formula for a Pure Real Tone in a DFT blog article. From (3) it is clear that $ \alpha $ represents the frequency of the tone signal.

Once $ \alpha $ is known, $ A $ and $ B $ can be constructed. This can be done by generating the signal values and applying the DFT to them. It makes sense to limit $ Z $, $ A $, and $ B $ to the three bin set surrounding the peak value in the DFT. The nature of the DFT is to concentrate the information about a pure tone to the bins near its frequency. Even so, having to calculate the DFT values for just three bins for the two basis vectors from generated signals means a lot of calculations. The bigger the number of sample points ($ N $) is, the more calculations that are needed. The blog article DFT Bin Value Formulas for Pure Real Tones gives formulas to efficiently calculate the DFT values of pure real tones with the number of calculations independent of the number of sample points. These formulas greatly facilitate the construction of $ A $ and $ B $.

The only unknowns remaining in (10) are $ a $ and $ b $. Solving for them is a well known standard technique in Linear Algebra. The equation can be turned into two scalar equations by dotting it with two different vectors. With real valued vectors, one would use the basis vectors. With complex valued vectors, it is better to use the complex conjugates of the basis vectors. The complex conjugate is denoted by the bar across the top of the vector label.

$$ \bar A \cdot Z = a ( \bar A \cdot A ) + b ( \bar A \cdot B ) \tag {13} $$

$$ \bar B \cdot Z = a ( \bar B \cdot A ) + b ( \bar B \cdot B ) \tag {14} $$

The dot products are either real or complex valued scalars so (13) and (14) form a system of two linear equations with two unknowns. By using the conjugates of the basis vectors, $ \bar A \cdot A $ and $ \bar B \cdot B $ are guaranteed to be real valued.

Multiplying (13) by $ \bar B \cdot B $ and (14) by $ \bar A \cdot B $ yields the following two equations:

$$ ( \bar B \cdot B )( \bar A \cdot Z ) = a ( \bar B \cdot B )( \bar A \cdot A ) + b ( \bar B \cdot B )( \bar A \cdot B ) \tag {15} $$

$$ ( \bar A \cdot B )( \bar B \cdot Z ) = a ( \bar A \cdot B )( \bar B \cdot A ) + b ( \bar A \cdot B )( \bar B \cdot B ) \tag {16} $$

Notice that the $ b $ terms are now identical. Subtracting (16) from (15) eliminates the $ b $ terms and leaves $ a $ easy to solve for.

$$ a = { ( \bar B \cdot B )( \bar A \cdot Z ) - ( \bar A \cdot B )( \bar B \cdot Z ) \over ( \bar B \cdot B )( \bar A \cdot A ) - ( \bar A \cdot B )( \bar B \cdot A ) } \tag {17} $$

The value of $ b $ can be solved for in a similar manner, or using the symmetry between $ A $ and $ B $ in (10), derived directly from (17) by swapping the $ A $ and $ B $ occurrences.

$$ b = { ( \bar A \cdot A )( \bar B \cdot Z ) - ( \bar B \cdot A )( \bar A \cdot Z ) \over ( \bar A \cdot A )( \bar B \cdot B ) - ( \bar B \cdot A )( \bar A \cdot B ) } \tag {18} $$

These equations seem to require a lot of calculations, but a little inspection will show there are fewer than there appear to be. First off, the two denominators are the same so they need only be calculated once. Second $ \bar B \cdot A $ is the conjugate of $ \bar A \cdot B $ so only one dot product needs to be done and finding the other is trivial. Third, the denominator is real valued so a complex division is avoided. Fourth, the numerators share common terms that need only be calculated once. Fifth, calculating the product of conjugates means summing the squares of the real and imaginary parts and is simpler than an actual complex multiplication.

Once $ a $ and $ b $ have been solved for, calculating the amplitude and phase becomes straightforward. Using (6) and (7), the solution can be seen as being the equivalent of a cartesian to polar coordinates conversion.

$$ M = \sqrt { a^2 + b^2 } \tag {19} $$

$$ \phi = atan2( -b, a ) \tag {20} $$

Another way to look at it is $ a $ and $ b $ act like a rescaled "virtual DFT bin" at the tone's frequency.

The vectors $ A $, $ B $, and $ Z $ have to have at least two elements or (17) and (18) become indeterminate.  Similarly, the equations will also become indeterminate if the signal has a whole integer frequency.  In that case, of course, the phase and amplitude can be determined from the sole non-zero DFT bin the signal produces.

Here is a sample program to demonstrate the formulas. The signal and DFT values come from the blog article DFT Bin Value Formulas for Pure Real Tones. The value of $ \alpha $ is hard coded. In practice, it would be solved from the DFT values using the formulas in the Exact Frequency Formula for a Pure Real Tone in a DFT.

The source code for the ComplexDouble class is not shown. The functionality of the methods should be obvious from their names and the formulas being implemented.

#include <math.h>
#include <stdio.h>
#include "ComplexDouble.hpp"

//============================================================================
int main( int argc, char *argv[] )
{
        ComplexDouble Z[3];

        Z[0].setTo( -0.00032563186,   0.10802118551 );
        Z[1].setTo( -0.07619790924,   0.36944527683 );
        Z[2].setTo(  0.10202082457,  -0.23340312262 );

        double freq = 10.4;
        double N    = 32;
        double alpha = freq * 2.0 * M_PI / N;

        double UA = cos( alpha * N )         - 1.0;
        double VA = cos( alpha * N - alpha ) - cos( -alpha );

        double UB = sin( alpha * N );
        double VB = sin( alpha * N - alpha ) + sin( alpha );

        ComplexDouble A[3];
        ComplexDouble B[3];

        for( int s = 0; s < 3; s++ )
        {
            double k      = 9 + s;  
            double beta_k = k * 2.0 * M_PI / N;

            double factor =  1.0 / ( 2.0 * N * ( cos( alpha ) - cos( beta_k ) ) );

            A[s].setTo( factor * ( UA * cos( beta_k ) - VA ),
                        factor * ( UA * sin( beta_k ) ) );

            B[s].setTo( factor * ( UB * cos( beta_k ) - VB ),
                        factor * ( UB * sin( beta_k ) ) );
        }

        ComplexDouble A_B;
        ComplexDouble A_Z;
        ComplexDouble B_Z;

        A_B.setTo( 0.0, 0.0 );
        A_Z.setTo( 0.0, 0.0 );
        B_Z.setTo( 0.0, 0.0 );

        double A_A = 0.0;                
        double B_B = 0.0;                

        for( int s = 0; s < 3; s++ )
        {
            A_B.addConjugateLeftProductOf( &A[s], &B[s] );
            A_Z.addConjugateLeftProductOf( &A[s], &Z[s] );
            B_Z.addConjugateLeftProductOf( &B[s], &Z[s] );

            A_A += A[s].sumOfSquares();
            B_B += B[s].sumOfSquares();
        }

        ComplexDouble B_A;
        B_A.setToConjugateOf( &A_B );

        double DenomRecip = 1.0 / ( A_A * B_B - A_B.sumOfSquares() );

        ComplexDouble a;
        ComplexDouble b;

        a.setTo( &A_Z );
        a.rescaleBy( B_B );
        a.subtractProductOf( &A_B, &B_Z );
        a.rescaleBy( DenomRecip );

        b.setTo( &B_Z );
        b.rescaleBy( A_A );
        b.subtractProductOf( &B_A, &A_Z );
        b.rescaleBy( DenomRecip );

        double M   = sqrt( a.myReal * a.myReal + b.myReal * b.myReal );
        double phi = atan2( -b.myReal, a.myReal );

        printf( "        Real             Imaginary        Mag    Angle\n" );
        printf( "        --------------   --------------   ------ -------\n" );
        Z[0].debugPrint( "Z[0]" );
        Z[1].debugPrint( "Z[1]" );
        Z[2].debugPrint( "Z[2]" );
        printf( "\n" );
        A[0].debugPrint( "A[0]" );
        A[1].debugPrint( "A[1]" );
        A[2].debugPrint( "A[2]" );
        printf( "\n" );
        B[0].debugPrint( "B[0]" );
        B[1].debugPrint( "B[1]" );
        B[2].debugPrint( "B[2]" );
        printf( "\n" );
        printf( "  A_A = %14.11f\n", A_A );
        printf( "  B_B = %14.11f\n", B_B );
        A_B.debugPrint( "A_B" );
        B_A.debugPrint( "B_A" );
        A_Z.debugPrint( "A_Z" );
        B_Z.debugPrint( "B_Z" );
        printf( "\n" );
        a.debugPrint( "a" );
        b.debugPrint( "b" );
        printf( "\n" );
        printf( "    M = %14.11f\n", M );
        printf( "  phi = %14.11f\n", phi );
        return 0;
}
//============================================================================

This is the output from the sample program. All the primary values are shown so the calculations can be verified without having to run the program.

        Real             Imaginary        Mag    Angle
        --------------   --------------   ------ -------
 Z[0] = -0.00032563186    0.10802118551   0.1080  1.5738
 Z[1] = -0.07619790924    0.36944527683   0.3772  1.7742
 Z[2] =  0.10202082457   -0.23340312262   0.2547 -1.1587

 A[0] =  0.05987317949    0.10707898936   0.1227  1.0610
 A[1] =  0.14302496043    0.36622285416   0.3932  1.1985
 A[2] = -0.05229282392   -0.23136730416   0.2372 -1.7931

 B[0] =  0.08809308831   -0.03479207270   0.0947 -0.3761
 B[1] =  0.34400795557   -0.11899301853   0.3640 -0.3330
 B[2] = -0.25711837384    0.07517579419   0.2679  2.8571

  A_A =  0.22589139467
  B_B =  0.21323295068
  A_B =  0.00322489968   -0.21793851756   0.2180 -1.5560
  B_A =  0.00322489968    0.21793851756   0.2180  1.5560
  A_Z =  0.18461529779    0.12305734360   0.2219  0.5879
  B_Z = -0.11773875612    0.17987242041   0.2150  2.1504

    a =  0.82533561487    0.00000000001   0.8253  0.0000
    b = -0.56464247339    0.00000000005   0.5646  3.1416

    M =  0.99999999996
  phi =  0.60000000002

Notice that the results for the amplitude and phase (1.0 and .6) have small errors in them. These are from the truncated input values of $ Z $, not any inexactness of the formulas. In the presence of noise or other tones, the results will be good estimates.

The DFT of a pure real tone can be used to determine the exact amplitude and phase of the signal even if the sampling frame is not lined up on a whole number of cycles. Just like the frequency formula in the previous blog, these equations are exact.  

[Edit 2015-05-29: Changed 'Magnitude' to 'Amplitude']

[Edit 2015-06-01: Swap n and k in subscripts]

Sign in