package jp.hasc.hasctool.core.runtime.filter.frequency;

import jp.hasc.hasctool.core.data.Complex;


/**
 * FFTクラスです
 * inverse:true(逆フーリエ変換)
 * inverse:false(フーリエ変換)
 * @author hiro
 */
public class FFT {

	private Complex[] input;
	private Complex[] output;
	private boolean inverse;
	private int length;

	public void setdata(Complex[] data){
		length = data.length;
		input = data;
	}

	public Complex[] getdata(){
		return output;
	}

	public FFT(boolean inv){
		inverse = inv;
	}

	public void execute() {
		output = new Complex[length];

		if(!(inverse)){
			fft();
		}else{
			ifft();
		}
	}

	//FFT
	private void fft(){
		int N = length;
		if (N == 1){
			output = new Complex[] { input[0] };
		}else if (N % 2 != 0){
			throw new RuntimeException("N is not a power of 2");
		}else{
			Complex[] even = new Complex[N/2];
			for (int k = 0; k < N/2; k++) {
				even[k] = input[2*k];
			}
			FFT fft = new FFT(false);
			fft.setdata(even);
			fft.execute();
			Complex[] q = fft.getdata();

			Complex[] odd  = even;
			for (int k = 0; k < N/2; k++) {
				odd[k] = input[2*k + 1];
			}

			fft = new FFT(false);
			fft.setdata(odd);
			fft.execute();
			Complex[] r = fft.getdata();

			for (int k = 0; k < N/2; k++) {
				double kth = -2 * k * Math.PI / N;
				Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
				output[k]       = q[k].plus(wk.times(r[k]));
				output[k + N/2] = q[k].minus(wk.times(r[k]));
			}
		}
	}

	//inverseFFT
	private void ifft(){
		int N = length;
		Complex[] val = new Complex[N];
		for (int i = 0; i < N; i++) {
			val[i] = input[i].conjugate();
		}

		FFT fft = new FFT(false);
		fft.setdata(val);
		fft.execute();
		val = fft.getdata();

		for (int i = 0; i < N; i++) {
			val[i] = val[i].conjugate();
		}
		for (int i = 0; i < N; i++) {
			output[i] = val[i].times(1.0 / N);
		}
	}

}
