001package cnslab.cnsmath;
002public class FFT
003{
004        public static void SWAP(double [] data, int x, int y)
005        {
006                double tmp;
007                tmp = data[y];
008                data[y] = data[x];
009                data[x] = tmp;
010        }
011
012        public static void fourn(double [] data, int [] nn, int isign)
013        {
014                int idim,i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
015                int ibit,k1,k2,n,nprev,nrem,ntot;
016                double tempi,tempr,theta,wi,wpi,wpr,wr,wtemp;
017
018                int ndim=nn.length;
019                ntot=data.length/2;
020                nprev=1;
021                for (idim=ndim-1;idim>=0;idim--) {
022                        n=nn[idim];
023                        nrem=ntot/(n*nprev);
024                        ip1=nprev << 1;
025                        ip2=ip1*n;
026                        ip3=ip2*nrem;
027                        i2rev=0;
028                        for (i2=0;i2<ip2;i2+=ip1) {
029                                if (i2 < i2rev) {
030                                        for (i1=i2;i1<i2+ip1-1;i1+=2) {
031                                                for (i3=i1;i3<ip3;i3+=ip2) {
032                                                        i3rev=i2rev+i3-i2;
033                                                        SWAP(data,i3,i3rev);
034                                                        SWAP(data,i3+1,i3rev+1);
035                                                }
036                                        }
037                                }
038                                ibit=ip2 >> 1;
039                                while (ibit >= ip1 && i2rev+1 > ibit) {
040                                        i2rev -= ibit;
041                                        ibit >>= 1;
042                                }
043                                i2rev += ibit;
044                        }
045                        ifp1=ip1;
046                        while (ifp1 < ip2) {
047                                ifp2=ifp1 << 1;
048                                theta=isign*6.28318530717959/(ifp2/ip1);
049                                wtemp=Math.sin(0.5*theta);
050                                wpr= -2.0*wtemp*wtemp;
051                                wpi=Math.sin(theta);
052                                wr=1.0;
053                                wi=0.0;
054                                for (i3=0;i3<ifp1;i3+=ip1) {
055                                        for (i1=i3;i1<i3+ip1-1;i1+=2) {
056                                                for (i2=i1;i2<ip3;i2+=ifp2) {
057                                                        k1=i2;
058                                                        k2=k1+ifp1;
059                                                        tempr=wr*data[k2]-wi*data[k2+1];
060                                                        tempi=wr*data[k2+1]+wi*data[k2];
061                                                        data[k2]=data[k1]-tempr;
062                                                        data[k2+1]=data[k1+1]-tempi;
063                                                        data[k1] += tempr;
064                                                        data[k1+1] += tempi;
065                                                }
066                                        }
067                                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
068                                        wi=wi*wpr+wtemp*wpi+wi;
069                                }
070                                ifp1=ifp2;
071                        }
072                        nprev *= n;
073                }
074        }
075}