tags:

views:

45

answers:

0

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.