debug fft and add new tests

1 files changed, 200 insertions, 317 deletions

diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c
index 60129b1c..74cfcb4a 100644
--- a/src/signalProcessing/fft/dfftmx.c
+++ b/src/signalProcessing/fft/dfftmx.c
@@ -97,7 +97,6 @@ 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 ;
@@ -119,41 +118,23 @@ int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan,
    nt = inc*ntot ;
    ks = inc*nspan;
    rad = atan ( 1 );
-
    c72 = cos (rad/0.6250);
    s72 = sin (rad/0.6250);
    s120= sqrt(0.750);
 
 
 
-
-
-
-
    preliminaryWork () ;
 
 
-   while ( retVal == 0 )
-      {
-
-        retVal = factorTransform ( ) ;
-      }
-
+   while ( retVal == 0 ) retVal = factorTransform ( ) ;
 
    np[0] = ks ;
 
-   if ( kt != 0)
-    {
-     permute_stage1 ( ) ;
-    }
-
-   if ( 2*kt + 1 < m )
-    {
-
-      permute_stage2 ( ) ;
-    }
-
+   if ( kt != 0) permute_stage1 ( ) ;
+    
 
+   if ( 2*kt + 1 < m ) permute_stage2 ( ) ;
 
    _pdblA = a ;
    _pdblB = b ;
@@ -170,14 +151,9 @@ Sous-Fonctions
 /* this function only set the value of variable */
 void preliminaryWork (void)
 {
-
-    int lim ;
-
-
-
-       s72 = -s72 ;
-       s120= -s120;
-       rad = -rad ;
+   s72 = -s72 ;
+   s120= -s120;
+   rad = -rad ;
 
    kspan = ks ;
    nn = nt -inc ;
@@ -185,19 +161,12 @@ void preliminaryWork (void)
 
    /* sin , cos values are re-initialized each lim steps */
 
-   lim = 32 ;
-   klim = lim*jc ;
+   klim = 32*jc ;
    i = 0;
    jf = 0 ;
 
-   maxf = m -kt ;
-   maxf = nfac[maxf-1] ;
-
-
-
-   if ( kt > 0 )
-      maxf = max ( nfac[kt-1] , maxf );
-
+   if ( kt > 0 )  maxf = max ( nfac[kt-1] , nfac[m-kt-1] );
+   else 	  maxf = nfac[m-kt-1] ;
 
 }
 
@@ -214,8 +183,6 @@ 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) ;
@@ -227,38 +194,32 @@ int  factorTransform (void)
 
 switch ( nfac[i-1] )
    {
-      case 2 :
+      case 2 :printf("passge en 2\n");
          /*transform for factor of 2 (including rotation factor)*/
 
-         retVal = pre_fOf2Trans( ) ;
-         if ( retVal == 0 )
-            {
-
-            factorOf2Transform ( ) ;
+         retVal = pre_fOf2Trans() ;
+         if ( retVal == 0 )  factorOf2Transform () ;
 
-            }
          break ;
 
-      case 4 :
+      case 4 :printf("passge en 4\n");
          /*transform for factor of 4 */
          kspnn = kspan ;
          kspan = kspan >> 2 ; /*kspan /= 4 */
 
-
-         retVal = factorOf4Transform ( ) ;
+         retVal = factorOf4Transform () ;
          break ;
 
-      case 3 :
-
+      case 3 :printf("passge en 3\n");
+printf("valeurs init : %d and %d and %d\n",nfac[i-1],kspan,nn);
          k = nfac[i-1] ;
          kspnn = kspan ;
          kspan = kspan / k ;
-
-
+         
          factorOf3Transform ( ) ;
          break ;
 
-      case 5 :
+      case 5 :printf("passge en 5\n");
 
          k = nfac[i-1] ;
          kspnn = kspan ;
@@ -267,19 +228,14 @@ switch ( nfac[i-1] )
          factorOf5Transform ( ) ;
          break ;
 
-      default :
+      default :printf("passge en autre\n");
 
          k = nfac[i-1] ;
          kspnn = kspan ;
          kspan = kspan / k ;
 
-
-         if ( nfac[i-1] != jf)
-            {
-
-             preFOtherTransform ( ) ;
-            }
-
+         if ( nfac[i-1] != jf) preFOtherTransform ( ) ;
+         
          factorOfOtherTransform ( ) ;
          break ;
     }
@@ -288,24 +244,12 @@ switch ( nfac[i-1] )
 
    if ( retVal == 42 )
     {
-
-
-      if ( i !=  m)
-        {
-
-         retVal = mulByRotationFactor ( ) ;
-        }
-      else
-        retVal = 1 ;
+      if ( i !=  m)  retVal = mulByRotationFactor ( ) ;
+      else           retVal = 1 ;
     }
 
-
-   if ( retVal == 1 )
-      return 1 ; /*goto permute */
-   else
-      return 0 ; /*goto factor_transform => once again*/
-
-
+   if ( retVal == 1 ) return 1 ; /*goto permute */
+   else 	      return 0 ; /*goto factor_transform => once again*/
 
 }
 
@@ -315,59 +259,45 @@ void permute_stage1 (void)
 
  int retVal = 1 ;
 
-   pre_sqFactor2NormlOrder ( ) ;
+   pre_sqFactor2NormlOrder () ;
 
    if ( n == ntot )
       /*permutation for single-variate transform (optional code)*/
       while ( retVal == 1)
         {
-         single_sqFactor2NormlOrder ( ) ;
+         single_sqFactor2NormlOrder () ;
          retVal = post_sqFactor2NormlOrder () ;
         }
    else
       /*permutation for multivariate transform*/
-      while ( retVal == 1)
-        {
-         retVal = multi_sqFactor2NormlOrder ( );
-        }
-
+      while ( retVal == 1) retVal = multi_sqFactor2NormlOrder ();
 
 }
 
 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 ;
+	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
 ******************************************/
 
@@ -377,42 +307,30 @@ Sous-Sous-Fonctions
 
 int pre_fOf2Trans (void)
 {
-   int ktemp = 0 ;
+	kspan /= 2;
+	k1 = kspan + 2 ;
+	/*50*/
+	do{
+		do{
+			k2 = kk + kspan ;
+			ak = a[k2-1] ;
+			bk = b[k2-1] ;
 
-   kspan /= 2;
-   k1 = kspan + 2 ;
-/*50*/
-   do
-     {
+			a[k2-1] = a[kk-1] - ak;
+			b[k2-1] = b[kk-1] - bk;
 
-      k2 = kk + kspan ;
-      ak = a[k2-1] ;
-      bk = b[k2-1] ;
+			a[kk-1] = a[kk-1] + ak;
+			b[kk-1] = b[kk-1] + bk;
 
-      a[k2-1] = a[kk-1] - ak;
-      b[k2-1] = b[kk-1] - bk;
+			kk = k2 + kspan ;
+		}while (kk <= nn);
 
-      a[kk-1] = a[kk-1] + ak;
-      b[kk-1] = b[kk-1] + bk;
-
-      kk = k2 + kspan ;
-      ktemp = kk ;
-
-      if ( kk > nn )
-         {
-           kk -= nn ;
-         }
-      }while (ktemp <= nn ||( kk <= jc && ktemp <= nn));
-
-
-
-   if ( kk > kspan )
-      return 1 ; /*goto350*/
-   else
-
-      return 0 ;/*goto60*/
+		kk -= nn ;
+	}while (kk <= jc);
 
 
+	if ( kk > kspan )  return 1 ; /*goto350*/
+	else               return 0 ; /*goto60*/
 
 
 }
@@ -421,79 +339,66 @@ int pre_fOf2Trans (void)
 
 int factorOf2Transform (void)
 {
-int ktemp = 0 ;
-
-
-
-
-do /*60*/
-  {
-   c1 = 1 - cd ;
-   s1 = sd ;
-   mm = min( k1/2 , klim);
-
-    do/* do 80 */
-       {
-        do
-            {
-              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;
-              ktemp = kk ;
-              if (kk >= nt)
-                {
-                  k2 = kk - nt;
-                  c1 = -c1;
-                  kk = k1 - k2;
-
-                }
-
-            }while (  ktemp < nt || (kk > k2 &&  ( ktemp >= nt))  );
-
-        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*/
+	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*/
 }
 
 
@@ -501,7 +406,6 @@ do /*60*/
 
 int factorOf4Transform (void)
 {
-
    int return_value = 0 ;
 
    /*120*/
@@ -537,8 +441,6 @@ int factorOf4Transform (void)
 void f4t_150 (void)
 {
 
-   int sign = 1 ;
-
    do{
       k1 = kk + kspan ;
       k2 = k1 + kspan ;
@@ -562,17 +464,11 @@ void f4t_150 (void)
       b[kk-1] = bkp + bjp ;
       bjp = bkp - bjp ;
 
-      if ( isn <  0 )
-         sign =  1 ;
-      else
-         sign = -1 ;
-
+      akp = akm + bjm ;
+      akm = akm - bjm ;
 
-      akp = akm +(sign * bjm );
-      akm = akm -(sign * bjm );
-
-      bkp = bkm +(sign *  ajm) ;
-      bkm = bkm -(sign *  ajm) ;
+      bkp = bkm - ajm ;
+      bkm = bkm + ajm ;
 
       if ( s1 == 0 )/*190*/
          {
@@ -597,14 +493,14 @@ void f4t_150 (void)
          a[k2-1] = bjp*c2 + ajp*s2 ;
          a[k3-1] = bkm*c3 + akm*s3 ;
          }
-    }while  ( kk < nt ) ;
+       kk=k3+kspan;
+    }while  ( kk <= nt ) ;
 
 
 }
 
 int  f4t_170 (void)
 {
-
    kk += ( jc - nt ) ;
 
    if ( kk <= mm )
@@ -620,12 +516,13 @@ int  f4t_170 (void)
 
       c1 = 0.5/(c2*c2+s1*s1) + 0.5 ;
       s1 *= c1 ;
-      c2 *= c1 ;
+      c1 *= c2 ;
 
       /*140*/
 
       c2 = c1*c1 - s1*s1 ;
       s2 = c1*s1*2 ;
+      c3 = c2*c1 - s2*s1 ;
       s3 = c2*s1 + s2*c1 ;
 
 
@@ -637,7 +534,7 @@ int  f4t_170 (void)
       if ( kk <= kspan )
          {
           s1 = dr*rad * (kk-1)/jc ;
-          c2 = cos (s1) ;
+          c1 = cos (s1) ;
           s1 = sin (s1) ;
           mm = min ( kspan , mm  + klim );
 
@@ -645,6 +542,7 @@ int  f4t_170 (void)
 
          c2 = c1*c1 - s1*s1 ;
          s2 = c1*s1*2 ;
+      	 c3 = c2*c1 - s2*s1 ;
          s3 = c2*s1 + s2*c1 ;
 
          return 0 ;
@@ -659,114 +557,99 @@ int  f4t_170 (void)
 
 void factorOf3Transform (void)
 {
-int ktemp = 0 ;
-do
-   {
-
-   k1 = kk + kspan ;
-   k2 = k1 + kspan ;
+	do{
+		do{
+			k1 = kk + kspan ;
+			k2 = k1 + kspan ;
 
-   ak = a[kk-1] ;
-   bk = b[kk-1] ;
+			ak = a[kk-1] ;
+			bk = b[kk-1] ;
 
-   aj = a[k1-1] + a[k2-1] ;
-   bj = b[k1-1] + b[k2-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 ;
+			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 ;
-   ktemp = kk ;
+			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 ;	
 
-   if ( kk >= nn )
-      kk -= nn ;
-
-   }while (  ktemp < nn || (kk <= kspan &&  ( ktemp >= nn))  );
+			kk = k2 + kspan ;
+		} while (kk < nn);
+		
+		kk -= nn ;
+	}while (kk <= kspan);
 
 }
 
 void factorOf5Transform (void)
 {
-    int ktemp ;
-
-    c2 = c72*c72 - s72 *s72 ;
-    s2 = 2 * c72*s72;
-
-   do
-      {
-      k1 = kk + kspan ;
-      k2 = k1 + kspan ;
-      k3 = k2 + kspan ;
-      k4 = k3 + kspan ;
+	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] ;
+			akp = a[k1-1] + a[k4-1] ;
+			akm = a[k1-1] - a[k4-1] ;
 
-      ajp = a[k2-1] + a[k3-1] ;
-      ajm = a[k2-1] - a[k3-1] ;
+			bkp = b[k1-1] + b[k4-1] ;
+			bkm = b[k1-1] - b[k4-1] ;
 
-      bjp = b[k2-1] + b[k3-1] ;
-      bjm = b[k2-1] - b[k3-1] ;
+			ajp = a[k2-1] + a[k3-1] ;
+			ajm = a[k2-1] - a[k3-1] ;
 
-      aa = a[kk-1] ;
-      bb = b[kk-1] ;
+			bjp = b[k2-1] + b[k3-1] ;
+			bjm = b[k2-1] - b[k3-1] ;
 
-      a[kk-1] = aa + akp + ajp;
-      b[kk-1] = bb + bkp + bjp;
+			aa = a[kk-1] ;
+			bb = b[kk-1] ;
 
-      ak = akp*c72 + ajp*c2 + aa ;
-      bk = bkp*c72 + bjp*c2 + bb ;
+			a[kk-1] = aa + akp + ajp;
+			b[kk-1] = bb + bkp + bjp;
 
-      aj = akm*s72 + ajm*s2 ;
-      bj = bkm*s72 + bjm*s2 ;
+			ak = akp*c72 + ajp*c2 + aa ;
+			bk = bkp*c72 + bjp*c2 + bb ;
 
-      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 ;
+			aj = akm*s72 + ajm*s2 ;
+			bj = bkm*s72 + bjm*s2 ;
 
-      bj = bkm*s2 - bjm*s72 ;
+			a[k1-1] = ak - bj ;
+			a[k4-1] = ak + bj ;
+			b[k1-1] = bk + aj ;
+			b[k4-1] = bk - aj ;
 
-      a[k2-1] = ak - bj ;
-      a[k3-1] = ak + bj ;
-      b[k2-1] = bk + aj ;
-      b[k3-1] = bk - aj ;
+			ak = akp*c2 + ajp*c72 + aa ;
+			bk = bkp*c2 + bjp*c72 + bb ;
 
-      kk = k4 + kspan;
-      ktemp = kk ;
+			aj = akm*s2 - ajm*s72 ;
 
+			bj = bkm*s2 - bjm*s72 ;
 
-      if ( kk >= nn )
-        kk -= nn ;
-
-
-   }while (ktemp < nn || (  kk <= kspan && ktemp >= nn));
-
+			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
@@ -775,7 +658,7 @@ special case of this one */
 
 void preFOtherTransform (void)
 {
-
+printf("0.k=%d \n",k);
 
    jf = k ;
    s1 = (rad*8)/k ;