Is there a way to make Fourier's Re and Im components separately track the results of FourierTransform? Is there a similar way to do the same for InverseFourier and InverseFourierTransform? I have noticed other questions with answers for the absolute value of the transforms, but they oscillate wildly when dealing with the real or imaginary part individually.
$\begingroup$
$\endgroup$
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
Yes, this is possible. Code for the transforms is given after the graphs and the residual checks.
Here's a demonstration in the frequency domain with the real part as blue and imaginary as orange. The lines are from the continuous Fourier transform, F(w), and the dots are from the discrete transform.
f@x_=Sign@x Sqrt@x Exp@-Abs@x; (*Sinc[x-1] is also good*)
F@w_=FourierTransform[f@x,x,w];
SetOptions[Plot,ImageSize->600,PlotRange->{{-5,5},All}];
SetOptions[ListPlot,PlotStyle->PointSize->Medium];
Show[Plot[{Re@F@w,Im@F@w},{w,-5,5}],
ListPlot@Thread[{{#,Re@#2},{#,Im@#2}}&@@@
First[dftMost@Table[{x,f@x},{x,aa=-100,bb=100,.1}]]]]
Here's a demonstration in the space (or time) domain. Colors are the same as above. The lines are from (what would be) the continuous inverse Fourier transform of F[w] to get f[x], and the points are from the discrete version.
Show[Plot[{Re@f@x,Im@f@x},{x,-5,5}],
ListPlot@Thread[{{#,Re@#2},{#,Im@#2}}&@@@
First@idft[Table[{w,F@w},{w,-32,32,.03}],{-32,32,aa}]]]
Showing that IDFT of the DFT of f has near-zero residual vs f at all points:
Table[{x,f@x},{x,aa,bb,.1}];
First[idft@@dftMost@%]-Most@%//Flatten//Abs//Max
(*2.9708966550739*^-13*)
Showing that DFT of the IDFT of F has low residual vs F at all points:
Table[{w,F@w},{w,-31.38,31.42,3/100}];
First[dft@@idft[%,{-31.38,31.42,aa}]]-%//Flatten//Abs//Max//N
(*0.005002388915432476*)
Note: Because we directly typed in the frequency limits of -31.38 to 31.42 for the sake of clarity, there is not a direct match with aa. aa controls the phase offset of the spatial domain, a.k.a. the shift of the dots along the x axis. I have separately verified that if I pull the frequency limits from the 2nd list output of the DFT, or just use the relevant formulas, then the match with aa is almost exact and then this residual is again on the order of 10^-13 like the one above.
Here is the code for the discrete Fourier transform:
Options@dft=Options@dftMost={FourierParameters->{0,1}};
(*main dft algorithm*)
dft[f_?VectorQ,{xini_?NumericQ,xfin_?NumericQ},
opts:OptionsPattern[]]:=
Module[{a,b,dF,ini,NN=Length@f,k,dk,L=xfin-xini,dx},
{a,b}=OptionValue@FourierParameters;
dk=2Pi/L;dx=L/NN;ini=1-Ceiling[NN/2];
dF=Fourier[f Exp[I b Range[0,NN-1] dx dk ini],opts];
{Transpose@{k=Range[ini,Floor[NN/2]]dk,
Sqrt[(2Pi/NN)^(a-1) Abs@b] dF Exp[I b k xini]dx},
{First@k,Last@k,xini},opts}
]
(*consume output of idft*)
dft[m_?MatrixQ,l:(_?VectorQ):{},opts:OptionsPattern[]]:=
With[{len=Length@l},dft[m[[All,2]],
If[len==2,l,{m[[1,1]],m[[-1,1]]}],opts]]
(*consume output of a table of data, forcing peroidicity*)
dftMost[m_?MatrixQ,opts:OptionsPattern[]]:=
dft[Most@m[[All,2]],{m[[1,1]],m[[-1,1]]},opts]
Here is the code for the inverse discrete Fourier transform:
Options@idft={FourierParameters->{0,1}};
(*main idft algorithm*)
idft[f_?VectorQ,{kini_?NumericQ,kfin_?NumericQ,
xini:(_?NumericQ):0},opts:OptionsPattern[]]:=
Module[{a,b,dF,ini,NN=Length@f,dk,L,x,dx},
{a,b}=OptionValue@FourierParameters;
dk=(kfin-kini)/(NN-1);L=2Pi/dk;dx=L/NN;
dF=InverseFourier[f Exp[-I b Range[0,NN-1]dk xini],opts];
{Transpose@{x=Range[0,NN-1]dx+xini,
Sqrt[(2Pi/NN)^(-a-1) Abs@b] dF Exp[-I b x kini]dk},
{xini,xini+L},opts}
]
(*consume output of dft or of a table of data*)
idft[m_?MatrixQ,l:_List:{},opts:OptionsPattern[]]:=
With[{len=Length@l},idft[m[[All,2]],
If[len==3||len==2,l,{m[[1,1]],m[[-1,1]],0}],opts]]
As a bonus, here is an example of checking with Parseval's theorem to make sure the same energy is present before and after Fourier transformation. Most is only needed on temporal or spatial series data to prevent double counting: if a function is truly periodic, the answers at zero radians and 2 Pi radians of the fundamental waveform would be the same, so the last point should not be counted. Frequency data, however, is usually already clipped appropriately.
Table[{x,f@x},{x,aa,bb,.1}];
parseval/@{Most@%,First@dftMost@%}
(*{0.49833666138305915,0.4983366613830593}*)
Parseval check code:
parseval[l_List/;MatrixQ@l]:=Total[Abs[l[[All,2]]]^2]*
Subtract@@MinMax@l[[All,1]]/(1-Length@l)
For an extra bonus, these utility functions may be used to specify appropriate min, max, and fundamental step frequencies so that the endpoints of the range and the step fit together in exact integer multiples.
num@exp_=2^exp;
wFund[num_,wMax_]=wMax/Floor[num/2];
period[num_,wMax_]=2Pi/wFund[num,wMax];
wSampling[num_,wMax_]=num wFund[num,wMax];
wIni[num_,wMax_]=(1-Ceiling[num/2]) wFund[num,wMax];
wFin[_,wMax_]=wMax;
wMinMax[num_,wMax_]={wIni[num,wMax],wFin[num,wMax]};
wMinMaxStep[num_,wMax_]=Append[
wMinMax[num,wMax],wFund[num,wMax]]
(*{(wMax*(1-Ceiling[num/2]))/Floor[num/2],
wMax,wMax/Floor[num/2]}*)
Illustrating exact multiples of fundamental frequency step. Here I've used powers of two for the number of points, but it isn't required.
TableForm[Table@@{Join[{e,num@e},wMinMaxStep[num@e,wm]/
wFund[num@e,wm]],{e,0,10}},TableHeadings->{None,
{e,num ==2^e,wIni/wFund,wFin/wFund,"wFund"/wFund}}]
Illustrating near-zero residual for DFT of IDFT with the utility functions. Compare with the earlier DFT of IDFT transformation above. Note that this corresponds to starting the direct space integration at x==0.
Table@@{{w,F@w},Prepend[wMinMaxStep[2000,32],w]};
First[dft@@idft[%,wMinMax[2000,32]]]-%//Flatten//Abs//Max//N
(*2.2887833992611187*^-16*)