Forums

FFTW and Correct Result Interpretation

Started by jdgonko April 14, 2006
I'm using FFTW 3.1.1 and VC++ Pro 2005.  The input signal is a 27Hz,
16-bit, 2 channel, Fs=44100, fullscale sinewave (using -32768 to 32767 in
amplitude which is the size of type-short).  No window function used (a
blackman is selectable however it leaks outside of the 27Hz bin).  Prior
to the fftw_execute, I convert all the samples to an amplitude between
-1.0..+1.0 by dividing by 32767 and storing the type-double results in a
new array.  Type double is perfect for the fftw_execute input arrays (left
and right).  The last data input read from file is padded with zeroes to
make it a full 44,100 samples.

The resulting data is output to a comma-delimited file.  First column is
the FFT bin (which exactly corresponds to the frequency as the FFT is the
same as the Fs of 44100 so I get 0 to 22,050 bins where 0 = DC).  2nd
column is the average magnitude from all the FFT iterations (i.e. add-up
all the abs(mag) values from Re[x] and divide by the iterations needed to
process the file).  3rd column is the largest abs(mag) found during all
the FFT iterations.  The 4th and 5th columns are for the right channel
(same as the 2nd and 3rd columns for the left channel).

From 1 to 26Hz, the largest magnitude is around 0.00105 (or less).  Then,
the 84.77523 spike right at 27Hz.  Then around 0.00105 again for bins 28
through 22,050.  Seems to work great!

But, how do I convert these results to dB?  The 16-bit input data should
provide for a mathematical "-dB floor" value right?

And, I'm only using the FFT results from the Re[x] type-double memory cell
to learn everything on the magnitudes.  I'm ignoring the resulting Im[x]
cell values.  Is that correct?

Thanks!
Joe



Ignoring the imaginary values is not correct, for a given complex
frequency, a+bi, the  magnitude is equal a^2 + b^2.  Converting it to
dB is then just 10 * log10(a^2 + b^2).

If you want real power levels, you will have to figure out how your
count values should be converted to volts and how the sensitivity of
your sensor converts what you are trying to measure into volts.

The noise floor will be limited by the resolution of your A/D, probally
electrical noise will make it higher than that.

Thank you very much for providing the assistance!  I'll plug-in the formula
and give it some data to crunch.

Have a great weekend.
Joe

>Ignoring the imaginary values is not correct, for a given complex >frequency, a+bi, the magnitude is equal a^2 + b^2. Converting it to >dB is then just 10 * log10(a^2 + b^2). > >If you want real power levels, you will have to figure out how your >count values should be converted to volts and how the sensitivity of >your sensor converts what you are trying to measure into volts. > >The noise floor will be limited by the resolution of your A/D, probally >electrical noise will make it higher than that. > >
"jdgonko" <josef.gonko@co.anoka.mn.us> wrote in message 
news:YMOdneNz0LJMFaLZnZ2dnUVZ_tednZ2d@giganews.com...
> I'm using FFTW 3.1.1 and VC++ Pro 2005. The input signal is a 27Hz, > 16-bit, 2 channel, Fs=44100, fullscale sinewave (using -32768 to 32767 in > amplitude which is the size of type-short). No window function used (a > blackman is selectable however it leaks outside of the 27Hz bin). Prior > to the fftw_execute, I convert all the samples to an amplitude between > -1.0..+1.0 by dividing by 32767 and storing the type-double results in a > new array. Type double is perfect for the fftw_execute input arrays (left > and right). The last data input read from file is padded with zeroes to > make it a full 44,100 samples. > > The resulting data is output to a comma-delimited file. First column is > the FFT bin (which exactly corresponds to the frequency as the FFT is the > same as the Fs of 44100 so I get 0 to 22,050 bins where 0 = DC). 2nd > column is the average magnitude from all the FFT iterations (i.e. add-up > all the abs(mag) values from Re[x] and divide by the iterations needed to > process the file). 3rd column is the largest abs(mag) found during all > the FFT iterations. The 4th and 5th columns are for the right channel > (same as the 2nd and 3rd columns for the left channel). > > From 1 to 26Hz, the largest magnitude is around 0.00105 (or less). Then, > the 84.77523 spike right at 27Hz. Then around 0.00105 again for bins 28 > through 22,050. Seems to work great! > > But, how do I convert these results to dB? The 16-bit input data should > provide for a mathematical "-dB floor" value right? > > And, I'm only using the FFT results from the Re[x] type-double memory cell > to learn everything on the magnitudes. I'm ignoring the resulting Im[x] > cell values. Is that correct? > > Thanks! > Joe
There is going to be leakage in the process anyway. The only thing that stops you from seeing leakage, and therefore assuming there is none?, is that you select a time epoch that covers an integer number of cycles of the particular signal (27Hz) you're analyzing right now. This is a somewhat interesting case in that you are working with a sinusoid that has a frequency that's a multiple of 1Hz. 1Hz is the sample spacing if the sample rate is 44100 and if the number of samples is 44100 - i.e. the analysis epoch is exactly 1 second yielding a frequency sample interval of 1Hz. This would span exactly 27 cycles of the signal. But, you appear to be taking a smaller number of samples of the 27Hz signal because you have to zero pad to get 44100 samples. If you see no leakage then you must be sampling an integer number of cycles of the 27Hz signal as well. Since 44100 is not an integer multiple of 27, but it is an integer multiple of 3*27, then you must be taking some integral number of *3 cycles* of 27Hz such as 24 cycles which would be exactly 39200 samples at 44100Hz. That's the only way to get no leakage. If you take 24 cycles, the actual temporal epoch is 24/27 seconds and the spectral sample interval would be 27/24 or 1.125Hz. That's what the data supports. Then if you zero pad to 1 second the underlying data sequence is interpolated in frequency by a factor of 1.125 - adding one sample for every 8. This happens to generate a sequence in frequency that is exact at multiples of 1.125Hz and 27 is an integer multiple of 1.125 - so all of the energy remains at 27Hz in the transform I do believe. If you take a non-integer number of cycles of the 27Hz signal then you should see leakage in the transform. As has already been pointed out, you have to compute the abs() of the complex values to get the magnitude. dB are relative so you can compute dB relative to anything you like and then shift the result by a constant according to your system parameters to get "absolute" dBv, dB relative to a micropascal, dBm, etc. Fred
Fred Marshall wrote:
[...]
> There is going to be leakage in the process anyway. The only thing that > stops you from seeing leakage, and therefore assuming there is none?, is > that you select a time epoch that covers an integer number of cycles of the > particular signal (27Hz) you're analyzing right now. > > This is a somewhat interesting case in that you are working with a sinusoid > that has a frequency that's a multiple of 1Hz. 1Hz is the sample spacing if > the sample rate is 44100 and if the number of samples is 44100 - i.e. the > analysis epoch is exactly 1 second yielding a frequency sample interval of > 1Hz. This would span exactly 27 cycles of the signal. > > But, you appear to be taking a smaller number of samples of the 27Hz signal > because you have to zero pad to get 44100 samples. If you see no leakage > then you must be sampling an integer number of cycles of the 27Hz signal as > well. Since 44100 is not an integer multiple of 27, but it is an integer > multiple of 3*27, then you must be taking some integral number of *3 cycles* > of 27Hz such as 24 cycles which would be exactly 39200 samples at 44100Hz. > That's the only way to get no leakage.
Fred, it's obvious you've been doing this for quite some time but for newbie, am I led to belive that one would have to throw away ( based on the example provided ) 3 samples in order not to experience leakage? Apologies for being just a dumb software guy. Alot of this is new to me so I'm trying to get my head around your discussion on leakage etc.
> If you take 24 cycles, the actual temporal epoch is 24/27 seconds and the > spectral sample interval would be 27/24 or 1.125Hz. That's what the data > supports. Then if you zero pad to 1 second the underlying data sequence is > interpolated in frequency by a factor of 1.125 - adding one sample for every > 8. This happens to generate a sequence in frequency that is exact at > multiples of 1.125Hz and 27 is an integer multiple of 1.125 - so all of the > energy remains at 27Hz in the transform I do believe. > > If you take a non-integer number of cycles of the 27Hz signal then you > should see leakage in the transform.
Not sure I follow you here.
> As has already been pointed out, you have to compute the abs() of the > complex values to get the magnitude. > > dB are relative so you can compute dB relative to anything you like and then > shift the result by a constant according to your system parameters to get > "absolute" dBv, dB relative to a micropascal, dBm, etc. > > Fred
<forums_mp@hotmail.com> wrote in message 
news:1145083238.243203.157170@t31g2000cwb.googlegroups.com...
> > Fred Marshall wrote: > [...] >> There is going to be leakage in the process anyway. The only thing that >> stops you from seeing leakage, and therefore assuming there is none?, is >> that you select a time epoch that covers an integer number of cycles of >> the >> particular signal (27Hz) you're analyzing right now. >> >> This is a somewhat interesting case in that you are working with a >> sinusoid >> that has a frequency that's a multiple of 1Hz. 1Hz is the sample spacing >> if >> the sample rate is 44100 and if the number of samples is 44100 - i.e. the >> analysis epoch is exactly 1 second yielding a frequency sample interval >> of >> 1Hz. This would span exactly 27 cycles of the signal. >> >> But, you appear to be taking a smaller number of samples of the 27Hz >> signal >> because you have to zero pad to get 44100 samples. If you see no leakage >> then you must be sampling an integer number of cycles of the 27Hz signal >> as >> well. Since 44100 is not an integer multiple of 27, but it is an integer >> multiple of 3*27, then you must be taking some integral number of *3 >> cycles* >> of 27Hz such as 24 cycles which would be exactly 39200 samples at >> 44100Hz. >> That's the only way to get no leakage.
> Fred, it's obvious you've been doing this for quite some time but for > newbie, am I led to belive that one would have to throw away ( based on > the example provided ) 3 samples in order not to experience leakage? > Apologies for being just a dumb software guy. Alot of this is new to > me so I'm trying to get my head around your discussion on leakage etc. > > >> If you take 24 cycles, the actual temporal epoch is 24/27 seconds and the >> spectral sample interval would be 27/24 or 1.125Hz. That's what the data >> supports. Then if you zero pad to 1 second the underlying data sequence >> is >> interpolated in frequency by a factor of 1.125 - adding one sample for >> every >> 8. This happens to generate a sequence in frequency that is exact at >> multiples of 1.125Hz and 27 is an integer multiple of 1.125 - so all of >> the >> energy remains at 27Hz in the transform I do believe. >> >> If you take a non-integer number of cycles of the 27Hz signal then you >> should see leakage in the transform.
> Not sure I follow you here.
>> >> Fred
Well, I did this in a bit of a hurry and may have made a mistake but let's see: Set the sample rate to be 44100Hz. Let's start with exactly 24 cycles of a sinusoid at 27Hz. The temporal epoch is 24/27 seconds and there will be 39200 samples - exactly. If you compute the DFT of this sequence, there will be two non-zero values that correspond to 27Hz and (44100-27)Hz. All the other samples in frequency will be zero or in the numeric noise. I've confirmed this by simulation and it's as it must be. Now, zero pad the sequence out to 1 second or 44100 samples exactly. This doesn't do what I said earlier .... In the first case it's as though the sequence continues indefinitely (is periodic on 24/27 seconds). In the second case it's as though there's a tone burst of length 24/27 seconds that repeats indefinitely with period 1 second. The tone burst has a sinx/x spectrum or, if you like, lots of leakage. I now see no way that this "leakage" can be eliminated. This belies what I thought was the OP's assertion about leakage. In rereading the post it appears Blackman wasn't used because the main lobe is wider than 1 bin and not because there's no leakage as I'd assumed was meant. Fred
"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message 
news:E42dnUxU88hTo9zZnZ2dnUVZ_u-dnZ2d@centurytel.net...
>
> I now see no way that this "leakage" can be eliminated. > This belies what I thought was the OP's assertion about leakage. > In rereading the post it appears Blackman wasn't used because the main > lobe is wider than 1 bin and not because there's no leakage as I'd assumed > was meant. > > Fred
I said the leakage can't be eliminated .... well, if one windows the nonzero part of the sinusoid then the leakage will be reduced at the expense of a wider main lobe / slightly reduced resolution. Fred