00001
00002
00068 #include "_xfft.h"
00069
00095 PUBLIC pCFFT XAPI cfft_construct( SPL_INT np, SPL_INT fftnp,
00096 SPL_BOOL inv, SPL_FLOAT ufactor )
00097 {
00098 pCFFT c;
00099 SPL_INT n, n2;
00100 assert(np>0);
00101 assert(fftnp>=0);
00102
00103 c = (pCFFT)xfft_malloc(sizeof(CFFT));
00104 if (c!=NULL) {
00105 c->_nbits = fft_n_bits((fftnp>0) ? fftnp : np);
00106 n2 = ( (n = (SPL_INT)1 << c->_nbits) >> 1 );
00107
00108 c->_revec = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n);
00109 c->_imvec = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n);
00110 c->_tinv = (SPL_pINT)xfft_malloc(sizeof(SPL_INT)*n);
00111 c->_tsin = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n2);
00112 c->_tcos = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n2);
00113
00114 if ((((n!=0)&&((c->_revec==NULL)||(c->_imvec==NULL)||
00115 (c->_tinv==NULL)))|| ((n2!=0)&&
00116 ((c->_tsin==NULL)||(c->_tcos==NULL))))==0) {
00117 fft_fill_tinv(c->_tinv,n);
00118 if (n2)
00119 fft_fill_tsin_tcos(c->_tsin,c->_tcos,n2);
00120 c->_inv = inv;
00121 c->_ufactor = ufactor;
00122 c->_np = ( (np>n) ? n : np );
00123 }
00124 else {
00125 cfft_destruct(c);
00126 c=NULL;
00127 }
00128 }
00129
00130 return c;
00131 }
00132
00144 PUBLIC SPL_VOID XAPI cfft_destruct( pCFFT c )
00145 {
00146 assert(c!=NULL);
00147
00148 if (c->_revec!=NULL)
00149 xfft_free(c->_revec);
00150 if (c->_imvec!=NULL)
00151 xfft_free(c->_imvec);
00152 if (c->_tinv!=NULL)
00153 xfft_free(c->_tinv);
00154 if (c->_tsin!=NULL)
00155 xfft_free(c->_tsin);
00156 if (c->_tcos!=NULL)
00157 xfft_free(c->_tcos);
00158 xfft_free(c);
00159 }
00160
00177 PUBLIC SPL_VOID XAPI cfft_setnp( pCFFT c, SPL_INT np )
00178 {
00179 SPL_INT n;
00180 assert(c!=NULL);
00181 assert(np>0);
00182
00183 n = ((SPL_INT)1 << (c->_nbits));
00184
00185 c->_np = ((np>n) ? n : np );
00186 }
00187
00201 PUBLIC SPL_VOID XAPI cfft_setinv( pCFFT c, SPL_BOOL inv )
00202 {
00203 assert(c!=NULL);
00204
00205 c->_inv = inv;
00206 }
00207
00220 PUBLIC SPL_VOID XAPI cfft_setufactor( pCFFT c, SPL_FLOAT ufactor )
00221 {
00222 assert(c!=NULL);
00223
00224 c->_ufactor = ufactor;
00225 }
00226
00237 PUBLIC SPL_INT XAPI cfft_getnp( pCFFT c )
00238 {
00239 assert(c!=NULL);
00240
00241 return c->_np;
00242 }
00243
00254 PUBLIC SPL_INT XAPI cfft_getfftnp( pCFFT c )
00255 {
00256 assert(c!=NULL);
00257
00258 return ((SPL_INT)1 << (c->_nbits));
00259 }
00260
00271 PUBLIC SPL_INT XAPI cfft_getvecnp( pCFFT c )
00272 {
00273 assert(c!=NULL);
00274
00275 return ((SPL_INT)1 << (c->_nbits));
00276 }
00277
00287 PUBLIC SPL_BOOL XAPI cfft_getinv( pCFFT c )
00288 {
00289 assert(c!=NULL);
00290
00291 return c->_inv;
00292 }
00293
00301 PUBLIC SPL_FLOAT XAPI cfft_getufactor( pCFFT c )
00302 {
00303 assert(c!=NULL);
00304
00305 return c->_ufactor;
00306 }
00307
00325 PUBLIC SPL_pFLOAT XAPI cfft_getrevec( pCFFT c )
00326 {
00327 assert(c!=NULL);
00328
00329 return c->_revec;
00330 }
00331
00349 PUBLIC SPL_pFLOAT XAPI cfft_getimvec( pCFFT c )
00350 {
00351 assert(c!=NULL);
00352
00353 return c->_imvec;
00354 }
00355
00359 #define __cfft_cfft(c,re_invec,im_invec) \
00360 { \
00361 SPL_INT i, n; \
00362 SPL_pINT tinv; \
00363 \
00364 n = ((SPL_INT)1 << c->_nbits); \
00365 tinv = c->_tinv; \
00366 \
00367 for (i=0; i<c->_np; i++) { \
00368 c->_revec[*tinv] = (*(re_invec++)); \
00369 c->_imvec[*(tinv++)] = (*(im_invec++)); \
00370 } \
00371 \
00372 for (i=c->_np; i<n; i++) { \
00373 c->_revec[*tinv] = 0.0; \
00374 c->_imvec[*(tinv++)] = 0.0; \
00375 } \
00376 fft_fft(c->_nbits,c->_inv,c->_revec,c->_imvec,c->_tsin, \
00377 c->_tcos,c->_ufactor); \
00378 }
00379
00396 PUBLIC SPL_VOID XAPI cfft_cfft( pCFFT c,
00397 SPL_pFLOAT re_invec, SPL_pFLOAT im_invec )
00398 {
00399 assert(c!=NULL);
00400
00401 __cfft_cfft(c,re_invec,im_invec);
00402 }
00403
00418 PUBLIC SPL_VOID XAPI cfft_mag( pCFFT c )
00419 {
00420 SPL_INT i, n;
00421 assert(c!=NULL);
00422
00423 n = ( (SPL_INT)1 << c->_nbits );
00424 for (i = 0; i<n; i++)
00425 c->_revec[i] = __fft_cx_mag(c->_revec[i],c->_imvec[i]);
00426 }
00427
00442 PUBLIC SPL_VOID XAPI cfft_norm( pCFFT c )
00443 {
00444 SPL_INT i, n;
00445 assert(c!=NULL);
00446
00447 n = ( (SPL_INT)1 << c->_nbits );
00448 for (i = 0; i<n; i++)
00449 c->_revec[i] = __fft_cx_norm(c->_revec[i],c->_imvec[i]);
00450 }
00451
00466 PUBLIC SPL_VOID XAPI cfft_arg( pCFFT c )
00467 {
00468 SPL_INT i, n;
00469 assert(c!=NULL);
00470
00471 n = ( (SPL_INT)1 << c->_nbits );
00472 for (i = 0; i<n; i++)
00473 c->_imvec[i] = fft_zcx_arg(c->_revec[i],c->_imvec[i]);
00474 }
00475
00488 PUBLIC SPL_VOID XAPI cfft_magarg( pCFFT c )
00489 {
00490 SPL_INT i, n;
00491 SPL_FLOAT temp;
00492 assert(c!=NULL);
00493
00494 n = ( (SPL_INT)1 << c->_nbits );
00495 for (i = 0; i<n; i++) {
00496 temp = c->_revec[i];
00497 c->_revec[i] = __fft_cx_mag(temp,c->_imvec[i]);
00498 c->_imvec[i] = fft_zcx_arg(temp,c->_imvec[i]);
00499 }
00500 }
00501
00514 PUBLIC SPL_VOID XAPI cfft_normarg( pCFFT c )
00515 {
00516 SPL_INT i, n;
00517 SPL_FLOAT temp;
00518 assert(c!=NULL);
00519
00520 n = ( (SPL_INT)1 << c->_nbits );
00521 for (i = 0; i<n; i++) {
00522 temp = c->_revec[i];
00523 c->_revec[i] = __fft_cx_norm(temp,c->_imvec[i]);
00524 c->_imvec[i] = fft_zcx_arg(temp,c->_imvec[i]);
00525 }
00526 }
00527
00550 PUBLIC SPL_VOID XAPI cfft_trefmove_reim( pCFFT c, SPL_FLOAT nTs )
00551 {
00552 SPL_INT i, n;
00553 SPL_FLOAT d, a, b;
00554 assert(c!=NULL);
00555
00556 n = ( (SPL_INT)1 << c->_nbits );
00557 for (i = 0; i<n; i++) {
00558 d = (2*M_PI*i*nTs)/n;
00559 a = cos(d);
00560 b = sin(d);
00561 d = c->_revec[i]*a - c->_imvec[i]*b;
00562 c->_imvec[i] = c->_imvec[i]*a+c->_revec[i]*b;
00563 c->_revec[i] = d;
00564 }
00565 }
00566
00589 PUBLIC SPL_VOID XAPI cfft_trefmove_arg( pCFFT c, SPL_FLOAT nTs )
00590 {
00591 SPL_INT i, n;
00592 assert(c!=NULL);
00593
00594 n = ( (SPL_INT)1 << c->_nbits );
00595 for (i = 0; i<n; i++)
00596 c->_imvec[i] += ((2*M_PI*i*nTs)/n);
00597 }
00598
00621 PUBLIC SPL_VOID XAPI cfft_trefmove_argm( pCFFT c, SPL_FLOAT nTs )
00622 {
00623 SPL_INT i, n;
00624 assert(c!=NULL);
00625
00626 n = ( (SPL_INT)1 << c->_nbits );
00627 for (i = 0; i<n; i++) {
00628 SPL_FLOAT d=fmod(c->_imvec[i] + ((2*M_PI)*i*nTs)/n, (2*M_PI));
00629 if (d>M_PI) d-= (2*M_PI);
00630 else if (d<-M_PI) d+= (2*M_PI);
00631 c->_imvec[i] = d;
00632 }
00633 }
00634
00635
00636
