Re: Gsl FFT <-> FFTW compatibility

[Stefan Westerfeld <stefan space twc de>]

   Hi!

On Mon, May 16, 2011 at 11:35:32PM +0200, timj gnu org wrote:
> thanks for working on this.
> Attached you will find 4 patches which you can apply with git-am.
> These are your fft branch patches rebased onto current master plus
> documentation fixes where recent changes introduced conflicts.
> 
> Unfortunately the code can't be merged as is (thus this email)
> because of the following issues:
> 
> * testfft.c - buffer overrun in line 200: memcpy (dft_in, ref_fft_in, MAX_FFT_SIZE * sizeof (dft_in[0]));
>   from glibc:
> 	Checking reference fft for size 2: [*** buffer overflow detected ***: /opt/src/beast/bse/tests/.libs/lt-testfft terminated
>   This is possibly because MAX_DFT_SIZE != MAX_FFT_SIZE.

Right, I changed it to MAX_DFT_SIZE now.

>   Why's there a MAX_DFT_SIZE in the first place?

Because the DFT implementation is a lot slower than your
reference_power2_fftc implementation. So we only check small FFT sizes with
the DFT implementation (which is easy to check for correctness, because the
code is very simple), and all FFT sizes with your implementation.

> * Why does gsl_power2_fftsc() contain this line?
>   /* { unsigned int i; for (i = 0; i < n_values << 1; i++) rivalues_out[i] *= (double) n_values; } */

It was there before my changes, so I kept it as it is.

> * Why does gsl_power2_fftar(), an AnalysisReal function now contain
>   a call to a *synthesis* function? gsl_power2_fftsc
> 
> * Why does gsl_power2_fftsr(), a SynthesisReal function now contain a
>   call to an *analysis* function? gsl_power2_fft*analysis_skip2
> 
> For the latter two, should there possibly be some a<->s or analysis<->synthesis
> renames in place, to fix the obvious consistency discrepancy?

To make gslfft behave like FFTW, we change gsl_power2_fftac to do what was
previously done by gsl_power2_fftsc. Same thing for gsl_power2_fftac. We also
rename gsl_power2_fft*analysis_skip2 to gsl_power2_fft*synthesis_skip2. So for
complex ffts, everything is consistent.

However, there are the real variants, which used to work like that

fftar: real data -> fftac -> post processing -> complex output
fftsr: complex data -> pre processing -> fftsc -> real output

I changed these to FFTW behaviour by adjusting the (pre/post) processing steps,
so we have

fftar: real data -> old fftac -> post processing* -> complex output
fftsr: complex data -> pre processing* -> old fftsc -> real output

which is

fftar: real data -> fftsc -> post processing* -> complex output
fftsr: complex data -> pre processing* -> fftac -> real output

It would be possible to rederive the post/pre processing from numerical
recipies to do the right thing without still using the old style functions.
However that would be quite a bit of work, so I preferred to do it my way
which relies on the "old fftac" step and the "old fftsc" step.

I've put your version + dft size changes in a new branch, which is:

git fetch http://space.twc.de/public/git/stwbeast.git gslfft-fftw-compat2:gslfft-fftw-compat2

   Cu... Stefan
-- 
Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan