/**Calculates the DFT using a split-radix Fast Fourier Transform algorithm.
* @author Junan Zhang
*/
public final class FFT {
public int logm;
final int MAXLOGM=20; /* max FFT length 2^MAXLOGM */
final double TWOPI=6.28318530717958647692;
final double SQHALF=0.707106781186547524401;
int brseed[]= new int[4048];
float tab[][];
public FFT(int nlength ) {
double dtemp = Math.log(nlength) / Math.log(2);
if ( (dtemp - (int) dtemp) != 0.0) {
throw new Error("FFT length must be a power of 2.");
} else {
this.logm = (int) dtemp;
}
if (logm>=4) {
creattab(logm);
}
}
/** Calculates the magnitude spectrum of a real signal.
* The returned vector contains only the positive frequencies.
*/
public float[] calculateFFTMagnitude(float x[]) {
int i,n;
n=1<<this.logm;
if (x.length > n) {
throw new Error("Tried to use a " + n + "-points FFT for a vector with " +
x.length + " samples!");
}
rsfft(x);
float[] mag = new float[n/2 + 1];
mag[0] = x[0]; //DC frequency must be positive always
if (n==1) {
return mag;
}
mag[n/2] = (float) Math.abs(x[n/2]); //pi (meaning: fs / 2)
//System.out.println("FFT before magnitude");
//IO.DisplayVector(x);
for (i=1;i<n/2;i++) {
mag[i] = (float) Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
//System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
}
//IO.DisplayVector(mag);
return mag;
}
/** Calculates the magnitude spectrum of a real signal.
* The returned vector contains only the positive frequencies.
*/
public double[] calculateFFTMagnitude(double inputData[]) {
int i,n;
n=1<<this.logm;
if (inputData.length > n) {
throw new Error("Tried to use a " + n + "-points FFT for a vector with " +
inputData.length + " samples!");
}
//System.out.println("magnitude test");
//double[] dtest = DSP.DFTMagnitude(inputData);
//IO.DisplayVector(dtest);
float[] x = new float[n];
for (i=0; i<inputData.length; i++) {
x[i] = (float) inputData[i];
}
rsfft(x);
//System.out.println("FFT before magnitude");
//IO.DisplayVector(x);
double[] mag = new double[n/2 + 1];
mag[0] = x[0]; //DC frequency must be positive always
if (n==1) {
return mag;
}
mag[n/2] = Math.abs(x[n/2]); //pi (meaning: fs / 2)
for (i=1;i<n/2;i++) {
mag[i] = Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
//System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
}
//IO.DisplayVector(mag);
return mag;
}
/** Calculates the power (magnitude squared) spectrum of a real signal.
* The returned vector contains only the positive frequencies.
*/
public double[] calculateFFTPower(double inputData[]) {
int i,n;
n=1<<this.logm;
//System.out.println("magnitude test");
//double[] dtest = DSP.DFTMagnitude(inputData);
//IO.DisplayVector(dtest);
float[] x = new float[n];
for (i=0; i<inputData.length; i++) {
x[i] = (float) inputData[i];
}
rsfft(x);
//System.out.println("FFT before magnitude");
//IO.DisplayVector(x);
double[] mag = new double[n/2 + 1];
mag[0] = x[0]; //DC frequency must be positive always
if (n==1) {
return mag;
}
mag[n/2] = Math.abs(x[n/2]); //pi (meaning: fs / 2)
for (i=1;i<n/2;i++) {
mag[i] = x[i]*x[i]+x[n-i]*x[n-i];
//mag[i] = Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
//System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
}
//IO.DisplayVector(mag);
return mag;
}
/**In place calculation of FFT magnitude.
*/
public void FFTMagnitude(float x[])
{ int i,n;
rsfft(x);
n=1<<this.logm;
if (n==1) return;
for (i=1;i<n/2;i++)
{x[i]=(float)Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
x[n-i]=x[i];
}
//Aldebaro modification:
x[n/2] = Math.abs(x[n/2]);
}
void rsfft(float x[])
{
/*creat table*/
// if(logm>=4) creattab(logm);
/* Call recursive routine */
rsrec(x,logm);
/* Output array unshuffling using bit-reversed indices */
if (logm > 1) {
BR_permute(x, logm);
return ;
}
}
/* -------------------------------------------------------------------- *
* Inverse transform for real inputs *
*-------------------------------------------------------------------- */
void rsifft(float x[])
{
int i, m;
float fac;
int xp;
/* Output array unshuffling using bit-reversed indices */
if (logm > 1) {
BR_permute(x, logm);
}
x[0] *= 0.5;
if (logm > 0) x[1] *= 0.5;
/*creat table*/
//if(logm>=4) creattab(logm);
/* Call recursive routine */
rsirec(x, logm);
/* Normalization */
m = 1 << logm;
fac = (float)2.0 / m;
xp = 0;
for (i = 0; i < m; i++) {
x[xp++] *= fac;
}
}
/* -------------------------------------------------------------------- *
* Creat multiple fator table *
* -------------------------------------------------------------------- */
void creattab(int logm)
{ int m, m2, m4, m8, nel, n, rlogm;
int cn, spcn, smcn, c3n, spc3n, smc3n;
double ang, s, c;
tab=new float [logm-4+1][6*((1<<logm)/4-2)];
for(rlogm=logm; rlogm>=4;rlogm--)
{m = 1 << rlogm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2; nel=m4-2;
/* Initialize pointers */
cn =0; spcn = cn + nel; smcn = spcn + nel;c3n = smcn + nel; spc3n = c3n + nel; smc3n = spc3n + nel;
/* Compute tables */
for (n = 1; n < m4; n++) {
if (n == m8) continue;
ang = n * TWOPI / m;
c = Math.cos(ang); s = Math.sin(ang);
tab[rlogm-4][cn++] = (float)c; tab[rlogm-4][spcn++] = (float)(- (s + c)); tab[rlogm-4][smcn++] =(float)( s - c);
ang = 3 * n * TWOPI / m;
c = Math.cos(ang); s = Math. sin(ang);
tab[rlogm-4][c3n++] = (float)c; tab[rlogm-4][spc3n++] = (float)(- (s + c)); tab[rlogm-4][smc3n++] = (float)(s - c);
}
}
}
/* -------------------------------------------------------------------- *
* Recursive part of the RSFFT algorithm. Not externally *
* callable. *
* -------------------------------------------------------------------- */
void rsrec(float x[],int logm)
{
int m, m2, m4, m8, nel, n;
int x0=0;
int xr1, xr2, xi1;
int cn=0;
int spcn=0;
int smcn=0;
float tmp1, tmp2;
double ang, c, s;
/* Check range of logm */
try{ if ((logm < 0) || (logm > MAXLOGM)) {
System.err.println("FFT length m is too big: log2(m) = "+logm+"is out of bounds ["+0+","+MAXLOGM+"]");
throw new OutofborderException(logm);
}}
catch( OutofborderException e)
{throw new OutOfMemoryError();}
/* Compute trivial cases */
if (logm < 2) {
if (logm == 1) { /* length m = 2 */
xr2 = x0 + 1;
tmp1 = x[x0] + x[xr2];
x[xr2] = x[x0] - x[xr2];
x[x0] = tmp1;
return;
}
else if (logm == 0) return; /* length m = 1 */
}
/* Compute a few constants */
m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;
/* Build tables of butterfly coefficients, if necessary */
//if ((logm >= 4) && (tab[logm-4][0] == 0)) {
/* Allocate memory for tables */
// nel = m4 - 2;
/*if ((tab[logm-4] = (float *) calloc(3 * nel, sizeof(float)))
== NULL) {
printf("Error : RSFFT : not enough memory for cosine tables.\n");
error_exit();
}*/
/* Initialize pointers */
//tabi=logm-4;
// cn =0; spcn = cn + nel; smcn = spcn + nel;
/* Compute tables */
/*for (n = 1; n < m4; n++) {
if (n == m8) continue;
ang = n * (float)TWOPI / m;
c = Math.cos(ang); s = Math.sin(ang);
tab[tabi][cn++] = (float)c; tab[tabi][spcn++] = (float)(- (s + c)); tab[tabi][smcn++] =(float)( s - c);
}
}
/* Step 1 */
xr1 = x
