fft

I am using an FFT routine for the first time and this @$%& code is not
explained well.

The number of samples must be a power of 2, so I set this number to
256.

There is only one buffer used for input and output, which is an array
of complex.

I feed 256 double-s into the real part of the array.

The values I have to process were sampled at 8000 bytes per second.

I need to calculate the values of frequencies 800Kz and 1000KHz.

As I got it, After processing, the buffer contains all the frequency
spectrum, starting with frequency 0 at element 0 of the array.

Each subsequent element of the array contains the coefficient of an
higher frequency with step=4000/128.

So 
   element         contains the coefficient of frequency 
    128                    4000
     25                     800 
     32                    1000

It is not clear to me what is the use of the complex part of the
buffer. I thought it would contain the value of the frequency itself
(expressed with some black magic formula), but if the frequency raise
is constant at every element of the array, why not to drop this part
at all ?

Also, I found many negative values in the real part of the array: is
it correct ? (I had imagined I would get only positive or zero
values).

Thank you

Peppe Polpo

Peppe Polpo wrote:

> I am using an FFT routine for the first time and this @$%& code is not
> explained well.
> 
> The number of samples must be a power of 2, so I set this number to
> 256.
> 
> There is only one buffer used for input and output, which is an array
> of complex.
> 
> I feed 256 double-s into the real part of the array.
> 
> The values I have to process were sampled at 8000 bytes per second.
> 
> I need to calculate the values of frequencies 800Kz and 1000KHz.
> 
> As I got it, After processing, the buffer contains all the frequency
> spectrum, starting with frequency 0 at element 0 of the array.
> 
> Each subsequent element of the array contains the coefficient of an
> higher frequency with step=4000/128.
> 
> So 
>    element         contains the coefficient of frequency 
>     128                    4000
>      25                     800 
>      32                    1000
> 
> It is not clear to me what is the use of the complex part of the
> buffer. I thought it would contain the value of the frequency itself
> (expressed with some black magic formula), but if the frequency raise
> is constant at every element of the array, why not to drop this part
> at all ?
> 
> Also, I found many negative values in the real part of the array: is
> it correct ? (I had imagined I would get only positive or zero
> values).
> 
> Thank you
> 
> Peppe Polpo

The FFT finds how well the sample correlates to cosines and sines.  So a 
negative real part corresponds to a negative cosine wave, and the 
imaginary part corresponds to a sine wave (just trust me on this, 
without going through all the theory I can't make it clearer).

If you want a number proportional to the total _power_ existent at one 
frequency then you need to square the real and the imaginary part and 
add the result.  This will automatically be a positive number, which is 
what you want.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

On Sun, 20 Jun 2004 18:01:08 -0700, Tim Wescott wrote:

> The FFT finds how well the sample correlates to cosines and sines.  So a
> negative real part corresponds to a negative cosine wave, and the
> imaginary part corresponds to a sine wave (just trust me on this, without
> going through all the theory I can't make it clearer).
> 
> If you want a number proportional to the total _power_ existent at one
> frequency then you need to square the real and the imaginary part and add
> the result.  This will automatically be a positive number, which is what
> you want.

hmmm, don't you have to sqareroot the result again to find the magnitude?
[img, re are imaginary and real component, sqrt is squareroot] 

sqrt(img*img + re*re)

i'm a bloody fft beginner myself though so i might be wrong..

flo

-- 
Palimm Palimm!

florian schmidt wrote:
> On Sun, 20 Jun 2004 18:01:08 -0700, Tim Wescott wrote:
>>If you want a number proportional to the total _power_ existent at one
>>frequency then you need to square the real and the imaginary part and add
>>the result.  This will automatically be a positive number, which is what
>>you want.
> 
> 
> hmmm, don't you have to sqareroot the result again to find the magnitude?
> [img, re are imaginary and real component, sqrt is squareroot] 
> 
> sqrt(img*img + re*re)
> 
> i'm a bloody fft beginner myself though so i might be wrong..
> 
> flo
> 

The square root is required to find the magnitude, but Tim showed how to 
find the power (and he underlined it).  For power, you don't need the sqrt.

-- 
Jim Thomas            Principal Applications Engineer  Bittware, Inc
jthomas@bittware.com  http://www.bittware.com          (703) 779-7770
To mess up a Linux box, you need to work at it; to mess up your Windows
box, you just need to work on it. - Scott Granneman

Peppe Polpo wrote:
> There is only one buffer used for input and output, which is an array
> of complex.
> 

This is known as an "in-place" fft.  It uses a little less memory and 
(maybe) time.  This is mostly important for huge transforms or embedded 
systems.

> I feed 256 double-s into the real part of the array.
> 
> The values I have to process were sampled at 8000 bytes per second.
> 

Do you mean 8000 samples/s?  If not, then you haven't told us how wide 
your data samples are.  Are they 64 bit doubles?  16 bit integers? ....


> I need to calculate the values of frequencies 800Kz and 1000KHz.
> 

This violates the Nyquist criteria.  If you sample a real signal at 
8kHz, only frequencies less than 4kHz can be reliably detected.


Do you mean 800Hz and 1000Hz???

If that is the case then you won't be able detect 800Hz exactly ( i.e. 
without spectral leakage).  To do that for 800 and 1000 with a 8K 
sampling rate, you need an FFT size with a factor of 40 in it.

That would force you to use a mixed radix FFT library.


If you decide to go that route,  two good choices for mixed-radix fft 
libraries are


FFTW, which is large, complicated, and fast as hell.
http://www.fftw.org/

KISSFFT, which is tiny, simple as possible, half as fast.
(and authored by me)
http://sourceforge.net/projects/kissfft/

-- Mark Borgerding

"Jim Thomas" <jthomas@bittware.com> wrote in message
news:10ddm0pn6dk2624@corp.supernews.com...
> florian schmidt wrote:
> > On Sun, 20 Jun 2004 18:01:08 -0700, Tim Wescott wrote:
> >>If you want a number proportional to the total _power_ existent at one
> >>frequency then you need to square the real and the imaginary part and
add
> >>the result.  This will automatically be a positive number, which is what
> >>you want.
> >
> >
> > hmmm, don't you have to sqareroot the result again to find the
magnitude?
> > [img, re are imaginary and real component, sqrt is squareroot]
> >
> > sqrt(img*img + re*re)
> >
> > i'm a bloody fft beginner myself though so i might be wrong..
> >
> > flo
> >
>
> The square root is required to find the magnitude, but Tim showed how to
> find the power (and he underlined it).  For power, you don't need the
sqrt.
>

To take this even further, you very rarely would WANT to do the square root.
Typically you'd be comparing the magnitude to some kind of threshold or
something.  In that case it's much more efficient to just square your
threshold.  A square root can be a fairly "expensive" operation.

Brad

Brad Griffis wrote:
> A square root can be a fairly "expensive" operation.

If you're interested in square roots, Paul Hsieh has a cool site 
dedicated to various implementations and techniques.

http://www.azillionmonkeys.com/qed/sqroot.html

-- Mark Borgerding

Mark Borgerding <mark@borgerding.net> wrote in message news:<BONBc.140684$DG4.96988@fe2.columbus.rr.com>...

sorry for my lack of precision.

I must have been very tired when I wrote the message ... 

> Do you mean 8000 samples/s?  

yes

> your data samples are.  Are they 64 bit doubles?  16 bit integers? ....

signed 16 bit

> Do you mean 800Hz and 1000Hz???

er...yes

>To do that for 800 and 1000 with a 8K sampling rate, you need an FFT
size with a factor of 40 in it.

could you please expand on this ?

Will I have to use a buffer size multiple of 40, like 400 ?

> That would force you to use a mixed radix FFT library.

I read about the algorithm in Wikipedia

> http://sourceforge.net/projects/kissfft/

I need to work in Delphi. Can your code be compiled into a DLL ?

Thank you

Peppe Polpo

Peppe Polpo wrote:
> Mark Borgerding <mark@borgerding.net> wrote in message news:<BONBc.140684$DG4.96988@fe2.columbus.rr.com>...
>>To do that for 800 and 1000 with a 8K sampling rate, you need an FFT
 >> size with a factor of 40 in it.

I think I said in order to *exactly* detect those frequencies with a 256 
point FFT, you need a factor of 40.  Some applications can live with a 
little frequency aliasing.  You must decide for yourself.

> 
> could you please expand on this ?


Why 256 won't *exactly* work:

If you do a 256 point fft on a 8KS/s sequence, then the frequency 
resolution will be 8000/256 == 31.25 Hz.  Any sinusoids that are 
multiples of 31.25 can be exactly captured by the fft. Their periods 
coincide with the sampling window length.

1000/31.25 == 32, so there's no problem with 1000 Hz.
800/31.25 == 25.6, so an 800 Hz tone will have its energy spread across 
multiple bins.  With the bins at 26*31.25Hz and 25*31.25Hz receiving the 
lion's share.

For more info, look for "spectral leakage", "frequency leakage", and 
"frequency resolution" in a DSP book, or online.


Why 40?

8000/1000 == 8
8000/800 == 10

So if you want both 800 and 1000 to be multiples of the frequency 
resolution, then the fft length should be a multiple of both 8 and 10.
The LCM of 8 & 10 is 40.


>>http://sourceforge.net/projects/kissfft/
> 
> I need to work in Delphi. Can your code be compiled into a DLL ?

I don't see why not.  If you do it, send me your DLL project file.  I 
can put it in the package.

> Thank you
> 
> Peppe Polpo

Prego,
Mark

> The FFT finds how well the sample correlates to cosines and sines.  So a 
> negative real part corresponds to a negative cosine wave, and the 
> imaginary part corresponds to a sine wave (just trust me on this, 
> without going through all the theory I can't make it clearer).

A simple way of explaining this is that the use of the sine and cosine
allows you to represent the phase of the signal. Equal amounts of both
result in a 45 degree phase shift.