I learned a few things in the meantime: In my installation, NumPy uses fftpack_lite while SciPy uses FFTW3. There are more test cases in SciPy which all pass. So I am confirmed my problem is a pure usage problem. One thing I was confused about is the fact that even if I calculate the function over a certain interval, I cannot tell FFT which interval this is, it will instead assume [0...n]. So actually I did not transform a Lorentz function centered at zero but rather centered at 500. Unfortunately, this solves only half of my problem, because I still cannot reproduce the exact FT. I'll ask for that on the SciPy list, this now seems more appropriate.

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

I quickly copy-pasted and ran your code; it looks to me like the results you calculated analytically oscillate too fast to be represented discretely. Did you try to transform different, simpler signals? (e.g. a Gaussian?) Yes, I run into the same problem.

Since the oscillation frequency is given by the point around which the function is centered, it would be good to have it centered around zero. The FFT assumes the x axis to be [0..n], so how should I do this? The functions I have to transform later won't be symmetrical, so the trick abs(fftdata) is not possible. Felix

Felix Richter <felix <at> physik3.uni-rostock.de> writes:

...

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!

This is not voodoo, this is signal processing, which is itself harmonic analysis. Just because the Fast Fourier Transform is fast doesn't mean that this stuff is easy. You seem to be looking for a simple relationship between the Fourier Transform (an integral transform from L^2(R) -> L^2(R)) of a function f and the Discrete Fourier Transform (a linear transformation from R^n to R^n) of the vector of f sampled at regularly-spaced points. Such a simple relationship does not exist. That is why you found no such examples on the internet. The closest you might come is to study the surprisingly cogent explanation at http://en.wikipedia.org/wiki/Fourier_analysis of the differences between the various types of Fourier analysis. Remember that the DFT (as implemented by an FFT algorithm) is *not* an approximation to the Fourier transform, but rather a streamlined way of computing the coefficients of a Fourier series of a particular periodic function (that contains a finite number of Fourier modes). Rather than look for errors in the scaling factors or errors in your code, I think that you should try to expand your understanding of the (subtly) different types of Fourier representations. -Neil

Thanks for all your comments. It's definitely time to read a good book now. My original problem is a convolution of two complex functions given as samples over quite different intervals with different n. The imaginary part of one of these functions is Lorentz-shaped. I thought it might be good to resample them in the frequency domain, then multiply and transform back. For the resampling I have to make sure the two resulting frequency axises are equivalent/physically meaningful. Also, of course, I wanted to understand what the NumPy/SciPy routines do and how to use them correctly. Now I'll try to just resample in the time domain, transform without looking at the result, then blindly multiply and transform back. This seems to work, but I'll have to find a different testcase so I can make sure the results are trustworthy. Felix