# Fourier transform in three dimensions

I have a four-dimensional array of the form $$A_{ijkl}=\{x_i,y_j,z_k,f(x_i,y_j,z_k)\}$$, where $$f$$ is some known function. I would like to compute the Fourier transform of $$f(x,y,z)$$ numerically. In order to do this, I apply

B = Fourier[A[[All,All,All,-1]]];


Infamously, Fourier returns the data with the axes organised in a "funny" way. So my question is, given the array $$\bf A$$ and $$\bf B$$, how do I reconstruct the four dimensional array $$\bf C$$, with $$C_{ijkl}=\{k^x_i,k^y_j,k^z_k,\mathcal{F}(k^x_i,k^y_j,k^z_k)\}$$ where $$\mathcal{F}$$ is the Fourier transform of $$f$$?

Example:

Take $$f(x,y,z)=e^{i x}+\frac{1}{4}e^{i y}+\frac{1}{2}e^{i z}$$. Take as sampling points \begin{align} x_i = -3.2+\frac{i}{10}\,,\quad i=1,\ldots,63 \\ y_j = -3.2+\frac{j}{10}\,,\quad j=1,\ldots,63 \\ z_k = -3.2+\frac{k}{10}\,,\quad k=1,\ldots,63 \end{align} How do I apply Fourier, and arrange the axis, so that my array $$\bf C$$ only has non-zero elements (within numerical error) $$\{1,0,0,1\}\,,\quad\{0,1,0,\frac{1}{4}\}\,,\quad\{0,0,1,\frac{1}{2}\}\,.$$

To get an idea about the axes in a multidimensionafourier transform, we make some simple examples:

For a 2D Fouriertransform where only the x or y components vary:

d = Table[{Sin[1.1 x]}, {x, 0, 2 Pi, Pi/2}, {y, 0, 2 Pi, Pi/2}];
Fourier[d] // Chop // MatrixForm
d = Table[{Sin[1.1 y]}, {x, 0, 2 Pi, Pi/2}, {y, 0, 2 Pi, Pi/2}];
Fourier[d] // Chop // MatrixForm It is obvious that the first axes belongs to x variations and the second axes to y variations. If we have x and y variations:

d = Table[{Sin[1.1 x y]}, {x, 0, 2 Pi, Pi/2}, {y, 0, 2 Pi, Pi/2}];
Fourier[d] // Chop // MatrixForm Now e.g. the element with indices {1,1} is the dc value. The element {1,2} belongs to x/y frequencies={0,1}. The part {2,2} belongs to frequencies: {1,1}. The part {3,2} belongs to frequencies: {2,1} e.t.c.

The 3 dimensional case is similar. The part with indices {i1,i2,i3} belongs to frequencies: {i1-1,i2-1,i3-1}

To get only the function value, you would write:

data[[All, All, All, -1]]


Here is an example:

n = 3;
data = Table[{x, y, z, x + y + z}, {x, n}, {y, n}, {z, n}];
Fourier[data[[All, All, All, -1]]] • I just amended my edits. My concern are the axes. How do I rebuild the actual Fourier transform from the output of Fourier? In one dimension, this is is rather trivial (we cut the output of Fourier in half, and do some shifts). In two dimensions, it gets less trivial but I found a post in Stack Exchange that does the trick . For three-dimensional data, I have not managed. Jun 30 at 9:31
• Look at my second answer Jun 30 at 10:26

In the end, the solution was not that hard. I do think, however, that Mathematica could make this easier for us and provide a function that also returns the axis in frequency space.

f[x_, y_, z_] := Exp[I x] + 1/4 Exp[I y] + 1/2 Exp[I z]
R = \[Pi];
d = 0.1;
x = Join[-Reverse@#, Rest@#] &@Range[0., R, d];
y = Join[-Reverse@#, Rest@#] &@Range[0., R, d];
z = Join[-Reverse@#, Rest@#] &@Range[0., R, d];
n = Length@x;
m = (n - 1)/2;
kx = \[Pi]/(d m) Range[0. - m, m];
ky = \[Pi]/(d m) Range[0. - m, m];
kz = \[Pi]/(d m) Range[0. - m, m];
F = Outer[f, x, y, z];
fF = 1/(2 \[Pi])^3 d^3 RotateRight[#, {m, m, m}] &@
Fourier[#, FourierParameters -> {1, -1}] &@
RotateLeft[F, {m, m, m}] // Chop;

data = Flatten[
Table[{kx[[i1]], ky[[i2]], kz[[i3]], fF[[i1, i2, i3]]}, {i1,
1, n}, {i2, 1, n}, {i3, 1, n}], 2];

DeleteCases[data, _?(Abs[#[[-1]]] < 10^-2 &)]

(*{{0., 0., 1.01342, 0.504019}, {0., 1.01342, 0., 0.25201}, {1.01342,
0., 0., 1.00804}}*)
$$$$
`