/***************************************************************************** The functions in this file are adaptations of csfft.c that should be thread-safe for use with OpenMP. As such, this file is intended to be #include-d into a source file that uses OpenMP directives. In such a case, the following inclusions are not really necessary, but included here for reference. When using these functions, csfft_cox_OMP_SETUP() must be called once first, in order to setup some thread specific static storage. These variable names all start with 'CFO_'. -- RWCox -- Dec 2015 ***** I'm not actually sure these functions work -- BEWARE !! :( ***** ***** Instead ... try using the contents of file fftn_OMP.c ***** ******************************************************************************/ #include "mrilib.h" /* AFNI package library header */ #ifdef USE_OMP # include #endif /***=====================================================================** *** Prototypes for externally callable functions: ** *** Complex-to-complex FFT in place: ** *** mode = -1 or +1 (for inverse scaling, csfft_scale_inverse_OMP ** *** idim = dimension (power of 2, maybe with factors of 3 and 5) ** *** xc = input/output array ** *** Re-initializes itself only when idim changes from previous call. ** ***=====================================================================**/ void csfft_cox_OMP( int mode , int idim , complex *xc ) ; void csfft_scale_inverse_OMP( int scl ) ; /* scl=1 ==> force 1/N for mode=+1 **/ /* scl=0 ==> no 1/N scaling **/ /*** Aug 1999: ** *** idim can now contain factors of 3 and/or 5, up to and including ** *** 3^3 * 5^3. Some examples: ** *** 135 144 150 160 180 192 200 216 225 240 250 ** *** The routine csfft_nextup(n) returns the smallest size FFT ** *** >= n that csfft_cox_OMP() knows how to do. ** *** Note that efficiency of the different lengths can be quite ** *** different. In general, powers of 2 are fastest, but a ** *** single factor of 3 and/or 5 doesn't cause too much slowdown. ** *** The routine csfft_nextup_one35(n) returns the smallest size FFT ** *** >= n that contains at most one power of 3, at most one power ** *** of 5, and least one power of 2. In trials, these are the most ** *** time efficient. **/ /*-- Aug 1999: internal routines to do FFTs by decimation by 3, 4, 5 --*/ static void fft_3dec_OMP( int , int , complex * ) ; static void fft_4dec_OMP( int , int , complex * ) ; static void fft_5dec_OMP( int , int , complex * ) ; #undef RMAX #define RMAX 7 /* maximum levels of recursion for decimation by 3, 4, 5 */ #undef PI #define PI (3.141592653589793238462643) /* that should be good enough */ /*---------------------------------------------------------------------------*/ /* When using these OpenMP-ized functions, csfft_cox_OMP_SETUP() must be called once first! [Dec 2015] *//*-------------------------------------------------------------------------*/ static int CFO_nthr = 1 ; /* max number of threads to allow for */ static complex **CFO_csplus =NULL ; /* per-thread arrays of trig consts */ static complex **CFO_csminus=NULL ; static int *CFO_nold =NULL ; /* per-thread last dimension used */ static int *CFO_sclinv =NULL ; /* per-thread inverse scaling flag */ static int *CFO_recC =NULL ; /* per-thread recursion level for csfft_cox_OMP */ static int *CFO_recX =NULL ; /* per-thread recursion level for decimaion funcs */ static int *CFO_rmold3[RMAX] ; /* per-thread, per-recursion last M value used */ static int *CFO_rmold4[RMAX] ; /* for decimation by 3, 4, 5 functions */ static int *CFO_rmold5[RMAX] ; static int *nCFO_aa[RMAX] ; /* per-thread, per-recursion temp arrays */ static int *nCFO_bb[RMAX] ; /* nCFO_aa[r][i] is the size of the 'aa' */ static int *nCFO_cc[RMAX] ; /* array for recursion level r, thread i */ static int *nCFO_dd[RMAX] ; /* CFO_aa[r][i] is the corresponding array */ static int *nCFO_ee[RMAX] ; static int *nCFO_cs[RMAX] ; static complex **CFO_aa[RMAX] ; static complex **CFO_bb[RMAX] ; static complex **CFO_cc[RMAX] ; static complex **CFO_dd[RMAX] ; static complex **CFO_ee[RMAX] ; static complex **CFO_cs[RMAX] ; /* macro to declare and set the ithr = thread number (local variable) */ #undef DECLARE_ithr #ifdef USE_OMP # define DECLARE_ithr int ithr=omp_get_thread_num() #else # define DECLARE_ithr int ithr=0 #endif /*----------------------------------------------------------------------------*/ /* Allocate the CFO arrays, just one time. */ void csfft_cox_OMP_SETUP(void) { int rr ; if( CFO_nold != NULL ) return ; /* was called before! */ #ifdef USE_OMP CFO_nthr = omp_get_max_threads() ; #endif /* arrays that don't depend on recursion level, just thread number */ CFO_nold = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_sclinv = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_recC = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_recX = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_csplus = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_csminus = (complex **)calloc(sizeof(complex *),CFO_nthr) ; /* arrays that depend on recursion level and on thread number; note that we are just allocating pointer space now; the actual arrays are allocated by the CFO_assign macro (infra) */ for( rr=0 ; rr < RMAX ; rr++ ){ nCFO_aa[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; nCFO_bb[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; nCFO_cc[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; nCFO_dd[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; nCFO_ee[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; nCFO_cs[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_aa[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_bb[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_cc[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_dd[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_ee[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_cs[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ; CFO_rmold3[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_rmold4[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; CFO_rmold5[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ; } return ; } /** this macro assigns a workspace of length ll (complex), at recursion level rec (0..RMAX-1), to pointer named 'ab', from static array CFO_'ab'[rr][ithr], in thread number ithr. ithr and rec are local to the function whence this is called! The name 'ab' can be one of 'aa' 'bb' 'cc' 'dd' 'ee' 'cs'. The C preprocessor token catenation operator ## is used here. */ #undef CFO_assign #define CFO_assign(ab,ll) \ do{ if( nCFO_##ab[rec][ithr] < (ll) ){ \ CFO_##ab[rec][ithr] = (complex *)realloc(CFO_##ab[rec][ithr], \ sizeof(complex)*(ll)); \ nCFO_##ab[rec][ithr] = (ll) ; \ } \ ab = CFO_##ab[rec][ithr] ; \ } while(0) /*-------------- To order the 1/N scaling on inversion: Aug 1999 ------------*/ void csfft_scale_inverse_OMP( int scl ){ DECLARE_ithr ; CFO_sclinv[ithr] = scl; return; } /*-------------- For the unrolled FFT routines: November 1998 --------------*/ static void fft8 ( int mode , complex *xc ) ; /* completely */ static void fft16( int mode , complex *xc ) ; /* unrolled */ static void fft32( int mode , complex *xc ) ; /* routines */ static void fft64( int mode , complex *xc ) ; /* convenience macros for decimation by 4 */ #undef DONT_USE_4DEC #ifndef DONT_USE_4DEC #define fft128(m,x) fft_4dec_OMP(m, 128,x) /* via fft32 (non-recursive) */ #define fft256(m,x) fft_4dec_OMP(m, 256,x) /* via fft64 (non-recursive) */ #define fft512(m,x) fft_4dec_OMP(m, 512,x) /* via fft128 (recursive X 1) */ #define fft1024(m,x) fft_4dec_OMP(m, 1024,x) /* via fft256 (recursive X 1) */ #if 1 #define fft2048(m,x) fft_4dec_OMP(m, 2048,x) /* via fft512 (recursive X 2) */ #define fft4096(m,x) fft_4dec_OMP(m, 4096,x) /* via fft1024 (recursive X 2) */ #define fft8192(m,x) fft_4dec_OMP(m, 8192,x) /* via fft2048 (recursive X 3) */ #define fft16384(m,x) fft_4dec_OMP(m,16384,x) /* via fft4096 (recursive X 3) */ #define fft32768(m,x) fft_4dec_OMP(m,32768,x) /* via fft8192 (recursive X 4) */ #define fft65536(m,x) fft_4dec_OMP(m,65536,x) /* via fft16384 (recursive X 4) */ #endif #else #define fft128(m,x) csfft_cox_OMP(m, 128,x) #define fft256(m,x) csfft_cox_OMP(m, 256,x) #define fft512(m,x) csfft_cox_OMP(m, 512,x) #define fft1024(m,x) csfft_cox_OMP(m, 1024,x) #if 1 #define fft2048(m,x) csfft_cox_OMP(m, 2048,x) #define fft4096(m,x) csfft_cox_OMP(m, 4096,x) #define fft8192(m,x) csfft_cox_OMP(m, 8192,x) #define fft16384(m,x) csfft_cox_OMP(m,16384,x) #define fft32768(m,x) csfft_cox_OMP(m,32768,x) #define fft65536(m,x) csfft_cox_OMP(m,65536,x) #endif #endif /*-- do FFTs of size 2 and 4 with macros, not functions --*/ #define fft2(m,x) do{ float a,b,c,d ; \ a = x[0].r + x[1].r ; b = x[0].i + x[1].i ; \ c = x[0].r - x[1].r ; d = x[0].i - x[1].i ; \ x[0].r = a ; x[0].i = b ; \ x[1].r = c ; x[1].i = d ; } while(0) #define fft4(m,x) do{ float acpr,acmr,bdpr,bdmr, acpi,acmi,bdpi,bdmi; \ acpr = x[0].r + x[2].r; acmr = x[0].r - x[2].r; \ bdpr = x[1].r + x[3].r; bdmr = x[1].r - x[3].r; \ acpi = x[0].i + x[2].i; acmi = x[0].i - x[2].i; \ bdpi = x[1].i + x[3].i; bdmi = x[1].i - x[3].i; \ x[0].r = acpr+bdpr ; x[0].i = acpi+bdpi ; \ x[2].r = acpr-bdpr ; x[2].i = acpi-bdpi ; \ if(m > 0){ \ x[1].r = acmr-bdmi ; x[1].i = acmi+bdmr ; \ x[3].r = acmr+bdmi ; x[3].i = acmi-bdmr ; \ } else { \ x[1].r = acmr+bdmi ; x[1].i = acmi-bdmr ; \ x[3].r = acmr-bdmi ; x[3].i = acmi+bdmr ; \ } } while(0) /*-------------------------------------------------------------------- Initialize csfft trig constants table. Adapted from AJ's code. Does both mode=1 and mode=-1 tables. (mode=+1 == inversion) ----------------------------------------------------------------------*/ static void csfft_trigconsts( int idim ) /* internal function */ { DECLARE_ithr ; register unsigned int m, n, i0, i1, i2, i3, k; register complex *r0, *csp; register float co, si, f0, f1, f2, f3, f4; double al; complex *my_csp , *my_csm ; if( idim == CFO_nold[ithr] ) return ; /* previously done */ { if( idim > CFO_nold[ithr] ){ CFO_csplus[ithr] = (complex *)realloc( CFO_csplus[ithr] , sizeof(complex)*idim ) ; CFO_csminus[ithr] = (complex *)realloc( CFO_csminus[ithr], sizeof(complex)*idim ) ; if( CFO_csplus[ithr] == NULL || CFO_csminus[ithr] == NULL ){ fprintf(stderr,"\n*** csfft_cox_OMP cannot realloc space idim=%d! ***\n",idim); EXIT(1) ; } }} CFO_nold[ithr] = n = idim ; my_csp = CFO_csplus[ithr] ; my_csm = CFO_csminus[ithr] ; f1 = 1.0 ; /* csplus init */ m = 1; k = 0; while (n > m) { i3 = m << 1; f2 = m; al = f1*PI/f2; co = cos(al); si = sin(al); (my_csp + k)->r = 1.; (my_csp + k)->i = 0.; for (i0=0; i0 < m; i0++) { k++; csp = my_csp + k; r0 = csp - 1; csp->r = r0->r * co - r0->i * si; csp->i = r0->i * co + r0->r * si; } m = i3; } f1 = -1.0 ; /* csminus init */ m = 1; k = 0; while (n > m) { i3 = m << 1; f2 = m; al = f1*PI/f2; co = cos(al); si = sin(al); (my_csm + k)->r = 1.; (my_csm + k)->i = 0.; for (i0=0; i0 < m; i0++) { k++; csp = my_csm + k; r0 = csp - 1; csp->r = r0->r * co - r0->i * si; csp->i = r0->i * co + r0->r * si; } m = i3; } return ; } /*-------------------------------------------------------------------- Complex-to-complex FFT in place: mode = -1 or +1 (1/idim scaling on inverse [+1] if CFO_sclinv was set in the call to csfft_scale_inverse_OMP) idim = dimension (power of 2, possibly with a factor of 3^p * 5^q also, for p,q=0..RMAX) xc = input/output array Automagically re-initializes itself when idim changes from previous call. By AJ (Andrzej Jesmanowicz), modified by RWCox. ----------------------------------------------------------------------*/ /*- Macro to do 1/N scaling on inverse: if it is ordered by the user, and if mode is positive, and if this is not a recursive call. -*/ #define SCLINV \ if( CFO_sclinv[ithr] && mode > 0 && rec == 0 ){ \ register int qq ; register float ff = 1.0 / (float) idim ; \ for( qq=0 ; qq < idim ; qq++ ){ xc[qq].r *= ff ; xc[qq].i *= ff ; } } void csfft_cox_OMP( int mode , int idim , complex *xc ) { DECLARE_ithr ; register unsigned int m, n, i0, i1, i2, i3, k; register complex *r0, *r1, *csp; register float co, si, f0, f1, f2, f3, f4; int rec = CFO_recC[ithr] ; if( idim <= 1 ) return ; /* stoopid inpoot */ /*-- November 1998: maybe use the unrolled FFT routines --*/ /* INFO_message("csfft_cox_OMP(%d)",idim) ; */ switch( idim ){ /* none of these */ case 1: return; /* routines will */ case 2: fft2 (mode,xc); SCLINV; return; /* call csfft_cox */ case 4: fft4 (mode,xc); SCLINV; return; /* so don't need */ case 8: fft8 (mode,xc); SCLINV; return; /* to do rec++/-- */ case 16: fft16 (mode,xc); SCLINV; return; case 32: fft32 (mode,xc); SCLINV; return; case 64: fft64 (mode,xc); SCLINV; return; #ifndef DONT_USE_4DEC case 128: fft128 (mode,xc); SCLINV; return; case 256: fft256 (mode,xc); SCLINV; return; case 512: fft512 (mode,xc); SCLINV; return; case 1024: fft1024 (mode,xc); SCLINV; return; #if 1 case 2048: fft2048 (mode,xc); SCLINV; return; case 4096: fft4096 (mode,xc); SCLINV; return; case 8192: fft8192 (mode,xc); SCLINV; return; case 16384: fft16384(mode,xc); SCLINV; return; case 32768: fft32768(mode,xc); SCLINV; return; case 65536: fft65536(mode,xc); SCLINV; return; #endif #endif } /* fall thru to here ==> not one of the shorter 2^p cases above! */ /*-- Aug 1999: deal with factors of 3 or 5 [might call csfft_cox_OMP again] --*/ if( idim%3 == 0 ){CFO_recC[ithr]++; fft_3dec_OMP(mode,idim,xc); CFO_recC[ithr]--; SCLINV; return;} if( idim%5 == 0 ){CFO_recC[ithr]++; fft_5dec_OMP(mode,idim,xc); CFO_recC[ithr]--; SCLINV; return;} /*-- below here: the old general power of 2 routine (i.e., for idim > 65536) --*/ /**-- perhaps initialize trig consts (for this thread) --**/ if( CFO_nold[ithr] != idim ) csfft_trigconsts( idim ) ; /** Main loop starts here **/ n = idim; i2 = idim >> 1; i1 = 0; csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */ /*-- swap data --*/ for (i0=0; i0 < n; i0 ++) { if ( i1 > i0 ) { r0 = xc + i0; r1 = xc + i1; f1 = r0->r; f2 = r0->i; r0->r = r1->r; r0->i = r1->i; r1->r = f1; r1->i = f2; } m = i2; while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; } i1 += m; } /*-- compute compute compute --*/ m = 1; k = 0; while (n > m) { i3 = m << 1; for (i0=0; i0 < m; i0 ++) { co = (csp + k)->r; si = (csp + k)->i; for (i1=i0; i1 < n; i1 += i3) { r0 = xc + i1; r1 = r0 + m; f1 = r1->r * co; f2 = r1->i * si; f3 = r1->r * si; f4 = r1->i * co; f1 -= f2; f3 += f4; f2 = r0->r; f4 = r0->i; r1->r = f2 - f1; r1->i = f4 - f3; r0->r = f2 + f1; r0->i = f4 + f3; } k++; } m = i3; } SCLINV ; return ; } /****************************************************************************** Routines below for unrolled FFT routines; generated with program fftprint.c. *******************************************************************************/ /**************************************/ /* FFT routine unrolled of length 8 */ static void fft8( int mode , complex *xc ) { DECLARE_ithr ; register complex *csp , *xcx=xc; register float f1,f2,f3,f4 ; /* ININFO_message("fft8") ; */ /** perhaps initialize **/ if( CFO_nold[ithr] != 8 ) csfft_trigconsts( 8 ) ; csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */ /** data swapping part **/ f1 = xcx[1].r ; f2 = xcx[1].i ; xcx[1].r = xcx[4].r ; xcx[1].i = xcx[4].i ; xcx[4].r = f1 ; xcx[4].i = f2 ; f1 = xcx[3].r ; f2 = xcx[3].i ; xcx[3].r = xcx[6].r ; xcx[3].i = xcx[6].i ; xcx[6].r = f1 ; xcx[6].i = f2 ; /** butterflying part **/ f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */ f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */ f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */ f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; return ; } /**************************************/ /* FFT routine unrolled of length 16 */ static void fft16( int mode , complex *xc ) { DECLARE_ithr ; register complex *csp , *xcx=xc; register float f1,f2,f3,f4 ; /* ININFO_message("fft16") ; */ /** perhaps initialize **/ if( CFO_nold[ithr] != 16 ) csfft_trigconsts( 16 ) ; csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */ /** data swapping part **/ f1 = xcx[1].r ; f2 = xcx[1].i ; xcx[1].r = xcx[8].r ; xcx[1].i = xcx[8].i ; xcx[8].r = f1 ; xcx[8].i = f2 ; f1 = xcx[2].r ; f2 = xcx[2].i ; xcx[2].r = xcx[4].r ; xcx[2].i = xcx[4].i ; xcx[4].r = f1 ; xcx[4].i = f2 ; f1 = xcx[3].r ; f2 = xcx[3].i ; xcx[3].r = xcx[12].r ; xcx[3].i = xcx[12].i ; xcx[12].r = f1 ; xcx[12].i = f2 ; f1 = xcx[5].r ; f2 = xcx[5].i ; xcx[5].r = xcx[10].r ; xcx[5].i = xcx[10].i ; xcx[10].r = f1 ; xcx[10].i = f2 ; f1 = xcx[7].r ; f2 = xcx[7].i ; xcx[7].r = xcx[14].r ; xcx[7].i = xcx[14].i ; xcx[14].r = f1 ; xcx[14].i = f2 ; f1 = xcx[11].r ; f2 = xcx[11].i ; xcx[11].r = xcx[13].r ; xcx[11].i = xcx[13].i ; xcx[13].r = f1 ; xcx[13].i = f2 ; /** butterflying part **/ f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */ f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */ f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */ f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */ f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[11].r * csp[2].r - xcx[11].i * csp[2].i ; /* twiddles */ f3 = xcx[11].r * csp[2].i + xcx[11].i * csp[2].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[15].r * csp[2].r - xcx[15].i * csp[2].i ; /* twiddles */ f3 = xcx[15].r * csp[2].i + xcx[15].i * csp[2].r ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */ f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */ f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[14].r * csp[5].r - xcx[14].i * csp[5].i ; /* twiddles */ f3 = xcx[14].r * csp[5].i + xcx[14].i * csp[5].r ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */ f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */ f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */ f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */ f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[12].r * csp[11].r - xcx[12].i * csp[11].i ; /* twiddles */ f3 = xcx[12].r * csp[11].i + xcx[12].i * csp[11].r ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */ f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */ f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */ f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; return ; } /**************************************/ /* FFT routine unrolled of length 32 */ static void fft32( int mode , complex *xc ) { DECLARE_ithr ; register complex *csp , *xcx=xc; register float f1,f2,f3,f4 ; /* ININFO_message("fft32") ; */ /** perhaps initialize **/ if( CFO_nold[ithr] != 32 ) csfft_trigconsts( 32 ) ; csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */ /** data swapping part **/ f1 = xcx[1].r ; f2 = xcx[1].i ; xcx[1].r = xcx[16].r ; xcx[1].i = xcx[16].i ; xcx[16].r = f1 ; xcx[16].i = f2 ; f1 = xcx[2].r ; f2 = xcx[2].i ; xcx[2].r = xcx[8].r ; xcx[2].i = xcx[8].i ; xcx[8].r = f1 ; xcx[8].i = f2 ; f1 = xcx[3].r ; f2 = xcx[3].i ; xcx[3].r = xcx[24].r ; xcx[3].i = xcx[24].i ; xcx[24].r = f1 ; xcx[24].i = f2 ; f1 = xcx[5].r ; f2 = xcx[5].i ; xcx[5].r = xcx[20].r ; xcx[5].i = xcx[20].i ; xcx[20].r = f1 ; xcx[20].i = f2 ; f1 = xcx[6].r ; f2 = xcx[6].i ; xcx[6].r = xcx[12].r ; xcx[6].i = xcx[12].i ; xcx[12].r = f1 ; xcx[12].i = f2 ; f1 = xcx[7].r ; f2 = xcx[7].i ; xcx[7].r = xcx[28].r ; xcx[7].i = xcx[28].i ; xcx[28].r = f1 ; xcx[28].i = f2 ; f1 = xcx[9].r ; f2 = xcx[9].i ; xcx[9].r = xcx[18].r ; xcx[9].i = xcx[18].i ; xcx[18].r = f1 ; xcx[18].i = f2 ; f1 = xcx[11].r ; f2 = xcx[11].i ; xcx[11].r = xcx[26].r ; xcx[11].i = xcx[26].i ; xcx[26].r = f1 ; xcx[26].i = f2 ; f1 = xcx[13].r ; f2 = xcx[13].i ; xcx[13].r = xcx[22].r ; xcx[13].i = xcx[22].i ; xcx[22].r = f1 ; xcx[22].i = f2 ; f1 = xcx[15].r ; f2 = xcx[15].i ; xcx[15].r = xcx[30].r ; xcx[15].i = xcx[30].i ; xcx[30].r = f1 ; xcx[30].i = f2 ; f1 = xcx[19].r ; f2 = xcx[19].i ; xcx[19].r = xcx[25].r ; xcx[19].i = xcx[25].i ; xcx[25].r = f1 ; xcx[25].i = f2 ; f1 = xcx[23].r ; f2 = xcx[23].i ; xcx[23].r = xcx[29].r ; xcx[23].i = xcx[29].i ; xcx[29].r = f1 ; xcx[29].i = f2 ; /** butterflying part **/ f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */ f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */ f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[17].r ; f3 = xcx[17].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[19].r ; f3 = xcx[19].i ; /* cos=1 sin=0 */ f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = xcx[21].r ; f3 = xcx[21].i ; /* cos=1 sin=0 */ f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[23].r ; f3 = xcx[23].i ; /* cos=1 sin=0 */ f2 = xcx[22].r ; f4 = xcx[22].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; f1 = xcx[25].r ; f3 = xcx[25].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[27].r ; f3 = xcx[27].i ; /* cos=1 sin=0 */ f2 = xcx[26].r ; f4 = xcx[26].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; f1 = xcx[29].r ; f3 = xcx[29].i ; /* cos=1 sin=0 */ f2 = xcx[28].r ; f4 = xcx[28].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; f1 = xcx[31].r ; f3 = xcx[31].i ; /* cos=1 sin=0 */ f2 = xcx[30].r ; f4 = xcx[30].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ; f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[18].r ; f3 = xcx[18].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[22].r ; f3 = xcx[22].i ; /* cos=1 sin=0 */ f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[26].r ; f3 = xcx[26].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[30].r ; f3 = xcx[30].i ; /* cos=1 sin=0 */ f2 = xcx[28].r ; f4 = xcx[28].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[3].r * csp[2].i ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[7].r * csp[2].i ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[11].r * csp[2].i ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[15].r * csp[2].i ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[19].r * csp[2].i ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[23].r * csp[2].i ; f2 = xcx[21].r ; f4 = xcx[21].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[27].r * csp[2].i ; f2 = xcx[25].r ; f4 = xcx[25].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[31].r * csp[2].i ; f2 = xcx[29].r ; f4 = xcx[29].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ; f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[20].r ; f3 = xcx[20].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[28].r ; f3 = xcx[28].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */ f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */ f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */ f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ; f2 = xcx[25].r ; f4 = xcx[25].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[6].r * csp[5].i ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[14].r * csp[5].i ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[22].r * csp[5].i ; f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[30].r * csp[5].i ; f2 = xcx[26].r ; f4 = xcx[26].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */ f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */ f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ; f2 = xcx[19].r ; f4 = xcx[19].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */ f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ; f2 = xcx[27].r ; f4 = xcx[27].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ; f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[24].r ; f3 = xcx[24].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */ f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */ f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */ f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */ f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ; f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */ f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */ f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ; f2 = xcx[19].r ; f4 = xcx[19].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[12].r * csp[11].i ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[28].r * csp[11].i ; f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */ f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */ f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ; f2 = xcx[21].r ; f4 = xcx[21].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */ f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */ f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ; f2 = xcx[22].r ; f4 = xcx[22].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */ f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */ f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ; f2 = xcx[23].r ; f4 = xcx[23].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ; f1 = xcx[16].r ; f3 = xcx[16].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */ f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */ f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */ f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */ f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */ f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */ f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */ f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */ f3 = xcx[24].r * csp[23].i ; f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */ f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */ f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */ f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */ f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ; f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */ f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */ f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ; f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */ f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ; f2 = xcx[15].r ; f4 = xcx[15].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ; return ; } /**************************************/ /* FFT routine unrolled of length 64 */ void fft64( int mode , complex * xc ) { DECLARE_ithr ; register complex *csp , *xcx=xc; register float f1,f2,f3,f4 ; /* ININFO_message("fft64") ; */ /** perhaps initialize **/ if( CFO_nold[ithr] != 64 ) csfft_trigconsts( 64 ) ; csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */ /** data swapping part **/ f1 = xcx[1].r ; f2 = xcx[1].i ; xcx[1].r = xcx[32].r ; xcx[1].i = xcx[32].i ; xcx[32].r = f1 ; xcx[32].i = f2 ; f1 = xcx[2].r ; f2 = xcx[2].i ; xcx[2].r = xcx[16].r ; xcx[2].i = xcx[16].i ; xcx[16].r = f1 ; xcx[16].i = f2 ; f1 = xcx[3].r ; f2 = xcx[3].i ; xcx[3].r = xcx[48].r ; xcx[3].i = xcx[48].i ; xcx[48].r = f1 ; xcx[48].i = f2 ; f1 = xcx[4].r ; f2 = xcx[4].i ; xcx[4].r = xcx[8].r ; xcx[4].i = xcx[8].i ; xcx[8].r = f1 ; xcx[8].i = f2 ; f1 = xcx[5].r ; f2 = xcx[5].i ; xcx[5].r = xcx[40].r ; xcx[5].i = xcx[40].i ; xcx[40].r = f1 ; xcx[40].i = f2 ; f1 = xcx[6].r ; f2 = xcx[6].i ; xcx[6].r = xcx[24].r ; xcx[6].i = xcx[24].i ; xcx[24].r = f1 ; xcx[24].i = f2 ; f1 = xcx[7].r ; f2 = xcx[7].i ; xcx[7].r = xcx[56].r ; xcx[7].i = xcx[56].i ; xcx[56].r = f1 ; xcx[56].i = f2 ; f1 = xcx[9].r ; f2 = xcx[9].i ; xcx[9].r = xcx[36].r ; xcx[9].i = xcx[36].i ; xcx[36].r = f1 ; xcx[36].i = f2 ; f1 = xcx[10].r ; f2 = xcx[10].i ; xcx[10].r = xcx[20].r ; xcx[10].i = xcx[20].i ; xcx[20].r = f1 ; xcx[20].i = f2 ; f1 = xcx[11].r ; f2 = xcx[11].i ; xcx[11].r = xcx[52].r ; xcx[11].i = xcx[52].i ; xcx[52].r = f1 ; xcx[52].i = f2 ; f1 = xcx[13].r ; f2 = xcx[13].i ; xcx[13].r = xcx[44].r ; xcx[13].i = xcx[44].i ; xcx[44].r = f1 ; xcx[44].i = f2 ; f1 = xcx[14].r ; f2 = xcx[14].i ; xcx[14].r = xcx[28].r ; xcx[14].i = xcx[28].i ; xcx[28].r = f1 ; xcx[28].i = f2 ; f1 = xcx[15].r ; f2 = xcx[15].i ; xcx[15].r = xcx[60].r ; xcx[15].i = xcx[60].i ; xcx[60].r = f1 ; xcx[60].i = f2 ; f1 = xcx[17].r ; f2 = xcx[17].i ; xcx[17].r = xcx[34].r ; xcx[17].i = xcx[34].i ; xcx[34].r = f1 ; xcx[34].i = f2 ; f1 = xcx[19].r ; f2 = xcx[19].i ; xcx[19].r = xcx[50].r ; xcx[19].i = xcx[50].i ; xcx[50].r = f1 ; xcx[50].i = f2 ; f1 = xcx[21].r ; f2 = xcx[21].i ; xcx[21].r = xcx[42].r ; xcx[21].i = xcx[42].i ; xcx[42].r = f1 ; xcx[42].i = f2 ; f1 = xcx[22].r ; f2 = xcx[22].i ; xcx[22].r = xcx[26].r ; xcx[22].i = xcx[26].i ; xcx[26].r = f1 ; xcx[26].i = f2 ; f1 = xcx[23].r ; f2 = xcx[23].i ; xcx[23].r = xcx[58].r ; xcx[23].i = xcx[58].i ; xcx[58].r = f1 ; xcx[58].i = f2 ; f1 = xcx[25].r ; f2 = xcx[25].i ; xcx[25].r = xcx[38].r ; xcx[25].i = xcx[38].i ; xcx[38].r = f1 ; xcx[38].i = f2 ; f1 = xcx[27].r ; f2 = xcx[27].i ; xcx[27].r = xcx[54].r ; xcx[27].i = xcx[54].i ; xcx[54].r = f1 ; xcx[54].i = f2 ; f1 = xcx[29].r ; f2 = xcx[29].i ; xcx[29].r = xcx[46].r ; xcx[29].i = xcx[46].i ; xcx[46].r = f1 ; xcx[46].i = f2 ; f1 = xcx[31].r ; f2 = xcx[31].i ; xcx[31].r = xcx[62].r ; xcx[31].i = xcx[62].i ; xcx[62].r = f1 ; xcx[62].i = f2 ; f1 = xcx[35].r ; f2 = xcx[35].i ; xcx[35].r = xcx[49].r ; xcx[35].i = xcx[49].i ; xcx[49].r = f1 ; xcx[49].i = f2 ; f1 = xcx[37].r ; f2 = xcx[37].i ; xcx[37].r = xcx[41].r ; xcx[37].i = xcx[41].i ; xcx[41].r = f1 ; xcx[41].i = f2 ; f1 = xcx[39].r ; f2 = xcx[39].i ; xcx[39].r = xcx[57].r ; xcx[39].i = xcx[57].i ; xcx[57].r = f1 ; xcx[57].i = f2 ; f1 = xcx[43].r ; f2 = xcx[43].i ; xcx[43].r = xcx[53].r ; xcx[43].i = xcx[53].i ; xcx[53].r = f1 ; xcx[53].i = f2 ; f1 = xcx[47].r ; f2 = xcx[47].i ; xcx[47].r = xcx[61].r ; xcx[47].i = xcx[61].i ; xcx[61].r = f1 ; xcx[61].i = f2 ; f1 = xcx[55].r ; f2 = xcx[55].i ; xcx[55].r = xcx[59].r ; xcx[55].i = xcx[59].i ; xcx[59].r = f1 ; xcx[59].i = f2 ; /** butterflying part **/ f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */ f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */ f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[17].r ; f3 = xcx[17].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[19].r ; f3 = xcx[19].i ; /* cos=1 sin=0 */ f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = xcx[21].r ; f3 = xcx[21].i ; /* cos=1 sin=0 */ f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[23].r ; f3 = xcx[23].i ; /* cos=1 sin=0 */ f2 = xcx[22].r ; f4 = xcx[22].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; f1 = xcx[25].r ; f3 = xcx[25].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[27].r ; f3 = xcx[27].i ; /* cos=1 sin=0 */ f2 = xcx[26].r ; f4 = xcx[26].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; f1 = xcx[29].r ; f3 = xcx[29].i ; /* cos=1 sin=0 */ f2 = xcx[28].r ; f4 = xcx[28].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; f1 = xcx[31].r ; f3 = xcx[31].i ; /* cos=1 sin=0 */ f2 = xcx[30].r ; f4 = xcx[30].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ; f1 = xcx[33].r ; f3 = xcx[33].i ; /* cos=1 sin=0 */ f2 = xcx[32].r ; f4 = xcx[32].i ; xcx[33].r = f2-f1 ; xcx[33].i = f4-f3 ; xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ; f1 = xcx[35].r ; f3 = xcx[35].i ; /* cos=1 sin=0 */ f2 = xcx[34].r ; f4 = xcx[34].i ; xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ; xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ; f1 = xcx[37].r ; f3 = xcx[37].i ; /* cos=1 sin=0 */ f2 = xcx[36].r ; f4 = xcx[36].i ; xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ; xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ; f1 = xcx[39].r ; f3 = xcx[39].i ; /* cos=1 sin=0 */ f2 = xcx[38].r ; f4 = xcx[38].i ; xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ; xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ; f1 = xcx[41].r ; f3 = xcx[41].i ; /* cos=1 sin=0 */ f2 = xcx[40].r ; f4 = xcx[40].i ; xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ; xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ; f1 = xcx[43].r ; f3 = xcx[43].i ; /* cos=1 sin=0 */ f2 = xcx[42].r ; f4 = xcx[42].i ; xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ; xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ; f1 = xcx[45].r ; f3 = xcx[45].i ; /* cos=1 sin=0 */ f2 = xcx[44].r ; f4 = xcx[44].i ; xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ; xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ; f1 = xcx[47].r ; f3 = xcx[47].i ; /* cos=1 sin=0 */ f2 = xcx[46].r ; f4 = xcx[46].i ; xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ; xcx[46].r = f2+f1 ; xcx[46].i = f4+f3 ; f1 = xcx[49].r ; f3 = xcx[49].i ; /* cos=1 sin=0 */ f2 = xcx[48].r ; f4 = xcx[48].i ; xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ; xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ; f1 = xcx[51].r ; f3 = xcx[51].i ; /* cos=1 sin=0 */ f2 = xcx[50].r ; f4 = xcx[50].i ; xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ; xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ; f1 = xcx[53].r ; f3 = xcx[53].i ; /* cos=1 sin=0 */ f2 = xcx[52].r ; f4 = xcx[52].i ; xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ; xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ; f1 = xcx[55].r ; f3 = xcx[55].i ; /* cos=1 sin=0 */ f2 = xcx[54].r ; f4 = xcx[54].i ; xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ; xcx[54].r = f2+f1 ; xcx[54].i = f4+f3 ; f1 = xcx[57].r ; f3 = xcx[57].i ; /* cos=1 sin=0 */ f2 = xcx[56].r ; f4 = xcx[56].i ; xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ; xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ; f1 = xcx[59].r ; f3 = xcx[59].i ; /* cos=1 sin=0 */ f2 = xcx[58].r ; f4 = xcx[58].i ; xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ; xcx[58].r = f2+f1 ; xcx[58].i = f4+f3 ; f1 = xcx[61].r ; f3 = xcx[61].i ; /* cos=1 sin=0 */ f2 = xcx[60].r ; f4 = xcx[60].i ; xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ; xcx[60].r = f2+f1 ; xcx[60].i = f4+f3 ; f1 = xcx[63].r ; f3 = xcx[63].i ; /* cos=1 sin=0 */ f2 = xcx[62].r ; f4 = xcx[62].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[62].r = f2+f1 ; xcx[62].i = f4+f3 ; f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */ f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[18].r ; f3 = xcx[18].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[22].r ; f3 = xcx[22].i ; /* cos=1 sin=0 */ f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[26].r ; f3 = xcx[26].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[30].r ; f3 = xcx[30].i ; /* cos=1 sin=0 */ f2 = xcx[28].r ; f4 = xcx[28].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; f1 = xcx[34].r ; f3 = xcx[34].i ; /* cos=1 sin=0 */ f2 = xcx[32].r ; f4 = xcx[32].i ; xcx[34].r = f2-f1 ; xcx[34].i = f4-f3 ; xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ; f1 = xcx[38].r ; f3 = xcx[38].i ; /* cos=1 sin=0 */ f2 = xcx[36].r ; f4 = xcx[36].i ; xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ; xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ; f1 = xcx[42].r ; f3 = xcx[42].i ; /* cos=1 sin=0 */ f2 = xcx[40].r ; f4 = xcx[40].i ; xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ; xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ; f1 = xcx[46].r ; f3 = xcx[46].i ; /* cos=1 sin=0 */ f2 = xcx[44].r ; f4 = xcx[44].i ; xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ; xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ; f1 = xcx[50].r ; f3 = xcx[50].i ; /* cos=1 sin=0 */ f2 = xcx[48].r ; f4 = xcx[48].i ; xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ; xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ; f1 = xcx[54].r ; f3 = xcx[54].i ; /* cos=1 sin=0 */ f2 = xcx[52].r ; f4 = xcx[52].i ; xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ; xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ; f1 = xcx[58].r ; f3 = xcx[58].i ; /* cos=1 sin=0 */ f2 = xcx[56].r ; f4 = xcx[56].i ; xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ; xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ; f1 = xcx[62].r ; f3 = xcx[62].i ; /* cos=1 sin=0 */ f2 = xcx[60].r ; f4 = xcx[60].i ; xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ; xcx[60].r = f2+f1 ; xcx[60].i = f4+f3 ; f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[3].r * csp[2].i ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[7].r * csp[2].i ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[11].r * csp[2].i ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[15].r * csp[2].i ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[19].r * csp[2].i ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[23].r * csp[2].i ; f2 = xcx[21].r ; f4 = xcx[21].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[27].r * csp[2].i ; f2 = xcx[25].r ; f4 = xcx[25].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[31].r * csp[2].i ; f2 = xcx[29].r ; f4 = xcx[29].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ; f1 = - xcx[35].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[35].r * csp[2].i ; f2 = xcx[33].r ; f4 = xcx[33].i ; xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ; xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ; f1 = - xcx[39].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[39].r * csp[2].i ; f2 = xcx[37].r ; f4 = xcx[37].i ; xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ; xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ; f1 = - xcx[43].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[43].r * csp[2].i ; f2 = xcx[41].r ; f4 = xcx[41].i ; xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ; xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ; f1 = - xcx[47].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[47].r * csp[2].i ; f2 = xcx[45].r ; f4 = xcx[45].i ; xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ; xcx[45].r = f2+f1 ; xcx[45].i = f4+f3 ; f1 = - xcx[51].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[51].r * csp[2].i ; f2 = xcx[49].r ; f4 = xcx[49].i ; xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ; xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ; f1 = - xcx[55].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[55].r * csp[2].i ; f2 = xcx[53].r ; f4 = xcx[53].i ; xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ; xcx[53].r = f2+f1 ; xcx[53].i = f4+f3 ; f1 = - xcx[59].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[59].r * csp[2].i ; f2 = xcx[57].r ; f4 = xcx[57].i ; xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ; xcx[57].r = f2+f1 ; xcx[57].i = f4+f3 ; f1 = - xcx[63].i * csp[2].i ; /* cos=0 twiddles */ f3 = xcx[63].r * csp[2].i ; f2 = xcx[61].r ; f4 = xcx[61].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[61].r = f2+f1 ; xcx[61].i = f4+f3 ; f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */ f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[20].r ; f3 = xcx[20].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[28].r ; f3 = xcx[28].i ; /* cos=1 sin=0 */ f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[36].r ; f3 = xcx[36].i ; /* cos=1 sin=0 */ f2 = xcx[32].r ; f4 = xcx[32].i ; xcx[36].r = f2-f1 ; xcx[36].i = f4-f3 ; xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ; f1 = xcx[44].r ; f3 = xcx[44].i ; /* cos=1 sin=0 */ f2 = xcx[40].r ; f4 = xcx[40].i ; xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ; xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ; f1 = xcx[52].r ; f3 = xcx[52].i ; /* cos=1 sin=0 */ f2 = xcx[48].r ; f4 = xcx[48].i ; xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ; xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ; f1 = xcx[60].r ; f3 = xcx[60].i ; /* cos=1 sin=0 */ f2 = xcx[56].r ; f4 = xcx[56].i ; xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ; xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ; f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */ f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */ f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */ f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ; f2 = xcx[25].r ; f4 = xcx[25].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; f1 = xcx[37].r * csp[4].r - xcx[37].i * csp[4].i ; /* twiddles */ f3 = xcx[37].r * csp[4].i + xcx[37].i * csp[4].r ; f2 = xcx[33].r ; f4 = xcx[33].i ; xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ; xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ; f1 = xcx[45].r * csp[4].r - xcx[45].i * csp[4].i ; /* twiddles */ f3 = xcx[45].r * csp[4].i + xcx[45].i * csp[4].r ; f2 = xcx[41].r ; f4 = xcx[41].i ; xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ; xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ; f1 = xcx[53].r * csp[4].r - xcx[53].i * csp[4].i ; /* twiddles */ f3 = xcx[53].r * csp[4].i + xcx[53].i * csp[4].r ; f2 = xcx[49].r ; f4 = xcx[49].i ; xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ; xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ; f1 = xcx[61].r * csp[4].r - xcx[61].i * csp[4].i ; /* twiddles */ f3 = xcx[61].r * csp[4].i + xcx[61].i * csp[4].r ; f2 = xcx[57].r ; f4 = xcx[57].i ; xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ; xcx[57].r = f2+f1 ; xcx[57].i = f4+f3 ; f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[6].r * csp[5].i ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[14].r * csp[5].i ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[22].r * csp[5].i ; f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[30].r * csp[5].i ; f2 = xcx[26].r ; f4 = xcx[26].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; f1 = - xcx[38].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[38].r * csp[5].i ; f2 = xcx[34].r ; f4 = xcx[34].i ; xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ; xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ; f1 = - xcx[46].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[46].r * csp[5].i ; f2 = xcx[42].r ; f4 = xcx[42].i ; xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ; xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ; f1 = - xcx[54].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[54].r * csp[5].i ; f2 = xcx[50].r ; f4 = xcx[50].i ; xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ; xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ; f1 = - xcx[62].i * csp[5].i ; /* cos=0 twiddles */ f3 = xcx[62].r * csp[5].i ; f2 = xcx[58].r ; f4 = xcx[58].i ; xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ; xcx[58].r = f2+f1 ; xcx[58].i = f4+f3 ; f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */ f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */ f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ; f2 = xcx[19].r ; f4 = xcx[19].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */ f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ; f2 = xcx[27].r ; f4 = xcx[27].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ; f1 = xcx[39].r * csp[6].r - xcx[39].i * csp[6].i ; /* twiddles */ f3 = xcx[39].r * csp[6].i + xcx[39].i * csp[6].r ; f2 = xcx[35].r ; f4 = xcx[35].i ; xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ; xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ; f1 = xcx[47].r * csp[6].r - xcx[47].i * csp[6].i ; /* twiddles */ f3 = xcx[47].r * csp[6].i + xcx[47].i * csp[6].r ; f2 = xcx[43].r ; f4 = xcx[43].i ; xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ; xcx[43].r = f2+f1 ; xcx[43].i = f4+f3 ; f1 = xcx[55].r * csp[6].r - xcx[55].i * csp[6].i ; /* twiddles */ f3 = xcx[55].r * csp[6].i + xcx[55].i * csp[6].r ; f2 = xcx[51].r ; f4 = xcx[51].i ; xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ; xcx[51].r = f2+f1 ; xcx[51].i = f4+f3 ; f1 = xcx[63].r * csp[6].r - xcx[63].i * csp[6].i ; /* twiddles */ f3 = xcx[63].r * csp[6].i + xcx[63].i * csp[6].r ; f2 = xcx[59].r ; f4 = xcx[59].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[59].r = f2+f1 ; xcx[59].i = f4+f3 ; f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[24].r ; f3 = xcx[24].i ; /* cos=1 sin=0 */ f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[40].r ; f3 = xcx[40].i ; /* cos=1 sin=0 */ f2 = xcx[32].r ; f4 = xcx[32].i ; xcx[40].r = f2-f1 ; xcx[40].i = f4-f3 ; xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ; f1 = xcx[56].r ; f3 = xcx[56].i ; /* cos=1 sin=0 */ f2 = xcx[48].r ; f4 = xcx[48].i ; xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ; xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ; f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */ f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */ f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = xcx[41].r * csp[8].r - xcx[41].i * csp[8].i ; /* twiddles */ f3 = xcx[41].r * csp[8].i + xcx[41].i * csp[8].r ; f2 = xcx[33].r ; f4 = xcx[33].i ; xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ; xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ; f1 = xcx[57].r * csp[8].r - xcx[57].i * csp[8].i ; /* twiddles */ f3 = xcx[57].r * csp[8].i + xcx[57].i * csp[8].r ; f2 = xcx[49].r ; f4 = xcx[49].i ; xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ; xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ; f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */ f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */ f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ; f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = xcx[42].r * csp[9].r - xcx[42].i * csp[9].i ; /* twiddles */ f3 = xcx[42].r * csp[9].i + xcx[42].i * csp[9].r ; f2 = xcx[34].r ; f4 = xcx[34].i ; xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ; xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ; f1 = xcx[58].r * csp[9].r - xcx[58].i * csp[9].i ; /* twiddles */ f3 = xcx[58].r * csp[9].i + xcx[58].i * csp[9].r ; f2 = xcx[50].r ; f4 = xcx[50].i ; xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ; xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ; f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */ f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */ f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ; f2 = xcx[19].r ; f4 = xcx[19].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; f1 = xcx[43].r * csp[10].r - xcx[43].i * csp[10].i ; /* twiddles */ f3 = xcx[43].r * csp[10].i + xcx[43].i * csp[10].r ; f2 = xcx[35].r ; f4 = xcx[35].i ; xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ; xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ; f1 = xcx[59].r * csp[10].r - xcx[59].i * csp[10].i ; /* twiddles */ f3 = xcx[59].r * csp[10].i + xcx[59].i * csp[10].r ; f2 = xcx[51].r ; f4 = xcx[51].i ; xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ; xcx[51].r = f2+f1 ; xcx[51].i = f4+f3 ; f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[12].r * csp[11].i ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[28].r * csp[11].i ; f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = - xcx[44].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[44].r * csp[11].i ; f2 = xcx[36].r ; f4 = xcx[36].i ; xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ; xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ; f1 = - xcx[60].i * csp[11].i ; /* cos=0 twiddles */ f3 = xcx[60].r * csp[11].i ; f2 = xcx[52].r ; f4 = xcx[52].i ; xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ; xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ; f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */ f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */ f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ; f2 = xcx[21].r ; f4 = xcx[21].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; f1 = xcx[45].r * csp[12].r - xcx[45].i * csp[12].i ; /* twiddles */ f3 = xcx[45].r * csp[12].i + xcx[45].i * csp[12].r ; f2 = xcx[37].r ; f4 = xcx[37].i ; xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ; xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ; f1 = xcx[61].r * csp[12].r - xcx[61].i * csp[12].i ; /* twiddles */ f3 = xcx[61].r * csp[12].i + xcx[61].i * csp[12].r ; f2 = xcx[53].r ; f4 = xcx[53].i ; xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ; xcx[53].r = f2+f1 ; xcx[53].i = f4+f3 ; f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */ f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */ f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ; f2 = xcx[22].r ; f4 = xcx[22].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; f1 = xcx[46].r * csp[13].r - xcx[46].i * csp[13].i ; /* twiddles */ f3 = xcx[46].r * csp[13].i + xcx[46].i * csp[13].r ; f2 = xcx[38].r ; f4 = xcx[38].i ; xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ; xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ; f1 = xcx[62].r * csp[13].r - xcx[62].i * csp[13].i ; /* twiddles */ f3 = xcx[62].r * csp[13].i + xcx[62].i * csp[13].r ; f2 = xcx[54].r ; f4 = xcx[54].i ; xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ; xcx[54].r = f2+f1 ; xcx[54].i = f4+f3 ; f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */ f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */ f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ; f2 = xcx[23].r ; f4 = xcx[23].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ; f1 = xcx[47].r * csp[14].r - xcx[47].i * csp[14].i ; /* twiddles */ f3 = xcx[47].r * csp[14].i + xcx[47].i * csp[14].r ; f2 = xcx[39].r ; f4 = xcx[39].i ; xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ; xcx[39].r = f2+f1 ; xcx[39].i = f4+f3 ; f1 = xcx[63].r * csp[14].r - xcx[63].i * csp[14].i ; /* twiddles */ f3 = xcx[63].r * csp[14].i + xcx[63].i * csp[14].r ; f2 = xcx[55].r ; f4 = xcx[55].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[55].r = f2+f1 ; xcx[55].i = f4+f3 ; f1 = xcx[16].r ; f3 = xcx[16].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[48].r ; f3 = xcx[48].i ; /* cos=1 sin=0 */ f2 = xcx[32].r ; f4 = xcx[32].i ; xcx[48].r = f2-f1 ; xcx[48].i = f4-f3 ; xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ; f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */ f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[49].r * csp[16].r - xcx[49].i * csp[16].i ; /* twiddles */ f3 = xcx[49].r * csp[16].i + xcx[49].i * csp[16].r ; f2 = xcx[33].r ; f4 = xcx[33].i ; xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ; xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ; f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */ f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[50].r * csp[17].r - xcx[50].i * csp[17].i ; /* twiddles */ f3 = xcx[50].r * csp[17].i + xcx[50].i * csp[17].r ; f2 = xcx[34].r ; f4 = xcx[34].i ; xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ; xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ; f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */ f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[51].r * csp[18].r - xcx[51].i * csp[18].i ; /* twiddles */ f3 = xcx[51].r * csp[18].i + xcx[51].i * csp[18].r ; f2 = xcx[35].r ; f4 = xcx[35].i ; xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ; xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ; f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */ f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[52].r * csp[19].r - xcx[52].i * csp[19].i ; /* twiddles */ f3 = xcx[52].r * csp[19].i + xcx[52].i * csp[19].r ; f2 = xcx[36].r ; f4 = xcx[36].i ; xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ; xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ; f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */ f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[53].r * csp[20].r - xcx[53].i * csp[20].i ; /* twiddles */ f3 = xcx[53].r * csp[20].i + xcx[53].i * csp[20].r ; f2 = xcx[37].r ; f4 = xcx[37].i ; xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ; xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ; f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */ f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[54].r * csp[21].r - xcx[54].i * csp[21].i ; /* twiddles */ f3 = xcx[54].r * csp[21].i + xcx[54].i * csp[21].r ; f2 = xcx[38].r ; f4 = xcx[38].i ; xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ; xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ; f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */ f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; f1 = xcx[55].r * csp[22].r - xcx[55].i * csp[22].i ; /* twiddles */ f3 = xcx[55].r * csp[22].i + xcx[55].i * csp[22].r ; f2 = xcx[39].r ; f4 = xcx[39].i ; xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ; xcx[39].r = f2+f1 ; xcx[39].i = f4+f3 ; f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */ f3 = xcx[24].r * csp[23].i ; f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = - xcx[56].i * csp[23].i ; /* cos=0 twiddles */ f3 = xcx[56].r * csp[23].i ; f2 = xcx[40].r ; f4 = xcx[40].i ; xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ; xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ; f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */ f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[57].r * csp[24].r - xcx[57].i * csp[24].i ; /* twiddles */ f3 = xcx[57].r * csp[24].i + xcx[57].i * csp[24].r ; f2 = xcx[41].r ; f4 = xcx[41].i ; xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ; xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ; f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */ f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[58].r * csp[25].r - xcx[58].i * csp[25].i ; /* twiddles */ f3 = xcx[58].r * csp[25].i + xcx[58].i * csp[25].r ; f2 = xcx[42].r ; f4 = xcx[42].i ; xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ; xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ; f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */ f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[59].r * csp[26].r - xcx[59].i * csp[26].i ; /* twiddles */ f3 = xcx[59].r * csp[26].i + xcx[59].i * csp[26].r ; f2 = xcx[43].r ; f4 = xcx[43].i ; xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ; xcx[43].r = f2+f1 ; xcx[43].i = f4+f3 ; f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */ f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ; f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[60].r * csp[27].r - xcx[60].i * csp[27].i ; /* twiddles */ f3 = xcx[60].r * csp[27].i + xcx[60].i * csp[27].r ; f2 = xcx[44].r ; f4 = xcx[44].i ; xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ; xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ; f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */ f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = xcx[61].r * csp[28].r - xcx[61].i * csp[28].i ; /* twiddles */ f3 = xcx[61].r * csp[28].i + xcx[61].i * csp[28].r ; f2 = xcx[45].r ; f4 = xcx[45].i ; xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ; xcx[45].r = f2+f1 ; xcx[45].i = f4+f3 ; f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */ f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ; f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[62].r * csp[29].r - xcx[62].i * csp[29].i ; /* twiddles */ f3 = xcx[62].r * csp[29].i + xcx[62].i * csp[29].r ; f2 = xcx[46].r ; f4 = xcx[46].i ; xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ; xcx[46].r = f2+f1 ; xcx[46].i = f4+f3 ; f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */ f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ; f2 = xcx[15].r ; f4 = xcx[15].i ; xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ; f1 = xcx[63].r * csp[30].r - xcx[63].i * csp[30].i ; /* twiddles */ f3 = xcx[63].r * csp[30].i + xcx[63].i * csp[30].r ; f2 = xcx[47].r ; f4 = xcx[47].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[47].r = f2+f1 ; xcx[47].i = f4+f3 ; f1 = xcx[32].r ; f3 = xcx[32].i ; /* cos=1 sin=0 */ f2 = xcx[0].r ; f4 = xcx[0].i ; xcx[32].r = f2-f1 ; xcx[32].i = f4-f3 ; xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; f1 = xcx[33].r * csp[32].r - xcx[33].i * csp[32].i ; /* twiddles */ f3 = xcx[33].r * csp[32].i + xcx[33].i * csp[32].r ; f2 = xcx[1].r ; f4 = xcx[1].i ; xcx[33].r = f2-f1 ; xcx[33].i = f4-f3 ; xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; f1 = xcx[34].r * csp[33].r - xcx[34].i * csp[33].i ; /* twiddles */ f3 = xcx[34].r * csp[33].i + xcx[34].i * csp[33].r ; f2 = xcx[2].r ; f4 = xcx[2].i ; xcx[34].r = f2-f1 ; xcx[34].i = f4-f3 ; xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; f1 = xcx[35].r * csp[34].r - xcx[35].i * csp[34].i ; /* twiddles */ f3 = xcx[35].r * csp[34].i + xcx[35].i * csp[34].r ; f2 = xcx[3].r ; f4 = xcx[3].i ; xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ; xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; f1 = xcx[36].r * csp[35].r - xcx[36].i * csp[35].i ; /* twiddles */ f3 = xcx[36].r * csp[35].i + xcx[36].i * csp[35].r ; f2 = xcx[4].r ; f4 = xcx[4].i ; xcx[36].r = f2-f1 ; xcx[36].i = f4-f3 ; xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; f1 = xcx[37].r * csp[36].r - xcx[37].i * csp[36].i ; /* twiddles */ f3 = xcx[37].r * csp[36].i + xcx[37].i * csp[36].r ; f2 = xcx[5].r ; f4 = xcx[5].i ; xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ; xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; f1 = xcx[38].r * csp[37].r - xcx[38].i * csp[37].i ; /* twiddles */ f3 = xcx[38].r * csp[37].i + xcx[38].i * csp[37].r ; f2 = xcx[6].r ; f4 = xcx[6].i ; xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ; xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; f1 = xcx[39].r * csp[38].r - xcx[39].i * csp[38].i ; /* twiddles */ f3 = xcx[39].r * csp[38].i + xcx[39].i * csp[38].r ; f2 = xcx[7].r ; f4 = xcx[7].i ; xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ; xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; f1 = xcx[40].r * csp[39].r - xcx[40].i * csp[39].i ; /* twiddles */ f3 = xcx[40].r * csp[39].i + xcx[40].i * csp[39].r ; f2 = xcx[8].r ; f4 = xcx[8].i ; xcx[40].r = f2-f1 ; xcx[40].i = f4-f3 ; xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; f1 = xcx[41].r * csp[40].r - xcx[41].i * csp[40].i ; /* twiddles */ f3 = xcx[41].r * csp[40].i + xcx[41].i * csp[40].r ; f2 = xcx[9].r ; f4 = xcx[9].i ; xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ; xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; f1 = xcx[42].r * csp[41].r - xcx[42].i * csp[41].i ; /* twiddles */ f3 = xcx[42].r * csp[41].i + xcx[42].i * csp[41].r ; f2 = xcx[10].r ; f4 = xcx[10].i ; xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ; xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; f1 = xcx[43].r * csp[42].r - xcx[43].i * csp[42].i ; /* twiddles */ f3 = xcx[43].r * csp[42].i + xcx[43].i * csp[42].r ; f2 = xcx[11].r ; f4 = xcx[11].i ; xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ; xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; f1 = xcx[44].r * csp[43].r - xcx[44].i * csp[43].i ; /* twiddles */ f3 = xcx[44].r * csp[43].i + xcx[44].i * csp[43].r ; f2 = xcx[12].r ; f4 = xcx[12].i ; xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ; xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; f1 = xcx[45].r * csp[44].r - xcx[45].i * csp[44].i ; /* twiddles */ f3 = xcx[45].r * csp[44].i + xcx[45].i * csp[44].r ; f2 = xcx[13].r ; f4 = xcx[13].i ; xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ; xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; f1 = xcx[46].r * csp[45].r - xcx[46].i * csp[45].i ; /* twiddles */ f3 = xcx[46].r * csp[45].i + xcx[46].i * csp[45].r ; f2 = xcx[14].r ; f4 = xcx[14].i ; xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ; xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; f1 = xcx[47].r * csp[46].r - xcx[47].i * csp[46].i ; /* twiddles */ f3 = xcx[47].r * csp[46].i + xcx[47].i * csp[46].r ; f2 = xcx[15].r ; f4 = xcx[15].i ; xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ; xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ; f1 = - xcx[48].i * csp[47].i ; /* cos=0 twiddles */ f3 = xcx[48].r * csp[47].i ; f2 = xcx[16].r ; f4 = xcx[16].i ; xcx[48].r = f2-f1 ; xcx[48].i = f4-f3 ; xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; f1 = xcx[49].r * csp[48].r - xcx[49].i * csp[48].i ; /* twiddles */ f3 = xcx[49].r * csp[48].i + xcx[49].i * csp[48].r ; f2 = xcx[17].r ; f4 = xcx[17].i ; xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ; xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; f1 = xcx[50].r * csp[49].r - xcx[50].i * csp[49].i ; /* twiddles */ f3 = xcx[50].r * csp[49].i + xcx[50].i * csp[49].r ; f2 = xcx[18].r ; f4 = xcx[18].i ; xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ; xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; f1 = xcx[51].r * csp[50].r - xcx[51].i * csp[50].i ; /* twiddles */ f3 = xcx[51].r * csp[50].i + xcx[51].i * csp[50].r ; f2 = xcx[19].r ; f4 = xcx[19].i ; xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ; xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; f1 = xcx[52].r * csp[51].r - xcx[52].i * csp[51].i ; /* twiddles */ f3 = xcx[52].r * csp[51].i + xcx[52].i * csp[51].r ; f2 = xcx[20].r ; f4 = xcx[20].i ; xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ; xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; f1 = xcx[53].r * csp[52].r - xcx[53].i * csp[52].i ; /* twiddles */ f3 = xcx[53].r * csp[52].i + xcx[53].i * csp[52].r ; f2 = xcx[21].r ; f4 = xcx[21].i ; xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ; xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; f1 = xcx[54].r * csp[53].r - xcx[54].i * csp[53].i ; /* twiddles */ f3 = xcx[54].r * csp[53].i + xcx[54].i * csp[53].r ; f2 = xcx[22].r ; f4 = xcx[22].i ; xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ; xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; f1 = xcx[55].r * csp[54].r - xcx[55].i * csp[54].i ; /* twiddles */ f3 = xcx[55].r * csp[54].i + xcx[55].i * csp[54].r ; f2 = xcx[23].r ; f4 = xcx[23].i ; xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ; xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ; f1 = xcx[56].r * csp[55].r - xcx[56].i * csp[55].i ; /* twiddles */ f3 = xcx[56].r * csp[55].i + xcx[56].i * csp[55].r ; f2 = xcx[24].r ; f4 = xcx[24].i ; xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ; xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; f1 = xcx[57].r * csp[56].r - xcx[57].i * csp[56].i ; /* twiddles */ f3 = xcx[57].r * csp[56].i + xcx[57].i * csp[56].r ; f2 = xcx[25].r ; f4 = xcx[25].i ; xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ; xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; f1 = xcx[58].r * csp[57].r - xcx[58].i * csp[57].i ; /* twiddles */ f3 = xcx[58].r * csp[57].i + xcx[58].i * csp[57].r ; f2 = xcx[26].r ; f4 = xcx[26].i ; xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ; xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; f1 = xcx[59].r * csp[58].r - xcx[59].i * csp[58].i ; /* twiddles */ f3 = xcx[59].r * csp[58].i + xcx[59].i * csp[58].r ; f2 = xcx[27].r ; f4 = xcx[27].i ; xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ; xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ; f1 = xcx[60].r * csp[59].r - xcx[60].i * csp[59].i ; /* twiddles */ f3 = xcx[60].r * csp[59].i + xcx[60].i * csp[59].r ; f2 = xcx[28].r ; f4 = xcx[28].i ; xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ; xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; f1 = xcx[61].r * csp[60].r - xcx[61].i * csp[60].i ; /* twiddles */ f3 = xcx[61].r * csp[60].i + xcx[61].i * csp[60].r ; f2 = xcx[29].r ; f4 = xcx[29].i ; xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ; xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ; f1 = xcx[62].r * csp[61].r - xcx[62].i * csp[61].i ; /* twiddles */ f3 = xcx[62].r * csp[61].i + xcx[62].i * csp[61].r ; f2 = xcx[30].r ; f4 = xcx[30].i ; xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ; xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ; f1 = xcx[63].r * csp[62].r - xcx[63].i * csp[62].i ; /* twiddles */ f3 = xcx[63].r * csp[62].i + xcx[63].i * csp[62].r ; f2 = xcx[31].r ; f4 = xcx[31].i ; xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ; xcx[31].r = f2+f1 ; xcx[31].i = f4+f3 ; return ; } /*================================================================*/ /*------------------------------------------------------------------ fft_4dec: do a decimation by 4. At most RMAX levels of recursion are allowed. --------------------------------------------------------------------*/ static void fft_4dec_OMP( int mode , int idim , complex *xc ) { DECLARE_ithr ; int rec , mold ; int N=idim , M=idim/4 , M2=2*M , M3=3*M ; complex *cs , *aa , *bb , *cc , *dd ; int k , i ; float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 , acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ; /* ININFO_message("fft_4dec_OMP(%d)",idim) ; */ /*-- initialize cosine/sine table and memory space --*/ rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */ if( rec >= RMAX ){ fprintf(stderr,"\n*** fft_4dec_OMP: too many recursions with idim=%d\n",idim); EXIT(1); } CFO_assign(cs,M3) ; /* assign workspace of length M3 to cs pointer */ CFO_assign(aa,M) ; /* at this recursion level for this thread (etc) */ CFO_assign(bb,M) ; CFO_assign(cc,M) ; CFO_assign(dd,M) ; mold = CFO_rmold4[rec][ithr] ; /* what was done before at this recursion level */ { if( M != mold ){ double th = (2.0*PI/N) ; cs[0].r = 1.0 ; cs[0].i = 0.0 ; for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } CFO_rmold4[rec][ithr] = M ; }} /*-- load subarrays, and FFT each one --*/ for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; } CFO_recX[ithr]++ ; /* bump recursion level */ switch( N ){ default: fprintf(stderr,"\n*** Illegal call to fft_4dec: N=%d\n",N); EXIT(1) ; case 128: fft32(mode,aa); fft32(mode,bb); /* no recursion */ fft32(mode,cc); fft32(mode,dd); break; case 256: fft64(mode,aa); fft64(mode,bb); /* no recursion */ fft64(mode,cc); fft64(mode,dd); break; case 512: fft128(mode,aa); fft128(mode,bb); /* recurse once */ fft128(mode,cc); fft128(mode,dd); break; case 1024: fft256(mode,aa); fft256(mode,bb); /* recurse once */ fft256(mode,cc); fft256(mode,dd); break; #if 1 case 2048: fft512(mode,aa); fft512(mode,bb); /* recurse twice */ fft512(mode,cc); fft512(mode,dd); break; case 4096: fft1024(mode,aa); fft1024(mode,bb); /* recurse twice */ fft1024(mode,cc); fft1024(mode,dd); break; case 8192: fft2048(mode,aa); fft2048(mode,bb); /* recurse 3 times */ fft2048(mode,cc); fft2048(mode,dd); break; case 16384: fft4096(mode,aa); fft4096(mode,bb); /* recurse 3 times */ fft4096(mode,cc); fft4096(mode,dd); break; case 32768: fft8192(mode,aa); fft8192(mode,bb); /* recurse 4 times */ fft8192(mode,cc); fft8192(mode,dd); break; case 65536: fft16384(mode,aa); fft16384(mode,bb); /* recurse 4 times */ fft16384(mode,cc); fft16384(mode,dd); break; #endif } CFO_recX[ithr]-- ; /* reduce recursion level */ /*-- recombination --*/ if( mode > 0 ){ for( k=0 ; k < M ; k++ ){ t1 = bb[k].r; t2 = bb[k].i; tr = cs[k].r; ti = cs[k].i; bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ t1 = cc[ k].r; t2 = cc[ k].i; tr = cs[2*k].r; ti = cs[2*k].i; ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ t1 = dd[ k].r; t2 = dd[ k].i; tr = cs[3*k].r; ti = cs[3*k].i; ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */ aar = aa[k].r; aai = aa[k].i; /* a[k] */ acpr = aar + ccr ; acmr = aar - ccr ; bdpr = bbr + ddr ; bdmr = bbr - ddr ; acpi = aai + cci ; acmi = aai - cci ; bdpi = bbi + ddi ; bdmi = bbi - ddi ; xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ; xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ; } } else { for( k=0 ; k < M ; k++ ){ t1 = bb[k].r; t2 = bb[k].i; tr = cs[k].r; ti = cs[k].i; bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ t1 = cc[ k].r; t2 = cc[ k].i; tr = cs[2*k].r; ti = cs[2*k].i; ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ t1 = dd[ k].r; t2 = dd[ k].i; tr = cs[3*k].r; ti = cs[3*k].i; ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */ aar = aa[k].r; aai = aa[k].i; /* a[k] */ acpr = aar + ccr ; acmr = aar - ccr ; bdpr = bbr + ddr ; bdmr = bbr - ddr ; acpi = aai + cci ; acmi = aai - cci ; bdpi = bbi + ddi ; bdmi = bbi - ddi ; xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ; xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ; } } return ; } /*======================================================================= The following radix-3 and radix-5 routines are by RWCox -- Aug 1999. They are not inside unrolled code. =========================================================================*/ /*------------------------------------------------------------------ fft_3dec: do a decimation-by-3, plus a power-of-2. At most RMAX levels of recursion are allowed, which means that at most 3**RMAX can be a factor of idim. --------------------------------------------------------------------*/ #undef CC3 #undef SS3 #define CC3 (-0.5) /* cos(2*Pi/3) */ #define SS3 (0.8660254038) /* sin(2*Pi/3) */ static void fft_3dec_OMP( int mode , int idim , complex *xc ) { DECLARE_ithr ; int rec , mold ; int N=idim , M=idim/3 , M2=2*M ; complex *cs=NULL , *aa , *bb , *cc ; register int k , i ; register float aar,aai, tr,ti, bbr,bbi, ccr,cci, t1,t2,t4,t5,t6,t8 ; /* ININFO_message("fft_3dec_OMP(%d)",idim) ; */ /*-- initialize cosine/sine table and memory space --*/ rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */ if( rec >= RMAX ){ fprintf(stderr,"\n*** fft_3dec_OMP: too many recursions at idim=%d\n",idim); EXIT(1); } CFO_assign(cs,M2) ; CFO_assign(aa,M) ; CFO_assign(bb,M) ; CFO_assign(cc,M) ; mold = CFO_rmold3[rec][ithr] ; /* what was done before at this recursion level */ { if( M != mold ){ double th = (2.0*PI/N) ; cs[0].r = 1.0 ; cs[0].i = 0.0 ; for( k=1 ; k < M2 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } CFO_rmold3[rec][ithr] = M ; }} /*-- load subarrays, and FFT each one --*/ for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; } CFO_recX[ithr]++ ; switch( M ){ case 1: break ; case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ; break; case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ; break; case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ; break; case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; break; case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; break; case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; break; case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; break; case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; break; case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; break; case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; break; #if 1 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; break; case 4096: fft4096(mode,aa) ; fft4096(mode,bb) ; fft4096(mode,cc) ; break; case 8192: fft8192(mode,aa) ; fft8192(mode,bb) ; fft8192(mode,cc) ; break; case 16384: fft16384(mode,aa); fft16384(mode,bb); fft16384(mode,cc); break; #endif default: csfft_cox_OMP(mode,M,aa); csfft_cox_OMP(mode,M,bb); csfft_cox_OMP(mode,M,cc); break ; } CFO_recX[ithr]-- ; /*-- recombination --*/ if( mode > 0 ){ for( k=0 ; k < M ; k++ ){ t1 = bb[k].r; t2 = bb[k].i; tr = cs[k].r; ti = cs[k].i; bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ t1 = cc[ k].r; t2 = cc[ k].i; tr = cs[2*k].r; ti = cs[2*k].i; ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ t4 = bbr+ccr ; t1 = t4*CC3 ; t8 = bbi+cci ; t6 = t8*CC3 ; t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3; aar = aa[k].r; aai = aa[k].i; /* a[k] */ xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ; xc[k+M ].r = (aar+t1)-t2 ; xc[k+M ].i = (aai+t6)+t5 ; xc[k+M2].r = (aar+t1)+t2 ; xc[k+M2].i = (aai+t6)-t5 ; } } else { for( k=0 ; k < M ; k++ ){ t1 = bb[k].r; t2 = bb[k].i; tr = cs[k].r; ti = cs[k].i; bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ t1 = cc[ k].r; t2 = cc[ k].i; tr = cs[2*k].r; ti = cs[2*k].i; ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ t4 = bbr+ccr ; t1 = t4*CC3 ; t8 = bbi+cci ; t6 = t8*CC3 ; t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3; aar = aa[k].r; aai = aa[k].i; /* a[k] */ xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ; xc[k+M ].r = (aar+t1)+t2 ; xc[k+M ].i = (aai+t6)-t5 ; xc[k+M2].r = (aar+t1)-t2 ; xc[k+M2].i = (aai+t6)+t5 ; } } return ; } /*------------------------------------------------------------------ fft_5dec: do a decimation-by-5, plus a power-of-2. At most RMAX levels of recursion are allowed, which means that at most 5**RMAX can be a factor of idim. --------------------------------------------------------------------*/ #undef COS72 #undef SIN72 #define COS72 0.30901699 /* cos(72 deg) */ #define SIN72 0.95105652 /* sin(72 deg) */ static void fft_5dec_OMP( int mode , int idim , complex *xc ) { DECLARE_ithr ; int rec , mold ; int N=idim , M=idim/5 , M2=2*M , M3=3*M , M4=4*M ; complex *cs, *aa, *bb, *cc, *dd, *ee ; register int k , i ; register float aar,aai,bbr,bbi,ccr,cci,ddr,ddi,eer,eei ; register float akp,akm,bkp,bkm , ajp,ajm,bjp,bjm ; register float ak,bk,aj,bj ; float c72 , s72 , c2,s2 ; int ss ; /* ININFO_message("fft_5dec_OMP(%d)",idim) ; */ /*-- initialize cosine/sine table and memory space --*/ rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */ if( rec >= RMAX ){ fprintf(stderr,"\n*** fft_5dec_OMP: too many recursions at idim=%d\n",idim); EXIT(1); } CFO_assign(cs,M4) ; CFO_assign(aa,M) ; CFO_assign(bb,M) ; CFO_assign(cc,M) ; CFO_assign(dd,M) ; CFO_assign(ee,M) ; mold = CFO_rmold5[rec][ithr] ; /* what was done before at this recursion level */ { if( M != mold ){ /* new M? */ double th = (2.0*PI/N) ; cs[0].r = 1.0 ; cs[0].i = 0.0 ; for( k=1 ; k < M4 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } CFO_rmold5[rec][ithr] = M ; }} /*-- load subarrays, and FFT each one --*/ for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; ee[k]=xc[i++]; } CFO_recX[ithr]++ ; switch( M ){ case 1: break ; case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ; fft2 (mode,dd) ; fft2 (mode,ee) ; break; case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ; fft4 (mode,dd) ; fft4 (mode,ee) ; break; case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ; fft8 (mode,dd) ; fft8 (mode,ee) ; break; case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; fft16 (mode,dd) ; fft16 (mode,ee) ; break; case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; fft32 (mode,dd) ; fft32 (mode,ee) ; break; case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; fft64 (mode,dd) ; fft64 (mode,ee) ; break; case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; fft128(mode,dd) ; fft128(mode,ee) ; break; case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; fft256(mode,dd) ; fft256(mode,ee) ; break; case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; fft512(mode,dd) ; fft512(mode,ee) ; break; case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; fft1024(mode,dd) ; fft1024(mode,ee) ; break; #if 1 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; fft2048(mode,dd) ; fft2048(mode,ee) ; break; case 4096: fft4096(mode,aa) ; fft4096(mode,bb) ; fft4096(mode,cc) ; fft4096(mode,dd) ; fft4096(mode,ee) ; break; case 8192: fft8192(mode,aa) ; fft8192(mode,bb) ; fft8192(mode,cc) ; fft8192(mode,dd) ; fft8192(mode,ee) ; break; case 16384: fft16384(mode,aa) ; fft16384(mode,bb) ; fft16384(mode,cc) ; fft16384(mode,dd) ; fft16384(mode,ee) ; break; #endif default: csfft_cox_OMP(mode,M,aa) ; csfft_cox_OMP(mode,M,bb) ; csfft_cox_OMP(mode,M,cc) ; csfft_cox_OMP(mode,M,dd) ; csfft_cox_OMP(mode,M,ee) ; break ; } CFO_recX[ithr]-- ; /*-- recombination --*/ if( mode > 0 ){ s72 = SIN72 ; ss = 1 ; } else { s72 = -SIN72 ; ss = -1 ; } c72 = COS72 ; c2 = c72 * c72 - s72 * s72 ; s2 = 2.0 * c72 * s72 ; for( k=0 ; k < M ; k++ ){ ak = bb[k].r ; bk = bb[k].i ; aj = cs[k].r ; bj = cs[k].i * ss ; bbr = aj*ak - bj*bk ; bbi = aj*bk + bj*ak ; /* b[k]*exp(ss*2*Pi*k/N) */ ak = cc[ k].r ; bk = cc[ k].i ; aj = cs[2*k].r ; bj = cs[2*k].i * ss ; ccr = aj*ak - bj*bk ; cci = aj*bk + bj*ak ; /* c[k]*exp(ss*4*Pi*k/N) */ ak = dd[ k].r ; bk = dd[ k].i ; aj = cs[3*k].r ; bj = cs[3*k].i * ss ; ddr = aj*ak - bj*bk ; ddi = aj*bk + bj*ak ; /* d[k]*exp(ss*6*Pi*k/N) */ ak = ee[ k].r ; bk = ee[ k].i ; aj = cs[4*k].r ; bj = cs[4*k].i * ss ; eer = aj*ak - bj*bk ; eei = aj*bk + bj*ak ; /* e[k]*exp(ss*8*Pi*k/N) */ aar = aa[k].r ; aai = aa[k].i ; /* a[k] */ /* the code below is (heavily) adapted from fftn.c */ akp = bbr + eer ; akm = bbr - eer ; bkp = bbi + eei ; bkm = bbi - eei ; ajp = ccr + ddr ; ajm = ccr - ddr ; bjp = cci + ddi ; bjm = cci - ddi ; xc[k].r = aar + akp + ajp ; xc[k].i = aai + bkp + bjp ; ak = akp * c72 + ajp * c2 + aar ; bk = bkp * c72 + bjp * c2 + aai ; aj = akm * s72 + ajm * s2 ; bj = bkm * s72 + bjm * s2 ; xc[k+M ].r = ak - bj ; xc[k+M ].i = bk + aj ; xc[k+M4].r = ak + bj ; xc[k+M4].i = bk - aj ; ak = akp * c2 + ajp * c72 + aar ; bk = bkp * c2 + bjp * c72 + aai ; aj = akm * s2 - ajm * s72 ; bj = bkm * s2 - bjm * s72 ; xc[k+M2].r = ak - bj ; xc[k+M2].i = bk + aj ; xc[k+M3].r = ak + bj ; xc[k+M3].i = bk - aj ; } return ; }