C++

Forum

 
1D FFTs of an Eigen tensor

JamieAl (236)
So I am trying to take an FFT along a specific direction of an Eigen tensor, similarly to my code in MATLAB.

In MATLAB:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 
Nx = 8;  
Ny = 6; 
Nz= 8;
Lx =6;
Ly = 6;
xi_x =  (2*pi)/Lx;
yi_y =  (2*pi)/Ly;
xi  = ((0:Nx-1)/Nx)*(2*pi);
yi  = ((0:Ny-1)/Ny)*(2*pi);
x =  xi/xi_x;
y =  yi/yi_y;

zlow = 0;%a
zupp =6; %b
Lz = (zupp-zlow);
zgl = (1/2)*(((zupp-zlow)*zgl) + (zupp+zlow));
[X,Z,Y] = meshgrid(x,zgl,y);
%ICs
A = 2*pi / Lx;
B = 2*pi / Ly;
u = (Z-zlow) .* (Z-zupp) .* sin(A*X).* sin(B*Y);
%Take FFT along only X and Y direction
uh =fft(fft(u,[],2),[],3);

The meshgrid function in MATLAB works in such a way that this line [X,Z,Y] = meshgrid(x,zgl,y); returns a 3d grid with a size of z-by-x-by-y. Then when taking my FFT along X which is the second element:
 
 
uhx =fft(u,[],2);

and along Y which is 3rd element:
 
 
uhx =fft(u,[],3);


Now, I am trying to do something similar in C++. I can already perform the following in C++:
 
 
uhx =fft(u,[],2);

where I will be taking the FFT along the columns which is nx. But, is there a way for me to do this uh =fft(fft(u,[],2),[],3); in C++ without having to take the FFTs separately twice? Or I am not sure maybe something similar.

For example, I am using Eigen tensor here for simplicity (prints 3D matrices nicely):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 
static const int nx = 8; 
static const int ny = 6;  
static const int nz = 8; 
using namespace std;
using namespace Eigen;

double SpatialMesh3DTrans(double dx,double dy, double dz, double zlow, double zupp,Eigen::Tensor<double, 3>& xarr, Eigen::Tensor<double, 3>& yarr, Eigen::Tensor<double, 3>& zarr);
int main(){

double Lx = 6.;
double zlow = 0.;
double zupp = 6.;   
double Ly = 6.;
double Lz = (zupp-zlow);
double dx = Lx / nx;
double dy = Ly / ny;
double dz = Lz / nz;

Eigen::Tensor<double, 3> eXX((nz+1),nx,ny); 
eXX.setZero(); 

Eigen::Tensor<double, 3> eYY((nz+1),nx,ny);   //tensor(row,col,matrix)
eYY.setZero(); 
Eigen::Tensor<double, 3> eZZ((nz+1),nx,ny);   
eZZ.setZero(); 

double kmax = SpatialMesh3DTrans(dx,dy,dz,zlow,zupp,eXX,eYY,eZZ);

Eigen::Tensor<double, 3> uFun((nz+1),nx,ny);   //tensor(row,col,matrix)
uFun.setZero();
 
    double A = 2*EIGEN_PI/Lx;
    double A1 = 2*EIGEN_PI/Ly;
   
    	for(int x = 0; x< nx; x++){//cols
		    for(int y = 0; y< ny; y++){
                for(int z = 0; z< nz+1; z++){
                uFun(z,x,y) = ( (eZZ(z,x,y)-zlow) * (eZZ(z,x,y)-zupp) * sin(A*eXX(z,x,y) ) * sin(A1 *eYY(z,x,y) ) );
                
            }
        }
    }
}

double SpatialMesh3DTrans(double dx,double dy, double dz, double zlow, double zupp,Eigen::Tensor<double, 3>& xarr, Eigen::Tensor<double, 3>& yarr, Eigen::Tensor<double, 3>& zarr){
	for(int x = 0; x< nx; x++){
		for(int y = 0; y< ny; y++){
            for(int z = 0; z< nz+1; z++){
                xarr(z,x,y)  = x*dx;
                yarr(z,x,y)  = y*dy;
			    zarr(z,x,y)  = 1. * cos(((z) * EIGEN_PI )/nz); 
                zarr(z,x,y) = (1./2)*(((zupp-zlow)*zarr(z,x,y) ) + (zupp+zlow));
            }
		}		
	}
	
		if (dx < dy){
			return 1/(2*dx);
		}else {
			return 1/(2*dy);
		}
}

Then I can print the "first" matrix from that 3D tensor which is equivalent to uFun(:,:,1) in MATLAB:
1
2
3
4
 
std::array<long,3> offsetA = {0,0,0};        //Starting point:(row,column,matrix)
std::array<long,3> extentA = {(nz+1),nx,1}; //Finish point:(row,column,matrix)
std::array<long,2> shapeA = {(nz+1),(nx)};
std::cout << uFun.slice(offsetA, extentA).reshape(shapeA) << std::endl;  //Extract slice and reshape it into a nz+1xnx matrix.

This returns:

     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0


So, for the given uFun tensor how can I take the FFT along the X and Y direction without having something too messy I guess?
Last edited on
JamieAl (236)
So, in the past I have been able to take the full 3D FFT of a tensor and worked perfectly with fftw_plan_dft_3d. I noticed from my very old posts there were some good suggestions that I will finally attempt to get working. For example, @mbozzi recommended using the fftw_plan_dft_r2c to do a 2D transforms of size nx*ny only of the Eigen matrix and I think I will test that.
JamieAl (236)
This is the attempt I made which is not really working:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
void r2cfft2d(Eigen::Tensor<double, 3>& rArr, Eigen::Tensor<std::complex<double>, 3>& cArr){	
	
			double* input_array = fftw_alloc_real(nx*ny*(nz+1));
          
			memcpy(input_array, rArr.data(), (nx*ny*(nz+1))* sizeof(double));

			fftw_complex *output_array;
			output_array = (fftw_complex*) fftw_malloc((nx*ny*(nz+1)) * sizeof(fftw_complex));

			int constexpr n[]{ny,nx};

			fftw_plan forward = fftw_plan_many_dft_r2c(2,n, (nz+1), input_array,n,1,(ny*nx), output_array,n,1,(ny*nx), FFTW_MEASURE);
			
            fftw_execute(forward);
		memcpy(cArr.data(),output_array, (nx*ny*(nz+1)) * sizeof(fftw_complex));
		fftw_free(input_array);
		fftw_free(output_array);
	
	

}


I think maybe I have to pass the tensor differently and not just use memcpy or use a different plan.
JamieAl (236)
Ok, my second attempt looks better. Basically, I read somewhere that one of the ways to deal with my issue is to extract a 2D tensor out of a 3D tensor and directly apply the FFT to it. What I have been able to do is extract the first 2D matrix out of my 3D tensor with nz+1 rows and nx columns and perform a 2D FFT on that and get a correct result. Now, two issues: 1) how can I easily iterate over all the matrices inside of the 3D tensor. For example, I have Tensor(row,col,matrix) and I have been able to take 2D FFT on Tensor(row,col,1) i.e. the first matrix only. (This can be dealt with easily I think)
2. For some reason applying fftw_plan_many_dft_r2c does not work and returns zeros? Does that mean I have to use a dummy variable to store my matrix so it's not overwritten?

Example 1: 2D FFT of 2D matrix extracted out of 3D tensor:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
//Extract 2D matrix out of 3D tensor 
Eigen::Tensor<double, 2> uFun2D((nz+1),nx);   //tensor(row,col,matrix)
uFun2D.setZero();
std::array<long,3> offsetA = {0,0,1};        //Starting point:(row,column,matrix)
std::array<long,3> extentA = {(nz+1),nx,1}; //Finish point:(row,column,matrix)
std::array<long,2> shapeA = {(nz+1),(nx)};
uFun2D = uFun.slice(offsetA, extentA).reshape(shapeA);
std::cout<<uFun2D<<endl;//correct
//Apply 2D FFT directly on 2D matrix
Eigen::Tensor<std::complex<double>, 2> uFun2Dk((nz+1),nx);   //tensor(row,col,matrix)
uFun2Dk.setZero();
	
r2cfft2d(uFun2D,uFun2Dk);
std::cout<<uFun2Dk<<endl;//weird output does NOT match MATLAB

Eigen::Tensor<double, 2> Out((nz+1),nx);   //tensor(row,col,matrix)
Out.setZero();

c2rfft2d(LujK,LujOut);
std::cout<<Out<<endl;//correct inverse

the FFT functions look like:
r2cfft2d:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
void r2cfft2d(Eigen::Tensor<double, 2>& rArr, Eigen::Tensor<std::complex<double>, 2>& cArr){	

fftw_complex *input_array;
input_array = (fftw_complex*) fftw_malloc((nx*(nz+1)) * sizeof(fftw_complex));

memcpy(input_array, rArr.data(), (nx*(nz+1))* sizeof(double));
fftw_complex *output_array;
output_array = (fftw_complex*) fftw_malloc((nx*(nz+1)) * sizeof(fftw_complex));
			
fftw_plan forward = fftw_plan_dft_2d(nx, (nz+1), input_array, output_array, FFTW_FORWARD, FFTW_ESTIMATE);
			
fftw_execute(forward);
	
memcpy(cArr.data(),output_array, (nx*(nz+1))* sizeof(fftw_complex));
fftw_free(input_array);
fftw_free(output_array);
		

}

and c2rfft2d:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
void c2rfft2d(Eigen::Tensor<std::complex<double>, 2>& cArr, Eigen::Tensor<double, 2>& rArr){	
		fftw_complex *input_array;
		input_array = (fftw_complex*) fftw_malloc((nx*(nz+1))* sizeof(fftw_complex));
		memcpy(input_array, cArr.data(), (nx*(nz+1))* sizeof(fftw_complex));
		fftw_complex *output_array;
		output_array = (fftw_complex*) fftw_malloc((nx*(nz+1)) * sizeof(fftw_complex));
		fftw_plan backward = fftw_plan_dft_2d(nx, (nz+1), input_array, output_array, FFTW_BACKWARD, FFTW_ESTIMATE);
		fftw_execute(backward);
		
		memcpy(rArr.data(),output_array, (nx*(nz+1)) * sizeof(double)); 
		Eigen::Tensor<double, 2> dummy((nz+1),nx); 
		dummy =  1.0/((nx*(nz+1))) * rArr; 
		
		fftw_free(input_array);
		fftw_free(output_array);

}

The above code returns an inverse function correctly but the 2D FFT does not necessarily return uh =fft(fft(u,[],2),[],3); obviously.
mbozzi (3933)
The above code returns an inverse function correctly but the 2D FFT does not necessarily return uh =fft(fft(u,[],2),[],3); obviously.

To make sure I understand what you're trying to compute, each matrix is like sheet of paper that spans the z and x axes. The tensor is a stack of papers that grows up along the y axis. And uh should be the FFT through the x and y axes -- i.e., the FFT of each tall slice formed by the individual columns of different matrices. In other words, I think you're considering the indices (z, x, y) to be row z, column x, and "sheet of paper"/matrix y. Am I right?

The above code returns an inverse function correctly but the 2D FFT does not necessarily return uh = fft(fft(u,[],2),[],3); obviously.
Can you give an example input/output for u and uh? So I can tell if the program matches the MATLAB output?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
 
#include "Tensor"
#include "fftw3.h"

#include <complex>
#include <array>
#include <iostream>
#include <iomanip>

static const int nx = 8; 
static const int ny = 6;  
static const int nz = 8; 

using namespace std;
using namespace Eigen;

double SpatialMesh3DTrans(
  double dx,
  double dy,
  double dz,
  double zlow,
  double zupp,
  Tensor<double, 3>& xarr,
  Tensor<double, 3>& yarr,
  Tensor<double, 3>& zarr)
{
  for (int x = 0; x < nx; x++)
  {
    for (int y = 0; y < ny; y++)
    {
      for(int z = 0; z < nz+1; z++)
      {
        xarr(z,x,y) = x*dx;
        yarr(z,x,y) = y*dy;
        zarr(z,x,y) = 1. * cos(((z) * EIGEN_PI) / nz); 
        zarr(z,x,y) = (1./2)*(((zupp-zlow)*zarr(z,x,y) ) + (zupp+zlow));
      }
    }		
  }
	
  if (dx < dy){
    return 1/(2*dx);
  }else {
    return 1/(2*dy);
  }
}

int main()
{
  double const Lx   = 6.;
  double const zlow = 0.;
  double const zupp = 6.;   
  double const Ly   = 6.;
  double const Lz   = zupp - zlow;
  double const dx   = Lx / nx;
  double const dy   = Ly / ny;
  double const dz   = Lz / nz;

  Tensor<double, 3> eXX(nz + 1, nx, ny); 
  Tensor<double, 3> eYY(nz + 1, nx, ny);
  Tensor<double, 3> eZZ(nz + 1, nx, ny);

  Tensor<double, 3> uFun(nz + 1, nx, ny);
  Tensor<std::complex<double>, 3> uFunk(nz + 1, nx, ny);
  Tensor<double, 3> out(nz + 1, nx, ny);
  
  int fft_size[] { ny, nx }; // NOTE(mbozzi): Indices are in reverse order because the tensor is column-major.
  fftw_plan fft_through_each_zx_slice = fftw_plan_many_dft_r2c(
      2,             // int rank,
      fft_size,      // const int *n,
      ny,            // int howmany,
      uFun.data(),   // double *in,
      nullptr,       // const int *inembed,
      1,             // int istride, // adjacent elements within the same zx-slice are adjacent in memory
      (nz + 1) * ny, // int idist,   // adjacent elements within separate zx-slices are separate in memory.
      reinterpret_cast<fftw_complex*>(uFunk.data()), // fftw_complex *out,
      nullptr,       // const int *onembed,
      1,             // int ostride,
      (nz + 1) * ny, // int odist,
      FFTW_ESTIMATE  // unsigned flags
      );
  
  fftw_plan inverse_fft_through_each_zx_slice = fftw_plan_many_dft_c2r(
    2,             // int rank,
    fft_size,      // const int *n,
    ny,            // int howmany,
    reinterpret_cast<fftw_complex*>(uFunk.data()), // fftw_complex *in,
    nullptr,       // const int *inembed,
    1,             // int istride, // adjacent elements within the same zx-slice are adjacent in memory
    (nz + 1) * ny, // int idist,   // adjacent elements within separate zx-slices are separate in memory.
    out.data(),    // fftw_complex *out,
    nullptr,       // const int *onembed,
    1,             // int ostride,
    (nz + 1) * ny, // int odist,
    FFTW_ESTIMATE  // unsigned flags
    );

  double const kmax = SpatialMesh3DTrans(dx,dy,dz,zlow,zupp,eXX,eYY,eZZ);
  
  double const A  = 2*EIGEN_PI / Lx;
  double const A1 = 2*EIGEN_PI / Ly;
   
  for (int x = 0; x < nx; x++)
  { 
    for (int y = 0; y < ny; y++)
    {
      for (int z = 0; z < nz + 1; z++)
      {
        uFun(z,x,y) = (eZZ(z,x,y)-zlow) * (eZZ(z,x,y)-zupp) * sin(A*eXX(z,x,y)) * sin(A1 *eYY(z,x,y));
      }
    }
  }

  std::cout << std::setprecision(2) << std::fixed;
  
  std::cout << uFun.slice(std::array{0, 0, 1}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) << "\n\n"; 
  fftw_execute(fft_through_each_zx_slice);            // uFun --> uFunk
  fftw_execute(inverse_fft_through_each_zx_slice);    // uFunk --> out
  out = (1.0 / (nx * ny)) * out; // normalize
  std::cout << out.slice(std::array{0, 0, 1}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) << '\n';
}

Last edited on
JamieAl (236)
Can you give an example input/output for u and uh? So I can tell if the program matches the MATLAB output?


Yes, so in MATLAB:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
clearvars; clc; close all;

Nx = 8;  
Ny = 6; 
Nz= 8;
Lx =6;
Ly = 6;

%1. Exact Case vs Approximation Case
dx = Lx/Nx; % Need to calculate dx
dy = Ly/Ny;
xi_x =  (2*pi)/Lx;
yi_y =  (2*pi)/Ly;
xi  = ((0:Nx-1)/Nx)*(2*pi);
yi  = ((0:Ny-1)/Ny)*(2*pi);
x =  xi/xi_x;
y =  yi/yi_y;

zlow = 0;%a
zupp =6; %b
Lz = (zupp-zlow);
zgl = cos ( pi * ( 0 : Nz ) / Nz )';

zgl = (1/2)*(((zupp-zlow)*zgl) + (zupp+zlow));
[X,Z,Y] = meshgrid(x,zgl,y);
A = 2*pi / Lx;
B = 2*pi / Ly;
u = (Z-zlow) .* (Z-zupp) .* sin(A*X).* sin(B*Y);
uh =fft(fft(u,[],2),[],3);

So, in MATLAB the input looks like:

val(:,:,1) =

     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0

val(:,:,2) =

         0         0         0         0         0         0         0         0
         0   -0.8071   -1.1414   -0.8071   -0.0000    0.8071    1.1414    0.8071
         0   -2.7557   -3.8971   -2.7557   -0.0000    2.7557    3.8971    2.7557
         0   -4.7042   -6.6528   -4.7042   -0.0000    4.7042    6.6528    4.7042
         0   -5.5114   -7.7942   -5.5114   -0.0000    5.5114    7.7942    5.5114
         0   -4.7042   -6.6528   -4.7042   -0.0000    4.7042    6.6528    4.7042
         0   -2.7557   -3.8971   -2.7557   -0.0000    2.7557    3.8971    2.7557
         0   -0.8071   -1.1414   -0.8071   -0.0000    0.8071    1.1414    0.8071
         0         0         0         0         0         0         0         0

Which matches my C++ input using Eigen tensor.
Now this is MATLAB's output :

val(:,:,1) =

   1.0e-13 *

   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.0016 + 0.0242i  -0.0000 + 0.0000i  -0.0016 - 0.0044i  -0.0000 + 0.0000i  -0.0016 + 0.0044i  -0.0000 - 0.0000i   0.0016 - 0.0242i
  -0.0000 + 0.0000i   0.0000 + 0.0931i  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0931i
  -0.0000 + 0.0000i  -0.0251 + 0.1087i  -0.0000 + 0.0000i   0.0251 - 0.0178i  -0.0000 + 0.0000i   0.0251 + 0.0178i  -0.0000 - 0.0000i  -0.0251 - 0.1087i
  -0.0000 + 0.0000i   0.0000 + 0.1507i  -0.0000 + 0.0000i   0.0000 + 0.0178i  -0.0000 + 0.0000i   0.0000 - 0.0178i  -0.0000 - 0.0000i   0.0000 - 0.1507i
  -0.0000 + 0.0000i  -0.0251 + 0.1087i  -0.0000 + 0.0000i   0.0251 - 0.0178i  -0.0000 + 0.0000i   0.0251 + 0.0178i  -0.0000 - 0.0000i  -0.0251 - 0.1087i
  -0.0000 + 0.0000i   0.0000 + 0.0753i  -0.0000 + 0.0000i   0.0000 + 0.0089i  -0.0000 + 0.0000i   0.0000 - 0.0089i  -0.0000 - 0.0000i   0.0000 - 0.0753i
  -0.0000 + 0.0000i   0.0016 + 0.0242i  -0.0000 + 0.0000i  -0.0016 - 0.0044i  -0.0000 + 0.0000i  -0.0016 + 0.0044i  -0.0000 - 0.0000i   0.0016 - 0.0242i
   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


val(:,:,2) =

   1.0e+02 *

   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.1582 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.1582 + 0.0000i
   0.0000 + 0.0000i   0.5400 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.5400 + 0.0000i
   0.0000 - 0.0000i   0.9218 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.9218 + 0.0000i
   0.0000 + 0.0000i   1.0800 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -1.0800 - 0.0000i
   0.0000 - 0.0000i   0.9218 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.9218 + 0.0000i
   0.0000 + 0.0000i   0.5400 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.5400 - 0.0000i
   0.0000 - 0.0000i   0.1582 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.1582 + 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

I am getting a 2D FFT to work on 3D Eigen and able to recover the inverse FFT correctly, however the forward FFT itself does not match MATLAB so I am not sure what does the 2D FFT of 3D Eigen is supposed to represent :/

Also, I only included the first two matrices of input/output as it will be too much code for this post. But, technically I have ny number of matrices inside the tensor.
Last edited on
JamieAl (236)
Sorry I thought I responded to this!
In other words, I think you're considering the indices (z, x, y) to be row z, column x, and "sheet of paper"/matrix y. Am I right?

Yes, I think this is correct! we have nz+1 rows, nx columns, ny matrices or number of "sheets".
JamieAl (236)
@mbozzi
The code you provided sort of works. It gives me the same issue I had with my code. So, when passing a full tensor, taking FFT and IFFT will return the correct first matrix of tensor but sort of missing elements in the next matrices. That tells me something wrong with passing the data as a tensor. For example, from your code:
1
2
3
4
5
6
7
 
std::cout << std::setprecision(3) << std::fixed; 
std::cout << uFun.slice(std::array{0, 0, 1}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) << "\n\n"; 
fftw_execute(fft_through_each_zx_slice);            // uFun --> uFunk
fftw_execute(inverse_fft_through_each_zx_slice);    // uFunk --> out
out = (1.0 / (nx * ny)) * out; // normalize
std::cout << out.slice(std::array{0, 0, 1}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) << '\n';

This returns:

0.000  0.000  0.000  0.000 -0.000 -0.000 -0.000 -0.000
-0.000 -0.807 -1.141 -0.807  0.000  0.807  1.141  0.807
-0.000 -2.756 -3.897 -2.756  0.000  2.756  3.897  2.756
-0.000 -4.704 -6.653 -4.704  0.000  4.704  6.653  4.704
-0.000 -5.511 -7.794 -5.511  0.000  5.511  7.794  5.511
-0.000 -4.704 -6.653 -4.704  0.000  4.704  6.653  4.704
-0.000 -2.756 -3.897 -2.756  0.000  2.756  3.897  2.756
-0.000 -0.807 -1.141 -0.807  0.000  0.807  1.141  0.807
-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.807 -1.141 -0.807  0.000  0.807  1.141  0.807
 0.000 -2.756 -3.897 -2.756  0.000  2.756  3.897  2.756
-0.000 -4.704 -6.653  0.000  0.000  4.704  6.653  4.704
-0.000 -5.511 -7.794  0.000  0.000  5.511  7.794  5.511
-0.000 -4.704 -6.653  0.000  0.000  4.704  6.653  4.704
 0.000 -2.756 -3.897  0.000  0.000  2.756  3.897  2.756
 0.000 -0.807 -1.141  0.000  0.000  0.807  1.141  0.807
 0.000  0.000 -0.000  0.000  0.000 -0.000  0.000  0.000

Notice the output isn't matching like starting from fourth column and 4th row. This only gets worse as I print other matrices. I am not sure why is this happening? As I recall this wasn't an issue if I pass an extracted 2D matrix out of 3D tensor instead.
mbozzi (3933)
You said that val(:,:,2) was one of the input blocks, and we programmed the computer to take the FFT of it. But that's a view through dimensions #1 and #2, while dimension #3 is constant. I thought we were trying to take the FFT through dimensions #2 and #3.

If I'm right, the slice should hold z constant while allowing dimensions x and y to vary. In other words, we should be taking the FFT of val(2,:,:) and not val(:,:,2). Because the axes are supposed to be (z, x, y), right?

This would explain the weird output that didn't match MATLAB from before.
Last edited on
JamieAl (236)
@mbozzi
Okay, sorry I am definitely confusing you and myself here. So my MATLAB code uses built-in meshgrid(X,Y,Z) function that rerturns a meshgrid of the size Y-by-X-by-Z meaning when taking FFTs dimension #1 refers to Y, #2 refers to X, and #3 refers to Z. However, I changed my meshgrid function inputs as the following: meshgrid(X,Z,Y) which means now my dimensions are: #1 for Z, #2 for X and #3 for Y. So, the equivalent of taking the 1D FFT of input uFun along dimension #2 (i.e X) looks like the following in C++:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
	fftw_plan const the_plan = fftw_plan_many_dft_r2c(
		1,         // compute rank-one transforms.
		&nx,      // the rank-one transforms are this long in their first (only) dimension
		(nz+1),        // do one transform per row
		Tensor2D.data(),         // the input data is here
		nullptr,   // the input data is not a row-major subarray of a larger array
		(nz+1),        // within this input array the next element is offset ri elements from prior
		1,         // address of next input array is offset 1 element from start of prior
		reinterpret_cast<fftw_complex*>(Tensor2DK.data()),         // the output data is here
		nullptr,   // the output data is not a row-major subarray of a larger array
		(nz+1),         // within this output array the next element is offset ro elements from prior
		1,         // address of next output array is offset 1 element from start of prior
		FFTW_MEASURE
	);

which is (fft(u,[],2)) in MATLAB.

what you're saying makes sense but I don't know why I am so confused :/
JamieAl (236)
I was trying to take the 1D FFT along each direction separately to understand how things work. But, I keep getting the equivalent of (fft(u,[],3)) wrong I guess.
Topic archived. No new replies allowed.