[Help-gsl] Requesting feedback on the FFT interface

[Dimitris Papavasiliou]

Hello again,

Upon trying to implement a dct transform with the GSL fft routines I 
noticed a couple of things regarding the fft interface:

1) There are no routines accepting gsl_vector arguments instead of 
(data, stride, n) triplets. Unless I'm missing something this seems 
almost a bug to me.

2) There's a weird assymetry in the interface.  Complex routines have 
*_forward, *_transform, *_inverse and *_backward flavors where 
*_transform accept a direction argument. Real routines on the other hand 
only *_transform and *_inverse flavors and *_transform simply does an 
forward transform.

3) Well I suppose there's a good reason for this but the halfcomplex 
arrays produced by radix2 and mixed-radix routines are different.  This 
tends to make postprocessing of the results of FFTs harder if you need 
to support both radix2 and mixed-radix routines.  For example I thought 
I'd implement DCTs for both routines and it was enough of a headache 
figuring out how these halfcomplex sequences where stored, how to 
multiply the real and imaginary parts with certain sines and cosines 
(in-place) and where to put the results to end up with the real result 
of the DCT without having to do it twice.

I'd therefore like to propose the following changes, either as a patch 
or as an extension, whichever you prefer:

1) Adding gsl_fft_*_*_vector routines.  I have code for that already.

2) Changing both real and complex transforms to a single scheme, either 
*_forward, *_transform, *_inverse, or *_transform, *_inverse.

3) In order to avoid rewriting the final modulation step of the DCT for 
radix2 and mixed transforms I have written a small function that 
permutes the mixed-radix version of the halfcomplex array (the FFTPACK 
version) to the radix2 version (and back).  The radix2 version is 
actually (at least in my case) better since you can process the array a 
pair at a time, possibly overwriting each pair without destroying data 
from other coefficients (well this probably doesn't make sense but I 
can't put it much better).  It is technically O(n^2) but it only 
performs element swaps and it probably could be optimized.
   So I think it might be better to use such a function to permute the 
output of the real mixed-radix functions so that the end-user won't have 
to a) learn and get accustomed to two halfcomplex formats and b) 
potentially write code to do the same thing with both formats.  We can 
use my function or a more optimized one but in any case, I doubt the 
overhead will be a problem.

This one is not a problem with the current interface but:

4) We can add a dct transform pair.  It (hopefully) works for arbitrary 
length sequences.

If the library maintainers and users agree to any of these then I can 
handle 1) and 4) I suppose but I'd need some help from the original 
implementor of the FFT routines for 2) and 3) as I don't have a lot of 
spare time right now.  I'd also need some help testing this stuff.

Sorry for the endless blabbering,
Dimitris P.