00001
00002
00029 #include "_xfft.h"
00030
00073 PUBLIC pRFFT XAPI rfft_construct( SPL_INT np, SPL_INT fftnp,
00074 SPL_BOOL inv, SPL_FLOAT ufactor,
00075 XFFT_WIN win, XFFT_WIN_FUNC winfunc )
00076 {
00077
00078
00079
00080
00081
00082
00083 pRFFT r;
00084 SPL_INT n, n2, n3;
00085 assert(np>0);
00086 assert(fftnp>=0);
00087
00088 r = (pRFFT)xfft_malloc(sizeof(RFFT));
00089 if (r!=NULL) {
00090 n = fft_n_bits((fftnp>0) ? fftnp : np);
00091 r->_nbits = ((n) ? n-1 : n);
00092 n2 = ( (n = (SPL_INT)1 << r->_nbits) >> 1 );
00093 n3 = n<<1;
00094 r->_np = ( (np>n3) ? n3 : np );
00095
00096 r->_revec = (SPL_pFLOAT)xfft_malloc(
00097 sizeof(SPL_FLOAT)*(n+1));
00098 r->_imvec = (SPL_pFLOAT)xfft_malloc(
00099 sizeof(SPL_FLOAT)*(n+1));
00100 r->_tinv = (SPL_pINT)xfft_malloc(sizeof(SPL_INT)*n);
00101 r->_half_tsin = (SPL_pFLOAT)xfft_malloc(
00102 sizeof(SPL_FLOAT)*n2);
00103 r->_tsin = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n2);
00104 r->_tcos = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*n2);
00105 if ((win!=XFFT_WIN_NONE)&&(r->_np))
00106 r->_winvec = (SPL_pFLOAT)xfft_malloc(
00107 sizeof(SPL_FLOAT)*r->_np);
00108 else
00109 r->_winvec = NULL;
00110
00111 if ( ((r->_revec==NULL)||(r->_imvec==NULL)||
00112 ((n)&&(r->_tinv==NULL))||
00113 ((n2)&&((r->_half_tsin==NULL)||
00114 (r->_tsin==NULL)||(r->_tcos==NULL)))||
00115 ((win!=XFFT_WIN_NONE)&&(r->_np)&&
00116 (r->_winvec==NULL))) == 0 ) {
00117
00118 if (n2) {
00119 fft_fill_half_tsin(r->_half_tsin,n2);
00120 fft_htsin_fill_tsin_tcos(r->_tsin,
00121 r->_tcos,r->_half_tsin,n2);
00122 }
00123 fft_fill_tinv(r->_tinv,n);
00124 r->_inv = inv;
00125 r->_ufactor = ufactor;
00126 r->_weight0pi=FALSE;
00127
00128 r->_owin = XFFT_WIN_NONE;
00129 rfft_setwin(r,win,winfunc);
00130 }
00131 else {
00132 rfft_destruct(r);
00133 r=NULL;
00134 }
00135 }
00136
00137 return r;
00138 }
00139
00150 PUBLIC SPL_VOID XAPI rfft_destruct( pRFFT r )
00151 {
00152 assert(r!=NULL);
00153
00154 if (r->_revec)
00155 xfft_free(r->_revec);
00156 if (r->_imvec)
00157 xfft_free(r->_imvec);
00158 if (r->_tinv)
00159 xfft_free(r->_tinv);
00160 if (r->_half_tsin)
00161 xfft_free(r->_half_tsin);
00162 if (r->_tsin)
00163 xfft_free(r->_tsin);
00164 if (r->_tcos)
00165 xfft_free(r->_tcos);
00166 if (r->_winvec)
00167 xfft_free(r->_winvec);
00168 xfft_free(r);
00169 }
00170
00198 PUBLIC SPL_BOOL XAPI rfft_setwin( pRFFT r, XFFT_WIN win,
00199 XFFT_WIN_FUNC winfunc )
00200 {
00201 INT i;
00202 assert(r!=NULL);
00203
00204 r->_win = win;
00205 r->_winfunc=winfunc;
00206 r->_ufactor0 = (r->_np) ? 2./r->_np : 1;
00207
00208 if ((win!=XFFT_WIN_NONE)&&(r->_np)) {
00209 if (r->_winvec==NULL) {
00210 r->_owin=XFFT_WIN_NONE;
00211 r->_winvec = (SPL_pFLOAT)xfft_malloc(
00212 sizeof(SPL_FLOAT)*r->_np);
00213 }
00214 if (r->_winvec!=NULL) {
00215 if (r->_owin!=r->_win) {
00216 switch (r->_win) {
00217 case XFFT_WIN_HAMM :
00218 win_hamm(r->_winvec,r->_np);
00219 break;
00220 case XFFT_WIN_HANN :
00221 win_hann(r->_winvec,r->_np);
00222 break;
00223 case XFFT_WIN_BLACK :
00224 win_black(r->_winvec,r->_np);
00225 break;
00226 case XFFT_WIN_BART :
00227 win_bart(r->_winvec,r->_np);
00228 break;
00229 case XFFT_WIN_KAIS5 :
00230 win_kais(r->_winvec,r->_np,5.0);
00231 break;
00232 case XFFT_WIN_KAIS6 :
00233 win_kais(r->_winvec,r->_np,6.0);
00234 break;
00235 case XFFT_WIN_KAIS7 :
00236 win_kais(r->_winvec,r->_np,7.0);
00237 break;
00238 case XFFT_WIN_KAIS8 :
00239 win_kais(r->_winvec,r->_np,8.0);
00240 break;
00241 case XFFT_WIN_KAIS9 :
00242 win_kais(r->_winvec,r->_np,9.0);
00243 break;
00244 case XFFT_WIN_KAIS10 :
00245 win_kais(r->_winvec,r->_np,10.0);
00246 break;
00247 case XFFT_WIN_USER :
00248 if ((r->_winfunc)!=NULL)
00249 (r->_winfunc)(r->_winvec,
00250 r->_np);
00251 else
00252 r->_win=XFFT_WIN_NONE;
00253 break;
00254 default :
00255
00256 assert(1==0);
00257 r->_win=XFFT_WIN_NONE;
00258 break;
00259 }
00260 r->_owin=r->_win;
00261 if (r->_win!=XFFT_WIN_NONE) {
00262 r->_winmean = 0;
00263 for (i=0; i<r->_np; i++) r->_winmean += r->_winvec[i];
00264 r->_winmean /= r->_np;
00265 }
00266 else
00267 r->_winmean = 1;
00268 }
00269 r->_ufactor0 /= r->_winmean;
00270 }
00271 else {
00272 r->_win=XFFT_WIN_NONE;
00273 return SPL_FALSE;
00274 }
00275 }
00276
00277 return SPL_TRUE;
00278 }
00279
00296 PUBLIC SPL_BOOL XAPI rfft_setnp( pRFFT r, SPL_INT np )
00297 {
00298 SPL_INT n;
00299 assert(r!=NULL);
00300 assert(np>0);
00301
00302 n = ((SPL_INT)1 << (r->_nbits+1));
00303 if (np>n)
00304 np = n;
00305 if (np!=r->_np) {
00306 r->_np = np;
00307 if (r->_winvec!=NULL) {
00308 xfft_free(r->_winvec);
00309 r->_winvec = NULL;
00310 }
00311 return rfft_setwin(r,r->_win,r->_winfunc);
00312 }
00313 return SPL_TRUE;
00314 }
00315
00349 PUBLIC SPL_BOOL XAPI rfft_setnpwin( pRFFT r, SPL_INT np,
00350 XFFT_WIN win, XFFT_WIN_FUNC winfunc )
00351 {
00352 assert(r!=NULL);
00353 assert(np>0);
00354
00355 r->_win = XFFT_WIN_NONE;
00356 rfft_setnp(r,np);
00357 return rfft_setwin(r,win,winfunc);
00358 }
00359
00373 PUBLIC SPL_VOID XAPI rfft_setinv( pRFFT r, SPL_BOOL inv )
00374 {
00375 assert(r!=NULL);
00376
00377 r->_inv = inv;
00378 }
00379
00392 PUBLIC SPL_VOID XAPI rfft_setufactor( pRFFT r, SPL_FLOAT ufactor )
00393 {
00394 assert(r!=NULL);
00395
00396 r->_ufactor = ufactor;
00397 }
00398
00411 PUBLIC SPL_VOID XAPI rfft_set0piweight( pRFFT r, SPL_BOOL weight )
00412 {
00413 assert(r!=NULL);
00414
00415 r->_weight0pi=weight;
00416 }
00417
00428 PUBLIC SPL_BOOL XAPI rfft_get0piweight( pRFFT r )
00429 {
00430 assert(r!=NULL);
00431
00432 return r->_weight0pi;
00433 }
00434
00445 PUBLIC SPL_INT XAPI rfft_getnp( pRFFT r )
00446 {
00447 assert(r!=NULL);
00448
00449 return r->_np;
00450 }
00451
00463
00464
00465
00466
00467 PUBLIC SPL_INT XAPI rfft_getfftnp( pRFFT r )
00468 {
00469 assert(r!=NULL);
00470
00471 return ((SPL_INT)1 << (r->_nbits+1));
00472 }
00473
00489 PUBLIC SPL_INT XAPI rfft_getvecnp( pRFFT r )
00490 {
00491 assert(r!=NULL);
00492
00493 return ((SPL_INT)1 << r->_nbits)+1;
00494 }
00495
00503 PUBLIC XFFT_WIN XAPI rfft_getwin( pRFFT r )
00504 {
00505 assert(r!=NULL);
00506
00507 return r->_win;
00508 }
00509
00521 PUBLIC SPL_pFLOAT XAPI rfft_getwinvec( pRFFT r )
00522 {
00523 assert(r!=NULL);
00524
00525 if (r->_win==XFFT_WIN_NONE)
00526 return NULL;
00527 return r->_winvec;
00528 }
00529
00537 PUBLIC XFFT_WIN_FUNC XAPI rfft_getwinfunc( pRFFT r )
00538 {
00539 assert(r!=NULL);
00540
00541 return r->_winfunc;
00542 }
00543
00553 PUBLIC SPL_BOOL XAPI rfft_getinv( pRFFT r )
00554 {
00555 assert(r!=NULL);
00556
00557 return r->_inv;
00558 }
00559
00567 PUBLIC SPL_FLOAT XAPI rfft_getufactor( pRFFT r )
00568 {
00569 assert(r!=NULL);
00570
00571 return r->_ufactor;
00572 }
00573
00591 PUBLIC SPL_pFLOAT XAPI rfft_getrevec( pRFFT r )
00592 {
00593 assert(r!=NULL);
00594
00595 return r->_revec;
00596 }
00597
00615 PUBLIC SPL_pFLOAT XAPI rfft_getimvec( pRFFT r )
00616 {
00617 assert(r!=NULL);
00618
00619 return r->_imvec;
00620 }
00621
00625 PRIVATE SPL_VOID rfft_dofft( pRFFT r )
00626 {
00627 SPL_INT n;
00628 SPL_FLOAT factor;
00629 assert(r!=NULL);
00630
00631 factor = r->_ufactor ? r->_ufactor : r->_ufactor0;
00632 n = ((SPL_INT)1 << r->_nbits);
00633
00634
00635 if ((r->_weight0pi) && (r->_inv!=FALSE) && (r->_ufactor==0)) {
00636 r->_revec[r->_tinv[0]] *= 2;
00637 r->_revec[r->_tinv[(SPL_INT)1 << (r->_nbits-1)]] *= 2;
00638
00639
00640 }
00641
00642 fft_fft( r->_nbits,r->_inv,r->_revec,r->_imvec,r->_tsin,r->_tcos,
00643 ((r->_inv)?factor*2.0:factor));
00644 fft_unscramble_cx_out(r->_revec,r->_imvec,n,r->_half_tsin,r->_inv);
00645 r->_revec[n] = r->_imvec[0];
00646 r->_imvec[0] = r->_imvec[n] = 0.0;
00647
00648
00649 if ((r->_weight0pi) && (r->_inv==FALSE) && (r->_ufactor==0)) {
00650 r->_revec[0] *= 0.5;
00651 r->_revec[n] *= 0.5;
00652 }
00653 }
00654
00658 #define __rfft_rfft(r,invec) \
00659 { \
00660 SPL_INT i, n2; \
00661 SPL_pINT tinv; \
00662 SPL_pFLOAT winv; \
00663 SPL_INT n = ((SPL_INT)1 << r->_nbits); \
00664 \
00665 \
00666 n2 = (r->_np) >> 1; \
00667 tinv = r->_tinv; \
00668 winv = r->_winvec; \
00669 if (r->_win!=XFFT_WIN_NONE) { \
00670 for (i=0; i<n2; i++) { \
00671 r->_revec[*tinv] = (*(invec++))*(*(winv++)); \
00672 r->_imvec[*(tinv++)] = (*(invec++))*(*(winv++)); \
00673 } \
00674 \
00675 if ((r->_np)&1) { \
00676 r->_revec[*tinv] = (*invec)*(*winv); \
00677 r->_imvec[*(tinv++)] = 0.0; \
00678 n2++; \
00679 } \
00680 } \
00681 else { \
00682 for (i=0; i<n2; i++) { \
00683 r->_revec[*tinv] = (*(invec++)); \
00684 r->_imvec[*(tinv++)] = (*(invec++)); \
00685 } \
00686 \
00687 if ((r->_np)&1) { \
00688 r->_revec[*tinv] = (*invec); \
00689 r->_imvec[*(tinv++)] = 0.0; \
00690 n2++; \
00691 } \
00692 } \
00693 for (i=n2; i<n; i++) { \
00694 r->_revec[*tinv] = 0.0; \
00695 r->_imvec[*(tinv++)] = 0.0; \
00696 } \
00697 rfft_dofft(r); \
00698 }
00699
00713 PUBLIC SPL_VOID XAPI rfft_rfft( pRFFT r, SPL_pFLOAT invec )
00714 {
00715 assert(r!=NULL);
00716
00717 __rfft_rfft(r,invec);
00718 }
00719
00733 PUBLIC SPL_VOID XAPI rfft_rfft_i( pRFFT r, SPL_pINT invec )
00734 {
00735 assert(r!=NULL);
00736
00737 __rfft_rfft(r,invec);
00738 }
00739
00753 PUBLIC SPL_VOID XAPI rfft_rfft_i16( pRFFT r, pINT16 invec )
00754 {
00755 assert(r!=NULL);
00756
00757 __rfft_rfft(r,invec);
00758 }
00759
00760
00774 PUBLIC SPL_VOID XAPI rfft_rfft_i32( pRFFT r, pINT32 invec )
00775 {
00776 assert(r!=NULL);
00777
00778 __rfft_rfft(r,invec);
00779 }
00780
00794 PUBLIC SPL_VOID XAPI rfft_rfft_u32( pRFFT r, pUINT32 invec )
00795 {
00796 assert(r!=NULL);
00797
00798 __rfft_rfft(r,invec);
00799 }
00800
00814 PUBLIC SPL_VOID XAPI rfft_rfft_f32( pRFFT r, pFLOAT32 invec )
00815 {
00816 assert(r!=NULL);
00817
00818 __rfft_rfft(r,invec);
00819 }
00820
00835 PUBLIC SPL_VOID XAPI rfft_mag( pRFFT r )
00836 {
00837 SPL_INT i, n;
00838 assert(r!=NULL);
00839
00840 n = ( (SPL_INT)1 << r->_nbits );
00841 for (i = 0; i<=n; i++)
00842 r->_revec[i] = __fft_cx_mag(r->_revec[i],r->_imvec[i]);
00843 }
00844
00859 PUBLIC SPL_VOID XAPI rfft_norm( pRFFT r )
00860 {
00861 SPL_INT i, n;
00862 assert(r!=NULL);
00863
00864 n = ( (SPL_INT)1 << r->_nbits );
00865 for (i = 0; i<=n; i++)
00866 r->_revec[i] = __fft_cx_norm(r->_revec[i],r->_imvec[i]);
00867 }
00868
00883 PUBLIC SPL_VOID XAPI rfft_arg( pRFFT r )
00884 {
00885 SPL_INT i, n;
00886 assert(r!=NULL);
00887
00888 n = ( (SPL_INT)1 << r->_nbits );
00889 for (i = 0; i<=n; i++)
00890 r->_imvec[i] = fft_zcx_arg(r->_revec[i],r->_imvec[i]);
00891 }
00892
00905 PUBLIC SPL_VOID XAPI rfft_magarg( pRFFT r )
00906 {
00907 SPL_INT i, n;
00908 SPL_FLOAT temp;
00909 assert(r!=NULL);
00910
00911 n = ( (SPL_INT)1 << r->_nbits );
00912 for (i = 0; i<=n; i++) {
00913 temp = r->_revec[i];
00914 r->_revec[i] = __fft_cx_mag(temp,r->_imvec[i]);
00915 r->_imvec[i] = fft_zcx_arg(temp,r->_imvec[i]);
00916 }
00917 }
00918
00931 PUBLIC SPL_VOID XAPI rfft_normarg( pRFFT r )
00932 {
00933 SPL_INT i, n;
00934 SPL_FLOAT temp;
00935 assert(r!=NULL);
00936
00937 n = ( (SPL_INT)1 << r->_nbits );
00938 for (i = 0; i<=n; i++) {
00939 temp = r->_revec[i];
00940 r->_revec[i] = __fft_cx_norm(temp,r->_imvec[i]);
00941 r->_imvec[i] = fft_zcx_arg(temp,r->_imvec[i]);
00942 }
00943 }
00944
00967 PUBLIC SPL_VOID XAPI rfft_trefmove_reim( pRFFT r, SPL_FLOAT nTs )
00968 {
00969 SPL_INT i, n;
00970 SPL_FLOAT d, a, b;
00971 assert(r!=NULL);
00972
00973 n = ( (SPL_INT)1 << r->_nbits );
00974 for (i = 0; i<=n; i++) {
00975 d = (M_PI*i*nTs)/n;
00976 a = cos(d);
00977 b = sin(d);
00978 d = r->_revec[i]*a - r->_imvec[i]*b;
00979 r->_imvec[i] = r->_imvec[i]*a+r->_revec[i]*b;
00980 r->_revec[i] = d;
00981 }
00982 }
00983
01006 PUBLIC SPL_VOID XAPI rfft_trefmove_arg( pRFFT r, SPL_FLOAT nTs )
01007 {
01008 SPL_INT i, n;
01009 assert(r!=NULL);
01010
01011 n = ( (SPL_INT)1 << r->_nbits );
01012 for (i = 0; i<=n; i++)
01013 r->_imvec[i] += ((M_PI*i*nTs)/n);
01014 }
01015
01038 PUBLIC SPL_VOID XAPI rfft_trefmove_argm( pRFFT r, SPL_FLOAT nTs )
01039 {
01040 SPL_INT i, n;
01041 assert(r!=NULL);
01042
01043 n = ( (SPL_INT)1 << r->_nbits );
01044 for (i = 0; i<=n; i++) {
01045 SPL_FLOAT d=fmod(r->_imvec[i] + (M_PI*i*nTs)/n, (2*M_PI));
01046 if (d>M_PI) d-= (2*M_PI);
01047 else if (d<-M_PI) d+= (2*M_PI);
01048 r->_imvec[i] = d;
01049 }
01050 }
01051
01052
01053
01054
