Jamie
Reputation: 519
Taking FFTW in MATLAB vs C/C++: Different outputs
I am exploring the differences between fftw in MATLAB and C/C++. I think this maybe a straightforward question but for some reason I don't see it.
Assume the following MATLAB code:
Nx = 8;
Ny = 8;
Lx =16;
%-----------
xi_x = (2*pi)/Lx;
xi = ((0:Nx-1)/Nx)*(2*pi);
x = xi/xi_x;
ylow = 0;
yupp =6;
ygl = cos ( pi * ( 0 : Ny ) / Ny )';
ygl = (1/2)*(((yupp-ylow)*ygl) + (yupp+ylow));
[X,Y] = meshgrid(x,ygl);
%define function u:
u = (Y-ylow) .* (Y-yupp) .* sin( (2*pi / Lx) * X);
uh = fft(u,[],2);
%pad uh with zeros
uhp=[uh(:,1:Nx/2) zeros(Ny,(m-Nx)) uh(:,Nx/2+1:Nx)];
%take ifft
up=real(ifft(uhp,[],2));
u1 = real(ifft(uh,[],2)); %same as u
My question above since both uh and uhp are the same values with extra zero padding, why is the inverse of fft of these two different? Meaning if uh is:
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 5.2721i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 5.2721i
-0.0000 + 0.0000i 0.0000 +18.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -18.0000i
-0.0000 + 0.0000i 0.0000 +30.7279i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -30.7279i
-0.0000 + 0.0000i 0.0000 +36.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -36.0000i
-0.0000 + 0.0000i 0.0000 +30.7279i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -30.7279i
-0.0000 + 0.0000i 0.0000 +18.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -18.0000i
0.0000 + 0.0000i 0.0000 + 5.2721i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 5.2721i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
and uhp is:
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 5.2721i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 5.2721i
-0.0000 + 0.0000i 0.0000 +18.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -18.0000i
-0.0000 + 0.0000i 0.0000 +30.7279i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -30.7279i
-0.0000 + 0.0000i 0.0000 +36.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -36.0000i
-0.0000 + 0.0000i 0.0000 +30.7279i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -30.7279i
-0.0000 + 0.0000i 0.0000 +18.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 -18.0000i
0.0000 + 0.0000i 0.0000 + 5.2721i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i 0.0000 - 5.2721i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
Then why the real(ifft(uh,[],2)) equals:
0 0 0 0 0 0 0 0
0 -0.9320 -1.3180 -0.9320 -0.0000 0.9320 1.3180 0.9320
-0.0000 -3.1820 -4.5000 -3.1820 -0.0000 3.1820 4.5000 3.1820
0 -5.4320 -7.6820 -5.4320 -0.0000 5.4320 7.6820 5.4320
-0.0000 -6.3640 -9.0000 -6.3640 -0.0000 6.3640 9.0000 6.3640
0 -5.4320 -7.6820 -5.4320 -0.0000 5.4320 7.6820 5.4320
-0.0000 -3.1820 -4.5000 -3.1820 -0.0000 3.1820 4.5000 3.1820
0 -0.9320 -1.3180 -0.9320 -0.0000 0.9320 1.3180 0.9320
0 0 0 0 0 0 0 0
and up=ifft(uhp,[],2); is:
0 0 0 0 0 0 0 0 0 0 0 0
0 -0.4393 -0.7610 -0.8787 -0.7610 -0.4393 -0.0000 0.4393 0.7610 0.8787 0.7610 0.4393
-0.0000 -1.5000 -2.5981 -3.0000 -2.5981 -1.5000 -0.0000 1.5000 2.5981 3.0000 2.5981 1.5000
0 -2.5607 -4.4352 -5.1213 -4.4352 -2.5607 -0.0000 2.5607 4.4352 5.1213 4.4352 2.5607
-0.0000 -3.0000 -5.1962 -6.0000 -5.1962 -3.0000 -0.0000 3.0000 5.1962 6.0000 5.1962 3.0000
0 -2.5607 -4.4352 -5.1213 -4.4352 -2.5607 -0.0000 2.5607 4.4352 5.1213 4.4352 2.5607
-0.0000 -1.5000 -2.5981 -3.0000 -2.5981 -1.5000 -0.0000 1.5000 2.5981 3.0000 2.5981 1.5000
0 -0.4393 -0.7610 -0.8787 -0.7610 -0.4393 -0.0000 0.4393 0.7610 0.8787 0.7610 0.4393
0 0 0 0 0 0 0 0 0 0 0 0
The reason I ask b/c in C/C++ when padding my arrays similarly and taking the ifft my code returns the same output, so real(ifft(uh,[],2)) = real(ifft(uhp,[],2))!
Here's my C/C++ code:
static const int nx = 8;
static const int ny = 8;
static const int ncomp = 2;
static const int nxk = nx/2 + 1;
static const int Mx = 3*nx/2, My= 3*ny/2;
double makeSpatialMesh2DTrans(double dx, double dy, double ylow, double yupp,double xarr[], double yarr[]);
int xy(int x, int y, int number_of_cols);
class MyFftwClass {
public:
MyFftwClass(void);
~MyFftwClass(void);
void r2cexecute(double rArr[], double cArr[][ncomp]);
private:
static fftw_plan s_plan;
double *m_buffer_in;
fftw_complex *m_buffer_out;
};
class MyiFftwClass { //inverse
public:
MyiFftwClass(void);
~MyiFftwClass(void);
void c2rexecute(double cArr[][ncomp],double rArr[]);
private:
static fftw_plan s_plan1;
fftw_complex *m_buffer_in1;
double *m_buffer_out1;
};
MyFftwClass r2c;
MyiFftwClass c2r;
int main(){
double Lx = 16.;
double ylow = 0.;
double yupp = 6.;
double Ly = (yupp-ylow);
double dx = Lx / nx;
double dy = Ly / ny;
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 kmax = makeSpatialMesh2DTrans(dx, dy, ylow, yupp, XX, YY);
double *uFun;
uFun = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
double A = 2*M_PI/Lx;
for(int y = 0; y< ny+1; y++){
for(int x = 0; x< nx; x++){
uFun[xy(x, y, nx)] = ( (YY[xy(x, y, nx)]-ylow) * (YY[xy(x, y, nx)]-yupp) * sin(A*XX[xy(x, y, nx)]));
}
}
}
fftw_complex *uFunk;
uFunk = (fftw_complex*) fftw_malloc((nx*(ny + 1)) * sizeof(fftw_complex));
memset(uFunk, 42, (nx*(ny + 1))* sizeof(fftw_complex));
r2c.r2cexecute(uFun,uFunk);
double N[2] = {nx,ny};
double M[2] = {Mx,My};
fftw_complex *uFunk_pad;
uFunk_pad = (fftw_complex*) fftw_malloc((Mx*My)* sizeof(fftw_complex));
//Add zero padding to the function u
std::cout << std::fixed << std::setprecision(3);
for(int y=0; y<ny+1; y++){
for(int x=0; x<Mx; x++){
if(x<nxk && y<ny){
uFunk_pad[xy(x, y, nx)][0] = uFunk[xy(x, y, nx)][0];
uFunk_pad[xy(x, y, nx)][1] = uFunk[xy(x, y, nx)][1];
}
else{
uFunk_pad[xy(x, y, nx)][0] = 0.0;
uFunk_pad[xy(x, y, nx)][1] = 0.0;
}
std::cout << std::setw(6) << std::complex{uFunk_pad[xy(x, y, nx)][0], uFunk_pad[xy(x, y, nx)][1]} << ' ';
}
std::cout << "\n";
}
//Take ifft of padded inputs
double *f_pad;
f_pad = (double*) fftw_malloc((Mx*My)*sizeof(double));
memset(f_pad, 42, (Mx*My)* sizeof(double));
c2r.c2rexecute(uFunk_pad,f_pad);
std::cout << std::fixed << std::setprecision(3);
for(int y=0; y<ny+1; y++){
for(int x=0; x<Mx; x++){
std::cout << std::setw(6) << f_pad[xy(x, y, nx)] << ' ';
}
std::cout << "\n";
//free stuff
}
int xy(int x, int y, int number_of_cols) { return x + y * number_of_cols; }
double makeSpatialMesh2DTrans(double dx, double dy, double ylow, double yupp,double xarr[], double yarr[]){
for(int y = 0; y< ny+1; y++){
for(int x = 0; x< nx; x++){
xarr[xy(x, y, nx)] = x*dx;
yarr[xy(x, y, nx)] = 1. * cos(((y) * M_PI )/ny);
yarr[xy(x, y, nx)] = (1./2)*(((yupp-ylow)*yarr[xy(x, y, nx)]) + (yupp+ylow));
}
}
if (dx < dy){
return 1/(2*dx);
}else {
return 1/(2*dy);
}
}
MyiFftwClass::MyiFftwClass(void)
{
m_buffer_in1 = fftw_alloc_complex(nx * (ny + 1));
m_buffer_out1 = fftw_alloc_real((nx * (ny + 1)));
if (!(m_buffer_in1 && m_buffer_out1))
{
throw new std::runtime_error("Failed to allocate memory!");
}
if (!s_plan1)
{
s_plan1 = fftw_plan_many_dft_c2r(1,&nx, (ny+1),m_buffer_in1,nullptr,1,(nx),m_buffer_out1,nullptr,1, (nx), FFTW_ESTIMATE);
if (!s_plan1)
{
throw new std::runtime_error("Failed to create plan!");
}
}
}
MyiFftwClass::~MyiFftwClass(void)
{
fftw_free(m_buffer_in1);
fftw_free(m_buffer_out1);
}
void MyiFftwClass::c2rexecute(double cArr[][ncomp],double rArr[])
{
memcpy(m_buffer_in1,cArr, sizeof(fftw_complex) * (nx * (ny + 1)));
fftw_execute_dft_c2r(s_plan1, m_buffer_in1, m_buffer_out1);
memcpy(rArr, m_buffer_out1,sizeof(double) * (nx * (ny + 1)));
// renormalize:
for(int y = 0; y < ny + 1; y++)
for(int x = 0; x < nx; x++)
rArr[xy(x, y, nx)] /= nx;
}
fftw_plan MyiFftwClass::s_plan1 = NULL;
MyFftwClass::MyFftwClass(void)
{
m_buffer_in = fftw_alloc_real((nx * (ny + 1)));
m_buffer_out = fftw_alloc_complex(nx * (ny + 1));
if (!(m_buffer_in && m_buffer_out))
{
throw new std::runtime_error("Failed to allocate memory!");
}
if (!s_plan)
{
s_plan = fftw_plan_many_dft_r2c(1,&nx, (ny+1),m_buffer_in,nullptr,1,(nx),m_buffer_out,nullptr,1, (nx), FFTW_ESTIMATE);
if (!s_plan)
{
throw new std::runtime_error("Failed to create plan!");
}
}
}
MyFftwClass::~MyFftwClass(void)
{
fftw_free(m_buffer_in);
fftw_free(m_buffer_out);
}
void MyFftwClass::r2cexecute(double rArr[], double cArr[][ncomp])
MyFftwClass::execute(double *const in, fftw_complex *const out)
{
memcpy(m_buffer_in, rArr, sizeof(double) * (nx * (ny + 1)));
fftw_execute_dft_r2c(s_plan, m_buffer_in, m_buffer_out);
memcpy(cArr, m_buffer_out, sizeof(fftw_complex) * (nx * (ny + 1)));
}
fftw_plan MyFftwClass::s_plan = NULL;
Outputs of the C/C++ are:
For uFunk the output:
(0.000,0.000) (0.000,-0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(0.000,0.000) (0.000,5.272) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,18.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,30.728) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,36.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,30.728) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,18.000) (-0.000,0.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,5.272) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(0.000,0.000) (0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
and output of the padded function uFunk_pad:
(0.000,0.000) (0.000,-0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(0.000,0.000) (0.000,5.272) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,18.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,30.728) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,36.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,30.728) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,18.000) (-0.000,0.000) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(-0.000,0.000) (0.000,5.272) (-0.000,0.000) (-0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
(0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
Then when taking the inverse of padded function uFunk_paddedit returns the same value as real(ifft(uh,[],2)) in MATLAB function and NOT real(ifft(uhp,[],2)):
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.932 -1.318 -0.932
0.000 -0.932 -1.318 -0.932 -0.000 0.932 1.318 0.932 -0.000 -3.182 -4.500 -3.182
-0.000 -3.182 -4.500 -3.182 -0.000 3.182 4.500 3.182 -0.000 -5.432 -7.682 -5.432
-0.000 -5.432 -7.682 -5.432 -0.000 5.432 7.682 5.432 0.000 -6.364 -9.000 -6.364
0.000 -6.364 -9.000 -6.364 -0.000 6.364 9.000 6.364 0.000 -5.432 -7.682 -5.432
0.000 -5.432 -7.682 -5.432 -0.000 5.432 7.682 5.432 0.000 -3.182 -4.500 -3.182
0.000 -3.182 -4.500 -3.182 -0.000 3.182 4.500 3.182 0.000 -0.932 -1.318 -0.932
0.000 -0.932 -1.318 -0.932 -0.000 0.932 1.318 0.932 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Why is my C/C++ output not matching the one from MATLAB? Am I misunderstanding something fundamental? Thanks
Upvotes: 3
Views: 65
Answers (0)
Related Questions
- 1D FFT MATLAB results are different from FFTW in C++
- Can code that is valid in both C and C++ produce different behavior when compiled in each language?
- Different results with fft2 (matlab) and fftw (C)
- fftw/c++ computes fft wrong, compared to matlab
- FFTW vs Matlab FFT
- Matlab FFT and FFTW
- Incorrect result after serializing and deserializing time_t variable