FFTW 最后用零填充

Question

你能帮我找出为什么 FFTW 的一个计划在输出数组的末尾给出零吗？当我用 Matlab 检查时，"fftw_plan_dft_1d" 产生了正确的结果。 Real to Complex plan "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_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 个 unitilialized（那些在调用 fftw_plan_dft_r2c_1d 之前，值恰好为零，不要依赖它们被设置为 0）。如果需要，那些其他值当然可以通过对称获得：

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

FFTW 最后用零填充

FFTW filled with zeros at the end

fft

complex-numbers

fftw