I am writing a Python program which should do some simulation. I have a working code example in C, which uses the FFTW library (version 2.1.5) to do Fourier transforms.
A simplified version of the C code is as follows:
#include <rfftw.h>
#include <stdio.h>
#include <math.h>
int main(int argc, char** argv)
{
size_t i, n = 2;
fftw_real *vx = (fftw_real*) malloc(8 * sizeof(fftw_real));
for (i=0; i<8; ++i)
vx[i] = i%2 ? 1.0f : 0.5f;
rfftwnd_plan plan_rc;
plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE);
rfftwnd_one_real_to_complex(plan_rc,(fftw_real*)vx,(fftw_complex*)vx);
for (i=0; i<8; ++i)
printf("Value: %f\n", vx[i]);
return 0;
}
I am now trying to translate this into equivalent Python code using Python FFTW which is a wrapper around FFTW 3. At the moment I have the following (also simplified):
import fftw3
import numpy as np
vxr = np.zeros(8, np.float64)
vxc = np.zeros(8, np.complex)
for i in range(8):
vxr[i] = 1.0 if i%2 else 0.5
plan_rc = fftw3.Plan(vxr, vxc, direction = 'forward')
plan_rc()
print vxc
Both versions generate output, however this isn't equal and in the real program it also gives very different results. I think I don't really understand what happens during the in-place real-to-complex conversion in the C code and thus cannot reproduce this in Python.