FFTW 结果为零答案

【问题标题】：FFTW results in zerosFFTW 结果为零
【发布时间】：2012-03-02 15:43:26
【问题描述】：

我正试图让 FFTW 在 C 中工作。它曾经为另一个项目（在 JNI 中）工作，我或多或少地从那个项目中复制代码，遗憾的是没有结果。

首先我生成一个正弦信号，如下所示：

double* generateSignal() {
    int fs=44100;
    double fsd = 44100.0; // fs in double format
    double f1=1000.0;
    int i;
    double PI = 3.141592653589793238462643;

    double t[fs];
    double value = 0.0;
    for (i = 0; i < fs; i++) {
        t[i] = value;
        value += 1.0/fsd;
    }

    double* signal = (double*) malloc(sizeof(double) * fs);
    for (i = 0; i < fs; i++) {
        signal[i] = sqrt(2) * sin(2 * PI * f1 * t[i]);
    }

    return signal;
}

这行得通，我只是为了完整性而发布它。

接下来，我想使用 FFTW 转换信号，我使用以下方法（基于FFTW documentation）：

void processSignal(double* signal) {
    int size = 44100;
    int i;

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

    for (i = 0; i < size; i++) {
        double* ptr = in[i];
        *ptr = signal[i];       // set first double, real part
        *(ptr + 1) = 0.0;       // set second double, imaginary part
    }

    fftw_plan p = fftw_plan_dft_1d(size, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 
    fftw_execute(p);

    for (i = 0; i < size; i++) {
        double* ptr = out[i];
        signal[i] = *ptr;       // get real part
    }

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

请从 FFTW 文档中注意这一点：typedef double fftw_complex[2];

现在，这导致signal-array 的所有值都为 -0.000000。我真的看不出这段代码的问题，所以请任何人指出我做错了什么？

谢谢。

PS：为了清楚起见，我从我的代码中删除了打印语句。

【问题讨论】：

您是否收到任何警告/错误？此行不应该编译：double* ptr = in[i];
没有，即使设置了 -Wall 标志（我意识到，这很奇怪）。我明白你为什么它不能编译的观点，但如果是这种情况，我会通过转换为双指针来解决它。无论如何，我想这不是导致问题的原因。
你没有转换 fftw_malloc(sizeof(fftw_complex) * size); 的输出；
您的预期输出是什么？你知道它应该是什么吗？编译带有警告并查看它的内容。我看不出有什么问题，但是如果您将 ptr 用作 ptr[0] 和 ptr[1] 会更清楚。
@Mystical 这会编译，因为 in 本质上是一个二维数组。所以 in[i] 是一个指针类型。

标签： c signal-processing fftw

【解决方案1】：

我认为问题可能在于您正在用当前迭代的实部覆盖上一次迭代的虚部。

for (i = 0; i < size; i++) {
        double* ptr = in[i];    //  <-- here's a problem
        *ptr = signal[i];       // set first double, real part
        *(ptr + 1) = 0.0;       // set second double, imaginary part
    }

i 在每次迭代中递增，所以在第一次迭代中，ptr 指向输入复数[0] 的实部，ptr + 1 指向复数[0] 的虚部]，但在第二次迭代中，ptr 指向复数[0] 的虚部，ptr + 1 指向复数[1] 的实部。

解决此问题的一些建议可能是：

for (i = 0, j = 0; i < size; i++, j+= 2) {
    double* ptr = in[j];    // j increments by 2 making ptr alays point to real part 
    *ptr = signal[i];       // set first double, real part
    *(ptr + 1) = 0.0;       // set second double, imaginary part
}

或

double* ptr = in[0]
for (i = 0; i < size; i++) {
        *ptr++ = signal[i]; // set first double, real part
        *ptr++ = 0.0;       // set second double, imaginary part
    }

【讨论】：

【解决方案2】：

天哪，这太愚蠢了。

这些值的大小类似于 -2.0E-9。结果四舍五入，因此打印 -0.00000。

非常感谢弗雷德。他的问题“你的预期输出是什么？”让我仔细查看了我之前提到的 JNI 项目的输出，并意识到我的错误。

也感谢其他答案！

【讨论】：