Search code examples
ccomplex-numbersfftw

FFT of an array of complex numbers


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?


Solution

  • 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;
    }