FFT usage

Hi all,
I'm using an FFT routine as given at
http://www.yov408.com/html/codespot.php?gg=36
- in fact, its been used by my group a lot, so it has a good reputation
;)

I have a few questions to ask about it though.
I've modified the code to allow the user to specifiy whether they want
the complex conjugate of an FFT.
The logic seems fairly simply - just a simple flag and a reversal of
the imaginary part output value's sign.

I've extracted the part that checks whether the size of the input
arrays are a power of two and if so, works out a new array size. The
reason for this is as follows: I have two arrays, the first being quite
big, the second small.
By calculating the minimum power of two size of the input arrays, this
will resize both arrays, padding out both of them with 0's until they
are the requisite size.

However, I read (and I can't find it now!) that there are problems with
the FFT when the NumSamples = NumBits.
NumSamples is the size of the new power of two arrays large enough to
hold both arrays, and I have set NumSamples = NumBits in the code. Can
I do this? My output arrays are both zero, which is obviously wrong!

To give an example, I have Input array1, size = 240 x320; the
NumSamples array (which is an array of ints) is then recalculated to be
256x512. The second input array2, size=30x30, is also set to 256x512. I
then do various things like multiply the elements of FFT(array1) by the
FFT*(array2) , but again the results are zero.

Can anyone please help?

Thanks a lot

Paul

Use an industry standard FFT implementation:  http://www.fftw.org/


paul@paullee.com wrote:
> Hi all,
> I'm using an FFT routine as given at
> http://www.yov408.com/html/codespot.php?gg=36
> - in fact, its been used by my group a lot, so it has a good reputation
> ;)
>
> I have a few questions to ask about it though.
> I've modified the code to allow the user to specifiy whether they want
> the complex conjugate of an FFT.
> The logic seems fairly simply - just a simple flag and a reversal of
> the imaginary part output value's sign.
>
> I've extracted the part that checks whether the size of the input
> arrays are a power of two and if so, works out a new array size. The
> reason for this is as follows: I have two arrays, the first being quite
> big, the second small.
> By calculating the minimum power of two size of the input arrays, this
> will resize both arrays, padding out both of them with 0's until they
> are the requisite size.
>
> However, I read (and I can't find it now!) that there are problems with
> the FFT when the NumSamples = NumBits.
> NumSamples is the size of the new power of two arrays large enough to
> hold both arrays, and I have set NumSamples = NumBits in the code. Can
> I do this? My output arrays are both zero, which is obviously wrong!
>
> To give an example, I have Input array1, size = 240 x320; the
> NumSamples array (which is an array of ints) is then recalculated to be
> 256x512. The second input array2, size=30x30, is also set to 256x512. I
> then do various things like multiply the elements of FFT(array1) by the
> FFT*(array2) , but again the results are zero.
> 
> Can anyone please help?
> 
> Thanks a lot
> 
> Paul

Well, a few points:

Firstly, although the implementation looks good, my superiors don't
want anything to do with the GPL.

Also (and this is a niggle but is solveable), my input data array would
be a linked list of floating point values. I can convert to a simply
array, but this would add overheads in other parts of the program (and
I can't use a fixed sized array as I don't know how much data I will be
loading in at run time)

Cheers

Paul
--
http://www.paullee.com

paul@paullee.com wrote:
> Well, a few points:
>
> Firstly, although the implementation looks good, my superiors don't
> want anything to do with the GPL.

I *think* there is a commercial lisence available, too. It's kind
of a marketing trick; hook the students with the GPL and get
them addicted. Then they will  be compelled to buy the commercial
lisence when they are in positions with power and $$$. That's at
least part of the explanation of how matlab gained the position they
have today. (The difference being that Matlab sold/sells very cheap
student lisences, not GPL freeware)

> Also (and this is a niggle but is solveable), my input data array would
> be a linked list of floating point values. I can convert to a simply
> array, but this would add overheads in other parts of the program (and
> I can't use a fixed sized array as I don't know how much data I will be
> loading in at run time)

These are two very bad solutions that hit wide on each side
of the sensible one: a dynamically allocated array of floats.
A linked list of floats is a ridiculous idea if you want any sort
of performance. A statically allocated array is fast, but requires
you to be clairvoyant as to the size of the larges data set your
program will ever be used on.

The sensible way that preserves performance and also only
requires mere human skills from the programmer, is to
dynamically allocate just enough *continuous* memory at
run-time. 

Rune

Rune Allnor wrote:

> paul@paullee.com wrote:

(snip)

>>Also (and this is a niggle but is solveable), my input data array would
>>be a linked list of floating point values. 
(snip)

> These are two very bad solutions that hit wide on each side
> of the sensible one: a dynamically allocated array of floats.
> A linked list of floats is a ridiculous idea if you want any sort
> of performance. 

I haven't thought of it in too much detail, but it might be that a 
linked list wouldn't be so bad for an FFT.

If you don't need to keep the input list, it might be that it can
be done by rearranging the links instead of exchanging array elements
in the butterfly pattern.

If you do need to keep the input, find the length, allocate the 
appropriate array, FFT, generate the output linked list.

-- glen

Shlomo Kashani wrote:
> Use an industry standard FFT implementation:  http://www.fftw.org/

i think, under some conditions (like we don't need the absolute best
optimization for speed), there are surviving reasons for using a simple
and fully unencumbered FFT routine.  problem for the OP, besides not
having an idea on how to help him (other than to suggest to stay away
from non-power of two FFTs), N=2^p is likely another of those
conditions that one might choose something other than an industry
standard.

r b-j

I'll look into the FFTW license, but I think I can guess what the
response will be from those "up high"!

I'm now looking at ways to simplify the data set which will hopefully
lead to a faster FFT.

The input data are real values, between 0 and 256. There are no complex
input values.

Does this mean that I can limit the size of the input/output arrays
(use integers, rather than longs or high precision data types)? Also,
will all the imaginary data arrays be zero too?

I've had an experiment with the code given in the link in my original
posting, but it produces some odd results. I thought that if you took a
data set, did an FFT on it, and then did an inverse FFT on this output,
you'd get the original data set back? I must be doing something
wrong....

paul@paullee.com wrote:
> I'll look into the FFTW license, but I think I can guess what the
> response will be from those "up high"!

If you can't use FFTW under the GPL (note that the GPL has essentially
no effect for "in-house" software since it only affects distribution)
and don't want to pay MIT for a less restrictive license, there are
plenty of other decent FFT codes out there with unrestrictive licenses.
 e.g. for power-of-two sizes the code by Ooura is very good, and a
compact mixed-radix code with reasonable performance is kissfft.

> The input data are real values, between 0 and 256. There are no complex
> input values.
>
> Does this mean that I can limit the size of the input/output arrays
> (use integers, rather than longs or high precision data types)? Also,
> will all the imaginary data arrays be zero too?

The output will still be complex, but half of the output will be the
complex-conjugate of the other half.  Because of this you can gain
roughly a factor of two in time and storage (you don't need imaginary
parts for the input, and only have to store half of the output) by
using a specialized FFT for real data.  FFTW, Ooura's code, and kissfft
all include such specialized code.

You're not going to be able to get away with integer (fixed-point) data
types unless you are willing to sacrifice substantial accuracy...the
DFT of integer input does not yield integer output.  If you have an
embedded CPU without floating-point hardware this might be worth it,
but on a general-purpose CPU you should just use floating-point, which
will be just as fast and far more accurate.

Don't use linked lists if you care about FFT performance.

Cordially,
Steven G. Johnson