/*******************************************************************************
 * Copyright (c) 2010 Attila Kovacs <attila_kovacs@post.harvard.edu>.
 * All rights reserved. 
 * 
 * This file is part of crush.
 * 
 *     crush is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by
 *     the Free Software Foundation, either version 3 of the License, or
 *     (at your option) any later version.
 * 
 *     crush 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 for more details.
 * 
 *     You should have received a copy of the GNU General Public License
 *     along with crush.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Contributors:
 *     Attila Kovacs <attila_kovacs@post.harvard.edu> - initial API and implementation
 ******************************************************************************/
// Copyright (c) 2010 Attila Kovacs 

//Copyright (c) 2007 Attila Kovacs 

package util;

//package crush.util;

import java.io.*;
import java.text.*;
import java.util.*;


public final class ArrayUtil {

	// A bunch of FFT related stuff here based on Numerical Recipes FFT routines.
	//double[][] getChunks(double[] data, long chunkSize) {}

	public static void print(double[] data) { print(data, System.out, new DecimalFormat("0.000000E0")); }

	public static void print(double[] data, DecimalFormat df) { print(data, System.out, df); }

	public static void print(double[] data, PrintStream out, DecimalFormat df) {
		for(int i=0; i<data.length; i++) out.println(df.format(data[i]));
	}

	public static double median(double[] data) { return median(data, 0, data.length); }

	public static double median(double[] data, int fromIndex, int toIndex) {
		Arrays.sort(data, fromIndex, toIndex);
		int n = toIndex - fromIndex;
		return n % 2 == 0 ? 0.5 * (data[fromIndex + n/2-1] + data[fromIndex + n/2]) : data[fromIndex + (n-1)/2];
	}

	public static float median(float[] data) { return median(data, 0, data.length); }

	public static float median(float[] data, int fromIndex, int toIndex) {
		Arrays.sort(data, fromIndex, toIndex);
		int n = toIndex - fromIndex;
		return n % 2 == 0 ? 0.5F * (data[fromIndex + n/2-1] + data[fromIndex + n/2]) : data[fromIndex + (n-1)/2];
	}

	public static double median(WeightedPoint[] data) { return median(data, 0, data.length); }
	
	public static double median(WeightedPoint[] data, int fromIndex, int toIndex) {
		return smartMedian(data, fromIndex, toIndex, 1.0);		
	}
	
	public static double smartMedian(final WeightedPoint[] data, final int fromIndex, final int toIndex, final double maxDependence) {
		if(toIndex - fromIndex == 1) return data[fromIndex].value;

		Arrays.sort(data, fromIndex, toIndex);

		// wt is the sum of all weights
		// wi is the integral sum including the current point.
		double wt = 0.0, wmax = 0.0;
		
		for(int i=fromIndex; i<toIndex; i++) {
			final double w = data[i].weight;
			if(w > wmax) wmax = w;
			if(w > 0.0) wt += w;
		}

		// If a single datum dominates the data, then return the weighted mean...
		if(wmax >= maxDependence * wt) {
			double sum=0.0, sumw=0.0;
			for(int i=fromIndex; i<toIndex; i++) {
				final double w = data[i].weight;
				if(w > 0.0) {
					sum += w * data[i].value;
					sumw += w;
				}
			}
			return sum/sumw;
		}
		
		// If all weights are zero return the arithmetic median...
		// This should never happen, but just in case...
		if(wt == 0.0) {
			final int n = toIndex - fromIndex;
			return n % 2 == 0 ? 0.5F * (data[fromIndex + n/2-1].value + data[fromIndex + n/2].value) : data[fromIndex + (n-1)/2].value;
		}


		final double midw = wt / 2.0; 
		int ig = fromIndex; 
		
		WeightedPoint last = WeightedPoint.NaN;
		WeightedPoint point = data[fromIndex];

		double wi = point.weight;
		
		
		while(wi < midw) if(data[++ig].weight > 0.0) {
			last = point;
			point = data[ig];	    
			wi += 0.5 * (last.weight + point.weight);    
		}
		
		double wplus = wi;
		double wminus = wi - 0.5 * (last.weight + point.weight);
		
		double w1 = (wplus - midw) / (wplus + wminus);
		return w1 * last.value + (1.0-w1) * point.value;
	}

	// Interpret integral normalized spectrum as signal amplitudes
	public static double[] amplitude(Complex[] data) {
		double[] amp = new double[data.length];
		for(int i=0; i<data.length; i++) amp[i] = data[i].length();
		amp[0] /= 2.0;
		amp[amp.length - 1] /= 2.0;
		return amp;
	}

	// Interpret integral normalized spectrum as noise amplitudes
	public static double[] noiseAmplitude(Complex[] data) {
		double[] amp = new double[data.length];
		for(int i=0; i<data.length; i++) amp[i] = data[i].length();
		return amp;
	}

	// Interpret integral normalized spectrum as signal power
	public static double[] norm(Complex[] data) {
		double[] norm = new double[data.length];
		for(int i=0; i<data.length; i++) norm[i] = data[i].norm();
		norm[0] /= 4.0;
		norm[norm.length - 1] /= 4.0;	
		return norm;
	}

	// Interpret integral normalized spectrum as noise power
	public static double[] noiseNorm(Complex[] data) {
		double[] norm = new double[data.length];
		for(int i=0; i<data.length; i++) norm[i] = data[i].norm();
		return norm;
	}

	public static double[] real(Complex[] data) {
		double[] real = new double[data.length];
		for(int i=0; i<data.length; i++) real[i] = data[i].x;
		return real;
	}

	public static double[] re(Complex[] data) { return real(data); }

	public static double[] imaginary(Complex[] data) {
		double[] imaginary = new double[data.length];
		for(int i=0; i<data.length; i++) imaginary[i] = data[i].y;
		return imaginary;
	}


	public static double[] im(Complex[] data) { return imaginary(data); }

	public static int getPaddedFFTSize(int size) {
		int n=1;	
		while(n < size) n<<=1;
		return n;
	}

	public static int getTruncatedFFTSize(int size) {
		int n=1;	
		while(n <= size) n<<=1;
		return n>>1;
	}

	public static double[] truncateForFFT(double[] data) {
		return padWithZeros(data, getTruncatedFFTSize(data.length));
	}

	public static Vector2D[] truncateForFFT(Complex[] data) {
		int n = getTruncatedFFTSize(data.length);	
		return data.length == n ? data : padWithZeros(data, n);
	}

	public static double[] padForFFT(double[] data) {
		int n = getPaddedFFTSize(data.length);
		return data.length == n ? data : padWithZeros(data, n);
	}

	public static Complex[] padForFFT(Complex[] data) {
		return padWithZeros(data, getPaddedFFTSize(data.length));
	}

	public static double[] padWithZeros(double[] data, int n) {
		double[] padded = new double[n];

		int N = Math.min(data.length, n);

		System.arraycopy(data, 0, padded, 0, N);

		return padded;
	}

	public static Complex[] padWithZeros(Complex[] data, int n) {
		Complex[] padded = new Complex[n];

		int N = Math.min(data.length, n);

		System.arraycopy(data, 0, padded, 0, N);
		//for(int i=0; i<N; i++) padded[i] = (Complex) data[i].clone();

		return padded;
	}

	public static String[] getWindowFunctionNames() {
		return new String[] { "Rectangular", "Hamming", "Hann", "Blackman", "Nutall", 
				"Blackman-Harris", "Blackman-Nutall", "Flat Top"};
	}

	public static double getEquivalentWindowWidth(String name) {
		if(name.equalsIgnoreCase("Rectangular") || name.equalsIgnoreCase("Box")) return 1.0;
		if(name.equalsIgnoreCase("Hamming")) return 1.37;
		if(name.equalsIgnoreCase("Hann")) return 1.50;
		if(name.equalsIgnoreCase("Blackman")) return 1.73;
		if(name.equalsIgnoreCase("Nutall")) return 2.02;
		if(name.equalsIgnoreCase("Blackman-Harris")) return 2.01;
		if(name.equalsIgnoreCase("Blackman-Nutall")) return 1.98;
		if(name.equalsIgnoreCase("Flat Top")) return 3.77;
		return -1.0;
	}

	public static double[] getWindowFunction(String name, int n) {
		if(name.equalsIgnoreCase("Rectangular") || name.equalsIgnoreCase("Box")) return getRectangularWindow(n);
		if(name.equalsIgnoreCase("Hamming")) return getHammingWindow(n);
		if(name.equalsIgnoreCase("Hann")) return getHannWindow(n);
		if(name.equalsIgnoreCase("Blackman")) return getBlackmanWindow(n);
		if(name.equalsIgnoreCase("Nutall")) return getNutallWindow(n);
		if(name.equalsIgnoreCase("Blackman-Harris")) return getBlackmanHarrisWindow(n);
		if(name.equalsIgnoreCase("Blackman-Nutall")) return getBlackmanNutallWindow(n);
		if(name.equalsIgnoreCase("Flat Top")) return getFlatTopWindow(n);
		return null;
	}

	// W = 1.00
	// -15dB
	public static double[] getRectangularWindow(int n) {
		return getWindowFunction(n, new double[] { 1.0 });
	}

	// W = 1.37
	// -40dB
	public static double[] getHammingWindow(int n) {
		return getWindowFunction(n, new double[] { 0.53836, -0.46164 });
	}

	// W = 1.50
	// -30dB
	public static double[] getHannWindow(int n) {
		return getWindowFunction(n, new double[] { 0.5, -0.5 });
	}

	// W = 1.73
	// -60dB
	public static double[] getBlackmanWindow(int n) {
		return getWindowFunction(n, new double[] { 0.5, -0.5 });
	}

	// W = 2.02
	// -100dB
	public static double[] getNutallWindow(int n) {
		return getWindowFunction(n, new double[] { 0.355768, -0.487396, 0.144232, -0.012604 });
	}

	// W = 2.01
	// -100dB
	public static double[] getBlackmanHarrisWindow(int n) {
		return getWindowFunction(n, new double[] { 0.35875, -0.48829, 0.14128, -0.01168 });
	}

	// W = 1.98
	// -100dB
	public static double[] getBlackmanNutallWindow(int n) {
		return getWindowFunction(n, new double[] { 0.3635819, -0.4891775, 0.1365995, -0.0106411 });
	}

	// W = 3.77
	// -70dB
	public static double[] getFlatTopWindow(int n) {
		return getWindowFunction(n, new double[] { 1.0, -1.93, 1.29, -0.388, 0.032 });
	}

	// normalize to sum w^2
	// Needed for provididng true PSD estimates...
	public static double[] getWindowFunction(int n, double[] coeff) {
		double[] w = new double[n];
		double norm = 0.0;
		for(int i=0; i<n; i++) {
			for(int k=0; k<coeff.length; k++) w[i] += coeff[k] * Math.cos(2.0*k*Math.PI * i/(n-1));
			norm += w[i]*w[i];
		}
		for(int i=0; i<n; i++) w[i] /= norm;

		return w;
	}
	
	// Rewritten to skip costly intermediate Complex storage...
	public static double[] convolveViaFFT(double[] data, double[] beam) {
		return convolveViaFFT(data, beam, false);		
	}
	
	public static double[] convolveViaFFT(double[] data, double[] beam, boolean inPlace) {
		final int n = inPlace ? data.length : getPaddedFFTSize(data.length + beam.length);
		final int nF = n >> 1;
		
		final double[] a = inPlace ? data : new double[n];
		final double[] b = inPlace ? beam : new double[n];
		if(!inPlace) {
			System.arraycopy(data, 0, a, 0, data.length);
			System.arraycopy(beam, 0, b, 0, data.length);
		}
			
		inplaceRealFFT(a, true);
		inplaceRealFFT(b, true);
	
		final double A = 1.0 / (nF * nF);
		
		a[0] *= b[0];
		a[1] *= b[1];
		
		for(int i=2; i<n; i+=2) {
			final double temp = a[i];
			final int j = i+1;
			a[i] = A * (a[i] * a[i] - b[i] * b[i]) ;
			a[j] = A * (temp * b[j] + a[j] * b[i]);
		}
		
		inplaceRealFFT(a, false);
		
		return a;
	}

	public static double[] autoCorrelateViaFFT(double[] data) { return convolveViaFFT(data, data); }

	public static double[] averagePower(double[] data, int windowSize) {
		if(data.length < windowSize) return averagePower(padWithZeros(data, windowSize), windowSize);
		return averagePower(data, getHammingWindow(windowSize));
	}

	// Rewritten to skip costly intermediate Complex storage...
	public static double[] averagePower(double[] data, double[] w) {
		int windowSize = w.length;
		int stepSize = windowSize >> 1;
		final double[] block = new double[getPaddedFFTSize(w.length)];
		final int nF = block.length >> 1;
		
		// Create the accumulated spectrum array
		double[] spectrum = null;

		int start = 0, N = 0;
		while(start + windowSize <= data.length) {

			for(int i=0,j=start; i<windowSize; i++,j++) block[i] = w[i] * data[j];	
			Arrays.fill(block, windowSize, block.length, 0.0);
			inplaceRealFFT(block, true);
			
			if(spectrum == null) spectrum = new double[nF + 1];

			spectrum[0] += block[0] * block[0];
			spectrum[nF] += block[1] * block[1];
			
			for(int i=nF; --i>=1; ) {
				int j = i << 1;
				spectrum[i] += block[j] * block[j];
				spectrum[i] += block[++j] * block[j];
			}

			start += stepSize;

			N++;
		}

		// Should not use amplitude normalization here but power...
		// The spectal power per frequency component.
		double norm = (double) 1.0 / N;

		for(int i=0; i<spectrum.length; i++) spectrum[i] *= norm;

		return spectrum;	
	}

	/*
	public static Complex[] overlapFFT(double[] data, int windowSize) {
		if(data.length < windowSize) return overlapFFT(padWithZeros(data, windowSize), windowSize);

		// figure out the FFT size
		int n=1;	
		while(n < windowSize) n<<=1;
		double[] block = new double[n];

		// Create the cumulated spectrum array
		Complex[] spectrum = null;

		int start = 0, N=0;
		while(start + windowSize <= data.length) {	    

			for(int i=0,j=start; i<windowSize; i++,j++) block[i] = data[j];	    

			Complex[] blockSpectrum = FFT(block);

			if(spectrum == null) spectrum = blockSpectrum;
			else for(int i=0; i<spectrum.length; i++) spectrum[i].add(blockSpectrum[i]);

			start += windowSize/2;

			N++;
		}

		double norm = 1.0 / N;

		for(int i=0; i<spectrum.length; i++) spectrum[i].scale(norm);

		return spectrum;	
	}
	*/

	public static Complex[] FFT(double[] data) { return FFT(data, getPaddedFFTSize(data.length)); }

	public static double[][] backFFT(Complex[][] spectrum) {
		return null;
	}

	public static Complex[] FFT(double[] data, int n) {
		double[] tdata = padWithZeros(data, n);
		int N = n/2;

		Complex[] fdata = new Complex[N+1];
		for(int i=N+1; --i >= 0; ) fdata[i] = new Complex();

		uncheckedFFT(tdata, fdata);

		return fdata;
	}

	// This performs the fft on tdata, and puts the
	// apmlitude normalized frequency data into the supplied Complex array.
	// but the tdata is also destroyed in the process...
	public static void uncheckedFFT(final double[] data, final Complex[] fdata) {
		final int N = data.length>>1;
		final double norm = 1.0 / N;

		inplaceRealFFT(data, true);

		fdata[0].x = norm * data[0];
		fdata[N].x = norm * data[1];

		for(int i=1,ii=2; i<N; i++) {
			fdata[i].x = norm * data[ii++];
			fdata[i].y = norm * data[ii++];
		}
	}
	
	
	// This performs the fft on tdata, and puts the
	// apmlitude normalized frequency data into the supplied Complex array.
	// but the tdata is also destroyed in the process...
	public static void inplaceRealFFT(final double[] data) {
		inplaceRealFFT(data, true);
		
		final double norm = 2.0 / data.length;
		for(int i=0; i<data.length; i++) data[i] *= norm;
	}

	public static void uncheckedBackFFT(final Complex[] fdata, final double[] data) {
		final int N = data.length>>1;

		data[0] = fdata[0].x;
		data[1] = fdata[N].x;

		for(int i=1,ii=2; i<N; i++) {
			data[ii++] = fdata[i].x;
			data[ii++] = fdata[i].y;
		}

		inplaceRealFFT(data, false);
	}

	public static void inplaceRealBackFFT(double[] data) { inplaceRealFFT(data, false); }


	public static double[] backFFT(Complex[] spectrum) {
		return backFFT(spectrum, getPaddedFFTSize(spectrum.length-1));
	}

	public static double[] backFFT(Complex[] spectrum, int n) {
		double[] tdata = new double[2 * n];
		uncheckedBackFFT(spectrum, tdata);	
		return tdata;
	}

	public static void uncheckedFFT(final Complex[] data, final boolean isForward, double[] temp) {
		if(temp == null) temp = new double[2*data.length];
		else if(temp.length != 2*data.length) temp = new double[2*data.length];

		for(int i=0, ii=0; i<data.length; i++) {
			temp[ii++] = data[i].x;
			temp[ii++] = data[i].y;
		}

		inplaceFFT(temp, isForward);

		double norm = isForward ? 1.0 / data.length : 1.0;

		for(int i=0, ii=0; i<data.length; i++) {
			data[i].x = norm * temp[ii++];
			data[i].y = norm * temp[ii++];
		}	    
	}

	public static void uncheckedFFT(Complex[][] data, boolean isForward) {
		Complex[] ctemp = null;
		double[] dtemp = null;
		uncheckedFFT(data, isForward, ctemp, dtemp);
	}

	public static void uncheckedFFT(Complex[][] data, boolean isForward, Complex[] temp, double[] dtemp) {
		int nx = data.length;
		int ny = data[0].length;

		if(temp == null) temp = new Complex[nx];
		else if(temp.length != nx) temp = new Complex[nx];

		for(int i=0; i<nx; i++) {
			if(temp[i] == null) temp[i] = new Complex();
			uncheckedFFT(data[i], isForward, dtemp);
		}

		for(int j=0; j<ny; j++) {
			for(int i=0; i<nx; i++) temp[i] = data[i][j];
			uncheckedFFT(temp, isForward, dtemp);
			for(int i=0; i<nx; i++) data[i][j] = temp[i];
		}
	}

	public static Complex[][] loadFFT(double[][] data) {
		Complex[][] cdata = null;
		return loadFFT(data, cdata);
	}

	public static Complex[][] loadFFT(double[][] data, Complex[][] cdata) {
		int nx = data.length;
		int ny = data[0].length;

		int nu = getPaddedFFTSize(nx);
		int nv = getPaddedFFTSize(ny);

		if(cdata == null) cdata = new Complex[nu][nv];
		else if(cdata.length != nu || cdata[0].length != nv) cdata = new Complex[nu][nv];

		for(int i=0; i<nu; i++) for(int j=0; j<nv; j++) {
			if(cdata[i][j] == null) cdata[i][j] = new Complex();
			cdata[i][j].x = (i<nx && j<ny) ? data[i][j] : 0.0;
			cdata[i][j].y = 0.0;
		} 

		return cdata;
	}

	public static void unloadFFT(Complex[][] cdata, double[][] data) {
		int nx = data.length;
		int ny = data[0].length;

		for(int i=0; i<nx; i++) for(int j=0; j<ny; j++) data[i][j] = cdata[i][j].x;
	}

	public static double[][] amplitude(Complex[][] data) {
		int nx = data.length;
		int ny = data[0].length;

		double[][] a = new double[nx][ny];

		for(int i=0; i<nx; i++) for(int j=0; j<ny; j++) a[i][j] = data[i][j].length();

		// Change from integral to amplitide normalization
		for(int i=0; i<nx; i++) {
			a[i][0] /= 2.0;
			a[i][ny-1] /= 2.0;
		}

		for(int j=0; j<ny; j++) {
			a[0][j] /= 2.0;
			a[nx-1][j] /= 2.0;
		}

		return a;
	}

	public static double[][] norm(Complex[][] data) {
		int nx = data.length;
		int ny = data[0].length;

		double[][] a = new double[nx][ny];

		for(int i=0; i<nx; i++) for(int j=0; j<ny; j++) a[i][j] = data[i][j].norm();

		return a;
	}


	// Loosely based on the Numberical Recipes routine realft.c
	public static void inplaceRealFFT(final double data[], final boolean forward) {
		if(forward) inplaceFFT(data, true);
		
		final int N = data.length>>1;
		
		final double uh = 0.5;
		final double sh = forward ? 0.5 : -0.5;
		final double theta = (forward ? 1 : -1) * Math.PI / N;
		final double c = Math.cos(theta);
		final double s = Math.sin(theta);
		
		double wr = c;
		double wi = s;
	
		for(int r1=2, r2=data.length-2; r1<N; r1+=2, r2-=2) {
			final int i1 = r1 + 1;
			final int i2 = r2 + 1;

			double hr = sh * (data[i1] + data[i2]);
			double hi = sh * (data[r2] - data[r1]);
			final double r = wr * hr - wi * hi;
			final double i = wr * hi + wi * hr;

			hr = uh * (data[r1] + data[r2]);
			hi = uh * (data[i1] - data[i2]);
			
			data[r1] = hr + r;
			data[i1] = hi + i;
			data[r2] = hr - r;
			data[i2] = i - hi;
			
			final double temp = wr;
			wr = temp * c - wi * s;
			wi = wi * c + temp * s;
		}
		
		final double d0 = data[0];

		if(forward) {
			data[0] = d0 + data[1];
			data[1] = d0 - data[1];
		} 
		else {
			data[0] = uh * (d0 + data[1]);
			data[1] = uh * (d0 - data[1]);
			inplaceFFT(data, false);
		}
	}
	
	public static void bitReverse(final double[] data) {
		for(int i=0, j=0; i<data.length; i+=2) {
			if(j > i) {	
				double temp = data[i]; data[i] = data[j]; data[j] = temp;
				temp = data[i+1]; data[i+1] = data[j+1]; data[j+1] = temp;
			}

			int m = data.length >> 1;
			while(m >= 2 && j >= m) {
				j -= m;
				m >>= 1;
			}
			j += m;
		}
	}
	
	// Loosely based on the Numerical Recipes routine four1.c
	public static void inplaceFFT(final double data[], final boolean isForward) {
		bitReverse(data);
		
		int blockSize = 2;
		final double sg2Pi = (isForward ? 1: -1) * 2.0 * Math.PI;		
		
		while(blockSize < data.length) {
			final int step = (blockSize << 1) - 1;
			final double theta = sg2Pi / blockSize;
			final double c = Math.cos(theta);
			final double s = Math.sin(theta);
			
			double wr = 1.0;
			double wi = 0.0;
				
			for(int m=0; m<blockSize; m+=2) {
				for(int i1=m; i1<data.length; i1+=step) {
					int i2 = i1 + blockSize;
					int i2p = i2 + 1;
					
					final double d2r = data[i2];
					final double d2i = data[i2p];
					
					double x = wr * d2r - wi * d2i;
					data[i2] = data[i1] - x;
					data[i1] += x;
					
					x = wr * d2i + wi * d2r;
					data[i2p] = data[++i1] - x;
					data[i1] += x;
				}
				final double temp = wr;
				wr = temp * c - wi * s;
				wi = wi * c + temp * s;
			}
			blockSize <<= 1;
			Thread.yield();
		}
	} 
	
	public static void bitReverse(final Complex[] data) {
		Complex temp = null;
		final int halfLength = data.length >> 1;
		
		for(int i=0, j=0; i<data.length; i++) {
			if (j > i) { temp = data[i]; data[i]=data[j]; data[j] = temp; }
			int k = halfLength;	
			while(k > 0 && j >= k) {
				j -= k;
				k >>= 1;
			}
			j += k;
		}	
		
	}
	
	

	// First element is both the zero freq and the cutoff freq!
	public static void FFT(final Complex data[], boolean isForward) {
		bitReverse(data);
		
		final double sgPi = (isForward ? 1: -1) * Math.PI;
		int blockSize = 1;
		final Complex w = new Complex();
		final Complex wp = new Complex();
		//final Complex wXj = new Complex();

		// Now the FFT using the Cooley-Tukey algorithm...
		while(blockSize < data.length) {
			final int iStep = blockSize << 1;
			final double theta = sgPi / blockSize;

			wp.set(Math.cos(theta), Math.sin(theta));
			w.set(1.0, 0.0);

			for(int k=0; k<blockSize; k++) {
				for(int i=k; i<data.length; i+=iStep) {
					w.omegaFFT(data[i], data[i+blockSize]);
					
					/*
					final Complex di = data[i];
					final Complex dj = data[i+blockSize];
					
					// X_j = X_i - w X_j
					dj.multiplyBy(w);
					wXj.copy(dj);
					dj.isubtract(di);
					// X_i = X_i + w X_j
					di.add(wXj);
					*/
				}

				w.multiplyBy(wp);
			}
			blockSize <<= 1;
			Thread.yield();
		}
	} 
	
	/*
	public static void FFT(final Complex data[], boolean isForward) {
		bitReverse(data);
		
		final double sgPi = (isForward ? 1: -1) * Math.PI;
		int blockSize = 1;
		//final Complex wXj = new Complex();

		// Now the FFT using the Cooley-Tukey algorithm...
		while(blockSize < data.length) {
			final int iStep = blockSize << 1;
			final double theta = sgPi / blockSize;
			final double c = Math.cos(theta);
			final double s = Math.sin(theta);
			
			double wr = 1.0;
			double wi = 0.0;

			for(int k=0; k<blockSize; k++) {
				for(int i=k; i<data.length; i+=iStep) {
					final Complex c1 = data[i];
					final Complex c2 = data[i+blockSize];
					
					final double d2r = c2.x;
					final double d2i = c2.y;
					
					double x = wr * d2r - wi * d2i;
					c2.x = c1.x - x;
					c1.x += x;
					
					x = wr * d2i + wi * d2r;
					c2.y = c1.y - x;
					c1.y += x;
					
				}
				final double temp = wr;
				wr = temp * c - wi * s;
				wi = wi * c + temp * s;
			}
			blockSize <<= 1;
			Thread.yield();
		}
	} 
	*/
	
}
