Continuous Fourier Transform VS FFT

[evad1089]

I have about 40 tabs open on this right now and something important is slipping my grasp. I know this has been covered a million and a half times, but for some reason I cannot seem to find a straight answer (or more probably realize and understand it when I see it).

When I take the Continuous Fourier transform of a function, I get back a very nice magnitude and phase.

When I try the exact same thing with a FFT I get back a very nice magnitude, but my phase is highly erratic. Most of the time, it is just oscillates in between two values every point (eg -\pi and \pi or does this on a slope. I am using the atan2 function.

If I unwrap this, the phase will go to ridiculously huge numbers.

Can someone explain if why this happens and more importantly if and how I can make the discrete phase match the continuous phase?

 

[torquil]

I'm not sure exactly what you are doing, so this is just a shot in the dark. But if I recall correctly, sometimes the data structure used in the output of a FFT algorithm arranges the Fourier coefficients in a different order than what you would expect. In Matlab, I think you have to use the "fftshift" or "ifftshift" function to correct for this? Which FFT software are you using?

 

[evad1089]

I have tried MatLab, Maxima, and Mathematica. Currently I am using Octave (free MatLab style program). I also have been using fftshift.

Here is what I mean with the phase. The below two picutures are of the phase of a Gaussian Beam. The first is the continuous transform and the second is the FFT.

[PLAIN]http://blackbricksoftware.com/custom/Screenshot-1.png

[PLAIN]http://blackbricksoftware.com/custom/Screenshot.png

There is a highly erratic behivour in the FFT phase.

 

[torquil]

I'm just guessing since I haven't seen your code, but perhaps most of those values on the x-axis correspond to Fourier coefficient magnitudes that are extremely small? If so, those phase values might just represent numerical inaccuracies?

I mean, it doesn't make a whole lot sense to plot the phase of a complex Fourier coefficient if the magnitude of that coefficient is comparable to the numerical accuracy.

 

[AlephZero]

I think you either misunderstood what the input to your FFT routine should be (for example, did you tell the FFT routine correctly whether the input was real or complex data), and/or you misunderstood how the the output is formatted (for example the data for the first and last Fourier coefficients are often stored differently to the rest, because their imaginary parts are always zero, and omitting the zeros means you need 2N real numbers to represent the FFT of N data points, not 2N+2 real numbers two of which are always 0).

Start with a simpler test data set where you know EXACTLY what the output should be. For example, generate a single sine or cosine wave, so all the output should be 0 except for one Fourier coefficient and one phase angle. Then try something like a square or triangle wave, where the exact FFT is well known.

 

[evad1089]

For anyone who stumbles across this in the future, my problem turned out to be two fold.

First, floating point rounding errors in calculating the phase cause a lot of jumps of pi.

Second, the FFT is in forward time, ie t>0. This is a phase shift. There are charts and such, depending on the type of function transformed, on how to correct this and obtain a centered transformation. This is not as useful as one would think because the rounding errors, the phase will still be off from the continuous transformation. Take a rectangular pulse centered around zero. Its transform is a sinc function. The phase of a sinc function jumps from zero to pi (or vice versa) every time the magnitude is zero. Unless the sinc is totally centered around zero, the phase does not have this behavior. It is almost impossible to transform a dataset and divide out the time shift term without having some rounding errors.

The solution I found for approximating the Fourier transform of a centered functions was to use the discrete time Fourier transform. Basically, instead of an integral, it used a Riemann sum. of each of the points multiplied by the e^{-i t w} to create a continuous function that will approximate the Fourier transform. I used Maxima CAS to do this, but Maple or Mathematica would work.

Here is my code:
dftp(points,k) := block([n:1,N:length(points),dftpoints:[]],
while n<=N do (
dftpoints: endcons(points[n]*exp(-%i*2*%pi*n*k/N),dftpoints),
n: n+1
),
dftpoints/N*2)$
dftps(points,k) := block([n:1,N:length(points),sum:0],
while n<=N do (
sum : sum+points[n]*exp(-%i*2*%pi*(n)*k/N),
n: n+1
),
sum/N*2)$
dftpso(points,k,o) := block([n:1,N:length(points),sum:0],
while n<=N do (
sum : sum+points[n]*exp(-%i*2*%pi*(n-o)*k/N),
n: n+1
),
sum/N*2)$