How to take 2D FFT of 3D matrix?

JamieAl (236)
I am using eigen tensor to define my 3D matrix and would like to perform a 2d fft on it. Originally, I was trying to follow with my other code where I take 1d fft of 2D eigen matrix, but I am not sure how that would work with a tensor.

My tensor is defined as the following:
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
 
static const int nx = 4;
static const int ny = 4; 
static const int nz = 4;

double Lx = 2*EIGEN_PI; 
double Ly = 2*EIGEN_PI; 
double A = (2 * EIGEN_PI)/Lx;
double B = (2 * EIGEN_PI)/ Ly;

//create X,Y,Z meshes 
Eigen::Tensor<double, 3> eXX((nz+1),nx,ny);
eXX.setZero(); 
Eigen::Tensor<double, 3> eYY((nz+1),nx,ny); 
eYY.setZero(); 
Eigen::Tensor<double, 3> eZZ((nz+1),nx,ny); 
eZZ.setZero(); 
for(int i = 0; i< nx; i++){
	for(int j = 0; j< ny; j++){
		for(int k = 0; k< nz+1; k++){ 
			eXX(k,i,j) = i*dx;
			eYY(k,i,j) = j*dy;
			eZZ(k,i,j) = -1. * cos(((k) * EIGEN_PI )/nz); 
		}

	}		
}

Eigen::Tensor<double, 3> uFun((nz+1),nx,ny); 
uFun.setZero(); 
	
for(int i = 0; i< nx; i++){
   for(int j = 0; j< ny; j++){
	for(int k = 0; k< nz+1; k++){ 
		uFun(k,i,j) = pow(eZZ(k,i,j),2.0) * sin(A * eXX(k,i,j)) * sin(B * eYY(k,i,j));
			
	}

   }		
}

In MATLAB this would be as the following: 
 
 
uh =fft(fft(uFun,[],2),[],3);

I was trying to copy using memcpy uFun to an array and then take the 1d fft of that array twice, but that didn't really work. Similarly to what I do when taking 1D FFT of 2D eigen matrix.
helios (17574)
I believe you could apply the 2D FFT three times. First on dimensions 0 and 1, then on dimensions 1 and 2, then on dimensions 0 and 3.
JamieAl (236)
@helios
why three times?
mbozzi (3932)
You can still use the fftw_plan_many_dft_r2c family of functions like here
https://www.fftw.org/fftw3_doc/Advanced-Real_002ddata-DFTs.html

fftw_plan_many_dft_r2c is fairly hard to call but it is flexible enough to set up the whole operation in one step.

I didn't test it this time but the call should be something like
1
2
 
int constexpr n[]{ ny, nx };
fftw_plan_many_dft_r2c(2, n, nz, in, n, 1, ny*nx, out, n, 1, ny*nx, flags);

That should do nz independent 2d transforms of size ny*nx.
When you execute the returned fftw_plan of course.
Last edited on
helios (17574)
On second thought, I think applying the 2D transform three times like that would effectively transform each dimension to the frequency domain's frequency domain. You should apply the 2D transform on any two dimensions, and then the corresponding 1D transform on the remaining orthogonal dimension. Or just do what mbozzi suggests.
mbozzi (3932)
I think you need to leave off the 1D transform at the end, because that final operation would give you the 3D FFT of the whole tensor.
Last edited on
helios (17574)
Oh, I misunderstood. I thought that was what OP wanted.
JamieAl (236)
Thanks everyone!
@mbozzi
I guess my main issue is passing the Eigen tensor. In the past I used to memcpy my Eigen matrix to a C++ array and perform my FFT on that array. However, doing the same thing with a tensor or a 3D matrix is making things difficult. For example this code returns strange output:
1
2
3
4
5
6
7
8
9
10
11
12
13
 
double *Mytest;
Mytest = (double*) fftw_malloc(nx*ny*nz*sizeof(double));
memset(Mytest, 42, nx*ny*nz* sizeof(double)); 

memcpy(Mytest,uFun.data(),uFun.size() * sizeof(double));

for(int i = 0; i< nx; i++){
	for(int j = 0; j< ny; j++){
		for(int k = 0; k< nz+1; k++){ 
			std::cout<<Mytest[k,i,j]<<endl;
		}
	}
}

The output here "Mytest" isn't the same as uFun. I am trying to pass this to the 2D FFT plan to test it.
Last edited on
JamieAl (236)
This is what my failed attempt look like:

Attempt1:
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
 
double A = (2 * EIGEN_PI)/Lx;
double B = (2 * EIGEN_PI)/ Ly;

Eigen::Tensor<double, 3> uFun((nz+1),nx,ny); 
uFun.setZero(); 
	
for(int i = 0; i< nx; i++){
	for(int j = 0; j< ny; j++){
		for(int k = 0; k< nz+1; k++){ 
			uFun(k,i,j) = pow(eZZ(k,i,j),2.0) * sin(A * eXX(k,i,j)) * sin(B * eYY(k,i,j));
		}

	}		
}
double *Mytest;
Mytest = (double*) fftw_malloc(nx*ny*nz*sizeof(double));
memset(Mytest, 42, nx * ny * nz * sizeof(double));

memcpy(Mytest,uFun.data(),uFun.size() * sizeof(double));
fftw_complex *dummyO; 

dummyO = (fftw_complex*) fftw_malloc((((nx)*(nz+1)))* sizeof(fftw_complex));  
memset(dummyO, 42, (((nx)*(nz+1)))* sizeof(fftw_complex)); 
	
int constexpr n[]{ ny, nx };
fftw_execute(fftw_plan_many_dft_r2c(2, n, nz, Mytest, n, 1, ny*nx, dummyO, n, 1, ny*nx, FFTW_ESTIMATE));

This returns an error: 
1
2
 
 malloc(): corrupted top size
Aborted


JamieAl (236)
Well, few cases to work, but I have questions if that's okay. The first case is a full 3D FFT of a 3D matrix that works.

Case 1 (3D FFT):
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
 
#define IMAG 1
#define REAL 0

Eigen::Tensor<double, 3> Aa((nz),nx,ny); 
Aa.setZero();
    for(int i=0; i < nx; ++i) {

        for(int j=0; j < ny; ++j) {

            for(int k=0; k < nz; ++k) {
                Aa(i, j, k) = 1.0; 

        }
    }
}
fftw_complex *input_array;
fftw_complex *output_array;

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

for (int i = 0; i < nx; ++i) {

        for (int j = 0; j < ny; ++j) {

            for (int k = 0; k < nz; ++k) {
                {
                    input_array[k + nz * (j + ny * i)][REAL] = Aa(i, j, k);
                    input_array[k + nz * (j + ny * i)][IMAG] = 0;
                }
            }
        }
    }
fftw_plan forward = fftw_plan_dft_3d(nx, ny, nz, input_array, output_array, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(forward);
fftw_destroy_plan(forward);
fftw_cleanup();
//print output
for(int i=0; i < nx; ++i){

        for(int j=0; j < ny; ++j) {

            for(int k=0; k < nz; ++k) {
                if (output_array[k + nz * (j + ny * i)][IMAG] < 0) {
                    cout << output_array[k + nz * (j + ny * i)][REAL] <<
                    " - " << abs(output_array[k + nz * (j + ny * i)][IMAG])
                    << "i" << endl;
                }
                else{
                    cout << output_array[k + nz * (j + ny * i)][REAL] << " + "
                    << abs(output_array[k + nz * (j + ny * i)][IMAG]) << "i"
                    << endl;
          }
      }
   }
}

Above example gives similar output to MATLAB: fft(fft(fft(Aa,[],1),[],2),[],3);

Case 2 (2D FFT of 3D matrix):
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
 
#define IMAG 1
#define REAL 0
	
double *Mytest;
Mytest = (double*) fftw_malloc(nx*ny*nz *sizeof(double));
memset(Mytest, 42, nx*ny*nz * sizeof(double)); 
	
	for(int i = 0; i< nx; i++){
		for(int j = 0; j< ny; j++){
			for(int k = 0; k< nz+1; k++){ 
				Mytest[k,i,j] = uFun(k,i,j); 
				
				
		}

	}		
}
fftw_complex *output_array1;
output_array1 = (fftw_complex*) fftw_malloc(nx*ny*nz * sizeof(fftw_complex));

int constexpr n[]{ ny, nx };
fftw_execute(fftw_plan_many_dft_r2c(2, n, nz, Mytest, n, 1, ny*nx, output_array1, n, 1, ny*nx, FFTW_ESTIMATE));
fftw_destroy_plan(fftw_plan_many_dft_r2c(2, n, nz, Mytest, n, 1, ny*nx, output_array1, n, 1, ny*nx, FFTW_ESTIMATE));
fftw_cleanup();
//output
for(int i=0; i < nx; ++i){

        for(int j=0; j < ny; ++j) {

            for(int k=0; k < nz; ++k) {
                if (output_array1[k + nz * (j + ny * i)][IMAG] < 0) {
                    cout << output_array1[k + nz * (j + ny * i)][REAL] <<
                    " - " << abs(output_array1[k + nz * (j + ny * i)][IMAG])
                    << "i" << endl;
                }
                else{
                    cout << output_array1[k + nz * (j + ny * i)][REAL] << " + "
                    << abs(output_array1[k + nz * (j + ny * i)][IMAG]) << "i"
                    << endl;
             }
         }
    }
}

The second case I am not sure of really. I am trying to compare it to the MATLAB case: fft(fft(uFun,[],2),[],3); . Is there a better way to print this as a 3D matrix similarly to a MATLAB output?
Also, I found out the error: 
1
2
3
 
 
malloc(): corrupted top size
Aborted

is coming from the memcpy line which I am not sure why:
 
 
memcpy(Mytest,uFun.data(),uFun.size() * sizeof(double));

Thanks again everyone!
lastchance (6980)
I doubt that you can confirm that your first code is working when all you are doing is taking the Fourier transform of a constant. That has zero frequency content.
JamieAl (236)
@lastchance
I doubt that you can confirm that your first code is working when all you are doing is taking the Fourier transform of a constant. That has zero frequency content.

You're right! I am going to test a sine function as well and report back. Thanks!
JamieAl (236)
So, this works for 3D FFT for a sine test function and I will be double checking the second example.

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
 
Eigen::Tensor<double, 3> Aa((nz+1),nx,ny); 
Aa.setZero();

	
 for(int i=0; i < nx; ++i) {

        for(int j=0; j < ny; ++j) {

            for(int k=0; k < nz+1; ++k) {
                Aa(k,i,j) = sin(A * eXX(k,i,j));

         }
     }
 }
	
	
fftw_complex *input_array;
fftw_complex *output_array;

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


for (int i = 0; i < nx; ++i) {

     for (int j = 0; j < ny; ++j) {

            for (int k = 0; k < nz+1; ++k) {
                {
                    input_array[k + (nz+1) * (j + ny * i)][REAL] = Aa(k,i,j); 
                    input_array[k + (nz+1) * (j + ny * i)][IMAG] = 0;
            }
        }
    }
}
	
fftw_plan forward = fftw_plan_dft_3d(nx, ny, (nz+1), input_array, output_array, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(forward);
fftw_destroy_plan(forward);
fftw_cleanup();
	
for(int i=0; i < nx; ++i){

    for(int j=0; j < ny; ++j) {

        for(int k=0; k < nz+1; ++k) {
                if (output_array[k + (nz+1) * (j + ny * i)][IMAG] < 0) {
                    cout << output_array[k + (nz+1) * (j + ny * i)][REAL] <<
                    " - " << abs(output_array[k + (nz+1) * (j + ny * i)][IMAG])
                    << "i" << endl;
                }
                else{
                    cout << output_array[k + (nz+1) * (j + ny * i)][REAL] << " + "
                    << abs(output_array[k + (nz+1) * (j + ny * i)][IMAG]) << "i"
                    << endl;
          }
       }
    }
}
Topic archived. No new replies allowed.