Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newfft.cpp

Go to the documentation of this file.
00001 //$ newfft.cpp
00002 
00003 // This is originally by Sande and Gentleman in 1967! I have translated from
00004 // Fortran into C and a little bit of C++.
00005 
00006 // It takes about twice as long as fftw
00007 // (http://theory.lcs.mit.edu/~fftw/homepage.html)
00008 // but is much shorter than fftw  and so despite its age
00009 // might represent a reasonable
00010 // compromise between speed and complexity.
00011 // If you really need the speed get fftw.
00012 
00013 
00014 //    THIS SUBROUTINE WAS WRITTEN BY G.SANDE OF PRINCETON UNIVERSITY AND
00015 //    W.M.GENTLMAN OF THE BELL TELEPHONE LAB.  IT WAS BROUGHT TO LONDON
00016 //    BY DR. M.D. GODFREY AT THE IMPERIAL COLLEGE AND WAS ADAPTED FOR
00017 //    BURROUGHS 6700 BY D. R. BRILLINGER AND J. PEMBERTON
00018 //    IT REPRESENTS THE STATE OF THE ART OF COMPUTING COMPLETE FINITE
00019 //    DISCRETE FOURIER TRANSFORMS AS OF NOV.1967.
00020 //    OTHER PROGRAMS REQUIRED.
00021 //                                 ONLY THOSE SUBROUTINES INCLUDED HERE.
00022 //                      USAGE.
00023 //       CALL AR1DFT(N,X,Y)
00024 //            WHERE  N IS THE NUMBER OF POINTS IN THE SEQUENCE .
00025 //                   X - IS A ONE-DIMENSIONAL ARRAY CONTAINING THE REAL
00026 //                       PART OF THE SEQUENCE.
00027 //                   Y - IS A ONE-DIMENSIONAL ARRAY CONTAINING THE
00028 //                       IMAGINARY PART OF THE SEQUENCE.
00029 //    THE TRANSFORM IS RETURNED IN X AND Y.
00030 //            METHOD
00031 //               FOR A GENERAL DISCUSSION OF THESE TRANSFORMS AND OF
00032 //    THE FAST METHOD FOR COMPUTING THEM, SEE GENTLEMAN AND SANDE,
00033 //    @FAST FOURIER TRANSFORMS - FOR FUN AND PROFIT,@ 1966 FALL JOINT
00034 //    COMPUTER CONFERENCE.
00035 //    THIS PROGRAM COMPUTES THIS FOR A COMPLEX SEQUENCE Z(T) OF LENGTH
00036 //    N WHOSE ELEMENTS ARE STORED AT(X(I) , Y(I)) AND RETURNS THE
00037 //    TRANSFORM COEFFICIENTS AT (X(I), Y(I)).
00038 //        DESCRIPTION
00039 //    AR1DFT IS A HIGHLY MODULAR ROUTINE CAPABLE OF COMPUTING IN PLACE
00040 //    THE COMPLETE FINITE DISCRETE FOURIER TRANSFORM  OF A ONE-
00041 //    DIMENSIONAL SEQUENCE OF RATHER GENERAL LENGTH N.
00042 //       THE MAIN ROUTINE , AR1DFT ITSELF, FACTORS N. IT THEN CALLS ON
00043 //    ON GR 1D FT TO COMPUTE THE ACTUAL TRANSFORMS, USING THESE FACTORS.
00044 //    THIS GR 1D FT DOES, CALLING AT EACH STAGE ON THE APPROPRIATE KERN
00045 //    EL R2FTK, R4FTK, R8FTK, R16FTK, R3FTK, R5FTK, OR RPFTK TO PERFORM
00046 //    THE COMPUTATIONS FOR THIS PASS OVER THE SEQUENCE, DEPENDING ON
00047 //    WHETHER THE CORRESPONDING FACTOR IS 2, 4, 8, 16, 3, 5, OR SOME
00048 //    MORE GENERAL PRIME P. WHEN GR1DFT IS FINISHED THE TRANSFORM IS
00049 //    COMPUTED, HOWEVER, THE RESULTS ARE STORED IN "DIGITS REVERSED"
00050 //    ORDER. AR1DFT THEREFORE, CALLS UPON GR 1S FS TO SORT THEM OUT.
00051 //    TO RETURN TO THE FACTORIZATION, SINGLETON HAS POINTED OUT THAT
00052 //    THE TRANSFORMS ARE MORE EFFICIENT IF THE SAMPLE SIZE N, IS OF THE
00053 //    FORM B*A**2 AND B CONSISTS OF A SINGLE FACTOR.  IN SUCH A CASE
00054 //    IF WE PROCESS THE FACTORS IN THE ORDER ABA  THEN
00055 //    THE REORDERING CAN BE DONE AS FAST IN PLACE, AS WITH SCRATCH
00056 //    STORAGE.  BUT AS B BECOMES MORE COMPLICATED, THE COST OF THE DIGIT
00057 //    REVERSING DUE TO B PART BECOMES VERY EXPENSIVE IF WE TRY TO DO IT
00058 //    IN PLACE.  IN SUCH A CASE IT MIGHT BE BETTER TO USE EXTRA STORAGE
00059 //    A ROUTINE TO DO THIS IS, HOWEVER, NOT INCLUDED HERE.
00060 //    ANOTHER FEATURE INFLUENCING THE FACTORIZATION IS THAT FOR ANY FIXED
00061 //    FACTOR N WE CAN PREPARE A SPECIAL KERNEL WHICH WILL COMPUTE
00062 //    THAT STAGE OF THE TRANSFORM MORE EFFICIENTLY THAN WOULD A KERNEL
00063 //    FOR GENERAL FACTORS, ESPECIALLY IF THE GENERAL KERNEL HAD TO BE
00064 //    APPLIED SEVERAL TIMES. FOR EXAMPLE, FACTORS OF 4 ARE MORE
00065 //    EFFICIENT THAN FACTORS OF 2, FACTORS OF 8 MORE EFFICIENT THAN 4,ETC
00066 //    ON THE OTHER HAND DIMINISHING RETURNS RAPIDLY SET IN, ESPECIALLY
00067 //    SINCE THE LENGTH OF THE KERNEL FOR A SPECIAL CASE IS ROUGHLY
00068 //    PROPORTIONAL TO THE FACTOR IT DEALS WITH. HENCE THESE PROBABLY ARE
00069 //    ALL THE KERNELS WE WISH TO HAVE.
00070 //            RESTRICTIONS.
00071 //    AN UNFORTUNATE FEATURE OF THE SORTING PROBLEM IS THAT THE MOST
00072 //    EFFICIENT WAY TO DO IT IS WITH NESTED DO LOOPS, ONE FOR EACH
00073 //    FACTOR. THIS PUTS A RESTRICTION ON N AS TO HOW MANY FACTORS IT
00074 //    CAN HAVE.  CURRENTLY THE LIMIT IS 16, BUT THE LIMIT CAN BE READILY
00075 //    RAISED IF NECESSARY.
00076 //    A SECOND RESTRICTION OF THE PROGRAM IS THAT LOCAL STORAGE OF THE
00077 //    THE ORDER P**2 IS REQUIRED BY THE GENERAL KERNEL RPFTK, SO SOME
00078 //    LIMIT MUST BE SET ON P.  CURRENTLY THIS IS 19, BUT IT CAN BE INCRE
00079 //    INCREASED BY TRIVIAL CHANGES.
00080 //       OTHER COMMENTS.
00081 //(1) THE ROUTINE IS ADAPTED TO CHECK WHETHER A GIVEN N WILL MEET THE
00082 //    ABOVE FACTORING REQUIREMENTS AN, IF NOT, TO RETURN THE NEXT HIGHER
00083 //    NUMBER, NX, SAY, WHICH WILL MEET THESE REQUIREMENTS.
00084 //    THIS CAN BE ACCHIEVED BY   A STATEMENT OF THE FORM
00085 //            CALL FACTR(N,X,Y).
00086 //    IF A DIFFERENT N, SAY NX, IS RETURNED THEN THE TRANSFORMS COULD BE
00087 //    OBTAINED BY EXTENDING THE SIZE OF THE X-ARRAY AND Y-ARRAY TO NX,
00088 //    AND SETTING X(I) = Y(I) = 0., FOR I = N+1, NX.
00089 //(2) IF THE SEQUENCE Z IS ONLY A REAL SEQUENCE, THEN THE IMAGINARY PART
00090 //    Y(I)=0., THIS WILL RETURN THE COSINE TRANSFORM OF THE REAL SEQUENCE
00091 //    IN X, AND THE SINE TRANSFORM IN Y.
00092 
00093 
00094 #ifndef WANT_STREAM
00095 #define WANT_STREAM
00096 #endif
00097 
00098 #ifndef WANT_MATH
00099 #define WANT_MATH
00100 #endif
00101 
00102 #include "newmatap.h"
00103 
00104 #ifdef use_namespace
00105 namespace NEWMAT {
00106 #endif
00107 
00108 #ifdef DO_REPORT
00109 #define REPORT { static ExeCounter ExeCount(__LINE__,20); ++ExeCount; }
00110 #else
00111 #define REPORT {}
00112 #endif
00113 
00114 inline Real square(Real x) { return x*x; }
00115 inline int square(int x) { return x*x; }
00116 
00117 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00118    const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00119    Real* X, Real* Y);
00120 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00121    Real* X, Real* Y);
00122 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y);
00123 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1);
00124 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00125    Real* X2, Real* Y2);
00126 static void R_4_FTK (int N, int M,
00127    Real* X0, Real* Y0, Real* X1, Real* Y1,
00128    Real* X2, Real* Y2, Real* X3, Real* Y3);
00129 static void R_5_FTK (int N, int M,
00130    Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00131    Real* X3, Real* Y3, Real* X4, Real* Y4);
00132 static void R_8_FTK (int N, int M,
00133    Real* X0, Real* Y0, Real* X1, Real* Y1,
00134    Real* X2, Real* Y2, Real* X3, Real* Y3,
00135    Real* X4, Real* Y4, Real* X5, Real* Y5,
00136    Real* X6, Real* Y6, Real* X7, Real* Y7);
00137 static void R_16_FTK (int N, int M,
00138    Real* X0, Real* Y0, Real* X1, Real* Y1,
00139    Real* X2, Real* Y2, Real* X3, Real* Y3,
00140    Real* X4, Real* Y4, Real* X5, Real* Y5,
00141    Real* X6, Real* Y6, Real* X7, Real* Y7,
00142    Real* X8, Real* Y8, Real* X9, Real* Y9,
00143    Real* X10, Real* Y10, Real* X11, Real* Y11,
00144    Real* X12, Real* Y12, Real* X13, Real* Y13,
00145    Real* X14, Real* Y14, Real* X15, Real* Y15);
00146 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f);
00147 
00148 
00149 bool FFT_Controller::ar_1d_ft (int PTS, Real* X, Real *Y)
00150 {
00151 //    ARBITRARY RADIX ONE DIMENSIONAL FOURIER TRANSFORM
00152 
00153    REPORT
00154 
00155    int  F,J,N,NF,P,PMAX,P_SYM,P_TWO,Q,R,TWO_GRP;
00156 
00157    // NP is maximum number of squared factors allows PTS up to 2**32 at least
00158    // NQ is number of not-squared factors - increase if we increase PMAX
00159    const int NP = 16, NQ = 10;
00160    SimpleIntArray PP(NP), QQ(NQ);
00161 
00162    TWO_GRP=16; PMAX=19;
00163 
00164    // PMAX is the maximum factor size
00165    // TWO_GRP is the maximum power of 2 handled as a single factor
00166    // Doesn't take advantage of combining powers of 2 when calculating
00167    // number of factors
00168 
00169    if (PTS<=1) return true;
00170    N=PTS; P_SYM=1; F=2; P=0; Q=0;
00171 
00172    // P counts the number of squared factors
00173    // Q counts the number of the rest
00174    // R = 0 for no non-squared factors; 1 otherwise
00175 
00176    // FACTOR holds all the factors - non-squared ones in the middle
00177    //   - length is 2*P+Q
00178    // SYM also holds all the factors but with the non-squared ones
00179    //   multiplied together - length is 2*P+R
00180    // PP holds the values of the squared factors - length is P
00181    // QQ holds the values of the rest - length is Q
00182 
00183    // P_SYM holds the product of the squared factors
00184 
00185    // find the factors - load into PP and QQ
00186    while (N > 1)
00187    {
00188       bool fail = true;
00189       for (J=F; J<=PMAX; J++)
00190          if (N % J == 0) { fail = false; F=J; break; }
00191       if (fail || P >= NP || Q >= NQ) return false; // can't factor
00192       N /= F;
00193       if (N % F != 0) QQ[Q++] = F;
00194       else { N /= F; PP[P++] = F; P_SYM *= F; }
00195    }
00196 
00197    R = (Q == 0) ? 0 : 1;  // R = 0 if no not-squared factors, 1 otherwise
00198 
00199    NF = 2*P + Q;
00200    SimpleIntArray FACTOR(NF + 1), SYM(2*P + R);
00201    FACTOR[NF] = 0;                // we need this in the "combine powers of 2"
00202 
00203    // load into SYM and FACTOR
00204    for (J=0; J<P; J++)
00205       { SYM[J]=FACTOR[J]=PP[P-1-J]; FACTOR[P+Q+J]=SYM[P+R+J]=PP[J]; }
00206 
00207    if (Q>0)
00208    {
00209       REPORT
00210       for (J=0; J<Q; J++) FACTOR[P+J]=QQ[J];
00211       SYM[P]=PTS/square(P_SYM);
00212    }
00213 
00214    // combine powers of 2
00215    P_TWO = 1;
00216    for (J=0; J < NF; J++)
00217    {
00218       if (FACTOR[J]!=2) continue;
00219       P_TWO=P_TWO*2; FACTOR[J]=1;
00220       if (P_TWO<TWO_GRP && FACTOR[J+1]==2) continue;
00221       FACTOR[J]=P_TWO; P_TWO=1;
00222    }
00223 
00224    if (P==0) R=0;
00225    if (Q<=1) Q=0;
00226 
00227    // do the analysis
00228    GR_1D_FT(PTS,NF,FACTOR,X,Y);                 // the transform
00229    GR_1D_FS(PTS,2*P+R,Q,SYM,P_SYM,QQ,X,Y);      // the reshuffling
00230 
00231    return true;
00232 
00233 }
00234 
00235 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00236    const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00237    Real* X, Real* Y)
00238 {
00239 //    GENERAL RADIX ONE DIMENSIONAL FOURIER SORT
00240 
00241 // PTS = number of points
00242 // N_SYM = length of SYM
00243 // N_UN_SYM = length of UN_SYM
00244 // SYM: squared factors + product of non-squared factors + squared factors
00245 // P_SYM = product of squared factors (each included only once)
00246 // UN_SYM: not-squared factors
00247 
00248    REPORT
00249 
00250    Real T;
00251    int  JJ,KK,P_UN_SYM;
00252 
00253    // I have replaced the multiple for-loop used by Sande-Gentleman code
00254    // by the following code which does not limit the number of factors
00255 
00256    if (N_SYM > 0)
00257    {
00258       REPORT
00259       SimpleIntArray U(N_SYM);
00260       for(MultiRadixCounter MRC(N_SYM, SYM, U); !MRC.Finish(); ++MRC)
00261       {
00262          if (MRC.Swap())
00263          {
00264             int P = MRC.Reverse(); int JJ = MRC.Counter(); Real T;
00265             T=X[JJ]; X[JJ]=X[P]; X[P]=T; T=Y[JJ]; Y[JJ]=Y[P]; Y[P]=T;
00266          }
00267       }
00268    }
00269 
00270    int J,JL,K,L,M,MS;
00271 
00272    // UN_SYM contains the non-squared factors
00273    // I have replaced the Sande-Gentleman code as it runs into
00274    // integer overflow problems
00275    // My code (and theirs) would be improved by using a bit array
00276    // as suggested by Van Loan
00277 
00278    if (N_UN_SYM==0) { REPORT return; }
00279    P_UN_SYM=PTS/square(P_SYM); JL=(P_UN_SYM-3)*P_SYM; MS=P_UN_SYM*P_SYM;
00280 
00281    for (J = P_SYM; J<=JL; J+=P_SYM)
00282    {
00283       K=J;
00284       do K = P_SYM * BitReverse(K / P_SYM, P_UN_SYM, N_UN_SYM, UN_SYM);
00285       while (K<J);
00286 
00287       if (K!=J)
00288       {
00289          REPORT
00290          for (L=0; L<P_SYM; L++) for (M=L; M<PTS; M+=MS)
00291          {
00292             JJ=M+J; KK=M+K;
00293             T=X[JJ]; X[JJ]=X[KK]; X[KK]=T; T=Y[JJ]; Y[JJ]=Y[KK]; Y[KK]=T;
00294          }
00295       }
00296    }
00297 
00298    return;
00299 }
00300 
00301 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00302    Real* X, Real* Y)
00303 {
00304 //    GENERAL RADIX ONE DIMENSIONAL FOURIER TRANSFORM;
00305 
00306    REPORT
00307 
00308    int  M = N;
00309 
00310    for (int i = 0; i < N_FACTOR; i++)
00311    {
00312       int P = FACTOR[i]; M /= P;
00313 
00314       switch(P)
00315       {
00316       case 1: REPORT break;
00317       case 2: REPORT R_2_FTK (N,M,X,Y,X+M,Y+M); break;
00318       case 3: REPORT R_3_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M); break;
00319       case 4: REPORT R_4_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M); break;
00320       case 5:
00321          REPORT
00322          R_5_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M,X+4*M,Y+4*M);
00323          break;
00324       case 8:
00325          REPORT
00326          R_8_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00327             X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00328             X+6*M,Y+6*M,X+7*M,Y+7*M);
00329          break;
00330       case 16:
00331          REPORT
00332          R_16_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00333             X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00334             X+6*M,Y+6*M,X+7*M,Y+7*M,X+8*M,Y+8*M,
00335             X+9*M,Y+9*M,X+10*M,Y+10*M,X+11*M,Y+11*M,
00336             X+12*M,Y+12*M,X+13*M,Y+13*M,X+14*M,Y+14*M,
00337             X+15*M,Y+15*M);
00338          break;
00339       default: REPORT R_P_FTK (N,M,P,X,Y); break;
00340       }
00341    }
00342 
00343 }
00344 
00345 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y)
00346 //    RADIX PRIME FOURIER TRANSFORM KERNEL;
00347 // X and Y are treated as M * P matrices with Fortran storage
00348 {
00349    REPORT
00350    bool NO_FOLD,ZERO;
00351    Real ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT;
00352    int  J,JJ,K0,K,M_OVER_2,MP,PM,PP,U,V;
00353 
00354    Real AA [9][9], BB [9][9];
00355    Real A [18], B [18], C [18], S [18];
00356    Real IA [9], IB [9], RA [9], RB [9];
00357 
00358    TWOPI=8.0*atan(1.0);
00359    M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;
00360 
00361    for (U=0; U<PP; U++)
00362    {
00363       ANGLE=TWOPI*Real(U+1)/Real(P);
00364       JJ=P-U-2;
00365       A[U]=cos(ANGLE); B[U]=sin(ANGLE);
00366       A[JJ]=A[U]; B[JJ]= -B[U];
00367    }
00368 
00369    for (U=1; U<=PP; U++)
00370    {
00371       for (V=1; V<=PP; V++)
00372          { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
00373    }
00374 
00375    for (J=0; J<M_OVER_2; J++)
00376    {
00377       NO_FOLD = (J==0 || 2*J==M);
00378       K0=J;
00379       ANGLE=TWOPI*Real(J)/Real(MP); ZERO=ANGLE==0.0;
00380       C[0]=cos(ANGLE); S[0]=sin(ANGLE);
00381       for (U=1; U<PM; U++)
00382       {
00383          C[U]=C[U-1]*C[0]-S[U-1]*S[0];
00384          S[U]=S[U-1]*C[0]+C[U-1]*S[0];
00385       }
00386       goto L700;
00387    L500:
00388       REPORT
00389       if (NO_FOLD) { REPORT goto L1500; }
00390       REPORT
00391       NO_FOLD=true; K0=M-J;
00392       for (U=0; U<PM; U++)
00393          { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
00394    L700:
00395       REPORT
00396       for (K=K0; K<N; K+=MP)
00397       {
00398          XT=X[K]; YT=Y[K];
00399          for (U=1; U<=PP; U++)
00400          {
00401             RA[U-1]=XT; IA[U-1]=YT;
00402             RB[U-1]=0.0; IB[U-1]=0.0;
00403          }
00404          for (U=1; U<=PP; U++)
00405          {
00406             JJ=P-U;
00407             RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
00408             RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
00409             XT=XT+RS; YT=YT+IS;
00410             for (V=0; V<PP; V++)
00411             {
00412                RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
00413                RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
00414             }
00415          }
00416          X[K]=XT; Y[K]=YT;
00417          for (U=1; U<=PP; U++)
00418          {
00419             if (!ZERO)
00420             {
00421                REPORT
00422                XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
00423                X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
00424                JJ=P-U;
00425                XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
00426                X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
00427                Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
00428             }
00429             else
00430             {
00431                REPORT
00432                X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
00433                JJ=P-U;
00434                X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
00435             }
00436          }
00437       }
00438       goto L500;
00439 L1500: ;
00440    }
00441    return;
00442 }
00443 
00444 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1)
00445 //    RADIX TWO FOURIER TRANSFORM KERNEL;
00446 {
00447    REPORT
00448    bool NO_FOLD,ZERO;
00449    int  J,K,K0,M2,M_OVER_2;
00450    Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;
00451 
00452    M2=M*2; M_OVER_2=M/2+1;
00453    TWOPI=8.0*atan(1.0);
00454 
00455    for (J=0; J<M_OVER_2; J++)
00456    {
00457       NO_FOLD = (J==0 || 2*J==M);
00458       K0=J;
00459       ANGLE=TWOPI*Real(J)/Real(M2); ZERO=ANGLE==0.0;
00460       C=cos(ANGLE); S=sin(ANGLE);
00461       goto L200;
00462    L100:
00463       REPORT
00464       if (NO_FOLD) { REPORT goto L600; }
00465       REPORT
00466       NO_FOLD=true; K0=M-J; C= -C;
00467    L200:
00468       REPORT
00469       for (K=K0; K<N; K+=M2)
00470       {
00471          RS=X0[K]+X1[K]; IS=Y0[K]+Y1[K];
00472          RU=X0[K]-X1[K]; IU=Y0[K]-Y1[K];
00473          X0[K]=RS; Y0[K]=IS;
00474          if (!ZERO) { X1[K]=RU*C+IU*S; Y1[K]=IU*C-RU*S; }
00475          else { X1[K]=RU; Y1[K]=IU; }
00476       }
00477       goto L100;
00478    L600: ;
00479    }
00480 
00481    return;
00482 }
00483 
00484 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00485    Real* X2, Real* Y2)
00486 //    RADIX THREE FOURIER TRANSFORM KERNEL
00487 {
00488    REPORT
00489    bool NO_FOLD,ZERO;
00490    int  J,K,K0,M3,M_OVER_2;
00491    Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
00492    Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;
00493 
00494    M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
00495    A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);
00496 
00497    for (J=0; J<M_OVER_2; J++)
00498    {
00499       NO_FOLD = (J==0 || 2*J==M);
00500       K0=J;
00501       ANGLE=TWOPI*Real(J)/Real(M3); ZERO=ANGLE==0.0;
00502       C1=cos(ANGLE); S1=sin(ANGLE);
00503       C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00504       goto L200;
00505    L100:
00506       REPORT
00507       if (NO_FOLD) { REPORT goto L600; }
00508       REPORT
00509       NO_FOLD=true; K0=M-J;
00510       T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
00511       T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
00512    L200:
00513       REPORT
00514       for (K=K0; K<N; K+=M3)
00515       {
00516          R0 = X0[K]; I0 = Y0[K];
00517          RS=X1[K]+X2[K]; IS=Y1[K]+Y2[K];
00518          X0[K]=R0+RS; Y0[K]=I0+IS;
00519          RA=R0+RS*A; IA=I0+IS*A;
00520          RB=(X1[K]-X2[K])*B; IB=(Y1[K]-Y2[K])*B;
00521          if (!ZERO)
00522          {
00523             REPORT
00524             R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
00525             X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00526             X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00527          }
00528          else { REPORT X1[K]=RA+IB; Y1[K]=IA-RB; X2[K]=RA-IB; Y2[K]=IA+RB; }
00529       }
00530       goto L100;
00531    L600: ;
00532    }
00533 
00534    return;
00535 }
00536 
00537 static void R_4_FTK (int N, int M,
00538    Real* X0, Real* Y0, Real* X1, Real* Y1,
00539    Real* X2, Real* Y2, Real* X3, Real* Y3)
00540 //    RADIX FOUR FOURIER TRANSFORM KERNEL
00541 {
00542    REPORT
00543    bool NO_FOLD,ZERO;
00544    int  J,K,K0,M4,M_OVER_2;
00545    Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
00546    Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;
00547 
00548    M4=M*4; M_OVER_2=M/2+1;
00549    TWOPI=8.0*atan(1.0);
00550 
00551    for (J=0; J<M_OVER_2; J++)
00552    {
00553       NO_FOLD = (J==0 || 2*J==M);
00554       K0=J;
00555       ANGLE=TWOPI*Real(J)/Real(M4); ZERO=ANGLE==0.0;
00556       C1=cos(ANGLE); S1=sin(ANGLE);
00557       C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00558       C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00559       goto L200;
00560    L100:
00561       REPORT
00562       if (NO_FOLD) { REPORT goto L600; }
00563       REPORT
00564       NO_FOLD=true; K0=M-J;
00565       T=C1; C1=S1; S1=T;
00566       C2= -C2;
00567       T=C3; C3= -S3; S3= -T;
00568    L200:
00569       REPORT
00570       for (K=K0; K<N; K+=M4)
00571       {
00572          RS0=X0[K]+X2[K]; IS0=Y0[K]+Y2[K];
00573          RU0=X0[K]-X2[K]; IU0=Y0[K]-Y2[K];
00574          RS1=X1[K]+X3[K]; IS1=Y1[K]+Y3[K];
00575          RU1=X1[K]-X3[K]; IU1=Y1[K]-Y3[K];
00576          X0[K]=RS0+RS1; Y0[K]=IS0+IS1;
00577          if (!ZERO)
00578          {
00579             REPORT
00580             R1=RU0+IU1; I1=IU0-RU1;
00581             R2=RS0-RS1; I2=IS0-IS1;
00582             R3=RU0-IU1; I3=IU0+RU1;
00583             X2[K]=R1*C1+I1*S1; Y2[K]=I1*C1-R1*S1;
00584             X1[K]=R2*C2+I2*S2; Y1[K]=I2*C2-R2*S2;
00585             X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00586          }
00587          else
00588          {
00589             REPORT
00590             X2[K]=RU0+IU1; Y2[K]=IU0-RU1;
00591             X1[K]=RS0-RS1; Y1[K]=IS0-IS1;
00592             X3[K]=RU0-IU1; Y3[K]=IU0+RU1;
00593          }
00594       }
00595       goto L100;
00596    L600: ;
00597    }
00598 
00599    return;
00600 }
00601 
00602 static void R_5_FTK (int N, int M,
00603    Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00604    Real* X3, Real* Y3, Real* X4, Real* Y4)
00605 //    RADIX FIVE FOURIER TRANSFORM KERNEL
00606 
00607 {
00608    REPORT
00609    bool NO_FOLD,ZERO;
00610    int  J,K,K0,M5,M_OVER_2;
00611    Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
00612    Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
00613    Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;
00614 
00615    M5=M*5; M_OVER_2=M/2+1;
00616    TWOPI=8.0*atan(1.0);
00617    A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
00618    A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);
00619 
00620    for (J=0; J<M_OVER_2; J++)
00621    {
00622       NO_FOLD = (J==0 || 2*J==M);
00623       K0=J;
00624       ANGLE=TWOPI*Real(J)/Real(M5); ZERO=ANGLE==0.0;
00625       C1=cos(ANGLE); S1=sin(ANGLE);
00626       C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00627       C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00628       C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00629       goto L200;
00630    L100:
00631       REPORT
00632       if (NO_FOLD) { REPORT goto L600; }
00633       REPORT
00634       NO_FOLD=true; K0=M-J;
00635       T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
00636       T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
00637       T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
00638       T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
00639    L200:
00640       REPORT
00641       for (K=K0; K<N; K+=M5)
00642       {
00643          R0=X0[K]; I0=Y0[K];
00644          RS1=X1[K]+X4[K]; IS1=Y1[K]+Y4[K];
00645          RU1=X1[K]-X4[K]; IU1=Y1[K]-Y4[K];
00646          RS2=X2[K]+X3[K]; IS2=Y2[K]+Y3[K];
00647          RU2=X2[K]-X3[K]; IU2=Y2[K]-Y3[K];
00648          X0[K]=R0+RS1+RS2; Y0[K]=I0+IS1+IS2;
00649          RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
00650          RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
00651          RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
00652          RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
00653          if (!ZERO)
00654          {
00655             REPORT
00656             R1=RA1+IB1; I1=IA1-RB1;
00657             R2=RA2+IB2; I2=IA2-RB2;
00658             R3=RA2-IB2; I3=IA2+RB2;
00659             R4=RA1-IB1; I4=IA1+RB1;
00660             X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00661             X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00662             X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00663             X4[K]=R4*C4+I4*S4; Y4[K]=I4*C4-R4*S4;
00664          }
00665          else
00666          {
00667             REPORT
00668             X1[K]=RA1+IB1; Y1[K]=IA1-RB1;
00669             X2[K]=RA2+IB2; Y2[K]=IA2-RB2;
00670             X3[K]=RA2-IB2; Y3[K]=IA2+RB2;
00671             X4[K]=RA1-IB1; Y4[K]=IA1+RB1;
00672          }
00673       }
00674       goto L100;
00675    L600: ;
00676    }
00677 
00678    return;
00679 }
00680 
00681 static void R_8_FTK (int N, int M,
00682    Real* X0, Real* Y0, Real* X1, Real* Y1,
00683    Real* X2, Real* Y2, Real* X3, Real* Y3,
00684    Real* X4, Real* Y4, Real* X5, Real* Y5,
00685    Real* X6, Real* Y6, Real* X7, Real* Y7)
00686 //    RADIX EIGHT FOURIER TRANSFORM KERNEL
00687 {
00688    REPORT
00689    bool NO_FOLD,ZERO;
00690    int  J,K,K0,M8,M_OVER_2;
00691    Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
00692    Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
00693    Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
00694    Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
00695    Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;
00696 
00697    M8=M*8; M_OVER_2=M/2+1;
00698    TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);
00699 
00700    for (J=0;J<M_OVER_2;J++)
00701    {
00702       NO_FOLD= (J==0 || 2*J==M);
00703       K0=J;
00704       ANGLE=TWOPI*Real(J)/Real(M8); ZERO=ANGLE==0.0;
00705       C1=cos(ANGLE); S1=sin(ANGLE);
00706       C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00707       C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00708       C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00709       C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00710       C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00711       C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00712       goto L200;
00713    L100:
00714       REPORT
00715       if (NO_FOLD) { REPORT goto L600; }
00716       REPORT
00717       NO_FOLD=true; K0=M-J;
00718       T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
00719       T=S2; S2=C2; C2=T;
00720       T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
00721       C4= -C4;
00722       T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
00723       T= -S6; S6= -C6; C6=T;
00724       T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
00725    L200:
00726       REPORT
00727       for (K=K0; K<N; K+=M8)
00728       {
00729          RS0=X0[K]+X4[K]; IS0=Y0[K]+Y4[K];
00730          RU0=X0[K]-X4[K]; IU0=Y0[K]-Y4[K];
00731          RS1=X1[K]+X5[K]; IS1=Y1[K]+Y5[K];
00732          RU1=X1[K]-X5[K]; IU1=Y1[K]-Y5[K];
00733          RS2=X2[K]+X6[K]; IS2=Y2[K]+Y6[K];
00734          RU2=X2[K]-X6[K]; IU2=Y2[K]-Y6[K];
00735          RS3=X3[K]+X7[K]; IS3=Y3[K]+Y7[K];
00736          RU3=X3[K]-X7[K]; IU3=Y3[K]-Y7[K];
00737          RSS0=RS0+RS2; ISS0=IS0+IS2;
00738          RSU0=RS0-RS2; ISU0=IS0-IS2;
00739          RSS1=RS1+RS3; ISS1=IS1+IS3;
00740          RSU1=RS1-RS3; ISU1=IS1-IS3;
00741          RUS0=RU0-IU2; IUS0=IU0+RU2;
00742          RUU0=RU0+IU2; IUU0=IU0-RU2;
00743          RUS1=RU1-IU3; IUS1=IU1+RU3;
00744          RUU1=RU1+IU3; IUU1=IU1-RU3;
00745          T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
00746          T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
00747          X0[K]=RSS0+RSS1; Y0[K]=ISS0+ISS1;
00748          if (!ZERO)
00749          {
00750             REPORT
00751             R1=RUU0+RUU1; I1=IUU0+IUU1;
00752             R2=RSU0+ISU1; I2=ISU0-RSU1;
00753             R3=RUS0+IUS1; I3=IUS0-RUS1;
00754             R4=RSS0-RSS1; I4=ISS0-ISS1;
00755             R5=RUU0-RUU1; I5=IUU0-IUU1;
00756             R6=RSU0-ISU1; I6=ISU0+RSU1;
00757             R7=RUS0-IUS1; I7=IUS0+RUS1;
00758             X4[K]=R1*C1+I1*S1; Y4[K]=I1*C1-R1*S1;
00759             X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00760             X6[K]=R3*C3+I3*S3; Y6[K]=I3*C3-R3*S3;
00761             X1[K]=R4*C4+I4*S4; Y1[K]=I4*C4-R4*S4;
00762             X5[K]=R5*C5+I5*S5; Y5[K]=I5*C5-R5*S5;
00763             X3[K]=R6*C6+I6*S6; Y3[K]=I6*C6-R6*S6;
00764             X7[K]=R7*C7+I7*S7; Y7[K]=I7*C7-R7*S7;
00765          }
00766          else
00767          {
00768             REPORT
00769             X4[K]=RUU0+RUU1; Y4[K]=IUU0+IUU1;
00770             X2[K]=RSU0+ISU1; Y2[K]=ISU0-RSU1;
00771             X6[K]=RUS0+IUS1; Y6[K]=IUS0-RUS1;
00772             X1[K]=RSS0-RSS1; Y1[K]=ISS0-ISS1;
00773             X5[K]=RUU0-RUU1; Y5[K]=IUU0-IUU1;
00774             X3[K]=RSU0-ISU1; Y3[K]=ISU0+RSU1;
00775             X7[K]=RUS0-IUS1; Y7[K]=IUS0+RUS1;
00776          }
00777       }
00778       goto L100;
00779    L600: ;
00780    }
00781 
00782    return;
00783 }
00784 
00785 static void R_16_FTK (int N, int M,
00786    Real* X0, Real* Y0, Real* X1, Real* Y1,
00787    Real* X2, Real* Y2, Real* X3, Real* Y3,
00788    Real* X4, Real* Y4, Real* X5, Real* Y5,
00789    Real* X6, Real* Y6, Real* X7, Real* Y7,
00790    Real* X8, Real* Y8, Real* X9, Real* Y9,
00791    Real* X10, Real* Y10, Real* X11, Real* Y11,
00792    Real* X12, Real* Y12, Real* X13, Real* Y13,
00793    Real* X14, Real* Y14, Real* X15, Real* Y15)
00794 //    RADIX SIXTEEN FOURIER TRANSFORM KERNEL
00795 {
00796    REPORT
00797    bool NO_FOLD,ZERO;
00798    int  J,K,K0,M16,M_OVER_2;
00799    Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
00800    Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
00801    Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
00802    Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
00803    Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
00804    Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
00805    Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
00806    Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
00807    Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
00808    Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
00809    Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
00810    Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
00811    Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
00812    Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
00813    Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
00814    Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
00815    Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;
00816 
00817    M16=M*16; M_OVER_2=M/2+1;
00818    TWOPI=8.0*atan(1.0);
00819    ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
00820    E2=cos(TWOPI/8.0);
00821    ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
00822    ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);
00823 
00824    for (J=0; J<M_OVER_2; J++)
00825    {
00826       NO_FOLD = (J==0 || 2*J==M);
00827       K0=J;
00828       ANGLE=TWOPI*Real(J)/Real(M16);
00829       ZERO=ANGLE==0.0;
00830       C1=cos(ANGLE); S1=sin(ANGLE);
00831       C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00832       C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00833       C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00834       C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00835       C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00836       C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00837       C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
00838       C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
00839       C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
00840       C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
00841       C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
00842       C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
00843       C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
00844       C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
00845       goto L200;
00846    L100:
00847       REPORT
00848       if (NO_FOLD) { REPORT goto L600; }
00849       REPORT
00850       NO_FOLD=true; K0=M-J;
00851       T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
00852       T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
00853       T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
00854       T=S4; S4=C4; C4=T;
00855       T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
00856       T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
00857       T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
00858       C8= -C8;
00859       T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
00860       T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
00861       T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
00862       T= -S12; S12= -C12; C12=T;
00863       T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
00864       T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
00865       T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
00866    L200:
00867       REPORT
00868       for (K=K0; K<N; K+=M16)
00869       {
00870          RS0=X0[K]+X8[K]; IS0=Y0[K]+Y8[K];
00871          RU0=X0[K]-X8[K]; IU0=Y0[K]-Y8[K];
00872          RS1=X1[K]+X9[K]; IS1=Y1[K]+Y9[K];
00873          RU1=X1[K]-X9[K]; IU1=Y1[K]-Y9[K];
00874          RS2=X2[K]+X10[K]; IS2=Y2[K]+Y10[K];
00875          RU2=X2[K]-X10[K]; IU2=Y2[K]-Y10[K];
00876          RS3=X3[K]+X11[K]; IS3=Y3[K]+Y11[K];
00877          RU3=X3[K]-X11[K]; IU3=Y3[K]-Y11[K];
00878          RS4=X4[K]+X12[K]; IS4=Y4[K]+Y12[K];
00879          RU4=X4[K]-X12[K]; IU4=Y4[K]-Y12[K];
00880          RS5=X5[K]+X13[K]; IS5=Y5[K]+Y13[K];
00881          RU5=X5[K]-X13[K]; IU5=Y5[K]-Y13[K];
00882          RS6=X6[K]+X14[K]; IS6=Y6[K]+Y14[K];
00883          RU6=X6[K]-X14[K]; IU6=Y6[K]-Y14[K];
00884          RS7=X7[K]+X15[K]; IS7=Y7[K]+Y15[K];
00885          RU7=X7[K]-X15[K]; IU7=Y7[K]-Y15[K];
00886          RSS0=RS0+RS4; ISS0=IS0+IS4;
00887          RSS1=RS1+RS5; ISS1=IS1+IS5;
00888          RSS2=RS2+RS6; ISS2=IS2+IS6;
00889          RSS3=RS3+RS7; ISS3=IS3+IS7;
00890          RSU0=RS0-RS4; ISU0=IS0-IS4;
00891          RSU1=RS1-RS5; ISU1=IS1-IS5;
00892          RSU2=RS2-RS6; ISU2=IS2-IS6;
00893          RSU3=RS3-RS7; ISU3=IS3-IS7;
00894          RUS0=RU0-IU4; IUS0=IU0+RU4;
00895          RUS1=RU1-IU5; IUS1=IU1+RU5;
00896          RUS2=RU2-IU6; IUS2=IU2+RU6;
00897          RUS3=RU3-IU7; IUS3=IU3+RU7;
00898          RUU0=RU0+IU4; IUU0=IU0-RU4;
00899          RUU1=RU1+IU5; IUU1=IU1-RU5;
00900          RUU2=RU2+IU6; IUU2=IU2-RU6;
00901          RUU3=RU3+IU7; IUU3=IU3-RU7;
00902          T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
00903          T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
00904          T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
00905          T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
00906          T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
00907          T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
00908          T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
00909          T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
00910          RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
00911          RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
00912          RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
00913          RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
00914          RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
00915          RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
00916          RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
00917          RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
00918          RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
00919          RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
00920          RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
00921          RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
00922          RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
00923          RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
00924          RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
00925          RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
00926          X0[K]=RSSS0+RSSS1; Y0[K]=ISSS0+ISSS1;
00927          if (!ZERO)
00928          {
00929             REPORT
00930             R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
00931             R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
00932             R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
00933             R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
00934             R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
00935             R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
00936             R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
00937             R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
00938             R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
00939             R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
00940             R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
00941             R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
00942             R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
00943             R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
00944             R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
00945             X8[K]=R1*C1+I1*S1; Y8[K]=I1*C1-R1*S1;
00946             X4[K]=R2*C2+I2*S2; Y4[K]=I2*C2-R2*S2;
00947             X12[K]=R3*C3+I3*S3; Y12[K]=I3*C3-R3*S3;
00948             X2[K]=R4*C4+I4*S4; Y2[K]=I4*C4-R4*S4;
00949             X10[K]=R5*C5+I5*S5; Y10[K]=I5*C5-R5*S5;
00950             X6[K]=R6*C6+I6*S6; Y6[K]=I6*C6-R6*S6;
00951             X14[K]=R7*C7+I7*S7; Y14[K]=I7*C7-R7*S7;
00952             X1[K]=R8*C8+I8*S8; Y1[K]=I8*C8-R8*S8;
00953             X9[K]=R9*C9+I9*S9; Y9[K]=I9*C9-R9*S9;
00954             X5[K]=R10*C10+I10*S10; Y5[K]=I10*C10-R10*S10;
00955             X13[K]=R11*C11+I11*S11; Y13[K]=I11*C11-R11*S11;
00956             X3[K]=R12*C12+I12*S12; Y3[K]=I12*C12-R12*S12;
00957             X11[K]=R13*C13+I13*S13; Y11[K]=I13*C13-R13*S13;
00958             X7[K]=R14*C14+I14*S14; Y7[K]=I14*C14-R14*S14;
00959             X15[K]=R15*C15+I15*S15; Y15[K]=I15*C15-R15*S15;
00960          }
00961          else
00962          {
00963             REPORT
00964             X8[K]=RUUS0+RUUS1; Y8[K]=IUUS0+IUUS1;
00965             X4[K]=RSUU0+RSUU1; Y4[K]=ISUU0+ISUU1;
00966             X12[K]=RUSU0+RUSU1; Y12[K]=IUSU0+IUSU1;
00967             X2[K]=RSSU0+ISSU1; Y2[K]=ISSU0-RSSU1;
00968             X10[K]=RUUU0+IUUU1; Y10[K]=IUUU0-RUUU1;
00969             X6[K]=RSUS0+ISUS1; Y6[K]=ISUS0-RSUS1;
00970             X14[K]=RUSS0+IUSS1; Y14[K]=IUSS0-RUSS1;
00971             X1[K]=RSSS0-RSSS1; Y1[K]=ISSS0-ISSS1;
00972             X9[K]=RUUS0-RUUS1; Y9[K]=IUUS0-IUUS1;
00973             X5[K]=RSUU0-RSUU1; Y5[K]=ISUU0-ISUU1;
00974             X13[K]=RUSU0-RUSU1; Y13[K]=IUSU0-IUSU1;
00975             X3[K]=RSSU0-ISSU1; Y3[K]=ISSU0+RSSU1;
00976             X11[K]=RUUU0-IUUU1; Y11[K]=IUUU0+RUUU1;
00977             X7[K]=RSUS0-ISUS1; Y7[K]=ISUS0+RSUS1;
00978             X15[K]=RUSS0-IUSS1; Y15[K]=IUSS0+RUSS1;
00979          }
00980       }
00981       goto L100;
00982    L600: ;
00983    }
00984 
00985    return;
00986 }
00987 
00988 // can the number of points be factorised sufficiently
00989 // for the fft to run
00990 
00991 bool FFT_Controller::CanFactor(int PTS)
00992 {
00993    REPORT
00994    const int NP = 16, NQ = 10, PMAX=19;
00995 
00996    if (PTS<=1) { REPORT return true; }
00997 
00998    int N = PTS, F = 2, P = 0, Q = 0;
00999 
01000    while (N > 1)
01001    {
01002       bool fail = true;
01003       for (int J = F; J <= PMAX; J++)
01004          if (N % J == 0) { fail = false; F=J; break; }
01005       if (fail || P >= NP || Q >= NQ) { REPORT return false; }
01006       N /= F;
01007       if (N % F != 0) Q++; else { N /= F; P++; }
01008    }
01009 
01010    return true;    // can factorise
01011 
01012 }
01013 
01014 bool FFT_Controller::OnlyOldFFT;         // static variable
01015 
01016 // **************************** multi radix counter **********************
01017 
01018 MultiRadixCounter::MultiRadixCounter(int nx, const SimpleIntArray& rx,
01019    SimpleIntArray& vx)
01020    : Radix(rx), Value(vx), n(nx), reverse(0),
01021       product(1), counter(0), finish(false)
01022 {
01023    REPORT for (int k = 0; k < n; k++) { Value[k] = 0; product *= Radix[k]; }
01024 }
01025 
01026 void MultiRadixCounter::operator++()
01027 {
01028    REPORT
01029    counter++; int p = product;
01030    for (int k = 0; k < n; k++)
01031    {
01032       Value[k]++; int p1 = p / Radix[k]; reverse += p1;
01033       if (Value[k] == Radix[k]) { REPORT Value[k] = 0; reverse -= p; p = p1; }
01034       else { REPORT return; }
01035    }
01036    finish = true;
01037 }
01038 
01039 
01040 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f)
01041 {
01042    // x = c[0]+f[0]*(c[1]+f[1]*(c[2]+...
01043    // return c[n-1]+f[n-1]*(c[n-2]+f[n-2]*(c[n-3]+...
01044    // prod is the product of the f[i]
01045    // n is the number of f[i] (don't assume f has the correct length)
01046 
01047    REPORT
01048    const int* d = f.Data() + n; int sum = 0; int q = 1;
01049    while (n--)
01050    {
01051       prod /= *(--d);
01052       int c = x / prod; x-= c * prod;
01053       sum += q * c; q *= *d;
01054    }
01055    return sum;
01056 }
01057 
01058 
01059 #ifdef use_namespace
01060 }
01061 #endif
01062 
01063 

newmat11b
Generated Mon May 9 04:54:17 2016 by Doxygen 1.6.3