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