1

FFTW3 ライブラリを使用して C で Matlab fft2() 関数を実装しようとしています。

しかし、私は異なる結果を得ました。

行列 [1 2; 3 4]、Matlab fft2 での結果は [10 -2; -4 0] と FFTW3 [10 -2; (3.05455e-314 + 2.122e-314i) (9.31763e-315 + 9.32558e-315i)]

次のコードを使用しました。

fftw_plan planG;
fftw_complex *inG, *outG;

inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2 * 2);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2 * 2);

inG[0][0]=1.0;
inG[1][0]=2.0;
inG[2][0]=3.0;
inG[3][0]=4.0;

inG[0][1]=0.0;
inG[1][1]=0.0;
inG[2][1]=0.0;
inG[3][1]=0.0;

planG = fftw_plan_dft_2d(2, 2, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(planG);

私は何を間違っていますか？

前もって感謝します。

4

2 に答える 2

1

あなたのコードは正しいようです。64 ビット マシンで GCC gcc-4_7-branch リビジョン 189773 で試したところ、戻ってきました

10 -2
-4 0

あるべきように。

printfまた、フォーマット (%f、%g、%lg) を使用してさまざまな出力手法を試しcreal、適切な方法でデータにアクセスしました。常に期待どおりの結果が返されています。

GCC をquad_precision有効にして使用していますか? それは私が今テストできないことの1つです。

于 2012-10-29T13:10:31.490 に答える
1

それは私にとってはうまくいくようです：

// fftw_test.c

#include <stdio.h>
#include <fftw3.h>

int main(int argc, char * argv[])
{
    fftw_plan planG;
    fftw_complex *inG, *outG;

    inG = fftw_malloc(sizeof(fftw_complex) * 2 * 2);
    outG = fftw_malloc(sizeof(fftw_complex) * 2 * 2);

    inG[0][0] = 1.0; // real input
    inG[1][0] = 2.0;
    inG[2][0] = 3.0;
    inG[3][0] = 4.0;

    inG[0][1] = 0.0; // imag input
    inG[1][1] = 0.0;
    inG[2][1] = 0.0;
    inG[3][1] = 0.0;

    planG = fftw_plan_dft_2d(2, 2, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(planG);

    printf("outg[0] = { %g, %g }\n", outG[0][0], outG[0][1]);
    printf("outg[1] = { %g, %g }\n", outG[1][0], outG[1][1]);
    printf("outg[2] = { %g, %g }\n", outG[2][0], outG[2][1]);
    printf("outg[3] = { %g, %g }\n", outG[3][0], outG[3][1]);

    return 0;
}

コンパイルして実行します。

$ gcc -Wall -lfftw3 fftw_test.c -o fftw_test
$ ./fftw_test 
outg[0] = { 10, 0 }
outg[1] = { -2, 0 }
outg[2] = { -4, 0 }
outg[3] = { 0, 0 }
$

結果を表示するために使用しているコードが間違っているだけなのでしょうか?

3x3の場合、私も一致しています：

// fftw_test_3_x_3.c

#include <stdio.h>
#include <fftw3.h>

#define M 3
#define N 3

int main(int argc, char * argv[])
{
    fftw_plan planG;
    fftw_complex *inG, *outG;
    int i;

    inG = fftw_malloc(sizeof(fftw_complex) * M * N);
    outG = fftw_malloc(sizeof(fftw_complex) * M * N);

    for (i = 0; i < M * N; ++i)
    {
        inG[i][0] = (double)(i + 1);
        inG[i][1] = 0.0; // imag input
    }

    planG = fftw_plan_dft_2d(M, N, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(planG);

    for (i = 0; i < M * N; ++i)
    {
        printf("outg[%d] = { %g, %g }\n", i, outG[i][0], outG[i][1]);
    }

    return 0;
}

コンパイルして実行します。

$ gcc -Wall -lfftw3 fftw_test_3_x_3.c -o fftw_test_3_x_3
$ ./fftw_test_3_x_3 
outg[0] = { 45, 0 }
outg[1] = { -4.5, 2.59808 }
outg[2] = { -4.5, -2.59808 }
outg[3] = { -13.5, 7.79423 }
outg[4] = { 0, 0 }
outg[5] = { 0, 0 }
outg[6] = { -13.5, -7.79423 }
outg[7] = { 0, 0 }
outg[8] = { 0, 0 }
$

Octave (MATLAB クローン) と比較します。

octave:1> a = [1 2 3 ; 4 5 6 ; 7 8 9]
a =

   1   2   3
   4   5   6
   7   8   9

octave:2> b = fft2(a, 3, 3)
b =

   45.00000 +  0.00000i   -4.50000 +  2.59808i   -4.50000 -  2.59808i
  -13.50000 +  7.79423i    0.00000 +  0.00000i    0.00000 +  0.00000i
  -13.50000 -  7.79423i    0.00000 -  0.00000i    0.00000 -  0.00000i

octave:2>
于 2012-10-29T12:50:53.000 に答える