fft and listconvolve: 
exercise fft functions in mathematica, determine proper scaling.

program put together by Chip Durfee, CSM
latest revision: 4 Feb 2014	fixed factor of 2π in normalization of FFT’s, normalized convolution.

t-ω domain

gridding definitions and notes

The FFT algorithm itself does not make use of the labels for the grids used to sample the function. But given a sampling in one domain, we can construct the grid in the other domain. 

To use the FFT algorithm, we need the number of points in the array to be a power of 2. 
With an even number of points, the pulse can’t be exactly centered on the grid if one the points is at the peak of the pulse. Here, I’ve chosen the gridding so that tLst[[ npts/2 +1]] = 0

The choice of dω=((2π)/(npts dt)) ensures the scaling of the transform is correct. The sampling in the t-domain often needs to be adjusted so that the transform in ω-domain fits within the overall range of the grid without spilling over, or dω might need to be decreased to get better spectral resolution. A larger span in ω space is obtained by finer sampling in t space: npts dω = 2 π/dt. More resolution in ω space requires a larger temporal span. 
 
In this case, the center frequency is 0 and ωLst[[ npts/2+1 ]] = 0
If there is a carrier, e.g. Exp[-i (Subscript[ω, 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 (Subscript[ω, 0]): 
ωLst=Table[ω0 +i dω,{i,-(npts/2-1),npts/2}];

dt=1;
npts=256 4;
tLst=Table[i dt,{i,-npts/2,npts/2-1}];
dω=((2 π)/(npts dt));
ωLst=Table[i dω,{i,-npts/2,npts/2-1}];

tLst[[npts/2+1]]

create a test Gaussian pulse profile, τp is 1/e half-width
Later, we want to compare the numeric and analytic result, so we can also calculate the FT of the function using our convention for the FT.

Clear[τp]

f[t_]:=(E^(-((((t/(τp))))^2)));
F[ω_]=FourierTransform[f[t],t,ω,FourierParameters->{1,1}]

We will also want to check normalization to ensure energy is conserved during the transform. Notice that Parseval’s theorem says 
(∫_(-∞)^∞)(((|f(t)|))^2)dt = (1/(2π))(∫_(-∞)^∞)(((|F(ω)|))^2)dω
Notice the factor of 1/2π on the integral in ω space.

Integrate[((f[t])^2),{t,-∞,∞},Assumptions->τp>0]
Integrate[(1/(2π))((F[ω])^2),{ω,-∞,∞},Assumptions->τp>0]

do the sampling of the function
note that for most functions, Mathematica allows you to put the list directly into the argument of the function to produce the sampled list.

τp=5;
fLst=f[tLst];
Print["direct list plot of points in array"];
ListPlot[fLst,PlotRange->All]
Print["list plot with t-scale on x-axis"];
ListLinePlot[Transpose[{tLst,fLst}],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;
energyT=Total[((Abs[fLst])^2)]dt; (* shortcut *)
Print["energy in time domain: ", N[energyT]];

FFT calculation: list rotation for correct phase shifts

Although result should be real, there are imaginary parts that show up. If the gridding is chosen so that tLst[[ npts/2 +1 ]] = 0, then these imaginary components are small. Also note that alternating values of the output array have opposite sign. 

The complex character of the transform has to do with the shift theorem: a change of position of the signal in the grid by N dt is like a shift in time. This in turn results in a phase shift in ω space: exp[ i ω N dt ]

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

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]

The FFT algorithm internally assumes the 0 point of the array is positioned at the end point. We can prerotate the array, then do the transform, and the phase shift goes away.

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

What this rotation has done is to place the t=0 point as the first element of the array for the signal. Here we do the same rotate on tLst to show this:

tLstRotated=RotateRight[tLst,npts/2];
tLstRotated[[1]]

Now we have zero phase (a real spectrum) when we should, but the resulting signal is also rotated on the grid in ω space.

ωLstRotated=RotateRight[ωLst,npts/2];
ωLstRotated[[1]]

so we can rotate the result back to recenter the spectrum

FLst=RotateLeft[Fourier[RotateRight[fLst,(npts/2)]],(npts/2)];
Print["direct list plot of points in array"];
ListPlot[FLst,Joined->True,PlotRange->All]
Print["list plot with ω-scale on x-axis"];
ListPlot[Transpose[{ωLst,FLst}],Joined->True,PlotRange->All]

check grid scaling and energy conservation

If we compare the output of the FFT and compare it to the analytic result, we see that the results are the same, but are scaled differently.

FωSampledLst=F[ωLst];
Print["red=numeric, blue=analytic"];
ListLinePlot[{Transpose[{ωLst,FLst}],Transpose[{ωLst,FωSampledLst}]},PlotRange->All,PlotStyle->{Red,Blue}]

it looks like these are the same width, but the vertical scaling is off. 
Show this explicitly by normalizing each plot:

ListLinePlot[{Transpose[{ωLst,(FLst/(Max[Abs[FLst]]))}],Transpose[{ωLst,((FωSampledLst)/(Max[FωSampledLst]))}]},PlotRange->All,PlotStyle->{Red,Blue}]

Remember that the FFT algorithm doesn’t know or care about what the grids are. The Fourier[ ] and InverseFourier[ ] commands internally conserve energy. 
To check energy conservation, note that without dt and dω, the total energy in the arrays are the same

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

To get a result for the energy that accounts for the time grid, we define total energy in time domain with dt so energy is independent of sampling rate. 
the summing is the numeric equivalent of integrating over these spaces.

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

A more compact way to obtain this is:

energyT=Total[((Abs[fLst])^2)]dt //N

since we know that the sum without dt and dω is the same, we will get the same result by adding in extra factors as shown. Remember Parseval’s theorem has 1/2π in front of the frequency domain integral.

energyω=(1/(2π))(2π(dt/(dω)))Total[((Abs[FLst])^2)]dω

We can finally set up a definition of the transform function that accounts for the energy scaling factor.

FLst=(Sqrt[2π(dt/(dω))])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
energyω=(1/(2π))Total[((Abs[FLst])^2)]dω;
Print["energy in ω domain: ", energyω];

recheck the comparison with analytic

Print["red=numeric, blue=analytic"];
ListLinePlot[{Transpose[{ωLst,FLst}],Transpose[{ωLst,FωSampledLst}]},PlotRange->All,PlotStyle->{{Red,Thick,Dashed},Blue}]

now check inverse transform: do exact opposite

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

Note that the amplitude scaling factors (Sqrt[2π(dt/(dω))]) and (Sqrt[((dω)/(2π dt))]) cancel out when you do (FFT^(-1)){ FFT{ f[t] } }, so if we don’t care about having the proper amplitude scaling in ω space, these factors aren’t necessary. 

Also note that the factor of (Sqrt[2π]) in the amplitude scaling factor is directly connected to our convention for the FT and Parseval’s theorem.

finally, define new function to take all this into acct

make use of the connection between dt and dω:
dω=((2π)/(npts dt)), so (Sqrt[2π(dt/(dω))])=(Sqrt[2π(dt/((2π)/(npts dt)))])=(Sqrt[npts (dt^2)])

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

Print["red=numeric, blue=analytic"];
ListLinePlot[{Transpose[{ωLst,FFT[fLst,npts,dt]}],Transpose[{ωLst,FωSampledLst}]},PlotRange->All,PlotStyle->{{Red,Thick,Dashed},Blue}]

Print["red=numeric, blue=analytic"];
ListLinePlot[{Transpose[{ωLst,Log[((Abs[FFT[fLst,npts,dt]])^2)]}]},PlotRange->All,PlotStyle->{{Red},Blue}]

ListConvolve

product in time domain results in a convolution in ω domain

convolution theorem: 
transform of product in the time domain:
FT{f(t) g(t)}=(1/(2 π))(∫_(-∞)^∞)F(ω')G(ω-ω')dω'

Clear[f,F,tau1,tau2,g,G,h1,h2]

f[t_]:=(E^(-((((t/tau1)))^2)));
F[ω_]=FourierTransform[f[t],t,ω,FourierParameters->{1,1}];
g[t_]:=((((t/tau2)))^2)(E^(-((((t/tau2)))^2)));
G[ω_]=FourierTransform[g[t],t,ω,FourierParameters->{1,1}];
h1[ω]=FourierTransform[f[t]g[t],t,ω,FourierParameters->{1,1}]
h2[ω]=(1/(2 π))Integrate[F[ωp]G[ω-ωp],{ωp,-∞,∞},Assumptions->{tau1>0,tau2>0}]
FullSimplify[((h1[ω])/(h2[ω])),Assumptions->{tau1>0,tau2>0}]

these are the same

Now do numeric convolution.  The npts/2+1 makes the rotation of the output arrays line up. The final 0 in the ListConvolve[ ] pads the second array with 0's rather than doing the convolution cyclically. 

The scale factor of dω/2π with the ListConvolve are there to make sure the convolution integral is done w.r.t. ω, and to have the 1/2π for the convolution theorem.

dt=1;
npts=256*8;
tLst=Table[i dt,{i,-npts/2,npts/2-1}];
dω=2 π/(npts dt);
ωLst=Table[i dω,{i,-npts/2,npts/2-1}];
ωMin=Min[ωLst];ωMax=Max[ωLst];

tau1=10;tau2=40;
fLst=f[tLst];
gLst=g[tLst];
FLst=FFT[fLst,npts,dt];
GLst=FFT[gLst,npts,dt]; 
H1Lst=FFT[fLst gLst,npts,dt];
H2Lst=ListConvolve[FLst,GLst,npts/2+1,0]dω/(2π);

Print["red = FFT{ f g}, blue = F⊗G"];
ListPlot[{Transpose[{ωLst,Abs[H1Lst]}],Transpose[{ωLst,Abs[H2Lst]}]},Joined->True,PlotRange->{{ωMin,ωMax},All},PlotStyle->{{Red,Thick,Dashed},Blue}]

product in frequency domain

convolution theorem: 
transform of product in the frequency domain:
(FT^(-1)){F(ω) G(ω)}=(∫_(-∞)^∞)F(t')G(t-t')dt'

F[ω_]:=(E^(-((((ω/(dω1))))^2)));
f[t_]=InverseFourierTransform[F[ω],ω,t,FourierParameters->{1,1}]
G[ω_]:=(E^(-((((ω/(dω2))))^2)));
g[t_]=InverseFourierTransform[G[ω],ω,t,FourierParameters->{1,1}];
h1[t]=InverseFourierTransform[F[ω]G[ω],ω,t,FourierParameters->{1,1}]
h2[t]=Integrate[f[tp]g[t-tp],{tp,-∞,∞},Assumptions->{dω1>0,dω2>0}]
FullSimplify[((h1[t])/(h2[t])),Assumptions->{dω1>0,dω2>0}]

these are the same

Now do numeric convolution. The scale factors in front of the ListConvolve are there empirically. The npts/2+1 makes the rotation of the output arrays line up. The final 0 in the ListConvolve[ ] pads the second array with 0's rather than doing the convolution cyclically. 

May want to do the convolution manually in a loop to check where the scale factors come from. 
Everything seems the same in this version except you pass in the freq spacing dw instead of dt. Same normalization factor.

dt=.05;
npts=256*2;
tLst=Table[i dt,{i,-npts/2,npts/2-1}];
dw=2 π/(npts dt);
wLst=Table[i dw,{i,-npts/2,npts/2-1}];
tMin=Min[tLst];tMax=Max[tLst];
wMin=Min[wLst];wMax=Max[wLst];
Print["tMin: ",tMin//N," tMax: ",tMax//N];
Print["wMin: ",wMin//N," wMax: ",wMax//N];

dω1=1;dω2=4;
FLst=F[wLst];GLst=g[wLst];
fLst=InverseFFT[FLst,npts,dt];gLst=InverseFFT[GLst,npts,dt]; 
H1Lst=InverseFFT[FLst GLst,npts,dt];
H2Lst=dt ListConvolve[fLst,gLst,npts/2+1,0];

ListPlot[{Transpose[{tLst,Abs[H1Lst]}],Transpose[{tLst,Abs[H2Lst]}]},Joined->True,PlotRange->{{tMin,tMax},All},PlotStyle->{{Red,Thick,Dashed},Blue}]

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}];

create a test Gaussian beam profile, wx is 1/e half-width
Later, we want to compare the numeric and analytic result, so we can also calculate the FT of the function using our convention for the FT.

Clear[wx]

g[x_]:=(E^(-((((x/wx)))^2)));
G[fx_]=FourierTransform[g[x],x,fx,FourierParameters->{0,-2π}];
FullSimplify[G[fx],Assumptions->wx>0]

We will also want to check normalization to ensure energy is conserved during the transform. Notice that Parseval’s theorem says 
(∫_(-∞)^∞)(((|g(x)|))^2)dx = (∫_(-∞)^∞)(((|G((Subscript[f, x]))|))^2)d(Subscript[f, x])
Notice there is no factor of 1/2π on the integral in (Subscript[f, x]) space like there is in the t, ω domains

Integrate[((g[x])^2),{x,-∞,∞},Assumptions->wx>0]
Integrate[((G[fx])^2),{fx,-∞,∞},Assumptions->wx>0]

sample the function on the grid

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

FullSimplify[G[0],Assumptions->wx>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,G[fxLst]}],Transpose[{fxLst,Chop[GLst]}]},Joined->True,PlotRange->{{-1,1},All},PlotStyle->{{Red,Thick,Dashed},Blue}]

both width and height look right