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; }
2D FFTW3 vs. MATLAB difference only near center place
Started by ●February 17, 2012
Reply by ●February 17, 20122012-02-17
On Feb 17, 9:38�am, "hyun.ha" <hyuncheong.ha@n_o_s_p_a_m.gmail.com> wrote:> 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-06iThese 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.)
Reply by ●February 18, 20122012-02-18
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. �����������������������������������������������������������������������
Reply by ●February 23, 20122012-02-23