The DFT Magnitude of a Real-valued Cosine Sequence

Rick LyonsJune 17, 20148 comments

This blog may seem a bit trivial to some readers here but, then again, it might be of some value to DSP beginners. It presents a mathematical proof of what is the magnitude of an N-point discrete Fourier transform (DFT) when the DFT's input is a real-valued sinusoidal sequence.

To be specific, if we perform an N-point DFT on N real-valued time-domain samples of a discrete cosine wave, having exactly integer k cycles over N time samples, the peak magnitude of the cosine wave's positive-frequency spectral component will be

This article is available in PDF format for easy printing

where A is the peak amplitude of the discrete cosine sequence. Graphically, Eq. (1) is shown in Figure 1.

I presented the equivalent of Eq. (1) in the 3rd edition of my DSP textbook. However, I recently received an e-mail from a reader who asked me:

"(I can see that) |X(k)|=AN/2 indeed correct, but
I cannot see mathematically where it comes from!
Such a simple, intuitive thing, but believe it or
not, even after having gone through your book,
I just can't figure out where this is derived from!
Sometimes the obvious things are made difficult by
the reader himself!

After reading that e-mail I thought, "Didn't I give the derivation of Eq. (1) in my book?" The answer is, shame on me, "No." Out of curiosity I checked four other DSP textbooks from my bookshelf and they chose not to provide Eq. (1) at all. How strange! I believe Eq. (1) is important and later in this blog I explain why. In any case, below is the mathematical derivation of Eq. (1).

Deriving Equation (1)

We start with a discrete time-domain cosine sequence x(n) as:

where A is the peak value of the cosine wave, k is the integer number of complete sinusoidal cycles occurring in the N samples, and variable n is the time index. We can show the desired N-point DFT of x(n), X(m), as:

where m is the frequency-domain index.

OK. Remembering Euler's relationship of cos(α) = (e + e-jα)/2, we can write X(m) as:

The left summation in Eq. (4) represents a positive-frequency spectral component, and the right summation represents a negative-frequency spectral component. In deriving Eq. (1) we need concentrate only on the positivefrequency summation which is:

Equation (5) is a geometric series of the form:

where variable q = -j2π(m-k)/N is a simple dummy variable we’ll use to keep the next equation simple. For mk Eq. (6) can be written in closed-form as:

(Equation (7) is a very important equation in the field of DSP, however it's typically presented as a simple unjustified fact in DSP books. If you want to see the full derivation of Eq. (7), see reference [1] or reference [2].)

Substituting q = -j2π(m-k)/N into the right side of Eq. (7), we rewrite Eq. (5) as:

OK, we have what we want. Equation (8) is a closed-form expression for the positive-frequency DFT of a real-valued input cosine sequence. (We could perform the algebraic acrobatics to convert Eq. (8) into a familiar sin(x)/x form, but we need not do that here.)

With the original DFT input being exactly integer k cycles of a cosine sequence, to verify Eq. (1) we evaluate Eq. (8) when the frequency index is m = k. Doing that we run into trouble because spectral sample X(m=k) is:

My original solution to the Eq. (9) problem was to use l'Hôpital's Rule to evaluate Eq. (8) when m = k.

That solution, while valid, was overly complicated and unnecessary. As DSP guru Prof. Dilip Sarwate [3] showed me, a simpler solution is to merely evaluate Eq. (6) for q = 0 (m = k in Eq. (5)) as

That's it. X(k) is real-valued and its magnitude is AN/2 which verifies Eq. (1).

Why Equation (1) is Important

Beyond the facts that knowing Eq. (1) will make you live longer, improve your sex life, and reduce global warming, Eq. (1) is important when computing fast Fourier transforms (FFTs) on fixed-point binary hardware. When you're computing FFT magnitude samples using binary two's-complement number representation, you want to avoid numerical overflow errors. So this means, for example, if you're using 16-bit hardware you want no computed FFT magnitude sample value to be greater than 215.

Exploring this idea, let's assume that the input data to an N-point FFT, computed using B-bit hardware (B = 16 in the above paragraph), is b-bits in word width. That means the maximum positive peak amplitude value of an input cosine wave is Amax = 2b-1. From our Eq. (1) we can write:

However, to avoid overflow using two's-complement arithmetic:

Rearranging Eq. (14) we have:

Sorry for all the algebra rigmarole. What Eq. (14) is telling is, based on Eq. (1), there is a maximum-sized FFT you can perform on fixed-point hardware while avoiding numerical overflow. For example, on B = 16-bit hardware and b = 8-bit input samples, the largest sized FFT you can compute to determine spectral magnitude values while avoid numerical overflow is

That Nmax = 512 value is shown in Figure 2 above the b = 8 x-axis value.

So this means if we perform 1024- or 2014-point FFTs, with a 16-bit machine accepting 8-bit input samples, it's possible that numerical overflow will occur.

Now, during a System Design Review meeting, if someone suggests performing 2048-point FFTs using 16-bit hardware on 8-bit data, you can pipe up with Eq. (16), derived from Eq. (1), to warn your colleagues of potential numerical overflow. And then tell them that some sort of FFT-internal data-checking and conditional truncation, or multi-word processing, will be necessary to perform those 2048-point FFTs. (Think of the increased software development cost. Yikes!) Or else you could just recommend using floating-point hardware.

A Few Final Thoughts

To determine the negative-frequency X(N-k) DFT sample's magnitude, we keep in mind that due to the circular nature of the DFT, X(N-k) = X*(k) = AN/2, where the '*' symbol means complex conjugate. (See reference [4].)

If our N-point DFT's input, in Eq. (3), had been a sine wave sequence,

the above derivation method, using Euler's relationship of sin(α) = (e - e-jα)/j2, would produce the same positive-frequency result of X(k) = AN/2.

One last thought from me, and it's a criticism. I stated that I couldn't find a derivation of Eq. (1) in the DSP books on my bookshelf. Well, the closest match I found was in a popular college DSP textbook. In that book the authors have the homework problem:

"Compute the DFT of x(n) = cos(2πnko/N), where 0 ≤ nN-1"

In their book's Solution Manual (I own a legal copy), the authors' solution consisted of nothing more than the following two lines:

Now you tell me, ...how can anyone expect a first-semester DSP student to go from Eq. (17), by inspection in one step, to Eq. (18)? I found no discussion in that book's DFT chapter that would show the student how to go from Eq. (17) to the delta-function expression in Eq. (18). It's no wonder why DSP can be so painful to learn in college classrooms!

References
[1] R. Lyons, “Understanding Digital Signal Processing”, 3rd Edition, Prentice Hall Publishing, 2011, Appendix B.
[2] http://en.wikipedia.org/wiki/Geometric_series
[3] http://en.wikipedia.org/wiki/L%27H%C3%B4pital%27s_rule
[4] J. Proakis and D. Manolakis, Digital Signal Processing- Principles, Algorithms, and Applications, 3rd Ed., Prentice Hall, 1996, pp. 413.


Previous post by Rick Lyons:
   Specifying the Maximum Amplifier Noise When Driving an ADC
Next post by Rick Lyons:
   Sum of Two Equal-Frequency Sinusoids

Comments:

[ - ]
Comment by Dilip SarwateAugust 19, 2014
Rick:

A little over ten years ago, on comp.dsp you had posted a similar (and completely unnecessary) application of L'Hopital's rule in evaluating the sum (6) via the geometric series formula (7), and I had pointed out that (7) DOES NOT HOLD when m = k giving that q = 0. In this case, the summands in (6) all have value e^{nq} = e^0 = 1, and so, since we are summing N terms, the sum is N.

Please consider editing your otherwise very nice write-up to exclude the long-winded discussion of calculus and L' Hopital's rule.
[ - ]
Comment by MRU13December 4, 2014
Hi. You note that, "We could perform the algebraic acrobatics to convert Eq. (8) into a familiar sin(x)/x form, but we need not do that here". Is the conversion of Eq. (8) into the sin(x)/x available somewhere? I have been unable to do this conversion and unable to find it anywhere. Thanks.
[ - ]
Comment by Rick LyonsMarch 28, 2017

Good Gosh!  MRU13, I just today saw your post! Here's the information you requested:

eq 8_78312.jpg

[ - ]
Comment by Christoph59March 22, 2015
Great post!

"Now you tell me, ...how can anyone expect a first-semester DSP student to go from Eq. (18), by inspection in one step, to Eq. (19)? I found no discussion in that book's DFT chapter that would show the student how to go from Eq. (18) to the delta-function expression in Eq. (19). It's no wonder why DSP can be so painful to learn in college classrooms!"

And in our day and age where you expect every answer to be no more than a click away, you find that the internet is practically empty of explanations. Thanks for this enlightening post!
[ - ]
Comment by Dilip SarwateNovember 13, 2015
"Now you tell me, ...how can anyone expect a first-semester DSP student to go from Eq. (18), by inspection in one step, to Eq. (19)?"

Rick:

Perhaps you did not notice that in your derivation of (1) you actually proved (19)? You showed in (8) that

X(m) = (1/2) [1 - exp(j2 pi (m-k))]/[1 - exp(j2 pi (m-k)/N)]

provided that m was not equal to k. But, in that case, exp(j2 pi (m-k)) = 1 and so the numerator is zero!. That is, you have shown that X(m) = 0 unless m = k (or N-k) and also (later, after much computation that I suggested in an earlier comment was not really necessary) that X(k) = X(N-k) = N/2; and you are done! You have proved (19)!! In passing, I note that while you set out to prove that the MAXIMUM DFT coefficient has value N/2 and occurs at k and N-k, you never really stated (or proved) that the N/2 is indeed the largest value of the DFT coefficient, that is, that the magnitude of (8) is always smaller than N/2. But just EVALUATING (8) instead of suggesting that it could be transformed to sin x / x form does the trick.

[ - ]
Comment by CedronMarch 28, 2017
Hi Rick,

I just found this blog article and have a few comments.

When you say "we need concentrate only on the positivefrequency summation" and then proceed to discard the negative frequency component, you are actually transforming the analysis from being on a real valued pure tone to a complex valued pure tone.  That's not a big deal and rather subtle, but I think you should have said that the negative frequency (N-m) bin will be zero valued because for integer frequencies of complex tones all bins except the bin matching the frequency are zero.  You have the equation in (8) to show this, but you never make it clear.

Notice that in (8), since $m$ is an integer that $ e^{-j2 \pi m} = 1 $ and therefore the numerator can be simplified to $ 1 - e^{j2 \pi k}  $.  This simplification allows for the simple derivation of an exact frequency formula for complex tones.

Your statement "We could perform the algebraic acrobatics to convert Eq. (8) into a familiar sin(x)/x form, ..." is simply not true.  MRU13 should not feel bad.


For readers who want further reading on this topic, I have several blog articles on this site that cover the same material.

For a derivation of the geometric summation formula and an alternative understanding of the meaning of a DFT:

DFT Graphical Interpretation: Centroids of Weighted Roots of Unity
https://www.dsprelated.com/showarticle/768.php

For a more comprehensive analysis of the material presented in this blog article (up to the overflow discussion):

DFT Bin Value Formulas for Pure Real Tones
https://www.dsprelated.com/showarticle/771.php

For the corresponding treatment of a complex tone, including the relationship between the sinc function (sin(x)/x) and the Dirichlet kernel:

DFT Bin Value Formulas for Pure Complex Tones
https://www.dsprelated.com/showarticle/1038.php

For the derivation of the exact frequency formula I mentioned previously:

A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT
https://www.dsprelated.com/showarticle/1039.php

Ced
[ - ]
Comment by Rick LyonsMarch 28, 2017

Cedron, you wrote:

   "Your statement "We could perform the algebraic acrobatics

    to convert Eq. (8) into a familiar sin(x)/x form, ..." is

    simply not true."

Why do you think my statement is not true?

[ - ]
Comment by CedronMarch 28, 2017

Hi Rick,

Same reason you put a squiggly equal sign in your equation (8b).  When you say something "is" when it should truly be "is like" you are making an incorrect assertion.  Converting an equation from one form to another is quite different than finding an approximation for it.  You could say an infinite number of decimal places worth.  To not make the distinction is to be misleading and simply not true.

In a recent thread on comp.dsp I pinpoint to the day, less than a year ago, when E.J. learned that the Sinc function is an approximation in this case, so this misconception is quite prevalent.

The Sinc function (sin(x)/x) pertains to the continuous FT.  The Dirichlet kernel (sin(Pi*d)/sin(Pi*d/N)) pertains to the discrete (normalization factor = 1) FT.  You don't make learning DSP any easier by conflating the two.

To be clear, this pertains to complex tones, not real valued sinusoids.  In this case you are talking about the positive frequency component which is a complex tone.  

The real tone case is a lot more complicated.  There is a version of my bin value formulas that has the Dirichlet kernel as a factor.  Therefore it is also a factor in the magnitude.  A blog article about this is in my to do list.

If anybody wants to figure it out on their own, start with equation (25) in "DFT Bin Value Formulas for Pure Real Tones" (linked above) and recognize that the difference of two Cosines can be expressed as the product of two Sines. 

Ced

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.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.

Sign up
or Sign in