Partial Derivative Using FFTW in 2D (C++)

[negarahani]

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 NxNy (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 kappa 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!