/*
  * Copyright 2007 Sun Microsystems, Inc.  All Rights Reserved.
  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
  *
  * This code is free software; you can redistribute it and/or modify it
  * under the terms of the GNU General Public License version 2 only, as
  * published by the Free Software Foundation.  Sun designates this
  * particular file as subject to the "Classpath" exception as provided
  * by Sun in the LICENSE file that accompanied this code.
  *
  * This code is distributed in the hope that it will be useful, but WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  * version 2 for more details (a copy is included in the LICENSE file that
  * accompanied this code).
  *
  * You should have received a copy of the GNU General Public License version
  * 2 along with this work; if not, write to the Free Software Foundation,
  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  *
  * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
  * CA 95054 USA or visit www.sun.com if you need additional information or
  * have any questions.
  */ 

package com.sun.media.sound;

/**
 * Fast Fourier Transformer.
 * 
 * @version %I%, %E%
 * @author Karl Helgason
 */ 
public final class FFT {

	private double[] w;

	private int fftFrameSize;

	private int sign;

	private int[] bitm_array;

	private int fftFrameSize2;

	// Sign = -1 is FFT, 1 is IFFT (inverse FFT)
	// Data = Interlaced double array to be transformed.
	// The order is: real (sin), complex (cos)
	// Framesize must be power of 2
	
	public FFT(int fftFrameSize, int sign) {
		w = computeTwiddleFactors(fftFrameSize, sign);

		this.fftFrameSize = fftFrameSize;
		this.sign = sign;
		fftFrameSize2 = fftFrameSize << 1;

		// Pre-process Bit-Reversal
		bitm_array = new int[fftFrameSize2];
		for (int i = 2; i < fftFrameSize2; i += 2) {
			int j;
			int bitm;
			for (bitm = 2, j = 0; bitm < fftFrameSize2; bitm <<= 1) {
				if ((i & bitm) != 0)
					j++;
				j <<= 1;
			}
			bitm_array[i] = j;
		}

	}

	public void transform(double[] data) {
		bitreversal(data);
		calc(fftFrameSize, data, sign, w);
	}

	private final static double[] computeTwiddleFactors(int fftFrameSize,
			int sign) {

		int imax = (int) (Math.log(fftFrameSize) / Math.log(2.));

		double[] warray = new double[(fftFrameSize - 1) * 4];
		int w_index = 0;

		for (int i = 0, nstep = 2; i < imax; i++) {

			int jmax = nstep;
			nstep <<= 1;

			double wr = 1.0;
			double wi = 0.0;

			double arg = Math.PI / (jmax >> 1);
			double wfr = Math.cos(arg);
			double wfi = sign * Math.sin(arg);

			for (int j = 0; j < jmax; j += 2) {

				warray[w_index++] = wr;
				warray[w_index++] = wi;

				double tempr = wr;
				wr = tempr * wfr - wi * wfi;
				wi = tempr * wfi + wi * wfr;
			}
		}

		// PRECOMPUTATION of wwr1, wwi1 for factor 4 Decomposition (3 * complex
		// operators and 8 +/- complex operators)
		{

			w_index = 0;
			int w_index2 = warray.length >> 1;
			for (int i = 0, nstep = 2; i < (imax - 1); i++) {

				int jmax = nstep;
				nstep *= 2;

				int ii = w_index + jmax;
				for (int j = 0; j < jmax; j += 2) {

					double wr = warray[w_index++];
					double wi = warray[w_index++];
					double wr1 = warray[ii++];
					double wi1 = warray[ii++];
					warray[w_index2++] = wr * wr1 - wi * wi1;
					warray[w_index2++] = wr * wi1 + wi * wr1;

				}
			}

		}

		return warray;

	}

	private final static void calc(int fftFrameSize, double[] data, int sign,
			double[] w)

	{

		final int fftFrameSize2 = fftFrameSize << 1;

		int nstep = 2;

		if (nstep >= fftFrameSize2)
			return;
		int i = nstep - 2;
		if (sign == -1)
			calcF4F(fftFrameSize, data, i, nstep, w);
		else
			calcF4I(fftFrameSize, data, i, nstep, w);

	}

	private final static void calcF2E(int fftFrameSize, double[] data, int i,
			int nstep, double[] w) {
		int jmax = nstep;
		nstep <<= 1;
		for (int n = 0; n < jmax; n += 2) {
			double wr = w[i++];
			double wi = w[i++];
			int m = n + jmax;
			double datam_r = data[m];
			double datam_i = data[m + 1];
			double datan_r = data[n];
			double datan_i = data[n + 1];
			double tempr = datam_r * wr - datam_i * wi;
			double tempi = datam_r * wi + datam_i * wr;
			data[m] = datan_r - tempr;
			data[m + 1] = datan_i - tempi;
			data[n] = datan_r + tempr;
			data[n + 1] = datan_i + tempi;
		}
		return;

	}

	// Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
	// complex operators

	private final static void calcF4F(int fftFrameSize, double[] data, int i,
			int nstep, double[] w) {
		final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
		// Factor-4 Decomposition

		int w_len = w.length >> 1;
		while (nstep < fftFrameSize2) {

			if (nstep << 2 == fftFrameSize2) {
				// Goto Factor-4 Final Decomposition
				// calcF4E(data, i, nstep, -1, w);
				calcF4FE(fftFrameSize, data, i, nstep, w);
				return;
			}
			int jmax = nstep;
			int nnstep = nstep << 1;
			if (nnstep == fftFrameSize2) {
				// Factor-4 Decomposition not possible
				calcF2E(fftFrameSize, data, i, nstep, w);
				return;
			}
			nstep <<= 2;
			int ii = i + jmax;
			int iii = i + w_len;

			{
				i += 2;
				ii += 2;
				iii += 2;

				for (int n = 0; n < fftFrameSize2; n += nstep) {
					int m = n + jmax;

					double datam1_r = data[m];
					double datam1_i = data[m + 1];
					double datan1_r = data[n];
					double datan1_i = data[n + 1];

					n += nnstep;
					m += nnstep;
					double datam2_r = data[m];
					double datam2_i = data[m + 1];
					double datan2_r = data[n];
					double datan2_i = data[n + 1];

					double tempr = datam1_r;
					double tempi = datam1_i;

					datam1_r = datan1_r - tempr;
					datam1_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					double n2w1r = datan2_r;
					double n2w1i = datan2_i;
					double m2ww1r = datam2_r;
					double m2ww1i = datam2_i;

					tempr = m2ww1r - n2w1r;
					tempi = m2ww1i - n2w1i;

					datam2_r = datam1_r + tempi;
					datam2_i = datam1_i - tempr;
					datam1_r = datam1_r - tempi;
					datam1_i = datam1_i + tempr;

					tempr = n2w1r + m2ww1r;
					tempi = n2w1i + m2ww1i;

					datan2_r = datan1_r - tempr;
					datan2_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					data[m] = datam2_r;
					data[m + 1] = datam2_i;
					data[n] = datan2_r;
					data[n + 1] = datan2_i;

					n -= nnstep;
					m -= nnstep;
					data[m] = datam1_r;
					data[m + 1] = datam1_i;
					data[n] = datan1_r;
					data[n + 1] = datan1_i;

				}

			}
			for (int j = 2; j < jmax; j += 2) {
				double wr = w[i++];
				double wi = w[i++];
				double wr1 = w[ii++];
				double wi1 = w[ii++];
				double wwr1 = w[iii++];
				double wwi1 = w[iii++];
				// double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
				// precomputed!!!
				// double wwi1 = wr * wi1 + wi * wr1;

				for (int n = j; n < fftFrameSize2; n += nstep) {
					int m = n + jmax;

					double datam1_r = data[m];
					double datam1_i = data[m + 1];
					double datan1_r = data[n];
					double datan1_i = data[n + 1];

					n += nnstep;
					m += nnstep;
					double datam2_r = data[m];
					double datam2_i = data[m + 1];
					double datan2_r = data[n];
					double datan2_i = data[n + 1];

					double tempr = datam1_r * wr - datam1_i * wi;
					double tempi = datam1_r * wi + datam1_i * wr;

					datam1_r = datan1_r - tempr;
					datam1_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					double n2w1r = datan2_r * wr1 - datan2_i * wi1;
					double n2w1i = datan2_r * wi1 + datan2_i * wr1;
					double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
					double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;

					tempr = m2ww1r - n2w1r;
					tempi = m2ww1i - n2w1i;

					datam2_r = datam1_r + tempi;
					datam2_i = datam1_i - tempr;
					datam1_r = datam1_r - tempi;
					datam1_i = datam1_i + tempr;

					tempr = n2w1r + m2ww1r;
					tempi = n2w1i + m2ww1i;

					datan2_r = datan1_r - tempr;
					datan2_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					data[m] = datam2_r;
					data[m + 1] = datam2_i;
					data[n] = datan2_r;
					data[n + 1] = datan2_i;

					n -= nnstep;
					m -= nnstep;
					data[m] = datam1_r;
					data[m + 1] = datam1_i;
					data[n] = datan1_r;
					data[n + 1] = datan1_i;

				}
			}

			i += jmax << 1;

		}

		calcF2E(fftFrameSize, data, i, nstep, w);

	}

	// Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
	// complex operators

	private final static void calcF4I(int fftFrameSize, double[] data, int i,
			int nstep, double[] w) {
		final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
		// Factor-4 Decomposition

		int w_len = w.length >> 1;
		while (nstep < fftFrameSize2) {

			if (nstep << 2 == fftFrameSize2) {
				// Goto Factor-4 Final Decomposition
				// calcF4E(data, i, nstep, 1, w);
				calcF4IE(fftFrameSize, data, i, nstep, w);
				return;
			}
			int jmax = nstep;
			int nnstep = nstep << 1;
			if (nnstep == fftFrameSize2) {
				// Factor-4 Decomposition not possible
				calcF2E(fftFrameSize, data, i, nstep, w);
				return;
			}
			nstep <<= 2;
			int ii = i + jmax;
			int iii = i + w_len;
			{
				i += 2;
				ii += 2;
				iii += 2;

				for (int n = 0; n < fftFrameSize2; n += nstep) {
					int m = n + jmax;

					double datam1_r = data[m];
					double datam1_i = data[m + 1];
					double datan1_r = data[n];
					double datan1_i = data[n + 1];

					n += nnstep;
					m += nnstep;
					double datam2_r = data[m];
					double datam2_i = data[m + 1];
					double datan2_r = data[n];
					double datan2_i = data[n + 1];

					double tempr = datam1_r;
					double tempi = datam1_i;

					datam1_r = datan1_r - tempr;
					datam1_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					double n2w1r = datan2_r;
					double n2w1i = datan2_i;
					double m2ww1r = datam2_r;
					double m2ww1i = datam2_i;

					tempr = n2w1r - m2ww1r;
					tempi = n2w1i - m2ww1i;

					datam2_r = datam1_r + tempi;
					datam2_i = datam1_i - tempr;
					datam1_r = datam1_r - tempi;
					datam1_i = datam1_i + tempr;

					tempr = n2w1r + m2ww1r;
					tempi = n2w1i + m2ww1i;

					datan2_r = datan1_r - tempr;
					datan2_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					data[m] = datam2_r;
					data[m + 1] = datam2_i;
					data[n] = datan2_r;
					data[n + 1] = datan2_i;

					n -= nnstep;
					m -= nnstep;
					data[m] = datam1_r;
					data[m + 1] = datam1_i;
					data[n] = datan1_r;
					data[n + 1] = datan1_i;

				}

			}
			for (int j = 2; j < jmax; j += 2) {
				double wr = w[i++];
				double wi = w[i++];
				double wr1 = w[ii++];
				double wi1 = w[ii++];
				double wwr1 = w[iii++];
				double wwi1 = w[iii++];
				// double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
				// precomputed!!!
				// double wwi1 = wr * wi1 + wi * wr1;

				for (int n = j; n < fftFrameSize2; n += nstep) {
					int m = n + jmax;

					double datam1_r = data[m];
					double datam1_i = data[m + 1];
					double datan1_r = data[n];
					double datan1_i = data[n + 1];

					n += nnstep;
					m += nnstep;
					double datam2_r = data[m];
					double datam2_i = data[m + 1];
					double datan2_r = data[n];
					double datan2_i = data[n + 1];

					double tempr = datam1_r * wr - datam1_i * wi;
					double tempi = datam1_r * wi + datam1_i * wr;

					datam1_r = datan1_r - tempr;
					datam1_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					double n2w1r = datan2_r * wr1 - datan2_i * wi1;
					double n2w1i = datan2_r * wi1 + datan2_i * wr1;
					double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
					double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;

					tempr = n2w1r - m2ww1r;
					tempi = n2w1i - m2ww1i;

					datam2_r = datam1_r + tempi;
					datam2_i = datam1_i - tempr;
					datam1_r = datam1_r - tempi;
					datam1_i = datam1_i + tempr;

					tempr = n2w1r + m2ww1r;
					tempi = n2w1i + m2ww1i;

					datan2_r = datan1_r - tempr;
					datan2_i = datan1_i - tempi;
					datan1_r = datan1_r + tempr;
					datan1_i = datan1_i + tempi;

					data[m] = datam2_r;
					data[m + 1] = datam2_i;
					data[n] = datan2_r;
					data[n + 1] = datan2_i;

					n -= nnstep;
					m -= nnstep;
					data[m] = datam1_r;
					data[m + 1] = datam1_i;
					data[n] = datan1_r;
					data[n + 1] = datan1_i;

				}
			}

			i += jmax << 1;

		}

		calcF2E(fftFrameSize, data, i, nstep, w);

	}

	// Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
	// complex operators

	private final static void calcF4FE(int fftFrameSize, double[] data, int i,
			int nstep, double[] w) {
		final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
		// Factor-4 Decomposition

		int w_len = w.length >> 1;
		while (nstep < fftFrameSize2) {

			int jmax = nstep;
			int nnstep = nstep << 1;
			if (nnstep == fftFrameSize2) {
				// Factor-4 Decomposition not possible
				calcF2E(fftFrameSize, data, i, nstep, w);
				return;
			}
			nstep <<= 2;
			int ii = i + jmax;
			int iii = i + w_len;
			for (int n = 0; n < jmax; n += 2) {
				double wr = w[i++];
				double wi = w[i++];
				double wr1 = w[ii++];
				double wi1 = w[ii++];
				double wwr1 = w[iii++];
				double wwi1 = w[iii++];
				// double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
				// precomputed!!!
				// double wwi1 = wr * wi1 + wi * wr1;

				int m = n + jmax;

				double datam1_r = data[m];
				double datam1_i = data[m + 1];
				double datan1_r = data[n];
				double datan1_i = data[n + 1];

				n += nnstep;
				m += nnstep;
				double datam2_r = data[m];
				double datam2_i = data[m + 1];
				double datan2_r = data[n];
				double datan2_i = data[n + 1];

				double tempr = datam1_r * wr - datam1_i * wi;
				double tempi = datam1_r * wi + datam1_i * wr;

				datam1_r = datan1_r - tempr;
				datam1_i = datan1_i - tempi;
				datan1_r = datan1_r + tempr;
				datan1_i = datan1_i + tempi;

				double n2w1r = datan2_r * wr1 - datan2_i * wi1;
				double n2w1i = datan2_r * wi1 + datan2_i * wr1;
				double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
				double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;

				tempr = m2ww1r - n2w1r;
				tempi = m2ww1i - n2w1i;

				datam2_r = datam1_r + tempi;
				datam2_i = datam1_i - tempr;
				datam1_r = datam1_r - tempi;
				datam1_i = datam1_i + tempr;

				tempr = n2w1r + m2ww1r;
				tempi = n2w1i + m2ww1i;

				datan2_r = datan1_r - tempr;
				datan2_i = datan1_i - tempi;
				datan1_r = datan1_r + tempr;
				datan1_i = datan1_i + tempi;

				data[m] = datam2_r;
				data[m + 1] = datam2_i;
				data[n] = datan2_r;
				data[n + 1] = datan2_i;

				n -= nnstep;
				m -= nnstep;
				data[m] = datam1_r;
				data[m + 1] = datam1_i;
				data[n] = datan1_r;
				data[n + 1] = datan1_i;

			}

			i += jmax << 1;

		}

	}

	// Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
	// complex operators

	private final static void calcF4IE(int fftFrameSize, double[] data, int i,
			int nstep, double[] w) {
		final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
		// Factor-4 Decomposition

		int w_len = w.length >> 1;
		while (nstep < fftFrameSize2) {

			int jmax = nstep;
			int nnstep = nstep << 1;
			if (nnstep == fftFrameSize2) {
				// Factor-4 Decomposition not possible
				calcF2E(fftFrameSize, data, i, nstep, w);
				return;
			}
			nstep <<= 2;
			int ii = i + jmax;
			int iii = i + w_len;
			for (int n = 0; n < jmax; n += 2) {
				double wr = w[i++];
				double wi = w[i++];
				double wr1 = w[ii++];
				double wi1 = w[ii++];
				double wwr1 = w[iii++];
				double wwi1 = w[iii++];
				// double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
				// precomputed!!!
				// double wwi1 = wr * wi1 + wi * wr1;

				int m = n + jmax;

				double datam1_r = data[m];
				double datam1_i = data[m + 1];
				double datan1_r = data[n];
				double datan1_i = data[n + 1];

				n += nnstep;
				m += nnstep;
				double datam2_r = data[m];
				double datam2_i = data[m + 1];
				double datan2_r = data[n];
				double datan2_i = data[n + 1];

				double tempr = datam1_r * wr - datam1_i * wi;
				double tempi = datam1_r * wi + datam1_i * wr;

				datam1_r = datan1_r - tempr;
				datam1_i = datan1_i - tempi;
				datan1_r = datan1_r + tempr;
				datan1_i = datan1_i + tempi;

				double n2w1r = datan2_r * wr1 - datan2_i * wi1;
				double n2w1i = datan2_r * wi1 + datan2_i * wr1;
				double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
				double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;

				tempr = n2w1r - m2ww1r;
				tempi = n2w1i - m2ww1i;

				datam2_r = datam1_r + tempi;
				datam2_i = datam1_i - tempr;
				datam1_r = datam1_r - tempi;
				datam1_i = datam1_i + tempr;

				tempr = n2w1r + m2ww1r;
				tempi = n2w1i + m2ww1i;

				datan2_r = datan1_r - tempr;
				datan2_i = datan1_i - tempi;
				datan1_r = datan1_r + tempr;
				datan1_i = datan1_i + tempi;

				data[m] = datam2_r;
				data[m + 1] = datam2_i;
				data[n] = datan2_r;
				data[n + 1] = datan2_i;

				n -= nnstep;
				m -= nnstep;
				data[m] = datam1_r;
				data[m + 1] = datam1_i;
				data[n] = datan1_r;
				data[n + 1] = datan1_i;

			}

			i += jmax << 1;

		}

	}

	private final void bitreversal(double[] data) {
		if (fftFrameSize < 4)
			return;

		int inverse = fftFrameSize2 - 2;
		for (int i = 0; i < fftFrameSize; i += 4) {
			int j = bitm_array[i];

			// Performing Bit-Reversal, even v.s. even, O(2N)
			if (i < j) {

				int n = i;
				int m = j;

				// COMPLEX: SWAP(data[n], data[m])
				// Real Part
				double tempr = data[n];
				data[n] = data[m];
				data[m] = tempr;
				// Imagery Part
				n++;
				m++;
				double tempi = data[n];
				data[n] = data[m];
				data[m] = tempi;

				n = inverse - i;
				m = inverse - j;

				// COMPLEX: SWAP(data[n], data[m])
				// Real Part
				tempr = data[n];
				data[n] = data[m];
				data[m] = tempr;
				// Imagery Part
				n++;
				m++;
				tempi = data[n];
				data[n] = data[m];
				data[m] = tempi;
			}

			// Performing Bit-Reversal, odd v.s. even, O(N)

			int m = j + fftFrameSize; // bitm_array[i+2];
			// COMPLEX: SWAP(data[n], data[m])
			// Real Part
			int n = i + 2;
			double tempr = data[n];
			data[n] = data[m];
			data[m] = tempr;
			// Imagery Part
			n++;
			m++;
			double tempi = data[n];
			data[n] = data[m];
			data[m] = tempi;
		}

	}

}