I have the MATLAB code:
Nx = 10;
Ny = 10;
Lx = 2*pi;
ygl = -cos(pi*(0:Ny)/Ny)'; %Gauss-Lobatto chebyshev points
x = (0:Nx-1)/Nx*2*pi;
%make mesh
[X,Y] = meshgrid(x,ygl);
A = 2*pi / Lx;
u = sin( (2*pi / Lx) * X);
uh = fft(u)
This gives the following output:
0 6.4656 10.4616 10.4616 6.4656 0.0000 -6.4656 -10.4616 -10.4616 -6.4656
0 0.0000 0 0 0.0000 0.0000 -0.0000 0 0 0
0 0.0000 0 0 0.0000 0.0000 0 0 0 -0.0000
0 0 0 0 0.0000 0.0000 0 0 0 0
0 0.0000 0 0 0.0000 0.0000 -0.0000 0 0 0
0 0.0000 0 0 0.0000 0.0000 0 0 0 -0.0000
0 0.0000 0 0 0.0000 0.0000 0 0 0 -0.0000
0 0.0000 0 0 0.0000 0.0000 -0.0000 0 0 0
0 0 0 0 0.0000 0.0000 0 0 0 0
0 0.0000 0 0 0.0000 0.0000 0 0 0 -0.0000
0 0.0000 0 0 0.0000 0.0000 -0.0000 0 0 0
And the C code using FFTW:
static const int nx = 10;
static const int ny = 10;
static const int nyk = ny/2 1;
double Lx = 2 * M_PI;
double A = (2 * M_PI)/Lx;
double *XX;
XX = (double*) fftw_malloc((nx*(ny 1))*sizeof(double));
memset(XX, 42, (nx*(ny 1))* sizeof(double));
double *YY;
YY = (double*) fftw_malloc((nx*(ny 1))*sizeof(double));
memset(YY, 42, (nx*(ny 1))* sizeof(double));
double *u;
u = (double*) fftw_malloc((((ny 1)*nx))*sizeof(double));
fftw_complex *uh;
uh = (fftw_complex*) fftw_malloc(((ny 1)*nyk)*sizeof(fftw_complex));
memset(uh, 42, ((ny 1)*nyk)* sizeof(fftw_complex));
for(int i = 0; i< nx 1; i ){
for(int j = 0; j< ny; j ){
XX[i (ny 1)*j] = (j)*2*M_PI/nx;
YY[i (ny 1)*j] = -1. * cos(((i) * M_PI )/ny);
u[i (ny 1)*j] = sin(A * XX[i (ny 1)*j]);
}
}
fftw_plan r2c1;
r2c1 = fftw_plan_dft_r2c_1d(nx , &u[0], &uh[0], FFTW_ESTIMATE);
fftw_execute(r2c1);
fftw_destroy_plan(r2c1);
fftw_cleanup();
//print uh
This give the output:
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
(0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000), (0.0000,0.0000),
Why are these two results are NOT matching?? I understand there are some differences between FFT in MATLAB and C but this is a simple 1D FFT and it's WAY off than the MATLAB results.
NOTE: I added the results of MATLAB and C for comparison ad requested by someone in the comments.
CodePudding user response:
Was able to fix this by taking the FFT of each columns as someone suggested in the comments:
{
cplx_buffer out = my_fftw_allocate_cplx(nyk, nx);
fftw_execute(fftw_plan_many_dft_r2c(1, &in.rows, in.cols, in.a, &in.rows, in.cols, 1, out.a, &in.rows, out.cols, 1, FFTW_ESTIMATE));
std::cout << "DFT of input, column-wise\n"; print_matrix(out);
}