Reputation: 87
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
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