复数数组的 FFT

Question

使用 FFTW，我想执行复数数组的 IFFT，从文件 "numbers.txt" 中读取数据（它们以 a+b*I 格式存储，每行一个复数).我尝试了以下方法，但它不起作用：

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

样本输入是（每行一个数字）：

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

特别是，我在从文件中读取数据时遇到了问题。有人有什么建议吗？提前致谢

编辑：

经过一些改进，我现在使用这个修改后的代码：

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

现在的问题是我总是得到一个实值转换数组。我该如何解决？

Answer 1

这个

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

必须

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

因为 fftw_complex 只是

typedef double fftw_complex[2];

Answer 2

代码没有为输入分配内存（这是主要问题）

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

While 循环应该检查 == 2 而不是 > 0。还添加了空格以允许复数中的空格。如果这不起作用，建议 OP 提供示例输入。
```
// while(fscanf(file, "%f+%fi", &re, &im) > 0) {
while(fscanf(file, "%f +%f i", &re, &im) == 2) { 
```
代码在循环 i == size.

建议更容易编码和维护 malloc() 风格：

// 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);

其他格式可能提供更多信息

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

健壮的代码会检查 scanf("%d", &size); 和 malloc() 的 return 值。
重新格式化代码会更容易理解。

Answer 3

你有几个问题。 fftw_plan_dft_c2r_1d的调用顺序如下：

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);

因此，您遇到了以下问题：

您永远不会分配 fftw_complex *in 数组。应该是

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

不需要in_cpx，因为in很复杂。
建议对 out 也使用 fftw_malloc。
fftw_complex 是双精度，所以 re 和 im 也应该是双精度。
你的printfwill not work for C99 complex numbers
检查 scanf 和 fscanf 的 returns。
Update 正如chux上面说的，当i >= size时跳出循环。另外，检查以确保循环完成时 i == size。

我没有可用的 c99，所以我无法访问新的 built-in complex type，但修复可能看起来像（未测试）：

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

复数数组的 FFT

FFT of an array of complex numbers

c

complex-numbers

fftw