// This code computes Discrete Fourier Transform and Inverse Discrete
// Fourier Transform by two ways:
//  * non-recursively by performing a double-loop
//  * recursively by FFT algorithm
// It generates a random array of complex numbers of size 2^n where
// n is the first command line parameter, applies to it DFT followed by an
// application of IDFT to the result.
//
// 1. Save this file as FFT.java
// 2. Compile it with javac FFT.java
// 3. Run it with java FFT 10 (or n of your choice in the range 4..15)
// For n<5 the program outputs the original complex numbers and the ones
// obtained by applying DFT and IDFT. Match of the numbers indicates that
// the program works correctly.


import java.util.*;

public class FFT
{
  private static class Complex  // data structure for storing a complex number
  {
    double re;                  // real part
    double im;                  // imaginary part
    public Complex(double a, double b)
    {
      re = a;
      im = b;
    }
    public Complex(Complex c)   // copy constructor
    {
      re = c.re;
      im = c.im;
    }  
  }

  public static void main(String[] args)
  {
    // generating a random array of complex numbers of the length given as
    // the power of 2 of the command line parameter
    int n = Integer.parseInt(args[0]);
    int method = Integer.parseInt(args[1]);
    if (n >= 15)
    {
      System.out.println("Too big number");
      System.exit(0);
    }
    int n2 = 1;
    for (int i=0; i<n; i++)              // computing n2=2^n
      n2 *= 2;
    System.out.println("Generating a random array of size " + n2);
    Complex[] poly = randomComplexArray(n2, 1000);
    System.out.println("Array is ready");

    // NON-RECURSIVE DFT
    if (method == 0)
    {
      System.out.print("\nApplying non-recursive DFT ... ");
      Complex[] fourierCoeff1 = DFT1(poly);
      System.out.print("IDFT ...");
      Complex[] poly1 = IDFT1(fourierCoeff1);
      if (n<=4)
      {
        System.out.println("\nOriginal #s\t\tDFT/IDFT #s");
        for (int i=0; i<poly.length; i++)
          System.out.println(poly[i].re + "\t" + poly[i].im + "\t\t" + 
                           poly1[i].re + "\t" + poly1[i].im);
      }
      System.out.println("done");
    }

    // RECURSIVE FFT
    if (method == 1)
    {
      System.out.print("\nApplying FFT ... ");
      Complex[] fourierCoeff2 = DFT2(poly);
      System.out.print("IFFT ...");
      Complex[] poly2 = IDFT2(fourierCoeff2);
      for (int i=0; i<poly2.length; i++)    // normalization
      {
        poly2[i].re /= poly2.length;
        poly2[i].im /= poly2.length;
      }
      if (n<=4)
      {
        System.out.println("\nOriginal #s\t\tDFT/IDFT #s");
        for (int i=0; i<poly.length; i++)
        System.out.println(poly[i].re + "\t" + poly[i].im + "\t\t" +
                           poly2[i].re + "\t" + poly2[i].im);
      }
      System.out.println("done");
    }
  }

  // this function computes a nonrecursive Discrete Fourier Transform 
  public static Complex[] DFT1(Complex[] poly)
  {
    int m = poly.length;
    Complex[] cc = new Complex[m];
    double t = (2*Math.PI)/m;
    double t2 = 0;
    for (int k=0; k<m; k++)
    { 
      double a = 0;     
      double b = 0;
      double theta = 0;
      for (int j=0; j<m; j++)
      {
        double cos = Math.cos(theta);
        double sin = Math.sin(theta);
        a += (poly[j].re)*cos + (poly[j].im)*sin;
        b += (poly[j].im)*cos - (poly[j].re)*sin;
        theta += t2;
      }
      cc[k] = new Complex(a, b);
      t2 += t;
    }
    return(cc);
  }

  // this function computer a nonrecursive Inverse Discrete Fourier Transform
  public static Complex[] IDFT1(Complex[] fd)
  {
    int m = fd.length;
    Complex[] cc = new Complex[m];
    double t = (2*Math.PI)/m;
    double t2 = 0;
    for (int j=0; j<m; j++)
    {
      double a = 0;
      double b = 0;
      double theta = 0;
      for (int k=0; k<m; k++)
      {
        double cos = Math.cos(theta);
        double sin = Math.sin(theta);
	a += (fd[k].re)*cos - (fd[k].im)*sin;
	b += (fd[k].im)*cos + (fd[k].re)*sin;
        theta += t2;
      }
      cc[j] = new Complex(a/m, b/m);
      t2 += t;
    }
    return(cc);
  }

  // this function computes recursively discrete Fourier transform of poly
  // IMPORTANT: the number of nodes in poly must be a power of 2
  public static Complex[] DFT2(Complex[] poly)
  {
    int m = poly.length;
    Complex[] cc = new Complex[m];
    if (m > 1)
    {
      int i = 0, m2 = m/2;
      Complex[] a0 = new Complex[m2];
      Complex[] a1 = new Complex[m2];
      for (int k=0; k<m; k=k+2)
      {
        a0[i] = new Complex(poly[k]);
        a1[i] = new Complex(poly[k+1]);
        i++;
      }
      Complex[] y0 = DFT2(a0);
      Complex[] y1 = DFT2(a1);

      double t = (2*Math.PI)/m, theta=0;
      for (int k=0; k<m2; k++)
      {
        double cos = Math.cos(theta);
        double sin = Math.sin(theta);
        double a = (y1[k].re)*cos + (y1[k].im)*sin;
        double b = (y1[k].im)*cos - (y1[k].re)*sin;
        cc[k] = new Complex(y0[k].re + a, y0[k].im + b);
        cc[k+m2] = new Complex(y0[k].re - a, y0[k].im - b);
        theta += t;
      }
    }
    else
      cc[0] = new Complex(poly[0]);
    return(cc);
  }

  // this function computes recursively inverse Discrete Fourier Transform 
  // IMPORTANT: the number of entries in fd must be a power of 2
  // The function returns a non-normalized result. 
  // Normalization (division of the computed values by the length of the
  // array) should be done in the calling program
  public static Complex[] IDFT2(Complex[] fd)
  {
    int m = fd.length;
    Complex[] cc = new Complex[m];
    if (m > 1)
    {
      int i = 0, m2 = m/2;
      Complex[] y0 = new Complex[m2];
      Complex[] y1 = new Complex[m2];
      for (int k=0; k<m; k=k+2)
      {
        y0[i] = new Complex(fd[k]);
        y1[i] = new Complex(fd[k+1]);
        i++;
      }

      Complex[] a0 = IDFT2(y0);
      Complex[] a1 = IDFT2(y1);
 
      double t = (2*Math.PI)/m, theta=0;
      for (int k=0; k<m2; k++)
      {
        double cos = Math.cos(theta);
        double sin = Math.sin(theta);
        double a = (a1[k].re)*cos - (a1[k].im)*sin;
        double b = (a1[k].im)*cos + (a1[k].re)*sin;
        cc[k] = new Complex(a0[k].re + a, a0[k].im + b);
        cc[k+m2] = new Complex(a0[k].re - a, a0[k].im - b);
        theta += t;
      }
    }
    else
      cc[0] = new Complex(fd[0]);
    return(cc);
  }

  // this method generates a random array of complex numbers of specified
  // length. the real and imaginary parts of the generated complex numbers 
  // are integers.
  public static Complex[] randomComplexArray(int length, int range)
  {
    Complex[] a = new Complex[length];
    Random generator = new Random();
    for (int i=0; i<a.length; i++)
      a[i] = new Complex(generator.nextInt(range), generator.nextInt(range));
    return a;
  }
}