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

这对我来说似乎工作正常：

// 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 回答
1

您的代码看起来正确。我在 64 位机器上使用 GCC gcc-4_7-branch 修订版 189773 进行了尝试，它正在返回

10 -2
-4 0

正如它应该。

还尝试了具有printf格式（%f、%g、%lg）的不同输出技术，并creal以良好的度量访问数据——它总是返回预期的结果。

您是否使用已quad_precision启用的 GCC？这是我现在无法测试的事情之一。

于 2012-10-29T13:10:31.490 回答