FFTW filled with zeros at the end

Question

Can you help me find out why one of the FFTW's plans gives zeroes at the end of an output array? The "fftw_plan_dft_1d" yields proper result as I checked it with Matlab. The Real to Complex plan "fftw_plan_dft_r2c_1d" makes some zeroes at the end. I don't understand why.

Here is the simple testing code using both plans.

#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;
}

Result:

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

As you can see there is this strange difference between both plans and the result should be the same.

Answer 1

As you have noted, the fftw_plan_dft_1d function computes the standard FFT Y_k of the complex input sequence X_n defined as

where j=sqrt(-1), for all values k=0,...,N-1 (thus generating N complex outputs in the array out), . You may notice that since the input happens to be real, the output exhibits Hermitian symmetry, that is for 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]);

where conj is the usual complex conjugate operator.

Or course, when using fftw_plan_dft_1d FFTW doesn't know the input just happens to be real, and thus does not take advantage of the symmetry.

The fftw_plan_dft_r2c_1d on the other hand takes advantage of that symmetry, and as indicated in "What FFTW Really Computes" section for "1d real data" of FFTW's documentation (emphasis mine):

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.

Thus in your case with N=8, only N/2+1 == 5 complex values are filled in out2, leaving the remaining 3 unitilialized (those values just happened to be zeros before the call to fftw_plan_dft_r2c_1d, do not rely on them being set to 0). If needed, those other values could of course be obtained from symmetry with:

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