How do I access the dynamically allocated 2D array of type fftw_complex in FFTW after a 2D r2c transformation?

Question

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

int main(){
    int N=2000;
    int i,j;
    FILE *filepointer;
    filepointer=fopen("2DDFT_spacetime.plt","w");
    //double in[N][N];
    double *in;
    fftw_complex *out;
    fftw_plan p;
    double fx=13.0;
    double fz=9.0;
    double x[N];
    double xstart=0.0;
    double xend=5.0/fx;
    double z[N];
    double zstart=0.0;
    double zend=5.0/fz;
    double dx=(xend-xstart)/(N-1);
    double dz=(zend-zstart)/(N-1);
    x[0]=xstart;
    z[0]=zstart;
    in = (double*) malloc(sizeof(double) * N * N); //allocates input array
    out = fftw_alloc_complex(N*((int)floor(N/2)+1)); //wrapper function ;allocates output array
    p = fftw_plan_dft_r2c_2d(N,N, in, out, FFTW_MEASURE);
    for(i=1;i<N;i++) {
        x[i]=x[i-1]+dx;
    }
    for(i=1;i<N;i++) {
        z[i]=z[i-1]+dz;
    }
    for(i=0;i<N;i++) {
        for(j=0;j<N;j++) {
            in[i*N+j]=cos(2*M_PI*fx*x[i]+2*M_PI*fz*z[j]);
        }
    }

    fftw_execute(p);

    fprintf(filepointer,"TITLE =\"FFTW\"\nVARIABLES=\"Wavenumber-x\"\n\"Wavenumber-z\"\n\"Amplitude\"\nZONE T=\"Amplitude\"\n I=%d, J=%d, K=1, ZONETYPE=Ordered\n DATAPACKING=POINT\n DT=(SINGLE SINGLE SINGLE)\n",N,(int)floor(N/2)+1);
    for(j=0;j<(int)floor(N/2)+1;j++) { 
        for(i=0;i<N;i++) { 
            fprintf(filepointer," %.9E %.9E %.9E\n",i/(xend-xstart),j/(zend-zstart),sqrt(pow(out[i*((int)floor(N/2)+1)+j][0],2)+pow(out[i*((int)floor(N/2)+1)+j][1],2)));
        }
    }
    fftw_destroy_plan(p);
    free(in); 
    fftw_free(out);
    fclose(filepointer);
    return(1);
}

I begin by allocating memory for a NxN double array and a NxN fftw_complex array, which is defined as typedef double fftw_complex[2] in the FFTW library. I assign real numbers to the array of doubles, do the real to complex FFT, and get the output in the array of fftw_complex. Should I access the real and imaginary parts of the output complex number as out[i*((int)floor(N/2)+1)+j][0] and out[i*((int)floor(N/2)+1)+j][1] respectively?

Hiding arrays behind typedef is a bad idea. But they probably did it for a reason. Perhaps they meant [0] to be real and [1] to be complex? Also, C has built-in support for complex numbers since forever, so the best solution is to drop the icky custom types in favour for standard C. — Lundin
As for how to use 2D arrays. - Change double *in; to double (*in)[N];. - Change malloc to in = malloc(sizeof(double) * N * N); (remove the cast) - Access array by in[i][j]. — Lundin
@Lundin Indeed, [0] is meant to be real and [1] is meant to be complex, as mentioned in the documentation and confirmed by me with my experiments with 1D FFTW. Using the above program, I am now almost certain that the output variables must be accessed in the way I have mentioned in my question. — mmrbest
@Lundin As for the complex numbers, if the header for the complex number is included before the fftw header, FFTW will work with the complex type as defined by C. This has been mentioned in the documentation. — mmrbest

francis francis · Accepted Answer · 2019-12-09T21:19:45

The most beautiful course of action is to include the standard native complex type #include <complex.h> before <fftw3.h> as signaled in the documentation of FFTW.

In particular, if you #include before , then fftw_complex is defined to be the native complex type and you can manipulate it with ordinary arithmetic (e.g. x = y * (3+4*I), where x and y are fftw_complex and I is the standard symbol for the imaginary unit);

The following example is compiled using gcc main.c -o main -Wall -lfftw3 -lm

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

int main(void){

    fftw_plan p;
    unsigned long int N = 10;
    fftw_complex *in=fftw_malloc(N*sizeof(fftw_complex));
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftw_complex *out=fftw_malloc(N*sizeof(fftw_complex));
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}

    printf("sizeof fftw complex %ld\n",sizeof(fftw_complex));
    p=fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    if (p==NULL){fprintf(stderr,"plan creation failed\n");exit(1);}

    unsigned int i;
    for(i=0;i<N;i++){
                in[i]=30.+12.*sin(2*3.1415926535*i/((double)N));
    }

    fftw_execute(p);

    for (i = 0; i < N; i++){
        printf("result: %d || %g %gI |\n", i, creal(out[i]), cimag(out[i]));
    }


    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

    return(0);
}

It shows how to retreive the imaginary and real part of a complex. Now any function operating on complex number can be used, as listed in https://en.cppreference.com/w/c/numeric/complex .

How do I access the dynamically allocated 2D array of type fftw_complex in FFTW after a 2D r2c transformation?

1 Answers