fft demo: exercise fft functions in mathematica, determine proper scaling.

program put together by Chip Durfee, CSM

t-ω domain

create the arrays for labeling the time and frequency axes
In this case, the center frequency is 0. If there is a carrier, e.g. Exp[-i ω0 t], that is suppressed so that the temporal profile just represents the envelope, then the resulting shift in the frequency spectrum is accounted for by centering the frequency axes at ω0.

dt=1;
npts=256;
tLst=Table[i dt,{i,-npts/2,npts/2-1}];
dw=2 π/(npts dt);
wLst=Table[i dw,{i,-npts/2,npts/2-1}];

create a test Gaussian pulse profile, tau is 1/e half-width

f[t_]:=(E^(-((((t/tau)))^2)));
tau=20;
fLst=f[tLst];
Print["direct list plot of points in array"];
ListPlot[fLst,Joined->True,PlotRange->All]
Print["list plot with t-scale on x-axis"];ListPlot[Transpose[{tLst,fLst}],Joined->True,PlotRange->All]

calculate total energy by just summing Abs value squared of all points, mult by dt

energyT=(∑_(i=1)^npts)((Abs[(Subscript[fLst, [[i]]])])^2)dt;
Print["energy in time domain: ", N[energyT]];

calculate FFT of list
Although result should be real, there are small imaginary parts that show up.
Also note that alternating values of the output array have opposite sign.

FLst=Fourier[fLst];
Print["sample elements from array: ",FLst[[1]]," ",FLst[[2]] ];

(* ListPlot[FLst,Joined->True,PlotRange->All]*)

we can get rid of small imag parts using Chop[ ]
we're still left with points that alternate in sign

ListPlot[Chop[FLst],Joined->True,PlotRange->All]

this sign alternation doesn't affect the amplitude.

ListPlot[Abs[FLst],Joined->True,PlotRange->All]

If we prerotate the array, then do the transform, the phase alternation goes away

FLst=Fourier[RotateRight[fLst,(npts/2)]];
ListPlot[Chop[FLst],Joined->True,PlotRange->All]

rotate the list to recenter the array

so the method is to rotate first, FFT, then rotate back
this way the arrays you use are centered in their arrays when you use them

FLst=RotateLeft[Fourier[RotateRight[fLst,(npts/2)]],(npts/2)];
ListPlot[Chop[FLst],Joined->True,PlotRange->All]
ListPlot[Transpose[{wLst,Chop[FLst]}],Joined->True,PlotRange->All]

check energy conservation
without dt and dw, sums are equal

energyT=(∑_(i=1)^npts)((Abs[(Subscript[fLst, [[i]]])])^2);
Print["energy in time domain: ", N[energyT]];energyW=(∑_(i=1)^npts)((Abs[(Subscript[FLst, [[i]]])])^2);
Print["energy in w domain: ", energyW];

define total energy in time domain with dt so energy is independent of sampling rate

energyT=(∑_(i=1)^npts)((Abs[(Subscript[fLst, [[i]]])])^2)dt //N
energyW=(dt/dw)(∑_(i=1)^npts)((Abs[(Subscript[FLst, [[i]]])])^2)dw

therefore, set up a definition of the transform function to account for this factor

FLst=(Sqrt[(dt/dw)])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
energyW=(∑_(i=1)^npts)((Abs[(Subscript[FLst, [[i]]])])^2)dw;
Print["energy in w domain: ", energyW];

now check inverse transform: do exact opposite

fLstInv=(Sqrt[(dw/dt)]) RotateLeft[InverseFourier[RotateRight[FLst,(npts/2)]],(npts/2)];
ListPlot[Transpose[{tLst,fLstInv}],Joined->True,PlotRange->All]

finally, define new function to take all this into acct

Clear[fLst,npts,dt,FLst];
FFT[fLst_,npts_,dt_]:=(Sqrt[((npts (dt^2))/(2 π))])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
InverseFFT[FLst_,npts_,dt_]:=(Sqrt[((2 π)/(npts (dt^2)))])RotateLeft[InverseFourier[RotateRight[FLst,npts/2]],npts/2];

x, fx version

this uses conventions found in Goodman for 2D or 1D spatial Fourier transforms.

dx=1;
npts=256;
xLst=Table[i dx,{i,-npts/2,npts/2-1}];
dfx=1/(npts dx);
fxLst=Table[i dfx,{i,-npts/2,npts/2-1}];

g[x_]:=(E^(-((((x/wx)))^2)));
wx=20;
gLst=g[xLst];
Print["list plot with x-scale on x-axis"];ListPlot[Transpose[{xLst,gLst}],Joined->True,PlotRange->All]

calculate total energy by just summing Abs value squared of all points, mult by dt

energyX=(∑_(i=1)^npts)((Abs[(Subscript[gLst, [[i]]])])^2)dx;
Print["energy in spatial domain: ", N[energyX]];

calculate FFT of list
Although result should be real, there are small imaginary parts that show up.
Also note that alternating values of the output array have opposite sign.

rotate the list to recenter the array

GLst=RotateLeft[Fourier[RotateRight[gLst,(npts/2)]],(npts/2)];
ListPlot[Transpose[{fxLst,Chop[GLst]}],Joined->True,PlotRange->All]

check energy conservation
without dx and dfx, sums are equal

energyX=(∑_(i=1)^npts)((Abs[(Subscript[gLst, [[i]]])])^2);
Print["energy in spatial domain: ", N[energyX]];energyFX=(∑_(i=1)^npts)((Abs[(Subscript[GLst, [[i]]])])^2);
Print["energy in fx domain: ", energyFX];

define total energy in time domain with dt so energy is independent of sampling rate

energyX=(∑_(i=1)^npts)((Abs[(Subscript[gLst, [[i]]])])^2)dx;
energyFX=(dx/dfx)(∑_(i=1)^npts)((Abs[(Subscript[GLst, [[i]]])])^2)dfx;

therefore, define transform to account for this factor

GLst=(Sqrt[(dx/dfx)])RotateLeft[Fourier[RotateRight[gLst,npts/2]],npts/2];
energyFX=(∑_(i=1)^npts)((Abs[(Subscript[GLst, [[i]]])])^2)dfx;
Print["energy in fx domain: ", energyFX];

now check inverse transform: do exact opposite

gLstInv=(Sqrt[(dfx/dx)]) RotateLeft[InverseFourier[RotateRight[GLst,(npts/2)]],(npts/2)];
ListPlot[Transpose[{xLst,gLstInv}],Joined->True,PlotRange->All]

finally, define new function to take all this into acct

Clear[gLst,npts,dx,GLst];
FFT[gLst_,npts_,dx_]:=(Sqrt[npts (dx^2)])RotateLeft[Fourier[RotateRight[gLst,npts/2]],npts/2];
InverseFFT[GLst_,npts_,dx_]:=(Sqrt[(1/(npts (dx^2)))])RotateLeft[InverseFourier[RotateRight[GLst,npts/2]],npts/2];

check amplitude in transform domain to see if it is consistent with integral

Clear[wx];
G[fx_]=FourierTransform[g[x],x,fx,FourierParameters->{0,-2 π}]

the 1/e width of this exponential is 1/(π wx) in fx space

G[0]

for this function, this evaluates to a peak of wx (Sqrt[π])

dx=.1;
npts=256;
xLst=Table[i dx,{i,-npts/2,npts/2-1}];
dfx=1/(npts dx);
fxLst=Table[i dfx,{i,-npts/2,npts/2-1}];

wx=1;
Print["peak height should be: ",wx (Sqrt[π])//N];
Print["1/e width should be:   ",1/(π wx)//N];

g[x_]:=Exp[-(((x/wx))^2)];
gLst=g[xLst];
GLst=FFT[gLst,npts,dx];

ListPlot[Transpose[{fxLst,((Chop[GLst])/(Max[Re[GLst]]))}],Joined->True,PlotRange->{{-0.6`,0.6`},All}]

Max[Re[GLst]]

1/E //N

both width and height look right