Doxygen Source Code Documentation
csfft.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | RMAX 3 |
#define | PI (3.141592653589793238462643) |
#define | fft1024(m, x) fft_4dec(m,1024,x) |
#define | fft2048(m, x) fft_4dec(m,2048,x) |
#define | fft2(m, x) |
#define | USE_FFT4_MACRO |
#define | fft4(m, x) |
#define | SCLINV |
#define | I00 |
#define | I0BOT 1 |
#define | N 128 |
#define | M 32 |
#define | M2 64 |
#define | M3 96 |
#define | N 256 |
#define | M 64 |
#define | M2 128 |
#define | M3 192 |
#define | N 512 |
#define | M 128 |
#define | M2 256 |
#define | M3 384 |
#define | CC3 (-0.5) |
#define | SS3 (0.8660254038) |
#define | COS72 0.30901699 |
#define | SIN72 0.95105652 |
#define | N35 ((RMAX+1)*(RMAX+1)) |
Functions | |
void | csfft_use_fftw (int uf) |
void | csfft_cox (int mode, int idim, complex *xc) |
int | csfft_nextup (int idim) |
int | csfft_nextup_one35 (int idim) |
void | csfft_scale_inverse (int scl) |
void | fft_3dec (int, int, complex *) |
void | fft_5dec (int, int, complex *) |
void | fft8 (int mode, complex *xc) |
void | fft16 (int mode, complex *xc) |
void | fft32 (int mode, complex *xc) |
void | fft64 (int mode, complex *xc) |
void | fft128 (int mode, complex *xc) |
void | fft256 (int mode, complex *xc) |
void | fft512 (int mode, complex *xc) |
void | fft_4dec (int, int, complex *) |
void | csfft_trigconsts (int idim) |
void | csfft_many (int mode, int idim, int nvec, complex *xc) |
int | csfft_nextup_even (int idim) |
Variables | |
int | use_fftw = 1 |
int | sclinv = 0 |
complex * | csplus = NULL |
complex * | csminus = NULL |
int | nold = -666 |
Define Documentation
|
Definition at line 1864 of file csfft.c. Referenced by fft_3dec(). |
|
Definition at line 1997 of file csfft.c. Referenced by fft_5dec(). |
|
Definition at line 90 of file csfft.c. Referenced by csfft_cox(), fft_3dec(), and fft_5dec(). |
|
Value: 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) Definition at line 95 of file csfft.c. Referenced by csfft_cox(), fft_3dec(), and fft_5dec(). |
|
Definition at line 91 of file csfft.c. Referenced by csfft_cox(), fft_3dec(), and fft_5dec(). |
|
Value: 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) Definition at line 105 of file csfft.c. Referenced by csfft_cox(), fft_3dec(), and fft_5dec(). |
|
|
|
|
|
|
|
|
|
Definition at line 1614 of file csfft.c. Referenced by fft128(), fft256(), fft512(), fft_3dec(), fft_4dec(), fft_5dec(), intrcall(), iocname(), mkpower(), mktmpn(), newentry(), oneof_stg(), opconv_fudge(), out_addr(), out_call(), p1_addr(), putcx1(), putio(), putop(), set_externs(), wr_equiv_init(), and yyparse(). |
|
|
|
|
|
Definition at line 1615 of file csfft.c. Referenced by fft128(), fft256(), fft512(), fft_3dec(), fft_4dec(), and fft_5dec(). |
|
|
|
|
|
Definition at line 1616 of file csfft.c. Referenced by fft128(), fft256(), fft512(), fft_4dec(), and fft_5dec(). |
|
|
|
|
|
Definition at line 1613 of file csfft.c. Referenced by fft128(), fft256(), fft512(), fft_3dec(), fft_4dec(), and fft_5dec(). |
|
Definition at line 2163 of file csfft.c. Referenced by csfft_nextup(). |
|
Definition at line 68 of file csfft.c. Referenced by csfft_trigconsts(), fft128(), fft256(), fft512(), fft64(), fft_3dec(), fft_4dec(), and fft_5dec(). |
|
Definition at line 65 of file csfft.c. Referenced by csfft_nextup(), fft_3dec(), fft_4dec(), and fft_5dec(). |
|
Value: if( sclinv && 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 ; } } Definition at line 206 of file csfft.c. Referenced by csfft_cox(). |
|
Definition at line 1998 of file csfft.c. Referenced by fft_5dec(). |
|
Definition at line 1865 of file csfft.c. Referenced by fft_3dec(). |
|
|
Function Documentation
|
Definition at line 211 of file csfft.c. References AFMALL, close(), csfft_trigconsts(), fd, fft1024, fft128(), fft16(), fft2, fft2048, fft256(), fft32(), fft4, fft512(), fft64(), fft8(), fft_3dec(), fft_4dec(), fft_5dec(), FFTW_BACKWARD, fftw_create_plan(), fftw_destroy_plan(), FFTW_ESTIMATE, FFTW_FORWARD, fftw_import_wisdom_from_string(), FFTW_IN_PLACE, fftw_one(), FFTW_SUCCESS, FFTW_USE_WISDOM, free, getenv(), complex::i, i1, i2, nold, complex::r, read(), SCLINV, THD_filesize(), and xc.
00212 { 00213 register unsigned int m, n, i0, i1, i2, i3, k; 00214 register complex *r0, *r1, *csp; 00215 register float co, si, f0, f1, f2, f3, f4; 00216 00217 static int rec=0 ; /* recursion level */ 00218 00219 /*-- 20 Oct 2000: maybe use FFTW instead --*/ 00220 00221 #ifdef USE_FFTW 00222 if( use_fftw ){ 00223 static int first=1 , oldmode=0, oldidim=0 ; 00224 static fftw_plan fpl ; 00225 if( first ){ 00226 char * env = getenv("AFNI_FFTW_WISDOM") ; 00227 int gotit=0 ; 00228 first = 0 ; 00229 if( env != NULL ){ 00230 int nw = THD_filesize(env) ; 00231 if( nw > 0 ){ 00232 int fd = open( env , O_RDONLY ) ; 00233 if( fd >= 0 ){ 00234 char * w = AFMALL( char, nw+4 ) ; 00235 int ii = read( fd , w , nw ) ; 00236 close(fd) ; 00237 if( ii > 0 ){ 00238 w[nw] = '\0' ; 00239 nw = (int) fftw_import_wisdom_from_string(w) ; 00240 gotit=(nw == FFTW_SUCCESS); 00241 } /* else fprintf(stderr,"csfft_cox: read fails\n"); */ 00242 free(w) ; 00243 } /* else fprintf(stderr,"csfft_cox: open fails\n"); */ 00244 } /* else fprintf(stderr,"csfft_cox: THD_filesize(%s) fails\n",env); */ 00245 } /* else fprintf(stderr,"csfft_cox: getenv fails\n"); */ 00246 /* if( gotit ) 00247 fprintf(stderr,"csfft_cox imported FFTW wisdom from file %s\n",env); 00248 else 00249 fprintf(stderr,"csfft_cox failed to import FFTW wisdom\n"); */ 00250 } 00251 00252 if( idim != oldidim || mode != oldmode ){ 00253 if( oldidim ) fftw_destroy_plan(fpl) ; 00254 fpl = fftw_create_plan( idim , 00255 (mode < 0) ? FFTW_FORWARD : FFTW_BACKWARD , 00256 FFTW_ESTIMATE|FFTW_USE_WISDOM|FFTW_IN_PLACE ); 00257 oldidim = idim ; oldmode = mode ; 00258 } 00259 00260 fftw_one( fpl , (fftw_complex *)xc , NULL ) ; 00261 SCLINV ; return ; 00262 } 00263 #endif /* USE_FFTW */ 00264 00265 /*-- November 1998: maybe use the unrolled FFT routines --*/ 00266 00267 #ifndef DONT_UNROLL_FFTS 00268 switch( idim ){ /* none of these */ 00269 case 1: return; /* routines will */ 00270 case 2: fft2 (mode,xc); SCLINV; return; /* call csfft_cox */ 00271 case 4: fft4 (mode,xc); SCLINV; return; /* so don't need */ 00272 case 8: fft8 (mode,xc); SCLINV; return; /* to do rec++/-- */ 00273 case 16: fft16 (mode,xc); SCLINV; return; 00274 case 32: fft32 (mode,xc); SCLINV; return; 00275 case 64: fft64 (mode,xc); SCLINV; return; 00276 case 128: fft128 (mode,xc); SCLINV; return; 00277 case 256: fft256 (mode,xc); SCLINV; return; 00278 case 512: fft512 (mode,xc); SCLINV; return; 00279 case 1024: fft1024(mode,xc); SCLINV; return; 00280 case 2048: fft2048(mode,xc); SCLINV; return; 00281 00282 case 4096: fft_4dec(mode, 4096,xc); SCLINV; return; /* 4dec -> 1024 */ 00283 case 8192: fft_4dec(mode, 8192,xc); SCLINV; return; /* 4dec -> 2048 */ 00284 case 16384: fft_4dec(mode,16384,xc); SCLINV; return; /* 4dec -> 4096 */ 00285 case 32768: fft_4dec(mode,32768,xc); SCLINV; return; /* 4dec -> 8192 */ 00286 } 00287 #endif /* end of unrollificationizing */ 00288 00289 /*-- Aug 1999: deal with factors of 3 or 5 [might call csfft_cox] --*/ 00290 00291 if( idim%3 == 0 ){rec++; fft_3dec(mode,idim,xc); rec--; SCLINV; return;} 00292 if( idim%5 == 0 ){rec++; fft_5dec(mode,idim,xc); rec--; SCLINV; return;} 00293 00294 /*-- below here: the old general power of 2 routine --*/ 00295 00296 /**-- perhaps initialize --**/ 00297 00298 if( nold != idim ) csfft_trigconsts( idim ) ; 00299 00300 /** Main loop starts here **/ 00301 00302 n = idim; 00303 i2 = idim >> 1; 00304 i1 = 0; 00305 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00306 00307 /*-- swap data --*/ 00308 00309 for (i0=0; i0 < n; i0 ++) { 00310 if ( i1 > i0 ) { 00311 r0 = xc + i0; r1 = xc + i1; 00312 f1 = r0->r; f2 = r0->i; 00313 r0->r = r1->r; r0->i = r1->i; 00314 r1->r = f1; r1->i = f2; 00315 } 00316 m = i2; 00317 while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; } 00318 i1 += m; 00319 } 00320 00321 /*-- compute compute compute --*/ 00322 00323 m = 1; k = 0; 00324 while (n > m) { 00325 i3 = m << 1; 00326 for (i0=0; i0 < m; i0 ++) { 00327 co = (csp + k)->r; si = (csp + k)->i; 00328 for (i1=i0; i1 < n; i1 += i3) { 00329 r0 = xc + i1; r1 = r0 + m; 00330 f1 = r1->r * co; f2 = r1->i * si; 00331 f3 = r1->r * si; f4 = r1->i * co; 00332 f1 -= f2; f3 += f4; 00333 f2 = r0->r; f4 = r0->i; 00334 r1->r = f2 - f1; r1->i = f4 - f3; 00335 r0->r = f2 + f1; r0->i = f4 + f3; 00336 } 00337 k++; 00338 } 00339 m = i3; 00340 } 00341 00342 SCLINV ; return ; 00343 } |
|
Main loop starts here * Definition at line 355 of file csfft.c. References csfft_cox(), csfft_trigconsts(), fft_3dec(), fft_5dec(), complex::i, i1, i2, nold, complex::r, and xc. Referenced by main().
00356 { 00357 register unsigned int m, n, i0, i1, i2, i3, k , iv ; 00358 register complex *r0, *r1, *csp , *xcx; 00359 register float co, si, f0, f1, f2, f3, f4; 00360 00361 if( nvec == 1 ){ csfft_cox( mode , idim , xc ) ; return ; } 00362 00363 if( idim % 3 == 0 ){ /* Aug 1999 */ 00364 for( m=0 ; m < nvec ; m++ ) 00365 fft_3dec( mode , idim , xc + m*idim ) ; 00366 return ; 00367 } else if( idim % 5 == 0 ){ 00368 for( m=0 ; m < nvec ; m++ ) 00369 fft_5dec( mode , idim , xc + m*idim ) ; 00370 return ; 00371 } 00372 00373 /** perhaps initialize **/ 00374 00375 if( nold != idim ) csfft_trigconsts( idim ) ; 00376 00377 n = idim; 00378 i2 = idim >> 1; 00379 i1 = 0; 00380 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00381 00382 for (i0=0; i0 < n; i0 ++) { 00383 if ( i1 > i0 ) { 00384 for( iv=0,xcx=xc ; iv < nvec ; iv++,xcx+=n ){ 00385 r0 = xcx + i0; r1 = xcx + i1; 00386 f1 = r0->r; f2 = r0->i; 00387 r0->r = r1->r; r0->i = r1->i; 00388 r1->r = f1; r1->i = f2; 00389 } 00390 } 00391 m = i2; 00392 while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; } 00393 i1 += m; 00394 } 00395 00396 #define I00 00397 #ifdef I00 00398 # define I0BOT 1 00399 #else 00400 # define I0BOT 0 00401 #endif 00402 00403 m = 1; 00404 k = 0; 00405 while (n > m) { 00406 i3 = m << 1; 00407 00408 #ifdef I00 00409 /* handle i0=0 case [co=1,si=0] in special code */ 00410 00411 for (i1=0; i1 < n; i1 += i3) { 00412 for( iv=0,r0=xc+i1 ; iv < nvec ; iv++,r0+=n ){ 00413 r1 = r0 + m; 00414 f1 = r1->r ; f3 = r1->i ; 00415 f2 = r0->r ; f4 = r0->i ; 00416 r1->r = f2 - f1 ; r1->i = f4 - f3 ; 00417 r0->r = f2 + f1 ; r0->i = f4 + f3 ; 00418 } 00419 } 00420 k++; 00421 #endif 00422 00423 for (i0=I0BOT; i0 < m; i0 ++) { 00424 co = (csp + k)->r; si = (csp + k)->i; 00425 for (i1=i0; i1 < n; i1 += i3) { 00426 for( iv=0,r0=xc+i1 ; iv < nvec ; iv++,r0+=n ){ 00427 r1 = r0 + m; 00428 #if 1 00429 f1 = r1->r * co - r1->i * si ; 00430 f3 = r1->r * si + r1->i * co ; 00431 #else 00432 f1 = r1->r * co ; f2 = r1->i * si ; 00433 f3 = r1->r * si ; f4 = r1->i * co ; 00434 f1 -= f2 ; f3 += f4 ; 00435 #endif 00436 f2 = r0->r ; f4 = r0->i ; 00437 r1->r = f2 - f1 ; r1->i = f4 - f3 ; 00438 r0->r = f2 + f1 ; r0->i = f4 + f3 ; 00439 } 00440 } 00441 k++; 00442 } 00443 m = i3; 00444 } 00445 return ; 00446 } |
|
Definition at line 2165 of file csfft.c. References free, i, malloc, N35, p, q, r, RMAX, and tt.
02166 { 02167 static int * tf = NULL , * dn = NULL ; 02168 int ibase , ii ; 02169 02170 /*-- build table of allowable powers of 3 and 5 [tf], 02171 the powers of 2 just less than them [dn], 02172 and their ratios tf/dn [rat]. 02173 Then sort tf and dn to be increasing in rat. --*/ 02174 02175 if( tf == NULL ){ 02176 int p , q , tt,ff , i=0 , j ; float * rat ; float r ; 02177 02178 tf = (int *) malloc(sizeof(int) *N35) ; 02179 dn = (int *) malloc(sizeof(int) *N35) ; 02180 rat = (float *) malloc(sizeof(float)*N35) ; 02181 02182 /* create tables */ 02183 02184 for( p=0,tt=1 ; p <= RMAX ; p++,tt*=3 ){ /* tt = 3^p */ 02185 for( q=0,ff=1 ; q <= RMAX ; q++,ff*=5,i++ ){ /* ff = 5^q */ 02186 tf[i] = tt * ff ; 02187 02188 j=2; while( j < tf[i] ){ j*=2; } /* j = power of 2 just >= tf */ 02189 dn[i] = j/2 ; /* dn = power of 2 just < tf */ 02190 02191 rat[i] = (float)(tf[i]) / (float)(dn[i]) ; 02192 } 02193 } 02194 02195 /* sort on rat (crude, but it works) */ 02196 02197 do{ 02198 for( i=1,p=0 ; i < N35 ; i++ ){ 02199 if( rat[i-1] > rat[i] ){ 02200 r = rat[i-1] ; rat[i-1] = rat[i] ; rat[i] = r ; 02201 q = tf [i-1] ; tf [i-1] = tf [i] ; tf [i] = q ; 02202 q = dn [i-1] ; dn [i-1] = dn [i] ; dn [i] = q ; 02203 p++ ; 02204 } 02205 } 02206 } while( p > 0 ) ; 02207 02208 free(rat) ; 02209 } 02210 02211 /*-- loop on powers of 2 (ibase); 02212 we can do FFTs of size = tf*ibase/dn (note 1 < tf/dn < 2); 02213 sinc tf/dn is sorted, we're scanning in increasing sizes --*/ 02214 02215 ibase = 1 ; 02216 while(1){ 02217 if( idim <= ibase ) return ibase ; 02218 02219 for( ii=0 ; ii < N35 ; ii++ ) 02220 if( dn[ii] <= ibase && idim <= tf[ii]*ibase/dn[ii] ) 02221 return tf[ii]*ibase/dn[ii] ; 02222 02223 ibase *= 2 ; 02224 } 02225 } |
|
Definition at line 2247 of file csfft.c. References csfft_nextup(). Referenced by apply_xshear(), apply_yshear(), and apply_zshear().
02248 { 02249 int jj = idim ; 02250 do{ 02251 jj = csfft_nextup(jj) ; 02252 if( jj%2 == 1 ) jj++ ; 02253 else return jj ; 02254 } while(1) ; 02255 return 0 ; /* cannot be reached */ 02256 } |
|
Definition at line 2232 of file csfft.c. References csfft_nextup().
02233 { 02234 int jj = idim ; 02235 do{ 02236 jj = csfft_nextup(jj) ; 02237 if( jj%9 == 0 || jj%25 == 0 || jj%2 == 1 ) jj++ ; 02238 else return jj ; 02239 } while(1) ; 02240 return 0 ; /* cannot be reached */ 02241 } |
|
Definition at line 74 of file csfft.c. References sclinv.
00074 { sclinv = scl; return; } |
|
Definition at line 139 of file csfft.c. References EXIT, free, complex::i, i1, i2, malloc, nold, PI, and complex::r.
00140 { 00141 register unsigned int m, n, i0, i1, i2, i3, k; 00142 register complex *r0, *csp; 00143 register float co, si, f0, f1, f2, f3, f4; 00144 double al; 00145 00146 if( idim == nold ) return ; 00147 00148 if( idim > nold ){ 00149 if( csplus != NULL ){ free(csplus) ; free(csminus) ; } 00150 csplus = (complex *) malloc( sizeof(complex) * idim ) ; 00151 csminus = (complex *) malloc( sizeof(complex) * idim ) ; 00152 if( csplus == NULL || csminus == NULL ){ 00153 fprintf(stderr,"\n*** csfft cannot malloc space! ***\n"); EXIT(1) ; 00154 } 00155 } 00156 nold = n = idim ; 00157 00158 f1 = 1.0 ; /* csplus init */ 00159 m = 1; k = 0; 00160 while (n > m) { 00161 i3 = m << 1; f2 = m; al = f1*PI/f2; 00162 co = cos(al); si = sin(al); 00163 (csplus + k)->r = 1.; (csplus + k)->i = 0.; 00164 for (i0=0; i0 < m; i0++) { 00165 k++; 00166 csp = csplus + k; r0 = csp - 1; 00167 csp->r = r0->r * co - r0->i * si; 00168 csp->i = r0->i * co + r0->r * si; 00169 } 00170 m = i3; 00171 } 00172 00173 f1 = -1.0 ; /* csminus init */ 00174 m = 1; k = 0; 00175 while (n > m) { 00176 i3 = m << 1; f2 = m; al = f1*PI/f2; 00177 co = cos(al); si = sin(al); 00178 (csminus + k)->r = 1.; (csminus + k)->i = 0.; 00179 for (i0=0; i0 < m; i0++) { 00180 k++; 00181 csp = csminus + k; r0 = csp - 1; 00182 csp->r = r0->r * co - r0->i * si; 00183 csp->i = r0->i * co + r0->r * si; 00184 } 00185 m = i3; 00186 } 00187 return ; 00188 } |
|
Definition at line 28 of file csfft.c. References use_fftw. Referenced by main().
00028 { use_fftw = uf; return; } |
|
Definition at line 1417 of file csfft.c. References fft32(), complex::i, i, M, M2, M3, malloc, N, PI, complex::r, and xc. Referenced by csfft_cox(), fft512(), fft_3dec(), and fft_5dec().
01418 { 01419 static complex * cs=NULL , * aa , * bb , * cc , * dd ; 01420 register int k , i ; 01421 register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 , 01422 acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ; 01423 01424 /*-- initialize cosine/sine table and memory space --*/ 01425 01426 if( cs == NULL ){ 01427 double th = (2.0*PI/N) ; 01428 cs = (complex *) malloc(sizeof(complex)*M3) ; 01429 aa = (complex *) malloc(sizeof(complex)*M ) ; 01430 bb = (complex *) malloc(sizeof(complex)*M ) ; 01431 cc = (complex *) malloc(sizeof(complex)*M ) ; 01432 dd = (complex *) malloc(sizeof(complex)*M ) ; 01433 01434 cs[0].r = 1.0 ; cs[0].i = 0.0 ; 01435 for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } 01436 } 01437 01438 /*-- load subarrays, and FFT each one --*/ 01439 01440 for( i=k=0; k < M; k++ ){ 01441 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; 01442 } 01443 01444 fft32(mode,aa); fft32(mode,bb); fft32(mode,cc); fft32(mode,dd); 01445 01446 /*-- recombination --*/ 01447 01448 if( mode > 0 ){ 01449 for( k=0 ; k < M ; k++ ){ 01450 t1 = bb[k].r; t2 = bb[k].i; 01451 tr = cs[k].r; ti = cs[k].i; 01452 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ 01453 01454 t1 = cc[ k].r; t2 = cc[ k].i; 01455 tr = cs[2*k].r; ti = cs[2*k].i; 01456 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ 01457 01458 t1 = dd[ k].r; t2 = dd[ k].i; 01459 tr = cs[3*k].r; ti = cs[3*k].i; 01460 ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */ 01461 01462 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01463 01464 acpr = aar + ccr ; acmr = aar - ccr ; 01465 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01466 acpi = aai + cci ; acmi = aai - cci ; 01467 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01468 01469 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01470 xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ; 01471 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01472 xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ; 01473 } 01474 } else { 01475 for( k=0 ; k < M ; k++ ){ 01476 t1 = bb[k].r; t2 = bb[k].i; 01477 tr = cs[k].r; ti = cs[k].i; 01478 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ 01479 01480 t1 = cc[ k].r; t2 = cc[ k].i; 01481 tr = cs[2*k].r; ti = cs[2*k].i; 01482 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ 01483 01484 t1 = dd[ k].r; t2 = dd[ k].i; 01485 tr = cs[3*k].r; ti = cs[3*k].i; 01486 ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */ 01487 01488 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01489 01490 acpr = aar + ccr ; acmr = aar - ccr ; 01491 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01492 acpi = aai + cci ; acmi = aai - cci ; 01493 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01494 01495 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01496 xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ; 01497 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01498 xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ; 01499 } 01500 } 01501 01502 return ; 01503 } |
|
butterflying part * Definition at line 600 of file csfft.c. References csfft_trigconsts(), complex::i, nold, complex::r, and xc. Referenced by csfft_cox(), fft_3dec(), and fft_5dec().
00601 { 00602 register complex * csp , * xcx=xc; 00603 register float f1,f2,f3,f4 ; 00604 00605 /** perhaps initialize **/ 00606 00607 if( nold != 16 ) csfft_trigconsts( 16 ) ; 00608 00609 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00610 00611 /** data swapping part **/ 00612 00613 f1 = xcx[1].r ; f2 = xcx[1].i ; 00614 xcx[1].r = xcx[8].r ; xcx[1].i = xcx[8].i ; 00615 xcx[8].r = f1 ; xcx[8].i = f2 ; 00616 00617 f1 = xcx[2].r ; f2 = xcx[2].i ; 00618 xcx[2].r = xcx[4].r ; xcx[2].i = xcx[4].i ; 00619 xcx[4].r = f1 ; xcx[4].i = f2 ; 00620 00621 f1 = xcx[3].r ; f2 = xcx[3].i ; 00622 xcx[3].r = xcx[12].r ; xcx[3].i = xcx[12].i ; 00623 xcx[12].r = f1 ; xcx[12].i = f2 ; 00624 00625 f1 = xcx[5].r ; f2 = xcx[5].i ; 00626 xcx[5].r = xcx[10].r ; xcx[5].i = xcx[10].i ; 00627 xcx[10].r = f1 ; xcx[10].i = f2 ; 00628 00629 f1 = xcx[7].r ; f2 = xcx[7].i ; 00630 xcx[7].r = xcx[14].r ; xcx[7].i = xcx[14].i ; 00631 xcx[14].r = f1 ; xcx[14].i = f2 ; 00632 00633 f1 = xcx[11].r ; f2 = xcx[11].i ; 00634 xcx[11].r = xcx[13].r ; xcx[11].i = xcx[13].i ; 00635 xcx[13].r = f1 ; xcx[13].i = f2 ; 00636 00637 /** butterflying part **/ 00638 00639 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ 00640 f2 = xcx[0].r ; f4 = xcx[0].i ; 00641 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; 00642 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00643 00644 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ 00645 f2 = xcx[2].r ; f4 = xcx[2].i ; 00646 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 00647 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00648 00649 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ 00650 f2 = xcx[4].r ; f4 = xcx[4].i ; 00651 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 00652 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00653 00654 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ 00655 f2 = xcx[6].r ; f4 = xcx[6].i ; 00656 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00657 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 00658 00659 f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */ 00660 f2 = xcx[8].r ; f4 = xcx[8].i ; 00661 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; 00662 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 00663 00664 f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */ 00665 f2 = xcx[10].r ; f4 = xcx[10].i ; 00666 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 00667 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; 00668 00669 f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */ 00670 f2 = xcx[12].r ; f4 = xcx[12].i ; 00671 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 00672 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; 00673 00674 f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */ 00675 f2 = xcx[14].r ; f4 = xcx[14].i ; 00676 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 00677 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; 00678 00679 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ 00680 f2 = xcx[0].r ; f4 = xcx[0].i ; 00681 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; 00682 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00683 00684 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ 00685 f2 = xcx[4].r ; f4 = xcx[4].i ; 00686 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 00687 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00688 00689 f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */ 00690 f2 = xcx[8].r ; f4 = xcx[8].i ; 00691 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; 00692 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 00693 00694 f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */ 00695 f2 = xcx[12].r ; f4 = xcx[12].i ; 00696 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 00697 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; 00698 00699 f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */ 00700 f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ; 00701 f2 = xcx[1].r ; f4 = xcx[1].i ; 00702 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 00703 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 00704 00705 f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */ 00706 f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ; 00707 f2 = xcx[5].r ; f4 = xcx[5].i ; 00708 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00709 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 00710 00711 f1 = xcx[11].r * csp[2].r - xcx[11].i * csp[2].i ; /* twiddles */ 00712 f3 = xcx[11].r * csp[2].i + xcx[11].i * csp[2].r ; 00713 f2 = xcx[9].r ; f4 = xcx[9].i ; 00714 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 00715 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; 00716 00717 f1 = xcx[15].r * csp[2].r - xcx[15].i * csp[2].i ; /* twiddles */ 00718 f3 = xcx[15].r * csp[2].i + xcx[15].i * csp[2].r ; 00719 f2 = xcx[13].r ; f4 = xcx[13].i ; 00720 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 00721 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; 00722 00723 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ 00724 f2 = xcx[0].r ; f4 = xcx[0].i ; 00725 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; 00726 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00727 00728 f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */ 00729 f2 = xcx[8].r ; f4 = xcx[8].i ; 00730 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; 00731 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 00732 00733 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ 00734 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; 00735 f2 = xcx[1].r ; f4 = xcx[1].i ; 00736 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 00737 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 00738 00739 f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */ 00740 f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ; 00741 f2 = xcx[9].r ; f4 = xcx[9].i ; 00742 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 00743 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; 00744 00745 f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */ 00746 f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ; 00747 f2 = xcx[2].r ; f4 = xcx[2].i ; 00748 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 00749 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00750 00751 f1 = xcx[14].r * csp[5].r - xcx[14].i * csp[5].i ; /* twiddles */ 00752 f3 = xcx[14].r * csp[5].i + xcx[14].i * csp[5].r ; 00753 f2 = xcx[10].r ; f4 = xcx[10].i ; 00754 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 00755 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; 00756 00757 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ 00758 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; 00759 f2 = xcx[3].r ; f4 = xcx[3].i ; 00760 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00761 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 00762 00763 f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */ 00764 f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ; 00765 f2 = xcx[11].r ; f4 = xcx[11].i ; 00766 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 00767 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; 00768 00769 f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */ 00770 f2 = xcx[0].r ; f4 = xcx[0].i ; 00771 xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ; 00772 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00773 00774 f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */ 00775 f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ; 00776 f2 = xcx[1].r ; f4 = xcx[1].i ; 00777 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; 00778 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 00779 00780 f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */ 00781 f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ; 00782 f2 = xcx[2].r ; f4 = xcx[2].i ; 00783 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; 00784 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00785 00786 f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */ 00787 f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ; 00788 f2 = xcx[3].r ; f4 = xcx[3].i ; 00789 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 00790 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 00791 00792 f1 = xcx[12].r * csp[11].r - xcx[12].i * csp[11].i ; /* twiddles */ 00793 f3 = xcx[12].r * csp[11].i + xcx[12].i * csp[11].r ; 00794 f2 = xcx[4].r ; f4 = xcx[4].i ; 00795 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; 00796 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00797 00798 f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */ 00799 f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ; 00800 f2 = xcx[5].r ; f4 = xcx[5].i ; 00801 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 00802 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 00803 00804 f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */ 00805 f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ; 00806 f2 = xcx[6].r ; f4 = xcx[6].i ; 00807 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 00808 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 00809 00810 f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */ 00811 f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ; 00812 f2 = xcx[7].r ; f4 = xcx[7].i ; 00813 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 00814 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; 00815 00816 return ; 00817 } |
|
Definition at line 1517 of file csfft.c. References fft64(), complex::i, i, M, M2, M3, malloc, N, PI, complex::r, and xc. Referenced by csfft_cox(), fft_3dec(), fft_4dec(), and fft_5dec().
01518 { 01519 static complex * cs=NULL , * aa , * bb , * cc , * dd ; 01520 register int k , i ; 01521 register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 , 01522 acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ; 01523 01524 /*-- initialize cosine/sine table and memory space --*/ 01525 01526 if( cs == NULL ){ 01527 double th = (2.0*PI/N) ; 01528 cs = (complex *) malloc(sizeof(complex)*M3) ; 01529 aa = (complex *) malloc(sizeof(complex)*M ) ; 01530 bb = (complex *) malloc(sizeof(complex)*M ) ; 01531 cc = (complex *) malloc(sizeof(complex)*M ) ; 01532 dd = (complex *) malloc(sizeof(complex)*M ) ; 01533 01534 cs[0].r = 1.0 ; cs[0].i = 0.0 ; 01535 for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } 01536 } 01537 01538 /*-- load subarrays, and FFT each one --*/ 01539 01540 for( i=k=0; k < M; k++ ){ 01541 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; 01542 } 01543 01544 fft64(mode,aa); fft64(mode,bb); fft64(mode,cc); fft64(mode,dd); 01545 01546 /*-- recombination --*/ 01547 01548 if( mode > 0 ){ 01549 for( k=0 ; k < M ; k++ ){ 01550 t1 = bb[k].r; t2 = bb[k].i; 01551 tr = cs[k].r; ti = cs[k].i; 01552 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ 01553 01554 t1 = cc[ k].r; t2 = cc[ k].i; 01555 tr = cs[2*k].r; ti = cs[2*k].i; 01556 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ 01557 01558 t1 = dd[ k].r; t2 = dd[ k].i; 01559 tr = cs[3*k].r; ti = cs[3*k].i; 01560 ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */ 01561 01562 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01563 01564 acpr = aar + ccr ; acmr = aar - ccr ; 01565 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01566 acpi = aai + cci ; acmi = aai - cci ; 01567 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01568 01569 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01570 xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ; 01571 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01572 xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ; 01573 } 01574 } else { 01575 for( k=0 ; k < M ; k++ ){ 01576 t1 = bb[k].r; t2 = bb[k].i; 01577 tr = cs[k].r; ti = cs[k].i; 01578 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ 01579 01580 t1 = cc[ k].r; t2 = cc[ k].i; 01581 tr = cs[2*k].r; ti = cs[2*k].i; 01582 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ 01583 01584 t1 = dd[ k].r; t2 = dd[ k].i; 01585 tr = cs[3*k].r; ti = cs[3*k].i; 01586 ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */ 01587 01588 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01589 01590 acpr = aar + ccr ; acmr = aar - ccr ; 01591 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01592 acpi = aai + cci ; acmi = aai - cci ; 01593 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01594 01595 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01596 xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ; 01597 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01598 xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ; 01599 } 01600 } 01601 01602 return ; 01603 } |
|
butterflying part * Definition at line 822 of file csfft.c. References csfft_trigconsts(), complex::i, nold, complex::r, and xc. Referenced by csfft_cox(), fft128(), fft64(), fft_3dec(), and fft_5dec().
00823 { 00824 register complex * csp , * xcx=xc; 00825 register float f1,f2,f3,f4 ; 00826 00827 /** perhaps initialize **/ 00828 00829 if( nold != 32 ) csfft_trigconsts( 32 ) ; 00830 00831 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00832 00833 /** data swapping part **/ 00834 00835 f1 = xcx[1].r ; f2 = xcx[1].i ; 00836 xcx[1].r = xcx[16].r ; xcx[1].i = xcx[16].i ; 00837 xcx[16].r = f1 ; xcx[16].i = f2 ; 00838 00839 f1 = xcx[2].r ; f2 = xcx[2].i ; 00840 xcx[2].r = xcx[8].r ; xcx[2].i = xcx[8].i ; 00841 xcx[8].r = f1 ; xcx[8].i = f2 ; 00842 00843 f1 = xcx[3].r ; f2 = xcx[3].i ; 00844 xcx[3].r = xcx[24].r ; xcx[3].i = xcx[24].i ; 00845 xcx[24].r = f1 ; xcx[24].i = f2 ; 00846 00847 f1 = xcx[5].r ; f2 = xcx[5].i ; 00848 xcx[5].r = xcx[20].r ; xcx[5].i = xcx[20].i ; 00849 xcx[20].r = f1 ; xcx[20].i = f2 ; 00850 00851 f1 = xcx[6].r ; f2 = xcx[6].i ; 00852 xcx[6].r = xcx[12].r ; xcx[6].i = xcx[12].i ; 00853 xcx[12].r = f1 ; xcx[12].i = f2 ; 00854 00855 f1 = xcx[7].r ; f2 = xcx[7].i ; 00856 xcx[7].r = xcx[28].r ; xcx[7].i = xcx[28].i ; 00857 xcx[28].r = f1 ; xcx[28].i = f2 ; 00858 00859 f1 = xcx[9].r ; f2 = xcx[9].i ; 00860 xcx[9].r = xcx[18].r ; xcx[9].i = xcx[18].i ; 00861 xcx[18].r = f1 ; xcx[18].i = f2 ; 00862 00863 f1 = xcx[11].r ; f2 = xcx[11].i ; 00864 xcx[11].r = xcx[26].r ; xcx[11].i = xcx[26].i ; 00865 xcx[26].r = f1 ; xcx[26].i = f2 ; 00866 00867 f1 = xcx[13].r ; f2 = xcx[13].i ; 00868 xcx[13].r = xcx[22].r ; xcx[13].i = xcx[22].i ; 00869 xcx[22].r = f1 ; xcx[22].i = f2 ; 00870 00871 f1 = xcx[15].r ; f2 = xcx[15].i ; 00872 xcx[15].r = xcx[30].r ; xcx[15].i = xcx[30].i ; 00873 xcx[30].r = f1 ; xcx[30].i = f2 ; 00874 00875 f1 = xcx[19].r ; f2 = xcx[19].i ; 00876 xcx[19].r = xcx[25].r ; xcx[19].i = xcx[25].i ; 00877 xcx[25].r = f1 ; xcx[25].i = f2 ; 00878 00879 f1 = xcx[23].r ; f2 = xcx[23].i ; 00880 xcx[23].r = xcx[29].r ; xcx[23].i = xcx[29].i ; 00881 xcx[29].r = f1 ; xcx[29].i = f2 ; 00882 00883 /** butterflying part **/ 00884 00885 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ 00886 f2 = xcx[0].r ; f4 = xcx[0].i ; 00887 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; 00888 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00889 00890 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ 00891 f2 = xcx[2].r ; f4 = xcx[2].i ; 00892 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 00893 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00894 00895 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ 00896 f2 = xcx[4].r ; f4 = xcx[4].i ; 00897 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 00898 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00899 00900 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ 00901 f2 = xcx[6].r ; f4 = xcx[6].i ; 00902 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00903 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 00904 00905 f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */ 00906 f2 = xcx[8].r ; f4 = xcx[8].i ; 00907 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; 00908 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 00909 00910 f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */ 00911 f2 = xcx[10].r ; f4 = xcx[10].i ; 00912 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 00913 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; 00914 00915 f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */ 00916 f2 = xcx[12].r ; f4 = xcx[12].i ; 00917 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 00918 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; 00919 00920 f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */ 00921 f2 = xcx[14].r ; f4 = xcx[14].i ; 00922 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 00923 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; 00924 00925 f1 = xcx[17].r ; f3 = xcx[17].i ; /* cos=1 sin=0 */ 00926 f2 = xcx[16].r ; f4 = xcx[16].i ; 00927 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; 00928 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; 00929 00930 f1 = xcx[19].r ; f3 = xcx[19].i ; /* cos=1 sin=0 */ 00931 f2 = xcx[18].r ; f4 = xcx[18].i ; 00932 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; 00933 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; 00934 00935 f1 = xcx[21].r ; f3 = xcx[21].i ; /* cos=1 sin=0 */ 00936 f2 = xcx[20].r ; f4 = xcx[20].i ; 00937 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; 00938 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; 00939 00940 f1 = xcx[23].r ; f3 = xcx[23].i ; /* cos=1 sin=0 */ 00941 f2 = xcx[22].r ; f4 = xcx[22].i ; 00942 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; 00943 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; 00944 00945 f1 = xcx[25].r ; f3 = xcx[25].i ; /* cos=1 sin=0 */ 00946 f2 = xcx[24].r ; f4 = xcx[24].i ; 00947 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; 00948 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; 00949 00950 f1 = xcx[27].r ; f3 = xcx[27].i ; /* cos=1 sin=0 */ 00951 f2 = xcx[26].r ; f4 = xcx[26].i ; 00952 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; 00953 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; 00954 00955 f1 = xcx[29].r ; f3 = xcx[29].i ; /* cos=1 sin=0 */ 00956 f2 = xcx[28].r ; f4 = xcx[28].i ; 00957 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; 00958 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; 00959 00960 f1 = xcx[31].r ; f3 = xcx[31].i ; /* cos=1 sin=0 */ 00961 f2 = xcx[30].r ; f4 = xcx[30].i ; 00962 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; 00963 xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ; 00964 00965 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ 00966 f2 = xcx[0].r ; f4 = xcx[0].i ; 00967 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; 00968 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00969 00970 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ 00971 f2 = xcx[4].r ; f4 = xcx[4].i ; 00972 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 00973 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00974 00975 f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */ 00976 f2 = xcx[8].r ; f4 = xcx[8].i ; 00977 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; 00978 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 00979 00980 f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */ 00981 f2 = xcx[12].r ; f4 = xcx[12].i ; 00982 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 00983 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; 00984 00985 f1 = xcx[18].r ; f3 = xcx[18].i ; /* cos=1 sin=0 */ 00986 f2 = xcx[16].r ; f4 = xcx[16].i ; 00987 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; 00988 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; 00989 00990 f1 = xcx[22].r ; f3 = xcx[22].i ; /* cos=1 sin=0 */ 00991 f2 = xcx[20].r ; f4 = xcx[20].i ; 00992 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; 00993 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; 00994 00995 f1 = xcx[26].r ; f3 = xcx[26].i ; /* cos=1 sin=0 */ 00996 f2 = xcx[24].r ; f4 = xcx[24].i ; 00997 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; 00998 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; 00999 01000 f1 = xcx[30].r ; f3 = xcx[30].i ; /* cos=1 sin=0 */ 01001 f2 = xcx[28].r ; f4 = xcx[28].i ; 01002 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; 01003 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ; 01004 01005 f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */ 01006 f3 = xcx[3].r * csp[2].i ; 01007 f2 = xcx[1].r ; f4 = xcx[1].i ; 01008 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 01009 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 01010 01011 f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */ 01012 f3 = xcx[7].r * csp[2].i ; 01013 f2 = xcx[5].r ; f4 = xcx[5].i ; 01014 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 01015 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 01016 01017 f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */ 01018 f3 = xcx[11].r * csp[2].i ; 01019 f2 = xcx[9].r ; f4 = xcx[9].i ; 01020 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 01021 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; 01022 01023 f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */ 01024 f3 = xcx[15].r * csp[2].i ; 01025 f2 = xcx[13].r ; f4 = xcx[13].i ; 01026 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 01027 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; 01028 01029 f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */ 01030 f3 = xcx[19].r * csp[2].i ; 01031 f2 = xcx[17].r ; f4 = xcx[17].i ; 01032 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; 01033 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; 01034 01035 f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */ 01036 f3 = xcx[23].r * csp[2].i ; 01037 f2 = xcx[21].r ; f4 = xcx[21].i ; 01038 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; 01039 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; 01040 01041 f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */ 01042 f3 = xcx[27].r * csp[2].i ; 01043 f2 = xcx[25].r ; f4 = xcx[25].i ; 01044 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; 01045 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; 01046 01047 f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */ 01048 f3 = xcx[31].r * csp[2].i ; 01049 f2 = xcx[29].r ; f4 = xcx[29].i ; 01050 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; 01051 xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ; 01052 01053 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ 01054 f2 = xcx[0].r ; f4 = xcx[0].i ; 01055 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; 01056 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 01057 01058 f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */ 01059 f2 = xcx[8].r ; f4 = xcx[8].i ; 01060 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; 01061 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 01062 01063 f1 = xcx[20].r ; f3 = xcx[20].i ; /* cos=1 sin=0 */ 01064 f2 = xcx[16].r ; f4 = xcx[16].i ; 01065 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; 01066 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; 01067 01068 f1 = xcx[28].r ; f3 = xcx[28].i ; /* cos=1 sin=0 */ 01069 f2 = xcx[24].r ; f4 = xcx[24].i ; 01070 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; 01071 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ; 01072 01073 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ 01074 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; 01075 f2 = xcx[1].r ; f4 = xcx[1].i ; 01076 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 01077 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 01078 01079 f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */ 01080 f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ; 01081 f2 = xcx[9].r ; f4 = xcx[9].i ; 01082 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 01083 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; 01084 01085 f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */ 01086 f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ; 01087 f2 = xcx[17].r ; f4 = xcx[17].i ; 01088 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; 01089 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; 01090 01091 f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */ 01092 f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ; 01093 f2 = xcx[25].r ; f4 = xcx[25].i ; 01094 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; 01095 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ; 01096 01097 f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */ 01098 f3 = xcx[6].r * csp[5].i ; 01099 f2 = xcx[2].r ; f4 = xcx[2].i ; 01100 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 01101 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 01102 01103 f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */ 01104 f3 = xcx[14].r * csp[5].i ; 01105 f2 = xcx[10].r ; f4 = xcx[10].i ; 01106 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 01107 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; 01108 01109 f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */ 01110 f3 = xcx[22].r * csp[5].i ; 01111 f2 = xcx[18].r ; f4 = xcx[18].i ; 01112 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; 01113 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; 01114 01115 f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */ 01116 f3 = xcx[30].r * csp[5].i ; 01117 f2 = xcx[26].r ; f4 = xcx[26].i ; 01118 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; 01119 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ; 01120 01121 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ 01122 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; 01123 f2 = xcx[3].r ; f4 = xcx[3].i ; 01124 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 01125 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 01126 01127 f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */ 01128 f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ; 01129 f2 = xcx[11].r ; f4 = xcx[11].i ; 01130 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 01131 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; 01132 01133 f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */ 01134 f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ; 01135 f2 = xcx[19].r ; f4 = xcx[19].i ; 01136 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; 01137 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; 01138 01139 f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */ 01140 f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ; 01141 f2 = xcx[27].r ; f4 = xcx[27].i ; 01142 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; 01143 xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ; 01144 01145 f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */ 01146 f2 = xcx[0].r ; f4 = xcx[0].i ; 01147 xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ; 01148 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 01149 01150 f1 = xcx[24].r ; f3 = xcx[24].i ; /* cos=1 sin=0 */ 01151 f2 = xcx[16].r ; f4 = xcx[16].i ; 01152 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; 01153 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ; 01154 01155 f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */ 01156 f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ; 01157 f2 = xcx[1].r ; f4 = xcx[1].i ; 01158 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ; 01159 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 01160 01161 f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */ 01162 f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ; 01163 f2 = xcx[17].r ; f4 = xcx[17].i ; 01164 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; 01165 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ; 01166 01167 f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */ 01168 f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ; 01169 f2 = xcx[2].r ; f4 = xcx[2].i ; 01170 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ; 01171 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 01172 01173 f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */ 01174 f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ; 01175 f2 = xcx[18].r ; f4 = xcx[18].i ; 01176 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; 01177 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ; 01178 01179 f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */ 01180 f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ; 01181 f2 = xcx[3].r ; f4 = xcx[3].i ; 01182 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ; 01183 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 01184 01185 f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */ 01186 f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ; 01187 f2 = xcx[19].r ; f4 = xcx[19].i ; 01188 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; 01189 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ; 01190 01191 f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */ 01192 f3 = xcx[12].r * csp[11].i ; 01193 f2 = xcx[4].r ; f4 = xcx[4].i ; 01194 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ; 01195 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 01196 01197 f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */ 01198 f3 = xcx[28].r * csp[11].i ; 01199 f2 = xcx[20].r ; f4 = xcx[20].i ; 01200 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; 01201 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ; 01202 01203 f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */ 01204 f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ; 01205 f2 = xcx[5].r ; f4 = xcx[5].i ; 01206 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ; 01207 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 01208 01209 f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */ 01210 f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ; 01211 f2 = xcx[21].r ; f4 = xcx[21].i ; 01212 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; 01213 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ; 01214 01215 f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */ 01216 f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ; 01217 f2 = xcx[6].r ; f4 = xcx[6].i ; 01218 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ; 01219 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 01220 01221 f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */ 01222 f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ; 01223 f2 = xcx[22].r ; f4 = xcx[22].i ; 01224 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; 01225 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ; 01226 01227 f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */ 01228 f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ; 01229 f2 = xcx[7].r ; f4 = xcx[7].i ; 01230 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ; 01231 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; 01232 01233 f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */ 01234 f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ; 01235 f2 = xcx[23].r ; f4 = xcx[23].i ; 01236 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; 01237 xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ; 01238 01239 f1 = xcx[16].r ; f3 = xcx[16].i ; /* cos=1 sin=0 */ 01240 f2 = xcx[0].r ; f4 = xcx[0].i ; 01241 xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ; 01242 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 01243 01244 f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */ 01245 f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ; 01246 f2 = xcx[1].r ; f4 = xcx[1].i ; 01247 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ; 01248 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 01249 01250 f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */ 01251 f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ; 01252 f2 = xcx[2].r ; f4 = xcx[2].i ; 01253 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ; 01254 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 01255 01256 f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */ 01257 f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ; 01258 f2 = xcx[3].r ; f4 = xcx[3].i ; 01259 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ; 01260 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 01261 01262 f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */ 01263 f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ; 01264 f2 = xcx[4].r ; f4 = xcx[4].i ; 01265 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ; 01266 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 01267 01268 f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */ 01269 f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ; 01270 f2 = xcx[5].r ; f4 = xcx[5].i ; 01271 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ; 01272 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 01273 01274 f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */ 01275 f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ; 01276 f2 = xcx[6].r ; f4 = xcx[6].i ; 01277 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ; 01278 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 01279 01280 f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */ 01281 f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ; 01282 f2 = xcx[7].r ; f4 = xcx[7].i ; 01283 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ; 01284 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ; 01285 01286 f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */ 01287 f3 = xcx[24].r * csp[23].i ; 01288 f2 = xcx[8].r ; f4 = xcx[8].i ; 01289 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ; 01290 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ; 01291 01292 f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */ 01293 f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ; 01294 f2 = xcx[9].r ; f4 = xcx[9].i ; 01295 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ; 01296 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ; 01297 01298 f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */ 01299 f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ; 01300 f2 = xcx[10].r ; f4 = xcx[10].i ; 01301 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ; 01302 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ; 01303 01304 f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */ 01305 f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ; 01306 f2 = xcx[11].r ; f4 = xcx[11].i ; 01307 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ; 01308 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ; 01309 01310 f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */ 01311 f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ; 01312 f2 = xcx[12].r ; f4 = xcx[12].i ; 01313 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ; 01314 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ; 01315 01316 f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */ 01317 f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ; 01318 f2 = xcx[13].r ; f4 = xcx[13].i ; 01319 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ; 01320 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ; 01321 01322 f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */ 01323 f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ; 01324 f2 = xcx[14].r ; f4 = xcx[14].i ; 01325 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ; 01326 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ; 01327 01328 f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */ 01329 f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ; 01330 f2 = xcx[15].r ; f4 = xcx[15].i ; 01331 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ; 01332 xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ; 01333 01334 return ; 01335 } |
|
Definition at line 1617 of file csfft.c. References fft128(), complex::i, i, M, M2, M3, malloc, N, PI, complex::r, and xc. Referenced by csfft_cox(), fft_3dec(), fft_4dec(), and fft_5dec().
01618 { 01619 static complex * cs=NULL , * aa , * bb , * cc , * dd ; 01620 register int k , i ; 01621 register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 , 01622 acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ; 01623 01624 /*-- initialize cosine/sine table and memory space --*/ 01625 01626 if( cs == NULL ){ 01627 double th = (2.0*PI/N) ; 01628 cs = (complex *) malloc(sizeof(complex)*M3) ; 01629 aa = (complex *) malloc(sizeof(complex)*M ) ; 01630 bb = (complex *) malloc(sizeof(complex)*M ) ; 01631 cc = (complex *) malloc(sizeof(complex)*M ) ; 01632 dd = (complex *) malloc(sizeof(complex)*M ) ; 01633 01634 cs[0].r = 1.0 ; cs[0].i = 0.0 ; 01635 for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } 01636 } 01637 01638 /*-- load subarrays, and FFT each one --*/ 01639 01640 for( i=k=0; k < M; k++ ){ 01641 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; 01642 } 01643 01644 fft128(mode,aa); fft128(mode,bb); fft128(mode,cc); fft128(mode,dd); 01645 01646 /*-- recombination --*/ 01647 01648 if( mode > 0 ){ 01649 for( k=0 ; k < M ; k++ ){ 01650 t1 = bb[k].r; t2 = bb[k].i; 01651 tr = cs[k].r; ti = cs[k].i; 01652 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ 01653 01654 t1 = cc[ k].r; t2 = cc[ k].i; 01655 tr = cs[2*k].r; ti = cs[2*k].i; 01656 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ 01657 01658 t1 = dd[ k].r; t2 = dd[ k].i; 01659 tr = cs[3*k].r; ti = cs[3*k].i; 01660 ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */ 01661 01662 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01663 01664 acpr = aar + ccr ; acmr = aar - ccr ; 01665 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01666 acpi = aai + cci ; acmi = aai - cci ; 01667 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01668 01669 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01670 xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ; 01671 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01672 xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ; 01673 } 01674 } else { 01675 for( k=0 ; k < M ; k++ ){ 01676 t1 = bb[k].r; t2 = bb[k].i; 01677 tr = cs[k].r; ti = cs[k].i; 01678 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ 01679 01680 t1 = cc[ k].r; t2 = cc[ k].i; 01681 tr = cs[2*k].r; ti = cs[2*k].i; 01682 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ 01683 01684 t1 = dd[ k].r; t2 = dd[ k].i; 01685 tr = cs[3*k].r; ti = cs[3*k].i; 01686 ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */ 01687 01688 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01689 01690 acpr = aar + ccr ; acmr = aar - ccr ; 01691 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01692 acpi = aai + cci ; acmi = aai - cci ; 01693 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01694 01695 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01696 xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ; 01697 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01698 xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ; 01699 } 01700 } 01701 01702 return ; 01703 } |
|
butterflying part * Definition at line 1341 of file csfft.c. References fft32(), complex::i, i, malloc, PI, complex::r, and xc. Referenced by csfft_cox(), fft256(), fft_3dec(), and fft_5dec().
01342 { 01343 static complex * cs=NULL , * aa , * bb ; 01344 register int k , i ; 01345 register float akr,aki , bkr,bki , tr,ti , t1,t2 ; 01346 01347 /*-- initialize cosine/sine table --*/ 01348 01349 if( cs == NULL ){ 01350 double th = (PI/32.0) ; 01351 cs = (complex *) malloc(sizeof(complex)*32) ; 01352 aa = (complex *) malloc(sizeof(complex)*32) ; 01353 bb = (complex *) malloc(sizeof(complex)*32) ; 01354 01355 cs[0].r = 1.0 ; cs[0].i = 0.0 ; /* never used */ 01356 for( k=1 ; k < 32 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); } 01357 } 01358 01359 /*-- load subarrays, and FFT each one --*/ 01360 01361 for( i=k=0 ; k < 32 ; k++ ){ aa[k] = xc[i++] ; bb[k] = xc[i++] ; } 01362 01363 fft32( mode , aa ) ; fft32( mode , bb ) ; 01364 01365 /*-- recombine: 0 and Nyquist --*/ 01366 01367 xc[ 0].r = aa[0].r + bb[0].r ; xc[ 0].i = aa[0].i + bb[0].i ; 01368 xc[32].r = aa[0].r - bb[0].r ; xc[32].i = aa[0].i - bb[0].i ; 01369 01370 /*-- recombine: all others --*/ 01371 01372 if( mode > 0 ){ 01373 for( k=1 ; k < 32 ; k++ ){ 01374 bkr = bb[k].r; bki = bb[k].i; 01375 tr = cs[k].r; ti = cs[k].i; 01376 01377 t1 = tr*bkr - ti*bki ; t2 = tr*bki + ti*bkr ; 01378 01379 akr = aa[k].r; aki = aa[k].i; 01380 01381 xc[k ].r = akr+t1; xc[k ].i = aki+t2; 01382 xc[k+32].r = akr-t1; xc[k+32].i = aki-t2; 01383 } 01384 } else { 01385 for( k=1 ; k < 32 ; k++ ){ 01386 bkr = bb[k].r; bki = bb[k].i; 01387 tr = cs[k].r; ti = cs[k].i; 01388 /* only change from */ 01389 t1 = tr*bkr + ti*bki ; t2 = tr*bki - ti*bkr ; /* above is with ti */ 01390 01391 akr = aa[k].r; aki = aa[k].i; 01392 01393 xc[k ].r = akr+t1; xc[k ].i = aki+t2; 01394 xc[k+32].r = akr-t1; xc[k+32].i = aki-t2; 01395 } 01396 } 01397 01398 return ; 01399 } |
|
perhaps initialize * Definition at line 506 of file csfft.c. References csfft_trigconsts(), complex::i, nold, complex::r, and xc. Referenced by csfft_cox(), fft_3dec(), and fft_5dec().
00507 { 00508 register complex * csp , * xcx=xc; 00509 register float f1,f2,f3,f4 ; 00510 00511 /** perhaps initialize **/ 00512 00513 if( nold != 8 ) csfft_trigconsts( 8 ) ; 00514 00515 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00516 00517 /** data swapping part **/ 00518 00519 f1 = xcx[1].r ; f2 = xcx[1].i ; 00520 xcx[1].r = xcx[4].r ; xcx[1].i = xcx[4].i ; 00521 xcx[4].r = f1 ; xcx[4].i = f2 ; 00522 00523 f1 = xcx[3].r ; f2 = xcx[3].i ; 00524 xcx[3].r = xcx[6].r ; xcx[3].i = xcx[6].i ; 00525 xcx[6].r = f1 ; xcx[6].i = f2 ; 00526 00527 /** butterflying part **/ 00528 00529 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */ 00530 f2 = xcx[0].r ; f4 = xcx[0].i ; 00531 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ; 00532 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00533 00534 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */ 00535 f2 = xcx[2].r ; f4 = xcx[2].i ; 00536 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 00537 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00538 00539 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */ 00540 f2 = xcx[4].r ; f4 = xcx[4].i ; 00541 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 00542 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00543 00544 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */ 00545 f2 = xcx[6].r ; f4 = xcx[6].i ; 00546 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00547 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ; 00548 00549 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */ 00550 f2 = xcx[0].r ; f4 = xcx[0].i ; 00551 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ; 00552 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00553 00554 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */ 00555 f2 = xcx[4].r ; f4 = xcx[4].i ; 00556 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 00557 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ; 00558 00559 f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */ 00560 f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ; 00561 f2 = xcx[1].r ; f4 = xcx[1].i ; 00562 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ; 00563 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 00564 00565 f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */ 00566 f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ; 00567 f2 = xcx[5].r ; f4 = xcx[5].i ; 00568 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00569 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ; 00570 00571 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */ 00572 f2 = xcx[0].r ; f4 = xcx[0].i ; 00573 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ; 00574 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ; 00575 00576 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */ 00577 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ; 00578 f2 = xcx[1].r ; f4 = xcx[1].i ; 00579 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ; 00580 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ; 00581 00582 f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */ 00583 f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ; 00584 f2 = xcx[2].r ; f4 = xcx[2].i ; 00585 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ; 00586 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ; 00587 00588 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */ 00589 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ; 00590 f2 = xcx[3].r ; f4 = xcx[3].i ; 00591 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ; 00592 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ; 00593 00594 return ; 00595 } |
|
Definition at line 1867 of file csfft.c. References CC3, csfft_cox(), EXIT, fft1024, fft128(), fft16(), fft2, fft2048, fft256(), fft32(), fft4, fft512(), fft64(), fft8(), free, complex::i, i, M, M2, malloc, N, PI, complex::r, RMAX, SS3, and xc. Referenced by csfft_cox(), and csfft_many().
01868 { 01869 static int rec=0 ; 01870 static int rmold[RMAX] = {-1,-1,-1} ; 01871 static complex *rcs[RMAX] = {NULL,NULL,NULL} ; 01872 static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX] ; 01873 01874 int N=idim , M=idim/3 , M2=2*M ; 01875 int mold ; 01876 complex * cs=NULL , * aa , * bb , * cc ; 01877 register int k , i ; 01878 register float aar,aai, tr,ti, bbr,bbi, ccr,cci, 01879 t1,t2,t4,t5,t6,t8 ; 01880 01881 /*-- initialize cosine/sine table and memory space --*/ 01882 01883 if( rec >= RMAX ){ 01884 fprintf(stderr,"\n*** fft_3dec: too many recursions!\n"); EXIT(1); 01885 } 01886 01887 mold = rmold[rec] ; /* what was done before at this recursion level */ 01888 01889 if( M != mold ){ 01890 double th = (2.0*PI/N) ; 01891 if( M > mold ){ 01892 if( rcs[rec] != NULL ){ 01893 free(rcs[rec]);free(raa[rec]);free(rbb[rec]);free(rcc[rec]); 01894 } 01895 rcs[rec] = (complex *) malloc(sizeof(complex)*M2) ; 01896 raa[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01897 rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01898 rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01899 } 01900 01901 rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ; 01902 for( k=1 ; k < M2 ; k++ ){ 01903 rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th); 01904 } 01905 rmold[rec] = M ; 01906 } 01907 01908 cs = rcs[rec] ; aa = raa[rec] ; bb = rbb[rec] ; cc = rcc[rec] ; 01909 01910 /*-- load subarrays, and FFT each one --*/ 01911 01912 for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; } 01913 01914 #ifndef DONT_UNROLL_FFTS 01915 switch( M ){ 01916 case 1: break ; 01917 01918 case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ; break; 01919 case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ; break; 01920 case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ; break; 01921 case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; break; 01922 case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; break; 01923 case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; break; 01924 case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; break; 01925 case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; break; 01926 case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; break; 01927 01928 case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; break; 01929 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; break; 01930 01931 default: rec++ ; 01932 csfft_cox(mode,M,aa) ; 01933 csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ; 01934 rec-- ; break ; 01935 } 01936 #else 01937 rec++ ; 01938 csfft_cox(mode,M,aa) ; csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ; 01939 rec-- ; 01940 #endif 01941 01942 /*-- recombination --*/ 01943 01944 if( mode > 0 ){ 01945 for( k=0 ; k < M ; k++ ){ 01946 t1 = bb[k].r; t2 = bb[k].i; 01947 tr = cs[k].r; ti = cs[k].i; 01948 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ 01949 01950 t1 = cc[ k].r; t2 = cc[ k].i; 01951 tr = cs[2*k].r; ti = cs[2*k].i; 01952 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ 01953 01954 t4 = bbr+ccr ; t1 = t4*CC3 ; 01955 t8 = bbi+cci ; t6 = t8*CC3 ; 01956 t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3; 01957 01958 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01959 01960 xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ; 01961 xc[k+M ].r = (aar+t1)-t2 ; xc[k+M ].i = (aai+t6)+t5 ; 01962 xc[k+M2].r = (aar+t1)+t2 ; xc[k+M2].i = (aai+t6)-t5 ; 01963 } 01964 } else { 01965 for( k=0 ; k < M ; k++ ){ 01966 t1 = bb[k].r; t2 = bb[k].i; 01967 tr = cs[k].r; ti = cs[k].i; 01968 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ 01969 01970 t1 = cc[ k].r; t2 = cc[ k].i; 01971 tr = cs[2*k].r; ti = cs[2*k].i; 01972 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ 01973 01974 t4 = bbr+ccr ; t1 = t4*CC3 ; 01975 t8 = bbi+cci ; t6 = t8*CC3 ; 01976 t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3; 01977 01978 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01979 01980 xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ; 01981 xc[k+M ].r = (aar+t1)+t2 ; xc[k+M ].i = (aai+t6)-t5 ; 01982 xc[k+M2].r = (aar+t1)-t2 ; xc[k+M2].i = (aai+t6)+t5 ; 01983 } 01984 } 01985 01986 return ; 01987 } |
|
Definition at line 1714 of file csfft.c. References EXIT, fft256(), fft512(), free, complex::i, i, M, M2, M3, malloc, N, PI, complex::r, RMAX, and xc. Referenced by csfft_cox().
01715 { 01716 static int rec=0 ; 01717 static int rmold[RMAX] = {-1,-1,-1} ; 01718 static complex *rcs[RMAX] = {NULL,NULL,NULL} ; 01719 static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX] , *rdd[RMAX] ; 01720 01721 int N=idim , M=idim/4 , M2=2*M , M3=3*M ; 01722 int mold=-1 ; 01723 complex * cs , * aa , * bb , * cc , * dd ; 01724 register int k , i ; 01725 register float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 , 01726 acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ; 01727 01728 /*-- initialize cosine/sine table and memory space --*/ 01729 01730 if( rec >= RMAX ){ 01731 fprintf(stderr,"\n*** fft_4dec: too many recursions!\n"); EXIT(1); 01732 } 01733 01734 mold = rmold[rec] ; /* what was done before at this recursion level */ 01735 01736 if( M != mold ){ 01737 double th = (2.0*PI/N) ; 01738 if( M > mold ){ 01739 if( rcs[rec] != NULL ){ 01740 free(rcs[rec]);free(raa[rec]); 01741 free(rbb[rec]);free(rcc[rec]);free(rdd[rec]); 01742 } 01743 rcs[rec] = (complex *) malloc(sizeof(complex)*M3) ; 01744 raa[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01745 rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01746 rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01747 rdd[rec] = (complex *) malloc(sizeof(complex)*M ) ; 01748 } 01749 01750 rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ; 01751 for( k=1 ; k < M3 ; k++ ){ 01752 rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th); 01753 } 01754 rmold[rec] = M ; 01755 } 01756 01757 cs = rcs[rec]; aa = raa[rec]; bb = rbb[rec]; cc = rcc[rec]; dd = rdd[rec]; 01758 01759 /*-- load subarrays, and FFT each one --*/ 01760 01761 for( i=k=0; k < M; k++ ){ 01762 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; 01763 } 01764 01765 rec++ ; 01766 switch( N ){ 01767 01768 default: fprintf(stderr,"\n*** Illegal call to fft_4dec: N=%d\n",N); 01769 EXIT(1) ; 01770 01771 case 1024: fft256(mode,aa); fft256(mode,bb); 01772 fft256(mode,cc); fft256(mode,dd); break; 01773 01774 case 2048: fft512(mode,aa); fft512(mode,bb); 01775 fft512(mode,cc); fft512(mode,dd); break; 01776 01777 case 4096: fft_4dec(mode,1024,aa); fft_4dec(mode,1024,bb); /* recurse once */ 01778 fft_4dec(mode,1024,cc); fft_4dec(mode,1024,dd); break ; 01779 01780 case 8192: fft_4dec(mode,2048,aa); fft_4dec(mode,2048,bb); /* recurse once */ 01781 fft_4dec(mode,2048,cc); fft_4dec(mode,2048,dd); break ; 01782 01783 case 16384: fft_4dec(mode,4096,aa); fft_4dec(mode,4096,bb); /* recurse twice */ 01784 fft_4dec(mode,4096,cc); fft_4dec(mode,4096,dd); break ; 01785 01786 case 32768: fft_4dec(mode,8192,aa); fft_4dec(mode,8192,bb); /* recurse twice */ 01787 fft_4dec(mode,8192,cc); fft_4dec(mode,8192,dd); break ; 01788 } 01789 rec-- ; 01790 01791 /*-- recombination --*/ 01792 01793 if( mode > 0 ){ 01794 for( k=0 ; k < M ; k++ ){ 01795 t1 = bb[k].r; t2 = bb[k].i; 01796 tr = cs[k].r; ti = cs[k].i; 01797 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */ 01798 01799 t1 = cc[ k].r; t2 = cc[ k].i; 01800 tr = cs[2*k].r; ti = cs[2*k].i; 01801 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */ 01802 01803 t1 = dd[ k].r; t2 = dd[ k].i; 01804 tr = cs[3*k].r; ti = cs[3*k].i; 01805 ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */ 01806 01807 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01808 01809 acpr = aar + ccr ; acmr = aar - ccr ; 01810 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01811 acpi = aai + cci ; acmi = aai - cci ; 01812 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01813 01814 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01815 xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ; 01816 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01817 xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ; 01818 } 01819 } else { 01820 for( k=0 ; k < M ; k++ ){ 01821 t1 = bb[k].r; t2 = bb[k].i; 01822 tr = cs[k].r; ti = cs[k].i; 01823 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */ 01824 01825 t1 = cc[ k].r; t2 = cc[ k].i; 01826 tr = cs[2*k].r; ti = cs[2*k].i; 01827 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */ 01828 01829 t1 = dd[ k].r; t2 = dd[ k].i; 01830 tr = cs[3*k].r; ti = cs[3*k].i; 01831 ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */ 01832 01833 aar = aa[k].r; aai = aa[k].i; /* a[k] */ 01834 01835 acpr = aar + ccr ; acmr = aar - ccr ; 01836 bdpr = bbr + ddr ; bdmr = bbr - ddr ; 01837 acpi = aai + cci ; acmi = aai - cci ; 01838 bdpi = bbi + ddi ; bdmi = bbi - ddi ; 01839 01840 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ; 01841 xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ; 01842 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ; 01843 xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ; 01844 } 01845 } 01846 01847 return ; 01848 } |
|
Definition at line 2000 of file csfft.c. References COS72, csfft_cox(), EXIT, fft1024, fft128(), fft16(), fft2, fft2048, fft256(), fft32(), fft4, fft512(), fft64(), fft8(), free, complex::i, i, M, M2, M3, malloc, N, PI, complex::r, RMAX, s2, SIN72, and xc. Referenced by csfft_cox(), and csfft_many().
02001 { 02002 static int rec=0 ; 02003 static int rmold[RMAX] = {-1,-1,-1} ; 02004 static complex *rcs[RMAX] = {NULL,NULL,NULL} ; 02005 static complex *raa[RMAX], *rbb[RMAX], *rcc[RMAX], *rdd[RMAX], *ree[RMAX] ; 02006 02007 int N=idim , M=idim/5 , M2=2*M , M3=3*M , M4=4*M ; 02008 int mold ; 02009 complex * cs, *aa, *bb, *cc, *dd, *ee ; 02010 register int k , i ; 02011 register float aar,aai,bbr,bbi,ccr,cci,ddr,ddi,eer,eei ; 02012 register float akp,akm,bkp,bkm , ajp,ajm,bjp,bjm ; 02013 register float ak,bk,aj,bj ; 02014 float c72 , s72 , c2,s2 ; 02015 int ss ; 02016 02017 /*-- initialize cosine/sine table and memory space --*/ 02018 02019 if( rec >= RMAX ){ 02020 fprintf(stderr,"\n*** fft_5dec: too many recursions!\n"); EXIT(1); 02021 } 02022 02023 mold = rmold[rec] ; /* what was done before at this recursion level */ 02024 02025 if( M != mold ){ /* new value of M? */ 02026 double th = (2.0*PI/N) ; 02027 if( M > mold ){ /* need more space */ 02028 if( rcs[rec] != NULL ){ 02029 free(rcs[rec]);free(raa[rec]);free(rbb[rec]); 02030 free(rcc[rec]);free(rdd[rec]);free(ree[rec]); 02031 } 02032 rcs[rec] = (complex *) malloc(sizeof(complex)*M4) ; 02033 raa[rec] = (complex *) malloc(sizeof(complex)*M ) ; 02034 rbb[rec] = (complex *) malloc(sizeof(complex)*M ) ; 02035 rcc[rec] = (complex *) malloc(sizeof(complex)*M ) ; 02036 rdd[rec] = (complex *) malloc(sizeof(complex)*M ) ; 02037 ree[rec] = (complex *) malloc(sizeof(complex)*M ) ; 02038 } 02039 02040 rcs[rec][0].r = 1.0 ; rcs[rec][0].i = 0.0 ; 02041 for( k=1 ; k < M4 ; k++ ){ 02042 rcs[rec][k].r = cos(k*th); rcs[rec][k].i = sin(k*th); 02043 } 02044 rmold[rec] = M ; /* you must remember this */ 02045 } 02046 02047 cs = rcs[rec] ; aa = raa[rec] ; bb = rbb[rec] ; 02048 cc = rcc[rec] ; dd = rdd[rec] ; ee = ree[rec] ; 02049 02050 /*-- load subarrays, and FFT each one --*/ 02051 02052 for( i=k=0; k < M; k++ ){ 02053 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; ee[k]=xc[i++]; 02054 } 02055 02056 #ifndef DONT_UNROLL_FFTS 02057 switch( M ){ 02058 case 1: break ; 02059 02060 case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ; 02061 fft2 (mode,dd) ; fft2 (mode,ee) ; break; 02062 02063 case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ; 02064 fft4 (mode,dd) ; fft4 (mode,ee) ; break; 02065 02066 case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ; 02067 fft8 (mode,dd) ; fft8 (mode,ee) ; break; 02068 02069 case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; 02070 fft16 (mode,dd) ; fft16 (mode,ee) ; break; 02071 02072 case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; 02073 fft32 (mode,dd) ; fft32 (mode,ee) ; break; 02074 02075 case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; 02076 fft64 (mode,dd) ; fft64 (mode,ee) ; break; 02077 02078 case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; 02079 fft128(mode,dd) ; fft128(mode,ee) ; break; 02080 02081 case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; 02082 fft256(mode,dd) ; fft256(mode,ee) ; break; 02083 02084 case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; 02085 fft512(mode,dd) ; fft512(mode,ee) ; break; 02086 02087 case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; 02088 fft1024(mode,dd) ; fft1024(mode,ee) ; break; 02089 02090 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; 02091 fft2048(mode,dd) ; fft2048(mode,ee) ; break; 02092 02093 default: rec++ ; 02094 csfft_cox(mode,M,aa) ; 02095 csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ; 02096 csfft_cox(mode,M,dd) ; csfft_cox(mode,M,ee) ; 02097 rec-- ; break ; 02098 } 02099 #else 02100 rec++ ; 02101 csfft_cox(mode,M,aa) ; csfft_cox(mode,M,bb) ; csfft_cox(mode,M,cc) ; 02102 csfft_cox(mode,M,dd) ; csfft_cox(mode,M,ee) ; 02103 rec-- ; 02104 #endif 02105 02106 /*-- recombination --*/ 02107 02108 if( mode > 0 ){ s72 = SIN72 ; ss = 1 ; } 02109 else { s72 = -SIN72 ; ss = -1 ; } 02110 02111 c72 = COS72 ; 02112 c2 = c72 * c72 - s72 * s72 ; 02113 s2 = 2.0 * c72 * s72 ; 02114 02115 for( k=0 ; k < M ; k++ ){ 02116 ak = bb[k].r ; bk = bb[k].i ; 02117 aj = cs[k].r ; bj = cs[k].i * ss ; 02118 bbr = aj*ak - bj*bk ; bbi = aj*bk + bj*ak ; /* b[k]*exp(ss*2*Pi*k/N) */ 02119 02120 ak = cc[ k].r ; bk = cc[ k].i ; 02121 aj = cs[2*k].r ; bj = cs[2*k].i * ss ; 02122 ccr = aj*ak - bj*bk ; cci = aj*bk + bj*ak ; /* c[k]*exp(ss*4*Pi*k/N) */ 02123 02124 ak = dd[ k].r ; bk = dd[ k].i ; 02125 aj = cs[3*k].r ; bj = cs[3*k].i * ss ; 02126 ddr = aj*ak - bj*bk ; ddi = aj*bk + bj*ak ; /* d[k]*exp(ss*6*Pi*k/N) */ 02127 02128 ak = ee[ k].r ; bk = ee[ k].i ; 02129 aj = cs[4*k].r ; bj = cs[4*k].i * ss ; 02130 eer = aj*ak - bj*bk ; eei = aj*bk + bj*ak ; /* e[k]*exp(ss*8*Pi*k/N) */ 02131 02132 aar = aa[k].r ; aai = aa[k].i ; /* a[k] */ 02133 02134 /* the code below is (heavily) adapted from fftn.c */ 02135 02136 akp = bbr + eer ; akm = bbr - eer ; 02137 bkp = bbi + eei ; bkm = bbi - eei ; 02138 ajp = ccr + ddr ; ajm = ccr - ddr ; 02139 bjp = cci + ddi ; bjm = cci - ddi ; 02140 02141 xc[k].r = aar + akp + ajp ; xc[k].i = aai + bkp + bjp ; 02142 02143 ak = akp * c72 + ajp * c2 + aar ; bk = bkp * c72 + bjp * c2 + aai ; 02144 aj = akm * s72 + ajm * s2 ; bj = bkm * s72 + bjm * s2 ; 02145 02146 xc[k+M ].r = ak - bj ; xc[k+M ].i = bk + aj ; 02147 xc[k+M4].r = ak + bj ; xc[k+M4].i = bk - aj ; 02148 02149 ak = akp * c2 + ajp * c72 + aar ; bk = bkp * c2 + bjp * c72 + aai ; 02150 aj = akm * s2 - ajm * s72 ; bj = bkm * s2 - bjm * s72 ; 02151 02152 xc[k+M2].r = ak - bj ; xc[k+M2].i = bk + aj ; 02153 xc[k+M3].r = ak + bj ; xc[k+M3].i = bk - aj ; 02154 } 02155 02156 return ; 02157 } |
Variable Documentation
|
|
|
|
|
Definition at line 132 of file csfft.c. Referenced by csfft_cox(), csfft_many(), csfft_trigconsts(), fft16(), fft32(), and fft8(). |
|
Definition at line 72 of file csfft.c. Referenced by csfft_scale_inverse(). |
|
Definition at line 27 of file csfft.c. Referenced by csfft_use_fftw(). |