[Matrix-SIG] Re: [PYTHON MATRIX-SIG] FFT
pcj+NOSPAM@primenet.com
pcj+NOSPAM@primenet.com
18 Mar 1998 21:43:20 -0800
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