OK, I've been hacking on the fft stuff (share/numeric/fft.lisp)
& I have some stuff I'd like to commit.
I'm airing this plan in case somebody wants to talk me out of it.
What's new:
* fft and fft_inverse accept lists as arguments
(in addition to named and unnamed arrays)
* return value type is same as whatever was supplied
* non-float numeric values OK
* test script
What's changed:
* renamed ift to fft_inverse
* fft and fft_inverse do not modify their arguments
* fft and fft_inverse take 1 argument (complex),
not 2 (separate real and imaginary parts)
* attempted to clarify documentation
If for whatever reason a user wants to supply only
Lisp arrays (real and imaginary parts) and have them
modified in place, those functions are accessible
from Maxima as ?fft and ?fft_inverse.
(Modifying arguments is a mess, so is separating
real and imaginary parts by hand. But it's still there
if that's your pleasure.)
There was some discussion about the name for the inverse
transform and from what I recall, fft_inverse wasn't any less
popular than the alternatives that were floated.
I reworked the documentation to try to be concise and accurate.
There was a suggestion to explain the result in terms of
frequencies; I didn't incorporate any of that as it didn't seem
make the situation any clearer.
I've appended the current version of the documentation
for fft for your enjoyment.
Hope this helps,
Robert Dodier
PS.
(%i1) ? fft
-- Function: fft (<x>)
Computes the complex fast Fourier transform. <x> is a list or
array (named or unnamed) which contains the data to transform.
The number of elements must be a power of 2. The elements must be
literal numbers (integers, rationals, floats, or bigfloats) or
symbolic constants, or expressions `a + b*%i' where `a' and `b'
are literal numbers or symbolic constants.
`fft' returns a new object of the same type as <x>, which is not
modified. Results are always computed as floats or expressions `a
+ b*%i' where `a' and `b' are floats.
The discrete Fourier transform is defined as follows. Let `y' be
the output of the transform. Then for `k' from 0 through `n - 1',
y[k] = (1/n) sum(x[j] exp(-2 %i %pi j k / n), j, 0, n - 1)
When the data <x> are real, real coefficients `a' and `b' can be
computed such that
x[j] = sum (a[k] * cos (2*%pi*j*k / n) + b[k] * sin
(2*%pi*j*k / n), k, 0, n/2)
with
a[0] = realpart (y[0])
b[0] = 0
and, for k from 1 through n/2 - 1,
a[k] = realpart (y[k] + y[n - k])
b[k] = imagpart (y[n - k] - y[k])
and
a[n/2] = realpart (y[n/2])
b[n/2] = 0
`load(fft)' loads this function.
See also `fft_inverse' (inverse transform), `recttopolar', and
`polartorect'.
Examples:
Real data.
(%i1) load (fft) $
(%i2) fpprintprec : 4 $
(%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $
(%i4) L1 : fft (L);
(%o4) [0.0, - 1.811 %i - .1036, 0.0, .6036 - .3107 %i, 0.0,
.3107 %i + .6036, 0.0, 1.811 %i - .1036]
(%i5) L2 : fft_inverse (L1);
(%o5) [1.0, 2.168L-19 %i + 2.0, 7.525L-20 %i + 3.0,
4.256L-19 %i + 4.0, - 1.0, - 2.168L-19 %i - 2.0,
- 7.525L-20 %i - 3.0, - 4.256L-19 %i - 4.0]
(%i6) lmax (abs (L2 - L));
(%o6) 3.545L-16
Complex data.
(%i1) load (fft) $
(%i2) fpprintprec : 4 $
(%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $
(%i4) L1 : fft (L);
(%o4) [0.5, .3536 %i + .3536, - 0.25 %i - 0.25,
0.5 - 6.776L-21 %i, 0.0, - .3536 %i - .3536, 0.25 %i - 0.25,
0.5 - 3.388L-20 %i]
(%i5) L2 : fft_inverse (L1);
(%o5) [1.0 - 4.066E-20 %i, 1.0 %i + 1.0, 1.0 - 1.0 %i,
- 1.008L-19 %i - 1.0, 4.066E-20 %i - 1.0, 1.0 - 1.0 %i,
1.0 %i + 1.0, 1.947L-20 %i + 1.0]
(%i6) lmax (abs (L2 - L));
(%o6) 6.83L-17
Computation of sine and cosine coefficients.
(%i1) load (fft) $
(%i2) fpprintprec : 4 $
(%i3) L : [1, 2, 3, 4, 5, 6, 7, 8] $
(%i4) n : length (L) $
(%i5) x : make_array (any, n) $
(%i6) fillarray (x, L) $
(%i7) y : fft (x) $
(%i8) a : make_array (any, n/2 + 1) $
(%i9) b : make_array (any, n/2 + 1) $
(%i10) a[0] : realpart (y[0]) $
(%i11) b[0] : 0 $
(%i12) for k : 1 thru n/2 - 1 do
(a[k] : realpart (y[k] + y[n - k]),
b[k] : imagpart (y[n - k] - y[k]));
(%o12) done
(%i13) a[n/2] : y[n/2] $
(%i14) b[n/2] : 0 $
(%i15) listarray (a);
(%o15) [4.5, - 1.0, - 1.0, - 1.0, - 0.5]
(%i16) listarray (b);
(%o16) [0, - 2.414, - 1.0, - .4142, 0]
(%i17) f(j) := sum (a[k] * cos (2*%pi*j*k / n) + b[k] * sin
(2*%pi*j*k / n), k, 0, n/2) $
(%i18) makelist (float (f (j)), j, 0, n - 1);
(%o18) [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]