On Aug 24, 2009, at 4:10 PM, Raymond Toy wrote:
> FWIW, I took at look at maxima's fft and bordeaux fft.
>
> I ran a few tests, after fixing two things in maxima's fft:
> o log-base2 needs work. (log-base2 8192) returns 12 for me, when it
> should return 13.
> o I changed the array declaration in fft-dif-internal to be
> (simple-array double-float (*)) to make a fair comparison.
>
> With these changes, here are some test results for computing fft's of
> various lengths. (All run on a 1.6 GHz sparc with cmucl.)
>
> Length Maxima Bordeaux
> 16384 0.02 0.01
> 32768 0.05 0.02
> 65536 0.10 0.09
>
> I suspect the reason for the difference in speed is that bordeaux fft
> caches the sin/cos twiddle factors. Maxima's fft doesn't do that, but
> it could be made to do that easily. Or we could just take the code,
> since it is GPL. But it uses complex arithmetic, so those Lisp's that
> don't have unboxed complex arithmetic and arrays will surely suffer.
>
> Ray
>
Yes, the difference between fft algorithms is not that big.
I have a working fft.lisp now (which uses Ray's algorithm). I don't
know how to prepare a patch and probably this is not the final
version of it. It supports lists and lisp arrays (the ones found in
Maxima - by calling make_array('any , 4); for instance). I know there
is another kind of Maxima array... $fft always outputs a list for
now. Extending this to previous functionality is a matter of writing
a few short functions if we want $fft to output the same kind of data
that it recieves.. It's a pain to work with Maxima arrays in the
code, though, and in Maxima also. Luckily one can use quite large
lists in Maxima with reasonable speed on todays computers...
I've attached fft.lisp - I hope the attachment gets through.
I get
(%i110) b: table( random(1.0), [i,1024*64])$
Evaluation took 0.3000 seconds (0.3010 elapsed) using 18.562 MB.
(%i111) fft(b)$
Evaluation took 0.7180 seconds (0.7210 elapsed) using 164.494 MB.
(%i112) myfft(b)$
Evaluation took 0.3860 seconds (0.3870 elapsed) using 29.254 MB.
Maxima's $fft is now only 2x slower than bordeaux-fft in SBCL. I'm
still all for using bordeaux-fft since it's free GPLd code (and has
been even maintained, whatever that means), but if complex lisp
numbers are slow on other lisps, then it's not a good idea. I only
have Clozure CL and SBCL on my computer. Timings from ECL and GCL
would matter - perhaps I'll even install ECL and compare...
I might have switched the 'forward'/'backward' fft routines by
accident..
I have a question: right now I'm constructing Maxima complex numbers
with pushing such forms into a list:
(push `((mplus) ,(aref realparts i) ((mtimes) $%i ,(aref imagparts i)))
ans)
realparts and imagparts are double-float arrays. The result of this
seems correct. Should it be done another way (via m+ ?). Should $fft
not know about Maxima internal representation? It seem to me that
most Maxima functions construct expressions directly, but this might
be bad practice.
fft.lisp can probably be further optimized with (declare (optimize
(speed 3) (safety 0))) in inside loops and perhaps some more type
declarations...
Regards,
Ziga
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fft.lisp
Type: application/octet-stream
Size: 8474 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20090824/0ad67e5a/attachment.obj
-------------- next part --------------