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$?


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


3 Answers 3


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

enter image description here

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

enter image description here

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

enter image description here

  • $\begingroup$ 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. $\endgroup$
    – user12588
    Jun 30, 2022 at 9:31
  • $\begingroup$ Look at my second answer $\endgroup$ Jun 30, 2022 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}}*)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.