FFTW plan rountine returns back null

Question

I need someone with experience using FFTW

I am writing a program which needs to do a real to complex transform, but my planning routine returns back null and I am not sure why. I am passing valid integer parameters for the size, and non-NULL pointers for the array. I have consulted the documentation but all it says is that if it cannot make the plan it returns NULL, but besides doing the opposite of what I mentioned above, it does not list any other reasons for my plan failing. Below is piece of sample code

int size ={64, 128, 256};

float *spatial = malloc(size[1]*size[1]*sizeof(float));
fftwf_complex *fourier = fftwf_alloc_complex(size[1]*size[1]);

for(int i=0; i< size[1]; i++){
    for(int j=0; j<size[1]; j++){
        spatial = //fill spatial array;
    }
}

for(int i=0; i< size[1]; i++){
    for(int j=0; j<size[1]; j++){
        fourier = //fill fourier array;
    }
}

fftwf_plan plan = fftwf_plan_dft_r2c_2d(size[1], size[1], spatial, fourier, FFTW_FORWARD);
fftwf_execute(plan);

the fftwf_plan_dft_r2c_2d() continually returns back NULL for plan and I am not sure why

Note: I am using Windows 8.1 with Visual Studios 2012

Answer 1

I was able to reproduce your problem on ubuntu... The problem comes from the flag FFTW_FORWARD used by fftwf_plan_dft_r2c_2d. I don't know where it is documented, but it is not working. The reason is that fftwf_plan_dft_r2c_2d is always a forward transform, the backward transform being fftwf_plan_dft_c2r_2d

Moreover, you can reduce the memory footprint of your program. Indeed, in case of r2c transforms, only half the complex coefficients are computed and stored. Hence, the array allocated by fftwf_alloc_complex(size[1]*size[1]); is about twice the needed size. Take a look at the documentation to know more about the extra padding and the layout of frequencies.

Here is a program compiled by gcc main.c -o main -lfftw3f -lm -std=c99 (may be different on your computer...). The flag FFTW_FORWARD is replaced by FFTW_ESTIMATE. And #include <complex.h> is placed before #include <fftw3.h> so that fftw can use the float complex type of complex.h

#include <complex.h>

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


int main(int argc, char* argv[]){
    int size [3] ={64, 128, 256};

    float *spatial = malloc(size[1]*size[1]*sizeof(float));
    if(spatial==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftwf_complex *fourier = fftwf_alloc_complex(size[1]*(size[1]/2+1));
    if(fourier==NULL){fprintf(stderr,"fftwf_alloc_complex failed\n");exit(1);}
    for(int i=0; i< size[1]; i++){
        for(int j=0; j<size[1]; j++){
            spatial[i*size[1]+j] =1; //fill spatial array;
        }
    }

    for(int i=0; i< size[1]; i++){
        for(int j=0; j<(size[1]/2+1); j++){
            fourier[i*(size[1]/2+1)+j]  =2.0+I ;//fill fourier array;
        }
    }

    fftwf_plan plan = fftwf_plan_dft_r2c_2d(size[1], size[1], spatial, fourier, FFTW_ESTIMATE);

    if(plan==NULL){printf("fftwf_plan_dft_r2c_2d : big problem %d\n",size[1]);fflush(stdout);}
    fftwf_execute(plan);


    printf("average is %g\n",creal(fourier[0])/(size[1]*size[1]));

    fftwf_destroy_plan(plan);
    free(fourier);
    free(spatial);

    return 0;
}