/* $Log: FFT4.java,v $
 * Revision 1.2  1999/08/31  03:58:46  stuart
 * check size argument
 *
 * Revision 1.1  1999/08/31  03:48:12  stuart
 * Initial revision
 *
 */

import java.io.*;

/** Radix 4 Fast Fourier Transform. 

@author Stuart D. Gathman, adapted from a posting by Jo Desmet 

Jo Desmet wrote:
<p>
I have attached a uuencoded file containing plain C code to make a radix
4 FFT (it could also be a radix 8, i am not sure anymore because its a while
ago that i have coded it). It is translated from machine code for the ADSP2101,
and i have optimized it for the use in C. It is possible to build a C++ wrapper
around it, but the main code must stay plain C (speed). I will not explain the
functions, I will only say that there is also a Reverse FFT.
<p>
Remember: radix 4
is not the fastest!! If you are only interested in real numbers,
there exist faster routines for it!! And I am sure that even my routine can be
more optimized, there are some places in my code that i coppy an array: try to
eliminate this.
<p>
I use this code tho prove that it is possible to use 1D algoritms on 2D Images
to do image recognision.
<p>
happy programming,<br>
<a href="jdesmet@unicall.be">Jo Desmet</a>
<a href="http://uc2.unicall.be/jdesmet">Jo's Web Page</a>
*/

public class FFT4 {

  private final double[] ReplaceXData;
  private final double[] ReplaceYData;
  private final double[] Sin;
  private final int[] UnscrambledIndex;
  private final int N,N_DIV_4,LOG2Ndiv2;

  /** Create an FFT engine.
    @param n	the data size this engine will handle - must be a power of 4
   */
  public FFT4(int n) {
    int i,j;
    for (i = 4; i != n; i <<= 2)
      if (i == 0)
        throw new IllegalArgumentException("N="+n);
    N=n;
    N_DIV_4=N/4;
    LOG2Ndiv2=((int)(Math.log(N)/Math.log(2))+1)>>1; //Juist afronden en snel
    System.err.println("LOG2Ndiv2="+ LOG2Ndiv2);

    ReplaceXData=new double[N];
    ReplaceYData=new double[N];
    Sin=new double[N+N/4];
    UnscrambledIndex=new int[N];

    for (i=0;i<N+N_DIV_4;i++)
      Sin[i]=Math.sin(i*(2*Math.PI/N));
    for (i=0;i<N;i++) {
      UnscrambledIndex[i]=0;
      for (j=0;j<LOG2Ndiv2;j++)
	UnscrambledIndex[i]|=((i>>((LOG2Ndiv2-1-j)<<1))&0x3)<<(j<<1);
    }
  }

  public void FFT(double[] Input,double[] OutputRe,double[] OutputIm) {
    int i;

    System.arraycopy(Input,0,ReplaceXData,0,N);
    for (i=0;i<N;i++) ReplaceYData[i]=0;

    butterfly();

    /* Unscramble */
    for (i=0;i<N;i++) OutputRe[i]=ReplaceXData[UnscrambledIndex[i]];
    for (i=0;i<N;i++) OutputIm[i]=ReplaceYData[UnscrambledIndex[i]];
  }   /* FFT */

  public void IFFT(double[] InputRe, double[] InputIm, double[] Output) {
    int i;

    // FFT(complex toegevoegde)=N*IFFT
    System.arraycopy(InputRe,0,ReplaceXData,0,N);
    System.arraycopy(InputIm,0,ReplaceYData,0,N);

    butterfly();
    /* Unscramble */
    for (i=0;i<N;i++) Output[N-i]=ReplaceXData[UnscrambledIndex[i]]/N;
  }   /* FFT */

  private void butterfly() {
    int BFlys=N_DIV_4;
    int BFlysBy3=(N*3)/4;
    int Groups=1,GroupsBy2=2,GroupsBy3=3;
    for (int Stage=LOG2Ndiv2;Stage > 0;Stage--) {
      int XA=0;
      int XB=XA+BFlys;
      int XC=XB+BFlys;
      int XD=XC+BFlys;
      int YA=0;
      int YB=YA+BFlys;
      int YC=YB+BFlys;
      int YD=YC+BFlys;
      for (int Group=Groups;Group > 0;Group--) {
	int Cb, Cc, Cd;
	int Sb, Sc, Sd;
	Cb=Cc=Cd=N_DIV_4;
	Sb=Sc=Sd=0;
	for (int BFly=BFlys;BFly > 0;BFly--) {
	  double xa=ReplaceXData[XA];
	  double xb=ReplaceXData[XB];
	  double xc=ReplaceXData[XC];
	  double xd=ReplaceXData[XD];
	  double ya=ReplaceYData[YA];
	  double yb=ReplaceYData[YB];
	  double yc=ReplaceYData[YC];
	  double yd=ReplaceYData[YD];

	  ReplaceXData[XA++]=xa+xb+xc+xd;
	  ReplaceYData[YA++]=ya+yb+yc+yd;
	  ReplaceXData[XB++]=(xa+yb-xc-yd)*Sin[Cb]+(ya-xb-yc+xd)*Sin[Sb];
	  ReplaceYData[YB++]=(ya-xb-yc+xd)*Sin[Cb]-(xa+yb-xc-yd)*Sin[Sb];
	  ReplaceXData[XC++]=(xa-xb+xc-xd)*Sin[Cc]+(ya-yb+yc-yd)*Sin[Sc];
	  ReplaceYData[YC++]=(ya-yb+yc-yd)*Sin[Cc]-(xa-xb+xc-xd)*Sin[Sc];
	  ReplaceXData[XD++]=(xa-yb-xc+yd)*Sin[Cd]+(ya+xb-yc-xd)*Sin[Sd];
	  ReplaceYData[YD++]=(ya+xb-yc-xd)*Sin[Cd]-(xa-yb-xc+yd)*Sin[Sd];

	  Cb+=Groups;
	  Cc+=GroupsBy2;
	  Cd+=GroupsBy3;

	  Sb+=Groups;
	  Sc+=GroupsBy2;
	  Sd+=GroupsBy3;
	}  /* ButterFly */

	XA+=BFlysBy3;
	XB+=BFlysBy3;
	XC+=BFlysBy3;
	XD+=BFlysBy3;

	YA+=BFlysBy3;
	YB+=BFlysBy3;
	YC+=BFlysBy3;
	YD+=BFlysBy3;

      }  /* Group */

      Groups<<=2;
      GroupsBy2<<=2;
      GroupsBy3<<=2;
      BFlys>>=2;
      BFlysBy3>>=2;
    }  /* Stage */
  }

  public static void main(String[] argv) throws IOException {
    StreamTokenizer data = new StreamTokenizer(new FileReader(argv[0]));
    data.resetSyntax();
    data.parseNumbers();
    data.whitespaceChars(0,' ');
    data.whitespaceChars(',',',');
    double[] v = new double[256];
    double[] ore = new double[256];
    double[] oim = new double[256];
    int vi = 0;
    for (;;) {
      int tok = data.nextToken();
      if (tok == data.TT_EOF) break;
      if (tok == data.TT_NUMBER)
        v[vi++] = data.nval;
    }
    FFT4 fft = new FFT4(vi);
    fft.FFT(v,ore,oim);
    for (int i = 0; i < vi; ++i)
      System.out.println("(" + ore[i] + ',' + oim[i] + ')');
    fft.IFFT(ore,oim,v);
    for (int i = 0; i < vi; ++i)
      System.out.println(v[i]);
  }

}