FFT of an array of complex numbers

Question

Using FFTW, I want to perform the IFFT of an array of complex numbers, reading data from file "numbers.txt" (where they are stored in the format a+b*I, one complex number per each row). I have tried with the following, but it doesn't work:

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

int main(void)
{
fftw_complex *in;
fftw_complex *out;
float re,im;
int size;

printf("Insert size");
scanf("%d", &size);

FILE *file = fopen("numbers.txt", "r");

int i=0;

while(fscanf(file, "%f+%fi", &re, &im) > 0) {
in[i]=re+im*I;
i++;
}
fclose(file);

 fftw_complex *out_cpx;

 fftw_plan ifft;
 out_cpx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*size);
 out = (double *) malloc(size*sizeof(double));

 ifft = fftw_plan_dft_c2r_1d(size, out_cpx, out, FFTW_ESTIMATE);   //Setup fftw plan 

 fftw_execute(ifft);

 printf("Input:    \tOutput:\n");
 for(i=0;i<size;i++)
 {
 printf("%f\t%f\n",(in[i]),out[i]);
 }

 fftw_destroy_plan(ifft);
 fftw_free(out_cpx);
 free(out);


 return 0;
 }

Sample input is (one number per each row):

0+0*I
0.01198+-0.0825632*I
1.32139e-05+1.94777e-05*I 
-1.18289e-11+7.30045e-12*I 
2.71773e-20+0*I 
-1.18289e-11+-7.30045e-12*I 
1.32139e-05+-1.94777e-05*I 
0.01198+0.0825632*I

In particular, I have problems in reading the data from the file. Does someone have any suggestion? Thanks in advance

Edit:

After some improvements, I'm now using this modified code:

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


 int main(void)
{
fftw_complex *in;
fftw_complex *out;
double re,im;
int size;
int i=0;
FILE *file;
fftw_plan ifft;

printf("Insert size");
if (scanf("%d", &size) != 1 || size < 1)
    return 1;

in = fftw_malloc(sizeof(*in)*size);
out = fftw_malloc(sizeof(*out)*size);

file = fopen("numbers.txt", "r");

for (i = 0; i < size && fscanf(file, "%lf+%lf*I\n", &re, &im) == 2; i++)
{
    in[i]= re+im*I;
}

fclose(file);
// Error if i != size?


ifft = fftw_plan_dft_1d(size, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(ifft);

printf("Input:    \tOutput:\n");
for(i=0; i<size; i++)
{
printf("%lf+%lfI\t%lf+%lfI\n", creal(in[i]), cimag(in[i]), creal(out[i]), cimag(out[i]));
}

fftw_destroy_plan(ifft);
fftw_free(in);
fftw_free(out);
//free(out);

return 0;
}

The problem now is that I always get a real-valued transformed array. How can I fix it?

Answer 1

Your have several problems. The calling sequence for fftw_plan_dft_c2r_1d is as follows:

The inverse transforms, taking complex input (storing the non-redundant half of a logically Hermitian array) to real output, are given by:

 fftw_plan fftw_plan_dft_c2r_1d(int n0,
                                fftw_complex *in, double *out,
                                unsigned flags);

Thus, you have the following problems:

1. You never allocate the fftw_complex *in array. It should be

   in = fftw_malloc(sizeof(*in)*size);
   // And later
   fftw_free(in);
   

2. There is no need for in_cpx, since in is complex.

3. Suggest also using fftw_malloc for out.

4. fftw_complex is double precision so re and im should be as well.

5. Your printf will not work for C99 complex numbers

6. Check the returns from scanf and fscanf.

7. Update as chux says above, break out of the loop when i >= size. Also, check to make sure that i == size when the loop is complete.

I don't have c99 available, so I don't have access to the new built-in complex type, but a fix might look like (not tested):

int main(void)
{
    fftw_complex *in;
    double *out;
    double re,im;
    int size;
    int i=0;
    FILE *file;
    fftw_plan ifft;

    printf("Insert size");
    if (scanf("%d", &size) != 1 || size < 1)
        return 1;

    in = fftw_malloc(sizeof(*in)*size);
    out = fftw_malloc(sizeof(*out)*size);

    file = fopen("numbers.txt", "r");

    for (i = 0; i < size && fscanf(file, "%lf+%lf*I\n", &re, &im) == 2; i++)
    {
        in[i]= re + im*I;
    }

    fclose(file);
    // Error if i != size?

    ifft = fftw_plan_dft_c2r_1d(size, in, out, FFTW_ESTIMATE);   //Setup fftw plan 

    fftw_execute(ifft);

    printf("Input:    \tOutput:\n");
    for(i=0; i<size; i++)
    {
        printf("%lf+i%lf\t%lf\n", creal(in[i]), cimag(in[i]), out[i]);
    }

    fftw_destroy_plan(ifft);
    fftw_free(in);
    fftw_free(out);

    return 0;
}

Answer 2

1. Code has no memory allocated for input (this is the major issue)

   fftw_complex *in;
   ...
   in = malloc(size * sizeof *in); 
   // or
   in = calloc(size, sizeof *in);
   ...
   in[i]=re+im*I;
   ...
   free(in);
   

2. While loop should check for == 2 not > 0. Also added spaces to allow spacing in complex number. Suggest OP provide sample input if this does not work.

   // while(fscanf(file, "%f+%fi", &re, &im) > 0) {
   while(fscanf(file, "%f +%f i", &re, &im) == 2) { 
   

3. Code does not check that after loop i == size.

4. Suggest easier to code and maintain malloc() style:

   // out_cpx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*size);
   out_cpx = fftw_malloc(size * sizeof *out_cpx);
   
   // out = (double *) malloc(size*sizeof(double));
   out = malloc(size * sizeof *out);
   

5. Alternate format may be more informative

   // printf("%f\t%f\n",(in[i]),out[i]);
   printf("%e\t%e\n",(in[i]),out[i]);
   

6. Robust code would check the return value from scanf("%d", &size); and malloc().

7. Re-formatting code would make it easier to understand.

Answer 3

this

printf("%f\t%f\n",(in[i]),out[i]);

has to be

printf("%f + %f j\t%f + %f j\n", in[i][0], in[i][1], out[i][0], out[i][1]);

since fftw_complex is just

typedef double fftw_complex[2];