- #1

- 5

- 0

Thank you

You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- Thread starter foton7
- Start date

- #1

- 5

- 0

Thank you

- #2

- 8

- 1

how about input a mathematical representation of your process into mathcad or mathematica?

both software packages will spit out C code

- #3

- 5

- 0

I will have a look to sourceforge.

I need to do it in C or C++. The Fourier transform is just a step in a much bigger software we are developing. We already have a prototype in Matlab and I do the fourier transforms there without any problem.

- #4

- 77

- 0

FFTW3 is very convenient for multidimensional transforms, and an FFT deals with complex numbers by definition (each discrete value is a sum of complex numbers in polar form, you can separate the real and imaginary parts with Euler's formula).

I encourage you to try this library, it is very efficient and is written in C++. The most complicated thing it usually needs is for you to change your multidimensional array into a contiguous 1-D array of complex numbers.

- #5

- 5

- 0

This is the piece of code that __it does not work__:

I have removed the arguments of the function just in case the problem was related with that and so that the function is self contained to post it here.

If I place a “cout” trace just before and just after the first call to fftw_plan_dft_2d then the cout trace after fftw_plan_dft_2d is never printed so I conclude that the program crashes inside that function.

Any idea?

Code:

```
...
#include <complex>//Includes C++ complex numbers (uses templates)
#include <fftw3.h> // Library to perform FFT transform (GPL license)
const int BLOCKSIZE= 10;
...
// Applies an MTF to an image. This filters the image. It is done in the frequency domain
// image_in[BLOCKSIZE][BLOCKSIZE] : Input image
// double MTF[BLOCKSIZE][BLOCKSIZE] : Modulation Transfer function (MTF)
// double image_out[BLOCKSIZE][BLOCKSIZE]: Output image after applying the MTF
//void applyMTF(const double image_in[][BLOCKSIZE], const double MTF[][BLOCKSIZE], double image_out[][BLOCKSIZE]){
void applyMTF(){
fftw_complex *image_v_spat; //Image vector in spatial domain
fftw_complex *image_v_freq; //Image vector in frequency domain
fftw_complex *MTF_v; //Modulation Transfer Function (in frequency domain by definition)
fftw_plan p; //Plan used by the fftw library
int i= 0, j= 0, counter= 0; //Loop counters
//Input parameters temporarily defined inside the function
double image_in[BLOCKSIZE][BLOCKSIZE];
double image_out[BLOCKSIZE][BLOCKSIZE];
double MTF[BLOCKSIZE][BLOCKSIZE];
//Fill image and MTF with dummy values:
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
image_in[i][j]= (double)(i*i - i + j);
MTF[i][j]= (double)(j);
}
}
//Allocate memory
image_v_spat = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
image_v_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
MTF_v = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
//Expand the image and MTF arrays (2-dimensional real) into vectors (1-dimensional complex array). Fill the imaginary part with zeros.
counter= 0;
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
counter++;
//Image
image_v_spat[counter][0]= image_in[i][j]; //Real part
image_v_spat[counter][1]= 0; //Imaginary part
//Same for MTF
MTF_v[counter][0]= MTF[i][j]; //Real part
MTF_v[counter][1]= 0; //Imaginary part
}
}
//Create the plan for the DFT
p= fftw_plan_dft_2d(BLOCKSIZE, BLOCKSIZE, image_v_spat, image_v_freq, FFTW_FORWARD, FFTW_ESTIMATE);
//Execute the plan (perform the Fourier transform of the image)
fftw_execute(p);
//Multiply the image in frequency domain times the MTF taking into account that they are complex numbers: (a0 + a1 i) * (b0 + b1 i) = (a0*b0 - a1*b1) + (a0*b1 + a1*b0)i
//The real part is stored in index 0 and the complex part in index 1
for(counter=0; counter < BLOCKSIZE*BLOCKSIZE; counter++){
image_v_freq[counter][0]= image_v_freq[counter][0] * MTF[counter][0] - image_v_freq[counter][1] * MTF[counter][1]; //Real part
image_v_freq[counter][1]= image_v_freq[counter][0] * MTF[counter][1] + image_v_freq[counter][1] * MTF[counter][0]; //Imaginary part
}
//Do the inverse Fourier transform
//First create the plan for the DFT
p= fftw_plan_dft_2d(BLOCKSIZE, BLOCKSIZE, image_v_freq, image_v_spat, FFTW_BACKWARD, FFTW_ESTIMATE);
//Then execute the plan (perform the inverse Fourier transform of the image)
fftw_execute(p);
//Do the absolute value and store the resulting image in a real 2D matrix instead of a complex 1D vector. This image is then returned by the function by ref.
counter= 0;
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
counter++;
image_out[i][j]= sqrt(image_v_spat[counter][0]*image_v_spat[counter][0] + image_v_spat[counter][1]*image_v_spat[counter][1]);
}
}
//Deallocate memory
fftw_destroy_plan(p);
fftw_free(image_v_spat);
fftw_free(image_v_freq);
fftw_free(MTF_v);
return;
}
```

If I place a “cout” trace just before and just after the first call to fftw_plan_dft_2d then the cout trace after fftw_plan_dft_2d is never printed so I conclude that the program crashes inside that function.

Any idea?

Last edited:

- #6

- 77

- 0

Your code looks good, as far as I can see. I will make one or two suggestions, and then I will post some code which definitely works for getting the transform of a 2-D function.

My first suggestion is to try making the plan before filling the 1-D contiguous arrays image_v_spat and MTF_v with values. Some people report strange things happening when the plan is executed afterwards.

My second suggestion is to use the row-major formula to convert your 2-D arrays into the needed 1-D arrays, rather than using a counter. Although your method will work, there is less chance of an error this way and in my opinion it is cleaner (one less variable to declare).

Here is my working code. Note that I use the same 2-D array to store my output values, since I discard my input:

This can be compiled with g++ -o FFT FFT.cc -I/<include path> -L/<library path> -lfftw3 (of course the details depend on your system or installation).

EDIT: What information is displayed when your program crashes (e.g. segmentation fault)?

My first suggestion is to try making the plan before filling the 1-D contiguous arrays image_v_spat and MTF_v with values. Some people report strange things happening when the plan is executed afterwards.

My second suggestion is to use the row-major formula to convert your 2-D arrays into the needed 1-D arrays, rather than using a counter. Although your method will work, there is less chance of an error this way and in my opinion it is cleaner (one less variable to declare).

Here is my working code. Note that I use the same 2-D array to store my output values, since I discard my input:

Code:

```
#include <iostream>
#include <cstdlib>
#include <fftw3.h>
#include <fstream>
using namespace std;
int main()
{
const int Lx = 4; // Array dimensions
const int Lt = 2;
int var_x;
int var_T;
double F[Lx][Lt][2];
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("2d_function");
// 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][0];
//input >> Fi[x][T][1]; Can get imaginary part if it exists
}
}
// 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][0];
// 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][0] = 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
ofstream output;
output.open("2d_transform_out");
for(int x = 0; x < Lx; x++)
{
for(int T = 0; T < Lt; T++)
{
output << x << "\t" << T << "\t" << F[x][T][0] << "\t" << "+ i " << F[x][T][1] << "\n";
}
output << "\n";
}
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
output.close();
return 0;
}
```

This can be compiled with g++ -o FFT FFT.cc -I/<include path> -L/<library path> -lfftw3 (of course the details depend on your system or installation).

EDIT: What information is displayed when your program crashes (e.g. segmentation fault)?

Last edited:

- #7

- 5

- 0

I haven't been able to figure out what is causing the problem in my code but based on yours I've been able to do functions that suit my needs.

I've encapsulated the direct (forward) transform in one function and the inverse (backwords) transform in other. In this way I do the fourier transform by calling the forward function, then I do the multiplication by the matrix (MTF in my code) and then call the inverse fourier transform function.

I remember using the FFW once before several years ago and I didn't have any problem back then. I still don't know what was the problem this time.

- #8

- 77

- 0

Share: