title: Port Matlab's FFT to native Java
body: <p>I want to port Matlab's Fast Fourier transform function fft() to native Java code.</p>
  
  <p>As a starting point I am using the code of <a href="http://www.jmathlib.de/" rel="nofollow">JMathLib</a> where the FFT is implemented as follows:</p>
  
  <pre><code>    // given double[] x as the input signal
  
      n = x.length;  // assume n is a power of 2
      nu = (int)(Math.log(n)/Math.log(2));
      int n2 = n/2;
      int nu1 = nu - 1;
      double[] xre = new double[n];
      double[] xim = new double[n];
      double[] mag = new double[n2];
      double tr, ti, p, arg, c, s;
      for (int i = 0; i &lt; n; i++) {
          xre[i] = x[i];
          xim[i] = 0.0;
      }
      int k = 0;
  
      for (int l = 1; l &lt;= nu; l++) {
          while (k &lt; n) {
              for (int i = 1; i &lt;= n2; i++) {
                  p = bitrev (k &gt;&gt; nu1);
                  arg = 2 * (double) Math.PI * p / n;
                  c = (double) Math.cos (arg);
                  s = (double) Math.sin (arg);
                  tr = xre[k+n2]*c + xim[k+n2]*s;
                  ti = xim[k+n2]*c - xre[k+n2]*s;
                  xre[k+n2] = xre[k] - tr;
                  xim[k+n2] = xim[k] - ti;
                  xre[k] += tr;
                  xim[k] += ti;
                  k++;
              }
              k += n2;
          }
          k = 0;
          nu1--;
          n2 = n2/2;
      }
      k = 0;
      int r;
      while (k &lt; n) {
          r = bitrev (k);
          if (r &gt; k) {
              tr = xre[k];
              ti = xim[k];
              xre[k] = xre[r];
              xim[k] = xim[r];
              xre[r] = tr;
              xim[r] = ti;
          }
          k++;
      }
      // The result 
      // -&gt; real part stored in xre
      // -&gt; imaginary part stored in xim
  </code></pre>
  
  <p>Unfortunately it doesn't give me the right results when I unit test it, for example with the array</p>
  
  <blockquote>
    <p>double[] x = { 1.0d, 5.0d, 9.0d, 13.0d };</p>
  </blockquote>
  
  <p>the result in Matlab: </p>
  
  <blockquote>
    <p>28.0<br>
     -8.0 - 8.0i<br>
      -8.0<br>
     -8.0 + 8.0i</p>
  </blockquote>
  
  <p>the result in my implementation:</p>
  
  <blockquote>
    <p>28.0<br>
     -8.0 + 8.0i<br>
      -8.0<br>
     -8.0 - 8.0i</p>
  </blockquote>
  
  <p>Note how the signs are wrong in the complex part.</p>
  
  <p>When I use longer, more complex signals the differences between the implementations affects also the numbers. So the implementation differences does not only relate to some sign-"error".</p>
  
  <p>My question: how can I adapt my implemenation to make it "equal" to the Matlab one? 
  Or: is there already a library that does exactly this?</p>