Hi Liam,

As I briefly indicated on #lisp and contrary to our earlier belief, there _is_ at least one case in which a user would like to supply the stride, vector size n and offset to the FFT functions. Namely, this is for the 2-dimensional (more generally of course, N-dimensional) Fourier transform.

GSL only provides one-dimensional FFTs. Naturally, an N-dimensional FFT can be easily implemented from a one-dimensional function. From saying this, you can probably already immediately see what I mean, but nevertheless, to be sure there's no mistake here, I'll give a more concrete example.

Let us approach the problem in C, and consider a 2D array of size a x b, where a is the number of rows and b is the number of columns. This is a linear, contiguous array in C (a "raw" C array). A 2D FFT can be performed by first FFT'ing one dimension and then the other.

Thus, we first FFT each row, and then each column. The order actually doesn't matter. In C, this would be handled by first performing a times the (forward) FFT, thus applying an FFT to each row, where the vector size n = b, the stride = 1 and each consecutive FFT call supplies the vector with an offset with a multiple of b.

Next up is the FFT of each column. In C we perform b times the (forward) FFT, thus applying an FFT to each column where the vector size n = a, the stride = b, and each consecutive FFT call supplies the vector with an offset equal to the column number.

For N-dimensional arrays, simply lather, rinse and repeat.

This shows the necessity of the optional (keyword?) argument that indicates the size of the vector to be FFT'd, the stride (which is already in there) and also the offset (which is missing at the moment). This also shows that the Fourier transform environment that is in fft-interface.lisp is very nice to have, when dealing with non-radix-2 array lengths. Each consecutive call when FFT-ing in one direction only needs to allocate the workspace and wavetable once and, more importantly, calculate the wavetable once, while the FFT is repeated on each row (or column).

I have to admit that I haven't looked at how multi-dimensional marrays or more correctly, marrays representing tensors of rank > 1, from the GSLL point of view. As I understand, the interface allows for a truly multi-rank approach (using maref with multiple indices), as with the normal CL arrays and aref.

In my opinion, a nice approach would be to make the whole FFT by default work on the entire marray, thus performing multi-dimensional FFTs when marrays have higher ranks. That way, we can hide the messing about with strides, sizes and offsets from the user. Just think, loading an image into a marray, and just calling forward-fourier-transform, and instantly (depending on the speed of your computer...) have a Fourier-space image!

On the other hand, while such an interface would be nice for the casual user (and probably for most normal users) I don't know if there are "power users" that (would want to use GSL's FFT functions, and) would require direct control over such parameters.

As I may have mentioned before on #lisp, I have used GSL's FFT some time ago, as well as FFTW briefly. FFTW has functions for higher-dimensional FFTs, and deducing from that interface, I would say that an FFT function that "does the right thing (TM)" depending on the rank of the marray would not be so bad.

Anyway, sorry for the long post. I hope my explanations/considerations/ideas are useful and I would not mind implementing some/most of this in an experimental interface, possibly with some help with respect to some of the fancy internals of GSLL.

Please let me know what you think (and any suggestions from users of GSLL and FTT, I'm looking at you, Mirko...) are greatly appreciated!

Thanks, Sumant