Forums

2D FFTW3 vs. MATLAB difference only near center place

Started by hyun.ha February 17, 2012
Dear friends,

Here I have a problem to re-produce same results with MATLAB and FFTW
during 2 dimensional FFT.

I'm trying to do FFT with Gaussian convolution which coordinates varies
from negative to positive, FFTW3 seems not to calculate the values
correctly. It happens only near center place.
And it affects IFFT and it's shifted results seriously in my calculation.

------------------------------------------------------------------------------------------------------------------
The different results is

from MATLAB FFT,

                        1                                                  
                                      2                                    
                                      3                  
49 :  6.51969993785790e-08 + 1.03627663030827e-06i       
-3.25984630124937e-08 - 1.03729991235754e-06i	        6.61744490042422e-22
+ 1.03627197122179e-06i
50 : -1.63152920904142e-08 - 5.19161012253523e-07i         
-4.35097002202893e-22 + 5.19160428141898e-07i      1.62830243797501e-08 -
5.18134237051512e-07i
51 :  1.69406589448395e-21 + 0.00000000000000i             
-1.69239446844913e-21 - 5.31856388157782e-23i	        1.68738678669962e-21
+ 1.06161378213412e-22i
52 : -1.63152920904142e-08 + 5.19161012253523e-07i        
3.25983530020901e-08 - 5.18135983627880e-07i	       -4.87848114810276e-08 +
5.16089397930887e-07i
53 :  6.51969993785790e-08 - 1.03627663030827e-06i        
-9.76667377967905e-08 + 1.03320616350072e-06i       1.29879317005283e-07 -
1.02810065720927e-06i

from FFTW3,

49 :  6.5197e-08 + 1.0363e-06i                                          
-3.2598e-08 + -1.0373e-06i                                      
-4.2352e-22 + 1.0363e-06i 
50 : -1.6315e-08 + - 5.1916e-07i                                         
-4.2352e-22 +  5.1916e-07i                                       1.6283e-08
+ -5.1813e-07i 
51 :  0.0000e+00 + 0.0000e+00i                                          
0.0000e+00 + 0.0000e+00i                                      0.0000e+00 +
0.0000e+00i 
52 : -1.6315e-08 + 5.1916e-07i                                           
3.2598e-08 + -5.1814e-07i                                       -4.8785e-08
+ 5.1609e-07i 
53 :  6.5197e-08 + -1.0363e-06i                                          
-9.7667e-08 +  1.0332e-06i                                       1.2988e-07
+ -1.0281e-06i

Could you let me have some idea to avoid this abnormal point like a
singularity?


Sincerely,
Hyuncheong

ps. For better communication, I'm putting a MATLAB (MS Windows) and my C
(gcc 4.1.2 under Linux) source codes here. Please take a look into it.

------------------------------------------------------------------------------------------------------------------
% MATLAB code
%spacer, mesh grid 
[x_ y_] = meshgrid(linspace(-10000,10000,100),linspace(-10000,10000,100));

%
stdev = 100; 
g_    = 1 / (3.141592 * stdev .^ 2) * exp(-(x_ .^ 2 + y_ .^ 2) ./ stdev
^2);

%
g_fft = fft2(g_);

%
g_ifft = ifft2(g_fft);
------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------
// C source code
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>

#include <stdio.h>
#include <fftw3.h>
#include <math.h>

#define MAXSIZE	 100

#define TRUE	 true
#define FALSE	 false 

using namespace std;

int main (int argc, char** argv) 
{

// Lineary space vector
double spacer[MAXSIZE];

for (int i = 0; i < MAXSIZE; i++) {

spacer[i] = (-10000. + 20000./(MAXSIZE-1)*i);
}

// Mesh grid
double x_[MAXSIZE][MAXSIZE];
double y_[MAXSIZE][MAXSIZE];

for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {

x_[i][j] = spacer[j];
y_[i][j] = spacer[i];
}
}

// Gaussian convulution
double stdev = 100.;
double g_[MAXSIZE][MAXSIZE];

for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {

g_[i][j] = (1./(M_PI*stdev*stdev))
   *exp(-(x_[i][j]*x_[i][j] + y_[i][j]*y_[i][j])/(stdev*stdev));
}
}

// FFT g (FWD)
 fftw_complex *fft_in_g;
 fftw_complex *fft_out_g;

 fft_in_g  =
(fftw_complex*)fftw_malloc(MAXSIZE*MAXSIZE*sizeof(fftw_complex));
 fft_out_g	=
(fftw_complex*)fftw_malloc(MAXSIZE*MAXSIZE*sizeof(fftw_complex));

for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {

fft_in_g[i*MAXSIZE+j][0] = g_[i][j];
fft_in_g[i*MAXSIZE+j][1] = 0;
}
}

fftw_plan plan_fft_g;

plan_fft_g = fftw_plan_dft_2d(MAXSIZE, MAXSIZE, 
      fft_in_g, fft_out_g, 
      FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(plan_fft_g);

double g_fft_Re[MAXSIZE][MAXSIZE];
double g_fft_Im[MAXSIZE][MAXSIZE];

for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {

g_fft_Re[i][j] = fft_out_g[i*MAXSIZE+j][0];
g_fft_Im[i][j] = fft_out_g[i*MAXSIZE+j][1];
printf("%4.4e+%4.4ei ", g_fft_Re[i][j], g_fft_Im[i][j]);
}
cout << endl;
}

//
fftw_destroy_plan(plan_fft_g);
fftw_free(fft_in_g);
fftw_free(fft_out_g);

// IFFT g_ (BWD)
 fftw_complex *fft_in_ig;
 fftw_complex *fft_out_ig;

 fft_in_ig  =
(fftw_complex*)fftw_malloc(MAXSIZE*MAXSIZE*sizeof(fftw_complex));
 fft_out_ig	=
(fftw_complex*)fftw_malloc(MAXSIZE*MAXSIZE*sizeof(fftw_complex));

for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {

fft_in_ig[i*MAXSIZE+j][0] = g_fft_Re[i][j];
fft_in_ig[i*MAXSIZE+j][1] = g_fft_Im[i][j];
}
}

fftw_plan plan_ifft_ig;

plan_ifft_ig = fftw_plan_dft_2d(MAXSIZE, MAXSIZE, 
fft_in_ig, fft_out_ig, 
FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(plan_ifft_ig);

double c_ifft_Re[MAXSIZE][MAXSIZE];
double c_ifft_Im[MAXSIZE][MAXSIZE];

cout << endl;
cout << endl;
for (int i = 0; i < MAXSIZE; i++) {
for (int j = 0; j < MAXSIZE; j++) {


c_ifft_Re[i][j] = fft_out_ig[i*MAXSIZE+j][0] / (MAXSIZE*MAXSIZE);
/*normalized*/
c_ifft_Im[i][j] = fft_out_ig[i*MAXSIZE+j][1] / (MAXSIZE*MAXSIZE);

printf("%4.4e ", c_ifft_Re[i][j]);
}
cout << endl;
}

//
fftw_destroy_plan(plan_ifft_ig);
fftw_free(fft_in_ig);
fftw_free(fft_out_ig);

//
return 0;
}





On Feb 17, 9:38&#2013266080;am, "hyun.ha" <hyuncheong.ha@n_o_s_p_a_m.gmail.com>
wrote:
> 49 : &#2013266080;6.51969993785790e-08 + 1.03627663030827e-06i > -3.25984630124937e-08 - 1.03729991235754e-06i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 6.61744490042422e-22 > + 1.03627197122179e-06i > 50 : -1.63152920904142e-08 - 5.19161012253523e-07i > -4.35097002202893e-22 + 5.19160428141898e-07i &#2013266080; &#2013266080; &#2013266080;1.62830243797501e-08 - > 5.18134237051512e-07i > 51 : &#2013266080;1.69406589448395e-21 + 0.00000000000000i > -1.69239446844913e-21 - 5.31856388157782e-23i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 1.68738678669962e-21 > + 1.06161378213412e-22i > 52 : -1.63152920904142e-08 + 5.19161012253523e-07i > 3.25983530020901e-08 - 5.18135983627880e-07i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; -4.87848114810276e-08 + > 5.16089397930887e-07i > 53 : &#2013266080;6.51969993785790e-08 - 1.03627663030827e-06i > -9.76667377967905e-08 + 1.03320616350072e-06i &#2013266080; &#2013266080; &#2013266080; 1.29879317005283e-07 - > 1.02810065720927e-06i > > from FFTW3, > > 49 : &#2013266080;6.5197e-08 + 1.0363e-06i > -3.2598e-08 + -1.0373e-06i > -4.2352e-22 + 1.0363e-06i > 50 : -1.6315e-08 + - 5.1916e-07i > -4.2352e-22 + &#2013266080;5.1916e-07i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 1.6283e-08 > + -5.1813e-07i > 51 : &#2013266080;0.0000e+00 + 0.0000e+00i > 0.0000e+00 + 0.0000e+00i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0.0000e+00 + > 0.0000e+00i > 52 : -1.6315e-08 + 5.1916e-07i > 3.2598e-08 + -5.1814e-07i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; -4.8785e-08 > + 5.1609e-07i > 53 : &#2013266080;6.5197e-08 + -1.0363e-06i > -9.7667e-08 + &#2013266080;1.0332e-06i &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 1.2988e-07 > + -1.0281e-06i
These look the same to me (at least, to the precision you are printing). What is the problem? (Note that the row labeled "51" in the Matlab output is the same as zero up to roundoff accuracy, so it matches your C call to FFTW.)
On 2/17/2012 9:38 AM, hyun.ha wrote:
> Dear friends, > > Here I have a problem to re-produce same results with MATLAB and FFTW > during 2 dimensional FFT. > > I'm trying to do FFT with Gaussian convolution which coordinates varies > from negative to positive, FFTW3 seems not to calculate the values > correctly. It happens only near center place. > And it affects IFFT and it's shifted results seriously in my calculation.
Matlab uses FFTW code. Any difference is in result is dur to your setup. Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
Dear Steven, Jerry

Thanks a lot. I found there were some mistake at freeing memory allocations
inside my full code, especially inverse FFT part. It was too fast to
preceed following calculation. :)