# 2-Dimensional NFourierTransform

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}]

• Are you sure in the function? I has sharp long non-decaying tails: img. It it very difficult to numerically integrate this function. Sep 30 '14 at 13:40
• 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/…
– bcp
Oct 1 '14 at 5:01
• 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.
– bcp
Oct 3 '14 at 10:06

\$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}]


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.

DeveloperPackedArrayQ /@ {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


The tails are reduced with the window

F W // ArrayPlot


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}]


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}]


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


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

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

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