DFT Bin Value Formulas for Pure Complex Tones

March 17, 2017

Introduction

This is an article to hopefully give a better understanding to the Discrete Fourier Transform (DFT) by deriving an analytical formula for the DFT of pure complex tones and an alternative variation. It is basically a parallel treatment to the real case given in DFT Bin Value Formulas for Pure Real Tones. In order to understand how a multiple tone signal acts in a DFT it is necessary to first understand how a single pure tone acts. Since a DFT is a linear transform, the DFT of a signal consisting of a sum of pure tones will be the sum of the DFTs of the constituent tones.

A Pure Complex Tone Sequence

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$. For complex tones negative frequencies make a lot more sense 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 can be rescaled to a distance along the circumference which is also an angular scale. Defining an angle variable makes the definition simpler. $$\alpha = f \cdot { 2\pi \over N } \tag 2$$ Substituting in this definition in (1) gives a simpler formula for the signal. $$S_n = M \cdot { e^{ i(\alpha n + \phi) } } \tag 3$$

The Discrete Fourier Transform

This is the definition of DFT. The same one used in the real case article: $$Z_k = { 1 \over N } \sum_{n=0}^{N-1} { S_n e^{ -i2 \pi { k \over N } n } } \tag 4$$ The summation can be made more comprehensible by making a frequency index to angular measurement substitution: $$\beta_k = k \cdot { 2\pi \over N } \tag 5$$

Making the substitution results in the following version of the DFT calculation: $$Z_k = { 1 \over N } \sum_{n=0}^{N-1} { S_n e^{ -i\beta_k n } } \tag 6$$

Roots of Unity

The $\beta_k$s mark the locations along the circumference of a set of N Roots of Unity. $$R_k = e^{ -i \beta_k } \tag {7}$$ The $R_k$s are the complex number values at the Root of Unity locations. Because of the negative sign in the exponent, the $\beta_k$s are indexed in a counter-clockwise manner starting with index zero at 1 (also known as Unity), and the Roots of Unity are indexed clockwise, also with index zero at 1.

Showing the each $R_k$ is an Nth Root of Unity is straightforward: $$R_k^N = (e^{ -i \beta_k })^N = (e^{ -i (k \cdot { 2\pi \over N }) })^N = (e^{ -i 2\pi })^k = 1^k = 1 \tag {8}$$

The DFT of a Pure Complex Tone

The definition of the pure complex tone sequence can now be plugged into the DFT definition to obtain a formula for the DFT values. $$Z_k = { 1 \over N } \sum_{n=0}^{N-1} { M \cdot { e^{ i(\alpha n + \phi) } } \cdot e^{ -i\beta_k n } } \tag {9}$$ This form isn't that useful. It can be made a little better by factoring out the constants. $$Z_k = { M \over { N } } e^{ i\phi } \sum_{n=0}^{N-1} { e^{ i(\alpha - \beta_k ) n } } \tag {10}$$

Applying the Summation Formula

Unlike the real signal case, the complex case consists of just one geometric summation. $$\sum_{n=0}^{N-1} {x^n} = { { 1 - x^N } \over { 1 - x } } \tag {11}$$ Applying the summation formula to the DFT definition results in this equation: $$Z_k = { M \over N } e^{ i \phi } \cdot { 1 - e^{ i(\alpha - \beta_k ) N } \over 1 - e^{ i(\alpha - \beta_k ) } } \tag {12}$$ It is useful to note, like in the real case, that $\beta_k$ times $N$ is a whole number multiple of $2 \pi$. $$\beta_k N = k \cdot 2\pi \tag {13}$$ This allows the $\beta_k$ to be removed from the numerator. $$Z_k = { M \over N } e^{ i \phi } \cdot { 1 - e^{ i\alpha N } \over 1 - e^{ i(\alpha - \beta_k ) } } \tag {14}$$ At this point the summation has been eliminated and a closed form analytic formula remains. The complex tone version of this formula is much simpler than the real signal version. The most notable difference is the role of the phase variable. Shifting the phase in a complex tone in the time domain causes the entire DFT bin set in the frequency domain to rotate in parallel. This does not happen in the real signal case.

Handling the Integer Frequency Case

When f is an integer frequency the value $\alpha N$ becomes a whole number multiple of $2 \pi$. This means the exponent in the numerator of (14) become 1 resulting in the numerator equalling zero. $$\alpha N = f \cdot 2\pi \tag {15}$$ For all cases of $k$, except when $k=f$, the result is that the bin value is zero. When $k=f$ the equation becomes 0 divided by 0 which is an indeterminate form. This case can then be evaluated by going back to (10). The summation becomes a sum of ones. There are $N$ of them so the sum is $N$. This $N$ will cancel out the $N$ in the denominator in the outside factor leaving a very simple equation: $$Z_k = M e^{ i \phi } \, \text{ when } \, k=f \tag {16}$$

Qualitative Analysis

The behavior of the bin values can now be understood from the equation. First off, as just explained, when the frequency is an integer value all the bin values will be zero except for when $k = f$. The answer is then (16) which is well known and well understood.

When the frequency is close to an integer value, $e^{ i\alpha N }$ will be near one and the numerator will be small. So all the bin values will be small except the case of the closest $k$ to $f$. In this case the denominator will also be small and the quotient will be larger. As the $k$s get further from $f$, the denominator gets larger and the bin values taper off.

When the $k$ value goes from below $f$ to above $f$, the sign of the denominator changes and the numerator stays the same. This explains why the bin values on either side of $f$ are nearly opposite each other.

When the frequency is further from being an integer value, the numerator becomes relatively larger so all the bin values are also larger. Since the denominator doesn't get nearly as small compared to the near integer case, the peak value is much smaller compared to the neighboring bins.

The amplitude is in the numerator of the outside factor so the bin values are directly proportional to it. This follows from the DFT being a linear operation.

Signal and Roots of Unity Formulation

As a matter of being interesting, rather than any practicality, the bin value formula can be expressed exclusively in terms of some signal values and the Roots of Unity.

Dividing the first signal point by the second give the value of $e^{ i\alpha }$. Any two sample points in a row would give the same results. $${ S_1 \over S_0 } = { Me^{ i(\alpha + \phi) } \over Me^{ i\phi } } = e^{ i\alpha } \tag {17}$$ Next, multiplying the amplitude and phase term across the numerator turns those terms into the definition of the first signal point and the first point just past the DFT frame. $$Z_k = { 1 \over N } \cdot { Me^{ i \phi } - Me^{ i(\alpha N + \phi) } \over 1 - e^{ i\alpha } e^{ -i\beta_k } } \tag {18}$$ Substitute in (17) into the denominator and the equation is complete. The bin values are now expressed as a transform of the roots of unity using only signal values as parameters. $$Z_k = { 1 \over N } \cdot { S_0 - S_N \over 1 - { S_1 \over S_0 } R_k } \tag {19}$$ Since the complex tone case is simpler than the real tone case, this equation is simpler than its real counterpart.

An Alternative Form

The traditional approach to studying the DFT treats it as a sampled case of the continuous case. In this regard, the DFT is thought of being points along an underlying continuous function. The relative scale has its zero at the frequency of the tone. Within the context of this article, it is as if the index values of k no longer have to be integers but can take on any real value on the same scale. $$\delta = \beta_k - \alpha = ( k - f ) { 2 \pi \over N } \tag {20}$$ Since the assumption that k is an integer has been relaxed, the simplification above that $\beta_k$ is a multiple of $2 \pi$ is no longer valid and (12) will be used as the base formula. Substituting in the definition of $\delta$ yields this equation: $$Z_k = { M \over N } e^{ i \phi } \cdot { 1 - e^{ -i \delta N } \over 1 - e^{ -i\delta } } \tag {21}$$ Factoring out the appropriate exponential values gets the formula ready for a trigonometric substitution. $$Z_k = { M \over N } e^{ i \phi } \cdot { e^{ -i \delta N/2 } \over e^{ -i\delta/2 } } \cdot { e^{ i \delta N/2 } - e^{ -i \delta N/2 } \over e^{ i\delta/2 } - e^{ -i\delta/2 } } \tag {22}$$ This is the definition of the Sine function in exponential terms. It follows straight from Euler's equation which you read about in my first blog article titled " The Exponential Nature of the Complex Unit Circle". $$sin( \theta ) = { e^{ i\theta } - e^{ -i\theta } \over 2i } \tag {23}$$ Dividing the numerator and denominator of the rightmost fraction in (22) by $2i$ makes the substitution possible. The other exponential terms can be combined into one. $$Z_k = { M \over N } e^{ i \left[ -\delta (N-1) / 2 + \phi \right] } \cdot { sin( \delta N / 2 ) \over sin( \delta / 2 ) } \tag {24}$$ With the exception of the phase value being included, this equation is also known as the Dirichlet Kernel[1].

The Magnitude

The bin value is a complex number. The magnitude of a complex number can be found by taking the square root of its product with its conjugate. $$\| Z_k \| = ( Z_k \cdot \overline{ Z_k } ) ^ { 1 \over 2 } \tag {25}$$ Finding the conjugate in this case means flipping the sign in the exponential because the fractional portion is real. $$\overline{ Z_k } = { M \over N } e^{ -i \left[ -\delta (N-1) / 2 + \phi \right] } \cdot { sin( \delta N / 2 ) \over sin( \delta / 2 ) } \tag {26}$$ Evaluating the magnitude eliminates the phase variable completely. This is consistent with the observation that shifting the phase rotates the DFT bin value but doesn't change its magnitude. $$\| Z_k \| = { M \over N } \left\lvert { sin( \delta N / 2 ) \over sin( \delta / 2 ) } \right\rvert \tag {27}$$ The absolute value bars are necessary because the Sines can sometimes take on negative values.

As N Gets Large

A comparison can be made between the continuous case of the FT and the DFT by studying what happens as the number of samples, and the size of the DFT, gets large for a given analysis interval. The first step is to change the scale back from being angular to being bin based. From (17) we can define the difference between the bin index and the frequency of the complex tone as $x$ and relate it to $\delta$. $$x = k - f = \delta \cdot { N \over 2 \pi } \tag {28}$$ Substitute this into (27). $$\| Z_k \| = { M \over N } \left\lvert { sin( x \pi ) \over sin( x \pi /N ) } \right\rvert \tag {29}$$ Bring N inside the absolute value bars and divide the numerator and denominator by $x \pi$. $$\| Z_k \| = M \left\lvert { { sin( x \pi ) \over x \pi } \over { sin( x \pi/N ) \over x \pi/N } } \right\rvert \tag {30}$$ The denominator now takes the form of this well known limit: $$\lim_{\theta \to 0} { sin( \theta ) \over \theta } = 1 \tag {31}$$ Explicitly it looks like this: $$\lim_{ N \to \infty } { { sin( x \pi/N ) \over x \pi/N } } = 1 \tag {32}$$ This simplifies (30). $$\lim_{ N \to \infty } \| Z_k \| = M \left\lvert { sin( x \pi ) \over x \pi } \right\rvert \tag {33}$$ The expression inside the absolute value bars is known as the Normalized Sinc function[2]. Since the limit operations were all done within the absolute value bars, this equation shows that the Normalized Sinc function is the limit of the Dirichlet Kernel function. The Sinc function pertains to the continuous FT and the Dirichlet pertains to the discrete FT. They are not the same, but they are approximations of each other.

Quantitative Example

Here is a test program to demonstrate that (14) is correct and a comparison between the true Dirichlet Kernel values and the Normalized Sinc funtion values. It is designed purposefully to demonstrate the equations, not for optimal performance.

#include <math.h>
#include <stdio.h>
struct complex
{
double Real;
double Imag;
};
//============================================================================
void DFT( double beta_k, int N, double M, double phi, double alpha, complex* result )
{
result->Real = 0.0;
result->Imag = 0.0;
for( int n = 0; n < N; n++ )
{
double theAngle = alpha * (double) n + phi;
double s_Real = M * cos( theAngle );
double s_Imag = M * sin( theAngle );
theAngle = beta_k * (double) n;
double b_Real =  cos( theAngle );
double b_Imag = -sin( theAngle );
result->Real += s_Real * b_Real - s_Imag * b_Imag;
result->Imag += s_Real * b_Imag + s_Imag * b_Real;
}
result->Real /= (double) N;
result->Imag /= (double) N;
}
//============================================================================
int main( int argc, char *argv[] )
{
double freq =  5.5;
double phi  =  1.0;
double M    =  1.0;
int    N    = 16;
double alpha = freq * 2.0 * M_PI / (double) N;
complex bin;
printf( "           DFT               Formula\n" );
printf( " k    Real     Imag       Real     Imag\n" );
printf( "--  -------- --------   --------  -------\n" );
double q = M / (double) N;
double a = alpha * (double) N + phi;
double num_real = q * ( cos( phi ) - cos( a ) );
double num_imag = q * ( sin( phi ) - sin( a ) );
for( int k = 0; k < N; k++ )
{
double beta_k = (double) k * 2.0 * M_PI / (double) N;
DFT( beta_k, N, M, phi, alpha, &bin );
double d = alpha - beta_k;
double den_real = 1 - cos( d );
double den_imag =    -sin( d );
double den2 = den_real * den_real + den_imag * den_imag;
double quot_real = ( num_real * den_real
+ num_imag * den_imag ) / den2;
double quot_imag = ( num_imag * den_real
- num_real * den_imag ) / den2;
printf( "%2d  %8.5f %8.5f   %8.5f %8.5f\n",
k,
bin.Real,
bin.Imag,
quot_real,
quot_imag );
}
printf( "\n\n" );
printf( " k    Real     Imag      Mag       True     Sinc     x\n" );
printf( "--  -------- --------  -------   -------  -------  ----\n" );
for( int k = 0; k < N; k++ )
{
double beta_k = (double) k * 2.0 * M_PI / (double) N;
DFT( beta_k, N, M, phi, alpha, &bin );
double v = sqrt( bin.Real * bin.Real + bin.Imag * bin.Imag );
double x = k - freq;
double s = sin( M_PI * x ) / ( M_PI * x );
double t = sin( M_PI * x ) / ( N * sin( M_PI * x / N ) );
printf( "%2d  %8.5f %8.5f %8.5f  %8.5f %8.5f  %4.1f\n",
k,
bin.Real,
bin.Imag,
v,
fabs( t ),
fabs( s ), x );
}
return 0;
}
//============================================================================

This is the output:
           DFT               Formula
k    Real     Imag       Real     Imag
--  -------- --------   --------  -------
0   0.00566  0.07064    0.00566  0.07064
1  -0.00939  0.08031   -0.00939  0.08031
2  -0.03031  0.09374   -0.03031  0.09374
3  -0.06462  0.11577   -0.06462  0.11577
4  -0.13960  0.16391   -0.13960  0.16391
5  -0.50021  0.39545   -0.50021  0.39545
6   0.56774 -0.29027    0.56774 -0.29027
7   0.20714 -0.05873    0.20714 -0.05873
8   0.13216 -0.01059    0.13216 -0.01059
9   0.09785  0.01144    0.09785  0.01144
10   0.07693  0.02488    0.07693  0.02488
11   0.06188  0.03454    0.06188  0.03454
12   0.04972  0.04235    0.04972  0.04235
13   0.03895  0.04927    0.03895  0.04927
14   0.02859  0.05592    0.02859  0.05592
15   0.01782  0.06284    0.01782  0.06284
k    Real     Imag      Mag       True     Sinc     x
--  -------- --------  -------   -------  -------  ----
0   0.00566  0.07064  0.07087   0.07087  0.05787  -5.5
1  -0.00939  0.08031  0.08085   0.08085  0.07074  -4.5
2  -0.03031  0.09374  0.09852   0.09852  0.09095  -3.5
3  -0.06462  0.11577  0.13258   0.13258  0.12732  -2.5
4  -0.13960  0.16391  0.21531   0.21531  0.21221  -1.5
5  -0.50021  0.39545  0.63764   0.63764  0.63662  -0.5
6   0.56774 -0.29027  0.63764   0.63764  0.63662   0.5
7   0.20714 -0.05873  0.21531   0.21531  0.21221   1.5
8   0.13216 -0.01059  0.13258   0.13258  0.12732   2.5
9   0.09785  0.01144  0.09852   0.09852  0.09095   3.5
10   0.07693  0.02488  0.08085   0.08085  0.07074   4.5
11   0.06188  0.03454  0.07087   0.07087  0.05787   5.5
12   0.04972  0.04235  0.06531   0.06531  0.04897   6.5
13   0.03895  0.04927  0.06280   0.06280  0.04244   7.5
14   0.02859  0.05592  0.06280   0.06280  0.03745   8.5
15   0.01782  0.06284  0.06531   0.06531  0.03351   9.5

The code can easily be copied and pasted into a local file for testing.

Conclusion

Just like the real tone case, the DFT bin values of a complex pure tone can be calculated without doing a summation. Having a mathematical statement of this value gives both a qualitative understanding to the behavior of a pure complex tone in the DFT as well as a quantitative method for efficient computation of the same. This article has given two different functions in several different forms which both give exact values for the DFT bin values at integer bin numbers, but differ in the values in between. These equations can also be used as part of analyzing signals to find frequency, amplitude, and phase of constituent tones which will be covered by future blog articles.

References

Previous post by Cedron Dawg:
Exponential Smoothing with a Wrinkle
Next post by Cedron Dawg:
A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT

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.