FFTW: size of output array for 2D Fourier Transform of 3D data

Question

I have a real 3D array of dimensions Nx*Ny*Nz and want to take a 2D Fourier transform for each z value using FFTW. Here the z index is the fastest varying in memory. Currently the following code works as expected:

int Nx = 16; int Ny = 8; int Nz = 3;

// allocate memory
const int dims = Nx * Ny * Nz;

// input data (pre Fourier transform)
double *input = fftw_alloc_real(dims);              

// why is this the required output size?
const int outdims = Nx * (Ny/2 + 1) * Nz;           
// we want to perform the transform out of place
// (so seperate array for output)
fftw_complex *output = fftw_alloc_complex(outdims);    


// setup "plans" for forward and backward transforms
const int rank = 2; const int howmany = Nz;
const int istride = Nz; const int ostride = Nz;
const int idist = 1; const int odist = 1;
int n[] = {Nx, Ny};
int *inembed = NULL, *onembed = NULL;


fftw_plan fp = fftw_plan_many_dft_r2c(rank, n, howmany,
        input, inembed, istride, idist,
        output, onembed, ostride, odist,
        FFTW_PATIENT);

fftw_plan bp = fftw_plan_many_dft_c2r(rank, n, howmany,
        output, onembed, ostride, odist,
        input, inembed, istride, idist,
        FFTW_PATIENT);

As I understand it, transforming a 1D sequence of length N requires (N/2 + 1) complex values so why does the above code crash if instead I set outdims = (Nx/2 + 1)*(Ny/2 + 1)*Nz as one might expect for a 2D transform?

Secondly am I right in thinking that I can access the real and imaginary parts of modes from qx = 0 to Nx/2 (inclusive) using the following:

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][1] )

EDIT: Full code and Makefile for those who want to play around. Assumes fftw and gsl installed.

EDIT2: If I understand correctly, the indexing (allowing positive and negative frequencies) should be (probably getting too messy for a macro!):

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][1] )

for (int qx = -Nx/2; qx < Nx/2; ++qx)
    for (int qy = 0; qy <= Ny/2; ++qy)
        outputRe(qx, qy, d) = ...

where outputRe(-Nx/2, qy, d) points to the same data as outputRe(Nx/2, qy, d). In practice I would probably just to loop over the first index and convert to a frequency, rather than the other way round!

Answer 1

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.