1
$\begingroup$
#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 am using the above program to initialize a 2D field of 2000x2000 elements of data (say in a x-z plane) and am using FFTW to transform the data from the real plane to the wavenumber plane, which, as per my best knowledge, has 2000x1001 elements as per FFTW documentation. The data is given by the function:

$$ in=\cos(2\pi f_x+2\pi f_z) $$

where $f_x$ is the frequency in x-direction, $f_z$ is the frequency in z-direction. I expect the maximum amplitude at my chosen wavenumbers to be 1.0. However, I get the number 1999959 as the maximum amplitude. I know that I have to normalize this to get the maximum amplitude, but how do I calculate the number that I have to divide the maximum amplitude with to arrive at Amplitude=1.0? If the number is 2000x1001, does this mean that I cannot obtain the amplitude exactly as 1.0 due to some numerical error?

$\endgroup$
2
  • $\begingroup$ Do you devide by array length one time when doing forward and backward transformation, or do you devide by the root each time? Also make sure that the FT definition you have in your head matches: fftw.org/fftw3_doc/… $\endgroup$ – MPIchael Nov 29 '19 at 14:28
  • $\begingroup$ I am doing only the forward transform as seen in the code. Also, I divide the amplitude of the Fourier transform by the number of elements in the Fourier plane and not sqrt(number of elements in the Fourier plane). $\endgroup$ – mmrbest Dec 2 '19 at 6:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.