Re: [Numpydiscussion] FFT usage / consistency
Hi all, Stefan, thank you very much for your quick answer. This was an obvious silly mistake. Now the first function does what it should. Still, me and my colleagues can't make any sense of what happens in the second example. The retransformed function is identical to the original one, but the Fouriertransform doesn't have anything to do with what it should be mathematically. Could you or someone else please have a look at it? Thanks again, Felix
I have to correct myself. The function test_fft1() still does not work, it just looked good in the plots, but at a closer look, Re(IFFT) is close to zero and is far from matching the exact IFT. So it seems FFT(IFFT(f)) == IFFT(FFT(f)) == 1 (if done right ;), but I just cannot reproduce the exact (I)FT. Felix
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.
2008/7/29 Felix Richter <felix@physik3.unirostock.de>:
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.
Felix, Do your answers differ from the theory by a constant factor, or are they completely unrelated? Stéfan
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: UTF8 """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/((xdx)**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 xaxis # 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 humanreadable 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()
Hi Felix, I quickly copypasted 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?) Ciao, / / .o. // ..o / / ANS ooo
I quickly copypasted 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
On Tue, Jul 29, 2008 at 8:56 AM, Felix Richter <felix@physik3.unirostock.de
wrote:
I quickly copypasted 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.
You can apply a linear phase shift to the transformed data, i.e., multiply by something of the form exp(ixn), where x depends on where you want the center and n is the index of the transformed data point. This effectively rotates the original data. Or you can just rotate the data. If the data is not symmetric you are always going to have complex components. What exactly are you trying to do? I mean, what is the original problem that you are trying to solve by this method? Chuck
Felix Richter <felix <at> physik3.unirostock.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 regularlyspaced 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
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.
I'd strongly recommend "The Fourier Transform and its Applications" by Bracewell, if that helps. James.
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 Lorentzshaped. 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
participants (6)

Charles R Harris

Felix Richter

Hans Meine

James Turner

Neil MartinsenBurrell

Stéfan van der Walt