[Matrix-SIG] FFT:real_fft/FFT:inverse_real_fft
Warren B. Focke
warren@pfiesteria.gsfc.nasa.gov
Wed, 14 Jul 1999 11:32:59 -0400 (EDT)
Some confusion has been generated by FFT:inverse_real_fft.
Normally I like confusion is a good thing, but someone asked me to look at
this, so I'll raise^H^H^H^H^Hmake an exception.
FFT:real_fft was coded with the intent that it return the
nonnegative frequency part of the Fourier transform of data which is real
in the time domain. Only the nonnegative frequencies are necessary
because this data has hermitian symmetry in the frequency domain.
Many users expect FFT:inverse_real_fft to perform the inverse of
that operation - that is, to transform the nonnegative frequency part of
data which has hermitian symmetry in the frequency domain to the time
domain, where it would be real. It was, however, coded with the intent
that, given data which is real in the frequency domain, it would return
the nonnegative time portion of the time-domain representation of that
data, which would have hermitian symmetry about t=0.
Unfortunately, neither FFT:inverse_real_fft nor FFT:real_fft
perform as advertised, and I know that RFFTB (the routine in FFTPACK which
is used by FFT:inverse_real_fft) has the intent of the first transform
described in the previous paragraph. Or, at least, I have written working
programs based on that assumption. So I took a look. The results:
real_fft: Fixed. Does what it was intended to.
real_inverse_fft: Does the inverse of real_fft.
hermite_fft: Takes the nonnegative time portion of data with
hermitian symmetry, returns the (real) frequency
domain representation of that data.
hermite_inverse_fft: Does the inverse of hermite_fft - this is what
inverse_real_fft was meant to do.
inverse_real_fft: Another name for hermite_inverse_fft.
Test 'em, if you would. The test code is designed to work with the
original versions, too.
All have the same calling semantics as the routines that existed
before: real_fft(array, npts, axis). npts always refers to the number of
points in the whole sequence, not the number in the nonnegative
frequency/time portion of data with hermitian symmetry (that is, it's N,
not N/2+1). If this number is not supplied, there are two possible values
for the real_inverse and hermite transforms: (array.shape[axis]-1)*2 and
(array.shape[axis]-1)*2+1. The routines always guess the former, so if
you want your result from these routines to have an odd number of points,
you must supply it.
real_inverse_fft and hermite_fft pad the data to twice the size it
needs to be (this happens in _raw_fft). I didn't see a way around this
without special-casing _raw_fft or losing the ability to get back an odd
number of points.
I chose names that did not interfere with the old ones. I think I
would prefer to name real_inverse_fft and hermite_inverse_fft
inverse_real_fft and inverse_hermite_fft instead. This should not break
old code since inverse_real_fft did not work before. Opinons?
I'll be happy to document these after they are tested and naming
etc. issues are resolved. (OK, maybe not happy, but I'll do it.)
Looking over things while packing this up just now, I noticed that
real_fft2d does not do the right thing. I'll fix that later, I wanna get
this out now.
Warren Focke
Eh? Wassat? Oh, the code! Here, apply with "patch -Np0" from the directory
where you unzipped LLNLDistribution11.tgz.
----------------------- cut here --------------------------------------
diff -Naur orig.LLNLDistribution11/Numerical/Lib/FFT.py LLNLDistribution11/Numerical/Lib/FFT.py
--- orig.LLNLDistribution11/Numerical/Lib/FFT.py Thu Apr 1 19:14:45 1999
+++ LLNLDistribution11/Numerical/Lib/FFT.py Tue Jul 13 19:37:36 1999
@@ -44,9 +44,19 @@
def real_fft(a, n=None, axis=-1):
return _raw_fft(a.astype(Numeric.Float), n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)
-def inverse_real_fft(a, n=None, axis=-1):
+def real_inverse_fft(a, n=None, axis=-1):
+ if n == None: n = (Numeric.shape(a)[axis]-1)*2
+ return _raw_fft(a.astype(Numeric.Complex), n, axis, fftpack.rffti, fftpack.rfftb, _real_fft_cache)/n
+
+def hermite_fft(a, n=None, axis=-1):
+ if n == None: n = (Numeric.shape(a)[axis]-1)*2
+ return real_inverse_fft(Numeric.conjugate(a), n, axis)*n
+
+def hermite_inverse_fft(a, n=None, axis=-1):
if n == None: n = Numeric.shape(a)[axis]
- return _raw_fft(a.astype(Numeric.Float), n, axis, fftpack.rffti, fftpack.rfftb, _real_fft_cache)/n
+ return Numeric.conjugate(real_fft(a, n, axis))/n
+
+inverse_real_fft = hermite_inverse_fft
def _raw_fft2d(a, s=None, axes=(-2,-1), function=fft):
a = Numeric.asarray(a)
@@ -58,7 +68,7 @@
return _raw_fft2d(a,s,axes,fft)
def inverse_fft2d(a, s=None, axes=(-2,-1)):
- return _raw_fft2d(a, s, axes, inverse_fft)
+ return _raw_fft2d(a, s, axes, inverse_fft)
def real_fft2d(a, s=None, axes=(-2,-1)):
return _raw_fft2d(a, s, axes, real_fft)
@@ -73,5 +83,134 @@
print inverse_fft2d (fft2d( [(0, 1), (1, 0)] ) )
print real_fft2d([(0,1),(1,0)] )
print real_fft2d([(1,1),(1,1)] )
+
+
+ import sys
+
+ oosq2 = 1.0/Numeric.sqrt(2.0)
+
+ toler = 1.e-10
+
+ p = Numeric.array(((1, 1, 1, 1, 1, 1, 1, 1),
+ (1, oosq2, 0, -oosq2, -1, -oosq2, 0, oosq2),
+ (1, 0, -1, 0, 1, 0, -1, 0),
+ (1, -oosq2, 0, oosq2, -1, oosq2, 0, -oosq2),
+ (1, -1, 1, -1, 1, -1, 1, -1),
+ (1, 0, 0, 0, 0, 0, 0, 0),
+ (0, 0, 0, 0, 0, 0, 0, 0)))
+
+ q = Numeric.array(((8,0,0,0,0,0,0,0),
+ (0,4,0,0,0,0,0,4),
+ (0,0,4,0,0,0,4,0),
+ (0,0,0,4,0,4,0,0),
+ (0,0,0,0,8,0,0,0),
+ (1,1,1,1,1,1,1,1),
+ (0,0,0,0,0,0,0,0)))
+
+ def cndns(m):
+ return Numeric.maximum.reduce(abs(m).flat)
+
+ try:
+ junk = hermite_fft
+ new = 1
+ except NameError:
+ new = 0
+
+ # Tests for correctness.
+## Somewhat limited since
+## p (and thus q also) is real and hermite, and the dimension we're
+## testing is a power of 2. If someone can cook up more general data
+## in their head or with another program/library, splice it in!
+
+ print "\nCorrectness testing dimension -1."
+
+ sys.stdout.write("fft: ")
+ sys.stdout.flush()
+ P = fft(p)
+ if cndns(P-q) / cndns(q) > toler:
+ print "inaccurate"
+ else:
+ print "OK"
+
+ sys.stdout.write("real_fft: ")
+ sys.stdout.flush()
+ RP = real_fft(p)
+ npt = p.shape[-1]
+ rpt = npt/2 + 1
+ qr = q[:,:rpt]
+ if cndns(RP-qr) / cndns(qr) > toler:
+ print "inaccurate"
+ else:
+ print "OK"
+
+ if new:
+ sys.stdout.write("real_inverse_fft: ")
+ sys.stdout.flush()
+ if cndns(real_inverse_fft(q, npt)-p) / cndns(p) > toler:
+ print "inaccurate"
+ else:
+ print "OK"
+ else:
+ sys.stdout.write("inverse_real_fft: ")
+ sys.stdout.flush()
+ hp = inverse_real_fft(q)
+ pr = p[:,:rpt]
+ if cndns(hp-pr) / cndns(pr) > toler:
+ print "inaccurate"
+ else:
+ print "OK"
+
+ # now just test consistency
+
+ for dim in range(len(p.shape)):
+ print "\nConsistency testing dimension %d, length %d."%(dim, p.shape[dim])
+
+ sys.stdout.write("fft/inverse_fft: ")
+ sys.stdout.flush()
+ P = fft(p, None, dim)
+ Q = inverse_fft(P, None, dim)
+ if cndns(Q-p) / cndns(p) > toler:
+ print "inconsistent"
+ else:
+ print "OK"
+
+ sys.stdout.write("fft/real_fft: ")
+ sys.stdout.flush()
+ RP = real_fft(p, None, dim)
+ npt = p.shape[dim]
+ rpt = npt/2 + 1
+
+ P = Numeric.take(P, range(rpt), dim)
+ if cndns(RP-P) / cndns(RP) > toler:
+ print "inconsistent"
+ else:
+ print "OK"
+
+ sys.stdout.write("inverse_fft/inverse_real_fft: ")
+ sys.stdout.flush()
+ hp = inverse_real_fft(q, npt, dim)
+ Q = inverse_fft(q, None, dim)
+ Q = Numeric.take(Q, range(rpt), dim)
+ if cndns(hp-Q) / cndns(hp) > toler:
+ print "inconsistent"
+ else:
+ print "OK"
+
+ if new:
+ sys.stdout.write("real_fft/real_inverse_fft: ")
+ sys.stdout.flush()
+ if cndns(real_inverse_fft(RP, npt, dim)-p) / cndns(p) > toler:
+ print "inconsistent"
+ else:
+ print "OK"
+
+ sys.stdout.write("inverse_real_fft/hermite_fft: ")
+ sys.stdout.flush()
+ if cndns(hermite_fft(hp, npt, dim)-q) / cndns(q) > toler:
+ print "inconsistent"
+ else:
+ print "OK"
+
+
if __name__ == '__main__': test()
diff -Naur orig.LLNLDistribution11/Numerical/Src/fftpackmodule.c LLNLDistribution11/Numerical/Src/fftpackmodule.c
--- orig.LLNLDistribution11/Numerical/Src/fftpackmodule.c Thu Apr 1 19:14:46 1999
+++ LLNLDistribution11/Numerical/Src/fftpackmodule.c Tue Jul 13 19:37:02 1999
@@ -108,7 +108,7 @@
PyObject *op1, *op2;
PyArrayObject *data, *ret;
double *wsave, *dptr, *rptr;
- int npts, nsave, nrepeats, i;
+ int npts, nsave, nrepeats, i, rstep;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, PyArray_DOUBLE, 1, 0);
@@ -117,6 +117,7 @@
data->dimensions[data->nd-1] = npts/2+1;
ret = (PyArrayObject *)PyArray_FromDims(data->nd, data->dimensions, PyArray_CDOUBLE);
data->dimensions[data->nd-1] = npts;
+ rstep = (ret->dimensions[ret->nd-1])*2;
if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1)
goto fail;
@@ -132,11 +133,11 @@
dptr = (double *)data->data;
for (i=0; i<nrepeats; i++) {
- memcpy((char *)rptr, dptr, npts*sizeof(double));
+ memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
rfftf(npts, rptr+1, wsave);
rptr[0] = rptr[1];
rptr[1] = 0.0;
- rptr += npts+2;
+ rptr += rstep;
dptr += npts;
}
PyArray_Free(op2, (char *)wsave);
@@ -161,12 +162,10 @@
int npts, nsave, nrepeats, i;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) return NULL;
- data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, PyArray_DOUBLE, 1, 0);
+ data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, PyArray_CDOUBLE, 1, 0);
if (data == NULL) return NULL;
npts = data->dimensions[data->nd-1];
- data->dimensions[data->nd-1] = npts/2+1;
- ret = (PyArrayObject *)PyArray_FromDims(data->nd, data->dimensions, PyArray_CDOUBLE);
- data->dimensions[data->nd-1] = npts;
+ ret = (PyArrayObject *)PyArray_FromDims(data->nd, data->dimensions, PyArray_DOUBLE);
if (PyArray_As1D(&op2, (char **)&wsave, &nsave, PyArray_DOUBLE) == -1)
goto fail;
@@ -177,17 +176,16 @@
goto fail;
}
- nrepeats = PyArray_SIZE(data)/npts;
+ nrepeats = PyArray_SIZE(ret)/npts;
rptr = (double *)ret->data;
dptr = (double *)data->data;
for (i=0; i<nrepeats; i++) {
- memcpy((char *)rptr, dptr, npts*sizeof(double));
- rfftb(npts, rptr+1, wsave);
- rptr[0] = rptr[1];
- rptr[1] = 0.0;
- rptr += npts+2;
- dptr += npts;
+ memcpy((char *)(rptr+1), (dptr+2), (npts-1)*sizeof(double));
+ rptr[0] = dptr[0];
+ rfftb(npts, rptr, wsave);
+ rptr += npts;
+ dptr += npts*2;
}
PyArray_Free(op2, (char *)wsave);
Py_DECREF(data);