5
$\begingroup$

This question is related to "Numerical Fourier transform of a complicated function" from 2012.

I have pretty much the same problem: I have a function $f(x)$ which is complicated to evaluate, and want to calculate $$ F(t) = \frac{\mathrm e^t}{\pi} \int_0^\infty f(x)\, \mathrm e^{itx}\; \mathrm dx . $$ I have sampled $f(x)$, for example data see here. From these points I want to obtain $F(t)$. The problem is: I have tried to do this in 3 different ways, and each gives slightly different results. I want to know which method is most reliable, and how to best improve accuracy.

If you want to try this yourself, save the example data in a file "samples.txt" and load it:

samples = ToExpression /@
    Import[FileNameJoin[{NotebookDirectory[], "samples.txt"}], "Data"];
samplePlot = ListPlot[{
    Cases[samples, {x_, f_} -> {x, Re[f]}],
    Cases[samples, {x_, f_} -> {x, Im[f]}]},
  PlotTheme -> "Detailed", PlotLegends -> None, PlotRange -> {{0, 200}, Full},
  AspectRatio -> 1/5, ImageSize -> Full, PlotMarkers -> {Automatic, Tiny}]

Method 1: NFourierTransform

My first attempt was to simply use NFourierTransform. For this, I interpolated my sample data, continuing with zero outside of the sampled domain. (Note that I convert to an integral $\int_{-\infty}^\infty$ by setting $f(-x)=f(x)^\ast$.)

maxX = 3000;

realParts = Cases[samples, {x_, f_} -> {x, Re[f]}];
imagParts = Cases[samples, {x_, f_} -> {x, Im[f]}];
intRe = Interpolation[realParts];
intIm = Interpolation[imagParts];
interpolation[x_] := 
    Which[Abs[x] >= maxX, 0, x >= 0, intRe[x] + I*intIm[x], True, intRe[-x] - I*intIm[-x]];

Needs["FourierSeries`"];
method1 = Range[-3, 0.75, 0.1] // 10^# & //
  Table[{t, E^t / (2*\[Pi]) * Re@NFourierTransform[
    interpolation[x], x, t, FourierParameters -> {1, 1}]}, {t, #}] &;

plot1 = ListLogLogPlot[method1,
  PlotTheme -> "Detailed", PlotLegends -> None, ImageSize -> Large, 
  PlotStyle -> Black, PlotMarkers -> {Automatic, Small}, PlotRange -> {10^-4, 20}]

Note that this gives negative values $F(t)$ for the largest $t$-values, which should be impossible.

Method 2: DFT

Next I thought, I have discrete data points, I should maybe try a discrete Fourier transformation. Taking some code from this answer, I did the following:

dx = 0.1;
numpoints = maxX/dx;
upsampled = Range[0, maxX, dx] // Table[interpolation[x], {x, #}] &;

method2 = Re@Fourier[upsampled, FourierParameters -> {1, 1}]*dx //
    RotateRight[#, Floor[(numpoints - 1)/2]] & //
    Transpose@{Range @@ ({-Floor[(numpoints - 1)/2], Ceiling[(numpoints + 1)/2], 1}/(maxX)), #} & //
    Cases[#, {t_, v_} -> {2*\[Pi]*t, E^(2*\[Pi]*t)/\[Pi]*v}] & //
    Interpolation;

plot2 = LogLogPlot[method2[t], {t, 0.001, 5}, PlotRange -> {10^-4, 20},
  PlotTheme -> "Detailed", PlotLegends -> None, ImageSize -> Large]

Note that the original samples are not evenly spaced, I am using the interpolation from above to create more sample points.

Method 3: Fitting

I have found a simpler function which fits the sample data reasonably well:

samplesMirrored = Cases[samples, {t_, f_} -> {-t - 1, Re[f]}] \[Union]
    Cases[samples, {t_, f_} -> {t, Im[f]}];
fitfun[x_] := A1/(1 + B1 + I*x) + A2/(1 + B2 + I*x) + (A3*2!)/(1 + B3 + I*x)^3;
fitfunMirrored[x_] := Evaluate@If[x < -1/2, 
  Evaluate@ComplexExpand@Re[fitfun[-x - 1]], Evaluate@ComplexExpand@Im[fitfun[x]]];
fit = FindFit[
  samplesMirrored,
  {fitfunMirrored[x], B1 > 0, B2 > 0, B3 > 0},
  {{A1, 6.6}, {B1, 15.3}, {A2, 2.7}, {B2, 5.2}, {A3, 1.2}, {B3, 3.5}},
  x];

Show[Plot[{Re@fitfun[x] /. fit, Im@fitfun[x] /. fit}, {x, 0, 20}, 
  PlotTheme -> "Detailed", PlotLegends -> None, ImageSize -> Large], samplePlot]
Show[Plot[{Re@fitfun[x] /. fit, Im@fitfun[x] /. fit}, {x, 1000, 3000},
  PlotTheme -> "Detailed", PlotLegends -> None, ImageSize -> Large], samplePlot]

The whole "mirrored"-business just deals with the fact that FindFit can't handle complex functions. The fit function can be Fourier transformed analytically:

method3[t_] := A1*E^(-B1*t) + A2*E^(-B2*t) + A3*t^2*E^(-B3*t) /. fit;
plot3 = LogLogPlot[method3[t], {t, 0.001, 5}, PlotRange -> {10^-4, 20},
  PlotTheme -> "Detailed", PlotLegends -> None, ImageSize -> Large, PlotStyle -> Orange]

The results look as follows (black dots: NFourierTransform, blue: DFT, orange: Fit): Plot of all F(t)

As you can see, the methods agree for intermediate values of $t$ and I would be certain enough that the result there is correct. However, I need $F(t)$ ideally for the whole range of $t$-values that is displayed here. What do you think: is one of these three methods more trustworthy than the others? If not, what would be the best way to increase the accuracy? (I did experiment a bit with taking more sample points, but it didn't immediately seem to help too much.)

$\endgroup$
5
$\begingroup$

First, notice that you have two types of discrepancies: the one at smaller values (left of the graph where the blue curve is under the dots) and the one at larger values (right of the graph where the blue curve is over the dots).

The discrepancy for large values of t is related to the assumptions that each method makes about the behavior in between points in samples. Think of t as being the frequency. A related concept is the Nyquist–Shannon sampling theorem. In any case, the main idea is that because the distance between points in your sampling is not infinitesimally small then you can reproduce the Fourier transform only up to certain frequencies.

If the transformation $f\to F$ did not have the factor $e^t$ in front of it, I would estimate the cutoff frequency using $.1\ t = π$. I might be off by a factor, but this is an estimate of the order of magnitude anyways. Because of the extra $e^t$ factor in front of the integral, I would do

ProductLog[π/(.1)]
(*2.52219*)

That seems to be the value where the blue curve has its minimum.

Since your data is compatible with all three solutions, you pick one using what you think is the behavior of $f(x)$ for $x$ between the sample points. It you think it should be boring (almost straight lines) then pick the answer with $F(t)$ going to zero as $t$ goes to infinity (orange line). If you think that the Fourier transform (that is the integral) does not go to zero faster than $e^{-t}$, then pick the dots.

Regarding Mathematica, all three method are giving you (basically) the same answer for the range of values of t where an answer can be found without further assumptions.

By the way, using

ProductLog[ π / Max[ samples[[All, 1]] ] ] // N

to estimate the onset of discrepancy for low frequencies, I got 0.0010461. The order of magnitude is correct but the blue curve began to diverge from the other two at values of t slightly larger.

$\endgroup$

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.