diff options
author | simon | 2008-10-15 09:08:28 +0000 |
---|---|---|
committer | simon | 2008-10-15 09:08:28 +0000 |
commit | 1e049292f41e3487dfc20ede6ab58a8241c3fb36 (patch) | |
tree | 4cefe40fefe612f18987195dd43cec8bade295c8 /src/signalProcessing/fft | |
parent | 8d2aba552a523cd95585bc492b75e9b70438f922 (diff) | |
download | scilab2c-1e049292f41e3487dfc20ede6ab58a8241c3fb36.tar.gz scilab2c-1e049292f41e3487dfc20ede6ab58a8241c3fb36.tar.bz2 scilab2c-1e049292f41e3487dfc20ede6ab58a8241c3fb36.zip |
fft work with an 7-elements' vector , still not working with a 6-elements
Diffstat (limited to 'src/signalProcessing/fft')
-rw-r--r-- | src/signalProcessing/fft/dfftbi.c | 5 | ||||
-rw-r--r-- | src/signalProcessing/fft/dfftmx.c | 132 | ||||
-rw-r--r-- | src/signalProcessing/fft/testDoubleFft.c | 77 |
3 files changed, 150 insertions, 64 deletions
diff --git a/src/signalProcessing/fft/dfftbi.c b/src/signalProcessing/fft/dfftbi.c index c840903e..d9bf1bd6 100644 --- a/src/signalProcessing/fft/dfftbi.c +++ b/src/signalProcessing/fft/dfftbi.c @@ -284,8 +284,9 @@ c ********************************************* , ntot , nf , nspan , isn , m , kt , j ,jj, j2,j3 , k ); */ for (i = 0 ; i < 15 ; i++ ) - //printf ( "\t nfac[%d] = %d\n" , i , nfac[i]) ; - + { + printf ( "\t nfac[%d] = %d\n" , i , nfac[i]) ; + } dfftmx( a , b , ntot , nf , nspan , isn , m , kt , &rstak[j-1] , &rstak[jj-1] , &rstak[j2-1] , &rstak[j3-1] , &istak[k-1] , nfac); diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c index 26c995a3..16514578 100644 --- a/src/signalProcessing/fft/dfftmx.c +++ b/src/signalProcessing/fft/dfftmx.c @@ -124,21 +124,21 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, s120= sqrt(0.750); int iii = 0 ; - //fprintf (stderr , "\n\n" ); + fprintf (stderr , "\n\n" ); for ( iii = 0 ; iii < 3 ; iii++) { - //fprintf (stderr , "\t\t %d tot : %f \t %f\n" , iii ,a[iii], b[iii]); + fprintf (stderr , "\t\t %d tot : %f \t %f\n" , iii ,a[iii], b[iii]); } - //fprintf (stderr , "preliminary\n" ); + fprintf (stderr , "preliminary\n" ); preliminaryWork () ; while ( retVal == 0 ) { - //fprintf (stderr , "factortransform\n" ); + fprintf (stderr , "factortransform\n" ); retVal = factorTransform ( ) ; } @@ -147,13 +147,13 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, if ( kt != 0) { - //fprintf (stderr , "permute stage 1\n" ); + fprintf (stderr , "permute stage 1\n" ); permute_stage1 ( ) ; } if ( 2*kt + 1 < m ) { - //fprintf (stderr , "permute stage 2\n" ); + fprintf (stderr , "permute stage 2\n" ); permute_stage2 ( ) ; } @@ -162,7 +162,7 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, for ( iii = 0 ; iii < 3 ; iii++) { - //fprintf (stderr , "\t\t out %d tot : %f \t %f\n" , iii ,a[iii], b[iii]); + fprintf (stderr , "\t\t out %d tot : %f \t %f\n" , iii ,a[iii], b[iii]); } @@ -216,7 +216,7 @@ void preliminaryWork (void) maxf = m -kt ; maxf = nfac[maxf-1] ; - //fprintf (stderr , "\tm - kt %d , nfac %d kt %d\n" , m - kt , maxf , kt); + fprintf (stderr , "\tm - kt %d , nfac %d kt %d\n" , m - kt , maxf , kt); if ( kt > 0 ) maxf = max ( nfac[kt-1] , maxf ); @@ -237,16 +237,16 @@ int factorTransform (void) sd = sin(dr*rad) ; kk = 1 ; i++ ; - //fprintf (stderr , "avant switch i %d ,nfac[i-1] %d\n" , i , nfac[i-1]); + fprintf (stderr , "avant switch i %d ,nfac[i-1] %d\n" , i , nfac[i-1]); switch ( nfac[i-1] ) { case 2 : /*transform for factor of 2 (including rotation factor)*/ - //fprintf (stderr , "\tpre_fOf2Trans\n" ); + fprintf (stderr , "\tpre_fOf2Trans\n" ); retVal = pre_fOf2Trans( ) ; if ( retVal == 0 ) { - //fprintf (stderr , "\tfactorOf2Transform\n" ); + fprintf (stderr , "\tfactorOf2Transform\n" ); factorOf2Transform ( ) ; } @@ -257,26 +257,26 @@ switch ( nfac[i-1] ) kspnn = kspan ; kspan = kspan >> 2 ; /*kspan /= 4 */ - //fprintf (stderr , "\tfactorOf4Transform\n" ); + fprintf (stderr , "\tfactorOf4Transform\n" ); retVal = factorOf4Transform ( ) ; break ; case 3 : - //fprintf (stderr , "\tfactorOf3Transform\n" ); + fprintf (stderr , "\tfactorOf3Transform\n" ); k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; - //fprintf (stderr , "\t\t k %d i %d\n" , k ,i); + fprintf (stderr , "\t\t k %d i %d\n" , k ,i); factorOf3Transform ( ) ; break ; case 5 : - //fprintf (stderr , "\tfactorOf5Transform\n" ); + fprintf (stderr , "\tfactorOf5Transform\n" ); k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; - //fprintf (stderr , "\t\t k %d\n" , k); + fprintf (stderr , "\t\t k %d\n" , k); factorOf5Transform ( ) ; break ; @@ -285,14 +285,14 @@ switch ( nfac[i-1] ) k = nfac[i-1] ; kspnn = kspan ; kspan = kspan / k ; - //fprintf (stderr , "\t\t k %d\n" , k); - //fprintf (stderr , "\t\t nfac[i-1] %d jf %d\n" , nfac[i-1] , jf ) ; + fprintf (stderr , "\t\t k %d\n" , k); + fprintf (stderr , "\t\t nfac[i-1] %d jf %d\n" , nfac[i-1] , jf ) ; if ( nfac[i-1] != jf) { - //fprintf (stderr , "\tpreFOtherTransform \n" ); + fprintf (stderr , "\tpreFOtherTransform \n" ); preFOtherTransform ( ) ; } - //fprintf (stderr , "\tfactorOfOtherTransform \n" ); + fprintf (stderr , "\tfactorOfOtherTransform \n" ); factorOfOtherTransform ( ) ; break ; } @@ -301,10 +301,10 @@ switch ( nfac[i-1] ) if ( retVal == 42 ) { - //fprintf (stderr , "\t\t i %d m %d\n" , i , m ); + fprintf (stderr , "\t\t i %d m %d\n" , i , m ); if ( i != m) { - //fprintf (stderr , "\tmulByRotationFactor \n" ); + fprintf (stderr , "\tmulByRotationFactor \n" ); retVal = mulByRotationFactor ( ) ; } else @@ -348,22 +348,22 @@ void permute_stage2 (void) int retVal ; kspnn = np[kt] ; - //fprintf (stderr , "\t kspnn %d , kt %d\n" , kspnn , kt); + fprintf (stderr , "\t kspnn %d , kt %d\n" , kspnn , kt); /*permutation for square-free facotrs of n */ - //fprintf (stderr , "\tnonSqFactor2NormOrder 2\n" ); + fprintf (stderr , "\tnonSqFactor2NormOrder 2\n" ); nonSqFactor2NormOrder ( ) ; /*determine the permutation cycles of length greater than 1*/ - //fprintf (stderr , "\tdetPermutCycles 2\n" ); + fprintf (stderr , "\tdetPermutCycles 2\n" ); detPermutCycles ( ); - //fprintf (stderr , "\tend ploplpop 2\n" ); + fprintf (stderr , "\tend ploplpop 2\n" ); retVal = end ( ) ; - //fprintf (stderr , "\t on n'est plus dans le end\n "); + fprintf (stderr , "\t on n'est plus dans le end\n "); while ( retVal == 1) { - //fprintf (stderr , "\reorderMatrix 2\n" ); + fprintf (stderr , "\reorderMatrix 2\n" ); reorderMatrix ( ) ; retVal = end ( ) ; } @@ -386,7 +386,7 @@ int pre_fOf2Trans (void) /*50*/ do { - //fprintf (stderr , "\t 50 \n"); + fprintf (stderr , "\t 50 \n"); k2 = kk + kspan ; ak = a[k2-1] ; bk = b[k2-1] ; @@ -658,7 +658,7 @@ void factorOf3Transform (void) int ktemp = 0 ; do { - //fprintf (stderr , "\t\t une boucle dans factor of 3\n"); + fprintf (stderr , "\t\t une boucle dans factor of 3\n"); k1 = kk + kspan ; k2 = k1 + kspan ; @@ -685,11 +685,11 @@ do kk = k2 + kspan ; ktemp = kk ; - //fprintf (stderr , "\t\t kk %d , nn %d , kspan %d\n" , kk , nn , kspan ); + fprintf (stderr , "\t\t kk %d , nn %d , kspan %d\n" , kk , nn , kspan ); if ( kk >= nn ) kk -= nn ; - //fprintf (stderr , "\t\t 2kk %d , nn %d , kspan %d\n" , kk , nn , kspan ); + fprintf (stderr , "\t\t 2kk %d , nn %d , kspan %d\n" , kk , nn , kspan ); }while ( ktemp < nn || (kk <= kspan && ( ktemp >= nn)) ); } @@ -698,7 +698,7 @@ void factorOf5Transform (void) { c2 = c72*c72 - s72 *s72 ; s2 = 2 * c72*s72; - //fprintf (stderr , "plop\n" ) ; + fprintf (stderr , "plop\n" ) ; do { k1 = kk + kspan ; @@ -706,7 +706,7 @@ void factorOf5Transform (void) k3 = k2 + kspan ; k4 = k3 + kspan ; - //fprintf (stderr , "kk %d \t k1 %d \t k2 %d \t k3 %d \t k4 %d\n", kk , k1 , k2 ,k3,k4 ); + fprintf (stderr , "kk %d \t k1 %d \t k2 %d \t k3 %d \t k4 %d\n", kk , k1 , k2 ,k3,k4 ); akp = a[k1-1] + a[k4-1] ; akm = a[k1-1] - a[k4-1] ; @@ -741,7 +741,7 @@ void factorOf5Transform (void) bk = bkp*c2 + bjp*c72 + bb ; aj = akm*s2 - ajm*s72 ; - //fprintf (stderr , "aj %f \takm %f \tajm %f\n" , aj , akm , ajm ); + fprintf (stderr , "aj %f \takm %f \tajm %f\n" , aj , akm , ajm ); bj = bkm*s2 - bjm*s72 ; a[k2-1] = ak - bj ; @@ -750,15 +750,15 @@ void factorOf5Transform (void) b[k3-1] = bk - aj ; kk = k4 + kspan; - //fprintf (stderr , "ak %f \tbk %f\naj %f \tbj %f\n" , ak , bk , aj ,bj ); - //fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); + fprintf (stderr , "ak %f \tbk %f\naj %f \tbj %f\n" , ak , bk , aj ,bj ); + fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); if ( kk >= nn ) kk -= nn ; - //fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); + fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); }while ( kk <= kspan && kk < nn); -//fprintf (stderr , "fin 5\n" ); +fprintf (stderr , "fin 5\n" ); } @@ -770,6 +770,7 @@ 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 ; @@ -777,7 +778,7 @@ void preFOtherTransform (void) do { ck[j-1] = ck[k-1] * c1 + sk[k-1]*s1 ; - sk[j-1] = ck[k-1] * c1 + sk[k-1]*s1 ; + sk[j-1] = ck[k-1] * s1 - sk[k-1]*c1 ; k -- ; @@ -785,12 +786,14 @@ void preFOtherTransform (void) sk[k-1] = - sk[j-1] ; j++ ; + fprintf (stderr , "je boucle\n" ); }while ( j < k ); } void factorOfOtherTransform (void) { +int ktemp = 0 ; do { @@ -801,6 +804,7 @@ do bb = b[kk-1] ; ak = aa ; + bk = bb ; j = 1 ; k1 += kspan ; @@ -811,7 +815,10 @@ do wt[j-1] = a[k1-1] + a[k2-1] ; ak = wt[j-1] + ak ; - bt[j-1] = b[k1-1] + bk ; + + bt[j-1] = b[k1-1] + b[k2-1] ; + bk = bt[j-1] + bk ; + printf ("\t 260 ak %f bk %f\n" , ak , bk ) ; j++ ; wt[j-1] = a[k1-1] - a[k2-1] ; @@ -843,16 +850,18 @@ do k++ ; ak += ( wt[k-1] * ck[jj-1] ) ; bk += ( bt[k-1] * ck[jj-1] ) ; + printf ("\t 280 ak %f bk %f\n" , ak , bk ) ; k++ ; - aj += wt[k-1] * sk[jj] ; - bj += bt[k-1] * sk[jj] ; + aj += (wt[k-1] * sk[jj-1]) ; + bj += (bt[k-1] * sk[jj-1]) ; + printf ("\t aj %f bj %f\n" , aj , bj ) ; jj += j ; if ( jj > jf ) jj -= jf ; - } while ( k < jf ) ; + k = jf - j ; a[k1-1] = ak - bj ; b[k1-1] = bk + aj ; @@ -860,15 +869,22 @@ do b[k2-1] = bk - aj ; j++ ; + fprintf (stderr , "je boucle3\n" ); }while ( j < k ) ; + + + + kk += kspnn ; + ktemp = kk ; + fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); if ( kk > nn ) { - kk -= n; + kk -= nn; } -}while ( kk < kspan || kk <= nn) ; - + fprintf (stderr , "kk %d \t nn %d \t kspan %d\n", kk , nn , kspan ); +}while ( ktemp <= nn || (kk <= kspan && ( ktemp > nn)) ); } @@ -1108,7 +1124,7 @@ void nonSqFactor2NormOrder (void) do { nfac[j-1] *= nfac[j] ; - //fprintf (stderr , "\t\t m %d j %d , nfac[j-1] %d\n", m , j , nfac[j-1] ); + fprintf (stderr , "\t\t m %d j %d , nfac[j-1] %d\n", m , j , nfac[j-1] ); j-- ; @@ -1116,7 +1132,7 @@ void nonSqFactor2NormOrder (void) kt ++ ; nn = nfac[kt-1] - 1; - //fprintf (stderr , "\t\t nn %d \n" , nn ) ; + fprintf (stderr , "\t\t nn %d \n" , nn ) ; jj = 0 ; j = 0; @@ -1126,11 +1142,11 @@ void nonSqFactor2NormOrder (void) k = kt + 1 ; kk = nfac[k-1] ; j ++ ; - //fprintf (stderr , "\t\tj %d\n" , j ) ; + fprintf (stderr , "\t\tj %d\n" , j ) ; while ( j <= nn ) { jj += kk ; - //fprintf (stderr , "\t\t 1jj %d , kk %d\n" , jj , kk ) ; + fprintf (stderr , "\t\t 1jj %d , kk %d\n" , jj , kk ) ; while ( jj >= k2 ) { @@ -1140,7 +1156,7 @@ void nonSqFactor2NormOrder (void) kk = nfac[k-1] ; jj += kk ; - //fprintf (stderr , "\t\t jj %d , kk %d\n" , jj , kk ) ; + fprintf (stderr , "\t\t jj %d , kk %d\n" , jj , kk ) ; } @@ -1149,7 +1165,7 @@ void nonSqFactor2NormOrder (void) k = kt + 1 ; kk = nfac[k-1] ; j ++ ; - //fprintf (stderr , "\t\t1j %d\n" , j ) ; + fprintf (stderr , "\t\t1j %d\n" , j ) ; } j = 0 ; @@ -1165,17 +1181,17 @@ void detPermutCycles (void) { do { - //fprintf (stderr , "\t\t j %d \tnp[j-1] %d\n" , j , np[j-1]); + fprintf (stderr , "\t\t j %d \tnp[j-1] %d\n" , j , np[j-1]); j++ ; kk = np[j-1] ; }while ( kk < 0 ) ; - //fprintf (stderr , "\t\t kk %d\t j %d\t\n" , kk , j ); + fprintf (stderr , "\t\t kk %d\t j %d\t\n" , kk , j ); if ( kk != j ) { do { - //fprintf (stderr , "\t\t 2boucle kk %d\n" , kk ); + fprintf (stderr , "\t\t 2boucle kk %d\n" , kk ); k = kk ; kk = np[k-1] ; np[k-1] = -kk ; @@ -1194,13 +1210,13 @@ int end (void) { - //fprintf (stderr , "\t\t\t end of the final end\n"); + fprintf (stderr , "\t\t\t end of the final end\n"); /* j = k3 + 1 ; nt -= kspnn ; i = nt - inc + 1 ; - //fprintf (stderr , "\t\t\t end of the final end\n"); + fprintf (stderr , "\t\t\t end of the final end\n"); if ( nt >= 0 ) return 1 ; diff --git a/src/signalProcessing/fft/testDoubleFft.c b/src/signalProcessing/fft/testDoubleFft.c index 6927594e..9bb41ce4 100644 --- a/src/signalProcessing/fft/testDoubleFft.c +++ b/src/signalProcessing/fft/testDoubleFft.c @@ -21,7 +21,10 @@ #define COLS4 4 #define COLS5 5 #define COLS6 6 +#define COLS7 7 #define COLS8 8 +#define COLS9 9 +#define COLS10 10 #define COLS16 16 #define COLS32 32 @@ -45,6 +48,12 @@ #define ZIMAG_IN6 { 0.21460078610107303, 0.31264199689030647, 0.36163610080257058, 0.2922266637906432,\ 0.56642488157376647, 0.59350947011262178 } +#define ZREAL_IN7 { 0.54425731627270579, 0.23207478970289230, 0.23122371966019273, 0.21646326314657927,\ + 0.65251349471509457, 0.88338878145441413, 0.30760907428339124 } +#define ZIMAG_IN7 { 0.21460078610107303, 0.31264199689030647, 0.36163610080257058, 0.2922266637906432,\ + 0.40948254754766822, 0.56642488157376647, 0.59350947011262178 } + + #define ZREAL_IN8 { 0.54425731627270579, 0.23207478970289230, 0.23122371966019273, 0.21646326314657927,\ 0.88338878145441413, 0.65251349471509457, 0.30760907428339124, 0.93296162132173777 } @@ -95,6 +104,13 @@ 0.69993670814631914} +#define ZREAL_RESULT7 { 3.06753043923527002,-0.62032167153569062,-0.13156333379499591, 0.48353341667797933,\ + 0.63567251139259018, 0.05503001802946385, 0.31991983390432432} +#define ZIMAG_RESULT7 { 2.75052244681864977, 0.82490994311348309,-0.93592353228518299,-0.23131444371235776,\ + -0.12732936894919694, 0.16455873200809046,-0.94321827428597393} + + + #define ZREAL_RESULT8 { 4.00049206055700779,-0.43357241280891956, 0.79836636409163475,-0.91119240848798977,\ -0.06753427721560001,-0.18576209864995416, 0.97926024347543716, 0.17400105922003017} #define ZIMAG_RESULT8 { 3.15585898794233799, 0.62132445165622818, 0.35205427557229996, 0.28289917172258683,\ @@ -338,6 +354,59 @@ static void zfftmaTest6 (void ) } + +static void zfftmaTest7 (void ) +{ + int i = 0 ; + + double tRealIn [] = ZREAL_IN7; + double tImagIn [] = ZIMAG_IN7 ; + + + + double tRealResult [] = ZREAL_RESULT7; + double tImagResult [] = ZIMAG_RESULT7; + + + + doubleComplex* out = (doubleComplex*) malloc ( sizeof(doubleComplex) * (unsigned int) (ROW*COLS7)); + doubleComplex* in = DoubleComplexMatrix ( tRealIn , tImagIn , ROW*COLS7 ); + doubleComplex* Result = DoubleComplexMatrix ( tRealResult , tImagResult ,ROW*COLS7) ; + + + + zfftma ( in , ROW , COLS7 , out ) ; + + /* if we don't add that test assert failed if result = 0 'cause then we have |(out - 0)|/|out| = 1*/ + + + for ( i = 0 ; i < (ROW*COLS7 ) ; i++ ) + { + printf ( "\t\t %d out : %e\t %e\t * i result : %e\t %e\t * i assert : : %e\t %e\t * i \n" , + i , + zreals(out[i]) , + zimags(out[i]), + zreals (Result[i]) , + zimags (Result[i]), + fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) , + fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i]))); + + if ( zreals(out[i]) < 1e-14 && zreals (Result[i]) < 1e-18 ) + assert ( 1 ) ; + else + assert ( fabs( zreals(out[i]) - zreals (Result[i]) ) / fabs (zreals (out[i])) < 1e-12 ); + + + if ( zimags(out[i]) < 1e-14 && zimags (Result[i]) < 1e-18 ) + assert ( 1 ) ; + else + assert ( fabs( zimags(out[i]) - zimags (Result[i]) ) / fabs (zimags (out[i])) < 1e-12 ) ; + + } + + +} + static void zfftmaTest4 (void ) { int i = 0 ; @@ -545,7 +614,7 @@ static int testFft(void) { printf("\n\n\n"); -/* + printf("\t>>>> Vector 2 Double Complex Tests\n"); zfftmaTest2(); printf("\t>>>> Vector 3 Double Complex Tests\n"); @@ -554,9 +623,9 @@ static int testFft(void) { zfftmaTest4(); printf("\t>>>> Vector 5 Double Complex Tests\n"); zfftmaTest5(); -*/ - printf("\t>>>> Vector 6 Double Complex Tests\n"); - zfftmaTest6(); + + printf("\t>>>> Vector 7 Double Complex Tests\n"); + zfftmaTest7(); /* printf("\t>>>> Vector 8 Double Complex Tests\n"); zfftmaTest8(); |