package toolbox;
import java.io.PrintStream;
public class FFT1D
{
private boolean radix2;
private double Rearg[];
private double Imarg[];
private double yReOut[];
private double yImOut[];
private int maxPrimeFactor;
private int maxPrimeFactorDiv2;
private int maxFactorCount;
private int n;
private int nFactor;
private final double c3_1 = -1.5D;
private final double c3_2 = 0.86602540378444004D;
private final double c5_1 = -1.25D;
private final double c5_2 = 0.55901699437495D;
private final double c5_3 = -0.95105651629514998D;
private final double c5_4 = -1.5388417685875999D;
private final double c5_5 = 0.36327126400268001D;
private final double c8 = 0.70710678118655002D;
private int groupOffset;
private int dataOffset;
private int adr;
private int groupNo;
private int dataNo;
private int blockNo;
private int twNo;
private double omega;
private double tw_re;
private double tw_im;
private double twiddleRe[];
private double twiddleIm[];
private double trigRe[];
private double trigIm[];
private double zRe[];
private double zIm[];
private double vRe[];
private double vIm[];
private double wRe[];
private double wIm[];
private int sofarRadix[];
private int actualRadix[];
private int remainRadix[];
public FFT1D(int size)
{
radix2 = true;
maxFactorCount = 20;
int m = 1;
for(int size1 = size; size1 > 2;)
{
size1 /= 2;
m++;
}
if((int)Math.round(Math.pow(2D, m)) == size) //round()方法取最靠近的其的数值
{
radix2 = true;
n = 1 << m;
double fact = 6.2831853071795862D / (double)n;
Imarg = new double[n];
Rearg = new double[n];
for(int i = 0; i < n; i++)
{
double arg = fact * (double)i;
Rearg[i] = Math.cos(arg);
Imarg[i] = -Math.sin(arg);
}
} else
{
radix2 = false;
maxPrimeFactor = (int)(double)(size + 1);
maxPrimeFactorDiv2 = (maxPrimeFactor + 1) / 2;
twiddleRe = new double[maxPrimeFactor];
twiddleIm = new double[maxPrimeFactor];
trigRe = new double[maxPrimeFactor];
trigIm = new double[maxPrimeFactor];
zRe = new double[maxPrimeFactor];
zIm = new double[maxPrimeFactor];
vRe = new double[maxPrimeFactorDiv2];
vIm = new double[maxPrimeFactorDiv2];
wRe = new double[maxPrimeFactorDiv2];
wIm = new double[maxPrimeFactorDiv2];
yReOut = new double[size];
yImOut = new double[size];
n = size;
sofarRadix = new int[maxFactorCount];
actualRadix = new int[maxFactorCount];
remainRadix = new int[maxFactorCount];
transTableSetup(sofarRadix, actualRadix, remainRadix);
}
}
public final void transform(double Re[], double Im[], int size, int shift)
{
n = size;
if(radix2)
{
doFFT1D_CooleyTukey(Re, Im, size, shift);
} else
if(shift == 0)
{
doFFT_Mix(Re, Im, size);
} else
{
doFFT_Mix(Re, Im, size, shift);
}
}
public final void inverse(double Re[], double Im[], int size, int shift)
{
n = size;
if(radix2)
{
doIFFT1D_CooleyTukey(Re, Im, size, shift);
} else
if(shift == 0)
{
doIFFT_Mix(Re, Im, size);
} else
{
doIFFT_Mix(Re, Im, size, shift);
}
}
private void doFFT1D_CooleyTukey(double Re[], double Im[], int size, int shift)
{
int m = 1;
for(int size1 = size; size1 > 2;)
{
size1 /= 2;
m++;
}
int j;
for(int i = j = shift; i < (shift + n) - 1; i++)
{
if(i < j)
{
double Retmp = Re[i];
double Imtmp = Im[i];
Re[i] = Re[j];
Im[i] = Im[j];
Re[j] = Retmp;
Im[j] = Imtmp;
}
int k;
for(k = n >> 1; k + shift <= j; k /= 2)
{
j -= k;
}
j += k;
}
int stepsize = 1;
for(int shifter = m - 1; stepsize < n; shifter--)
{
for(int j = shift; j < shift + n; j += stepsize << 1)
{
for(int i = 0; i < stepsize; i++)
{
int i_j = i + j;
int i_j_s = i_j + stepsize;
double Retmp;
if(i > 0)
{
Retmp = Rearg[i << shifter] * Re[i_j_s] - Imarg[i << shifter] * Im[i_j_s];
Im[i_j_s] = Rearg[i << shifter] * Im[i_j_s] + Imarg[i << shifter] * Re[i_j_s];
Re[i_j_s] = Retmp;
}
Retmp = Re[i_j] - Re[i_j_s];
double Imtmp = Im[i_j] - Im[i_j_s];
Re[i_j] += Re[i_j_s];
Im[i_j] += Im[i_j_s];
Re[i_j_s] = Retmp;
Im[i_j_s] = Imtmp;
}
}
stepsize <<= 1;
}
}
private void doIFFT1D_CooleyTukey(double Re[], double Im[], int size, int shift)
{
for(int i = shift; i < shift + size; i++)
{
Im[i] = -Im[i];
}
transform(Re, Im, size, shift);
for(int i = shift; i < shift + size; i++)
{
Re[i] = Re[i] / (double)size;
Im[i] = -Im[i] / (double)size;
}
}
private void factorize(int fact[], int num)
{
int radices[] = new int[7];
int factors[] = new int[maxFactorCount];
int nRadix = 6;
radices[1] = 2;
radices[2] = 3;
radices[3] = 4;
radices[4] = 5;
radices[5] = 8;
radices[6] = 10;
int j;
if(num == 1)
{
j = 1;
factors[1] = 1;
} else
{
j = 0;
}
for(int i = nRadix; num > 1 && i > 0;)
{
if(num % radices[i] == 0)
{
num /= radices[i];
j++;
factors[j] = radices[i];
} else
{
i--;
}
}
if(factors[j] == 2)
{
int i;
for(i = j - 1; i > 0 && factors[i] != 8; i--) { }
if(i > 0)
{
factors[j] = 4;
factors[i] = 4;
}
}
if(num > 1)
{
for(int k = 2; (double)k < Math.sqrt(num) + 1.0D; k++)
{
while(num % k == 0)
{
num /= k;
j++;
factors[j] = k;
}
}
if(num > 1)
{
j++;
factors[j] = num;
}
}
for(int i = 1; i <= j; i++)
{
fact[i] = factors[(j - i) + 1];
}
nFactor = j;
}
private final void transTableSetup(int sofar[], int actual[], int remain[])
{
factorize(actual, n);
if(actual[1] > maxPrimeFactor)
{
System.out.println((new StringBuilder("\nPrime factor of FFT length too large : %6d")).append(actual[1]).toString());
