1
$\begingroup$

I want to apply Mathematica's DFT procedure, Fourier[], to a Gaussian. Why a Gaussian? Because I know the analytical result.

So I define

a = 0.5;
h[t_] := Exp[-(t^2/a^2)];

The analytic FourierTransform of h using FourierParameters -> {0, -2 Pi} is:

HAnalytic[v] := Sqrt[Pi]*a*Exp[-a^2*Pi^2*v^2];

Then I create time values:

dt = 1/100;
tList = Table[t, {t, 0, 20000 dt, dt}];

and I create sample data:

inputData = h[tList]

I then create the frequency values:

dv = 1/((Length[tList] - 1)*dt);
vList = dv*Table[i, {i, 0, Length[tList] - 1}];

and carry out the DFT calculation:

HDFT = Re[Fourier[inputData]];

Now I plot it:

Show[
 Plot[HAnalytic[v], {v, 0, 1000*dv}, PlotRange -> {{0, 1.5}, {-1, 1}},
   PlotStyle -> Green],
 ListPlot[{Thread[{vList, HDFT}][[1 ;; 500]]}]
 ]

enter image description here

So far no problem. Only the scaling is different.

BUT NOW MY QUESTION!

Why do I get

enter image description here

i.e. the Fourier-transformed values are mirrored around the frequency axis, when using

tListSym = Table[t, {t, -10000 dt, 10000 dt, dt}]

I mean h[tListSym] gives really kind of the whole function to the DFT, as h is symmetric around $t = 0$, whereas h[tList] reproduces only the right hand side of h.

If I look at the values of HDFT, I can see that the values are alternatively positive / negative.

Can someone explain to me what is going on?

$\endgroup$
  • $\begingroup$ Does this help ? dsp.stackexchange.com/questions/431/… $\endgroup$ – image_doctor Jun 17 '15 at 11:07
  • $\begingroup$ Not really... as I am not wondering about negative frequenies. I am wondering why the amplitude of the frequencies is alternating for t =ListSym = Table[t, {t, -10000 dt, 10000 dt, dt}]. $\endgroup$ – newandlost Jun 17 '15 at 12:24
  • $\begingroup$ @Hugh has recently written up an excellent primer on DFT in Mathematica. Maybe something in there might be of help to you. In particular, I think you might want to plot the Abs[] of your DFT values, rather than the Re[]. $\endgroup$ – MarcoB Jun 17 '15 at 14:45
  • $\begingroup$ Does not matter, as the fouriertransform of a real even function is a real even function. $\endgroup$ – newandlost Jun 17 '15 at 15:15
3
$\begingroup$

The answer is that Fourier does not know your time values and does not see an even function. Thus if we generate your new input data and plot we get

inputData2 = h[tListSym];
ListLinePlot[inputData2, PlotRange -> All]

Mathematica graphics

The problem is thus how to tell Fourier that the first half of your input is for negative x values. We could change the phase of the values created by Fourier or we can put your input into the usual Fourier form of positive values then negative values - a shifted version. Thus

nn = Length[tListSym]
inputData3 = 
Join[inputData2[[(nn - 1)/2 + 1 ;; -1]], 
     inputData2[[1 ;; (nn - 1)/2]]];
ListLinePlot[inputData3, PlotRange -> All]

Mathematica graphics

Thus we are starting with the positive values and then putting the negative values after. Now when we do the Fourier transform we get what you want.

HDFT2 = Fourier[inputData2];
HDFT3 = Fourier[inputData3];

Here I have not taken the real parts as you did. We can compare the two versions

ListPlot[{Re[HDFT2][[1 ;; 500]], Re[HDFT3][[1 ;; 500]]}, 
 PlotRange -> All]

Mathematica graphics

ListPlot[{Im[HDFT2][[1 ;; 500]], Im[HDFT3][[1 ;; 500]]}, 
 PlotRange -> All]

Mathematica graphics

The imaginary part of the shifted version is zero because we have constructed an even function. We now get what you want.

Show[
Plot[HAnalytic[v], {v, 0, 1000*dv}, 
     PlotRange -> {{0, 1.5}, {-1, 1}}, PlotStyle -> Green],
ListPlot[{Thread[{vList, Re[HDFT2]}][[1 ;; 500]]}, PlotStyle -> Blue],
ListPlot[{Thread[{vList, HDFT3}][[1 ;; 500]]}, PlotStyle -> Brown]
 ]

Mathematica graphics

On a secondary point I think you are off by one point in defining your frequency values. The increment should be 1/(Length[] dt).

I suppose a deeper question is why does Fourier do positive values then negative values? There is an assumption of a periodic function so it does not matter. However, it is often a nuisance as you have found.

With regard to getting the scaling correct you need to work out your version of Parseval's Theorem using FourierParameters.

$\endgroup$
  • $\begingroup$ I was hoping you might be induced to weigh in on this one :-) Let me also add a link to a post in the DSP stackexchange site dealing explicitly with the math behind the even/odd function problem, just for reference. (+1) $\endgroup$ – MarcoB Jun 17 '15 at 18:09
  • $\begingroup$ @Hugh Thanks you :) "I suppose a deeper question is why does Fourier do positive values then negative values?" Yes ... But I think no I understand it. You have to think about it like this: It is a N periodic function, hence 0 = N, AND the first value given to the DFT defines 0. This is the same as sying the function lives on a ring. Therefore for getting a even function we have fist hand it the positive and then the negative values. I will make a small drawing: drive.google.com/file/d/0B1ecQW5sX1LJLWhfUUt6ckktY1U/edit $\endgroup$ – newandlost Jun 18 '15 at 11:15
  • $\begingroup$ @newandlost Unfortunately I can't download or see your video -I will try when I get to another computer. However, I think I can see what is happening and I agree with your comment. There is also as solution that would involve adjusting the phase of the Fourier transform after the calculation i.e. making a transformation in time. Thus by writing Exp[- I 2 Pi f (t-t0)] in the Fourier integral choosing a t0 and then taking this phase outside the integral it would be possible to correct the spectrum. Glad I could help. $\endgroup$ – Hugh Jun 18 '15 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.