I am trying to use the maxima fft routines to do some numerical
solutions of PDEs that I derived using maxima. I want to discretize the
space, and then use the fft to calculate derivatives and second
derivatives etc in space.
There are some irregularities that I'm experiencing in the whole process.
For example if I have data y in an array, and I take its FFT, and then
multiply each coefficient by 2*%pi*%i*k/N and do an inverse transform I
would expect to get the derivative... but I seem to get the derivative
multiplied by %pi * (1-epsilon) where epsilon isn't something like
(N-1)/N or anything that would be obvious if i had a bug.
So i eliminated the factor of %pi in my coefficients, and that gets me
very close, but I'm not sure why, it doesn't make much sense except that
the fft implementation has some normalization issues. In any case,
here's my fftderiv function, can anyone tell if there's something wrong
in my function or if FFT is perhaps just "losing" energy in the signal
due to numerical problems?
fftderiv(vals) :=
block([n:length(vals),coefs:[],pi:ev(%pi,numer)],
array(fftdReal,n-1),
array(fftdImag,n-1),
vals:float(vals),
coefs:float(makelist(2*k/(n),k,0,n-1)),
print(coefs),
fillarray(fftdReal,vals),fillarray(fftdImag,[0.0]),
fft(fftdReal,fftdImag),
fillarray(fftdReal,listarray(fftdReal) * coefs),
fillarray(fftdImag,listarray(fftdImag) * -coefs),
ift(fftdImag,fftdReal),
/* multiplying by %i is the same as reversing the role of real/imag
and multiplying by -1 for the i^2 */
listarray(fftdImag));