00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #undef STANDALONE  /* Define this if you want to use this code */
00008                    /* outside of the AFNI package (libmri.a).  */
00009 
00010 #ifdef STANDALONE
00011 #  include <stdlib.h>   /* include standard system headers */
00012 #  include <stddef.h>
00013 #  include <stdio.h>
00014 #  include <math.h>
00015    typedef struct complex { float r,i ; } complex ;  /* need this! */
00016 #  undef USE_FFTW
00017 #  define EXIT exit
00018 #else
00019 #  include "mrilib.h"   /* AFNI package library header */
00020 
00021 #  ifdef USE_FFTW         /* 20 Oct 2000 */
00022 #    include "fftw.h"
00023 #  endif
00024 
00025 #endif  /* STANDALONE */
00026 
00027 static int use_fftw=1 ;
00028 void csfft_use_fftw( int uf ){ use_fftw = uf; return; }
00029 
00030 /***=====================================================================**
00031  *** Prototypes for externally callable functions:                       **
00032  ***    Complex-to-complex FFT in place:                                 **
00033  ***      mode = -1 or +1 (for inverse scaling, cf. csfft_scale_inverse  **
00034  ***      idim = dimension (power of 2, maybe with factors of 3 and 5)   **
00035  ***      xc   = input/output array                                      **
00036  ***    Re-initializes itself only when idim changes from previous call. **
00037  ***=====================================================================**/
00038 
00039 void csfft_cox( int mode , int idim , complex * xc ) ;
00040 int  csfft_nextup( int idim ) ;
00041 int  csfft_nextup_one35( int idim ) ;
00042 void csfft_scale_inverse( int scl ) ; /* scl=1 ==> force 1/N for mode=+1 **/
00043                                       /* scl=0 ==> no 1/N scaling        **/
00044 
00045 /*** Aug 1999:                                                           **
00046  ***   idim can now contain factors of 3 and/or 5, up to and including   **
00047  ***   3^3 * 5^3.  Some examples:                                        **
00048  ***       135 144 150 160 180 192 200 216 225 240 250                   **
00049  ***   The routine csfft_nextup(n) returns the smallest size FFT         **
00050  ***    >= n that csfft_cox() knows how to do.                           **
00051  ***   Note that efficiency of the different lengths can be quite        **
00052  ***   different.  In general, powers of 2 are fastest, but a            **
00053  ***   single factor of 3 and/or 5 doesn't cause too much slowdown.      **
00054  ***   The routine csfft_nextup_one35(n) returns the smallest size FFT   **
00055  ***    >= n that contains at most one power of 3, at most one power     **
00056  ***   of 5, and least one power of 2.  In trials, these are the most    **
00057  ***   time efficient.                                                   **/
00058 
00059 /*-- Aug 1999: routines to do FFTs by decimation by 3 or 5 --*/
00060 
00061 static void fft_3dec( int , int , complex * ) ;
00062 static void fft_5dec( int , int , complex * ) ;
00063 
00064 #undef  RMAX
00065 #define RMAX 3   /* maximum levels of recursion for radix-3 and radix-5 */
00066 
00067 #undef  PI
00068 #define PI (3.141592653589793238462643)
00069 
00070 /*-------------- To order the 1/N scaling on inversion: Aug 1999 ------------*/
00071 
00072 static int sclinv = 0 ;  /* set by csfft_scale_inverse */
00073 
00074 void csfft_scale_inverse( int scl ){ sclinv = scl; return; }
00075 
00076 /*-------------- For the unrolled FFT routines: November 1998 --------------*/
00077 
00078 #ifndef DONT_UNROLL_FFTS
00079    static void fft8  ( int mode , complex * xc ) ; /* completely */
00080    static void fft16 ( int mode , complex * xc ) ; /* unrolled   */
00081    static void fft32 ( int mode , complex * xc ) ; /* routines   */
00082 
00083    static void fft64 ( int mode , complex * xc ) ; /* internal 2dec ->  32 */
00084    static void fft128( int mode , complex * xc ) ; /* internal 4dec ->  32 */
00085    static void fft256( int mode , complex * xc ) ; /* internal 4dec ->  64 */
00086    static void fft512( int mode , complex * xc ) ; /* internal 4dec -> 128 */
00087 
00088    static void fft_4dec( int , int , complex * ) ; /* general 4dec */
00089 
00090 #  define fft1024(m,x) fft_4dec(m,1024,x)          /* 4dec -> 256 */
00091 #  define fft2048(m,x) fft_4dec(m,2048,x)          /* 4dec -> 512 */
00092 
00093 /*-- do FFTs of size 2 and 4 with macros --*/
00094 
00095 #  define fft2(m,x) do{ float a,b,c,d ;                             \
00096                         a = x[0].r + x[1].r ; b = x[0].i + x[1].i ; \
00097                         c = x[0].r - x[1].r ; d = x[0].i - x[1].i ; \
00098                         x[0].r = a ; x[0].i = b ;                   \
00099                         x[1].r = c ; x[1].i = d ; } while(0)
00100 
00101 # define USE_FFT4_MACRO
00102 # ifndef USE_FFT4_MACRO
00103    static void fft4( int mode , complex * xc ) ; /* unrolled routine */
00104 # else
00105 #  define fft4(m,x) do{ float acpr,acmr,bdpr,bdmr, acpi,acmi,bdpi,bdmi; \
00106                         acpr = x[0].r + x[2].r; acmr = x[0].r - x[2].r; \
00107                         bdpr = x[1].r + x[3].r; bdmr = x[1].r - x[3].r; \
00108                         acpi = x[0].i + x[2].i; acmi = x[0].i - x[2].i; \
00109                         bdpi = x[1].i + x[3].i; bdmi = x[1].i - x[3].i; \
00110                         x[0].r = acpr+bdpr ; x[0].i = acpi+bdpi ;       \
00111                         x[2].r = acpr-bdpr ; x[2].i = acpi-bdpi ;       \
00112                         if(m > 0){                                      \
00113                           x[1].r = acmr-bdmi ; x[1].i = acmi+bdmr ;     \
00114                           x[3].r = acmr+bdmi ; x[3].i = acmi-bdmr ;     \
00115                         } else {                                        \
00116                           x[1].r = acmr+bdmi ; x[1].i = acmi-bdmr ;     \
00117                           x[3].r = acmr-bdmi ; x[3].i = acmi+bdmr ;     \
00118                         } } while(0)
00119 # endif
00120 #endif
00121 
00122 /*----------------------------------------------------------------------
00123    Speedups with unrolled FFTs [program fftest.c]:
00124    Pentium II 400 MHz:  58% (32) to 36% (256)  [gcc -O3 -ffast-math]
00125    SGI R10000 175 MHz:  47% (32) to 18% (256)  [cc -Ofast]
00126    HP PA-8000 200 MHz:  58% (32) to 40% (256)  [cc +O3 +Oaggressive]
00127 ------------------------------------------------------------------------*/
00128 
00129 /*----------------- the csfft trig constants tables ------------------*/
00130 
00131 static complex * csplus = NULL , * csminus = NULL ;
00132 static int nold = -666 ;
00133 
00134 /*--------------------------------------------------------------------
00135    Initialize csfft trig constants table.  Adapted from AJ's code.
00136    Does both mode=1 and mode=-1 tables.
00137 ----------------------------------------------------------------------*/
00138 
00139 static void csfft_trigconsts( int idim )  /* internal function */
00140 {
00141    register unsigned int  m, n, i0, i1, i2, i3, k;
00142    register complex       *r0, *csp;
00143    register float         co, si, f0, f1, f2, f3, f4;
00144    double                 al;
00145 
00146    if( idim == nold ) return ;
00147 
00148    if( idim > nold ){
00149       if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
00150       csplus  = (complex *) malloc( sizeof(complex) * idim ) ;
00151       csminus = (complex *) malloc( sizeof(complex) * idim ) ;
00152       if( csplus == NULL || csminus == NULL ){
00153          fprintf(stderr,"\n*** csfft cannot malloc space! ***\n"); EXIT(1) ;
00154       }
00155    }
00156    nold = n = idim ;
00157 
00158    f1 = 1.0 ;  /* csplus init */
00159    m  = 1; k  = 0;
00160    while (n > m) {
00161       i3 = m << 1; f2 = m; al = f1*PI/f2;
00162       co = cos(al); si = sin(al);
00163       (csplus + k)->r = 1.; (csplus + k)->i = 0.;
00164       for (i0=0; i0 < m; i0++) {
00165          k++;
00166          csp = csplus + k; r0 = csp - 1;
00167          csp->r = r0->r * co - r0->i * si;
00168          csp->i = r0->i * co + r0->r * si;
00169       }
00170       m = i3;
00171    }
00172 
00173    f1 = -1.0 ;  /* csminus init */
00174    m  = 1; k  = 0;
00175    while (n > m) {
00176       i3 = m << 1; f2 = m; al = f1*PI/f2;
00177       co = cos(al); si = sin(al);
00178       (csminus + k)->r = 1.; (csminus + k)->i = 0.;
00179       for (i0=0; i0 < m; i0++) {
00180          k++;
00181          csp = csminus + k; r0  = csp - 1;
00182          csp->r = r0->r * co - r0->i * si;
00183          csp->i = r0->i * co + r0->r * si;
00184       }
00185       m = i3;
00186    }
00187    return ;
00188 }
00189 
00190 /*--------------------------------------------------------------------
00191    Complex-to-complex FFT in place:
00192      mode = -1 or +1 (1/idim scaling on inverse [+1] if sclinv
00193                       was set in the call to csfft_scale_inverse)
00194      idim = dimension (power of 2, possibly with a factor
00195                        of 3^p * 5^q also, for p,q=0..RMAX)
00196      xc   = input/output array
00197    Automagically re-initializes itself when idim changes from
00198    previous call.  By AJ (Andrzej Jesmanowicz), modified by RWCox.
00199 ----------------------------------------------------------------------*/
00200 
00201 /*- Macro to do 1/N scaling on inverse:
00202       if it is ordered by the user, and
00203       if mode is positive, and
00204       if this is not a recursive call.  -*/
00205 
00206 #define SCLINV                                                 \
00207  if( sclinv && mode > 0 && rec == 0 ){                         \
00208    register int qq ; register float ff = 1.0 / (float) idim ;  \
00209    for( qq=0 ; qq < idim ; qq++ ){ xc[qq].r *= ff ; xc[qq].i *= ff ; } }
00210 
00211 void csfft_cox( int mode , int idim , complex * xc )
00212 {
00213    register unsigned int  m, n, i0, i1, i2, i3, k;
00214    register complex       *r0, *r1, *csp;
00215    register float         co, si, f0, f1, f2, f3, f4;
00216 
00217    static int rec=0 ;  /* recursion level */
00218 
00219    /*-- 20 Oct 2000: maybe use FFTW instead --*/
00220 
00221 #ifdef USE_FFTW
00222    if( use_fftw ){
00223      static int first=1 , oldmode=0, oldidim=0 ;
00224      static fftw_plan fpl ;
00225      if( first ){
00226         char * env = getenv("AFNI_FFTW_WISDOM") ;
00227         int gotit=0 ;
00228         first = 0 ;
00229         if( env != NULL ){
00230            int nw = THD_filesize(env) ;
00231            if( nw > 0 ){
00232               int fd = open( env , O_RDONLY ) ;
00233               if( fd >= 0 ){
00234                  char * w = AFMALL( char, nw+4 ) ;
00235                  int ii = read( fd , w , nw ) ;
00236                  close(fd) ;
00237                  if( ii > 0 ){
00238                     w[nw] = '\0' ;
00239                     nw = (int) fftw_import_wisdom_from_string(w) ;
00240                     gotit=(nw == FFTW_SUCCESS);
00241                  } /* else fprintf(stderr,"csfft_cox: read fails\n"); */
00242                  free(w) ;
00243               } /* else fprintf(stderr,"csfft_cox: open fails\n"); */
00244            } /* else fprintf(stderr,"csfft_cox: THD_filesize(%s) fails\n",env); */
00245         } /* else fprintf(stderr,"csfft_cox: getenv fails\n"); */
00246      /* if( gotit )
00247           fprintf(stderr,"csfft_cox imported FFTW wisdom from file %s\n",env);
00248         else
00249           fprintf(stderr,"csfft_cox failed to import FFTW wisdom\n"); */
00250      }
00251 
00252      if( idim != oldidim || mode != oldmode ){
00253         if( oldidim ) fftw_destroy_plan(fpl) ;
00254         fpl = fftw_create_plan( idim ,
00255                                 (mode < 0) ? FFTW_FORWARD : FFTW_BACKWARD ,
00256                                 FFTW_ESTIMATE|FFTW_USE_WISDOM|FFTW_IN_PLACE );
00257         oldidim = idim ; oldmode = mode ;
00258      }
00259 
00260      fftw_one( fpl , (fftw_complex *)xc , NULL ) ;
00261      SCLINV ; return ;
00262    }
00263 #endif /* USE_FFTW */
00264 
00265    /*-- November 1998: maybe use the unrolled FFT routines --*/
00266 
00267 #ifndef DONT_UNROLL_FFTS
00268    switch( idim ){                                  /* none of these  */
00269       case    1:                           return;  /* routines will  */
00270       case    2: fft2   (mode,xc); SCLINV; return;  /* call csfft_cox */
00271       case    4: fft4   (mode,xc); SCLINV; return;  /* so don't need  */
00272       case    8: fft8   (mode,xc); SCLINV; return;  /* to do rec++/-- */
00273       case   16: fft16  (mode,xc); SCLINV; return;
00274       case   32: fft32  (mode,xc); SCLINV; return;
00275       case   64: fft64  (mode,xc); SCLINV; return;
00276       case  128: fft128 (mode,xc); SCLINV; return;
00277       case  256: fft256 (mode,xc); SCLINV; return;
00278       case  512: fft512 (mode,xc); SCLINV; return;
00279       case 1024: fft1024(mode,xc); SCLINV; return;
00280       case 2048: fft2048(mode,xc); SCLINV; return;
00281 
00282       case  4096: fft_4dec(mode, 4096,xc); SCLINV; return; /* 4dec -> 1024 */
00283       case  8192: fft_4dec(mode, 8192,xc); SCLINV; return; /* 4dec -> 2048 */
00284       case 16384: fft_4dec(mode,16384,xc); SCLINV; return; /* 4dec -> 4096 */
00285       case 32768: fft_4dec(mode,32768,xc); SCLINV; return; /* 4dec -> 8192 */
00286    }
00287 #endif  /* end of unrollificationizing */
00288 
00289    /*-- Aug 1999: deal with factors of 3 or 5 [might call csfft_cox] --*/
00290 
00291    if( idim%3 == 0 ){rec++; fft_3dec(mode,idim,xc); rec--; SCLINV; return;}
00292    if( idim%5 == 0 ){rec++; fft_5dec(mode,idim,xc); rec--; SCLINV; return;}
00293 
00294    /*-- below here: the old general power of 2 routine --*/
00295 
00296    /**-- perhaps initialize --**/
00297 
00298    if( nold != idim ) csfft_trigconsts( idim ) ;
00299 
00300    /** Main loop starts here **/
00301 
00302    n   = idim;
00303    i2  = idim >> 1;
00304    i1  = 0;
00305    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00306 
00307    /*-- swap data --*/
00308 
00309    for (i0=0; i0 < n; i0 ++) {
00310       if ( i1 > i0 ) {
00311          r0    = xc + i0; r1    = xc + i1;
00312          f1    = r0->r;   f2    = r0->i;
00313          r0->r = r1->r;   r0->i = r1->i;
00314          r1->r = f1;      r1->i = f2;
00315       }
00316       m = i2;
00317       while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; }
00318      i1 += m;
00319    }
00320 
00321    /*-- compute compute compute --*/
00322 
00323    m = 1; k = 0;
00324    while (n > m) {
00325       i3 = m << 1;
00326       for (i0=0; i0 < m; i0 ++) {
00327          co = (csp + k)->r; si = (csp + k)->i;
00328          for (i1=i0; i1 < n; i1 += i3) {
00329             r0    = xc + i1;    r1    = r0 + m;
00330             f1    = r1->r * co; f2    = r1->i * si;
00331             f3    = r1->r * si; f4    = r1->i * co;
00332             f1   -= f2;         f3   += f4;
00333             f2    = r0->r;      f4    = r0->i;
00334             r1->r = f2 - f1;    r1->i = f4 - f3;
00335             r0->r = f2 + f1;    r0->i = f4 + f3;
00336          }
00337          k++;
00338       }
00339       m = i3;
00340    }
00341 
00342    SCLINV ; return ;
00343 }
00344 
00345 /*--------------------------------------------------------------------
00346    Many complex-to-complex FFTs in place [vectorized]:
00347      mode = -1 or +1 (NO SCALING ON INVERSE!)
00348      idim = dimension (power of 2)
00349      nvec = number of vectors to do (vectors are contiguous)
00350      xc   = input/output array
00351    Automagically re-initializes itself when idim changes from
00352    previous call.  By AJJesmanowicz, modified by RWCox.
00353 ----------------------------------------------------------------------*/
00354 
00355 void csfft_many( int mode , int idim , int nvec , complex * xc )
00356 {
00357    register unsigned int  m, n, i0, i1, i2, i3, k , iv ;
00358    register complex       *r0, *r1, *csp , *xcx;
00359    register float         co, si, f0, f1, f2, f3, f4;
00360 
00361    if( nvec == 1 ){ csfft_cox( mode , idim , xc ) ; return ; }
00362 
00363    if( idim % 3 == 0 ){                           /* Aug 1999 */
00364       for( m=0 ; m < nvec ; m++ )
00365          fft_3dec( mode , idim , xc + m*idim ) ;
00366       return ;
00367    } else if( idim % 5 == 0 ){
00368       for( m=0 ; m < nvec ; m++ )
00369          fft_5dec( mode , idim , xc + m*idim ) ;
00370       return ;
00371    }
00372 
00373    /** perhaps initialize **/
00374 
00375    if( nold != idim ) csfft_trigconsts( idim ) ;
00376 
00377    n   = idim;
00378    i2  = idim >> 1;
00379    i1  = 0;
00380    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00381 
00382    for (i0=0; i0 < n; i0 ++) {
00383       if ( i1 > i0 ) {
00384          for( iv=0,xcx=xc ; iv < nvec ; iv++,xcx+=n ){
00385             r0    = xcx + i0; r1    = xcx + i1;
00386             f1    = r0->r;    f2    = r0->i;
00387             r0->r = r1->r;    r0->i = r1->i;
00388             r1->r = f1;       r1->i = f2;
00389          }
00390       }
00391       m = i2;
00392       while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; }
00393      i1 += m;
00394    }
00395 
00396 #define I00
00397 #ifdef I00
00398 #  define I0BOT 1
00399 #else
00400 #  define I0BOT 0
00401 #endif
00402 
00403    m = 1;
00404    k = 0;
00405    while (n > m) {
00406       i3 = m << 1;
00407 
00408 #ifdef I00
00409       /* handle i0=0 case [co=1,si=0] in special code */
00410 
00411       for (i1=0; i1 < n; i1 += i3) {
00412          for( iv=0,r0=xc+i1 ; iv < nvec ; iv++,r0+=n ){
00413             r1    = r0 + m;
00414             f1    = r1->r ;   f3    = r1->i ;
00415             f2    = r0->r ;   f4    = r0->i ;
00416             r1->r = f2 - f1 ; r1->i = f4 - f3 ;
00417             r0->r = f2 + f1 ; r0->i = f4 + f3 ;
00418          }
00419       }
00420       k++;
00421 #endif
00422 
00423       for (i0=I0BOT; i0 < m; i0 ++) {
00424          co = (csp + k)->r; si = (csp + k)->i;
00425          for (i1=i0; i1 < n; i1 += i3) {
00426             for( iv=0,r0=xc+i1 ; iv < nvec ; iv++,r0+=n ){
00427                r1    = r0 + m;
00428 #if 1
00429                f1    = r1->r * co - r1->i * si ;
00430                f3    = r1->r * si + r1->i * co ;
00431 #else
00432                f1    = r1->r * co ; f2    = r1->i * si ;
00433                f3    = r1->r * si ; f4    = r1->i * co ;
00434                f1   -= f2 ;         f3   += f4 ;
00435 #endif
00436                f2    = r0->r ;   f4    = r0->i ;
00437                r1->r = f2 - f1 ; r1->i = f4 - f3 ;
00438                r0->r = f2 + f1 ; r0->i = f4 + f3 ;
00439             }
00440          }
00441          k++;
00442       }
00443       m = i3;
00444    }
00445    return ;
00446 }
00447 
00448 /****************************************************************************
00449    Everything from here on down is for the unrolled FFT routines.
00450    They were generated with program fftprint.c.
00451 *****************************************************************************/
00452 
00453 #ifndef DONT_UNROLL_FFTS
00454 
00455 #ifndef USE_FFT4_MACRO
00456 /**************************************/
00457 /* FFT routine unrolled of length   4 */
00458 
00459 static void fft4( int mode , complex * xc )
00460 {
00461    register complex * csp , * xcx=xc;
00462    register float f1,f2,f3,f4 ;
00463 
00464    /** perhaps initialize **/
00465 
00466    if( nold != 4 ) csfft_trigconsts( 4 ) ;
00467 
00468    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00469 
00470    /** data swapping part **/
00471 
00472    f1 = xcx[1].r ; f2 = xcx[1].i ;
00473    xcx[1].r = xcx[2].r ; xcx[1].i = xcx[2].i ;
00474    xcx[2].r = f1 ; xcx[2].i = f2 ;
00475 
00476    /** butterflying part **/
00477 
00478    f1 = xcx[1].r ; f3 = xcx[1].i ;  /* cos=1 sin=0 */
00479    f2 = xcx[0].r ; f4 = xcx[0].i ;
00480    xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
00481    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00482 
00483    f1 = xcx[3].r ; f3 = xcx[3].i ;  /* cos=1 sin=0 */
00484    f2 = xcx[2].r ; f4 = xcx[2].i ;
00485    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00486    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00487 
00488    f1 = xcx[2].r ; f3 = xcx[2].i ;  /* cos=1 sin=0 */
00489    f2 = xcx[0].r ; f4 = xcx[0].i ;
00490    xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
00491    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00492 
00493    f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */
00494    f3 = xcx[3].r * csp[2].i ;
00495    f2 = xcx[1].r ; f4 = xcx[1].i ;
00496    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00497    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00498 
00499    return ;
00500 }
00501 #endif /* USE_FFT4_MACRO */
00502 
00503 /**************************************/
00504 /* FFT routine unrolled of length   8 */
00505 
00506 static void fft8( int mode , complex * xc )
00507 {
00508    register complex * csp , * xcx=xc;
00509    register float f1,f2,f3,f4 ;
00510 
00511    /** perhaps initialize **/
00512 
00513    if( nold != 8 ) csfft_trigconsts( 8 ) ;
00514 
00515    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00516 
00517    /** data swapping part **/
00518 
00519    f1 = xcx[1].r ; f2 = xcx[1].i ;
00520    xcx[1].r = xcx[4].r ; xcx[1].i = xcx[4].i ;
00521    xcx[4].r = f1 ; xcx[4].i = f2 ;
00522 
00523    f1 = xcx[3].r ; f2 = xcx[3].i ;
00524    xcx[3].r = xcx[6].r ; xcx[3].i = xcx[6].i ;
00525    xcx[6].r = f1 ; xcx[6].i = f2 ;
00526 
00527    /** butterflying part **/
00528 
00529    f1 = xcx[1].r ; f3 = xcx[1].i ;  /* cos=1 sin=0 */
00530    f2 = xcx[0].r ; f4 = xcx[0].i ;
00531    xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
00532    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00533 
00534    f1 = xcx[3].r ; f3 = xcx[3].i ;  /* cos=1 sin=0 */
00535    f2 = xcx[2].r ; f4 = xcx[2].i ;
00536    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00537    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00538 
00539    f1 = xcx[5].r ; f3 = xcx[5].i ;  /* cos=1 sin=0 */
00540    f2 = xcx[4].r ; f4 = xcx[4].i ;
00541    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
00542    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00543 
00544    f1 = xcx[7].r ; f3 = xcx[7].i ;  /* cos=1 sin=0 */
00545    f2 = xcx[6].r ; f4 = xcx[6].i ;
00546    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00547    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
00548 
00549    f1 = xcx[2].r ; f3 = xcx[2].i ;  /* cos=1 sin=0 */
00550    f2 = xcx[0].r ; f4 = xcx[0].i ;
00551    xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
00552    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00553 
00554    f1 = xcx[6].r ; f3 = xcx[6].i ;  /* cos=1 sin=0 */
00555    f2 = xcx[4].r ; f4 = xcx[4].i ;
00556    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
00557    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00558 
00559    f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */
00560    f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ;
00561    f2 = xcx[1].r ; f4 = xcx[1].i ;
00562    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00563    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00564 
00565    f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */
00566    f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ;
00567    f2 = xcx[5].r ; f4 = xcx[5].i ;
00568    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00569    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
00570 
00571    f1 = xcx[4].r ; f3 = xcx[4].i ;  /* cos=1 sin=0 */
00572    f2 = xcx[0].r ; f4 = xcx[0].i ;
00573    xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
00574    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00575 
00576    f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
00577    f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
00578    f2 = xcx[1].r ; f4 = xcx[1].i ;
00579    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
00580    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00581 
00582    f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */
00583    f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ;
00584    f2 = xcx[2].r ; f4 = xcx[2].i ;
00585    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
00586    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00587 
00588    f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
00589    f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
00590    f2 = xcx[3].r ; f4 = xcx[3].i ;
00591    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00592    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
00593 
00594    return ;
00595 }
00596 
00597 /**************************************/
00598 /* FFT routine unrolled of length  16 */
00599 
00600 static void fft16( int mode , complex * xc )
00601 {
00602    register complex * csp , * xcx=xc;
00603    register float f1,f2,f3,f4 ;
00604 
00605    /** perhaps initialize **/
00606 
00607    if( nold != 16 ) csfft_trigconsts( 16 ) ;
00608 
00609    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00610 
00611    /** data swapping part **/
00612 
00613    f1 = xcx[1].r ; f2 = xcx[1].i ;
00614    xcx[1].r = xcx[8].r ; xcx[1].i = xcx[8].i ;
00615    xcx[8].r = f1 ; xcx[8].i = f2 ;
00616 
00617    f1 = xcx[2].r ; f2 = xcx[2].i ;
00618    xcx[2].r = xcx[4].r ; xcx[2].i = xcx[4].i ;
00619    xcx[4].r = f1 ; xcx[4].i = f2 ;
00620 
00621    f1 = xcx[3].r ; f2 = xcx[3].i ;
00622    xcx[3].r = xcx[12].r ; xcx[3].i = xcx[12].i ;
00623    xcx[12].r = f1 ; xcx[12].i = f2 ;
00624 
00625    f1 = xcx[5].r ; f2 = xcx[5].i ;
00626    xcx[5].r = xcx[10].r ; xcx[5].i = xcx[10].i ;
00627    xcx[10].r = f1 ; xcx[10].i = f2 ;
00628 
00629    f1 = xcx[7].r ; f2 = xcx[7].i ;
00630    xcx[7].r = xcx[14].r ; xcx[7].i = xcx[14].i ;
00631    xcx[14].r = f1 ; xcx[14].i = f2 ;
00632 
00633    f1 = xcx[11].r ; f2 = xcx[11].i ;
00634    xcx[11].r = xcx[13].r ; xcx[11].i = xcx[13].i ;
00635    xcx[13].r = f1 ; xcx[13].i = f2 ;
00636 
00637    /** butterflying part **/
00638 
00639    f1 = xcx[1].r ; f3 = xcx[1].i ;  /* cos=1 sin=0 */
00640    f2 = xcx[0].r ; f4 = xcx[0].i ;
00641    xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
00642    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00643 
00644    f1 = xcx[3].r ; f3 = xcx[3].i ;  /* cos=1 sin=0 */
00645    f2 = xcx[2].r ; f4 = xcx[2].i ;
00646    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00647    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00648 
00649    f1 = xcx[5].r ; f3 = xcx[5].i ;  /* cos=1 sin=0 */
00650    f2 = xcx[4].r ; f4 = xcx[4].i ;
00651    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
00652    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00653 
00654    f1 = xcx[7].r ; f3 = xcx[7].i ;  /* cos=1 sin=0 */
00655    f2 = xcx[6].r ; f4 = xcx[6].i ;
00656    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00657    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
00658 
00659    f1 = xcx[9].r ; f3 = xcx[9].i ;  /* cos=1 sin=0 */
00660    f2 = xcx[8].r ; f4 = xcx[8].i ;
00661    xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
00662    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
00663 
00664    f1 = xcx[11].r ; f3 = xcx[11].i ;  /* cos=1 sin=0 */
00665    f2 = xcx[10].r ; f4 = xcx[10].i ;
00666    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
00667    xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
00668 
00669    f1 = xcx[13].r ; f3 = xcx[13].i ;  /* cos=1 sin=0 */
00670    f2 = xcx[12].r ; f4 = xcx[12].i ;
00671    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
00672    xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
00673 
00674    f1 = xcx[15].r ; f3 = xcx[15].i ;  /* cos=1 sin=0 */
00675    f2 = xcx[14].r ; f4 = xcx[14].i ;
00676    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
00677    xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
00678 
00679    f1 = xcx[2].r ; f3 = xcx[2].i ;  /* cos=1 sin=0 */
00680    f2 = xcx[0].r ; f4 = xcx[0].i ;
00681    xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
00682    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00683 
00684    f1 = xcx[6].r ; f3 = xcx[6].i ;  /* cos=1 sin=0 */
00685    f2 = xcx[4].r ; f4 = xcx[4].i ;
00686    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
00687    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00688 
00689    f1 = xcx[10].r ; f3 = xcx[10].i ;  /* cos=1 sin=0 */
00690    f2 = xcx[8].r ; f4 = xcx[8].i ;
00691    xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
00692    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
00693 
00694    f1 = xcx[14].r ; f3 = xcx[14].i ;  /* cos=1 sin=0 */
00695    f2 = xcx[12].r ; f4 = xcx[12].i ;
00696    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
00697    xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
00698 
00699    f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */
00700    f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ;
00701    f2 = xcx[1].r ; f4 = xcx[1].i ;
00702    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00703    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00704 
00705    f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */
00706    f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ;
00707    f2 = xcx[5].r ; f4 = xcx[5].i ;
00708    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00709    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
00710 
00711    f1 = xcx[11].r * csp[2].r - xcx[11].i * csp[2].i ; /* twiddles */
00712    f3 = xcx[11].r * csp[2].i + xcx[11].i * csp[2].r ;
00713    f2 = xcx[9].r ; f4 = xcx[9].i ;
00714    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
00715    xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
00716 
00717    f1 = xcx[15].r * csp[2].r - xcx[15].i * csp[2].i ; /* twiddles */
00718    f3 = xcx[15].r * csp[2].i + xcx[15].i * csp[2].r ;
00719    f2 = xcx[13].r ; f4 = xcx[13].i ;
00720    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
00721    xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
00722 
00723    f1 = xcx[4].r ; f3 = xcx[4].i ;  /* cos=1 sin=0 */
00724    f2 = xcx[0].r ; f4 = xcx[0].i ;
00725    xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
00726    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00727 
00728    f1 = xcx[12].r ; f3 = xcx[12].i ;  /* cos=1 sin=0 */
00729    f2 = xcx[8].r ; f4 = xcx[8].i ;
00730    xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
00731    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
00732 
00733    f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
00734    f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
00735    f2 = xcx[1].r ; f4 = xcx[1].i ;
00736    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
00737    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00738 
00739    f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */
00740    f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ;
00741    f2 = xcx[9].r ; f4 = xcx[9].i ;
00742    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
00743    xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
00744 
00745    f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */
00746    f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ;
00747    f2 = xcx[2].r ; f4 = xcx[2].i ;
00748    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
00749    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00750 
00751    f1 = xcx[14].r * csp[5].r - xcx[14].i * csp[5].i ; /* twiddles */
00752    f3 = xcx[14].r * csp[5].i + xcx[14].i * csp[5].r ;
00753    f2 = xcx[10].r ; f4 = xcx[10].i ;
00754    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
00755    xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
00756 
00757    f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
00758    f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
00759    f2 = xcx[3].r ; f4 = xcx[3].i ;
00760    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00761    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
00762 
00763    f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */
00764    f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ;
00765    f2 = xcx[11].r ; f4 = xcx[11].i ;
00766    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
00767    xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
00768 
00769    f1 = xcx[8].r ; f3 = xcx[8].i ;  /* cos=1 sin=0 */
00770    f2 = xcx[0].r ; f4 = xcx[0].i ;
00771    xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ;
00772    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00773 
00774    f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */
00775    f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ;
00776    f2 = xcx[1].r ; f4 = xcx[1].i ;
00777    xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
00778    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
00779 
00780    f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */
00781    f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ;
00782    f2 = xcx[2].r ; f4 = xcx[2].i ;
00783    xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
00784    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00785 
00786    f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */
00787    f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ;
00788    f2 = xcx[3].r ; f4 = xcx[3].i ;
00789    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
00790    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
00791 
00792    f1 = xcx[12].r * csp[11].r - xcx[12].i * csp[11].i ; /* twiddles */
00793    f3 = xcx[12].r * csp[11].i + xcx[12].i * csp[11].r ;
00794    f2 = xcx[4].r ; f4 = xcx[4].i ;
00795    xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
00796    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00797 
00798    f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */
00799    f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ;
00800    f2 = xcx[5].r ; f4 = xcx[5].i ;
00801    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
00802    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
00803 
00804    f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */
00805    f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ;
00806    f2 = xcx[6].r ; f4 = xcx[6].i ;
00807    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
00808    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
00809 
00810    f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */
00811    f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ;
00812    f2 = xcx[7].r ; f4 = xcx[7].i ;
00813    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
00814    xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
00815 
00816    return ;
00817 }
00818 
00819 /**************************************/
00820 /* FFT routine unrolled of length  32 */
00821 
00822 static void fft32( int mode , complex * xc )
00823 {
00824    register complex * csp , * xcx=xc;
00825    register float f1,f2,f3,f4 ;
00826 
00827    /** perhaps initialize **/
00828 
00829    if( nold != 32 ) csfft_trigconsts( 32 ) ;
00830 
00831    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00832 
00833    /** data swapping part **/
00834 
00835    f1 = xcx[1].r ; f2 = xcx[1].i ;
00836    xcx[1].r = xcx[16].r ; xcx[1].i = xcx[16].i ;
00837    xcx[16].r = f1 ; xcx[16].i = f2 ;
00838 
00839    f1 = xcx[2].r ; f2 = xcx[2].i ;
00840    xcx[2].r = xcx[8].r ; xcx[2].i = xcx[8].i ;
00841    xcx[8].r = f1 ; xcx[8].i = f2 ;
00842 
00843    f1 = xcx[3].r ; f2 = xcx[3].i ;
00844    xcx[3].r = xcx[24].r ; xcx[3].i = xcx[24].i ;
00845    xcx[24].r = f1 ; xcx[24].i = f2 ;
00846 
00847    f1 = xcx[5].r ; f2 = xcx[5].i ;
00848    xcx[5].r = xcx[20].r ; xcx[5].i = xcx[20].i ;
00849    xcx[20].r = f1 ; xcx[20].i = f2 ;
00850 
00851    f1 = xcx[6].r ; f2 = xcx[6].i ;
00852    xcx[6].r = xcx[12].r ; xcx[6].i = xcx[12].i ;
00853    xcx[12].r = f1 ; xcx[12].i = f2 ;
00854 
00855    f1 = xcx[7].r ; f2 = xcx[7].i ;
00856    xcx[7].r = xcx[28].r ; xcx[7].i = xcx[28].i ;
00857    xcx[28].r = f1 ; xcx[28].i = f2 ;
00858 
00859    f1 = xcx[9].r ; f2 = xcx[9].i ;
00860    xcx[9].r = xcx[18].r ; xcx[9].i = xcx[18].i ;
00861    xcx[18].r = f1 ; xcx[18].i = f2 ;
00862 
00863    f1 = xcx[11].r ; f2 = xcx[11].i ;
00864    xcx[11].r = xcx[26].r ; xcx[11].i = xcx[26].i ;
00865    xcx[26].r = f1 ; xcx[26].i = f2 ;
00866 
00867    f1 = xcx[13].r ; f2 = xcx[13].i ;
00868    xcx[13].r = xcx[22].r ; xcx[13].i = xcx[22].i ;
00869    xcx[22].r = f1 ; xcx[22].i = f2 ;
00870 
00871    f1 = xcx[15].r ; f2 = xcx[15].i ;
00872    xcx[15].r = xcx[30].r ; xcx[15].i = xcx[30].i ;
00873    xcx[30].r = f1 ; xcx[30].i = f2 ;
00874 
00875    f1 = xcx[19].r ; f2 = xcx[19].i ;
00876    xcx[19].r = xcx[25].r ; xcx[19].i = xcx[25].i ;
00877    xcx[25].r = f1 ; xcx[25].i = f2 ;
00878 
00879    f1 = xcx[23].r ; f2 = xcx[23].i ;
00880    xcx[23].r = xcx[29].r ; xcx[23].i = xcx[29].i ;
00881    xcx[29].r = f1 ; xcx[29].i = f2 ;
00882 
00883    /** butterflying part **/
00884 
00885    f1 = xcx[1].r ; f3 = xcx[1].i ;  /* cos=1 sin=0 */
00886    f2 = xcx[0].r ; f4 = xcx[0].i ;
00887    xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
00888    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00889 
00890    f1 = xcx[3].r ; f3 = xcx[3].i ;  /* cos=1 sin=0 */
00891    f2 = xcx[2].r ; f4 = xcx[2].i ;
00892    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
00893    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
00894 
00895    f1 = xcx[5].r ; f3 = xcx[5].i ;  /* cos=1 sin=0 */
00896    f2 = xcx[4].r ; f4 = xcx[4].i ;
00897    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
00898    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00899 
00900    f1 = xcx[7].r ; f3 = xcx[7].i ;  /* cos=1 sin=0 */
00901    f2 = xcx[6].r ; f4 = xcx[6].i ;
00902    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
00903    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
00904 
00905    f1 = xcx[9].r ; f3 = xcx[9].i ;  /* cos=1 sin=0 */
00906    f2 = xcx[8].r ; f4 = xcx[8].i ;
00907    xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
00908    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
00909 
00910    f1 = xcx[11].r ; f3 = xcx[11].i ;  /* cos=1 sin=0 */
00911    f2 = xcx[10].r ; f4 = xcx[10].i ;
00912    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
00913    xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
00914 
00915    f1 = xcx[13].r ; f3 = xcx[13].i ;  /* cos=1 sin=0 */
00916    f2 = xcx[12].r ; f4 = xcx[12].i ;
00917    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
00918    xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
00919 
00920    f1 = xcx[15].r ; f3 = xcx[15].i ;  /* cos=1 sin=0 */
00921    f2 = xcx[14].r ; f4 = xcx[14].i ;
00922    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
00923    xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
00924 
00925    f1 = xcx[17].r ; f3 = xcx[17].i ;  /* cos=1 sin=0 */
00926    f2 = xcx[16].r ; f4 = xcx[16].i ;
00927    xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
00928    xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
00929 
00930    f1 = xcx[19].r ; f3 = xcx[19].i ;  /* cos=1 sin=0 */
00931    f2 = xcx[18].r ; f4 = xcx[18].i ;
00932    xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
00933    xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
00934 
00935    f1 = xcx[21].r ; f3 = xcx[21].i ;  /* cos=1 sin=0 */
00936    f2 = xcx[20].r ; f4 = xcx[20].i ;
00937    xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
00938    xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
00939 
00940    f1 = xcx[23].r ; f3 = xcx[23].i ;  /* cos=1 sin=0 */
00941    f2 = xcx[22].r ; f4 = xcx[22].i ;
00942    xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
00943    xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
00944 
00945    f1 = xcx[25].r ; f3 = xcx[25].i ;  /* cos=1 sin=0 */
00946    f2 = xcx[24].r ; f4 = xcx[24].i ;
00947    xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
00948    xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
00949 
00950    f1 = xcx[27].r ; f3 = xcx[27].i ;  /* cos=1 sin=0 */
00951    f2 = xcx[26].r ; f4 = xcx[26].i ;
00952    xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
00953    xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
00954 
00955    f1 = xcx[29].r ; f3 = xcx[29].i ;  /* cos=1 sin=0 */
00956    f2 = xcx[28].r ; f4 = xcx[28].i ;
00957    xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
00958    xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
00959 
00960    f1 = xcx[31].r ; f3 = xcx[31].i ;  /* cos=1 sin=0 */
00961    f2 = xcx[30].r ; f4 = xcx[30].i ;
00962    xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
00963    xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ;
00964 
00965    f1 = xcx[2].r ; f3 = xcx[2].i ;  /* cos=1 sin=0 */
00966    f2 = xcx[0].r ; f4 = xcx[0].i ;
00967    xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
00968    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
00969 
00970    f1 = xcx[6].r ; f3 = xcx[6].i ;  /* cos=1 sin=0 */
00971    f2 = xcx[4].r ; f4 = xcx[4].i ;
00972    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
00973    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
00974 
00975    f1 = xcx[10].r ; f3 = xcx[10].i ;  /* cos=1 sin=0 */
00976    f2 = xcx[8].r ; f4 = xcx[8].i ;
00977    xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
00978    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
00979 
00980    f1 = xcx[14].r ; f3 = xcx[14].i ;  /* cos=1 sin=0 */
00981    f2 = xcx[12].r ; f4 = xcx[12].i ;
00982    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
00983    xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
00984 
00985    f1 = xcx[18].r ; f3 = xcx[18].i ;  /* cos=1 sin=0 */
00986    f2 = xcx[16].r ; f4 = xcx[16].i ;
00987    xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
00988    xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
00989 
00990    f1 = xcx[22].r ; f3 = xcx[22].i ;  /* cos=1 sin=0 */
00991    f2 = xcx[20].r ; f4 = xcx[20].i ;
00992    xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
00993    xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
00994 
00995    f1 = xcx[26].r ; f3 = xcx[26].i ;  /* cos=1 sin=0 */
00996    f2 = xcx[24].r ; f4 = xcx[24].i ;
00997    xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
00998    xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
00999 
01000    f1 = xcx[30].r ; f3 = xcx[30].i ;  /* cos=1 sin=0 */
01001    f2 = xcx[28].r ; f4 = xcx[28].i ;
01002    xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
01003    xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
01004 
01005    f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */
01006    f3 = xcx[3].r * csp[2].i ;
01007    f2 = xcx[1].r ; f4 = xcx[1].i ;
01008    xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
01009    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
01010 
01011    f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */
01012    f3 = xcx[7].r * csp[2].i ;
01013    f2 = xcx[5].r ; f4 = xcx[5].i ;
01014    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
01015    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
01016 
01017    f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */
01018    f3 = xcx[11].r * csp[2].i ;
01019    f2 = xcx[9].r ; f4 = xcx[9].i ;
01020    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
01021    xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
01022 
01023    f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */
01024    f3 = xcx[15].r * csp[2].i ;
01025    f2 = xcx[13].r ; f4 = xcx[13].i ;
01026    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
01027    xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
01028 
01029    f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */
01030    f3 = xcx[19].r * csp[2].i ;
01031    f2 = xcx[17].r ; f4 = xcx[17].i ;
01032    xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
01033    xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
01034 
01035    f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */
01036    f3 = xcx[23].r * csp[2].i ;
01037    f2 = xcx[21].r ; f4 = xcx[21].i ;
01038    xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
01039    xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
01040 
01041    f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */
01042    f3 = xcx[27].r * csp[2].i ;
01043    f2 = xcx[25].r ; f4 = xcx[25].i ;
01044    xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
01045    xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
01046 
01047    f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */
01048    f3 = xcx[31].r * csp[2].i ;
01049    f2 = xcx[29].r ; f4 = xcx[29].i ;
01050    xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
01051    xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ;
01052 
01053    f1 = xcx[4].r ; f3 = xcx[4].i ;  /* cos=1 sin=0 */
01054    f2 = xcx[0].r ; f4 = xcx[0].i ;
01055    xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
01056    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
01057 
01058    f1 = xcx[12].r ; f3 = xcx[12].i ;  /* cos=1 sin=0 */
01059    f2 = xcx[8].r ; f4 = xcx[8].i ;
01060    xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
01061    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
01062 
01063    f1 = xcx[20].r ; f3 = xcx[20].i ;  /* cos=1 sin=0 */
01064    f2 = xcx[16].r ; f4 = xcx[16].i ;
01065    xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
01066    xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
01067 
01068    f1 = xcx[28].r ; f3 = xcx[28].i ;  /* cos=1 sin=0 */
01069    f2 = xcx[24].r ; f4 = xcx[24].i ;
01070    xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
01071    xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
01072 
01073    f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
01074    f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
01075    f2 = xcx[1].r ; f4 = xcx[1].i ;
01076    xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
01077    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
01078 
01079    f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */
01080    f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ;
01081    f2 = xcx[9].r ; f4 = xcx[9].i ;
01082    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
01083    xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
01084 
01085    f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */
01086    f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ;
01087    f2 = xcx[17].r ; f4 = xcx[17].i ;
01088    xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
01089    xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
01090 
01091    f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */
01092    f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ;
01093    f2 = xcx[25].r ; f4 = xcx[25].i ;
01094    xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
01095    xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
01096 
01097    f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */
01098    f3 = xcx[6].r * csp[5].i ;
01099    f2 = xcx[2].r ; f4 = xcx[2].i ;
01100    xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
01101    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
01102 
01103    f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */
01104    f3 = xcx[14].r * csp[5].i ;
01105    f2 = xcx[10].r ; f4 = xcx[10].i ;
01106    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
01107    xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
01108 
01109    f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */
01110    f3 = xcx[22].r * csp[5].i ;
01111    f2 = xcx[18].r ; f4 = xcx[18].i ;
01112    xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
01113    xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
01114 
01115    f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */
01116    f3 = xcx[30].r * csp[5].i ;
01117    f2 = xcx[26].r ; f4 = xcx[26].i ;
01118    xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
01119    xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
01120 
01121    f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
01122    f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
01123    f2 = xcx[3].r ; f4 = xcx[3].i ;
01124    xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
01125    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
01126 
01127    f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */
01128    f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ;
01129    f2 = xcx[11].r ; f4 = xcx[11].i ;
01130    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
01131    xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
01132 
01133    f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */
01134    f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ;
01135    f2 = xcx[19].r ; f4 = xcx[19].i ;
01136    xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
01137    xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
01138 
01139    f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */
01140    f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ;
01141    f2 = xcx[27].r ; f4 = xcx[27].i ;
01142    xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
01143    xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ;
01144 
01145    f1 = xcx[8].r ; f3 = xcx[8].i ;  /* cos=1 sin=0 */
01146    f2 = xcx[0].r ; f4 = xcx[0].i ;
01147    xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ;
01148    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
01149 
01150    f1 = xcx[24].r ; f3 = xcx[24].i ;  /* cos=1 sin=0 */
01151    f2 = xcx[16].r ; f4 = xcx[16].i ;
01152    xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
01153    xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
01154 
01155    f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */
01156    f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ;
01157    f2 = xcx[1].r ; f4 = xcx[1].i ;
01158    xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
01159    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
01160 
01161    f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */
01162    f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ;
01163    f2 = xcx[17].r ; f4 = xcx[17].i ;
01164    xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
01165    xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
01166 
01167    f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */
01168    f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ;
01169    f2 = xcx[2].r ; f4 = xcx[2].i ;
01170    xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
01171    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
01172 
01173    f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */
01174    f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ;
01175    f2 = xcx[18].r ; f4 = xcx[18].i ;
01176    xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
01177    xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
01178 
01179    f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */
01180    f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ;
01181    f2 = xcx[3].r ; f4 = xcx[3].i ;
01182    xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
01183    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
01184 
01185    f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */
01186    f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ;
01187    f2 = xcx[19].r ; f4 = xcx[19].i ;
01188    xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
01189    xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
01190 
01191    f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */
01192    f3 = xcx[12].r * csp[11].i ;
01193    f2 = xcx[4].r ; f4 = xcx[4].i ;
01194    xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
01195    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
01196 
01197    f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */
01198    f3 = xcx[28].r * csp[11].i ;
01199    f2 = xcx[20].r ; f4 = xcx[20].i ;
01200    xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
01201    xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
01202 
01203    f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */
01204    f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ;
01205    f2 = xcx[5].r ; f4 = xcx[5].i ;
01206    xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
01207    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
01208 
01209    f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */
01210    f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ;
01211    f2 = xcx[21].r ; f4 = xcx[21].i ;
01212    xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
01213    xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
01214 
01215    f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */
01216    f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ;
01217    f2 = xcx[6].r ; f4 = xcx[6].i ;
01218    xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
01219    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
01220 
01221    f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */
01222    f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ;
01223    f2 = xcx[22].r ; f4 = xcx[22].i ;
01224    xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
01225    xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
01226 
01227    f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */
01228    f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ;
01229    f2 = xcx[7].r ; f4 = xcx[7].i ;
01230    xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
01231    xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
01232 
01233    f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */
01234    f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ;
01235    f2 = xcx[23].r ; f4 = xcx[23].i ;
01236    xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
01237    xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ;
01238 
01239    f1 = xcx[16].r ; f3 = xcx[16].i ;  /* cos=1 sin=0 */
01240    f2 = xcx[0].r ; f4 = xcx[0].i ;
01241    xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ;
01242    xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
01243 
01244    f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */
01245    f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ;
01246    f2 = xcx[1].r ; f4 = xcx[1].i ;
01247    xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
01248    xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
01249 
01250    f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */
01251    f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ;
01252    f2 = xcx[2].r ; f4 = xcx[2].i ;
01253    xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
01254    xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
01255 
01256    f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */
01257    f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ;
01258    f2 = xcx[3].r ; f4 = xcx[3].i ;
01259    xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
01260    xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
01261 
01262    f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */
01263    f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ;
01264    f2 = xcx[4].r ; f4 = xcx[4].i ;
01265    xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
01266    xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
01267 
01268    f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */
01269    f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ;
01270    f2 = xcx[5].r ; f4 = xcx[5].i ;
01271    xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
01272    xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
01273 
01274    f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */
01275    f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ;
01276    f2 = xcx[6].r ; f4 = xcx[6].i ;
01277    xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
01278    xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
01279 
01280    f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */
01281    f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ;
01282    f2 = xcx[7].r ; f4 = xcx[7].i ;
01283    xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
01284    xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
01285 
01286    f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */
01287    f3 = xcx[24].r * csp[23].i ;
01288    f2 = xcx[8].r ; f4 = xcx[8].i ;
01289    xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
01290    xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
01291 
01292    f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */
01293    f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ;
01294    f2 = xcx[9].r ; f4 = xcx[9].i ;
01295    xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
01296    xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
01297 
01298    f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */
01299    f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ;
01300    f2 = xcx[10].r ; f4 = xcx[10].i ;
01301    xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
01302    xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
01303 
01304    f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */
01305    f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ;
01306    f2 = xcx[11].r ; f4 = xcx[11].i ;
01307    xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
01308    xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
01309 
01310    f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */
01311    f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ;
01312    f2 = xcx[12].r ; f4 = xcx[12].i ;
01313    xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
01314    xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
01315 
01316    f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */
01317    f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ;
01318    f2 = xcx[13].r ; f4 = xcx[13].i ;
01319    xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
01320    xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
01321 
01322    f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */
01323    f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ;
01324    f2 = xcx[14].r ; f4 = xcx[14].i ;
01325    xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
01326    xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
01327 
01328    f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */
01329    f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ;
01330    f2 = xcx[15].r ; f4 = xcx[15].i ;
01331    xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
01332    xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ;
01333 
01334    return ;
01335 }
01336 
01337 /*----------------------------------------------------------------
01338    Do a 64 FFT using fft32 and decimation-by-2
01339 ------------------------------------------------------------------*/
01340 
01341 static void fft64( int mode , complex * xc )
01342 {
01343    static complex * cs=NULL , * aa , * bb ;
01344    register int k , i ;
01345    register float akr,aki , bkr,bki , tr,ti , t1,t2 ;
01346 
01347    /*-- initialize cosine/sine table --*/
01348 
01349    if( cs == NULL ){
01350       double th = (PI/32.0) ;
01351       cs = (complex *) malloc(sizeof(complex)*32) ;
01352       aa = (complex *) malloc(sizeof(complex)*32) ;
01353       bb = (complex *) malloc(sizeof(complex)*32) ;
01354 
01355       cs[0].r = 1.0 ; cs[0].i = 0.0 ; /* never used */
01356       for( k=1 ; k < 32 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
01357    }
01358 
01359    /*-- load subarrays, and FFT each one --*/
01360 
01361    for( i=k=0 ; k < 32 ; k++ ){ aa[k] = xc[i++] ; bb[k] = xc[i++] ; }
01362 
01363    fft32( mode , aa ) ; fft32( mode , bb ) ;
01364 
01365    /*-- recombine: 0 and Nyquist --*/
01366 
01367    xc[ 0].r = aa[0].r + bb[0].r ; xc[ 0].i = aa[0].i + bb[0].i ;
01368    xc[32].r = aa[0].r - bb[0].r ; xc[32].i = aa[0].i - bb[0].i ;
01369 
01370    /*-- recombine: all others --*/
01371 
01372    if( mode > 0 ){
01373       for( k=1 ; k < 32 ; k++ ){
01374          bkr = bb[k].r; bki = bb[k].i;
01375          tr  = cs[k].r; ti  = cs[k].i;
01376 
01377          t1  = tr*bkr - ti*bki ; t2 = tr*bki + ti*bkr ;
01378 
01379          akr = aa[k].r; aki = aa[k].i;
01380 
01381          xc[k   ].r = akr+t1; xc[k   ].i = aki+t2;
01382          xc[k+32].r = akr-t1; xc[k+32].i = aki-t2;
01383       }
01384    } else {
01385       for( k=1 ; k < 32 ; k++ ){
01386          bkr = bb[k].r; bki = bb[k].i;
01387          tr  = cs[k].r; ti  = cs[k].i;
01388                                                          /* only change from */
01389          t1  = tr*bkr + ti*bki ; t2 = tr*bki - ti*bkr ;  /* above is with ti */
01390 
01391          akr = aa[k].r; aki = aa[k].i;
01392 
01393          xc[k   ].r = akr+t1; xc[k   ].i = aki+t2;
01394          xc[k+32].r = akr-t1; xc[k+32].i = aki-t2;
01395       }
01396    }
01397 
01398    return ;
01399 }
01400 
01401 /*================================================================
01402   Decimation by 4 routines below from 12 Aug 1999 -- RWCox
01403   On a 400 MHz PII, about 8% faster than 2 decimation by 2 steps.
01404   Added fft512, fft1024, and fft2048 also, which are 50-80%
01405   faster than the unrolled csfft_cox() routine on the PII --
01406   somewhat less speedup on SGI R10K, but still worth something.
01407 ==================================================================*/
01408 
01409 /*----------------------------------------------------------------
01410    Do a 128 FFT using fft32 and decimation by 4
01411 ------------------------------------------------------------------*/
01412 
01413 #define N  128
01414 #define M   32
01415 #define M2  64
01416 #define M3  96
01417 static void fft128( int mode , complex * xc )
01418 {
01419    static complex * cs=NULL , * aa , * bb , * cc , * dd ;
01420    register int k , i ;
01421    register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 ,
01422                   acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ;
01423 
01424    /*-- initialize cosine/sine table and memory space --*/
01425 
01426    if( cs == NULL ){
01427       double th = (2.0*PI/N) ;
01428       cs = (complex *) malloc(sizeof(complex)*M3) ;
01429       aa = (complex *) malloc(sizeof(complex)*M ) ;
01430       bb = (complex *) malloc(sizeof(complex)*M ) ;
01431       cc = (complex *) malloc(sizeof(complex)*M ) ;
01432       dd = (complex *) malloc(sizeof(complex)*M ) ;
01433 
01434       cs[0].r = 1.0 ; cs[0].i = 0.0 ;
01435       for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
01436    }
01437 
01438    /*-- load subarrays, and FFT each one --*/
01439 
01440    for( i=k=0; k < M; k++ ){
01441       aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++];
01442    }
01443 
01444    fft32(mode,aa); fft32(mode,bb); fft32(mode,cc); fft32(mode,dd);
01445 
01446    /*-- recombination --*/
01447 
01448    if( mode > 0 ){
01449       for( k=0 ; k < M ; k++ ){
01450          t1  = bb[k].r; t2  = bb[k].i;
01451          tr  = cs[k].r; ti  = cs[k].i;
01452          bbr = tr*t1 - ti*t2  ; bbi = tr*t2 + ti*t1  ;  /* b[k]*exp(+2*Pi*k/N) */
01453 
01454          t1  = cc[  k].r; t2  = cc[  k].i;
01455          tr  = cs[2*k].r; ti  = cs[2*k].i;
01456          ccr = tr*t1 - ti*t2  ; cci = tr*t2 + ti*t1  ;  /* c[k]*exp(+4*Pi*k/N) */
01457 
01458          t1  = dd[  k].r; t2  = dd[  k].i;
01459          tr  = cs[3*k].r; ti  = cs[3*k].i;
01460          ddr = tr*t1 - ti*t2  ; ddi = tr*t2 + ti*t1  ;  /* d[k]*exp(+6*Pi*k/N) */
01461 
01462          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01463 
01464          acpr = aar + ccr ; acmr = aar - ccr ;
01465          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01466          acpi = aai + cci ; acmi = aai - cci ;
01467          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01468 
01469          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01470          xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ;
01471          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01472          xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ;
01473       }
01474    } else {
01475       for( k=0 ; k < M ; k++ ){
01476          t1  = bb[k].r; t2  = bb[k].i;
01477          tr  = cs[k].r; ti  = cs[k].i;
01478          bbr = tr*t1 + ti*t2  ; bbi = tr*t2 - ti*t1  ;  /* b[k]*exp(-2*Pi*k/N) */
01479 
01480          t1  = cc[  k].r; t2  = cc[  k].i;
01481          tr  = cs[2*k].r; ti  = cs[2*k].i;
01482          ccr = tr*t1 + ti*t2  ; cci = tr*t2 - ti*t1  ;  /* c[k]*exp(-4*Pi*k/N) */
01483 
01484          t1  = dd[  k].r; t2  = dd[  k].i;
01485          tr  = cs[3*k].r; ti  = cs[3*k].i;
01486          ddr = tr*t1 + ti*t2  ; ddi = tr*t2 - ti*t1  ;  /* d[k]*exp(-6*Pi*k/N) */
01487 
01488          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01489 
01490          acpr = aar + ccr ; acmr = aar - ccr ;
01491          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01492          acpi = aai + cci ; acmi = aai - cci ;
01493          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01494 
01495          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01496          xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ;
01497          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01498          xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ;
01499       }
01500    }
01501 
01502    return ;
01503 }
01504 #undef N
01505 #undef M
01506 #undef M2
01507 #undef M3
01508 
01509 /*------------------------------------------------------------------
01510    Do a 256 FFT using fft64 and decimation by 4
01511 --------------------------------------------------------------------*/
01512 
01513 #define N  256
01514 #define M   64
01515 #define M2 128
01516 #define M3 192
01517 static void fft256( int mode , complex * xc )
01518 {
01519    static complex * cs=NULL , * aa , * bb , * cc , * dd ;
01520    register int k , i ;
01521    register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 ,
01522                   acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ;
01523 
01524    /*-- initialize cosine/sine table and memory space --*/
01525 
01526    if( cs == NULL ){
01527       double th = (2.0*PI/N) ;
01528       cs = (complex *) malloc(sizeof(complex)*M3) ;
01529       aa = (complex *) malloc(sizeof(complex)*M ) ;
01530       bb = (complex *) malloc(sizeof(complex)*M ) ;
01531       cc = (complex *) malloc(sizeof(complex)*M ) ;
01532       dd = (complex *) malloc(sizeof(complex)*M ) ;
01533 
01534       cs[0].r = 1.0 ; cs[0].i = 0.0 ;
01535       for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
01536    }
01537 
01538    /*-- load subarrays, and FFT each one --*/
01539 
01540    for( i=k=0; k < M; k++ ){
01541       aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++];
01542    }
01543 
01544    fft64(mode,aa); fft64(mode,bb); fft64(mode,cc); fft64(mode,dd);
01545 
01546    /*-- recombination --*/
01547 
01548    if( mode > 0 ){
01549       for( k=0 ; k < M ; k++ ){
01550          t1  = bb[k].r; t2  = bb[k].i;
01551          tr  = cs[k].r; ti  = cs[k].i;
01552          bbr = tr*t1 - ti*t2  ; bbi = tr*t2 + ti*t1  ;  /* b[k]*exp(+2*Pi*k/N) */
01553 
01554          t1  = cc[  k].r; t2  = cc[  k].i;
01555          tr  = cs[2*k].r; ti  = cs[2*k].i;
01556          ccr = tr*t1 - ti*t2  ; cci = tr*t2 + ti*t1  ;  /* c[k]*exp(+4*Pi*k/N) */
01557 
01558          t1  = dd[  k].r; t2  = dd[  k].i;
01559          tr  = cs[3*k].r; ti  = cs[3*k].i;
01560          ddr = tr*t1 - ti*t2  ; ddi = tr*t2 + ti*t1  ;  /* d[k]*exp(+6*Pi*k/N) */
01561 
01562          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01563 
01564          acpr = aar + ccr ; acmr = aar - ccr ;
01565          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01566          acpi = aai + cci ; acmi = aai - cci ;
01567          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01568 
01569          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01570          xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ;
01571          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01572          xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ;
01573       }
01574    } else {
01575       for( k=0 ; k < M ; k++ ){
01576          t1  = bb[k].r; t2  = bb[k].i;
01577          tr  = cs[k].r; ti  = cs[k].i;
01578          bbr = tr*t1 + ti*t2  ; bbi = tr*t2 - ti*t1  ;  /* b[k]*exp(-2*Pi*k/N) */
01579 
01580          t1  = cc[  k].r; t2  = cc[  k].i;
01581          tr  = cs[2*k].r; ti  = cs[2*k].i;
01582          ccr = tr*t1 + ti*t2  ; cci = tr*t2 - ti*t1  ;  /* c[k]*exp(-4*Pi*k/N) */
01583 
01584          t1  = dd[  k].r; t2  = dd[  k].i;
01585          tr  = cs[3*k].r; ti  = cs[3*k].i;
01586          ddr = tr*t1 + ti*t2  ; ddi = tr*t2 - ti*t1  ;  /* d[k]*exp(-6*Pi*k/N) */
01587 
01588          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01589 
01590          acpr = aar + ccr ; acmr = aar - ccr ;
01591          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01592          acpi = aai + cci ; acmi = aai - cci ;
01593          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01594 
01595          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01596          xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ;
01597          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01598          xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ;
01599       }
01600    }
01601 
01602    return ;
01603 }
01604 #undef N
01605 #undef M
01606 #undef M2
01607 #undef M3
01608 
01609 /*------------------------------------------------------------------
01610    Do a 512 FFT using fft128 and decimation by 4
01611 --------------------------------------------------------------------*/
01612 
01613 #define N  512
01614 #define M  128
01615 #define M2 256
01616 #define M3 384
01617 static void fft512( int mode , complex * xc )
01618 {
01619    static complex * cs=NULL , * aa , * bb , * cc , * dd ;
01620    register int k , i ;
01621    register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 ,
01622                   acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ;
01623 
01624    /*-- initialize cosine/sine table and memory space --*/
01625 
01626    if( cs == NULL ){
01627       double th = (2.0*PI/N) ;
01628       cs = (complex *) malloc(sizeof(complex)*M3) ;
01629       aa = (complex *) malloc(sizeof(complex)*M ) ;
01630       bb = (complex *) malloc(sizeof(complex)*M ) ;
01631       cc = (complex *) malloc(sizeof(complex)*M ) ;
01632       dd = (complex *) malloc(sizeof(complex)*M ) ;
01633 
01634       cs[0].r = 1.0 ; cs[0].i = 0.0 ;
01635       for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
01636    }
01637 
01638    /*-- load subarrays, and FFT each one --*/
01639 
01640    for( i=k=0; k < M; k++ ){
01641       aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++];
01642    }
01643 
01644    fft128(mode,aa); fft128(mode,bb); fft128(mode,cc); fft128(mode,dd);
01645 
01646    /*-- recombination --*/
01647 
01648    if( mode > 0 ){
01649       for( k=0 ; k < M ; k++ ){
01650          t1  = bb[k].r; t2  = bb[k].i;
01651          tr  = cs[k].r; ti  = cs[k].i;
01652          bbr = tr*t1 - ti*t2  ; bbi = tr*t2 + ti*t1  ;  /* b[k]*exp(+2*Pi*k/N) */
01653 
01654          t1  = cc[  k].r; t2  = cc[  k].i;
01655          tr  = cs[2*k].r; ti  = cs[2*k].i;
01656          ccr = tr*t1 - ti*t2  ; cci = tr*t2 + ti*t1  ;  /* c[k]*exp(+4*Pi*k/N) */
01657 
01658          t1  = dd[  k].r; t2  = dd[  k].i;
01659          tr  = cs[3*k].r; ti  = cs[3*k].i;
01660          ddr = tr*t1 - ti*t2  ; ddi = tr*t2 + ti*t1  ;  /* d[k]*exp(+6*Pi*k/N) */
01661 
01662          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01663 
01664          acpr = aar + ccr ; acmr = aar - ccr ;
01665          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01666          acpi = aai + cci ; acmi = aai - cci ;
01667          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01668 
01669          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01670          xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ;
01671          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01672          xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ;
01673       }
01674    } else {
01675       for( k=0 ; k < M ; k++ ){
01676          t1  = bb[k].r; t2  = bb[k].i;
01677          tr  = cs[k].r; ti  = cs[k].i;
01678          bbr = tr*t1 + ti*t2  ; bbi = tr*t2 - ti*t1  ;  /* b[k]*exp(-2*Pi*k/N) */
01679 
01680          t1  = cc[  k].r; t2  = cc[  k].i;
01681          tr  = cs[2*k].r; ti  = cs[2*k].i;
01682          ccr = tr*t1 + ti*t2  ; cci = tr*t2 - ti*t1  ;  /* c[k]*exp(-4*Pi*k/N) */
01683 
01684          t1  = dd[  k].r; t2  = dd[  k].i;
01685          tr  = cs[3*k].r; ti  = cs[3*k].i;
01686          ddr = tr*t1 + ti*t2  ; ddi = tr*t2 - ti*t1  ;  /* d[k]*exp(-6*Pi*k/N) */
01687 
01688          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01689 
01690          acpr = aar + ccr ; acmr = aar - ccr ;
01691          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01692          acpi = aai + cci ; acmi = aai - cci ;
01693          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01694 
01695          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01696          xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ;
01697          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01698          xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ;
01699       }
01700    }
01701 
01702    return ;
01703 }
01704 #undef N
01705 #undef M
01706 #undef M2
01707 #undef M3
01708 
01709 /*------------------------------------------------------------------
01710    fft_4dec: do a decimation by 4 for N=1024 and 2048.
01711    At most RMAX levels of recursion are allowed.
01712 --------------------------------------------------------------------*/
01713 
01714 static void fft_4dec( int mode , int idim , complex * xc )
01715 {
01716    static int rec=0 ;
01717    static int rmold[RMAX] = {-1,-1,-1} ;
01718    static complex *rcs[RMAX] = {NULL,NULL,NULL} ;
01719    static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX] , *rdd[RMAX] ;
01720 
01721    int N=idim , M=idim/4 , M2=2*M , M3=3*M ;
01722    int mold=-1 ;
01723    complex * cs , * aa , * bb , * cc , * dd ;
01724    register int k , i ;
01725    register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 ,
01726                   acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ;
01727 
01728    /*-- initialize cosine/sine table and memory space --*/
01729 
01730    if( rec >= RMAX ){
01731       fprintf(stderr,"\n*** fft_4dec: too many recursions!\n"); EXIT(1);
01732    }
01733 
01734    mold = rmold[rec] ;  /* what was done before at this recursion level */
01735 
01736    if( M != mold ){
01737       double th = (2.0*PI/N) ;
01738       if( M > mold ){
01739          if( rcs[rec] != NULL ){
01740             free(rcs[rec]);free(raa[rec]);
01741             free(rbb[rec]);free(rcc[rec]);free(rdd[rec]);
01742          }
01743          rcs[rec] = (complex *) malloc(sizeof(complex)*M3) ;
01744          raa[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01745          rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01746          rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01747          rdd[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01748       }
01749 
01750       rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ;
01751       for( k=1 ; k < M3 ; k++ ){
01752          rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th);
01753       }
01754       rmold[rec] = M ;
01755    }
01756 
01757    cs = rcs[rec]; aa = raa[rec]; bb = rbb[rec]; cc = rcc[rec]; dd = rdd[rec];
01758 
01759    /*-- load subarrays, and FFT each one --*/
01760 
01761    for( i=k=0; k < M; k++ ){
01762       aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++];
01763    }
01764 
01765    rec++ ;
01766    switch( N ){
01767 
01768       default: fprintf(stderr,"\n*** Illegal call to fft_4dec: N=%d\n",N);
01769                EXIT(1) ;
01770 
01771       case 1024: fft256(mode,aa); fft256(mode,bb);
01772                  fft256(mode,cc); fft256(mode,dd); break;
01773 
01774       case 2048: fft512(mode,aa); fft512(mode,bb);
01775                  fft512(mode,cc); fft512(mode,dd); break;
01776 
01777       case 4096: fft_4dec(mode,1024,aa); fft_4dec(mode,1024,bb);         /* recurse once */
01778                  fft_4dec(mode,1024,cc); fft_4dec(mode,1024,dd); break ;
01779 
01780       case 8192: fft_4dec(mode,2048,aa); fft_4dec(mode,2048,bb);         /* recurse once */
01781                  fft_4dec(mode,2048,cc); fft_4dec(mode,2048,dd); break ;
01782 
01783       case 16384: fft_4dec(mode,4096,aa); fft_4dec(mode,4096,bb);        /* recurse twice */
01784                   fft_4dec(mode,4096,cc); fft_4dec(mode,4096,dd); break ;
01785 
01786       case 32768: fft_4dec(mode,8192,aa); fft_4dec(mode,8192,bb);        /* recurse twice */
01787                   fft_4dec(mode,8192,cc); fft_4dec(mode,8192,dd); break ;
01788    }
01789    rec-- ;
01790 
01791    /*-- recombination --*/
01792 
01793    if( mode > 0 ){
01794       for( k=0 ; k < M ; k++ ){
01795          t1  = bb[k].r; t2  = bb[k].i;
01796          tr  = cs[k].r; ti  = cs[k].i;
01797          bbr = tr*t1 - ti*t2  ; bbi = tr*t2 + ti*t1  ;  /* b[k]*exp(+2*Pi*k/N) */
01798 
01799          t1  = cc[  k].r; t2  = cc[  k].i;
01800          tr  = cs[2*k].r; ti  = cs[2*k].i;
01801          ccr = tr*t1 - ti*t2  ; cci = tr*t2 + ti*t1  ;  /* c[k]*exp(+4*Pi*k/N) */
01802 
01803          t1  = dd[  k].r; t2  = dd[  k].i;
01804          tr  = cs[3*k].r; ti  = cs[3*k].i;
01805          ddr = tr*t1 - ti*t2  ; ddi = tr*t2 + ti*t1  ;  /* d[k]*exp(+6*Pi*k/N) */
01806 
01807          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01808 
01809          acpr = aar + ccr ; acmr = aar - ccr ;
01810          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01811          acpi = aai + cci ; acmi = aai - cci ;
01812          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01813 
01814          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01815          xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ;
01816          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01817          xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ;
01818       }
01819    } else {
01820       for( k=0 ; k < M ; k++ ){
01821          t1  = bb[k].r; t2  = bb[k].i;
01822          tr  = cs[k].r; ti  = cs[k].i;
01823          bbr = tr*t1 + ti*t2  ; bbi = tr*t2 - ti*t1  ;  /* b[k]*exp(-2*Pi*k/N) */
01824 
01825          t1  = cc[  k].r; t2  = cc[  k].i;
01826          tr  = cs[2*k].r; ti  = cs[2*k].i;
01827          ccr = tr*t1 + ti*t2  ; cci = tr*t2 - ti*t1  ;  /* c[k]*exp(-4*Pi*k/N) */
01828 
01829          t1  = dd[  k].r; t2  = dd[  k].i;
01830          tr  = cs[3*k].r; ti  = cs[3*k].i;
01831          ddr = tr*t1 + ti*t2  ; ddi = tr*t2 - ti*t1  ;  /* d[k]*exp(-6*Pi*k/N) */
01832 
01833          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01834 
01835          acpr = aar + ccr ; acmr = aar - ccr ;
01836          bdpr = bbr + ddr ; bdmr = bbr - ddr ;
01837          acpi = aai + cci ; acmi = aai - cci ;
01838          bdpi = bbi + ddi ; bdmi = bbi - ddi ;
01839 
01840          xc[k   ].r = acpr+bdpr ; xc[k   ].i = acpi+bdpi ;
01841          xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ;
01842          xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
01843          xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ;
01844       }
01845    }
01846 
01847    return ;
01848 }
01849 #endif /* DONT_UNROLL_FFTS */
01850 
01851 /*=======================================================================
01852   The following radix-3 and radix-5 routines are by RWCox -- Aug 1999.
01853   They are not inside unrolled code.
01854 =========================================================================*/
01855 
01856 /*------------------------------------------------------------------
01857    fft_3dec: do a decimation-by-3, plus a power-of-2.
01858    At most RMAX levels of recursion are allowed, which means
01859    that at most 3**RMAX can be a factor of idim.
01860 --------------------------------------------------------------------*/
01861 
01862 #undef  CC3
01863 #undef  SS3
01864 #define CC3 (-0.5)         /* cos(2*Pi/3) */
01865 #define SS3 (0.8660254038) /* sin(2*Pi/3) */
01866 
01867 static void fft_3dec( int mode , int idim , complex * xc )
01868 {
01869    static int rec=0 ;
01870    static int rmold[RMAX] = {-1,-1,-1} ;
01871    static complex *rcs[RMAX] = {NULL,NULL,NULL} ;
01872    static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX] ;
01873 
01874    int N=idim , M=idim/3 , M2=2*M ;
01875    int mold ;
01876    complex * cs=NULL , * aa , * bb , * cc ;
01877    register int k , i ;
01878    register float aar,aai, tr,ti, bbr,bbi, ccr,cci,
01879                   t1,t2,t4,t5,t6,t8 ;
01880 
01881    /*-- initialize cosine/sine table and memory space --*/
01882 
01883    if( rec >= RMAX ){
01884       fprintf(stderr,"\n*** fft_3dec: too many recursions!\n"); EXIT(1);
01885    }
01886 
01887    mold = rmold[rec] ;  /* what was done before at this recursion level */
01888 
01889    if( M != mold ){
01890       double th = (2.0*PI/N) ;
01891       if( M > mold ){
01892          if( rcs[rec] != NULL ){
01893             free(rcs[rec]);free(raa[rec]);free(rbb[rec]);free(rcc[rec]);
01894          }
01895          rcs[rec] = (complex *) malloc(sizeof(complex)*M2) ;
01896          raa[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01897          rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01898          rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ;
01899       }
01900 
01901       rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ;
01902       for( k=1 ; k < M2 ; k++ ){
01903          rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th);
01904       }
01905       rmold[rec] = M ;
01906    }
01907 
01908    cs = rcs[rec] ; aa = raa[rec] ; bb = rbb[rec] ; cc = rcc[rec] ;
01909 
01910    /*-- load subarrays, and FFT each one --*/
01911 
01912    for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; }
01913 
01914 #ifndef DONT_UNROLL_FFTS
01915    switch( M ){
01916       case   1: break ;
01917 
01918       case   2: fft2  (mode,aa) ; fft2  (mode,bb) ; fft2  (mode,cc) ; break;
01919       case   4: fft4  (mode,aa) ; fft4  (mode,bb) ; fft4  (mode,cc) ; break;
01920       case   8: fft8  (mode,aa) ; fft8  (mode,bb) ; fft8  (mode,cc) ; break;
01921       case  16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; break;
01922       case  32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; break;
01923       case  64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; break;
01924       case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; break;
01925       case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; break;
01926       case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; break;
01927 
01928       case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; break;
01929       case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; break;
01930 
01931       default:  rec++ ;
01932                 csfft_cox(mode,M,aa) ;
01933                 csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ;
01934                 rec-- ; break ;
01935    }
01936 #else
01937       rec++ ;
01938       csfft_cox(mode,M,aa) ; csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ;
01939       rec-- ;
01940 #endif
01941 
01942    /*-- recombination --*/
01943 
01944    if( mode > 0 ){
01945       for( k=0 ; k < M ; k++ ){
01946          t1  = bb[k].r; t2  = bb[k].i;
01947          tr  = cs[k].r; ti  = cs[k].i;
01948          bbr = tr*t1  - ti*t2  ; bbi = tr*t2  + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */
01949 
01950          t1  = cc[  k].r; t2  = cc[  k].i;
01951          tr  = cs[2*k].r; ti  = cs[2*k].i;
01952          ccr = tr*t1  - ti*t2  ; cci = tr*t2  + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */
01953 
01954          t4 = bbr+ccr      ; t1 = t4*CC3       ;
01955          t8 = bbi+cci      ; t6 = t8*CC3       ;
01956          t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3;
01957 
01958          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01959 
01960          xc[k   ].r = aar+t4      ; xc[k   ].i = aai+t8      ;
01961          xc[k+M ].r = (aar+t1)-t2 ; xc[k+M ].i = (aai+t6)+t5 ;
01962          xc[k+M2].r = (aar+t1)+t2 ; xc[k+M2].i = (aai+t6)-t5 ;
01963       }
01964    } else {
01965       for( k=0 ; k < M ; k++ ){
01966          t1  = bb[k].r; t2  = bb[k].i;
01967          tr  = cs[k].r; ti  = cs[k].i;
01968          bbr = tr*t1  + ti*t2  ; bbi = tr*t2  - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */
01969 
01970          t1  = cc[  k].r; t2  = cc[  k].i;
01971          tr  = cs[2*k].r; ti  = cs[2*k].i;
01972          ccr = tr*t1  + ti*t2  ; cci = tr*t2  - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */
01973 
01974          t4 = bbr+ccr      ; t1 = t4*CC3       ;
01975          t8 = bbi+cci      ; t6 = t8*CC3       ;
01976          t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3;
01977 
01978          aar = aa[k].r; aai = aa[k].i;                  /* a[k] */
01979 
01980          xc[k   ].r = aar+t4      ; xc[k   ].i = aai+t8      ;
01981          xc[k+M ].r = (aar+t1)+t2 ; xc[k+M ].i = (aai+t6)-t5 ;
01982          xc[k+M2].r = (aar+t1)-t2 ; xc[k+M2].i = (aai+t6)+t5 ;
01983       }
01984    }
01985 
01986    return ;
01987 }
01988 
01989 /*------------------------------------------------------------------
01990    fft_5dec: do a decimation-by-5, plus a power-of-2.
01991    At most RMAX levels of recursion are allowed, which means
01992    that at most 5**RMAX can be a factor of idim.
01993 --------------------------------------------------------------------*/
01994 
01995 #undef  COS72
01996 #undef  SIN72
01997 #define COS72  0.30901699   /* cos(72 deg) */
01998 #define SIN72  0.95105652   /* sin(72 deg) */
01999 
02000 static void fft_5dec( int mode , int idim , complex * xc )
02001 {
02002    static int rec=0 ;
02003    static int rmold[RMAX] = {-1,-1,-1} ;
02004    static complex *rcs[RMAX] = {NULL,NULL,NULL} ;
02005    static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX], *rdd[RMAX], *ree[RMAX] ;
02006 
02007    int N=idim , M=idim/5 , M2=2*M , M3=3*M , M4=4*M ;
02008    int mold ;
02009    complex * cs, *aa, *bb, *cc, *dd, *ee ;
02010    register int k , i ;
02011    register float aar,aai,bbr,bbi,ccr,cci,ddr,ddi,eer,eei ;
02012    register float akp,akm,bkp,bkm , ajp,ajm,bjp,bjm ;
02013    register float ak,bk,aj,bj ;
02014    float c72 , s72 , c2,s2 ;
02015    int ss ;
02016 
02017    /*-- initialize cosine/sine table and memory space --*/
02018 
02019    if( rec >= RMAX ){
02020       fprintf(stderr,"\n*** fft_5dec: too many recursions!\n"); EXIT(1);
02021    }
02022 
02023    mold = rmold[rec] ;  /* what was done before at this recursion level */
02024 
02025    if( M != mold ){     /* new value of M? */
02026       double th = (2.0*PI/N) ;
02027       if( M > mold ){             /* need more space */
02028          if( rcs[rec] != NULL ){
02029             free(rcs[rec]);free(raa[rec]);free(rbb[rec]);
02030             free(rcc[rec]);free(rdd[rec]);free(ree[rec]);
02031          }
02032          rcs[rec] = (complex *) malloc(sizeof(complex)*M4) ;
02033          raa[rec] = (complex *) malloc(sizeof(complex)*M ) ;
02034          rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ;
02035          rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ;
02036          rdd[rec] = (complex *) malloc(sizeof(complex)*M ) ;
02037          ree[rec] = (complex *) malloc(sizeof(complex)*M ) ;
02038       }
02039 
02040       rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ;
02041       for( k=1 ; k < M4 ; k++ ){
02042          rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th);
02043       }
02044       rmold[rec] = M ;  /* you must remember this */
02045    }
02046 
02047    cs = rcs[rec] ; aa = raa[rec] ; bb = rbb[rec] ;
02048                    cc = rcc[rec] ; dd = rdd[rec] ; ee = ree[rec] ;
02049 
02050    /*-- load subarrays, and FFT each one --*/
02051 
02052    for( i=k=0; k < M; k++ ){
02053       aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; ee[k]=xc[i++];
02054    }
02055 
02056 #ifndef DONT_UNROLL_FFTS
02057    switch( M ){
02058       case   1: break ;
02059 
02060       case   2: fft2  (mode,aa) ; fft2  (mode,bb) ; fft2  (mode,cc) ;
02061                 fft2  (mode,dd) ; fft2  (mode,ee) ;                   break;
02062 
02063       case   4: fft4  (mode,aa) ; fft4  (mode,bb) ; fft4  (mode,cc) ;
02064                 fft4  (mode,dd) ; fft4  (mode,ee) ;                   break;
02065 
02066       case   8: fft8  (mode,aa) ; fft8  (mode,bb) ; fft8  (mode,cc) ;
02067                 fft8  (mode,dd) ; fft8  (mode,ee) ;                   break;
02068 
02069       case  16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ;
02070                 fft16 (mode,dd) ; fft16 (mode,ee) ;                   break;
02071 
02072       case  32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ;
02073                 fft32 (mode,dd) ; fft32 (mode,ee) ;                   break;
02074 
02075       case  64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ;
02076                 fft64 (mode,dd) ; fft64 (mode,ee) ;                   break;
02077 
02078       case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ;
02079                 fft128(mode,dd) ; fft128(mode,ee) ;                   break;
02080 
02081       case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ;
02082                 fft256(mode,dd) ; fft256(mode,ee) ;                   break;
02083 
02084       case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ;
02085                 fft512(mode,dd) ; fft512(mode,ee) ;                   break;
02086 
02087       case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ;
02088                  fft1024(mode,dd) ; fft1024(mode,ee) ;                   break;
02089 
02090       case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ;
02091                  fft2048(mode,dd) ; fft2048(mode,ee) ;                   break;
02092 
02093       default:  rec++ ;
02094                 csfft_cox(mode,M,aa) ;
02095                 csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ;
02096                 csfft_cox(mode,M,dd) ; csfft_cox(mode,M,ee) ;
02097                 rec-- ; break ;
02098    }
02099 #else
02100       rec++ ;
02101       csfft_cox(mode,M,aa) ; csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ;
02102       csfft_cox(mode,M,dd) ; csfft_cox(mode,M,ee) ;
02103       rec-- ;
02104 #endif
02105 
02106    /*-- recombination --*/
02107 
02108    if( mode > 0 ){ s72 =  SIN72 ; ss =  1 ; }
02109    else          { s72 = -SIN72 ; ss = -1 ; }
02110 
02111    c72 = COS72 ;
02112    c2  = c72 * c72 - s72 * s72 ;
02113    s2  = 2.0 * c72 * s72 ;
02114 
02115    for( k=0 ; k < M ; k++ ){
02116       ak  = bb[k].r ; bk  = bb[k].i      ;
02117       aj  = cs[k].r ; bj  = cs[k].i * ss ;
02118       bbr = aj*ak - bj*bk ; bbi = aj*bk + bj*ak  ;  /* b[k]*exp(ss*2*Pi*k/N) */
02119 
02120       ak  = cc[  k].r ; bk  = cc[  k].i      ;
02121       aj  = cs[2*k].r ; bj  = cs[2*k].i * ss ;
02122       ccr = aj*ak - bj*bk ; cci = aj*bk + bj*ak  ;  /* c[k]*exp(ss*4*Pi*k/N) */
02123 
02124       ak  = dd[  k].r ; bk  = dd[  k].i      ;
02125       aj  = cs[3*k].r ; bj  = cs[3*k].i * ss ;
02126       ddr = aj*ak - bj*bk ; ddi = aj*bk + bj*ak  ;  /* d[k]*exp(ss*6*Pi*k/N) */
02127 
02128       ak  = ee[  k].r ; bk  = ee[  k].i      ;
02129       aj  = cs[4*k].r ; bj  = cs[4*k].i * ss ;
02130       eer = aj*ak - bj*bk ; eei = aj*bk + bj*ak  ;  /* e[k]*exp(ss*8*Pi*k/N) */
02131 
02132       aar = aa[k].r ; aai = aa[k].i ;               /* a[k] */
02133 
02134       /* the code below is (heavily) adapted from fftn.c */
02135 
02136       akp = bbr + eer ; akm = bbr - eer ;
02137       bkp = bbi + eei ; bkm = bbi - eei ;
02138       ajp = ccr + ddr ; ajm = ccr - ddr ;
02139       bjp = cci + ddi ; bjm = cci - ddi ;
02140 
02141       xc[k].r = aar + akp + ajp ; xc[k].i = aai + bkp + bjp ;
02142 
02143       ak = akp * c72 + ajp * c2 + aar ; bk = bkp * c72 + bjp * c2 + aai ;
02144       aj = akm * s72 + ajm * s2       ; bj = bkm * s72 + bjm * s2       ;
02145 
02146       xc[k+M ].r = ak - bj ; xc[k+M ].i = bk + aj ;
02147       xc[k+M4].r = ak + bj ; xc[k+M4].i = bk - aj ;
02148 
02149       ak = akp * c2 + ajp * c72 + aar ; bk = bkp * c2 + bjp * c72 + aai ;
02150       aj = akm * s2 - ajm * s72       ; bj = bkm * s2 - bjm * s72       ;
02151 
02152       xc[k+M2].r = ak - bj ; xc[k+M2].i = bk + aj ;
02153       xc[k+M3].r = ak + bj ; xc[k+M3].i = bk - aj ;
02154    }
02155 
02156    return ;
02157 }
02158 
02159 /*--------------------------------------------------------------------
02160    Return the smallest integer >= idim for which we can do an FFT.
02161 ----------------------------------------------------------------------*/
02162 
02163 #define N35 ((RMAX+1)*(RMAX+1))
02164 
02165 int csfft_nextup( int idim )
02166 {
02167    static int * tf = NULL , * dn = NULL ;
02168    int ibase , ii ;
02169 
02170    /*-- build table of allowable powers of 3 and 5 [tf],
02171         the powers of 2 just less than them        [dn],
02172         and their ratios tf/dn                     [rat].
02173         Then sort tf and dn to be increasing in rat.     --*/
02174 
02175    if( tf == NULL ){
02176       int p , q , tt,ff , i=0 , j ; float * rat ; float r ;
02177 
02178       tf  = (int *)   malloc(sizeof(int)  *N35) ;
02179       dn  = (int *)   malloc(sizeof(int)  *N35) ;
02180       rat = (float *) malloc(sizeof(float)*N35) ;
02181 
02182       /* create tables */
02183 
02184       for( p=0,tt=1 ; p <= RMAX ; p++,tt*=3 ){         /* tt = 3^p */
02185          for( q=0,ff=1 ; q <= RMAX ; q++,ff*=5,i++ ){  /* ff = 5^q */
02186             tf[i] = tt * ff ;
02187 
02188             j=2; while( j < tf[i] ){ j*=2; } /* j = power of 2 just >= tf */
02189             dn[i] = j/2 ;                    /* dn = power of 2 just < tf */
02190 
02191             rat[i] = (float)(tf[i]) / (float)(dn[i]) ;
02192          }
02193       }
02194 
02195       /* sort on rat (crude, but it works) */
02196 
02197       do{
02198          for( i=1,p=0 ; i < N35 ; i++ ){
02199             if( rat[i-1] > rat[i] ){
02200                r = rat[i-1] ; rat[i-1] = rat[i] ; rat[i] = r ;
02201                q = tf [i-1] ; tf [i-1] = tf [i] ; tf [i] = q ;
02202                q = dn [i-1] ; dn [i-1] = dn [i] ; dn [i] = q ;
02203                p++ ;
02204             }
02205          }
02206       } while( p > 0 ) ;
02207 
02208       free(rat) ;
02209    }
02210 
02211    /*-- loop on powers of 2 (ibase);
02212         we can do FFTs of size = tf*ibase/dn (note 1 < tf/dn < 2);
02213         sinc tf/dn is sorted, we're scanning in increasing sizes  --*/
02214 
02215    ibase = 1 ;
02216    while(1){
02217       if( idim <= ibase ) return ibase ;
02218 
02219       for( ii=0 ; ii < N35 ; ii++ )
02220          if( dn[ii] <= ibase && idim <= tf[ii]*ibase/dn[ii] )
02221             return tf[ii]*ibase/dn[ii] ;
02222 
02223       ibase *= 2 ;
02224    }
02225 }
02226 
02227 /*----------------------------------------------------------------------
02228    return the next legal value that has only 1 power of 3 and/or 5,
02229    and that is also even
02230 ------------------------------------------------------------------------*/
02231 
02232 int csfft_nextup_one35( int idim )
02233 {
02234    int jj = idim ;
02235    do{
02236       jj = csfft_nextup(jj) ;
02237       if( jj%9 == 0 || jj%25 == 0 || jj%2 == 1 ) jj++ ;
02238       else                                       return jj ;
02239    } while(1) ;
02240    return 0 ; /* cannot be reached */
02241 }
02242 
02243 /*----------------------------------------------------------------------
02244    return the next legal value that is even [13 May 2003 - RWCox]
02245 ------------------------------------------------------------------------*/
02246 
02247 int csfft_nextup_even( int idim )
02248 {
02249    int jj = idim ;
02250    do{
02251       jj = csfft_nextup(jj) ;
02252       if( jj%2 == 1 ) jj++ ;
02253       else            return jj ;
02254    } while(1) ;
02255    return 0 ; /* cannot be reached */
02256 }