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()