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 フラグ ページから:

重要: 保存された計画 (Wisdom を参照) がその問題に対して利用可能でない限り、プランナは計画中に入力配列を上書きするため、計画の作成後に入力データを初期化する必要があります。これに対する唯一の例外は、以下で説明するように、FFTW_ESTIMATE および FFTW_WISDOM_ONLY フラグです。

したがって、計画を設定して入力を上書きしていました。簡単な解決策は、データ配列の初期化の前に計画の作成を移動することです。

動作するサンプル コード (FT が前方に移動してから戻る)はこちら.

于 2013-06-11T20:54:21.410 に答える