implement fft and related functions

[ViralBShah]

We need fft before launch. I assume FFTW is the only realistic library we can use. If it is too much of a pain to integrate until our full C FFI is ready, we can go with FFTPACK.

[JeffBezanson]

How does this look:

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

Or is FFTPACK the kind of code that's really hard to beat?

[StefanKarpinski]

My understanding is that FFTW is the kind of code that's hard to beat, but I suspect Viral and others know much more about it than I do.

[ViralBShah]

FFTPACK is an old well known code that many projects use. I don't think it is particularly fast or anything.

It may also be possible to write a useful fft in julia, as a placeholder until we get fftw.

If kissfft is not used by at least a few projects, we will end up becoming beta testers.

-viral

On Jun 19, 2011, at 2:01 AM, JeffBezanson wrote:

How does this look:

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

Or is FFTPACK the kind of code that's really hard to beat?

Reply to this email directly or view it on GitHub:
#60 (comment)

[StefanKarpinski]

I just pushed bc774de which lets you download, compile and setup FFTW. Haven't tried actually using it yet, but it's a start.

[ViralBShah]

From the manual, this is what the code looks like. So, either, we need to write a simple thin wrapper that is callable from julia, or wait until julia can support a more general C FFI.

Another benefit of FFTW could be that we can use their matrix transpose routines too, which are quite fast.

#include <fftw3.h>
...
{
    fftw_complex *in, *out;
    fftw_plan p;
    ...
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    ...
    fftw_execute(p); /* repeat as needed */
    ...
    fftw_destroy_plan(p);
    fftw_free(in); fftw_free(out);
}

-viral

[ViralBShah]

Coming to think of it, we can just use the Fortran calling interface, and get fftw to work immediately.

http://fftw.org/fftw3_doc/Calling-FFTW-from-Fortran.html#Calling-FFTW-from-Fortran

-viral

On Jun 19, 2011, at 4:31 AM, StefanKarpinski wrote:

I just pushed bc774de which lets you download, compile and setup FFTW. Haven't tried actually using it yet, but it's a start.

Reply to this email directly or view it on GitHub:
#60 (comment)

[StefanKarpinski]

I don't get why we can't use the C interface. Because of the data layout or because of the plan bit? The data layout is actually not a problem since fftw_complex is just a pair of doubles, so a Julia array of Complex128's has the desired layout. Not sure about the plan bit admittedly. Why is the Fortran interface better?

[ViralBShah]

How do we create all the structs that fftw defines in fftw.h in julia and pass them through the C interface? If we can, then we can use the C interface - but I always thought this was not possible.

-viral

On Jun 19, 2011, at 5:14 AM, StefanKarpinski wrote:

I don't get why we can't use the C interface. Because of the data layout or because of the plan bit? The data layout is actually not a problem since fftw_complex is just a pair of doubles, so a Julia array of Complex128's has the desired layout. Not sure about the plan bit admittedly. Why is the Fortran interface better?

Reply to this email directly or view it on GitHub:
#60 (comment)

[ViralBShah]

The only reason to use the fortran interface is that it is likely to be free of passing the C structs around.

-viral

On Jun 19, 2011, at 5:14 AM, StefanKarpinski wrote:

I don't get why we can't use the C interface. Because of the data layout or because of the plan bit? The data layout is actually not a problem since fftw_complex is just a pair of doubles, so a Julia array of Complex128's has the desired layout. Not sure about the plan bit admittedly. Why is the Fortran interface better?

Reply to this email directly or view it on GitHub:
#60 (comment)

[ViralBShah]

Well, seems like I can just treat all fftw data structures as void * and pass them around. Keeping this bug open until we implement most of the fft related functions - FFT2, FFTN, FFTSHIFT, IFFT, IFFT2, IFFTN.. We probably also need to do some saving of the plans and such.

julia> x = complex(ones(8),ones(8))
[1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im,1.0 + 1.0im]

julia> fft(x)
[8.0 + 8.0im,0.0 + 0.0im,0.0 + 0.0im,0.0 + 0.0im,0.0 + 0.0im,0.0 + 0.0im,0.0 + 0.0im,0.0 + 0.0im]

[JeffBezanson]

Excellent!
Don't forget fftw_destroy_plan.

[ViralBShah]

Good catch - fixed.

Commit 021afbf also makes includes fft in the regular julia build and startup.

[StefanKarpinski]

This is awesome. Just what I hoped might happen if I setup the download and compilation bit :-)

[ViralBShah]

Yes - I wouldn't have gone through all this just to try out and see if it would work. However, with your download and compile done, I found myself sitting on a plane, and thought, what the heck, let's try it out.

-viral

On Jun 22, 2011, at 12:37 PM, StefanKarpinski wrote:

This is awesome. Just what I hoped might happen if I setup the download and compilation bit :-)

Reply to this email directly or view it on GitHub:
#60 (comment)

[StefanKarpinski]

Excellent :-D

[ViralBShah]

Discussion on subarray relevant here, for ability to compute fft on slices of arrays.

http://groups.google.com/group/julia-math/browse_thread/thread/f07b1f6047d98a8c/

[ViralBShah]

Good discussion on real valued FFT, and general stuff - including posts by Steve Johnson:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/39968

[StefanKarpinski]

Yeah, that is an awesome thread. Gives a lot of faith in the quality of the FFTW implementation. Glad we're using it :-)