/**
 * A simple textbook implementation of the 1-dimensional Fast Fourier
 * Transform.
 * 
 * ...I found this code on my hard drive... It's been so many year
 * since I touched this that I don't remember if this ever worked as
 * it should.. so beware. It's possible that this works better than
 * the other example codes that I've been hacking at a fast pace:) but
 * I didn't try this out. Looks like I gave it some thought back then,
 * though...
 * 
 * This is for your convenience if you want to try the rotation- and
 * translation invariant shape descriptor features discussed on
 * Tuesdays lecture.
 * 
 * (You'll have to implement the boundary tracking by yourself, and
 * distance-from-center computation. Something to learn by doing
 * instead of just watching ;)).
 */

/**
 * A simple textbook implementation of the Fast Fourier Transform.
 * 
 * @author nieminen
 */
public class SimpleFFT {
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		System.out.println("Could have tests here, but so far there's none.");
	}
	
	/**
	 * Bit-reverse operation. 
	 * 
	 * (reverses the lowest bits of an integer, and clears the rest.)
	 * 
	 * example: bitrev(5,''10010'') is ''01001''.
	 * example: bitrev(5,''100110010'') is ''01001''.
	 * 
	 * @param nbits  - interpret val to have this many bits
	 * @param val    - value to be bit-reversed
	 * @return       - bit-reversed value.
	 */
	private static int bitreverse(int nbits, int val) {
		int res = 0;
		for (int bit=0; bit<nbits; bit++){
		    res <<= 1;
			res += val & 1;
			val >>= 1;
		}
     	return res;
	}
	

	/**
	 * Fast Fourier transform for complex-valued vectors of size N=2^n. 
	 * 
	 * "By the book". Only for myself to finally learn how this is done.
	 * Use a mature library in real applications!
	 * 
	 * Normalizes the result by 1/sqrt(N).
	 * 
	 * All arrays must have the same size, and the size must be 2^n for 
	 * integral n. The algorithm does the FFT essentially "in place", so you 
	 * may choose to pass arrays such that reSrc==reFt && imSrc==rmFt.
	 * 
	 * @param reSrc   Real parts of input data
	 * @param imSrc   Imaginary parts of input data
	 * @param reFt    Storage array for Re(ft(x[i])). Will be rewritten here;
	 * @param imFt    Storage array for Im(ft(x[i])). Will be rewritten here;
	 */
	public static void fft(double[] reSrc, double[] imSrc, double[] reFt, double[] imFt) {
		
		final int n = checkAndGiveMaxPowerOfTwo(reSrc, imSrc, reFt, imFt);
		final int flen = 1<<n;

		final double twoPi = 2 * Math.PI;

		/* We can use the output arrays for also the intermediate results, if
		 * we use the special indexing scheme with bit-reversing.
		 */
		for(int i=0; i<flen; i++){
			int j = bitreverse(n, i); 
			if (j>i) continue; // Do each bit-reverse swap only once. 
			double reTmp = reSrc[j];
			double imTmp = imSrc[j];
			reFt[j] = reSrc[i];
			imFt[j] = imSrc[i];
			reFt[i] = reTmp;
			imFt[i] = imTmp;			
		}

		// The Cooley-Tukey FFT as found in many textbook examples:
		for (int m=1; m<=n; m++){
			final int twoToM = 1 << m;
			final int twoToMMinus1 = 1 << (m-1);
			
			for(int k=0; k<twoToMMinus1; k++){
				//Twiddle by "exp(i2Pi/N)^(-kN/2^m)" == "exp(-ik2Pi/2^m)"
				final double angle = -k*twoPi/(double)twoToM;
				final double realpha = Math.cos(angle);
				final double imalpha = Math.sin(angle);
				
				final int jmax = 1 << (n-m); // j=0, ..., 2^(n-m)-1
				for (int j=0; j<jmax; j++){
					int uind = k+twoToM*j;
					int vind = twoToMMinus1+uind;
					
					double reu = reFt[uind];
					double imu = imFt[uind];
					
					// Remember complex multiply rule: (a+ib)(c+id)=ac-bd+i(ad+bc)
					double rev = realpha * reFt[vind] - imalpha * imFt[vind];
					double imv = realpha * imFt[vind] + imalpha * reFt[vind];
	
					// Then add.
					reFt[uind] = reu + rev;
					imFt[uind] = imu + imv;
					
					reFt[vind] = reu - rev;
					imFt[vind] = imu - imv;
				}
			}
		}
		
		final double sqrtn = Math.sqrt(flen);
		for(int k=0; k<flen; k++){
			reFt[k] /= sqrtn;
			imFt[k] /= sqrtn;
		}
	}

	
	/**
	 * Returns the biggest n for which 2^n <= length of array tables. 
	 * 
	 * Also checks the array sizes: Warns to std err if array size is not 
	 * a power of two; throws an {@link IllegalArgumentException} if not
	 * all arrays are of the same size.
	 *   
	 * (A private helper method.)
	 * 
	 * @param reSrc
	 * @param imSrc
	 * @param reFt
	 * @param imFt
	 * @return
	 */
	private static int checkAndGiveMaxPowerOfTwo(double[] reSrc, double[] imSrc,
			double[] reFt, double[] imFt) {
		int slen = reSrc.length;
		
		// Validate array lengths:
		if ((imSrc.length != slen) || (reFt.length != slen) || (imFt.length != slen)){
			throw new IllegalArgumentException("Arrays must be of same length.");
		}
		
		// Find the biggest fitting n; warn to stderr if array sizes are wrong.:
		// {after this 2^n is biggest that is at most slen.}
		int n=0;
		while((2<<n) <= slen) n++;
	
		final int flen = 1<<n;  // N=2^n;
		
		// System.out.printf("slen %d  flen %d%n", slen, flen);
		
		if (slen != flen){
			System.err.printf("WARNING: (FFT) Sample length is %d, "
					+"but FFT is done only on samples 1 to %d!%n",
					slen, flen);
		}
		return n;
	}

	/**
	 * Fourier transform directly using the mathematical definition. This is 
	 * O(N^2), and therefore quite useless.
	 * 
	 * @param src        re of source data
	 * @param isrc       im of source data
	 * @param transf     re of transformation
	 * @param itransf    im of transformation
	 */
	public static void fourierAsDefined(double[] src, double[] isrc,
			double[] transf, double[] itransf) {
		int ndata = src.length;
		
		double twopi = 2.0 * Math.PI;
		for (int k=0; k<ndata; k++){
			transf[k] = 0;
			itransf[k] = 0;
			for (int j=0;j<ndata;j++){
				int kj = k*j;
				double pow = kj * twopi / ndata;
				double x2=Math.cos(-pow);
				double y2=Math.sin(-pow);
				// (x1+iy1)(x2+iy2) = x1x2-y1y2 + i(x1y2+x2y1)
				transf[k] += src[j]*x2 - isrc[j]*y2;
				itransf[k] += src[j]*y2 + x2*isrc[j];
			}
			transf[k] /= Math.sqrt(ndata);
			itransf[k] /= Math.sqrt(ndata);
		}
	}
}