[negar68]

Hello!

I am trying to calculate the partial derivative of the function Sin(x)Sin(y) along the x-direction using Fourier Transform (FFTW library in C++). I take my function from real to Fourier space (fftw_plan_dft_r2c_2d), multiply the result by -i * kappa array (wavenumber/spatial frequency), and then transform the result back into real space (fftw_plan_dft_c2r_2d). Finally, I divide the result by the size of my 2D array Nx*Ny (I need to do this based on the documentation). However, the resulting function is scaled up by some value and the amplitude is not 1.

 
    #define Nx 360
    #define Ny 360
    #define REAL 0
    #define IMAG 1
    #define Pi 3.14
    
    #include <stdio.h>
    #include <math.h>
    #include <fftw3.h>
    #include <iostream>
    #include <iostream>
    #include <fstream>
    #include <iomanip> 
    
    int main() {
    
    int i, j;
    int Nyh = Ny / 2 + 1;
    
    double k1;
    double k2;
    double *signal;
    double *result2;
    double *kappa1;
    double *kappa2;
    double df;
    fftw_complex *result1;
    fftw_complex *dfhat;
    
    signal = (double*)fftw_malloc(sizeof(double) * Nx * Ny );
    result2 = (double*)fftw_malloc(sizeof(double) * Nx * Ny );
    kappa1 = (double*)fftw_malloc(sizeof(double) * Nx * Nyh );
    result1 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * Nx * Nyh);
    dfhat = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * Nx * Nyh);
    
    for (int i = 0; i < Nx; i++){
            if ( i  < Nx/2 + 1)
    		k1 = 2 * Pi * (-Nx + i) / Nx;
            else 
    		k1 = 2 * Pi * i / Nx;
           for (int j = 0; j < Nyh; j++){
                    kappa1[j + Nyh * i] = k1;
    }
    }
    
    fftw_plan plan1 = fftw_plan_dft_r2c_2d(Nx,Ny,signal,result1,FFTW_ESTIMATE);
    fftw_plan plan2 = fftw_plan_dft_c2r_2d(Nx,Ny,dfhat,result2,FFTW_ESTIMATE);
    
            for (int i=0; i < Nx; i++){
            	for (int j=0; j < Ny; j++){
    	          double x = (double)i / 180 * (double) Pi;
                      double y = (double)j / 180 * (double) Pi;
                      signal[j + Ny * i] = sin(x) * sin(y);
            	}
                    }
    
    fftw_execute(plan1);
    
    for (int i = 0; i < Nx; i++) {
    	for (int j = 0; j < Nyh; j++) {
    		dfhat[j + Nyh * i][REAL] = 0;
    		dfhat[j + Nyh * i][IMAG] = 0;
    		dfhat[j + Nyh * i][IMAG] = -result1[j + Nyh * i][REAL] * kappa1[j + Nyh * i];
    		dfhat[j + Nyh * i][REAL] = result1[j + Nyh * i][IMAG] * kappa1[j + Nyh * i];
        }
    }
    
    
    		fftw_execute(plan2);
    
    		for (int i = 0; i < Nx; i++) {
    			for (int j = 0; j < Ny; j++) {
    				result2[j + Ny * i] /= (Nx*Ny);
    			}
    		}
    
                   for (int j = 0; j < Ny; j++) {
                   for (int i = 0; i < Nx; i++) {
                   std::cout <<std::setw(15)<< result2[j + Ny * i];
                            }
                   std::cout << std::endl;
                    }
    
    fftw_destroy_plan(plan1);
    fftw_destroy_plan(plan2);
    
    fftw_free(result1);
    fftw_free(dfhat);
    
    return 0;
    }

I have defined my spatial frequency/wavenumber array to have Nx * Ny/2+1 elements and have split the array at the center (as all references suggest). But it does not work. I would appreciate any help!