0

我有一个由结构表示的实矩阵double input[N][M]，我想沿 N 方向对 j = 0..M 的每一列进行 1D FT。我的尝试如下：

#include <fftw3.h>
#include <math.h>

#define M_PI 3.14159265358979323846264338327

int main(void)
{
    int N = 16; int M = 2;

    double input[N][M];             // input data
    double *in = (double *) input;

    fftw_complex output[N][M];      // output data
    fftw_complex *out = (fftw_complex *) output;

    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < M; ++j)
        {
            double x = (double) i / (double) N;
            input[i][j] = sin(2*(j+1) * M_PI * x);
            fprintf (stderr, "input[%d][%d] = %.9f, ", i, j, input[i][j]);
        }
        fprintf (stderr, "\n");
    }
    fprintf (stderr, "\n");

    // setup plans
    int rank = 1; int howmany = M;
    int istride = M; int ostride = M; int idist = 1; int odist = 1;

    fftw_plan fp = fftw_plan_many_dft_r2c(rank, &N, howmany,
                                          in, NULL, istride, idist,
                                          out, NULL, ostride, odist,
                                          FFTW_MEASURE);


    // take the forward transform
    fftw_execute(fp); fprintf (stderr, "FFT complete\n");

    for (int j = 0; j < M; ++j)
        for (int i = 0; i < N/2; ++i)
            fprintf (stderr, "OUT[%3d] = (%.4f + %.4fi)\n",
                    i, output[i][j][0], output[i][j][1]);

    fprintf (stderr, "\n");

    fftw_destroy_plan(fp);

    return 0;
}

我正在编译的gcc fft.c -std=c99 -g -lfftw3 -lm. 但是代码似乎不起作用： FT 输出全为零（请参见此处）

该函数的文档在此处。

编辑：更新，所以它似乎只适用于 FFTW_ESTIMATE 而不是任何其他标志。知道为什么会这样吗？

4

2 回答 2

1

啊破解了！从Planner Flags 页面：

重要提示：规划器会在规划期间覆盖输入数组，除非已保存的规划（参见 Wisdom）可用于该问题，因此您应该在创建规划后初始化输入数据。唯一的例外是 FFTW_ESTIMATE 和 FFTW_WISDOM_ONLY 标志，如下所述。

所以它正在设置计划并覆盖输入：简单的解决方案是在初始化数据数组之前移动计划创建。

工作示例代码（FT 向前，然后向后）在这里。

于 2013-06-11T20:54:21.410 回答
0

罪魁祸首可能是这个声明：

double *in = (double *) in;

这声明in为指向 a double（或双精度数组）的指针。然后你让它指向它自己，但它并没有指向任何特别的地方。这意味着当指针被取消引用时，它与取消引用未初始化的指针相同，这是未定义的行为，并且很可能会导致崩溃。

你可能的意思是让它指向input?

于 2013-06-11T17:22:03.723 回答