/* * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - INRIA - Allan SIMON * * This file must be used under the terms of the CeCILL. * This source file is licensed as described in the file COPYING, which * you should have received as part of this distribution. The terms * are also available at * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt * */ #include #include #include "max.h" #include "min.h" #include "fft_internal.h" /*c'est moche je sais */ static double* a ; static double* b ; static int ntot ; static int n ; static int nspan ; static int isn ; static int m ; static int kt ; static double* wt ; static double* ck ; static double* bt ; static double* sk ; static int* np ; static int* nfac ; static int inc ; static int nt ; static int ks ; static double rad ; static double c72 ; static double s72 ; static double s120 ; static double aa ; static double ak ; static double akm ; static double akp ; static double aj ; static double ajp ; static double ajm ; static double bb ; static double bk ; static double bkm ; static double bkp ; static double bj ; static double bjp ; static double bjm ; static double dr ; static double cd ; static double c1 ; static double c2 ; static double c3 ; static double sd ; static double s1 ; static double s2 ; static double s3 ; static int kspan ; static int nn ; static int jc ; static int klim ; static int jf ; static int maxf ; static int kk ; static int k ; static int k1 ; static int k2 ; static int k3 ; static int k4 ; static int mm ; static int kspnn ; static int i ; static int j ; static int jj; /* Prototypes */ static void preliminaryWork (void); static void permute_stage1 (void); static void permute_stage2 (void); static void f4t_150 (void); static void factorOf3Transform (void) ; static void factorOf5Transform (void) ; static void preFOtherTransform (void); static void factorOfOtherTransform (void); static void pre_sqFactor2NormlOrder (void); static void nonSqFactor2NormOrder (void) ; static void detPermutCycles (void); static void reorderMatrix (void ) ; static int f4t_170 (void); static int factorTransform (void); static int pre_fOf2Trans (void); static int factorOf2Transform (void); static int factorOf4Transform (void); static int mulByRotationFactor (void ); static int post_sqFactor2NormlOrder (void); static void single_sqFactor2NormlOrder (void); static int multi_sqFactor2NormlOrder (void); /* End Prototypes */ /*note on this code all numbers alone in comment is a reference to the corresponding goto in the original fotran code */ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, int _iIsn, int _iM, int _iKt, double* _pdblWt, double* _pdblCk, double* _pdblBt, double* _pdblSk, int* _piNp, int* _piNfac) { int retVal = 0 ; a = _pdblA ; b = _pdblB ; ntot = _iNtot ; n = _iN ; nspan= _iNspan ; isn = _iIsn; m = _iM ; kt = _iKt ; wt = _pdblWt ; ck = _pdblCk; bt = _pdblBt; sk = _pdblSk; np = _piNp; nfac = _piNfac; inc = abs ( isn ) ; nt = inc*ntot ; ks = inc*nspan; rad = atan((double)1); c72 = cos (rad/0.6250); s72 = sin (rad/0.6250); s120= sqrt(0.750); preliminaryWork() ; while ( retVal == 0 ) retVal = factorTransform ( ) ; np[0] = ks ; if ( kt != 0) permute_stage1 ( ) ; if ( 2*kt + 1 < m ) permute_stage2 ( ) ; _pdblA = a ; _pdblB = b ; return 0 ; } /** ************************************** Sous-Fonctions ******************************************/ /* this function only set the value of variable */ static void preliminaryWork (void) { s72 = -s72 ; s120= -s120; rad = -rad ; kspan = ks ; nn = nt -inc ; jc = ks/n ; /* sin , cos values are re-initialized each lim steps */ klim = 32*jc ; i = 0; jf = 0 ; if ( kt > 0 ) maxf = max ( nfac[kt-1] , nfac[m-kt-1] ); else maxf = nfac[m-kt-1] ; } /*40*/ /* this function is call as many time as dfftbi has determined factor for the size of the input vector each time we call a transform function for each kind of factor , we begin by the smallest factor are stored in nfac */ static int factorTransform (void) { int retVal = 42; dr = 8 * (double)jc/(double)kspan ; cd = 2 * sin(0.5*dr*rad)*sin(0.5*dr*rad); sd = sin(dr*rad) ; kk = 1 ; i++ ; switch ( nfac[i-1] ) { case 2 : /*transform for factor of 2 (including rotation factor)*/ retVal = pre_fOf2Trans() ; if ( retVal == 0 ) factorOf2Transform () ; break ; case 4 : kspnn = kspan ; kspan = kspan >> 2 ; /*kspan /= 4 */ retVal = factorOf4Transform () ; break ; case 3 : k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; factorOf3Transform ( ) ; break ; case 5 : k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; factorOf5Transform ( ) ; break ; default : k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; if ( nfac[i-1] != jf) preFOtherTransform ( ) ; factorOfOtherTransform ( ) ; break ; } if ( retVal == 42 ) { if ( i != m) retVal = mulByRotationFactor ( ) ; else retVal = 1 ; } if ( retVal == 1 ) return 1 ; /*goto permute */ else return 0 ; /*goto factor_transform => once again*/ } /* permutation for square factor of n */ static void permute_stage1 (void) { int retVal = 1 ; pre_sqFactor2NormlOrder () ; if ( n == ntot ) /*permutation for single-variate transform (optional code)*/ while ( retVal == 1) { single_sqFactor2NormlOrder () ; retVal = post_sqFactor2NormlOrder () ; } else /*permutation for multivariate transform*/ while ( retVal == 1) retVal = multi_sqFactor2NormlOrder (); } static void permute_stage2 (void) { kspnn = np[kt] ; /*permutation for square-free facotrs of n */ nonSqFactor2NormOrder () ; /*determine the permutation cycles of length greater than 1*/ detPermutCycles (); j = k3 + 1; nt -= kspnn ; i = nt - inc + 1 ; while ( nt >= 0 ) { reorderMatrix ( ) ; j = k3 + 1 ; nt -= kspnn ; i = nt - inc + 1 ; } } /***************************************** Sous-Sous-Fonctions ******************************************/ static int pre_fOf2Trans (void) { kspan /= 2; k1 = kspan + 2 ; /*50*/ do{ do{ k2 = kk + kspan ; ak = a[k2-1] ; bk = b[k2-1] ; a[k2-1] = a[kk-1] - ak; b[k2-1] = b[kk-1] - bk; a[kk-1] = a[kk-1] + ak; b[kk-1] = b[kk-1] + bk; kk = k2 + kspan ; }while (kk <= nn); kk -= nn ; }while (kk <= jc); if ( kk > kspan ) return 1 ; /*goto350*/ else return 0 ; /*goto60*/ } static int factorOf2Transform (void) { do /*60*/ {/*while ( kk <= jc*2 )*/ c1 = 1 - cd ; s1 = sd ; mm = min( k1/2 , klim); do/* do 80 */ {/*while ( kk <= mm || ( kk > mm && kk < k2 ))*/ do {/*while(kk > k2) */ do { /*while ( kk < nt )*/ k2 = kk + kspan; ak = a[kk-1] - a[k2-1]; bk = b[kk-1] - b[k2-1]; a[kk-1] = a[kk-1] + a[k2-1]; b[kk-1] = b[kk-1] + b[k2-1]; a[k2-1] = c1*ak - s1*bk; b[k2-1] = s1*ak + c1*bk; kk = k2 + kspan; }while ( kk < nt ); k2 = kk - nt; c1 = -c1; kk = k1 - k2; }while (kk > k2); kk += jc; if ( kk <= mm ) /* 70 */ { ak = c1 - ( cd*c1+sd*s1) ; s1 += (sd*c1-cd*s1) ; /*c the following three statements compensate for truncation c error. if rounded arithmetic is used, substitute c c1=ak*/ c1 = 0.5/(ak*ak+s1*s1) + 0.5 ; s1 *= c1 ; c1 *= ak ; } else { if ( kk < k2 ) /*90*/ { s1 = dr*rad*((double)(kk-1)/(double)jc); c1 = cos(s1) ; s1 = sin(s1) ; mm = min(k1/2,mm+klim); } } } while ( kk <= mm || ( kk > mm && kk < k2 )); k1 += (inc+inc) ; kk = (k1-kspan)/2 + jc; } while ( kk <= jc*2 ); return 0 ; /*goto40*/ } /* this one is just an optimisation of the factor of 2 transform , we compute more things each turn */ static int factorOf4Transform (void) { int return_value = 0 ; /*120*/ do { c1 = 1 ; s1 = 0 ; mm = min ( kspan , klim ) ; do { f4t_150 () ; return_value = f4t_170 () ; } while ( return_value == 0 ); kk += ( inc - kspan ) ; } while ( kk <= jc ) ; if ( kspan == jc ) return 1 ; /*goto350*/ else return 0 ;/*goto40*/ } /*this function and the following are just here for conveniance , they just do fourier transformation for factor of 4 but as the code was a bit long in factorof4transform , we've created two sub-functions */ static void f4t_150 (void) { do{ k1 = kk + kspan ; k2 = k1 + kspan ; k3 = k2 + kspan ; akp = a[kk-1] + a[k2-1] ; akm = a[kk-1] - a[k2-1] ; ajp = a[k1-1] + a[k3-1] ; ajm = a[k1-1] - a[k3-1] ; a[kk-1] = akp + ajp ; ajp = akp - ajp ; bkp = b[kk-1] + b[k2-1] ; bkm = b[kk-1] - b[k2-1] ; bjp = b[k1-1] + b[k3-1] ; bjm = b[k1-1] - b[k3-1] ; b[kk-1] = bkp + bjp ; bjp = bkp - bjp ; akp = akm + bjm ; akm = akm - bjm ; bkp = bkm - ajm ; bkm = bkm + ajm ; if ( s1 == 0 )/*190*/ { a[k1-1] = akp ; a[k2-1] = ajp ; a[k3-1] = akm ; b[k1-1] = bkp ; b[k2-1] = bjp ; b[k3-1] = bkm ; } else /*160*/ { a[k1-1] = akp*c1 - bkp*s1 ; a[k2-1] = ajp*c2 - bjp*s2 ; a[k3-1] = akm*c3 - bkm*s3 ; a[k1-1] = bkp*c1 + akp*s1 ; a[k2-1] = bjp*c2 + ajp*s2 ; a[k3-1] = bkm*c3 + akm*s3 ; } kk=k3+kspan; }while ( kk <= nt ) ; } static int f4t_170 (void) { kk += ( jc - nt ) ; if ( kk <= mm ) { c2 = c1 - (cd*c1 + sd*s1); s1 = s1 + (sd*c1 - cd*s1); /* the following three statements compensate for truncation error. if rounded arithmetic is used, substitute c1=c2 */ c1 = 0.5/(c2*c2+s1*s1) + 0.5 ; s1 *= c1 ; c1 *= c2 ; /*140*/ c2 = c1*c1 - s1*s1 ; s2 = c1*s1*2 ; c3 = c2*c1 - s2*s1 ; s3 = c2*s1 + s2*c1 ; return 0 ; } else { if ( kk <= kspan ) { s1 = dr*rad * (kk-1)/jc ; c1 = cos (s1) ; s1 = sin (s1) ; mm = min ( kspan , mm + klim ); /*140*/ c2 = c1*c1 - s1*s1 ; s2 = c1*s1*2 ; c3 = c2*c1 - s2*s1 ; s3 = c2*s1 + s2*c1 ; return 0 ; } } return 1 ; } static void factorOf3Transform (void) { do{ do{ k1 = kk + kspan ; k2 = k1 + kspan ; ak = a[kk-1] ; bk = b[kk-1] ; aj = a[k1-1] + a[k2-1] ; bj = b[k1-1] + b[k2-1] ; a[kk-1] = ak + aj ; b[kk-1] = bk + bj ; ak = -0.5*aj + ak ; bk = -0.5*bj + bk ; aj = (a[k1-1] - a[k2-1])*s120 ; bj = (b[k1-1] - b[k2-1])*s120 ; a[k1-1] = ak - bj ; b[k1-1] = bk + aj ; a[k2-1] = ak + bj ; b[k2-1] = bk - aj ; kk = k2 + kspan ; } while (kk < nn); kk -= nn ; }while (kk <= kspan); } static void factorOf5Transform (void) { c2 = c72*c72 - s72 *s72 ; s2 = 2 * c72*s72; do{ do{ k1 = kk + kspan ; k2 = k1 + kspan ; k3 = k2 + kspan ; k4 = k3 + kspan ; akp = a[k1-1] + a[k4-1] ; akm = a[k1-1] - a[k4-1] ; bkp = b[k1-1] + b[k4-1] ; bkm = b[k1-1] - b[k4-1] ; ajp = a[k2-1] + a[k3-1] ; ajm = a[k2-1] - a[k3-1] ; bjp = b[k2-1] + b[k3-1] ; bjm = b[k2-1] - b[k3-1] ; aa = a[kk-1] ; bb = b[kk-1] ; a[kk-1] = aa + akp + ajp; b[kk-1] = bb + bkp + bjp; ak = akp*c72 + ajp*c2 + aa ; bk = bkp*c72 + bjp*c2 + bb ; aj = akm*s72 + ajm*s2 ; bj = bkm*s72 + bjm*s2 ; a[k1-1] = ak - bj ; a[k4-1] = ak + bj ; b[k1-1] = bk + aj ; b[k4-1] = bk - aj ; ak = akp*c2 + ajp*c72 + aa ; bk = bkp*c2 + bjp*c72 + bb ; aj = akm*s2 - ajm*s72 ; bj = bkm*s2 - bjm*s72 ; a[k2-1] = ak - bj ; a[k3-1] = ak + bj ; b[k2-1] = bk + aj ; b[k3-1] = bk - aj ; kk = k4 + kspan; }while (kk < nn); kk -= nn ; }while (kk <= kspan); } /* this function is the general case of non factor of 2 factor , the factorof3transform and factorof5trandform are just special case of this one */ static void preFOtherTransform (void) { jf = k ; s1 = (rad*8)/k ; c1 = cos (s1) ; s1 = sin (s1) ; ck[jf-1] = 1 ; sk[jf-1] = 0 ; j = 1 ; do { ck[j-1] = ck[k-1] * c1 + sk[k-1]*s1 ; sk[j-1] = ck[k-1] * s1 - sk[k-1]*c1 ; k -- ; ck[k-1] = ck[j-1] ; sk[k-1] = - sk[j-1] ; j++ ; }while ( j < k ); } static void factorOfOtherTransform (void) { int ktemp = 0 ; do { k1 = kk ; k2 = kk + kspnn ; aa = a[kk-1] ; bb = b[kk-1] ; ak = aa ; bk = bb ; j = 1 ; k1 += kspan ; do { k2 -= kspan ; j++ ; wt[j-1] = a[k1-1] + a[k2-1] ; ak = wt[j-1] + ak ; bt[j-1] = b[k1-1] + b[k2-1] ; bk = bt[j-1] + bk ; j++ ; wt[j-1] = a[k1-1] - a[k2-1] ; bt[j-1] = b[k1-1] - b[k2-1] ; k1 += kspan; }while ( k1 < k2 ) ; a[kk-1] = ak ; b[kk-1] = bk ; k1 = kk ; k2 = kk + kspnn ; j = 1 ; do { k1 += kspan ; k2 -= kspan ; jj = j ; ak = aa ; bk = bb ; aj = 0 ; bj = 0 ; k = 1 ; do { k++ ; ak += ( wt[k-1] * ck[jj-1] ) ; bk += ( bt[k-1] * ck[jj-1] ) ; k++ ; aj += (wt[k-1] * sk[jj-1]) ; bj += (bt[k-1] * sk[jj-1]) ; jj += j ; if ( jj > jf ) jj -= jf ; } while ( k < jf ) ; k = jf - j ; a[k1-1] = ak - bj ; b[k1-1] = bk + aj ; a[k2-1] = ak + bj ; b[k2-1] = bk - aj ; j++ ; }while ( j < k ) ; kk += kspnn ; ktemp = kk ; if ( kk > nn ) { kk -= nn; } }while ( ktemp <= nn || (kk <= kspan && ( ktemp > nn)) ); } static int mulByRotationFactor (void ) { int ktemp = 0 ; if ( i != m ) { kk = jc + 1 ; /*300*/ do { c2 = 1 - cd ; s1 = sd ; mm = min ( kspan , klim ) ; /*320 */ do { c1 = c2 ; s2 = s1 ; kk += kspan ; do { ak = a[kk-1] ; a[kk-1] = c2*ak - s2*b[kk-1] ; b[kk-1] = s2*ak + c2*b[kk-1] ; kk += kspnn ; ktemp = kk ; if ( kk > nt ) { ak = s1*s2 ; s2 = s1*c2 + s2*c1 ; c2 = c1*c2 - ak ; kk += (kspan - nt ) ; } }while (ktemp <= nt || ( kk <= kspnn && ktemp > nt )) ; kk += ( jc - kspnn ); if ( kk <= mm ) { /* 310*/ c2 = c1 - ( cd*c1 + sd*s1 ) ; s1 += (sd*c1 - cd*s1 ) ; /* the following three statements compensate for truncation error. if rounded arithmetic is used, substitute c1=c2 */ c1 = 0.5/(c2*c2+s1*s1) + 0.5 ; s1 *= c1 ; c2 *= c1 ; } else { if ( kk <= kspan ) { s1 = dr*rad * (kk-1)/jc ; c2 = cos (s1) ; s1 = sin (s1) ; mm = min ( kspan , mm + klim ); } } }while ( kk <= mm ||( kk <= kspan && kk > mm ) ) ; kk += (jc + inc -kspan ); }while ( kk <= jc+jc); return 0 ; /* goto40 */ } return 1 ; /* goto350*/ } static void pre_sqFactor2NormlOrder (void) { k = kt + kt + 1 ; if ( m < k ) k -- ; j = 1 ; np[k] = jc ; do { np[j] = np[j-1]/nfac[j-1] ; np[k-1] = np[k]*nfac[j-1] ; j++ ; k-- ; }while ( j < k ) ; k3 = np[k] ; kspan = np[1] ; kk = jc + 1 ; k2 = kspan + 1 ; j = 1; } static int post_sqFactor2NormlOrder (void) { do { do { k2 -= np[j-1] ; j++ ; k2 += np[j] ; } while ( k2 > np[j-1]); j = 1 ; /* 390 */ do { if ( kk < k2 ) { return 1 ; } else { kk += inc ; k2 += kspan ; } }while( k2 < ks ); }while ( kk < ks ) ; jc = k3 ; return 0; } /* appeler cetter fonction dans un do while valeur_retour != 1)*/ static void single_sqFactor2NormlOrder (void) { do { ak = a[kk-1] ; a[kk-1] = a[k2-1] ; a[k2-1] = ak ; bk = b[kk-1] ; b[kk-1] = b[k2-1] ; b[k2-1] = bk ; kk += inc ; k2 += kspan ; } while ( k2 < ks ); /*380*/ } /*idem que single_ */ static int multi_sqFactor2NormlOrder (void) { k = kk + jc ; do /*410*/ { ak = a[kk-1] ; a[kk-1] = a[k2-1] ; a[k2-1] = ak ; bk = b[kk-1] ; b[kk-1] = b[k2-1] ; b[k2-1] = bk ; kk += inc ; k2 += kspan ; } while ( kk < k ); kk += (ks - jc ) ; k2 += (ks - jc ) ; if ( kk < nt ) return 1 ; k2 += ( kspan - nt ); kk += ( jc - nt ); if ( k2 < ks ) { return 1 ; } if( post_sqFactor2NormlOrder ( ) == 1 ) { return 1 ; } jc = k3 ; return 0; } static void nonSqFactor2NormOrder (void) { j = m - kt ; nfac[j] = 1 ; do { nfac[j-1] *= nfac[j] ; j-- ; }while ( j != kt ) ; kt ++ ; nn = nfac[kt-1] - 1; jj = 0 ; j = 0; /*480*/ k2 = nfac[kt-1] ; k = kt + 1 ; kk = nfac[k-1] ; j ++ ; while ( j <= nn ) { jj += kk ; while ( jj >= k2 ) { jj -= k2 ; k2 = kk ; k++ ; kk = nfac[k-1] ; jj += kk ; } np[j-1] = jj ; k2 = nfac[kt-1] ; k = kt + 1 ; kk = nfac[k-1] ; j ++ ; } j = 0 ; return ; } /* here we determine how many permutation cycles we need to do */ static void detPermutCycles (void) { do { do { j++ ; kk = np[j-1] ; }while ( kk < 0 ) ; if ( kk != j ) { do { k = kk ; kk = np[k-1] ; np[k-1] = -kk ; }while ( kk != j ) ; k3 = kk ; } else np[j-1] = -j ; }while ( j != nn ); maxf *= inc ; return ; } static void reorderMatrix (void) { do { do { j-- ; }while (np[j-1] < 0 ) ; jj = jc ; /*520*/ do { kspan = jj ; if ( jj > maxf ) kspan = maxf ; jj -= kspan ; k = np [j-1]; kk = jc*k + i + jj ; k1 = kk + kspan ; k2 = 0 ; do /*530*/ { k2 ++ ; wt[k2-1] = a[k1-1] ; bt[k2-1] = b[k1-1] ; k1 -= inc ; }while ( k1 != kk ); do { k1 = kk + kspan ; k2 = k1 - jc * (k + np[k-1]); k = -np[k-1]; do { a[k1-1] = a[k2-1] ; b[k1-1] = b[k2-1] ; k1 -= inc ; k2 -= inc ; }while ( k1 != kk ) ; kk = k2 ; }while ( k != j ); k1 = kk +kspan ; k2 = 0 ; /*560*/ do { k2 ++ ; a[k1-1] = wt[k2-1] ; b[k1-1] = bt[k2-1] ; k1 -= inc ; }while ( k1 != kk ) ; } while ( jj != 0 ) ; }while ( j != 1 ) ; return ; }