Home > Software design >  1D FFT MATLAB results are different from FFTW in C
1D FFT MATLAB results are different from FFTW in C

Time:06-29

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);
  }

  • Related