2D FFT (Fast Fourier Transform librerie)

[foton7]

Does anyone know a good free library to do Fourier Transforms (FFT or DFT). I know FFTW but I'm having some problems with it. I want an alternative that do FFT in two dimensions with complex numbers. The libraries I have found doesn't fulfill this requirements.

Thank you

 

[eric_in_chica]

have you tried sourceforge.net?

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

both software packages will spit out C code

 

[foton7]

Thank you Eric,

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.

 

[Lt_Dax]

What are the problems you are having with FFTW? Might I be able to help?

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.

 

[foton7]

This is the piece of code that it does not work:

...
#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;
}

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?

 

[Lt_Dax]

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:

#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)?

 

[foton7]

Many thanks!

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.

 

[Lt_Dax]

Yes, it's very common to solve a problem in programming without ever figuring out what the problem was! It's happened to me many, many times...