/* Fast Fourier transform methods for 1-d and 2-d data, using Rocca's
     column FT routine for more efficiency.
     Adapted from Jon Claerbout's Imaging the Earth's Interior (1985)
     by John N. Louie, 24 November 1998.
*/

/* Include the usual classes */
import java.util.*;
import java.io.*;
import java.lang.Math;

/* Declare a class to contain all the 2-d fft methods we will use. */
/* When you change the code in this file, make a copy to a new file with
    a new name like "ex6fttestc.java", and be sure to make the same change
    to the class name below - e.g. "class ex6fttestc" */
class fttestc
	{
	/* The main() method compiles into a file named "fttestc.class" that
	    you run in the Java virtual machine as a stand-alone application
	    with a command such as "java fttestc"; or whatever you have
	    renamed the class and file. */
	public static void main(String args[]) /* No runtime arguments needed */
		{
		System.out.println("Starting fttestc...");
		/* Declare variables */
		int it, nt, ix, nx;
		/* Assign sizes: 64x64 is always quick; 1024x1024 takes several
		    minutes and more than 40 Mb memory */
		nx = 64;
		nt = 64;

		/* Declare and initialize planes to zero values */
		FltPlane p = new FltPlane(nx, nt);
		ComplexPlane cp = new ComplexPlane(nx, nt);
		
		/* Insert a spike into the real part of the complex array.
		    The class ComplexPlane has the method setElemReal(ix,it,val)
		    as a shortcut for the statement: cp.vec[ix*cp.ne + it].re = val; */
		cp.setElemReal(3,16,1F); cp.setElemReal(4,16,4F);
		cp.setElemReal(5,16,6F); cp.setElemReal(6,16,4F);
		cp.setElemReal(7,16,1F); cp.setElemReal(3,17,1F);
		cp.setElemReal(4,17,4F); cp.setElemReal(5,17,6F);
		cp.setElemReal(6,17,4F); cp.setElemReal(7,17,1F);

		/* Transform plane into the 2-d Fourier domain */
		ft2dc(cp, 1F, 1F);

		/* Copy the real part of the transformed complex plane to a real plane
		   to allow plotting */
		cp.realtoFltPlane(p);
		
		/* Plot the real plane in a new window titled "fttestc", with
		    the fastest-varying direction pointing down, no vertical
		    exaggeration, using white color for negative amplitudes and
		    black color for positive amplitudes with near-zero amplitudes
		    gray (like looking at a WTVA wiggle-trace-variable-area plot from
		    a distance), scaling from the maximum down to minus the maximum. */
		p.viewinFrame("fttestc", true, 1F, RgCtab.vagray, false, -1F);
		
		/* All Done with fttestc's main() application method */
		}

	/* 2-d FFT in-place method for class ComplexPlane data structures,
	    using Rocca's columnwise FFT scheme */
	static void ft2dc(ComplexPlane cp, float sign1, float sign2)
		{
		/* Declare local variables */
		int n1, n2;
		int i1, i2;
		
		/* Obtain n1=nx and n2=nt from the input ComplexPlane class
		    data structure, using RGData interface methods */
		n1 = cp.getNMembofDim(1); /* Dim=1 is vectors, or rows */
		n2 = cp.getNMembofDim(0); /* Dim=0 is fastest-varying, the elements */
		
		/* Adding String objects concatenates them */
		System.out.print("ft2dc: starting " + n2 + "x" + n1 + " 2-d FFT... ");
		
		/* Transform over the fast dimension, the rows */
		/* Loop over n1, number of rows */
		for (i1=0; i1<n1; i1++)
			/* 1-d Fourier transform of row of plane cp, done n1 times.
			    Here use a form of the fft method that does not require
			    copying the row to a scratch vector, and back into cp */
			fftr(n2, cp, i1*n2, sign1, 1F);
		
		/* Transform over the slow dimension, the columns; assumes
		    the complex planes are small enough to all fit in memory */
		
		colmft(cp, sign2, 1F);
		
		/* Standard output terminated with a newline */
		System.out.println("Done.");
		
		/* Return to calling method */
		}

	/* Fabio Rocca's "row" Fourier transform, from Claerbout (1985) p. 73-74,
	    modified by J. Louie for the indexing scheme used by C (Dec. 1987)
	    and Java (Nov. 1998). */
	static void colmft(ComplexPlane cx, float sign2, float scale)
		{
		/* Declare all local variables */
		int i, i1, i2, j, m, lstep, istep, n1, n2;
		float arg;
		
		Complex cw = new Complex();
		Complex cdel = new Complex();
		/* Obtain n1=nx and n2=nt from the input ComplexPlane class
		    data structure, using RGData interface methods */
		n1 = cx.getNMembofDim(1); /* Dim=1 is vectors, or rows */
		n2 = cx.getNMembofDim(0); /* Dim=0 is fastest-varying, the elements */
		
		/* Multiply each element in matrix by scale factor */
		for (i1=0; i1<n1; i1++)
			for (i2=0; i2<n2; i2++)
				cx.vec[i1*n2+i2].scale(scale);
		
		j = 0;
		for (i=0; i<n1; i++)
			{
			if (i <= j)
				twid1(n2, cx, i*n2+0, cx, j*n2+0);
			m = n1/2;
			while (j > m-1)
				{
				j = j - m;
				m = m/2;
				if (m < 1)
					break;
				}
			j = j + m;
			}
		lstep = 1;
		do
			{
			istep = 2*lstep;
			cw.re = 1F;
			cw.im = 0F;
			arg = (float)(sign2*Math.PI/lstep);
			cdel.re = (float)(Math.cos((double)(arg)));
			cdel.im = (float)(Math.sin((double)(arg)));
			for (m=0; m<lstep; m++)
				{
				for (i=m; i<n1; i+=istep)
					twid2(n2, cw, cx, i*n2+0, cx, (i+lstep)*n2+0);
				cw.mul(cw, cdel);
				}
			lstep = istep;
			}
			while (lstep < n1);
		/* Return to calling method */
		}
	static void twid1(int n, ComplexPlane cx, int j0, ComplexPlane cy, int k0)
		{
		int i;
		Complex ct = new Complex();
		for (i=0; i<n; i++)
			{
			ct.re = cx.vec[j0+i].re;
			ct.im = cx.vec[j0+i].im;
			cx.vec[j0+i].re = cy.vec[k0+i].re;
			cx.vec[j0+i].im = cy.vec[k0+i].im;
			cy.vec[k0+i].re = ct.re;
			cy.vec[k0+i].im = ct.im;
			}
		}
	static void twid2(int n, Complex cw, ComplexPlane cx, int j0,
		ComplexPlane cy, int k0)
		{
		int i;
		Complex ctemp = new Complex();
		for (i=0; i<n; i++)
			{
			ctemp.mul(cw, cy.vec[k0+i]);
			cy.vec[k0+i].sub(cx.vec[j0+i], ctemp);
			cx.vec[j0+i].add(cx.vec[j0+i], ctemp);
			}
		}

	/* Simple in-place 1-d FFT method for a vector within a class ComplexPlane
	    data structure. The number of complex elements in each vector
	    MUST be a power of 2 */
	static void fftr(int lx, ComplexPlane cx, int j0, float signi, float sc)
		{
		/* Declare all local variables */
		int i, j, k, m, istep;
		float arg;
		Complex cw = new Complex();
		Complex ct = new Complex();
		/* Here the lx argument supplies the vector length, since we are
		    working on one vector within the cx class ComplexPlane */
		
		j = 0;
		k = 1;
		/* Left out for compatibility with fttestc.colmft
		sc = (float)(sqrt(1.0/(double)(lx)));
		*/
		for (i=0; i<lx; i++)
			{
			if (i <= j)
				{
				ct.re = sc * cx.vec[j0+j].re;
				ct.im = sc * cx.vec[j0+j].im;
				cx.vec[j0+j].re = sc * cx.vec[j0+i].re;
				cx.vec[j0+j].im = sc * cx.vec[j0+i].im;
				cx.vec[j0+i].re = ct.re;
				cx.vec[j0+i].im = ct.im;
				}
			m = lx/2;
			while (j > m-1)
				{
				j = j - m;
				m = m/2;
				if (m < 1)
					break;
				}
			j = j + m;
			}
		do
			{
			istep = 2*k;
			for (m=0; m<k; m++)
				{
				arg = (float)(Math.PI*signi*m/k);
				cw.re = (float)(Math.cos((double)(arg)));
				cw.im = (float)(Math.sin((double)(arg)));
				for (i=m; i<lx; i+=istep)
					{
					/* class Complex multiplication method ct = cw*cx[i+k] */
					ct.mul(cw, cx.vec[j0+i+k]);
					/* class Complex subtraction method cx[i+k] = cx[i] - ct */
					cx.vec[j0+i+k].sub(cx.vec[j0+i], ct);
					/* class Complex addition method cx[i] = cx[i] + ct */
					cx.vec[j0+i].add(cx.vec[j0+i], ct);
					}
				}
			k = istep;
			}
			while (k < lx);
		/* Return to calling method */
		}

	/* End of methods in class fttestc */
	}