To help clarify (focusing in 2D as it easily extends to 2D transform of 3D data):
storage
An Nx * Ny
array requires Nx * (Ny / 2 + 1)
complex elements after a Fourier Transform.
First, in the y-direction, the negative frequencies can be reconstructed from the complex conjugate symmetry (that comes from transforming a real sequence). The y modes ky
then run from 0 to Ny/2
inclusive. So for y we need Ny/2 + 1
complex values.
Next we transform in the x-direction where we cannot use the same symmetry assumption, as we are acting on complex-valued y values. Therefore we must include positive and negative frequencies, so x-modes kx
run from -Nx/2 to Nx/2
inclusive. However kx = -Nx/2
and kx = Nx/2
are equivalent so only one is stored (see here). So for x we need Nx
complex values.
access
As tir38 points out the x index post-transform runs from 0 to Nx-1, however this doesn't mean that modes kx
run from 0 to Nx-1. FFTW packs positive frequencies in the first half of the array, then negative frequencies in the second half (in reverse order), like:
kx = 0, 1, 2, ..., Nx/2, -Nx/2 + 1, ..., -2, -1
There are two ways we can think about accessing these elements. First as tir38 suggests we can loop through in order and work out the mode kx
from the index:
for (int i = 0; i < Nx; i++)
{
// produces the list of kxs above
int kx = (i <= Nx/2) ? i : i - Nx;
// here we index with i, but with the knowledge that the mode is kx
outputRe(i, ...) = some function of kx
}
or we can loop through the modes kx
and convert to an index:
for (int kx = -Nx/2; kx < Nx/2; kx++)
{
// work out index from mode kx
int i = (kx >= 0) ? i : Nx + i;
// here we index with i, but with the knowledge that the mode is kx
outputRe(i, ...) = some function of kx
}
The two types of indexing along with the rest of the code can found here.