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 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!

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)