Fourier transform in three dimensions

Question

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}\}\,. $$

Answer 1

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}

Answer 2

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

Answer 3

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}}*)
```