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. :)
Reply by Jerry Avins●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 stevenj●February 17, 20122012-02-17
On Feb 17, 9:38�am, "hyun.ha" <hyuncheong.ha@n_o_s_p_a_m.gmail.com>
wrote:
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.)
Reply by hyun.ha●February 17, 20122012-02-17
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;
}