Mathematica FourierSeries package contains the NFourierTransform function for calculating 1-D Fourier integral numerically.

What is the best approach for calculting 2D Fourier integral numerically? Should I simply use the NIntegrate? If so, what is the best Method?

The example of the transform I am interested in is the following:

FourierTransform[1/(x^2 - 2 y^2 - 5 + 5 I), {x, y}, {fx, fy}]
  • $\begingroup$ Are you sure in the function? I has sharp long non-decaying tails: img. It it very difficult to numerically integrate this function. $\endgroup$
    – ybeltukov
    Sep 30 '14 at 13:40
  • $\begingroup$ Unfortunately yes. This is the function I need to integrate. Besides, I am interested in the analytic solution as well: math.stackexchange.com/questions/952480/… $\endgroup$
    – bcp
    Oct 1 '14 at 5:01
  • $\begingroup$ Actually, I can reduce this integral to 1D Fourier transform by calculating the integral with respect to x analytically. However, the convergence of the y-integral will still be rather poor. $\endgroup$
    – bcp
    Oct 3 '14 at 10:06

Your function

$HistoryLength = 0;
a = -2;
b = 5 I - 5;
f[x_, y_] := 1/(x^2 + a y^2 + b);

have long and sharp tails

Plot3D[Re[f[x, y]], {x, -20, 20}, {y, -20, 20}, PlotPoints -> 200, 
 MaxRecursion -> 5, PlotRange -> All, AxesLabel -> {x, y}]

enter image description here

However there is a common procedure to calculate the Fourier transform numerically. It is tricky from the first sight but it is quite obvious if you apply this technique several times. Two main ideas:

  1. Use the discrete fast Fourier transform.
  2. Use a window function.

Let us discretize from -R to R with the step d over x and y

R = 100.;
d = 0.1;

x = Join[-Reverse@#, Rest@#] &@Range[0., R, d];
n = Length@x;
m = (n - 1)/2;
k = π/(d m) Range[0. - m, m];
X = ConstantArray[x, n];
Y = Transpose@X;

F = f[X, Y];
W = # Transpose@# &@ConstantArray[BlackmanWindow[x/R/2], n];

Capital letters mean 2D arrays of values. Here W is a window. It smoothly cuts of the function and reduces Gibbs oscillations. BlackmanWindow is a good choice, but you can choose many other window functions.

Developer`PackedArrayQ /@ {X, Y, F, W}

{True, True, True, True, True}

All our arrays are packed so everything will be calculated very fast.

Without the window our function have the long tails

F // ArrayPlot

enter image description here

The tails are reduced with the window

F W // ArrayPlot

enter image description here

The Fourier transform and the interpolation

fF = d^2 RotateRight[#, {m, m}] &@
     Fourier[#, FourierParameters -> {1, 1}] &@RotateLeft[F W, {m, m}];

res = ListInterpolation[Transpose[fF], {k, k}];

Plot3D[Im@res[kx, ky], {kx, -5, 5}, {ky, -5, 5}, PlotPoints -> 50, 
 MaxRecursion -> 3, PlotRange -> {-1, 1}, AxesLabel -> {kx, ky}]

enter image description here

Now we can compare our numerical result with the exact formula which I propose on math SE

exact[kx_, ky_] := 2 π BesselK[0, Sqrt[b (kx^2 + ky^2/a)]]/Sqrt[a]

Plot3D[Im@exact[kx, ky], {kx, -5, 5}, {ky, -5, 5}, PlotPoints -> 50, 
 MaxRecursion -> 3, PlotRange -> {-1, 1}, Exclusions -> None, AxesLabel -> {kx, ky}]

enter image description here

Plot[{Re@res[kx, 1], Re@exact[kx, 1]}, {kx, -3, 3}, 
 Exclusions -> None, PlotRange -> 2, ImageSize -> 500]

enter image description here

Plot[{Im@res[kx, 1], Im@exact[kx, 1]}, {kx, -3, 3}, 
 Exclusions -> None, PlotRange -> 2, ImageSize -> 500]

enter image description here

Everything are the same (except the sharp peaks due to the finite discretization)!

  • $\begingroup$ Very nice exposition, thanks! I seem to have however a problem using this for a couple of specific examples in Multidimensional DFT of function $\endgroup$
    – tks
    Jun 29 '17 at 7:55

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.