Audio processing library
//GNU General Public License, version 2
package ca.uol.aig.fftpack;
/**
* Construct a 1-D complex data sequence.
*/
public class Complex1D
{
/**
* <em>x</em>[<em>i</em>] is the real part of <em>i</em>-th complex data.
*/
public double x[];
/**
* <em>y</em>[<em>i</em>] is the imaginary part of <em>i</em>-th complex data.
*/
public double y[];
}
package ca.uol.aig.fftpack;
/**
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
class RealDoubleFFT_Mixed
{
// ******************************************************************** //
// Real-Valued FFT Initialization.
// ******************************************************************** //
/**
* Initialization of Real FFT.
*/
void rffti(int n, double wtable[]) /* length of wtable = 2*n + 15 */
{
if (n == 1)
return;
rffti1(n, wtable, 0);
}
/*---------------------------------------------------------
rffti1: further initialization of Real FFT
--------------------------------------------------------*/
void rffti1(int n, double wtable[], int offset)
{
double argh;
int ntry=0, i, j;
double argld;
int k1, l1, l2, ib;
double fi;
int ld, ii, nf, ip, nl, is, nq, nr;
double arg;
int ido, ipm;
int nfm1;
// Create a working array.
tempData = new double[n];
nl=n;
nf=0;
j=0;
factorize_loop:
while(true)
{
++j;
if(j<=4)
ntry=NTRY_H[j-1];
else
ntry+=2;
do
{
nq=nl / ntry;
nr=nl-ntry*nq;
if(nr !=0) continue factorize_loop;
++nf;
wtable[nf+1+2*n+offset]=ntry;
nl=nq;
if(ntry==2 && nf !=1)
{
for(i=2; i<=nf; i++)
{
ib=nf-i+2;
wtable[ib+1+2*n+offset]=wtable[ib+2*n+offset];
}
wtable[2+2*n+offset]=2;
}
}while(nl !=1);
break factorize_loop;
}
wtable[0+2*n+offset] = n;
wtable[1+2*n+offset] = nf;
argh=TWO_PI /(double)(n);
is=0;
nfm1=nf-1;
l1=1;
if(nfm1==0) return;
for(k1=1; k1<=nfm1; k1++)
{
ip=(int)wtable[k1+1+2*n+offset];
ld=0;
l2=l1*ip;
ido=n / l2;
ipm=ip-1;
for(j=1; j<=ipm;++j)
{
ld+=l1;
i=is;
argld=(double)ld*argh;
fi=0;
for(ii=3; ii<=ido; ii+=2)
{
i+=2;
fi+=1;
arg=fi*argld;
wtable[i-2+n+offset] = Math.cos(arg);
wtable[i-1+n+offset] = Math.sin(arg);
}
is+=ido;
}
l1=l2;
}
} /*rffti1*/
// ******************************************************************** //
// Real-Valued FFT -- Forward Transform.
// ******************************************************************** //
/*---------------------------------------------------------
rfftf: Real forward FFT
--------------------------------------------------------*/
void rfftf(int n, double r[], double wtable[])
{
if(n==1) return;
rfftf1(n, r, wtable, 0);
} /*rfftf*/
/*---------------------------------------------------------
rfftf1: further processing of Real forward FFT
--------------------------------------------------------*/
void rfftf1(int n, double[] c, final double[] wtable, int offset)
{
final double[] td = tempData;
System.arraycopy(wtable, offset, td, 0, n);
int nf = (int) wtable[1 + 2 * n + offset];
int na = 1;
int l2 = n;
int iw = n - 1 + n + offset;
for (int k1 = 1; k1 <= nf; ++k1) {
int kh = nf - k1;
int ip = (int) wtable[kh + 2 + 2 * n + offset];
int l1 = l2 / ip;
int ido = n / l2;
int idl1 = ido * l1;
iw -= (ip - 1) * ido;
na = 1 - na;
if (ip == 4) {
if (na == 0)
radf4(ido, l1, c, td, wtable, iw);
else
radf4(ido, l1, td, c, wtable, iw);
} else if (ip == 2) {
if (na == 0)
radf2(ido, l1, c, td, wtable, iw);
else
radf2(ido, l1, td, c, wtable, iw);
} else if (ip == 3) {
if (na == 0)
radf3(ido, l1, c, td, wtable, iw);
else
radf3(ido, l1, td, c, wtable, iw);
} else if (ip == 5) {
if (na == 0)
radf5(ido, l1, c, td, wtable, iw);
else
radf5(ido, l1, td, c, wtable, iw);
} else {
if (ido == 1)
na = 1 - na;
if (na == 0) {
radfg(ido, ip, l1, idl1, c, c, c, td, td, wtable, iw);
na = 1;
} else {
radfg(ido, ip, l1, idl1, td, td, td, c, c, wtable, iw);
na = 0;
}
}
l2 = l1;
}
// If na == 1, the results are in c. Otherwise they're in tempData.
if (na == 0)
for (int i = 0; i < n; i++)
c[i] = td[i];
}
// ******************************************************************** //
// Real-Valued FFT -- Reverse Transform.
// ******************************************************************** //
/*---------------------------------------------------------
rfftb: Real backward FFT
--------------------------------------------------------*/
void rfftb(int n, double r[], double wtable[])
{
if(n==1) return;
rfftb1(n, r, wtable, 0);
} /*rfftb*/
/*---------------------------------------------------------
rfftb1: further processing of Real backward FFT
--------------------------------------------------------*/
void rfftb1(int n, double c[], final double wtable[], int offset)
{
int k1, l1, l2, na, nf, ip, iw, ido, idl1;
final double[] td = tempData;
System.arraycopy(wtable, offset, td, 0, n);
nf=(int)wtable[1+2*n+offset];
na=0;
l1=1;
iw=n+offset;
for(k1=1; k1<=nf; k1++)
{
ip=(int)wtable[k1+1+2*n+offset];
l2=ip*l1;
ido=n / l2;
idl1=ido*l1;
if(ip==4)
{
if(na==0)
{
radb4(ido, l1, c, td, wtable, iw);
}
else
{
radb4(ido, l1, td, c, wtable, iw);
}
na=1-na;
}
else if(ip==2)
{
if(na==0)
{
radb2(ido, l1, c, td, wtable, iw);
}
else
{
radb2(ido, l1, td, c, wtable, iw);
}
na=1-na;
}
else if(ip==3)
{
if(na==0)
{
radb3(ido, l1, c, td, wtable, iw);
}
else
{
radb3(ido, l1, td, c, wtable, iw);
}
na=1-na;
}
else if(ip==5)
{
if(na==0)
{
radb5(ido, l1, c, td, wtable, iw);
}
else
{
radb5(ido, l1, td, c, wtable, iw);
}
na=1-na;
}
else
{
if(na==0)
{
radbg(ido, ip, l1, idl1, c, c, c, td, td, wtable, iw);
}
else
{
radbg(ido, ip, l1, idl1, td, td, td, c, c, wtable, iw);
}
if(ido==1) na=1-na;
}
l1=l2;
iw+=(ip-1)*ido;
}
if (na == 1)
for (int i = 0; i < n; i++)
c[i] = td[i];
}
// ******************************************************************** //
// Real-Valued FFT -- General Subroutines.
// ******************************************************************** //
/*---------------------------------------------------------
radfg: Real FFT's forward processing of general factor
--------------------------------------------------------*/
private void radfg(int ido, int ip, int l1, int idl1, double cc[],
double c1[], double c2[], double ch[], double ch2[],
final double wtable[], int offset)
{
int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
double dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
int iw1 = offset;
arg=TWO_PI / (double)ip;
dcp=Math.cos(arg);
dsp=Math.sin(arg);
ipph=(ip+1)/ 2;
nbd=(ido-1)/ 2;
if(ido !=1)
{
for(ik=0; ik<idl1; ik++) ch2[ik]=c2[ik];
for(j=1; j<ip; j++)
for(k=0; k<l1; k++)
ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido];
if(nbd<=l1)
{
is=-ido;
for(j=1; j<ip; j++)
{
is+=ido;
idij=is-1;
for(i=2; i<ido; i+=2)
{
idij+=2;
for(k=0; k<l1; k++)
{
ch[i-1+(k+j*l1)*ido]=
wtable[idij-1+iw1]*c1[i-1+(k+j*l1)*ido]
+wtable[idij+iw1]*c1[i+(k+j*l1)*ido];
ch[i+(k+j*l1)*ido]=
wtable[idij-1+iw1]*c1[i+(k+j*l1)*ido]
-wtable[idij+iw1]*c1[i-1+(k+j*l1)*ido];
}
}
}
}
else
{
is=-ido;
for(j=1; j<ip; j++)
{
is+=ido;
for(k=0; k<l1; k++)
{
idij=is-1;
for(i=2; i<ido; i+=2)
{
idij+=2;
ch[i-1+(k+j*l1)*ido]=
wtable[idij-1+iw1]*c1[i-1+(k+j*l1)*ido]
+wtable[idij+iw1]*c1[i+(k+j*l1)*ido];
ch[i+(k+j*l1)*ido]=
wtable[idij-1+iw1]*c1[i+(k+j*l1)*ido]
-wtable[idij+iw1]*c1[i-1+(k+j*l1)*ido];
}
}
}
}
if(nbd>=l1)
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
c1[i-1+(k+j*l1)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
}
}
}
}
else
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(i=2; i<ido; i+=2)
{
for(k=0; k<l1; k++)
{
c1[i-1+(k+j*l1)*ido]=
ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
c1[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
c1[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
c1[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
}
}
}
}
}
else
{
for(ik=0; ik<idl1; ik++) c2[ik]=ch2[ik];
}
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido]+ch[(k+jc*l1)*ido];
c1[(k+jc*l1)*ido]=ch[(k+jc*l1)*ido]-ch[(k+j*l1)*ido];
}
}
ar1=1;
ai1=0;
for(l=1; l<ipph; l++)
{
lc=ip-l;
ar1h=dcp*ar1-dsp*ai1;
ai1=dcp*ai1+dsp*ar1;
ar1=ar1h;
for(ik=0; ik<idl1; ik++)
{
ch2[ik+l*idl1]=c2[ik]+ar1*c2[ik+idl1];
ch2[ik+lc*idl1]=ai1*c2[ik+(ip-1)*idl1];
}
dc2=ar1;
ds2=ai1;
ar2=ar1;
ai2=ai1;
for(j=2; j<ipph; j++)
{
jc=ip-j;
ar2h=dc2*ar2-ds2*ai2;
ai2=dc2*ai2+ds2*ar2;
ar2=ar2h;
for(ik=0; ik<idl1; ik++)
{
ch2[ik+l*idl1]+=ar2*c2[ik+j*idl1];
ch2[ik+lc*idl1]+=ai2*c2[ik+jc*idl1];
}
}
}
for(j=1; j<ipph; j++)
for(ik=0; ik<idl1; ik++)
ch2[ik]+=c2[ik+j*idl1];
if(ido>=l1)
{
for(k=0; k<l1; k++)
{
for(i=0; i<ido; i++)
{
cc[i+k*ip*ido]=ch[i+k*ido];
}
}
}
else
{
for(i=0; i<ido; i++)
{
for(k=0; k<l1; k++)
{
cc[i+k*ip*ido]=ch[i+k*ido];
}
}
}
for(j=1; j<ipph; j++)
{
jc=ip-j;
j2=2*j;
for(k=0; k<l1; k++)
{
cc[ido-1+(j2-1+k*ip)*ido]=ch[(k+j*l1)*ido];
cc[(j2+k*ip)*ido]=ch[(k+jc*l1)*ido];
}
}
if(ido==1) return;
if(nbd>=l1)
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
j2=2*j;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
}
}
}
}
else
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
j2=2*j;
for(i=2; i<ido; i+=2)
{
ic=ido-i;
for(k=0; k<l1; k++)
{
cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
}
}
}
}
}
/*---------------------------------------------------------
radbg: Real FFT's backward processing of general factor
--------------------------------------------------------*/
private void radbg(int ido, int ip, int l1, int idl1, double cc[], double c1[],
double c2[], double ch[], double ch2[], final double wtable[], int offset)
{
int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
double dc2, ai1, ai2, ar1, ar2, ds2;
int nbd;
double dcp, arg, dsp, ar1h, ar2h;
int iw1 = offset;
arg=TWO_PI / (double)ip;
dcp=Math.cos(arg);
dsp=Math.sin(arg);
nbd=(ido-1)/ 2;
ipph=(ip+1)/ 2;
if(ido>=l1)
{
for(k=0; k<l1; k++)
{
for(i=0; i<ido; i++)
{
ch[i+k*ido]=cc[i+k*ip*ido];
}
}
}
else
{
for(i=0; i<ido; i++)
{
for(k=0; k<l1; k++)
{
ch[i+k*ido]=cc[i+k*ip*ido];
}
}
}
for(j=1; j<ipph; j++)
{
jc=ip-j;
j2=2*j;
for(k=0; k<l1; k++)
{
ch[(k+j*l1)*ido]=cc[ido-1+(j2-1+k*ip)*ido]+cc[ido-1+(j2-1+k*ip)*ido];
ch[(k+jc*l1)*ido]=cc[(j2+k*ip)*ido]+cc[(j2+k*ip)*ido];
}
}
if(ido !=1)
{
if(nbd>=l1)
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
}
}
}
}
else
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(i=2; i<ido; i+=2)
{
ic=ido-i;
for(k=0; k<l1; k++)
{
ch[i-1+(k+j*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
ch[i-1+(k+jc*l1)*ido]=cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
}
}
}
}
}
ar1=1;
ai1=0;
for(l=1; l<ipph; l++)
{
lc=ip-l;
ar1h=dcp*ar1-dsp*ai1;
ai1=dcp*ai1+dsp*ar1;
ar1=ar1h;
for(ik=0; ik<idl1; ik++)
{
c2[ik+l*idl1]=ch2[ik]+ar1*ch2[ik+idl1];
c2[ik+lc*idl1]=ai1*ch2[ik+(ip-1)*idl1];
}
dc2=ar1;
ds2=ai1;
ar2=ar1;
ai2=ai1;
for(j=2; j<ipph; j++)
{
jc=ip-j;
ar2h=dc2*ar2-ds2*ai2;
ai2=dc2*ai2+ds2*ar2;
ar2=ar2h;
for(ik=0; ik<idl1; ik++)
{
c2[ik+l*idl1]+=ar2*ch2[ik+j*idl1];
c2[ik+lc*idl1]+=ai2*ch2[ik+jc*idl1];
}
}
}
for(j=1; j<ipph; j++)
{
for(ik=0; ik<idl1; ik++)
{
ch2[ik]+=ch2[ik+j*idl1];
}
}
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
ch[(k+j*l1)*ido]=c1[(k+j*l1)*ido]-c1[(k+jc*l1)*ido];
ch[(k+jc*l1)*ido]=c1[(k+j*l1)*ido]+c1[(k+jc*l1)*ido];
}
}
if(ido==1) return;
if(nbd>=l1)
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
}
}
}
}
else
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(i=2; i<ido; i+=2)
{
for(k=0; k<l1; k++)
{
ch[i-1+(k+j*l1)*ido]=c1[i-1+(k+j*l1)*ido]-c1[i+(k+jc*l1)*ido];
ch[i-1+(k+jc*l1)*ido]=c1[i-1+(k+j*l1)*ido]+c1[i+(k+jc*l1)*ido];
ch[i+(k+j*l1)*ido]=c1[i+(k+j*l1)*ido]+c1[i-1+(k+jc*l1)*ido];
ch[i+(k+jc*l1)*ido]=c1[i+(k+j*l1)*ido]-c1[i-1+(k+jc*l1)*ido];
}
}
}
}
for(ik=0; ik<idl1; ik++) c2[ik]=ch2[ik];
for(j=1; j<ip; j++)
for(k=0; k<l1; k++)
c1[(k+j*l1)*ido]=ch[(k+j*l1)*ido];
if(nbd<=l1)
{
is=-ido;
for(j=1; j<ip; j++)
{
is+=ido;
idij=is-1;
for(i=2; i<ido; i+=2)
{
idij+=2;
for(k=0; k<l1; k++)
{
c1[i-1+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido]
-wtable[idij+iw1]*ch[i+(k+j*l1)*ido];
c1[i+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido]
+wtable[idij+iw1]*ch[i-1+(k+j*l1)*ido];
}
}
}
}
else
{
is=-ido;
for(j=1; j<ip; j++)
{
is+=ido;
for(k=0; k<l1; k++)
{
idij=is-1;
for(i=2; i<ido; i+=2)
{
idij+=2;
c1[i-1+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido]
-wtable[idij+iw1]*ch[i+(k+j*l1)*ido];
c1[i+(k+j*l1)*ido] = wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido]
+wtable[idij+iw1]*ch[i-1+(k+j*l1)*ido];
}
}
}
}
}
// ******************************************************************** //
// Real-Valued FFT -- Factor-Specific Optimized Subroutines.
// ******************************************************************** //
/*-------------------------------------------------
radf2: Real FFT's forward processing of factor 2
-------------------------------------------------*/
private void radf2(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
int i, k, ic;
double ti2, tr2;
int iw1;
iw1 = offset;
for(k=0; k<l1; k++)
{
ch[2*k*ido]=cc[k*ido]+cc[(k+l1)*ido];
ch[(2*k+1)*ido+ido-1]=cc[k*ido]-cc[(k+l1)*ido];
}
if(ido<2) return;
if(ido !=2)
{
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
tr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
+wtable[i-1+iw1]*cc[i+(k+l1)*ido];
ti2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
-wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
ch[i+2*k*ido]=cc[i+k*ido]+ti2;
ch[ic+(2*k+1)*ido]=ti2-cc[i+k*ido];
ch[i-1+2*k*ido]=cc[i-1+k*ido]+tr2;
ch[ic-1+(2*k+1)*ido]=cc[i-1+k*ido]-tr2;
}
}
if(ido%2==1)return;
}
for(k=0; k<l1; k++)
{
ch[(2*k+1)*ido]=-cc[ido-1+(k+l1)*ido];
ch[ido-1+2*k*ido]=cc[ido-1+k*ido];
}
}
/*-------------------------------------------------
radb2: Real FFT's backward processing of factor 2
-------------------------------------------------*/
private void radb2(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
int i, k, ic;
double ti2, tr2;
int iw1 = offset;
for(k=0; k<l1; k++)
{
ch[k*ido]=cc[2*k*ido]+cc[ido-1+(2*k+1)*ido];
ch[(k+l1)*ido]=cc[2*k*ido]-cc[ido-1+(2*k+1)*ido];
}
if(ido<2) return;
if(ido !=2)
{
for(k=0; k<l1;++k)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
ch[i-1+k*ido]=cc[i-1+2*k*ido]+cc[ic-1+(2*k+1)*ido];
tr2=cc[i-1+2*k*ido]-cc[ic-1+(2*k+1)*ido];
ch[i+k*ido]=cc[i+2*k*ido]-cc[ic+(2*k+1)*ido];
ti2=cc[i+(2*k)*ido]+cc[ic+(2*k+1)*ido];
ch[i-1+(k+l1)*ido]=wtable[i-2+iw1]*tr2-wtable[i-1+iw1]*ti2;
ch[i+(k+l1)*ido]=wtable[i-2+iw1]*ti2+wtable[i-1+iw1]*tr2;
}
}
if(ido%2==1) return;
}
for(k=0; k<l1; k++)
{
ch[ido-1+k*ido]=2*cc[ido-1+2*k*ido];
ch[ido-1+(k+l1)*ido]=-2*cc[(2*k+1)*ido];
}
}
/*-------------------------------------------------
radf3: Real FFT's forward processing of factor 3
-------------------------------------------------*/
private void radf3(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
int i, k, ic;
double ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
int iw1, iw2;
iw1 = offset;
iw2 = iw1 + ido;
for(k=0; k<l1; k++)
{
cr2=cc[(k+l1)*ido]+cc[(k+2*l1)*ido];
ch[3*k*ido]=cc[k*ido]+cr2;
ch[(3*k+2)*ido]=TAU_I*(cc[(k+l1*2)*ido]-cc[(k+l1)*ido]);
ch[ido-1+(3*k+1)*ido]=cc[k*ido]+TAU_R*cr2;
}
if(ido==1) return;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
dr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
+wtable[i-1+iw1]*cc[i+(k+l1)*ido];
di2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
-wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
dr3 = wtable[i-2+iw2]*cc[i-1+(k+l1*2)*ido]
+wtable[i-1+iw2]*cc[i+(k+l1*2)*ido];
di3 = wtable[i-2+iw2]*cc[i+(k+l1*2)*ido]
-wtable[i-1+iw2]*cc[i-1+(k+l1*2)*ido];
cr2 = dr2+dr3;
ci2 = di2+di3;
ch[i-1+3*k*ido]=cc[i-1+k*ido]+cr2;
ch[i+3*k*ido]=cc[i+k*ido]+ci2;
tr2=cc[i-1+k*ido]+TAU_R*cr2;
ti2=cc[i+k*ido]+TAU_R*ci2;
tr3=TAU_I*(di2-di3);
ti3=TAU_I*(dr3-dr2);
ch[i-1+(3*k+2)*ido]=tr2+tr3;
ch[ic-1+(3*k+1)*ido]=tr2-tr3;
ch[i+(3*k+2)*ido]=ti2+ti3;
ch[ic+(3*k+1)*ido]=ti3-ti2;
}
}
}
/*-------------------------------------------------
radb3: Real FFT's backward processing of factor 3
-------------------------------------------------*/
private void radb3(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
int i, k, ic;
double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
int iw1, iw2;
iw1 = offset;
iw2 = iw1 + ido;
for(k=0; k<l1; k++)
{
tr2=2*cc[ido-1+(3*k+1)*ido];
cr2=cc[3*k*ido]+TAU_R*tr2;
ch[k*ido]=cc[3*k*ido]+tr2;
ci3=2*TAU_I*cc[(3*k+2)*ido];
ch[(k+l1)*ido]=cr2-ci3;
ch[(k+2*l1)*ido]=cr2+ci3;
}
if(ido==1) return;
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
tr2=cc[i-1+(3*k+2)*ido]+cc[ic-1+(3*k+1)*ido];
cr2=cc[i-1+3*k*ido]+TAU_R*tr2;
ch[i-1+k*ido]=cc[i-1+3*k*ido]+tr2;
ti2=cc[i+(3*k+2)*ido]-cc[ic+(3*k+1)*ido];
ci2=cc[i+3*k*ido]+TAU_R*ti2;
ch[i+k*ido]=cc[i+3*k*ido]+ti2;
cr3=TAU_I*(cc[i-1+(3*k+2)*ido]-cc[ic-1+(3*k+1)*ido]);
ci3=TAU_I*(cc[i+(3*k+2)*ido]+cc[ic+(3*k+1)*ido]);
dr2=cr2-ci3;
dr3=cr2+ci3;
di2=ci2+cr3;
di3=ci2-cr3;
ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*dr2
-wtable[i-1+iw1]*di2;
ch[i+(k+l1)*ido] = wtable[i-2+iw1]*di2
+wtable[i-1+iw1]*dr2;
ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*dr3
-wtable[i-1+iw2]*di3;
ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*di3
+wtable[i-1+iw2]*dr3;
}
}
}
/*-------------------------------------------------
radf4: Real FFT's forward processing of factor 4
-------------------------------------------------*/
private void radf4(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
final double hsqt2=0.7071067811865475D;
int i, k, ic;
double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
int iw1, iw2, iw3;
iw1 = offset;
iw2 = offset + ido;
iw3 = iw2 + ido;
for(k=0; k<l1; k++)
{
tr1=cc[(k+l1)*ido]+cc[(k+3*l1)*ido];
tr2=cc[k*ido]+cc[(k+2*l1)*ido];
ch[4*k*ido]=tr1+tr2;
ch[ido-1+(4*k+3)*ido]=tr2-tr1;
ch[ido-1+(4*k+1)*ido]=cc[k*ido]-cc[(k+2*l1)*ido];
ch[(4*k+2)*ido]=cc[(k+3*l1)*ido]-cc[(k+l1)*ido];
}
if(ido<2) return;
if(ido !=2)
{
for(k=0; k<l1; k++)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
cr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
+wtable[i-1+iw1]*cc[i+(k+l1)*ido];
ci2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
-wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
cr3 = wtable[i-2+iw2]*cc[i-1+(k+2*l1)*ido]
+wtable[i-1+iw2]*cc[i+(k+2*l1)*ido];
ci3 = wtable[i-2+iw2]*cc[i+(k+2*l1)*ido]
-wtable[i-1+iw2]*cc[i-1+(k+2*l1)*ido];
cr4 = wtable[i-2+iw3]*cc[i-1+(k+3*l1)*ido]
+wtable[i-1+iw3]*cc[i+(k+3*l1)*ido];
ci4 = wtable[i-2+iw3]*cc[i+(k+3*l1)*ido]
-wtable[i-1+iw3]*cc[i-1+(k+3*l1)*ido];
tr1=cr2+cr4;
tr4=cr4-cr2;
ti1=ci2+ci4;
ti4=ci2-ci4;
ti2=cc[i+k*ido]+ci3;
ti3=cc[i+k*ido]-ci3;
tr2=cc[i-1+k*ido]+cr3;
tr3=cc[i-1+k*ido]-cr3;
ch[i-1+4*k*ido]=tr1+tr2;
ch[ic-1+(4*k+3)*ido]=tr2-tr1;
ch[i+4*k*ido]=ti1+ti2;
ch[ic+(4*k+3)*ido]=ti1-ti2;
ch[i-1+(4*k+2)*ido]=ti4+tr3;
ch[ic-1+(4*k+1)*ido]=tr3-ti4;
ch[i+(4*k+2)*ido]=tr4+ti3;
ch[ic+(4*k+1)*ido]=tr4-ti3;
}
}
if(ido%2==1) return;
}
for(k=0; k<l1; k++)
{
ti1=-hsqt2*(cc[ido-1+(k+l1)*ido]+cc[ido-1+(k+3*l1)*ido]);
tr1=hsqt2*(cc[ido-1+(k+l1)*ido]-cc[ido-1+(k+3*l1)*ido]);
ch[ido-1+4*k*ido]=tr1+cc[ido-1+k*ido];
ch[ido-1+(4*k+2)*ido]=cc[ido-1+k*ido]-tr1;
ch[(4*k+1)*ido]=ti1-cc[ido-1+(k+2*l1)*ido];
ch[(4*k+3)*ido]=ti1+cc[ido-1+(k+2*l1)*ido];
}
}
/*-------------------------------------------------
radb4: Real FFT's backward processing of factor 4
-------------------------------------------------*/
private void radb4(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
int i, k, ic;
double ci2, ci3, ci4, cr2, cr3, cr4;
double ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
int iw1, iw2, iw3;
iw1 = offset;
iw2 = iw1 + ido;
iw3 = iw2 + ido;
for(k=0; k<l1; k++)
{
tr1=cc[4*k*ido]-cc[ido-1+(4*k+3)*ido];
tr2=cc[4*k*ido]+cc[ido-1+(4*k+3)*ido];
tr3=cc[ido-1+(4*k+1)*ido]+cc[ido-1+(4*k+1)*ido];
tr4=cc[(4*k+2)*ido]+cc[(4*k+2)*ido];
ch[k*ido]=tr2+tr3;
ch[(k+l1)*ido]=tr1-tr4;
ch[(k+2*l1)*ido]=tr2-tr3;
ch[(k+3*l1)*ido]=tr1+tr4;
}
if(ido<2) return;
if(ido !=2)
{
for(k=0; k<l1;++k)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
ti1=cc[i+4*k*ido]+cc[ic+(4*k+3)*ido];
ti2=cc[i+4*k*ido]-cc[ic+(4*k+3)*ido];
ti3=cc[i+(4*k+2)*ido]-cc[ic+(4*k+1)*ido];
tr4=cc[i+(4*k+2)*ido]+cc[ic+(4*k+1)*ido];
tr1=cc[i-1+4*k*ido]-cc[ic-1+(4*k+3)*ido];
tr2=cc[i-1+4*k*ido]+cc[ic-1+(4*k+3)*ido];
ti4=cc[i-1+(4*k+2)*ido]-cc[ic-1+(4*k+1)*ido];
tr3=cc[i-1+(4*k+2)*ido]+cc[ic-1+(4*k+1)*ido];
ch[i-1+k*ido]=tr2+tr3;
cr3=tr2-tr3;
ch[i+k*ido]=ti2+ti3;
ci3=ti2-ti3;
cr2=tr1-tr4;
cr4=tr1+tr4;
ci2=ti1+ti4;
ci4=ti1-ti4;
ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*cr2
-wtable[i-1+iw1]*ci2;
ch[i+(k+l1)*ido] = wtable[i-2+iw1]*ci2
+wtable[i-1+iw1]*cr2;
ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*cr3
-wtable[i-1+iw2]*ci3;
ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*ci3
+wtable[i-1+iw2]*cr3;
ch[i-1+(k+3*l1)*ido] = wtable[i-2+iw3]*cr4
-wtable[i-1+iw3]*ci4;
ch[i+(k+3*l1)*ido] = wtable[i-2+iw3]*ci4
+wtable[i-1+iw3]*cr4;
}
}
if(ido%2==1) return;
}
for(k=0; k<l1; k++)
{
ti1=cc[(4*k+1)*ido]+cc[(4*k+3)*ido];
ti2=cc[(4*k+3)*ido]-cc[(4*k+1)*ido];
tr1=cc[ido-1+4*k*ido]-cc[ido-1+(4*k+2)*ido];
tr2=cc[ido-1+4*k*ido]+cc[ido-1+(4*k+2)*ido];
ch[ido-1+k*ido]=tr2+tr2;
ch[ido-1+(k+l1)*ido]=SQRT_2*(tr1-ti1);
ch[ido-1+(k+2*l1)*ido]=ti2+ti2;
ch[ido-1+(k+3*l1)*ido]=-SQRT_2*(tr1+ti1);
}
}
/*-------------------------------------------------
radf5: Real FFT's forward processing of factor 5
-------------------------------------------------*/
private void radf5(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
final double tr11=0.309016994374947D;
final double ti11=0.951056516295154D;
final double tr12=-0.809016994374947D;
final double ti12=0.587785252292473D;
int i, k, ic;
double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
int iw1, iw2, iw3, iw4;
iw1 = offset;
iw2 = iw1 + ido;
iw3 = iw2 + ido;
iw4 = iw3 + ido;
for(k=0; k<l1; k++)
{
cr2=cc[(k+4*l1)*ido]+cc[(k+l1)*ido];
ci5=cc[(k+4*l1)*ido]-cc[(k+l1)*ido];
cr3=cc[(k+3*l1)*ido]+cc[(k+2*l1)*ido];
ci4=cc[(k+3*l1)*ido]-cc[(k+2*l1)*ido];
ch[5*k*ido]=cc[k*ido]+cr2+cr3;
ch[ido-1+(5*k+1)*ido]=cc[k*ido]+tr11*cr2+tr12*cr3;
ch[(5*k+2)*ido]=ti11*ci5+ti12*ci4;
ch[ido-1+(5*k+3)*ido]=cc[k*ido]+tr12*cr2+tr11*cr3;
ch[(5*k+4)*ido]=ti12*ci5-ti11*ci4;
}
if(ido==1) return;
for(k=0; k<l1;++k)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
dr2 = wtable[i-2+iw1]*cc[i-1+(k+l1)*ido]
+wtable[i-1+iw1]*cc[i+(k+l1)*ido];
di2 = wtable[i-2+iw1]*cc[i+(k+l1)*ido]
-wtable[i-1+iw1]*cc[i-1+(k+l1)*ido];
dr3 = wtable[i-2+iw2]*cc[i-1+(k+2*l1)*ido]
+wtable[i-1+iw2]*cc[i+(k+2*l1)*ido];
di3 = wtable[i-2+iw2]*cc[i+(k+2*l1)*ido]
-wtable[i-1+iw2]*cc[i-1+(k+2*l1)*ido];
dr4 = wtable[i-2+iw3]*cc[i-1+(k+3*l1)*ido]
+wtable[i-1+iw3]*cc[i+(k+3*l1)*ido];
di4 = wtable[i-2+iw3]*cc[i+(k+3*l1)*ido]
-wtable[i-1+iw3]*cc[i-1+(k+3*l1)*ido];
dr5 = wtable[i-2+iw4]*cc[i-1+(k+4*l1)*ido]
+wtable[i-1+iw4]*cc[i+(k+4*l1)*ido];
di5 = wtable[i-2+iw4]*cc[i+(k+4*l1)*ido]
-wtable[i-1+iw4]*cc[i-1+(k+4*l1)*ido];
cr2=dr2+dr5;
ci5=dr5-dr2;
cr5=di2-di5;
ci2=di2+di5;
cr3=dr3+dr4;
ci4=dr4-dr3;
cr4=di3-di4;
ci3=di3+di4;
ch[i-1+5*k*ido]=cc[i-1+k*ido]+cr2+cr3;
ch[i+5*k*ido]=cc[i+k*ido]+ci2+ci3;
tr2=cc[i-1+k*ido]+tr11*cr2+tr12*cr3;
ti2=cc[i+k*ido]+tr11*ci2+tr12*ci3;
tr3=cc[i-1+k*ido]+tr12*cr2+tr11*cr3;
ti3=cc[i+k*ido]+tr12*ci2+tr11*ci3;
tr5=ti11*cr5+ti12*cr4;
ti5=ti11*ci5+ti12*ci4;
tr4=ti12*cr5-ti11*cr4;
ti4=ti12*ci5-ti11*ci4;
ch[i-1+(5*k+2)*ido]=tr2+tr5;
ch[ic-1+(5*k+1)*ido]=tr2-tr5;
ch[i+(5*k+2)*ido]=ti2+ti5;
ch[ic+(5*k+1)*ido]=ti5-ti2;
ch[i-1+(5*k+4)*ido]=tr3+tr4;
ch[ic-1+(5*k+3)*ido]=tr3-tr4;
ch[i+(5*k+4)*ido]=ti3+ti4;
ch[ic+(5*k+3)*ido]=ti4-ti3;
}
}
}
/*-------------------------------------------------
radb5: Real FFT's backward processing of factor 5
-------------------------------------------------*/
private void radb5(int ido, int l1, final double cc[], double ch[],
final double wtable[], int offset)
{
final double tr11=0.309016994374947D;
final double ti11=0.951056516295154D;
final double tr12=-0.809016994374947D;
final double ti12=0.587785252292473D;
int i, k, ic;
double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
int iw1, iw2, iw3, iw4;
iw1 = offset;
iw2 = iw1 + ido;
iw3 = iw2 + ido;
iw4 = iw3 + ido;
for(k=0; k<l1; k++)
{
ti5=2*cc[(5*k+2)*ido];
ti4=2*cc[(5*k+4)*ido];
tr2=2*cc[ido-1+(5*k+1)*ido];
tr3=2*cc[ido-1+(5*k+3)*ido];
ch[k*ido]=cc[5*k*ido]+tr2+tr3;
cr2=cc[5*k*ido]+tr11*tr2+tr12*tr3;
cr3=cc[5*k*ido]+tr12*tr2+tr11*tr3;
ci5=ti11*ti5+ti12*ti4;
ci4=ti12*ti5-ti11*ti4;
ch[(k+l1)*ido]=cr2-ci5;
ch[(k+2*l1)*ido]=cr3-ci4;
ch[(k+3*l1)*ido]=cr3+ci4;
ch[(k+4*l1)*ido]=cr2+ci5;
}
if(ido==1) return;
for(k=0; k<l1;++k)
{
for(i=2; i<ido; i+=2)
{
ic=ido-i;
ti5=cc[i+(5*k+2)*ido]+cc[ic+(5*k+1)*ido];
ti2=cc[i+(5*k+2)*ido]-cc[ic+(5*k+1)*ido];
ti4=cc[i+(5*k+4)*ido]+cc[ic+(5*k+3)*ido];
ti3=cc[i+(5*k+4)*ido]-cc[ic+(5*k+3)*ido];
tr5=cc[i-1+(5*k+2)*ido]-cc[ic-1+(5*k+1)*ido];
tr2=cc[i-1+(5*k+2)*ido]+cc[ic-1+(5*k+1)*ido];
tr4=cc[i-1+(5*k+4)*ido]-cc[ic-1+(5*k+3)*ido];
tr3=cc[i-1+(5*k+4)*ido]+cc[ic-1+(5*k+3)*ido];
ch[i-1+k*ido]=cc[i-1+5*k*ido]+tr2+tr3;
ch[i+k*ido]=cc[i+5*k*ido]+ti2+ti3;
cr2=cc[i-1+5*k*ido]+tr11*tr2+tr12*tr3;
ci2=cc[i+5*k*ido]+tr11*ti2+tr12*ti3;
cr3=cc[i-1+5*k*ido]+tr12*tr2+tr11*tr3;
ci3=cc[i+5*k*ido]+tr12*ti2+tr11*ti3;
cr5=ti11*tr5+ti12*tr4;
ci5=ti11*ti5+ti12*ti4;
cr4=ti12*tr5-ti11*tr4;
ci4=ti12*ti5-ti11*ti4;
dr3=cr3-ci4;
dr4=cr3+ci4;
di3=ci3+cr4;
di4=ci3-cr4;
dr5=cr2+ci5;
dr2=cr2-ci5;
di5=ci2-cr5;
di2=ci2+cr5;
ch[i-1+(k+l1)*ido] = wtable[i-2+iw1]*dr2
-wtable[i-1+iw1]*di2;
ch[i+(k+l1)*ido] = wtable[i-2+iw1]*di2
+wtable[i-1+iw1]*dr2;
ch[i-1+(k+2*l1)*ido] = wtable[i-2+iw2]*dr3
-wtable[i-1+iw2]*di3;
ch[i+(k+2*l1)*ido] = wtable[i-2+iw2]*di3
+wtable[i-1+iw2]*dr3;
ch[i-1+(k+3*l1)*ido] = wtable[i-2+iw3]*dr4
-wtable[i-1+iw3]*di4;
ch[i+(k+3*l1)*ido] = wtable[i-2+iw3]*di4
+wtable[i-1+iw3]*dr4;
ch[i-1+(k+4*l1)*ido] = wtable[i-2+iw4]*dr5
-wtable[i-1+iw4]*di5;
ch[i+(k+4*l1)*ido] = wtable[i-2+iw4]*di5
+wtable[i-1+iw4]*dr5;
}
}
}
// ******************************************************************** //
// Private Constants.
// ******************************************************************** //
private static final int[] NTRY_H = { 4, 2, 3, 5 };
private static final double TWO_PI = 2.0 * Math.PI;
private static final double SQRT_2 = 1.414213562373095;
private static final double TAU_R = -0.5;
private static final double TAU_I = 0.866025403784439;
// ******************************************************************** //
// Private Data.
// ******************************************************************** //
// Working data array for FFT.
private double[] tempData = null;
}
package ca.uol.aig.fftpack;
/**
* sine FFT transform of a real odd sequence.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class RealDoubleFFT_Odd extends RealDoubleFFT_Mixed
{
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward
* transform (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
public double norm_factor;
private double wavetable[];
private int ndim;
/**
* Construct a wavenumber table with size <em>n</em>.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a real data sequence. When (<em>n</em>+1) is a multiplication of small
* numbers (4, 2, 3, 5), this FFT transform is very efficient.
*/
public RealDoubleFFT_Odd(int n)
{
ndim = n;
norm_factor = 2*(n+1);
int wtable_length = 2*ndim + ndim/2 + 3 + 15;
if(wavetable == null || wavetable.length != wtable_length)
{
wavetable = new double[wtable_length];
}
sinti(ndim, wavetable);
}
/**
* Forward sine FFT transform. It computes the discrete sine transform of
* an odd sequence.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
public void ft(double x[])
{
sint(ndim, x, wavetable);
}
/**
* Backward sine FFT transform. It is the unnormalized inverse transform of <em>ft</em>.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
public void bt(double x[])
{
sint(ndim, x, wavetable);
}
/*---------------------------------------
sint1: further processing of sine FFT.
--------------------------------------*/
void sint1(int n, double war[], double wtable[])
{
final double sqrt3=1.73205080756888;
int modn, i, k;
double xhold, t1, t2;
int kc, np1, ns2;
int iw1, iw2, iw3;
double[] wtable_p1 = new double[2*(n+1)+15];
iw1=n / 2;
iw2=iw1+n+1;
iw3=iw2+n+1;
double[] x = new double[n+1];
for(i=0; i<n; i++)
{
wtable[i+iw1]=war[i];
war[i]=wtable[i+iw2];
}
if(n<2)
{
wtable[0+iw1]+=wtable[0+iw1];
}
else if(n==2)
{
xhold=sqrt3*(wtable[0+iw1]+wtable[1+iw1]);
wtable[1+iw1]=sqrt3*(wtable[0+iw1]-wtable[1+iw1]);
wtable[0+iw1]=xhold;
}
else
{
np1=n+1;
ns2=n / 2;
wtable[0+iw2]=0;
for(k=0; k<ns2; k++)
{
kc=n-k-1;
t1=wtable[k+iw1]-wtable[kc+iw1];
t2=wtable[k]*(wtable[k+iw1]+wtable[kc+iw1]);
wtable[k+1+iw2]=t1+t2;
wtable[kc+1+iw2]=t2-t1;
}
modn=n%2;
if(modn !=0)
wtable[ns2+1+iw2]=4*wtable[ns2+iw1];
System.arraycopy(wtable, iw1, wtable_p1, 0, n+1);
System.arraycopy(war, 0, wtable_p1, n+1, n);
System.arraycopy(wtable, iw3, wtable_p1, 2*(n+1), 15);
System.arraycopy(wtable, iw2, x, 0, n+1);
rfftf1(np1, x, wtable_p1, 0);
System.arraycopy(x, 0, wtable, iw2, n+1);
wtable[0+iw1]=0.5*wtable[0+iw2];
for(i=2; i<n; i+=2)
{
wtable[i-1+iw1]=-wtable[i+iw2];
wtable[i+iw1]=wtable[i-2+iw1]+wtable[i-1+iw2];
}
if(modn==0)
wtable[n-1+iw1]=-wtable[n+iw2];
}
for(i=0; i<n; i++)
{
wtable[i+iw2]=war[i];
war[i]=wtable[i+iw1];
}
}
/*----------------
sint: sine FFT
---------------*/
void sint(int n, double x[], double wtable[])
{
sint1(n, x, wtable);
}
/*----------------------------------------------------------------------
sinti: initialization of sin-FFT
----------------------------------------------------------------------*/
void sinti(int n, double wtable[])
{
final double pi=Math.PI; //3.14159265358979;
int k, ns2;
double dt;
if(n<=1) return;
ns2=n / 2;
dt=pi /(double)(n+1);
for(k=0; k<ns2; k++)
wtable[k]=2*Math.sin((k+1)*dt);
rffti1(n+1, wtable, ns2);
}
}
package ca.uol.aig.fftpack;
/**
* FFT transform of a complex periodic sequence.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class ComplexDoubleFFT extends ComplexDoubleFFT_Mixed
{
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward transform
* (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
public double norm_factor;
private double wavetable[];
private int ndim;
/**
* Construct a wavenumber table with size <em>n</em> for Complex FFT.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a complex data sequence. When <em>n</em> is a multiplication of small
* numbers (4, 2, 3, 5), this FFT transform is very efficient.
*/
public ComplexDoubleFFT(int n)
{
ndim = n;
norm_factor = n;
if(wavetable == null || wavetable.length !=(4*ndim+15))
{
wavetable = new double[4*ndim + 15];
}
cffti(ndim, wavetable);
}
/**
* Forward complex FFT transform.
*
* @param x 2*<em>n</em> real double data representing <em>n</em> complex double data.
* As an input parameter, <em>x</em> is an array of 2*<em>n</em> real
* data representing <em>n</em> complex data. As an output parameter, <em>x</em> represents <em>n</em>
* FFT'd complex data. Their relation as follows:
* <br>
* x[2*i] is the real part of <em>i</em>-th complex data;
* <br>
* x[2*i+1] is the imaginary part of <em>i</em>-the complex data.
*
*/
public void ft(double x[])
{
if(x.length != 2*ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
cfftf(ndim, x, wavetable);
}
/**
* Forward complex FFT transform.
*
* @param x an array of <em>n</em> Complex data
*/
public void ft(Complex1D x)
{
if(x.x.length != ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
double[] y = new double[2*ndim];
for(int i=0; i<ndim; i++)
{
y[2*i] = x.x[i];
y[2*i+1] = x.y[i];
}
cfftf(ndim, y, wavetable);
for(int i=0; i<ndim; i++)
{
x.x[i]=y[2*i];
x.y[i]=y[2*i+1];
}
}
/**
* Backward complex FFT transform. It is the unnormalized inverse transform of <em>ft</em>(double[]).
*
* @param x 2*<em>n</em> real double data representing <em>n</em> complex double data.
*
* As an input parameter, <em>x</em> is an array of 2*<em>n</em>
* real data representing <em>n</em> complex data. As an output parameter, <em>x</em> represents
* <em>n</em> FFT'd complex data. Their relation as follows:
* <br>
* x[2*<em>i</em>] is the real part of <em>i</em>-th complex data;
* <br>
* x[2*<em>i</em>+1] is the imaginary part of <em>i</em>-the complex data.
*
*/
public void bt(double x[])
{
if(x.length != 2*ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
cfftb(ndim, x, wavetable);
}
/**
* Backward complex FFT transform. It is the unnormalized inverse transform of <em>ft</em>(Complex1D[]).
*
*
* @param x an array of <em>n</em> Complex data
*/
public void bt(Complex1D x)
{
if(x.x.length != ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
double[] y = new double[2*ndim];
for(int i=0; i<ndim; i++)
{
y[2*i] = x.x[i];
y[2*i+1] = x.y[i];
}
cfftb(ndim, y, wavetable);
for(int i=0; i<ndim; i++)
{
x.x[i]=y[2*i];
x.y[i]=y[2*i+1];
}
}
}
package ca.uol.aig.fftpack;
/**
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
class ComplexDoubleFFT_Mixed
{
/*----------------------------------------------------------------------
passf2: Complex FFT's forward/backward processing of factor 2;
isign is +1 for backward and -1 for forward transforms
----------------------------------------------------------------------*/
void passf2(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset, int isign)
/*isign==+1 for backward transform*/
{
int i, k, ah, ac;
double ti2, tr2;
int iw1;
iw1 = offset;
if(ido<=2)
{
for(k=0; k<l1; k++)
{
ah=k*ido;
ac=2*k*ido;
ch[ah]=cc[ac]+cc[ac+ido];
ch[ah+ido*l1]=cc[ac]-cc[ac+ido];
ch[ah+1]=cc[ac+1]+cc[ac+ido+1];
ch[ah+ido*l1+1]=cc[ac+1]-cc[ac+ido+1];
}
}
else
{
for(k=0; k<l1; k++)
{
for(i=0; i<ido-1; i+=2)
{
ah=i+k*ido;
ac=i+2*k*ido;
ch[ah]=cc[ac]+cc[ac+ido];
tr2=cc[ac]-cc[ac+ido];
ch[ah+1]=cc[ac+1]+cc[ac+1+ido];
ti2=cc[ac+1]-cc[ac+1+ido];
ch[ah+l1*ido+1]=wtable[i+iw1]*ti2+isign*wtable[i+1+iw1]*tr2;
ch[ah+l1*ido]=wtable[i+iw1]*tr2-isign*wtable[i+1+iw1]*ti2;
}
}
}
}
/*----------------------------------------------------------------------
passf3: Complex FFT's forward/backward processing of factor 3;
isign is +1 for backward and -1 for forward transforms
----------------------------------------------------------------------*/
void passf3(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset, int isign)
{
final double taur=-0.5;
final double taui=0.866025403784439;
int i, k, ac, ah;
double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
int iw1, iw2;
iw1 = offset;
iw2 = iw1 + ido;
if(ido==2)
{
for(k=1; k<=l1; k++)
{
ac=(3*k-2)*ido;
tr2=cc[ac]+cc[ac+ido];
cr2=cc[ac-ido]+taur*tr2;
ah=(k-1)*ido;
ch[ah]=cc[ac-ido]+tr2;
ti2=cc[ac+1]+cc[ac+ido+1];
ci2=cc[ac-ido+1]+taur*ti2;
ch[ah+1]=cc[ac-ido+1]+ti2;
cr3=isign*taui*(cc[ac]-cc[ac+ido]);
ci3=isign*taui*(cc[ac+1]-cc[ac+ido+1]);
ch[ah+l1*ido]=cr2-ci3;
ch[ah+2*l1*ido]=cr2+ci3;
ch[ah+l1*ido+1]=ci2+cr3;
ch[ah+2*l1*ido+1]=ci2-cr3;
}
}
else
{
for(k=1; k<=l1; k++)
{
for(i=0; i<ido-1; i+=2)
{
ac=i+(3*k-2)*ido;
tr2=cc[ac]+cc[ac+ido];
cr2=cc[ac-ido]+taur*tr2;
ah=i+(k-1)*ido;
ch[ah]=cc[ac-ido]+tr2;
ti2=cc[ac+1]+cc[ac+ido+1];
ci2=cc[ac-ido+1]+taur*ti2;
ch[ah+1]=cc[ac-ido+1]+ti2;
cr3=isign*taui*(cc[ac]-cc[ac+ido]);
ci3=isign*taui*(cc[ac+1]-cc[ac+ido+1]);
dr2=cr2-ci3;
dr3=cr2+ci3;
di2=ci2+cr3;
di3=ci2-cr3;
ch[ah+l1*ido+1]=wtable[i+iw1]*di2+isign*wtable[i+1+iw1]*dr2;
ch[ah+l1*ido]=wtable[i+iw1]*dr2-isign*wtable[i+1+iw1]*di2;
ch[ah+2*l1*ido+1]=wtable[i+iw2]*di3+isign*wtable[i+1+iw2]*dr3;
ch[ah+2*l1*ido]=wtable[i+iw2]*dr3-isign*wtable[i+1+iw2]*di3;
}
}
}
}
/*----------------------------------------------------------------------
passf4: Complex FFT's forward/backward processing of factor 4;
isign is +1 for backward and -1 for forward transforms
----------------------------------------------------------------------*/
void passf4(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset, int isign)
{
int i, k, ac, ah;
double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
int iw1, iw2, iw3;
iw1 = offset;
iw2 = iw1 + ido;
iw3 = iw2 + ido;
if(ido==2)
{
for(k=0; k<l1; k++)
{
ac=4*k*ido+1;
ti1=cc[ac]-cc[ac+2*ido];
ti2=cc[ac]+cc[ac+2*ido];
tr4=cc[ac+3*ido]-cc[ac+ido];
ti3=cc[ac+ido]+cc[ac+3*ido];
tr1=cc[ac-1]-cc[ac+2*ido-1];
tr2=cc[ac-1]+cc[ac+2*ido-1];
ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
ah=k*ido;
ch[ah]=tr2+tr3;
ch[ah+2*l1*ido]=tr2-tr3;
ch[ah+1]=ti2+ti3;
ch[ah+2*l1*ido+1]=ti2-ti3;
ch[ah+l1*ido]=tr1+isign*tr4;
ch[ah+3*l1*ido]=tr1-isign*tr4;
ch[ah+l1*ido+1]=ti1+isign*ti4;
ch[ah+3*l1*ido+1]=ti1-isign*ti4;
}
}
else
{
for(k=0; k<l1; k++)
{
for(i=0; i<ido-1; i+=2)
{
ac=i+1+4*k*ido;
ti1=cc[ac]-cc[ac+2*ido];
ti2=cc[ac]+cc[ac+2*ido];
ti3=cc[ac+ido]+cc[ac+3*ido];
tr4=cc[ac+3*ido]-cc[ac+ido];
tr1=cc[ac-1]-cc[ac+2*ido-1];
tr2=cc[ac-1]+cc[ac+2*ido-1];
ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
ah=i+k*ido;
ch[ah]=tr2+tr3;
cr3=tr2-tr3;
ch[ah+1]=ti2+ti3;
ci3=ti2-ti3;
cr2=tr1+isign*tr4;
cr4=tr1-isign*tr4;
ci2=ti1+isign*ti4;
ci4=ti1-isign*ti4;
ch[ah+l1*ido]=wtable[i+iw1]*cr2-isign*wtable[i+1+iw1]*ci2;
ch[ah+l1*ido+1]=wtable[i+iw1]*ci2+isign*wtable[i+1+iw1]*cr2;
ch[ah+2*l1*ido]=wtable[i+iw2]*cr3-isign*wtable[i+1+iw2]*ci3;
ch[ah+2*l1*ido+1]=wtable[i+iw2]*ci3+isign*wtable[i+1+iw2]*cr3;
ch[ah+3*l1*ido]=wtable[i+iw3]*cr4-isign*wtable[i+1+iw3]*ci4;
ch[ah+3*l1*ido+1]=wtable[i+iw3]*ci4+isign*wtable[i+1+iw3]*cr4;
}
}
}
}
/*----------------------------------------------------------------------
passf5: Complex FFT's forward/backward processing of factor 5;
isign is +1 for backward and -1 for forward transforms
----------------------------------------------------------------------*/
void passf5(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset, int isign)
/*isign==-1 for forward transform and+1 for backward transform*/
{
final double tr11=0.309016994374947;
final double ti11=0.951056516295154;
final double tr12=-0.809016994374947;
final double ti12=0.587785252292473;
int i, k, ac, ah;
double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
int iw1, iw2, iw3, iw4;
iw1 = offset;
iw2 = iw1 + ido;
iw3 = iw2 + ido;
iw4 = iw3 + ido;
if(ido==2)
{
for(k=1; k<=l1;++k)
{
ac=(5*k-4)*ido+1;
ti5=cc[ac]-cc[ac+3*ido];
ti2=cc[ac]+cc[ac+3*ido];
ti4=cc[ac+ido]-cc[ac+2*ido];
ti3=cc[ac+ido]+cc[ac+2*ido];
tr5=cc[ac-1]-cc[ac+3*ido-1];
tr2=cc[ac-1]+cc[ac+3*ido-1];
tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
ah=(k-1)*ido;
ch[ah]=cc[ac-ido-1]+tr2+tr3;
ch[ah+1]=cc[ac-ido]+ti2+ti3;
cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
cr5=isign*(ti11*tr5+ti12*tr4);
ci5=isign*(ti11*ti5+ti12*ti4);
cr4=isign*(ti12*tr5-ti11*tr4);
ci4=isign*(ti12*ti5-ti11*ti4);
ch[ah+l1*ido]=cr2-ci5;
ch[ah+4*l1*ido]=cr2+ci5;
ch[ah+l1*ido+1]=ci2+cr5;
ch[ah+2*l1*ido+1]=ci3+cr4;
ch[ah+2*l1*ido]=cr3-ci4;
ch[ah+3*l1*ido]=cr3+ci4;
ch[ah+3*l1*ido+1]=ci3-cr4;
ch[ah+4*l1*ido+1]=ci2-cr5;
}
}
else
{
for(k=1; k<=l1; k++)
{
for(i=0; i<ido-1; i+=2)
{
ac=i+1+(k*5-4)*ido;
ti5=cc[ac]-cc[ac+3*ido];
ti2=cc[ac]+cc[ac+3*ido];
ti4=cc[ac+ido]-cc[ac+2*ido];
ti3=cc[ac+ido]+cc[ac+2*ido];
tr5=cc[ac-1]-cc[ac+3*ido-1];
tr2=cc[ac-1]+cc[ac+3*ido-1];
tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
ah=i+(k-1)*ido;
ch[ah]=cc[ac-ido-1]+tr2+tr3;
ch[ah+1]=cc[ac-ido]+ti2+ti3;
cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
cr5=isign*(ti11*tr5+ti12*tr4);
ci5=isign*(ti11*ti5+ti12*ti4);
cr4=isign*(ti12*tr5-ti11*tr4);
ci4=isign*(ti12*ti5-ti11*ti4);
dr3=cr3-ci4;
dr4=cr3+ci4;
di3=ci3+cr4;
di4=ci3-cr4;
dr5=cr2+ci5;
dr2=cr2-ci5;
di5=ci2-cr5;
di2=ci2+cr5;
ch[ah+l1*ido]=wtable[i+iw1]*dr2-isign*wtable[i+1+iw1]*di2;
ch[ah+l1*ido+1]=wtable[i+iw1]*di2+isign*wtable[i+1+iw1]*dr2;
ch[ah+2*l1*ido]=wtable[i+iw2]*dr3-isign*wtable[i+1+iw2]*di3;
ch[ah+2*l1*ido+1]=wtable[i+iw2]*di3+isign*wtable[i+1+iw2]*dr3;
ch[ah+3*l1*ido]=wtable[i+iw3]*dr4-isign*wtable[i+1+iw3]*di4;
ch[ah+3*l1*ido+1]=wtable[i+iw3]*di4+isign*wtable[i+1+iw3]*dr4;
ch[ah+4*l1*ido]=wtable[i+iw4]*dr5-isign*wtable[i+1+iw4]*di5;
ch[ah+4*l1*ido+1]=wtable[i+iw4]*di5+isign*wtable[i+1+iw4]*dr5;
}
}
}
}
/*----------------------------------------------------------------------
passfg: Complex FFT's forward/backward processing of general factor;
isign is +1 for backward and -1 for forward transforms
----------------------------------------------------------------------*/
void passfg(int nac[], int ido, int ip, int l1, int idl1,
final double cc[], double c1[], double c2[], double ch[], double ch2[],
final double wtable[], int offset, int isign)
{
int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp;
double wai, war;
int iw1;
iw1 = offset;
idot=ido / 2;
ipph=(ip+1)/ 2;
idp=ip*ido;
if(ido>=l1)
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(k=0; k<l1; k++)
{
for(i=0; i<ido; i++)
{
ch[i+(k+j*l1)*ido]=cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
ch[i+(k+jc*l1)*ido]=cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
}
}
}
for(k=0; k<l1; k++)
for(i=0; i<ido; i++)
ch[i+k*ido]=cc[i+k*ip*ido];
}
else
{
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(i=0; i<ido; i++)
{
for(k=0; k<l1; k++)
{
ch[i+(k+j*l1)*ido]=cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
ch[i+(k+jc*l1)*ido]=cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
}
}
}
for(i=0; i<ido; i++)
for(k=0; k<l1; k++)
ch[i+k*ido]=cc[i+k*ip*ido];
}
idl=2-ido;
inc=0;
for(l=1; l<ipph; l++)
{
lc=ip-l;
idl+=ido;
for(ik=0; ik<idl1; ik++)
{
c2[ik+l*idl1]=ch2[ik]+wtable[idl-2+iw1]*ch2[ik+idl1];
c2[ik+lc*idl1]=isign*wtable[idl-1+iw1]*ch2[ik+(ip-1)*idl1];
}
idlj=idl;
inc+=ido;
for(j=2; j<ipph; j++)
{
jc=ip-j;
idlj+=inc;
if(idlj>idp) idlj-=idp;
war=wtable[idlj-2+iw1];
wai=wtable[idlj-1+iw1];
for(ik=0; ik<idl1; ik++)
{
c2[ik+l*idl1]+=war*ch2[ik+j*idl1];
c2[ik+lc*idl1]+=isign*wai*ch2[ik+jc*idl1];
}
}
}
for(j=1; j<ipph; j++)
for(ik=0; ik<idl1; ik++)
ch2[ik]+=ch2[ik+j*idl1];
for(j=1; j<ipph; j++)
{
jc=ip-j;
for(ik=1; ik<idl1; ik+=2)
{
ch2[ik-1+j*idl1]=c2[ik-1+j*idl1]-c2[ik+jc*idl1];
ch2[ik-1+jc*idl1]=c2[ik-1+j*idl1]+c2[ik+jc*idl1];
ch2[ik+j*idl1]=c2[ik+j*idl1]+c2[ik-1+jc*idl1];
ch2[ik+jc*idl1]=c2[ik+j*idl1]-c2[ik-1+jc*idl1];
}
}
nac[0]=1;
if(ido==2) return;
nac[0]=0;
for(ik=0; ik<idl1; ik++) c2[ik]=ch2[ik];
for(j=1; j<ip; j++)
{
for(k=0; k<l1; k++)
{
c1[(k+j*l1)*ido+0]=ch[(k+j*l1)*ido+0];
c1[(k+j*l1)*ido+1]=ch[(k+j*l1)*ido+1];
}
}
if(idot<=l1)
{
idij=0;
for(j=1; j<ip; j++)
{
idij+=2;
for(i=3; i<ido; i+=2)
{
idij+=2;
for(k=0; k<l1; k++)
{
c1[i-1+(k+j*l1)*ido]=
wtable[idij-2+iw1]*ch[i-1+(k+j*l1)*ido]-
isign*wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido];
c1[i+(k+j*l1)*ido]=
wtable[idij-2+iw1]*ch[i+(k+j*l1)*ido]+
isign*wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido];
}
}
}
}
else
{
idj=2-ido;
for(j=1; j<ip; j++)
{
idj+=ido;
for(k=0; k<l1; k++)
{
idij=idj;
for(i=3; i<ido; i+=2)
{
idij+=2;
c1[i-1+(k+j*l1)*ido]=
wtable[idij-2+iw1]*ch[i-1+(k+j*l1)*ido]-
isign*wtable[idij-1+iw1]*ch[i+(k+j*l1)*ido];
c1[i+(k+j*l1)*ido]=
wtable[idij-2+iw1]*ch[i+(k+j*l1)*ido]+
isign*wtable[idij-1+iw1]*ch[i-1+(k+j*l1)*ido];
}
}
}
}
}
/*---------------------------------------------------------
cfftf1: further processing of Complex forward FFT
--------------------------------------------------------*/
void cfftf1(int n, double c[], final double wtable[], int isign)
{
int idot, i;
int k1, l1, l2;
int na, nf, ip, iw, ido, idl1;
int[] nac = new int[1];
int iw1, iw2;
double[] ch = new double[2*n];
iw1=2*n;
iw2=4*n;
System.arraycopy(wtable, 0, ch, 0, 2*n);
nac[0] = 0;
nf=(int)wtable[1+iw2];
na=0;
l1=1;
iw=iw1;
for(k1=2; k1<=nf+1; k1++)
{
ip=(int)wtable[k1+iw2];
l2=ip*l1;
ido=n / l2;
idot=ido+ido;
idl1=idot*l1;
if(ip==4)
{
if(na==0)
{
passf4(idot, l1, c, ch, wtable, iw, isign);
}
else
{
passf4(idot, l1, ch, c, wtable, iw, isign);
}
na=1-na;
}
else if(ip==2)
{
if(na==0)
{
passf2(idot, l1, c, ch, wtable, iw, isign);
}
else
{
passf2(idot, l1, ch, c, wtable, iw, isign);
}
na=1-na;
}
else if(ip==3)
{
if(na==0)
{
passf3(idot, l1, c, ch, wtable, iw, isign);
}
else
{
passf3(idot, l1, ch, c, wtable, iw, isign);
}
na=1-na;
}
else if(ip==5)
{
if(na==0)
{
passf5(idot, l1, c, ch, wtable, iw, isign);
}
else
{
passf5(idot, l1, ch, c, wtable, iw, isign);
}
na=1-na;
}
else
{
if(na==0)
{
passfg(nac, idot, ip, l1, idl1, c, c, c, ch, ch, wtable, iw, isign);
}
else
{
passfg(nac, idot, ip, l1, idl1, ch, ch, ch, c, c, wtable, iw, isign);
}
if(nac[0] !=0) na=1-na;
}
l1=l2;
iw+=(ip-1)*idot;
}
if(na==0) return;
for(i=0; i<2*n; i++) c[i]=ch[i];
}
/*---------------------------------------------------------
cfftf: Complex forward FFT
--------------------------------------------------------*/
void cfftf(int n, double c[], double wtable[])
{
cfftf1(n, c, wtable, -1);
}
/*---------------------------------------------------------
cfftb: Complex borward FFT
--------------------------------------------------------*/
void cfftb(int n, double c[], double wtable[])
{
cfftf1(n, c, wtable, +1);
}
/*---------------------------------------------------------
cffti1: further initialization of Complex FFT
--------------------------------------------------------*/
void cffti1(int n, double wtable[])
{
final int[] ntryh = {3, 4, 2, 5};
final double twopi=2.0D*Math.PI;
double argh;
int idot, ntry=0, i, j;
double argld;
int i1, k1, l1, l2, ib;
double fi;
int ld, ii, nf, ip, nl, nq, nr;
double arg;
int ido, ipm;
nl=n;
nf=0;
j=0;
factorize_loop:
while(true)
{
j++;
if(j<=4)
ntry=ntryh[j-1];
else
ntry+=2;
do
{
nq=nl / ntry;
nr=nl-ntry*nq;
if(nr !=0) continue factorize_loop;
nf++;
wtable[nf+1+4*n]=ntry;
nl=nq;
if(ntry==2 && nf !=1)
{
for(i=2; i<=nf; i++)
{
ib=nf-i+2;
wtable[ib+1+4*n]=wtable[ib+4*n];
}
wtable[2+4*n]=2;
}
} while(nl !=1);
break factorize_loop;
}
wtable[0+4*n]=n;
wtable[1+4*n]=nf;
argh=twopi /(double)n;
i=1;
l1=1;
for(k1=1; k1<=nf; k1++)
{
ip=(int)wtable[k1+1+4*n];
ld=0;
l2=l1*ip;
ido=n / l2;
idot=ido+ido+2;
ipm=ip-1;
for(j=1; j<=ipm; j++)
{
i1=i;
wtable[i-1+2*n]=1;
wtable[i+2*n]=0;
ld+=l1;
fi=0;
argld=ld*argh;
for(ii=4; ii<=idot; ii+=2)
{
i+=2;
fi+=1;
arg=fi*argld;
wtable[i-1+2*n]=Math.cos(arg);
wtable[i+2*n]=Math.sin(arg);
}
if(ip>5)
{
wtable[i1-1+2*n]=wtable[i-1+2*n];
wtable[i1+2*n]=wtable[i+2*n];
}
}
l1=l2;
}
}
/*---------------------------------------------------------
cffti: Initialization of Real forward FFT
--------------------------------------------------------*/
void cffti(int n, double wtable[])
{
if(n==1) return;
cffti1(n, wtable);
} /*cffti*/
}
package ca.uol.aig.fftpack;
/**
* FFT transform of a real periodic sequence.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class RealDoubleFFT
extends RealDoubleFFT_Mixed
{
/**
* Construct a wavenumber table with size <em>n</em>.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a real data sequence. When <em>n</em> is a multiplication of small
* numbers (4, 2, 3, 5), this FFT transform is very efficient.
*/
public RealDoubleFFT(int n)
{
ndim = n;
norm_factor = n;
if(wavetable == null || wavetable.length !=(2*ndim+15))
{
wavetable = new double[2*ndim + 15];
}
rffti(ndim, wavetable);
}
/**
* Forward real FFT transform. It computes the discrete transform
* of a real data sequence.
*
* <p>The x parameter is both input and output data. After the FFT, x
* contains the transform coefficients used to construct n
* complex FFT coefficients.
*
* <p>The real part of the first complex FFT coefficients is x[0];
* its imaginary part is 0. If n is even set m = n/2, if n is odd set
* m = (n+1)/2, then for k = 1, ..., m-1:
* <ul>
* <li>the real part of k-th complex FFT coefficient is x[2 * k];
* <li>the imaginary part of k-th complex FFT coefficient is x[2 * k - 1].
* </ul>
*
* <p>If n is even, the real of part of (n/2)-th complex FFT
* coefficient is x[n - 1]; its imaginary part is 0.
*
* <p>The remaining complex FFT coefficients can be obtained by the
* symmetry relation: the (n-k)-th complex FFT coefficient is the
* conjugate of n-th complex FFT coefficient.
*
* @param x An array which contains the sequence to be
* transformed.
*/
public void ft(double[] x) {
if (x.length != ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
rfftf(ndim, x, wavetable);
}
/**
* Forward real FFT transform. It computes the discrete transform of a real data sequence.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients used to construct <em>n</em> complex FFT coeffients.
* <br>
* @param y the first complex (<em>n</em>+1)/2 (when <em>n</em> is odd) or (<em>n</em>/2+1) (when
* <em>n</em> is even) FFT coefficients.
* The remaining complex FFT coefficients can be obtained by the symmetry relation:
* the (<em>n</em>-<em>k</em>)-th complex FFT coefficient is the conjugate of <em>n</em>-th complex FFT coeffient.
*
*/
public void ft(double x[], Complex1D y) {
if (x.length != ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
rfftf(ndim, x, wavetable);
if(ndim%2 == 0)
{
y.x = new double[ndim/2 + 1];
y.y = new double[ndim/2 + 1];
}
else
{
y.x = new double[(ndim+1)/2];
y.y = new double[(ndim+1)/2];
}
y.x[0] = x[0];
y.y[0] = 0.0D;
for(int i=1; i<(ndim+1)/2; i++)
{
y.x[i] = x[2*i-1];
y.y[i] = x[2*i];
}
if(ndim%2 == 0)
{
y.x[ndim/2] = x[ndim-1];
y.y[ndim/2] = 0.0D;
}
}
/**
* Backward real FFT transform. It is the unnormalized inverse transform of <em>ft</em>(double[]).
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients. Also see the comments of <em>ft</em>(double[])
* for the relation between <em>x</em> and complex FFT coeffients.
*/
public void bt(double x[])
{
if(x.length != ndim)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
rfftb(ndim, x, wavetable);
}
/**
* Backward real FFT transform. It is the unnormalized inverse transform of <em>ft</em>(Complex1D, double[]).
*
* @param x an array which contains the sequence to be transformed. When <em>n</em> is odd, it contains the first
* (<em>n</em>+1)/2 complex data; when <em>n</em> is even, it contains (<em>n</em>/2+1) complex data.
* @param y the real FFT coeffients.
* <br>
* Also see the comments of <em>ft</em>(double[]) for the relation
* between <em>x</em> and complex FFT coeffients.
*/
public void bt(Complex1D x, double y[])
{
if(ndim%2 == 0)
{
if(x.x.length != ndim/2+1)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
}
else
{
if(x.x.length != (ndim+1)/2)
throw new IllegalArgumentException("The length of data can not match that of the wavetable");
}
y[0] = x.x[0];
for(int i=1; i<(ndim+1)/2; i++)
{
y[2*i-1]=x.x[i];
y[2*i]=x.y[i];
}
if(ndim%2 == 0)
{
y[ndim-1]=x.x[ndim/2];
}
rfftb(ndim, y, wavetable);
}
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward transform
* (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
public double norm_factor;
private double wavetable[];
private int ndim;
}
package ca.uol.aig.fftpack;
/**
* cosine FFT transform of a real even sequence.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class RealDoubleFFT_Even extends RealDoubleFFT_Mixed
{
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward transform
* (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
public double norm_factor;
private double wavetable[];
private int ndim;
/**
* Construct a wavenumber table with size <em>n</em>.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a real data sequence. When (<em>n</em>-1) is a multiplication of small
* numbers (4, 2, 3, 5), this FFT transform is very efficient.
*/
public RealDoubleFFT_Even(int n)
{
ndim = n;
norm_factor = 2 * (n - 1);
if (wavetable == null || wavetable.length != (3 * ndim + 15))
wavetable = new double[3 * ndim + 15];
costi(ndim, wavetable);
}
/**
* Forward cosine FFT transform. It computes the discrete sine transform of
* an odd sequence.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
public void ft(double[] x) {
cost(ndim, x, wavetable);
}
/**
* Backward cosine FFT transform. It is the unnormalized inverse transform of <em>ft</em>.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
public void bt(double[] x) {
cost(ndim, x, wavetable);
}
/*-------------------------------------------------------------
cost: cosine FFT. Backward and forward cos-FFT are the same.
------------------------------------------------------------*/
void cost(int n, double x[], final double wtable[])
{
int modn, i, k;
double c1, t1, t2;
int kc;
double xi;
int nm1;
double x1h;
int ns2;
double tx2, x1p3, xim2;
nm1=n-1;
ns2=n / 2;
if(n-2<0) return;
else if(n==2)
{
x1h=x[0]+x[1];
x[1]=x[0]-x[1];
x[0]=x1h;
}
else if(n==3)
{
x1p3=x[0]+x[2];
tx2=x[1]+x[1];
x[1]=x[0]-x[2];
x[0]=x1p3+tx2;
x[2]=x1p3-tx2;
}
else
{
c1=x[0]-x[n-1];
x[0]+=x[n-1];
for(k=1; k<ns2; k++)
{
kc=nm1-k;
t1=x[k]+x[kc];
t2=x[k]-x[kc];
c1+=wtable[kc]*t2;
t2=wtable[k]*t2;
x[k]=t1-t2;
x[kc]=t1+t2;
}
modn=n%2;
if(modn !=0) x[ns2]+=x[ns2];
rfftf1(nm1, x, wtable, n);
xim2=x[1];
x[1]=c1;
for(i=3; i<n; i+=2)
{
xi=x[i];
x[i]=x[i-2]-x[i-1];
x[i-1]=xim2;
xim2=xi;
}
if(modn !=0) x[n-1]=xim2;
}
}
/*----------------------------------
costi: initialization of cos-FFT
---------------------------------*/
void costi(int n, double wtable[])
{
final double pi=Math.PI; //3.14159265358979;
int k, kc, ns2;
double dt;
if(n<=3) return;
ns2=n / 2;
dt=pi /(double)(n-1);
for(k=1; k<ns2; k++)
{
kc=n-k-1;
wtable[k]=2*Math.sin(k*dt);
wtable[kc]=2*Math.cos(k*dt);
}
rffti1(n-1, wtable, n);
}
}
package ca.uol.aig.fftpack;
/**
* cosine FFT transform with odd wave numbers.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class RealDoubleFFT_Even_Odd extends RealDoubleFFT_Mixed
{
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward transform
* (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
public double norm_factor;
double wavetable[];
int ndim;
/**
* Construct a wavenumber table with size <em>n</em>.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a real data sequence. When <em>n</em> is a multiplication of small
* numbers(4, 2, 3, 5), this FFT transform is very efficient.
*/
public RealDoubleFFT_Even_Odd(int n) {
ndim = n;
norm_factor = 4*n;
if(wavetable == null || wavetable.length !=(3*ndim+15))
{
wavetable = new double[3*ndim + 15];
}
cosqi(ndim, wavetable);
}
/**
* Forward FFT transform of quarter wave data. It computes the coeffients in
* cosine series representation with only odd wave numbers.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
public void ft(double x[]) {
cosqf(ndim, x, wavetable);
}
/**
* Backward FFT transform of quarter wave data. It is the unnormalized inverse transform
* of <em>ft</em>.
*
* @param x an array which contains the sequence to be tranformed. After FFT, <em>x</em> contains
* the transform coeffients.
*/
public void bt(double x[]) {
cosqb(ndim, x, wavetable);
}
/*----------------------------------------------------------------------
cosqf1: further processing of forward cos-FFT with odd wave numbers.
----------------------------------------------------------------------*/
void cosqf1(int n, double x[], double wtable[])
{
int modn, i, k;
int kc, ns2;
double xim1;
ns2=(n+1)/ 2;
for(k=1; k<ns2; k++)
{
kc=n-k;
wtable[k+n]=x[k]+x[kc];
wtable[kc+n]=x[k]-x[kc];
}
modn=n%2;
if(modn==0) wtable[ns2+n]=x[ns2]+x[ns2];
for(k=1; k<ns2; k++)
{
kc=n-k;
x[k]=wtable[k-1]*wtable[kc+n]+wtable[kc-1]*wtable[k+n];
x[kc]=wtable[k-1]*wtable[k+n]-wtable[kc-1]*wtable[kc+n];
}
if(modn==0) x[ns2]=wtable[ns2-1]*wtable[ns2+n];
rfftf1(n, x, wtable, n);
for(i=2; i<n; i+=2)
{
xim1=x[i-1]-x[i];
x[i]=x[i-1]+x[i];
x[i-1]=xim1;
}
}
/*----------------------------------------------------------------------
cosqb1: further processing of backward cos-FFT with odd wave numbers.
----------------------------------------------------------------------*/
void cosqb1(int n, double x[], double wtable[])
{
int modn, i, k;
int kc, ns2;
double xim1;
ns2=(n+1)/ 2;
for(i=2; i<n; i+=2)
{
xim1=x[i-1]+x[i];
x[i]-=x[i-1];
x[i-1]=xim1;
}
x[0]+=x[0];
modn=n%2;
if(modn==0) x[n-1]+=x[n-1];
rfftb1(n, x, wtable, n);
for(k=1; k<ns2; k++)
{
kc=n-k;
wtable[k+n]=wtable[k-1]*x[kc]+wtable[kc-1]*x[k];
wtable[kc+n]=wtable[k-1]*x[k]-wtable[kc-1]*x[kc];
}
if(modn==0) x[ns2]=wtable[ns2-1]*(x[ns2]+x[ns2]);
for(k=1; k<ns2; k++)
{
kc=n-k;
x[k]=wtable[k+n]+wtable[kc+n];
x[kc]=wtable[k+n]-wtable[kc+n];
}
x[0]+=x[0];
}
/*-----------------------------------------------
cosqf: forward cosine FFT with odd wave numbers.
----------------------------------------------*/
void cosqf(int n, double x[], final double wtable[])
{
final double sqrt2=1.4142135623731;
double tsqx;
if(n<2)
{
return;
}
else if(n==2)
{
tsqx=sqrt2*x[1];
x[1]=x[0]-tsqx;
x[0]+=tsqx;
}
else
{
cosqf1(n, x, wtable);
}
}
/*-----------------------------------------------
cosqb: backward cosine FFT with odd wave numbers.
----------------------------------------------*/
void cosqb(int n, double x[], double wtable[])
{
final double tsqrt2=2.82842712474619;
double x1;
if(n<2)
{
x[0]*=4;
}
else if(n==2)
{
x1=4*(x[0]+x[1]);
x[1]=tsqrt2*(x[0]-x[1]);
x[0]=x1;
}
else
{
cosqb1(n, x, wtable);
}
}
/*-----------------------------------------------------------
cosqi: initialization of cosine FFT with odd wave numbers.
----------------------------------------------------------*/
void cosqi(int n, double wtable[])
{
final double pih=Math.PI/2.0D; //1.57079632679491;
int k;
double dt;
dt=pih / (double)n;
for(k=0; k<n; k++) wtable[k]=Math.cos((k+1)*dt);
rffti1(n, wtable, n);
}
}
package ca.uol.aig.fftpack;
/**
* sine FFT transform with odd wave numbers.
* @author Baoshe Zhang
* @author Astronomical Instrument Group of University of Lethbridge.
*/
public class RealDoubleFFT_Odd_Odd extends RealDoubleFFT_Even_Odd
{
/**
* <em>norm_factor</em> can be used to normalize this FFT transform. This is because
* a call of forward transform (<em>ft</em>) followed by a call of backward transform
* (<em>bt</em>) will multiply the input sequence by <em>norm_factor</em>.
*/
/**
* Construct a wavenumber table with size n.
* The sequences with the same size can share a wavenumber table. The prime
* factorization of <em>n</em> together with a tabulation of the trigonometric functions
* are computed and stored.
*
* @param n the size of a real data sequence. When <em>n</em> is a multiplication of small
* numbers (4, 2, 3, 5), this FFT transform is very efficient.
*/
public RealDoubleFFT_Odd_Odd(int n)
{
super(n);
}
/**
* Forward FFT transform of quarter wave data. It computes the coeffients in
* sine series representation with only odd wave numbers.
*
* @param x an array which contains the sequence to be transformed. After FFT,
* <em>x</em> contains the transform coeffients.
*/
@Override
public void ft(double x[])
{
sinqf(ndim, x, wavetable);
}
/**
* Backward FFT transform of quarter wave data. It is the unnormalized inverse transform
* of <em>ft</em>.
*
* @param x an array which contains the sequence to be tranformed. After FFT, <em>x</em> contains
* the transform coeffients.
*/
@Override
public void bt(double x[])
{
sinqb(ndim, x, wavetable);
}
/*-----------------------------------------------
sinqf: forward sine FFT with odd wave numbers.
----------------------------------------------*/
void sinqf(int n, double x[], double wtable[])
{
int k;
double xhold;
int kc, ns2;
if(n==1) return;
ns2=n / 2;
for(k=0; k<ns2; k++)
{
kc=n-k-1;
xhold=x[k];
x[k]=x[kc];
x[kc]=xhold;
}
cosqf(n, x, wtable);
for(k=1; k<n; k+=2) x[k]=-x[k];
}
/*-----------------------------------------------
sinqb: backward sine FFT with odd wave numbers.
----------------------------------------------*/
void sinqb(int n, double x[], double wtable[])
{
int k;
double xhold;
int kc, ns2;
if(n<=1)
{
x[0]*=4;
return;
}
ns2=n / 2;
for(k=1; k<n; k+=2) x[k]=-x[k];
cosqb(n, x, wtable);
for(k=0; k<ns2; k++)
{
kc=n-k-1;
xhold=x[k];
x[k]=x[kc];
x[kc]=xhold;
}
}
/*
void sinqi(int n, double wtable[])
{
cosqi(n, wtable);
}
*/
}
Related examples in the same category