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（从而在数组out 中生成N 复杂输出），。 您可能会注意到，由于输入恰好是真实的，因此输出表现出厄米对称性，即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_1dFFTW 时不知道输入恰好是真实的，因此没有利用对称性。

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

由于这种对称性，输出 Y 的一半是冗余的（是另一半的复共轭），因此 1d r2c 仅变换 Y 的输出元素 0...n/2 ( n/2+1 复数），其中除以 2 向下舍入。

因此，在您使用N=8 的情况下，只有N/2+1 == 5 复杂值被填充到out2，剩下的3 个单元化（在调用@ 之前，这些值恰好为零） 987654337@，不要依赖它们被设置为 0）。如果需要，这些其他值当然可以通过以下对称性获得：

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