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.