title: fftw: Why is my 2D DFT output not the same as the output from taking 1D DFT of each row?
body: <p>For my project, I have to take the DFT of a large 2D input matrix, process it, then convert it back with an IDFT and compare the results to the input matrix. My problem is in the 2D DFT step. I wrote a test with a small simple data set, which I execute in <code>main()</code>. I use the <a href="https://eigen.tuxfamily.org/dox/index.html" rel="nofollow noreferrer">Eigen</a> library for matrices and vectors. The output is this: </p>
  
  <pre><code>Using testM in a myTransform object:
  1 2 3
  4 5 6
  7 8 9
  calculateDFT took 33 Microseconds
  DFT
            (6,0)
  (-1.5,0.866025)
  (-1.5,0.866025)
           (15,0)
  (-1.5,0.866025)
  (-1.5,0.866025)
           (24,0)
  (-1.5,0.866025)
  (-1.5,0.866025)
  IDFT
  1 2 3 4 5 6 7 8 9
  Using testM in a myTransform2D object:
  1 2 3
  4 5 6
  7 8 9
  Default myTransform object created
  DFT2D
            (45,0)  (-4.5,-2.59808)          (45,-0)
            (45,0) (-13.5,-7.79423)          (45,-0)
            (45,0)           (0,-0)          (45,-0)
  IDFT
  27.5 -0.5 -1.5 -8.5  8.5  7.5   -7   10    9
  </code></pre>
  
  <p>In the snippets below, <code>this-&gt;N = this-&gt;nRows * this-&gt;nCols</code>. The results of Test 1 and Test 2 should be the same, but they are obvisouly different. I have read the documentation over and over and still can't find why it's going wrong. fftw does row-major multi-dimensional transforms, <code>in</code> is filled per row of the matrix. The <code>transfer_output</code> function does not do anything with the values themselves, only converts the standard array to an <code>Eigen::Matrix</code>. Where am I going wrong? Any help would be very much appreciated. I have also tried to find similar posts here but none had my problem as far as I could find.</p>
  
  <pre><code>void test()
  {
      RowVectorXf test(9);
      test &lt;&lt; 1, 2, 3, 4, 5, 6, 7, 8, 9;
      // Prep matrix 
      Map&lt;Matrix&lt;float, 3, 3, RowMajor&gt;&gt; testM(test.data()); // convert test to rowmajor matrix
  
      // Test 1: feed the matrix to a myTransform object and take 1D DFTs and 1D IDFTs 
      std::cout &lt;&lt; "Using testM in a myTransform object:\n" &lt;&lt; testM &lt;&lt; std::endl;
      myTransform testX1D(testM, 0);
      testX1D.vectorise();
      testX1D.calculateDFT();
      testX1D.calculateIDFT();
      std::cout &lt;&lt; "DFT" &lt;&lt; std::endl &lt;&lt; testX1D.dft &lt;&lt; std::endl;
      std::cout &lt;&lt; "IDFT" &lt;&lt; std::endl &lt;&lt; testX1D.idft &lt;&lt; std::endl; // works, too.
  
      .. Test 2: Feed the matrix to a myTransform2D object and take the 2D DFT and IDFT.
      std::cout &lt;&lt; "Using testM in a myTransform2D object:\n" &lt;&lt; testM &lt;&lt; std::endl;
      myTransform2D testX(testM, 0); // 2D version 
      testX.vectorise(); // stored in testX.m which will hold the same as test but in a colmajor vector.
      testX.calculateDFT(); // where it goes wrong?
      std::cout &lt;&lt; "DFT2D" &lt;&lt; std::endl &lt;&lt; testX.dft2D &lt;&lt; std::endl;
      testX.calculateIDFT();
      std::cout &lt;&lt; "IDFT" &lt;&lt; std::endl &lt;&lt; testX.idft &lt;&lt; std::endl;
  }
  </code></pre>
  
  <p>This is how I calculate the DFTs in each case, using the fftw library (fftwf because I use single precision to save memory and the values of the non-test data are on the order of -10000 to 10000, so I think that is not an issue). </p>
  
  <pre><code>void myTransform::calculateDFT()
  /// Calculates discrete fourier transform of vectorised data `m`.
  /** uses the FFTW library (https://fftw.org). The dft is stored in myTransform::dft*/
  {
      //std::cout &lt;&lt; m &lt;&lt; std::endl;
      fftwf_complex* out;
      fftwf_plan p;
      out = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * this-&gt;nCols);
      float* in = new float[static_cast&lt;const float&amp;&gt;(this-&gt;nCols)];
      p = fftwf_plan_dft_r2c_1d(this-&gt;nCols, in, out, FFTW_ESTIMATE);
      // calculate DFT for each trace and assign it to a segment of this-&gt;dft
      unsigned int factor = 0;
      auto check = std::chrono::high_resolution_clock::now();
      for (int k = 0; k &lt; this-&gt;nRows; k++)
      {
          factor = k * this-&gt;nCols;
          //TODO: if possible, fix this slow in[i] = ... part. 
          for (int i = 0; i &lt; this-&gt;nCols; i++)
          {
              in[i] = this-&gt;m[factor + i];
          }
          p = fftwf_plan_dft_r2c_1d(this-&gt;nCols, in, out, FFTW_ESTIMATE);
          fftwf_execute(p);
          this-&gt;transfer_output(out, k); // does nothing but add the output to a vector dft. 
      }
      delete [] in;
      fftwf_free(out);
      fftwf_destroy_plan(p);
  }
  </code></pre>
  
  <p>And for the 2D DFT case: Here I use std::complex as specified on <a href="http://www.fftw.org/fftw3_doc/Complex-numbers.html#Complex-numbers" rel="nofollow noreferrer">fftw.org</a>. I allocate <code>nRows * (nCols/2 + 1)</code> single-precision floats as instructed <a href="http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data" rel="nofollow noreferrer">here</a>. For the 1D case, this is done in the 1D <code>transfer_output</code> function where <code>dft</code> is filled with <code>out[this-&gt;nCols - i]</code> for <code>i &gt; this-&gt;nCols/2</code></p>
  
  <pre><code>void myTransform2D::calculateDFT()
  /// Should calculate the DFT in 2D with fftwf_plan_dft_r2c_2d(n0, n1, *in, *out, flags).
  {
      std::complex&lt;float&gt;* out;
      fftwf_plan p;
      out = (std::complex&lt;float&gt;*)fftwf_malloc(sizeof(std::complex&lt;float&gt;) * this-&gt;nRows * (this-&gt;nCols/2+1)); // Hermitian symmetry for r2c transforms
      float* in = new float[this-&gt;N];
      in = (float*)fftwf_malloc(sizeof(float) * this-&gt;N);
      p = fftwf_plan_dft_r2c_2d(this-&gt;nRows, this-&gt;nCols, in, reinterpret_cast&lt;fftwf_complex*&gt;(out), FFTW_ESTIMATE);
      // Fill input array
      for (int i = 0; i &lt; this-&gt;nRows; i++)
      {
          int factor = i * this-&gt;nCols;
          for (int j = 0; j &lt; this-&gt;nCols; j++)
          {
              in[factor + j] = this-&gt;m[factor + j];
          }
      }
      fftwf_execute(p);
      transfer_output(out);
      fftwf_free(in);
      fftwf_free(out);
      fftwf_destroy_plan(p);
  }
  </code></pre>
  
  <p>I convert back to the time domain using the IDFT, again in both 1D and 2D. I don't know for sure if the 2D version works, as the DFT goes wrong. The 1D case works, so I show the 2D case only.</p>
  
  <pre><code>void myTransform2D::calculateIDFT()
  /// Calculates inverse fourier transform of `this-&gt;dft2D`.
  /** Also uses the FFTW library. Results might not be perfect as we use floats
      instead of doubles because of large data sizes. */
  {
      float* out = new float[this-&gt;N];
      std::complex&lt;float&gt;* in; 
      fftwf_plan pr;
      in = (std::complex&lt;float&gt;*)fftwf_malloc(sizeof(std::complex&lt;float&gt;) * this-&gt;N);
      out = (float*)fftwf_malloc(sizeof(float) * this-&gt;N);
      pr = fftwf_plan_dft_c2r_2d(this-&gt;nRows, this-&gt;nCols, reinterpret_cast&lt;fftwf_complex*&gt;(in), out, FFTW_ESTIMATE);
  
      for (int i = 0; i &lt; this-&gt;nRows; i++)
      {
          for (int j = 0; j &lt; this-&gt;nCols; j++)
          {
              in[i * this-&gt;nCols + j] = this-&gt;dft2D(i, j);
          }
      }
  
      fftwf_execute(pr);
  
      for (int i = 0; i &lt; this-&gt;N; i++)
      {
          this-&gt;idft[i] = out[i] / this-&gt;N; // fftw does unnormalized inverse transforms.
      }
  
      fftwf_free(out);
      fftwf_free(in);
      fftwf_destroy_plan(pr);
  }
  </code></pre>
  
  <p>EDIT: removed some code as suggested
  EDIT2: removed image, added contents as text.</p>