Mega Code Archive

Audio processing library

//GNU General Public License, version 2 package ca.uol.aig.fftpack; /** * Construct a 1-D complex data sequence. */ public class Complex1D { /** * x[i] is the real part of i-th complex data. */ public double x[]; /** * y[i] is the imaginary part of i-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 { /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward * transform (bt) will multiply the input sequence by norm_factor. */ public double norm_factor; private double wavetable[]; private int ndim; /** * Construct a wavenumber table with size n. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a real data sequence. When (n+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, * x contains the transform coeffients. */ public void ft(double x[]) { sint(ndim, x, wavetable); } /** * Backward sine FFT transform. It is the unnormalized inverse transform of ft. * * @param x an array which contains the sequence to be transformed. After FFT, * x 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 { /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward transform * (bt) will multiply the input sequence by norm_factor. */ public double norm_factor; private double wavetable[]; private int ndim; /** * Construct a wavenumber table with size n for Complex FFT. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a complex data sequence. When n 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*n real double data representing n complex double data. * As an input parameter, x is an array of 2*n real * data representing n complex data. As an output parameter, x represents n * FFT'd complex data. Their relation as follows: * * x[2*i] is the real part of i-th complex data; * * x[2*i+1] is the imaginary part of i-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 n 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 ft(double[]). * * @param x 2*n real double data representing n complex double data. * * As an input parameter, x is an array of 2*n * real data representing n complex data. As an output parameter, x represents * n FFT'd complex data. Their relation as follows: * * x[2*i] is the real part of i-th complex data; * * x[2*i+1] is the imaginary part of i-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 ft(Complex1D[]). * * * @param x an array of n 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 n. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a real data sequence. When n 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. * * The x parameter is both input and output data. After the FFT, x * contains the transform coefficients used to construct n * complex FFT coefficients. * * 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> * * If n is even, the real of part of (n/2)-th complex FFT * coefficient is x[n - 1]; its imaginary part is 0. * * 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, * x contains the transform coeffients used to construct n complex FFT coeffients. * * @param y the first complex (n+1)/2 (when n is odd) or (n/2+1) (when * n is even) FFT coefficients. * 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 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 ft(double[]). * * @param x an array which contains the sequence to be transformed. After FFT, * x contains the transform coeffients. Also see the comments of ft(double[]) * for the relation between x 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 ft(Complex1D, double[]). * * @param x an array which contains the sequence to be transformed. When n is odd, it contains the first * (n+1)/2 complex data; when n is even, it contains (n/2+1) complex data. * @param y the real FFT coeffients. * * Also see the comments of ft(double[]) for the relation * between x 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); } /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward transform * (bt) will multiply the input sequence by norm_factor. */ 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 { /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward transform * (bt) will multiply the input sequence by norm_factor. */ public double norm_factor; private double wavetable[]; private int ndim; /** * Construct a wavenumber table with size n. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a real data sequence. When (n-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, * x contains the transform coeffients. */ public void ft(double[] x) { cost(ndim, x, wavetable); } /** * Backward cosine FFT transform. It is the unnormalized inverse transform of ft. * * @param x an array which contains the sequence to be transformed. After FFT, * x 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 { /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward transform * (bt) will multiply the input sequence by norm_factor. */ public double norm_factor; double wavetable[]; int ndim; /** * Construct a wavenumber table with size n. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a real data sequence. When n 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, * x 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 ft. * * @param x an array which contains the sequence to be tranformed. After FFT, x 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 { /** * norm_factor can be used to normalize this FFT transform. This is because * a call of forward transform (ft) followed by a call of backward transform * (bt) will multiply the input sequence by norm_factor. */ /** * Construct a wavenumber table with size n. * The sequences with the same size can share a wavenumber table. The prime * factorization of n together with a tabulation of the trigonometric functions * are computed and stored. * * @param n the size of a real data sequence. When n 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, * x 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 ft. * * @param x an array which contains the sequence to be tranformed. After FFT, x 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); } */ }