DSPRelated.com
Forums

Calibrating FFT results, amplitude in to magnitude out

Started by Brian Willoughby March 25, 2011
Hello all,

There is one factor that seems to be missing from texts which describe 
the DFT/FFT (or perhaps I have missed it).  That is: The correlation 
between time domain signal amplitude and frequency domain bin magnitudes.

On the one hand, many DSP libraries are very meticulous about 
documenting the differences between their FFT or IFFT implementation 
results versus Matlab.  For example, one implementation of an N-point 
FFT might require that the results be divided by N to correlate with the 
Matlab results, while another implementation would produce the exact 
results without further calculations.  Even though I don't use Matlab, 
these sorts of issues make sense to me as an expected side effect of 
certain code optimizations.

What I haven't found in the general texts is mention of the values I 
should expect if I plug sin() at a constant frequency into a time domain 
array and calculate the FFT.  Should I expect one frequency magnitude 
bin to have a value of 1.0 because the sin() function varies between 
-1.0 and +1.0? ... or would the magnitude bin be expected to hold a 
value of 2.0 because the peak-to-peak amplitude of sin() is actually 2.0 
(+1.0 - -1.0)?  I am asking specifically about sin() frequencies which 
fall precisely into a frequency bin based on the FFT size.  The general 
case would obviously be more complicated due to potential frequency 
smearing.

The reason I ask is that Apple's Accelerate/vecLib Framework seems to be 
returning 2.0 when I carefully construct time domain test signals with 
'unity' amplitude, and for some reason I was expecting 1.0 ... so what's 
the 'right' answer?

By the way, my test code is taking the additional step of converting the 
complex FFT results from their real and imaginary components into 
magnitude (and phase - which should be optional).



I also have a slightly related question: I assume that if the real and 
imaginary components never exceed +/-1.0, then the magnitude could not 
possibly be greater than 1.4142, but I also have a feeling that the 
nature of the FFT is that the results should probably never produce a 
magnitude greater than 1.0, because the real and imaginary components 
are not completely random - i.e. a given pair of real and imaginary 
components would never both be +1.0 or -1.0 unless something terribly 
unnatural had occurred.  Again, assuming that the time-domain signal 
does not exceed +/-1.0 (and I realize that is not always the case with 
uncontrolled inputs).


Brian Willoughby
Sound Consulting

P.S.  I'm sure that I'll happen across one of my texts that does spell 
this out, now that I've gone to the trouble of writing the question. 
But I do not that I've gone looking for a reference on more than one 
occasion and came up with nothing.
When I was introduced to the continuous Fourier transform and its inverse, there was a 1/sqrt(N) in front of each integral. Later, I ran across others in which one integral was left "raw", and the other divided by N. Dividing either the direct or inverse transform works out to be simplest for a particular application, but it seems to me that dividing each by the square root has only symmetry and aesthetics to recommend it.

Personally, I usually prefer the form where the frequency-domain coefficients coincide with the peak amplitudes of the components. That is, for a square wave of amplitude +/-1 and frequency w=2*pi*f, the components are
   (4/pi)*[sin(w) + sin(3*w)/3 + sin(5w)/5 + ...].

Needless to say, this applies also to DFT summations.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
On Mar 25, 9:50�am, Brian Willoughby
<Sound_Consulting-...@Sounds.wa.com> wrote:
> Hello all, > > There is one factor that seems to be missing from texts which describe > the DFT/FFT (or perhaps I have missed it). &#4294967295;That is: The correlation > between time domain signal amplitude and frequency domain bin magnitudes. > > On the one hand, many DSP libraries are very meticulous about > documenting the differences between their FFT or IFFT implementation > results versus Matlab. &#4294967295;For example, one implementation of an N-point > FFT might require that the results be divided by N to correlate with the > Matlab results, while another implementation would produce the exact > results without further calculations. &#4294967295;Even though I don't use Matlab, > these sorts of issues make sense to me as an expected side effect of > certain code optimizations. > > What I haven't found in the general texts is mention of the values I > should expect if I plug sin() at a constant frequency into a time domain > array and calculate the FFT. &#4294967295;Should I expect one frequency magnitude > bin to have a value of 1.0 because the sin() function varies between > -1.0 and +1.0? ... or would the magnitude bin be expected to hold a > value of 2.0 because the peak-to-peak amplitude of sin() is actually 2.0 > (+1.0 - -1.0)? &#4294967295;I am asking specifically about sin() frequencies which > fall precisely into a frequency bin based on the FFT size. &#4294967295;The general > case would obviously be more complicated due to potential frequency > smearing. > > The reason I ask is that Apple's Accelerate/vecLib Framework seems to be > returning 2.0 when I carefully construct time domain test signals with > 'unity' amplitude, and for some reason I was expecting 1.0 ... so what's > the 'right' answer? > > By the way, my test code is taking the additional step of converting the > complex FFT results from their real and imaginary components into > magnitude (and phase - which should be optional). > > I also have a slightly related question: I assume that if the real and > imaginary components never exceed +/-1.0, then the magnitude could not > possibly be greater than 1.4142, but I also have a feeling that the > nature of the FFT is that the results should probably never produce a > magnitude greater than 1.0, because the real and imaginary components > are not completely random - i.e. a given pair of real and imaginary > components would never both be +1.0 or -1.0 unless something terribly > unnatural had occurred. &#4294967295;Again, assuming that the time-domain signal > does not exceed +/-1.0 (and I realize that is not always the case with > uncontrolled inputs). > > Brian Willoughby > Sound Consulting > > P.S. &#4294967295;I'm sure that I'll happen across one of my texts that does spell > this out, now that I've gone to the trouble of writing the question. > But I do not that I've gone looking for a reference on more than one > occasion and came up with nothing.
Hello Brian, You've ventured into an area where things are done in more than one way. Most commonly done is the DFT(FFT) has a gain of N/2 where is is the "size" of the DFT(FFT). This gain results when no scale factor is applied to the direct DFT(FFT) process. Then the inverse DFT(FFT) will end up being multiplied by 2/N to make the round trip gain equal to one. But I know of at least one case where both the forward and inverse DFTs(FFTs) are multiplied by sqrt(2/N) to impart a sort of symmetry to the computational process. A simple test of transforming a sinusoid whose frequency is one of the bin frequencies (avoids leakage) and seeing what the gain is will tell you what gain your process has. But from the defn of a DFT and applying it to a sinusoid on bin frequency will show you where the gain of N/2 comes from. IHTH, Clay
Whoops! Those multipliers are all wrong. Here's the fix:

When I was introduced to the continuous Fourier transform and its inverse, there was a sqrt(2*N) in front of each integral. Later, I ran across others in which one integral was left "raw", and the other multiplied by 2*N. Dividing either the direct or inverse transform works out to be simplest for a particular application, but it seems to me that dividing each by the square root has only symmetry and aesthetics to recommend it.

Some days, it doesn't pay to get out of bed.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
On 2011/03/25 07:26, Clay wrote:
> You've ventured into an area where things are done in more than one > way. Most commonly done is the DFT(FFT) has a gain of N/2 where is is > the "size" of the DFT(FFT). This gain results when no scale factor is > applied to the direct DFT(FFT) process. Then the inverse DFT(FFT) will > end up being multiplied by 2/N to make the round trip gain equal to > one. But I know of at least one case where both the forward and > inverse DFTs(FFTs) are multiplied by sqrt(2/N) to impart a sort of > symmetry to the computational process. A simple test of transforming a > sinusoid whose frequency is one of the bin frequencies (avoids > leakage) and seeing what the gain is will tell you what gain your > process has. > > > But from the defn of a DFT and applying it to a sinusoid on bin > frequency will show you where the gain of N/2 comes from.
Is that what Matlab produces? I must admit that I haven't been able to justify the cost of a Matlab license, and thus I've never used it. Yet Apple and Texas Instruments seem to define their FFT libraries in terms of how Matlab does things. I always assumed that Matlab was producing the results as defined by the DFT, regardless of whether it was 'efficient' for the given processor, where Apple and TI are interested in implementations which take the fewest cycles, and thus may end up with a different scaling. Yours comments leave me thinking that a sinusoid with an amplitude of 1 (or a range of +/-1) should have a different bin magnitude for each FFT size. P.S. Another related question is what to expect from impulse response and step response inputs. I'm inclined to synthesize an impulse as a single +1.0 sample followed by 0.0 samples to pad out the FFT size. Similarly, a step input would be a 0.0 sample followed by +1.0 samples to pad out the FFT size. In both cases, I would hope that the output is limited to +1.0 in each bin magnitude, unless the system under test has overshoot, ringing, resonance, or something like that. e.g., I should be able to test the frequency response of an unknown process by sending an impulse to the input and then calculating an FFT on the output. A 'flat' frequency response would then be expected as bins with +1.0 at every frequency. Have I got the calibration right? Brian Willoughby Sound Consulting
On Mar 27, 9:29&#4294967295;am, Brian Willoughby
<Sound_Consulting-...@Sounds.wa.com> wrote:

...
> &#4294967295;Another related question is what to expect from impulse response > and step response inputs. &#4294967295;I'm inclined to synthesize an impulse as a > single +1.0 sample followed by 0.0 samples to pad out the FFT size. > Similarly, a step input would be a 0.0 sample followed by +1.0 samples > to pad out the FFT size. &#4294967295;In both cases, I would hope that the output is > limited to +1.0 in each bin magnitude, unless the system under test has > overshoot, ringing, resonance, or something like that. &#4294967295;e.g., I should > be able to test the frequency response of an unknown process by sending > an impulse to the input and then calculating an FFT on the output. &#4294967295;A > 'flat' frequency response would then be expected as &#4294967295;bins with +1.0 at > every frequency. &#4294967295;Have I got the calibration right?
Brian, this is a very similar issue to the one posted by Andy365 in the thread: "Energy in a signal as a function of sampling". there is a well-defined answer for both of you. or maybe a couple of different answers that will get you to the same place. the way i dealt with this when i was in college and using the computer to do FFTs and such on data that (in simulation) came from a known source, was to compare the Discrete Fourier Transform (however it's defined, whether it has a 1, a 1/N, or a sqrt(1/N) in front) to a Riemann Summation of the Fourier Integral. the other way i can think of is to simply define a unit sinusoid, bust it into a positive frequency e^(jwn) and a negative frequency component, e(-jwn), and stuff that into the DFT and see what comes out. if w = 2*pi*k/N for integer k, it will come out real clean and you will know how large the number in the FFT bin is for a unit- amplitude sine. r b-j
On 2011/03/27 20:21, robert bristow-johnson wrote:
> the way i dealt with this when i was in college and using the computer > to do FFTs and such on data that (in simulation) came from a known > source, was to compare the Discrete Fourier Transform (however it's > defined, whether it has a 1, a 1/N, or a sqrt(1/N) in front) to a > Riemann Summation of the Fourier Integral. > > the other way i can think of is to simply define a unit sinusoid, bust > it into a positive frequency e^(jwn) and a negative frequency > component, e(-jwn), and stuff that into the DFT and see what comes > out. if w = 2*pi*k/N for integer k, it will come out real clean and > you will know how large the number in the FFT bin is for a unit- > amplitude sine.
Thanks for the response. I hope I don't seem lazy, but I was hoping to just get the answer. ;-) I don't mind also knowing how to independently derive said answer. In college, I remember (thinking) that I understood all of those formulae with 'e' and found them to be a beautiful expression of the symmetry of nature. These days, I have forgotten quite a bit of that, assuming I really ever had it right. Now I tend to only understand the direct and practical application rather than the theoretical proofs. Actually, I've never been able to just look at a formula like e^(j2&pi;f/T) and 'know' what its integral would be. I tend to depend upon the math types to tell me what it works out to be. That said, if I make an attempt to follow what you wrote above, I find myself asking why you translate a sinusoid into both e^(jwn) and e^(-jwn) ... wouldn't that put double the signal in? Do you really need both the negative and positive frequencies for a real-valued sinusoid input? Why wouldn't a pure real sinusoid simply translate to the positive frequency only? I cracked open Steiglitz' DPS Primer and the first example I came across translated sin(kwt) into e^(jkwt), without the negative frequency term. To be precise, he translates sin(k&pi;t/T) with period of 2T into e^(jk2&pi;t/T) with period of T, which is a little confusing in itself because of the change in period, but there you have it. In any event, lately I've been finding myself wishing that someone could write a book on DSP which uses plain English at a somewhat high level to introduce each concept, even if there is still a followup mathematical proof to back up the claims. But I'd like it organized so that I could easily skip the proof until I need to know it and understand it. Every text I've seen seems to start with the math, but that ends up stretching out the description long enough that it's harder to follow unless you already understand the math. Sort of counterproductive. I'm not saying that some understanding of the math is not useful or even necessary, but I also have a feeling that some of these concepts could easily be described to 90% of the necessary level of detail without too much math, at least by someone gifted in writing (yes, an admittedly rare trait among engineers and scientists, not to mention mathematicians). I think Hal Chamberlin's "Musical Applications of Microprocessors" was the last time I sat down for some light reading in plain English and found myself in an "a ha, this all makes sense on a deeper level, now" kind of moment. Of course, only an engineer would consider that book to be light reading, but I find it a lot lighter than the mathematical texts. Anyway, sorry for the rambling, but my original question remains: Can anyone just spell out the answer in simple terms? Two very helpful individuals have replied in this thread, but basically suggested that I derive my answer by manually calculating the DFT. That's an intriguing challenge, but my initial goal was to find a confirmed reference that I could use to test my code against (and the library code that I am using), and thus I feel like if I derive the answer myself then I would still have some doubt. Actually, I'm refining some old code where I thought I had all of these calibrations figured out, and now I'm having doubts, thus the request for a solid answer. P.S. Who is Riemann? I never heard of him before. Actually, scratch my literal question, because Wikipedia has the answer. They even mention Fourier in the article... Brian Willoughby Sound Consulting
On Mar 25, 3:50&#4294967295;pm, Brian Willoughby
<Sound_Consulting-...@Sounds.wa.com> wrote:

> What I haven't found in the general texts is mention of the values I > should expect if I plug sin() at a constant frequency into a time domain > array and calculate the FFT. &#4294967295;Should I expect one frequency magnitude > bin to have a value of 1.0 because the sin() function varies between > -1.0 and +1.0? ... or would the magnitude bin be expected to hold a > value of 2.0 because the peak-to-peak amplitude of sin() is actually 2.0 > (+1.0 - -1.0)? &#4294967295;I am asking specifically about sin() frequencies which > fall precisely into a frequency bin based on the FFT size. &#4294967295;The general > case would obviously be more complicated due to potential frequency > smearing.
Those are two questions, and thus they have two answers. First, assume the sinusoidal has a frequency that coincides with a DFT coefficient. In that case, it is common to define the DFT such that the DFT coefficient is N times as large as the amplitude of the sinusoidal, where N is the DFT length. There is an 1/N assymetry in the scaling coefficinets of the DFT and IDFT, such that (formally) IDFT(DFT(x[n])) = x[n]. The second case is when the frequency of the sinusoidal does *not* coincide with a DFT coefficnet. In this case all the issues regarding scaling as considered in case 1 still applies, along with the frequency smearing issues. Rune
On 2011/03/28 00:18, Rune Allnor wrote:
> On Mar 25, 3:50 pm, Brian Willoughby > <Sound_Consulting-...@Sounds.wa.com> wrote: > >> What I haven't found in the general texts is mention of the values I >> should expect if I plug sin() at a constant frequency into a time domain >> array and calculate the FFT. Should I expect one frequency magnitude >> bin to have a value of 1.0 because the sin() function varies between >> -1.0 and +1.0? ... or would the magnitude bin be expected to hold a >> value of 2.0 because the peak-to-peak amplitude of sin() is actually 2.0 >> (+1.0 - -1.0)? I am asking specifically about sin() frequencies which >> fall precisely into a frequency bin based on the FFT size. The general >> case would obviously be more complicated due to potential frequency >> smearing. > > Those are two questions, and thus they have two answers. > > First, assume the sinusoidal has a frequency that coincides > with a DFT coefficient. In that case, it is common to define > the DFT such that the DFT coefficient is N times as large as > the amplitude of the sinusoidal, where N is the DFT length. > There is an 1/N assymetry in the scaling coefficinets of the > DFT and IDFT, such that (formally) IDFT(DFT(x[n])) = x[n]. > > The second case is when the frequency of the sinusoidal does > *not* coincide with a DFT coefficnet. In this case all the > issues regarding scaling as considered in case 1 still applies, > along with the frequency smearing issues.
Thanks. I realize that it's basically two questions, and I was attempting to exclude the second question. Getting back to it, I think that part of what I am missing is experience with Matlab. Both Apple and Texas Instruments seem intent to describe their FFT library in terms of how the results compare to Matlab's results. That's where I get a little lost. In your answer, you refer to DFT coefficients. Would that be the same as the FFT bins in the result? You mention that the DFT coefficient is N times as large as the amplitude of the sinusoid, but I would interpret that to mean that if I send in a +/-1.0 sin() at a nice frequency for a 128-point FFT, then I would see +128 in the result bin for that frequency (or maybe +256). What's confusing, though, is that Texas Instruments claims their unscaled FFT produces 1/N results compared to Matlab. i.e. you must multiply the Matlab results by N to compare to the TI DSPLIB results. This seems entirely opposite of what you describe. Meanwhile, TI says that their scaled FFT produces the exact same results as Matlab, which doesn't really tell me anything since I'm unfamiliar with Matlab. The reason I'm confused is that your description implies (to me) that I would need to divide by N to get the calibrated results that I want, and yet TI only mentions multiplying by N or leaving the values alone. Some piece here doesn't jibe with the rest. P.S. I assume that, when the input frequency does not coincide exactly with a DFT coefficient, the magnitude of the results for unaligned signal frequencies will be less than or almost equal to the aligned frequency results. In other words, I assume that the energy of the signal is diminished by being spread out over multiple bins. Brian Willoughby Sound Consulting
On Mar 28, 9:42&#4294967295;am, Brian Willoughby
<Sound_Consulting-...@Sounds.wa.com> wrote:

> You mention that the DFT coefficient is N times as large as the > amplitude of the sinusoid,
Use x[n] = exp(j 2*pi*q/N), q and N integers, q < N/2. Then (view with fixed-width font) N-1 X[k] = sum exp(j2pi q/N) exp(-j2pi k/N) n=0 If q =/= k then X[k] == 0. When q == k then N-1 N-1 X[k] = sum exp(j2pi k/N) exp(-j2pi k/N) = sum 1 = N. n=0 n=0 So the DFT coefficient X[q] has N times as large amplitude as the sinusoidal x[n]. Do note the 1/N factor in the expression for the IDFT. Rune