[Numpy-discussion] FFT usage / consistency

Felix Richter felix@physik3.uni-rostock...
Tue Jul 29 08:04:41 CDT 2008


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


More information about the Numpy-discussion mailing list