> Do your answers differ from the theory by a constant factor, or are
> they completely unrelated?
No, it's more complicated. Below you'll find my most recent, more stripped
down code.
- I don't know how to scale in a way that works for any n.
- I don't know how to get the oscillations to match. I suppose its a problem
with the frequency scale, but usage of fftfreq() is straightforward...
- I don't know why the imaginary part of the FFT behaves so different from the
real part. It should just be a matter of sin vs. cos.
Is this voodoo? ;-)
And I didn't find any example on the internet which tries just to reproduce an
analytic FT with the FFT...
Thanks for your help!
# coding: UTF-8
"""Test for FFT against analytic results"""
from scipy import *
from scipy import fftpack as fft
import pylab
def expdecay(t, dx, a):
return exp(-a*abs(t))*exp(1j*dx*t) * sqrt(pi/2.0)
def lorentz(x, dx, a):
return a/((x-dx)**2+a**2)
origfunc = lorentz
exactfft = expdecay
xrange, dxrange = linspace(0, 100, 2**12, retstep=True)
n = len(xrange)
# calculate original function over positive half of x-axis
# this serves as input to fft, make sure datatype is complex
ftdata = zeros(xrange.shape, complex128)
ftdata += origfunc(xrange, 50, 1.0)
# do FFT
fftft = fft.fft(ftdata)
# normalize
# but how exactly?
fftft /= sqrt(n)
# shift frequencies into human-readable order
fftfts= fft.fftshift(fftft)
# determine frequency axis
fftscale = fft.fftfreq(n, dxrange)
fftscale = fft.fftshift(fftscale)
# calculate exact result of FT for comparison
exactres = exactfft(fftscale, 50, 1.0)
pylab.subplot(211)
pylab.plot(xrange, ftdata.real, 'x', label='Re data')
pylab.legend()
pylab.subplot(212)
pylab.plot(fftscale, fftfts.real, 'x', label='Re FFT(data)')
pylab.plot(fftscale, fftfts.imag, '.', label='Im FFT(data)')
pylab.plot(fftscale, exactres.real, label='exact Re FT')
pylab.plot(fftscale, exactres.imag, label='exact Im FT')
pylab.legend()
pylab.show()
pylab.close()
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
