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):

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.)