FFTW 末尾补零

Question

你能帮我找出为什么 FFTW 的计划之一在输出数组末尾给出零吗？当我用 Matlab 检查时，“fftw_plan_dft_1d”产生了正确的结果。实数到复数计划“fftw_plan_dft_r2c_1d”在最后做一些零。我不明白为什么。

这是使用这两个计划的简单测试代码。

#include <iostream>
#include <complex.h>
#include <fftw3.h>

using namespace std;

int main()
{
    fftw_complex *in, *out, *out2;
    double array[] = {1.0,2.0,3.0,4.0,5.0,6.0,0.0,0.0};
    fftw_plan p, p2;

    int N = 8;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    for (int i = 0; i < N; i++) {
        in[i] = i+1+0*I;
    }

    in[6] = 0+0*I;
    in[7] = 0+0*I;

    cout << "complex array" << endl;
    for (int i = 0; i < N; i++) {
        cout << "[" << i << "]: " << creal(in[i]) << " + " << cimag(in[i]) << "i" << endl;
    }
    cout << endl;

    cout << "real array" << endl;
    for (int i = 0; i < N; i++) {
        cout << "[" << i << "]: " << array[i] << endl;
    }
    cout << endl;

    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    p2 = fftw_plan_dft_r2c_1d(N, array, out2, FFTW_ESTIMATE);

    fftw_execute(p); /* repeat as needed */
    fftw_execute(p2);

    cout << "fftw_plan_dft_1d:" << endl;
    for (int i = 0; i < N; i++) {
        cout << "[" << i << "]: " << creal(out[i]) << " + " << cimag(out[i]) << "i" << endl;
    }
    cout << endl;

    cout << "fftw_plan_dft_r2c_1d:" << endl;
    for (int i = 0; i < N; i++) {
        cout << "[" << i << "]: " << creal(out2[i]) << " + " << cimag(out2[i]) << "i" << endl;
    }
    cout << endl;

    fftw_destroy_plan(p);
    fftw_destroy_plan(p2);
    fftw_free(in);
    fftw_free(out);
    fftw_free(out2);

    return 0;
}

结果:

complex array
[0]: 1 + 0i
[1]: 2 + 0i
[2]: 3 + 0i
[3]: 4 + 0i
[4]: 5 + 0i
[5]: 6 + 0i
[6]: 0 + 0i
[7]: 0 + 0i

real array
[0]: 1
[1]: 2
[2]: 3
[3]: 4
[4]: 5
[5]: 6
[6]: 0
[7]: 0

fftw_plan_dft_1d:
[0]: 21 + 0i
[1]: -9.65685 + -3i
[2]: 3 + -4i
[3]: 1.65685 + 3i
[4]: -3 + 0i
[5]: 1.65685 + -3i
[6]: 3 + 4i
[7]: -9.65685 + 3i

fftw_plan_dft_r2c_1d:
[0]: 21 + 0i
[1]: -9.65685 + -3i
[2]: 3 + -4i
[3]: 1.65685 + 3i
[4]: -3 + 0i
[5]: 0 + 0i
[6]: 0 + 0i
[7]: 0 + 0i

正如您所看到的，两个计划之间存在这种奇怪的差异，但结果应该是相同的。

Answer 1

正如您所注意到的，fftw_plan_dft_1d 函数计算复数输入序列 X_n 的标准 FFT Y_k，定义为

其中j=sqrt(-1)，对于所有值k=0,...,N-1(从而生成N 数组中的复杂输出 out), . 您可能会注意到，由于输入恰好是实数，因此输出呈现埃尔米特对称性，即 N=8:

out[4] == conj(out[4]); // the central one (out[4] for N=8) must be real
out[5] == conj(out[3]);
out[6] == conj(out[2]);
out[7] == conj(out[1]);

哪里 conj是通常的复共轭算子。

或者当然，当使用 fftw_plan_dft_1d FFTW 不知道输入恰好是真实的，因此不会利用对称性。

另一方面，fftw_plan_dft_r2c_1d 利用了这种对称性，如 "What FFTW Really Computes" section for "1d real data" of FFTW's documentation 中所示。 (强调我的):

As a result of this symmetry, half of the output Y is redundant (being the complex conjugate of the other half), and so the 1d r2c transforms only output elements 0...n/2 of Y (n/2+1 complex numbers), where the division by 2 is rounded down.

因此，在 N=8 的情况下，只有 N/2+1 == 5 复数值填充在 out2 中，留下其余 3 个未初始化(这些值在调用 fftw_plan_dft_r2c_1d 之前恰好为零，不要依赖将它们设置为 0)。如果需要，这些其他值当然可以通过对称性获得:

for (i = (N/2)+1; i<N; i++) {
  out2[i] = conj(out2[N-i]);
}