Mathew77
Mathew77

Reputation: 87

1D FFT along different dimensions of 3D array in MATLAB vs C++

This maybe a simple question, but unfortunately I am unable to find a simple answer to. I am taking 1D FFT of a 3D matrix (tensor) of size: # of rows by # of columns by # of matrices and so in MATLAB this translates to:

clearvars; clc; close all;
Nx = 8;  
Ny = 4; 
Nz= 6;
Lx =16;
Ly = 6;
dx = Lx/Nx; 
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); %this gives 3d grid with z-by-x-by-y size (i.e. ZXY)
%ICs
A = 2*pi / Lx;
B = 2*pi / Ly;
u = (Z-zlow) .* (Z-zupp) .* sin(A*X).* sin(B*Y);
uh1 =(fft(u,[],1));
uh2 =(fft(u,[],2));
uh3 =(fft(u,[],3));

This produces a 3D tensor with (Nz+1) rows and Nx columns and Ny matrices and taking the 1D FFT along each "column" of u is represented by uh1 =(fft(u,[],1)); and the 1D FFT along each "row" of us is represented by uh2 =(fft(u,[],2)); according to the MATLAB FFT documentation.

This works fine with me when using FFT and Eigen library in C++! However, I am not sure what

uh3 =(fft(u,[],3));

means? what is the 3rd dimension represents? I thought that meant how many matrices I had inside of the tensor u? I am unsure how to calculate this in my C++ code.

Example of C++ code:

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

int main(){
Eigen::Tensor<double, 3> uFun((nz+1),nx,ny);//tensor(row,col,matrix)
uFun.setZero();
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) = //define function u same as MATLAB
            }
        }
    }
//define matrix and copy the tensor into it
//map tensor to matrix

Eigen::MatrixXd dummy((nz+1),nx); //to match tensor slices
dummy.setZero();
Eigen::Tensor<double, 2> Tensor2D((nz+1),nx);   //tensor(row,col,matrix)
Tensor2D.setZero();

Tensor2D = uFun.slice(std::array{0, 0, 3}, std::array{nz + 1, nx, 1}).reshape(std::array{nz + 1, nx}) ;
dummy = static_cast<MatrixType<double>>(Eigen::Map<const MatrixType<double>>(Tensor2D.data(), (nz+1), nx));
Eigen::MatrixXcd tempOut1((nz+1),nx);
//uh1= fft(u,[],1);
    for (int k=0; k< dummy.cols(); k++){
        tempOut1.col(k) = fft.fwd(dummy.col(k));        
    }
//uh2= fft(u,[],2);
Eigen::MatrixXcd tempOut((nz+1),nx);
Eigen::VectorXcd row(nx);
    for (int k=0; k< (dummy.rows()); k++){
        fft.fwd(row, dummy.row(k));
        tempOut.row(k) = row;
    }

}

The above code works for uh1=fft(u,[],1) and uh2=fft(u,[],2) but not uh3=fft(u,[],3) which I am not sure what does it represent here?

Here's the MATLAB FFT documentation that I am following with: https://www.mathworks.com/help/matlab/ref/fft.html#buuutyt-5

NOTE: I know I can use FFTW library for easier use with Eigen tensors, but in the meantime this has been much more user-friendly to me "personally"

EDIT: Thanks to @CrisLuengo I have better idea of what each dimension represents here. What I am trying to do looks like the following in MATLAB:

numRows=Nz+1;numCols=Nx;
for i= 1:numRows
   for j= 1:numCols
      uk(i,j,:) = fft(u(i,j,:));
   end
end

where the above code is basically uh3 =(fft(u,[],3)); and now using EigenFFT maybe something like this (NOT TESTED yet):

for (int i = 0; i < dummy.rows(); ++i) {
        for (int j = 0; j < dummy.cols(); ++j) {
            int index =  i * dummy.cols() + j; // Flattened index
            VectorXcd slice = dummy.row(index);
            tempOut.row(index) = fft.fwd(slice);
        }
    }

Upvotes: 0

Views: 98

Answers (1)

Jamie
Jamie

Reputation: 519

This should work:

for (int i = 0; i < dummy.rows(); ++i){
        VectorXd slice = dummy.row(i);
        VectorXcd tmpOut(nx);
        fft.fwd(tmpOut, slice);
        for (int j = 0; j < dummy.cols(); ++j){
            tempOut(i, j) = tmpOut(j);
        }
    }

    std::cout<<tempOut<<endl;

Upvotes: -1

Related Questions