c++matlabeigen3fftw

taking 1D FFT of Eigen matrix along rows using Eigen FFT library returns incorrect results

I have seen few examples on Eigen FFT library that I wanted to implement before moving on FFTW library. For some reason I am unable to get the correct 1D FFT along each row of an Eigen matrix. I get the 1D FFT along each column correctly, but not rows. I am comparing a very very simple C++ code with MATLAB for clarity.

In MATLAB:

Nx=8;
Ny=8;
u = eye(Ny+1,Nx); %very simple choice for now
%take FFT along columns
ucol = fft(u,[],1);
urow = fft(u,[],2);

In C++:

#include "fftw3.h"
#include<unsupported/Eigen/CXX11/Tensor>

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

using namespace std;
using namespace Eigen;
int main{
Eigen::FFT<double> fft;
Eigen::MatrixXcd u(ny+1,nx);
u.setIdentity();

//1D FFT along cols
Eigen::MatrixXcd ucol(ny+1,nx);
for (int k=0; k<u.cols(); k++){
     ucol.col(k) = fft.fwd(u.col(k));
}
std::cout<<ucol<<endl; //correct compared to MATLAB: ucol = fft(u,[],1);

//1D FFT along rows 
Eigen::MatrixXcd urow(ny+1,nx);
for (int k=0; k<u.rows(); k++){
     urow.row(k) = fft.fwd(u.row(k));
}
std::cout<<urow<<endl; //**Incorrect** compared to MATLAB: urow = fft(u,[],2);
}

Why is the second 1D FFT incorrect compared to MATLAB??

MATLAB outputs:

ucol =

   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   0.7660 - 0.6428i   0.1736 - 0.9848i  -0.5000 - 0.8660i  -0.9397 - 0.3420i  -0.9397 + 0.3420i  -0.5000 + 0.8660i   0.1736 + 0.9848i
   1.0000 + 0.0000i   0.1736 - 0.9848i  -0.9397 - 0.3420i  -0.5000 + 0.8660i   0.7660 + 0.6428i   0.7660 - 0.6428i  -0.5000 - 0.8660i  -0.9397 + 0.3420i
   1.0000 + 0.0000i  -0.5000 - 0.8660i  -0.5000 + 0.8660i   1.0000 + 0.0000i  -0.5000 - 0.8660i  -0.5000 + 0.8660i   1.0000 + 0.0000i  -0.5000 - 0.8660i
   1.0000 + 0.0000i  -0.9397 - 0.3420i   0.7660 + 0.6428i  -0.5000 - 0.8660i   0.1736 + 0.9848i   0.1736 - 0.9848i  -0.5000 + 0.8660i   0.7660 - 0.6428i
   1.0000 + 0.0000i  -0.9397 + 0.3420i   0.7660 - 0.6428i  -0.5000 + 0.8660i   0.1736 - 0.9848i   0.1736 + 0.9848i  -0.5000 - 0.8660i   0.7660 + 0.6428i
   1.0000 + 0.0000i  -0.5000 + 0.8660i  -0.5000 - 0.8660i   1.0000 + 0.0000i  -0.5000 + 0.8660i  -0.5000 - 0.8660i   1.0000 + 0.0000i  -0.5000 + 0.8660i
   1.0000 + 0.0000i   0.1736 + 0.9848i  -0.9397 + 0.3420i  -0.5000 - 0.8660i   0.7660 - 0.6428i   0.7660 + 0.6428i  -0.5000 + 0.8660i  -0.9397 - 0.3420i
   1.0000 + 0.0000i   0.7660 + 0.6428i   0.1736 + 0.9848i  -0.5000 + 0.8660i  -0.9397 + 0.3420i  -0.9397 - 0.3420i  -0.5000 - 0.8660i   0.1736 - 0.9848i


urow =

   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i   1.0000 + 0.0000i
   1.0000 + 0.0000i   0.7071 - 0.7071i   0.0000 - 1.0000i  -0.7071 - 0.7071i  -1.0000 + 0.0000i  -0.7071 + 0.7071i   0.0000 + 1.0000i   0.7071 + 0.7071i
   1.0000 + 0.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i   1.0000 + 0.0000i   0.0000 - 1.0000i  -1.0000 + 0.0000i   0.0000 + 1.0000i
   1.0000 + 0.0000i  -0.7071 - 0.7071i   0.0000 + 1.0000i   0.7071 - 0.7071i  -1.0000 + 0.0000i   0.7071 + 0.7071i   0.0000 - 1.0000i  -0.7071 + 0.7071i
   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i  -1.0000 + 0.0000i
   1.0000 + 0.0000i  -0.7071 + 0.7071i   0.0000 - 1.0000i   0.7071 + 0.7071i  -1.0000 + 0.0000i   0.7071 - 0.7071i   0.0000 + 1.0000i  -0.7071 - 0.7071i
   1.0000 + 0.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i   1.0000 + 0.0000i   0.0000 + 1.0000i  -1.0000 + 0.0000i   0.0000 - 1.0000i
   1.0000 + 0.0000i   0.7071 + 0.7071i   0.0000 + 1.0000i  -0.7071 + 0.7071i  -1.0000 + 0.0000i  -0.7071 - 0.7071i   0.0000 - 1.0000i   0.7071 - 0.7071i
   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

C++ output:

ucol=

(1.000,0.000)   (1.000,0.000)   (1.000,0.000)   (1.000,0.000)   (1.000,0.000)   (1.000,0.000)   (1.000,0.000)   (1.000,0.000)
  (1.000,0.000)  (0.766,-0.643)  (0.174,-0.985) (-0.500,-0.866) (-0.940,-0.342)  (-0.940,0.342)  (-0.500,0.866)   (0.174,0.985)
  (1.000,0.000)  (0.174,-0.985) (-0.940,-0.342)  (-0.500,0.866)   (0.766,0.643)  (0.766,-0.643) (-0.500,-0.866)  (-0.940,0.342)
  (1.000,0.000) (-0.500,-0.866)  (-0.500,0.866)   (1.000,0.000) (-0.500,-0.866)  (-0.500,0.866)   (1.000,0.000) (-0.500,-0.866)
  (1.000,0.000) (-0.940,-0.342)   (0.766,0.643) (-0.500,-0.866)   (0.174,0.985)  (0.174,-0.985)  (-0.500,0.866)  (0.766,-0.643)
 (1.000,-0.000)  (-0.940,0.342)  (0.766,-0.643)  (-0.500,0.866)  (0.174,-0.985)   (0.174,0.985) (-0.500,-0.866)   (0.766,0.643)
 (1.000,-0.000)  (-0.500,0.866) (-0.500,-0.866)  (1.000,-0.000)  (-0.500,0.866) (-0.500,-0.866)  (1.000,-0.000)  (-0.500,0.866)
 (1.000,-0.000)   (0.174,0.985)  (-0.940,0.342) (-0.500,-0.866)  (0.766,-0.643)   (0.766,0.643)  (-0.500,0.866) (-0.940,-0.342)
 (1.000,-0.000)   (0.766,0.643)   (0.174,0.985)  (-0.500,0.866)  (-0.940,0.342) (-0.940,-0.342) (-0.500,-0.866)  (0.174,-0.985)

urow=

(1.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)
 (1.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)
 (1.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)
 (1.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)
 (1.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)
 (1.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)
 (1.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)
 (1.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)

why is urow from the C++ code not matching the MATLAB output?? I am pretty sure the second loop is iterating over all rows of u

Thanks

Solution

  • Your usage of fft.fwd is not correct. The documentation shows that the function expects the frequency vector and then the time vector. You are providing only one input.

    Here is godbolt: https://www.godbolt.org/z/Mecfr487f

    Here is just the code.

    #include <Eigen/Core>
#include<unsupported/Eigen/CXX11/Tensor>
#include <unsupported/Eigen/FFT>
#include <iostream>

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

using namespace std;
using namespace Eigen;

int main(int argc, char * argv[]){
Eigen::FFT<double> fft;
Eigen::MatrixXcd u(ny+1,nx);
u.setIdentity();

//1D FFT along cols
Eigen::MatrixXcd ucol(ny+1,nx);
Eigen::VectorXcd col(ny+1);
for (int k=0; k<u.cols(); k++){
     // ucol.col(k) = fft.fwd(u.col(k));
     fft.fwd(col, u.col(k));
     ucol.col(k) = col;
}
std::cout<<ucol<<endl;

//1D FFT along rows 
Eigen::MatrixXcd urow(ny+1,nx);
Eigen::VectorXcd row(nx);
for (int k=0; k<u.rows(); k++){
     // urow.row(k) = fft.fwd(u.row(k));
     fft.fwd(row, u.row(k));
     urow.row(k) = row;
}
std::cout<<urow<<endl;
}