[Matrix-SIG] Re: [PYTHON MATRIX-SIG] FFT

[pcj+NOSPAM@primenet.com]

A really, really long time ago, Jim Hugunin wrote:
> real_fft returns the left half of the fft of a real array (the fft of a 
> real array is conjugate symmetric, so only half of it is needed).

Except it doesn't.  fftpackmodule.c needs this patch:

*** fftpackmodule.c~	Fri Mar 21 13:10:00 1997
--- fftpackmodule.c	Tue Mar 17 22:25:00 1998
***************
*** 132,138 ****
    dptr = (double *)data->data;
    
    for (i=0; i<nrepeats; i++) {
! 	memcpy((char *)rptr, dptr, npts*sizeof(double));
      rfftf(npts, rptr+1, wsave);
  	rptr[0] = rptr[1];
  	rptr[1] = 0.0;
--- 132,138 ----
    dptr = (double *)data->data;
    
    for (i=0; i<nrepeats; i++) {
! 	memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
      rfftf(npts, rptr+1, wsave);
  	rptr[0] = rptr[1];
  	rptr[1] = 0.0;

As David Buscher correctly pointed out, real_fft currently ignores the
first point of the sequence.

>  inverse_real_fft is not the inverse function of real_fft, but it instead 
> returns the left half of the inverse fft of a real valued vector.  This 
> naming convention is from fftpack, and I don't think it's bad enough to 
> fix.

This is not my understanding of fftpack's "real backwards" routine.
(Admittedly the documentation is a little misleading.)  rfftb takes an
array of the left half of the real and imaginary components of the
fourier transform of a real sequence (in the same format as the output
of rfftf), and returns the real sequence.  May I suggest instead the
following semantics for inverse_real_fft:

  - accepts a complex sequence, assumed to be the left half of a
  conjugate symmetric sequence

  - returns a real sequence of length 2n-2 if the imaginary part of
  the last value is equal to 0 and length 2n-1 otherwise.

This implementation would have the (IMHO desirable) property that
inverse_real_fft(real_fft(seq)) would in general equal seq.  My only
concern is that if the imaginary part of the last point was not
precisely zero (due to roundoff error), you could get a sequence one
longer than the input.  Offhand, I can't think of a robust algorithm
to decide if imag(seq[n]) was sufficently close to 0 to shorten the
returned sequence by 1.

--
Paul C. Janzen                    Research is what I'm doing when I
http://www.acs.psu.edu/users/pcj  don't know what I'm doing.
                                     - Wernher von Braun