Fourier transform using FFTW3 library - absurd numbers result

May 23, 2010 at 8:34pm
I have written a small program which uses the most popular Fast Fourier Transform library FFTW, which has a helpful manual on its website.

My program takes a 2-dimensional plot as input (a graph where value of the function varies with the two axes x and T, a 2-d array), and the FFTW execute function should produce another 2-d function.

My input function has the appearance of exponentially falling towards the middle, or rising towards the edges. Its fourier transform should have a large 'bump' in the middle, but worse than this, my output numbers are of the order 10^-300.

Since FFTW requires some trickery to make sure the 2-d array is in 1-d format, C-major order, I assume it is something to do with that. I have carefully reveiwed my code, and I cannot see what the error is. Can anybody help? Here is the code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#include <iostream>
#include <cstdlib>
#include <fftw3.h>
#include <math.h>
#include <fstream>
using namespace std;

int main()
{
    const int Lx = 24;
    const int Lt = 48;
    int var_x;
    int var_T;

    double F[Lx][Lt];
    double pi = 4*atan(1);

    fftw_complex *in, *out;
    fftw_plan p;

    // Declare one-dimensional contiguous arrays of dimension   Lx*Lt
    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Lx*Lt);
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Lx*Lt);
    

    // create plan, an object containing necessary parameters
    p = fftw_plan_dft_2d(Lx, Lt, in, out, FFTW_FORWARD, FFTW_MEASURE);
    
    ifstream input;
    input.open("out_b5p40kp13625-24x48_pion_80-200_2d");
    // Put input data into F[x][T]
    // Careful - the loop order depends on how the data is recorded
    for(int x = 0; x < Lx; x++)
    {
	for(int T = 0; T < Lt; T++)
	{
	    input >> var_T;
	    input >> var_x;
	    input >> F[x][T];
	}
    }

    // Now turn the array F[x][T] into C-order (row major) 1-d array
    // It also happens to be an array of fftw_complex values
    // However, we only need to store each value in the real part
    // The order here is important. Rightmost index loops on the inside
    for(int x = 0; x < Lx; x++)
    {
	for(int T = 0; T < Lt; T++)
	{
	    in[T + Lt*x][0] = F[x][T];
	    // in[T + Lt*x][1] = F[x][T][1]  : Can enter imaginary part if needed
	}
    }
    input.close();

    // Perform FFT
    fftw_execute(p);

    // Get output data into 2-d form - use same array as for input
    for(int x = 0; x < Lx; x++)
    {
	for(int T = 0; T < Lt; T++)
	{
	    F[x][T] = out[T + Lt*x][0];
	    // F[x][T][1] = out[T + Lt*x][1] : Can get imaginary data too if it exists
	}
    }

    // Print output to file
    // Remember that the array indices x and y represent a different variable
    // e.g. momentum space, k = 2*pix/L
    ofstream output;
    output.open("2d_transform_out");
    for(int x = 0; x < Lx; x++)
    {
	for(int T = 0; T < Lt; T++)
	{
	    output << 2*pi*x/Lx << "\t" << 2*pi*T/Lt << "\t" << F[x][T] << "\n";
	}
	output << "\n";
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

    output.close();
}


The code only compiles on a system which has the FFTW libraries installed. On the Linux system I use, this means compiling with "g++ -o FFTW_2d -I<location of include> -L<location of library> -lfftw3 FFTW2d.cc".

I realize this is highly specific problem, but help for scientific or numeric programming seems hard to come by.
Last edited on May 23, 2010 at 8:41pm
Topic archived. No new replies allowed.