FFT

Preliminaries

(* F1 on this for a tutorial *)
tutorial/FourierTransforms

(* define function to plot several things side by side *)
Clear[plots];
plots[plotlist_]:=GraphicsRow[plotlist,0,ImageSize->250Length[plotlist]]

Ringing

(* length of periodic sequence *)
Clear[Nseq];
Nseq=(2^4)

(* top hat function at center of interval *)
Clear[centersquare];
centersquare=Table[If[Abs[i]≤(Nseq/8),1,0],{i,-Quotient[Nseq-1,2],Quotient[Nseq,2]}];

(* top hat function at edge of interval *)
Clear[edgesquare];
edgesquare=Table[If[i≤1+Quotient[Nseq,8]||i>Nseq-Quotient[Nseq,8],1,0],{i,1,Nseq}];

plots[{ListLinePlot[centersquare],ListLinePlot[edgesquare]}]

(* The Fourier transforms of the centered and edge top hats differ only in phase *)
plots[{ ListLinePlot[Abs[Fourier[centersquare]],PlotRange->All],ListLinePlot[Re[Fourier[edgesquare]],PlotRange->All]}]

Mathematica's convention for FFT

(* mathematica's discrete Fourier transform is the complex list (Subscript[(OverTilde[a]), 0]), (Subscript[(OverTilde[a]), 1]), ..., (Subscript[(OverTilde[a]), [(N/2)]]), (Subscript[(OverTilde[a]), [(N/2)]+1]) == (Subscript[(OverTilde[a]), -N+[(N/2)]+1]), ..., (Subscript[(OverTilde[a]), -2]), (Subscript[(OverTilde[a]), -1]) *)
(* if the function (Subscript[a, n]) is real in real space, then (Subscript[(OverTilde[a]), 0]) is real, (Subscript[(OverTilde[a]), [(N/2)]]) is real for even N, and (Subscript[(OverTilde[a]), -m]) = (SuperStar[(Subscript[(OverTilde[a]), m])]), so the second half of the list is redundant, being the complex conjugate of the first half *)
Print[Subscript[(OverTilde[a]),0]," = ",Fourier[centersquare][[1]]];Print[Subscript[(OverTilde[a]),Quotient[Nseq+1,2]-1]," to ",Subscript[(OverTilde[a]),Quotient[Nseq,2]+1]," = ",Table[Fourier[centersquare][[i]],{i,Quotient[Nseq+1,2],Quotient[Nseq,2]+2}]];
Print[Subscript[(OverTilde[a]),Nseq-1]," = ",Fourier[centersquare][[Nseq]]];

(* mathematica defines its discrete Fourier transform to include a factor of (Sqrt[N]), so that FT∘(FT^\[Conjugate]) is the identity *)
Chop[Fourier[Conjugate[Fourier[centersquare]]]-centersquare]

Windowing remedies ringing caused by sharp truncation

(* define a "broad" window of the form (((1-(x^2)))^(n ))in real space, that extends over frac of the interval, and goes to zero at each end *)window[n_,frac_]:=Module[{x},Table[x=(1/frac)((2i)/Nseq);If[Abs[x]≤1,(((1-(x^2)))^n),0],{i,-Quotient[Nseq-1,2],Quotient[Nseq,2]}]];

ListLinePlot[window[1,1],PlotRange->{0,All}]

(* the Fourier transform of the broad window is a convolution kernel spiked around zero frequency *)
ListLinePlot[Take[Abs[Fourier[window[1,1]]],Quotient[Nseq,2]+1],PlotRange->All]

(* same plot, plotted logarithmically *)
ListLinePlot[Take[Log[10,Abs[Fourier[window[1,1]]]],Quotient[Nseq,2]+1],PlotRange->All]

Power spectrum of random data, and its shot noise

(* set expected mean of Poisson distribution *)
Clear[ranmean];
ranmean=100;

(* list of data drawn from a Poisson distribution distribution *)
Clear[randat];
randat=Table[RandomInteger[PoissonDistribution[ranmean]],{i,1,(2^10)}];

ListLinePlot[randat,PlotRange->All]

(* actual shot noise of random data = actual mean of the data; expected shot noise = expected mean = ranmean *)
Clear[randatmean];
randatmean=(1/(Length[randat]))(∑_(i=1)^(Length[randat]))randat[[i]];
N[randatmean]

(* square root of the power spectrum of the random data, compared to the square root of the shot noise *)
ListLinePlot[{Take[Abs[Fourier[randat-randatmean]],Quotient[Length[randat],2]+1],Table[(Sqrt[randatmean]),{Quotient[Length[randat],2]+1}]},PlotRange->All]