FFT algorithm completely botches phase, help!
sleepingawake86
I am trying to use igor pro to analyze some interferograms to find the amplitude and the phase. I have done a few trials(sine wave, Gaussian*sine, etc) and igor completely botches the phase, and I don't know why, or how to fix it.
For instance, I take a function Sin(w_0*t+phi) and perform the FFT on it and it gives me a spike at w_0 for the amplitude, as expected, but when I look at the phase as a function of w it looks random, not constant at phi.
If I do the DFT by hand on paper I get the right answer for the phase, so its not an issue with the DFT being discrete, its something with the FFT implementation, and I have no clue about how igor goes about it, or how to fix it. Help?
FT[f(x-a)] -> F(s)*exp(-i*2*pi*a*s)
where FT denotes the Fourier Transform and F(s) is the FT of the function f(x).
Using the latest version of IGOR I executed the following:
Variable phi=pi/4
make/n=(1024) ddd=cos(w0*x+phi)
FFT/Dest=fff /out=1 ddd
Display fff
ModifyGraph cmplxMode=1
Display fff
ModifyGraph cmplxMode=2
At this point you can set various phi values and observe the changes in the second graph. Feel free to send me an example experiment (support@wavemetrics.com) if you can reproduce any randomness or results that are not consistent with theory.
A.G.
WaveMetrics, Inc.
November 2, 2011 at 03:15 pm - Permalink
Using the DFT formula and calculating by hand when I do the FT of Sin(w0*t+phi) and taking the phase I get the right answer for w0, and junk elsewhere, which I expected. When I use the FFT in igor, the phase graph has a jump discontinuity at w0 and does not reach the correct phi. what is going on?
November 3, 2011 at 02:46 pm - Permalink
In order to investigate we need to reproduce what you are doing. So either post some commands that illustrate the problem or post an experiment or send commands or an experiment to WaveMetrics support.
November 3, 2011 at 02:55 pm - Permalink
Using a set of model data for a two level atom interferogram I produced a graph of phase vs freq for my DFT, an FFT, and the actual model.
The DFT almost matches perfectly, the FFT looks terrible. Attached is the experiment and graph.
November 3, 2011 at 03:29 pm - Permalink
November 3, 2011 at 10:56 pm - Permalink
I'm glad you provided an experiment as it illustrates a number of problems not all of which are relevant to your question but they might be worth mentioning:
1. Your scaling is all over the place. The FFT requires some kind of wave scaling and it is not possible to provide it via an x-wave. Since your x-wave appeared to be monotonic I just provided scaling using the command:
2. Your DFT code should be performed in DP not SP so your Make commands should be adjusted accordingly.
3. I suspect that your x-wave already reflects an extra factor of 2; otherwise the arguments for the cos() and sin() in your DFT function should have the form i*2*pi/N. This may be connected with (1) above but I wanted to mention it for clarity.
4. For various reasons it is recommended to use atan2() as opposed to atan() in the last line of the loop.
5. It is very instructive to compute and plot the magnitude of the signal. Your graph has conveniently omitted the start of the data where it is actually significant. The portion of the data that you are displaying is between 4 and 5 orders of magnitude smaller which should make us pause and consider the consequences.
6. I do not know where you got your phasemodel wave so it is not obvious to me why this should be the "correct" result. This is precisely why we need an example where one can evaluate the Fourier Transform analytically and thereby determine if there is an issue with the phase.
A.G.
WaveMetrics, Inc.
November 4, 2011 at 05:22 pm - Permalink
If I look at your data (wave0 and wave1) I can approximately fit it with a product of a Gaussian and a sine function. Something along the line of
f(x)=a+b*exp(-(x-c)^2/d)*sin(pi*x*e) with (approximately)
a=-2.3e-5
b=4.9
c=0.5
d=19.5
e=0.749
The exact values are not important for the following argument as is the basic shape of the function. If you now apply the convolution theorem of Fourier Transform you'd expect the transform to be the convolution of a Gaussian with a delta function. The Gaussian in the transform domain should have a width that is a scaled reciprocal of d and the convolution should result in nothing more than this Gaussian centered at the frequency corresponding to e. Once you get outside this Gaussian envelope, you have to be extremely cautious in computing anything e.g., the atan() of the ratio of the imaginary to the real part of the signal.
Now to explain the difference between the FFT and DFT calculations I suggest that you run the following example that is somewhat simpler than your functional form but illustrates the issue well:
•display yyy
•fft /dest=foo/out=1 yyy
•matrixop/o refoo=real(foo)
•matrixop/o imfoo=imag(foo)
•display imfoo
The result of the phase here may appear to be random but is mostly an issue of sampling if one remembers the shift theorem of Fourier transforms. To see that:
•fft /dest=foo/out=1 yyy
•matrixop/o refoo=real(foo)
In other words, the shift introduced in the "x" domain results in an oscillating phase which is sampled at a frequency other than that of its oscillation thus giving you the impression of a random phase.
I hope this helps,
A.G.
WaveMetrics, Inc.
November 5, 2011 at 04:09 am - Permalink
Ok so I believe we are on the same page now, the question now is given that in real life we can't get perfect sampling rates, etc. How do I take what igors FFT gives me, and get the right answer out?
November 9, 2011 at 03:11 pm - Permalink
For the record, it should be clear that IGOR computes the correct FFT.
In my first response on Nov 2, I suggested that you look at consequences of the shift theorem but you did not think I answered your question. My reply of Nov 4 should have shown you exactly how the shift applies to your data and why you see a difference between the FFT and your own procedure.
I could be wrong but your last response quoted above suggests that this is still not clear so let me be more specific:
If you want to account for the precise x origin you can either rotate the input wave (check the documentation for the Rotate operation) prior to executing the FFT or you could apply the FFT to your unmodified input wave (using FFT/out=1) and post multiply the complex result of the FFT by the complex conjugate of the phase that you expect based on the shift theorem. Both methods should produce identical results and agree with your manual computation.
Before you try either method I suggest that you take a quick look at the "FFT Swapping Demo" which you can find under File Menu->Example Experiments->Analysis. The demo illustrates both "rotation" and multiplication by a complex phase.
I hope that after all these examples you can convince yourself that IGOR is not "botching the phase"
A.G.
WaveMetrics, Inc.
November 10, 2011 at 01:34 pm - Permalink