The DFT Magnitude of a Real-valued Cosine Sequence
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
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(α) = (ejα + 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 m ≠ k 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(α) = (ejα - 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 ≤ n ≤ N-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.
- Comments
- Write a Comment Select to add a comment
Good Gosh! MRU13, I just today saw your post! Here's the information you requested:
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.
"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!
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.
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 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
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?
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
Nice post! I would like to add a recent finding of mine which appears to "have been forgotten" in the DSP domain: The inverse DFT models the periodic continuation of a sequence x(n) for which you computed the X(k) as a sum of cosines. The amplitudes of the cosines are 2|X(k)| and the phases are exactly the phase angles of the X(k). It is straight forward to show this for a real sequence.
Then it follows automatically that |X(k)| = A N/2 is true.
Hello DFT_Enthusiast.
Your |X(k)| = A*N/2 equation is correct under the conditions that the sine wave corresponding to freq index 'k' [1] has a peak amplitude of A and [2] contains an integer number of cycles over N time-domain samples.
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.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: