Difficulty with zeros from 2D FFTW

Question

I have a real matrix represented by the structure double input[N][M] where I want to take a 1D FT for each column j = 0..M along the N direction. My attempt at this is as follows:

#include <fftw3.h>
#include <math.h>

#define M_PI 3.14159265358979323846264338327

int main(void)
{
    int N = 16; int M = 2;

    double input[N][M];             // input data
    double *in = (double *) input;

    fftw_complex output[N][M];      // output data
    fftw_complex *out = (fftw_complex *) output;

    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < M; ++j)
        {
            double x = (double) i / (double) N;
            input[i][j] = sin(2*(j+1) * M_PI * x);
            fprintf (stderr, "input[%d][%d] = %.9f, ", i, j, input[i][j]);
        }
        fprintf (stderr, "\n");
    }
    fprintf (stderr, "\n");

    // setup plans
    int rank = 1; int howmany = M;
    int istride = M; int ostride = M; int idist = 1; int odist = 1;

    fftw_plan fp = fftw_plan_many_dft_r2c(rank, &N, howmany,
                                          in, NULL, istride, idist,
                                          out, NULL, ostride, odist,
                                          FFTW_MEASURE);


    // take the forward transform
    fftw_execute(fp); fprintf (stderr, "FFT complete\n");

    for (int j = 0; j < M; ++j)
        for (int i = 0; i < N/2; ++i)
            fprintf (stderr, "OUT[%3d] = (%.4f + %.4fi)\n",
                    i, output[i][j][0], output[i][j][1]);

    fprintf (stderr, "\n");

    fftw_destroy_plan(fp);

    return 0;
}

which I'm compiling with gcc fft.c -std=c99 -g -lfftw3 -lm. However the code doesn't appear to be working: FT output is all zeros (see here)

Documentation for the function is here.

EDIT: update, so it only seems to work for FFTW_ESTIMATE and not any of the other flags. Any idea why this might be?

Answer 1

Ah cracked it! From the Planner Flags page:

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan. The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

So it was setting up the plan and overwriting the input: easy solution is to move plan creation before initialisation of data array.

Working example code (which FTs forward, then back) here.

Answer 2

The culprit is probably this declaration:

double *in = (double *) in;

This declares in as a pointer to a double (or an array of doubles). You then make it point to itself, but it doesn't point anywhere special to begin with. This means that when the pointer is dereferenced it's the same as dereferencing an uninitialized pointer, which is undefined behavior, and will most likely cause a crash.

You probably mean to make it point to input?