I have not been following this fft discussion because I thought it was
all about twiddling constant multipliers or similar stuff. If it
turns out that the program is really using n^2 or worse operations to
copy stuff into an array, that is really an unnecessary expense.
I have been using (power of 2) fft for experiments in multiplying
polynomials, but that can work on bigfloats, all written in Lisp. I
have not timed this in comparison with FFT-Bordeaux, which I am not
familiar with.
The code is in the program "four1" in this file
http://www.cs.berkeley.edu/~fateman/generic/mpfr.lisp
and in slightly other guises in other files in that directory.
I suspect the program was copied from Numerical Recipes, but maybe not.
It is possible that one could have a much faster FFT for double-floats
by exploiting better cache strategies or parallel computing, in case
this makes a big difference for some application.
I also found some other FFT implementations in lisp. Clearly a poor
lisp compiler will produce a slow program; some of the lisp compilers
around are pretty good though.
RJF
On Aug 23, 2009, at 9:00 AM, ?iga Lenar?i? wrote:
>
> On Aug 23, 2009, at 5:22 AM, Robert Dodier wrote:
>
>> On 8/22/09, ?iga Lenar?i? <ziga.lenarcic at gmail.com> wrote:
>>
>>> I've tried the 'new' fft in Maxima with support for lists. While all
>>> is nice and it 'works', after inspecting sourcecode it seem to be
>>> implemented in a very inefficient way. I haven't tested the speed of
>>> the internal fft algorithm, but since it's a lisp rewrite from a
>>> fortran book, probably there are faster algorithms.
>>
>> What faster algorithm would that be?
>>
>>> (dotimes (i ($length a1))
>>> (setf (nth (1+ i) a1) (m+ (nth (1+ i) a1) (m* '$%i (nth (1
>>> + i) a2)))))
>>>
>>> can't be efficient on lists of length perhaps 1024*8. Surely the
>>> list
>>> should be traversed in some other way (and constructing a complex
>>> number with m+ and m* is probably an overkill also).
>>
>> Surely indeed. If you would like to suggest an alternative
>> (i.e. post a patch) I am willing to consider it.
>>
>>> Maxima's 'fft' is terribly slow, but doesn't need to be.
>>
>> Well, a claim like this would be a lot stronger if accompanied
>> by some evidence or argument from first principals or something.
>>
>>> I would suggest also to take a look at existing implementations in
>>> CL
>>> like Bordeaux-FFT which also allows for arrays which are not of
>>> length of power of 2 (but uses a much slower algorithm in such
>>> cases).
>>
>> Yes, an function to handle cases other than length = 2^n would
>> be a welcome addition. If you'd like to spell out how to do it
>> (i.e. post patches & details about how this whole business with
>> Bordeaux-FFT is going to work) then I'm willing to consider it.
>>
>> FWIW
>>
>> Robert Dodier
>
> I've made a very quick test:
>
> (%i51) a:makelist(random(1.0),i,1,1024*64)$
> Evaluation took 1.1000 seconds (1.1070 elapsed) using 35.064 MB.
> (%i52) fft(a)$
> Evaluation took 30.6790 seconds (30.7180 elapsed) using 97.500 MB.
> (%i56) myfft(a)$
> Evaluation took 0.3490 seconds (0.3500 elapsed) using 21.499 MB.
>
> myfft uses bordeaux-fft. I wrote a thin layer to convert maxima list
> to a complex lisp array and back. The difference grows for bigger
> sample sizer (for which the reason must be inefficient list
> conversion of O(N^2) in current fft). Probably the internal fft
> routine of bordeaux-fft is faster too (for instance it reuses some
> data from previous calculation if the sample is of the same size -
> which is often the case when working with fft).
>
> Ofcourse I don't have a 'drop-in' replacement ready yet (no error
> checking, no support for maxima arrays yet and bordeaux-fft uses
> different scaling)...
>
> I think I could write a proper maxima interface, along with perhaps
> sine and cosine transform. I'm not very familiar with maxima arrays
> though.
>
> Any interest?
>
> For a smaller sample the difference is very significant too
> (%i57) a:makelist(random(1.0),i,1,1024*8)$
> Evaluation took 0.1330 seconds (0.1330 elapsed) using 4.384 MB.
> (%i58) fft(a)$
> Evaluation took 0.7940 seconds (0.7980 elapsed) using 10.315 MB.
> (%i59) myfft(a)$
> Evaluation took 0.0480 seconds (0.0490 elapsed) using 3.658 MB.
>
> This is on SBCL.
>
> I think bordeaux-fft only supports 2^n samples. Fast non 2^n samples
> code is probably very complicated (it's better to just interface to
> FFTW library like everybody does), but a 'naive' DFT implementation
> could be put together just to be able to handle a general case. It
> has O(N^2) which might be fine for small sample sizes and one-time
> computations...
>
> Regards,
> Ziga
>
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima