added comment on code

8 files changed, 150 insertions, 119 deletions

diff --git a/src/signalProcessing/fft/dfftbi.c b/src/signalProcessing/fft/dfftbi.c
index f22ed995..e98b3bdd 100644
--- a/src/signalProcessing/fft/dfftbi.c
+++ b/src/signalProcessing/fft/dfftbi.c
@@ -15,28 +15,48 @@
 #include "max.h"
 #include "fft_internal.h"
 
-
 /*
-   iIw[0] = 0 ;
-   iIw[1] = 10 ;
-   iIw[2] = 10 ;
-   iIw[3] = lw ;
-   iIw[4] = 10 ;
-
-   iw[0] = 0 ;
-   iw[1] = 10 ;
-   iw[2] = 10 ;
-   iw[3] = lw ;
-   iw[4] = 10 ;
-
-   dfftbi ( a   , b     , nseg  , n      , nspn ,
-            isn , ierr  , iIw[0], iIw[1] , iIw[2],
-            iIw[3], iIw[4], iw, iIw );
- void dfftbi ( double* a , double* b , int nseg , int n , int nspn ,
-              int isn , int ierr, int lout , int lnow ,
-              int lused ,int lmax , int lbook , double* rstak , int* istakk )
+c arrays a and b originally hold the real and imaginary
+c      components of the data, and return the real and
+c      imaginary components of the resulting fourier coefficients.
+c multivariate data is indexed according to the fortran
+c      array element successor function, without limit
+c      on the number of implied multiple subscripts.
+c      the subroutine is called once for each variate.
+c      the calls for a multivariate transform may be in any order.
+c
+c n is the dimension of the current variable.
+c nspn is the spacing of consecutive data values
+c      while indexing the current variable.
+c nseg*n*nspn is the total number of complex data values.
+c the sign of isn determines the sign of the complex
+c      exponential, and the magnitude of isn is normally one.
+c      the magnitude of isn determines the indexing increment for a&b.
+c
+c if fft is called twice, with opposite signs on isn, an
+c      identity transformation is done...calls can be in either order.
+c      the results are scaled by 1/n when the sign of isn is positive.
+c
+c a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
+c is computed by
+c        call fft(a,b,n2*n3,n1,1,-1)
+c        call fft(a,b,n3,n2,n1,-1)
+c        call fft(a,b,1,n3,n1*n2,-1)
+c
+c a single-variate transform of n complex data values is computed by
+c        call fft(a,b,1,n,1,-1)
+c
+c the data may alternatively be stored in a single complex
+c      array a, then the magnitude of isn changed to two to
+c      give the correct indexing increment and a(2) used to
+c      pass the initial address for the sequence of imaginary
+c      values, e.g.
+c
+c
+c array nfac is working storage for factoring n.  the smallest
+c      number exceeding the 15 locations provided is 12,754,584.
+c!
 */
-
 void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
               int isn   , int ierr)
 {
@@ -94,7 +114,7 @@ void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
       }
 
 
-
+/* we search as much 4 in the factor of vector's length as we can */
 
    while ( (k- (int)(k/16)*16 ) == 0 )
       {
@@ -104,7 +124,7 @@ void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
       }
 
 
-
+/* we search all square factor */
 
    do
       {
@@ -124,7 +144,8 @@ void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
 
 
 
-
+/* if the remaining size after all the previous division is less than 4
+   then it's the last factor  */
   if ( k <= 4)
      {
 
@@ -221,7 +242,7 @@ void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
  /*i = ( (istkgt - 1 + nitems) * isize[3] -1) + 3 ;*/
  i = 12 + nitems*2;
 
-
+/* this part is mainly to allocate size for workspace */
 
    istak = (int*) malloc ( sizeof (int) * (unsigned int) i);
 
@@ -269,18 +290,20 @@ c    *********************************************
 
 
 
-   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);
+   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);
 
    k =2 ;
 
    in = 2 ;
-
+/*
    if (!( lbook <= lnow &&  lnow <= lused  ))
       {
          ierr = 3 ;
          return ;
       }
-
+*/
    while ( in > 0)
       {
          if ( lbook > istak[lnow-1] || istak[lnow-1] >=  lnow-1)
diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c
index e9ab1504..60129b1c 100644
--- a/src/signalProcessing/fft/dfftmx.c
+++ b/src/signalProcessing/fft/dfftmx.c
@@ -204,7 +204,10 @@ void preliminaryWork (void)
 
 
 /*40*/
-/* this function is call as many time as the size of the input vector has odd factors */
+/* 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
+ */
 
 int  factorTransform (void)
 {
@@ -306,7 +309,7 @@ switch ( nfac[i-1] )
 
 }
 
-
+/* permutation for square factor of n */
 void permute_stage1 (void)
 {
 
@@ -494,6 +497,7 @@ do /*60*/
 }
 
 
+/* this one is just an optimisation of the factor of 2 transform , we compute more things each turn */
 
 int factorOf4Transform (void)
 {
@@ -527,6 +531,9 @@ int factorOf4Transform (void)
 
 }
 
+/*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 */
+
 void f4t_150 (void)
 {
 
@@ -762,6 +769,9 @@ void factorOf5Transform (void)
 
 }
 
+/* this function is the general case of non factor of 2 factor , the factorof3transform and factorof5trandform are just
+special case of this one */
+
 
 void preFOtherTransform (void)
 {
@@ -1182,7 +1192,7 @@ void nonSqFactor2NormOrder (void)
    return ;
 }
 
-
+/*  here we determine how many permutation cycles we need to do */
 void detPermutCycles (void)
 {
 
diff --git a/src/signalProcessing/fft/fft842.c b/src/signalProcessing/fft/fft842.c
index cb5d97bb..e110ba9c 100644
--- a/src/signalProcessing/fft/fft842.c
+++ b/src/signalProcessing/fft/fft842.c
@@ -159,7 +159,7 @@ void fft842 (doubleComplex* b, int size , int in)
 
 
 
-
+/* this code is two fix a problem of result order which was different from what scilab  give */
        for ( i = 0 ; i < size /2 ; i++)
         {
           temp =  b[i] ;
diff --git a/src/signalProcessing/fft/r4tx.c b/src/signalProcessing/fft/r4tx.c
index f7b6400d..a6a94110 100644
--- a/src/signalProcessing/fft/r4tx.c
+++ b/src/signalProcessing/fft/r4tx.c
@@ -17,6 +17,7 @@
 /*
 ** radix 4 iteration subroutine
 */
+/* this function do in one turn the same computation that do radix 2 in two turns  */
 void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, doubleComplex* c3)
 {
   int kk;
@@ -26,39 +27,23 @@ void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, d
     {
        /* real and imag parts alternate */
 
-
+    /* this first step is strictly equivalent than calling radix 2
+       except that radix would have needed 2 turns to compute what radix4 do in one */
 	temp1 = zadds ( c0[kk] , c2[kk] ) ;
 	temp2 = zdiffs( c0[kk] , c2[kk] ) ;
 	temp3 = zadds ( c1[kk] , c3[kk] ) ;
 	temp4 = zdiffs( c1[kk] , c3[kk] ) ;
-/*
-      r1 = cr0[kk] + cr2[kk];
-      r2 = cr0[kk] - cr2[kk];
-      r3 = cr1[kk] + cr3[kk];
-      r4 = cr1[kk] - cr3[kk];
 
 
-      i1 = ci0[kk] + ci2[kk];
-      i2 = ci0[kk] - ci2[kk];
-      i3 = ci1[kk] + ci3[kk];
-      i4 = ci1[kk] - ci3[kk];
-*/
+    /* strictly equivalent than calling radix2 with the temporary vector , but here also , radix4 do it in one turn
+    instead of two */
 	c0[kk] = zadds ( temp1 , temp3 );
 	c1[kk] = zdiffs( temp1 , temp3 );
 
-/*
-      cr0[kk] = r1 + r3;
-      ci0[kk] = i1 + i3;
-      cr1[kk] = r1 - r3;
-      ci1[kk] = i1 - i3;
-*/
+
 	c2[kk] = DoubleComplex ( zreals ( temp2 ) - zimags( temp4 ) , zimags ( temp2 ) + zreals( temp4 ) );
 	c3[kk] = DoubleComplex ( zreals ( temp2 ) + zimags( temp4 ) , zimags ( temp2 ) - zreals( temp4 ) );
-/*
-      cr2[kk] = r2 - i4;
-      ci2[kk] = i2 + r4;
-      cr3[kk] = r2 + i4;
-      ci3[kk] = i2 - r4;
-*/
+
+
     }
 }
diff --git a/src/signalProcessing/fft/r8tx.c b/src/signalProcessing/fft/r8tx.c
index 8c44488e..bae32639 100644
--- a/src/signalProcessing/fft/r8tx.c
+++ b/src/signalProcessing/fft/r8tx.c
@@ -20,6 +20,9 @@
 /*
 ** radix 8 iteration subroutine
 */
+
+/* this function do in one turn the same computation that do radix 2 in three turns  */
+
 void r8tx ( int nxtlt,int nthpo,int lengt,
             doubleComplex* cc0,doubleComplex* cc1,doubleComplex* cc2,doubleComplex* cc3,
             doubleComplex* cc4,doubleComplex* cc5,doubleComplex* cc6,doubleComplex* cc7)
@@ -65,9 +68,11 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
 
       for(kk=j;kk<nthpo;kk+=lengt)
 	{
-	  ;/* (k-1)*2*/ /* index by twos; re & im alternate */
+	   /* (k-1)*2*/ /* index by twos; re & im alternate */
 
 
+       /*  first turn the same as calling radix 2 with the input vector */
+       /* but radix2 will have do it in three turn , radix8 do it in one */
 		Atemp0 =  zadds ( cc0[kk] , cc4[kk] ) ;
 		Atemp1 =  zadds ( cc1[kk] , cc5[kk] ) ;
 		Atemp2 =  zadds ( cc2[kk] , cc6[kk] ) ;
@@ -79,26 +84,24 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
 		Atemp6 =  zdiffs ( cc2[kk] , cc6[kk] ) ;
 		Atemp7 =  zdiffs ( cc3[kk] , cc7[kk] ) ;
 
-
+        /* second turn the same as  calling radix 2 with the vector transformed by a previous call of radix2 */
+        /* the same here , three turns in one */
 		Btemp0 =  zadds  ( Atemp0 , Atemp2 ) ;
 		Btemp1 =  zadds  ( Atemp1 , Atemp3 ) ;
 		Btemp2 =  zdiffs ( Atemp0 , Atemp2 ) ;
 		Btemp3 =  zdiffs ( Atemp1 , Atemp3 ) ;
 
+        Btemp4 = DoubleComplex ( zreals ( Atemp4 ) - zimags( Atemp6 ) , zimags ( Atemp4 ) + zreals( Atemp6 ) );
+        Btemp5 = DoubleComplex ( zreals ( Atemp5 ) - zimags( Atemp7 ) , zimags ( Atemp5 ) + zreals( Atemp7 ) );
+        Btemp6 = DoubleComplex ( zreals ( Atemp4 ) + zimags( Atemp6 ) , zimags ( Atemp4 ) - zreals( Atemp6 ) );
+        Btemp7 = DoubleComplex ( zreals ( Atemp5 ) + zimags( Atemp7 ) , zimags ( Atemp5 ) - zreals( Atemp7 ) );
 
-
-
-	Btemp4 = DoubleComplex ( zreals ( Atemp4 ) - zimags( Atemp6 ) , zimags ( Atemp4 ) + zreals( Atemp6 ) );
-	Btemp5 = DoubleComplex ( zreals ( Atemp5 ) - zimags( Atemp7 ) , zimags ( Atemp5 ) + zreals( Atemp7 ) );
-	Btemp6 = DoubleComplex ( zreals ( Atemp4 ) + zimags( Atemp6 ) , zimags ( Atemp4 ) - zreals( Atemp6 ) );
-	Btemp7 = DoubleComplex ( zreals ( Atemp5 ) + zimags( Atemp7 ) , zimags ( Atemp5 ) - zreals( Atemp7 ) );
-
-
-	cc0[kk] = zadds ( Btemp0 , Btemp1 );
-
+        /*third turn the same as calling radix 2 with the vector transformed by two previous call of radix2 */
+	    cc0[kk] = zadds ( Btemp0 , Btemp1 );
 
 
 
+        /* if we are not in the first turn */
 
 	  if(j>0)
 	    {
diff --git a/src/signalProcessing/ifft/dfftbi.c b/src/signalProcessing/ifft/dfftbi.c
index f22ed995..1e09df92 100644
--- a/src/signalProcessing/ifft/dfftbi.c
+++ b/src/signalProcessing/ifft/dfftbi.c
@@ -17,24 +17,46 @@
 
 
 /*
-   iIw[0] = 0 ;
-   iIw[1] = 10 ;
-   iIw[2] = 10 ;
-   iIw[3] = lw ;
-   iIw[4] = 10 ;
-
-   iw[0] = 0 ;
-   iw[1] = 10 ;
-   iw[2] = 10 ;
-   iw[3] = lw ;
-   iw[4] = 10 ;
-
-   dfftbi ( a   , b     , nseg  , n      , nspn ,
-            isn , ierr  , iIw[0], iIw[1] , iIw[2],
-            iIw[3], iIw[4], iw, iIw );
- void dfftbi ( double* a , double* b , int nseg , int n , int nspn ,
-              int isn , int ierr, int lout , int lnow ,
-              int lused ,int lmax , int lbook , double* rstak , int* istakk )
+c arrays a and b originally hold the real and imaginary
+c      components of the data, and return the real and
+c      imaginary components of the resulting fourier coefficients.
+c multivariate data is indexed according to the fortran
+c      array element successor function, without limit
+c      on the number of implied multiple subscripts.
+c      the subroutine is called once for each variate.
+c      the calls for a multivariate transform may be in any order.
+c
+c n is the dimension of the current variable.
+c nspn is the spacing of consecutive data values
+c      while indexing the current variable.
+c nseg*n*nspn is the total number of complex data values.
+c the sign of isn determines the sign of the complex
+c      exponential, and the magnitude of isn is normally one.
+c      the magnitude of isn determines the indexing increment for a&b.
+c
+c if fft is called twice, with opposite signs on isn, an
+c      identity transformation is done...calls can be in either order.
+c      the results are scaled by 1/n when the sign of isn is positive.
+c
+c a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
+c is computed by
+c        call fft(a,b,n2*n3,n1,1,-1)
+c        call fft(a,b,n3,n2,n1,-1)
+c        call fft(a,b,1,n3,n1*n2,-1)
+c
+c a single-variate transform of n complex data values is computed by
+c        call fft(a,b,1,n,1,-1)
+c
+c the data may alternatively be stored in a single complex
+c      array a, then the magnitude of isn changed to two to
+c      give the correct indexing increment and a(2) used to
+c      pass the initial address for the sequence of imaginary
+c      values, e.g.
+c
+c
+c array nfac is working storage for factoring n.  the smallest
+c      number exceeding the 15 locations provided is 12,754,584.
+c!
 */
 
 void dfftbi ( double* a , double* b , int nseg , int n    , int nspn  ,
diff --git a/src/signalProcessing/ifft/r4tx.c b/src/signalProcessing/ifft/r4tx.c
index f7b6400d..a6a94110 100644
--- a/src/signalProcessing/ifft/r4tx.c
+++ b/src/signalProcessing/ifft/r4tx.c
@@ -17,6 +17,7 @@
 /*
 ** radix 4 iteration subroutine
 */
+/* this function do in one turn the same computation that do radix 2 in two turns  */
 void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, doubleComplex* c3)
 {
   int kk;
@@ -26,39 +27,23 @@ void r4tx( int nthpo, doubleComplex* c0, doubleComplex* c1, doubleComplex* c2, d
     {
        /* real and imag parts alternate */
 
-
+    /* this first step is strictly equivalent than calling radix 2
+       except that radix would have needed 2 turns to compute what radix4 do in one */
 	temp1 = zadds ( c0[kk] , c2[kk] ) ;
 	temp2 = zdiffs( c0[kk] , c2[kk] ) ;
 	temp3 = zadds ( c1[kk] , c3[kk] ) ;
 	temp4 = zdiffs( c1[kk] , c3[kk] ) ;
-/*
-      r1 = cr0[kk] + cr2[kk];
-      r2 = cr0[kk] - cr2[kk];
-      r3 = cr1[kk] + cr3[kk];
-      r4 = cr1[kk] - cr3[kk];
 
 
-      i1 = ci0[kk] + ci2[kk];
-      i2 = ci0[kk] - ci2[kk];
-      i3 = ci1[kk] + ci3[kk];
-      i4 = ci1[kk] - ci3[kk];
-*/
+    /* strictly equivalent than calling radix2 with the temporary vector , but here also , radix4 do it in one turn
+    instead of two */
 	c0[kk] = zadds ( temp1 , temp3 );
 	c1[kk] = zdiffs( temp1 , temp3 );
 
-/*
-      cr0[kk] = r1 + r3;
-      ci0[kk] = i1 + i3;
-      cr1[kk] = r1 - r3;
-      ci1[kk] = i1 - i3;
-*/
+
 	c2[kk] = DoubleComplex ( zreals ( temp2 ) - zimags( temp4 ) , zimags ( temp2 ) + zreals( temp4 ) );
 	c3[kk] = DoubleComplex ( zreals ( temp2 ) + zimags( temp4 ) , zimags ( temp2 ) - zreals( temp4 ) );
-/*
-      cr2[kk] = r2 - i4;
-      ci2[kk] = i2 + r4;
-      cr3[kk] = r2 + i4;
-      ci3[kk] = i2 - r4;
-*/
+
+
     }
 }
diff --git a/src/signalProcessing/ifft/r8tx.c b/src/signalProcessing/ifft/r8tx.c
index 8c44488e..bae32639 100644
--- a/src/signalProcessing/ifft/r8tx.c
+++ b/src/signalProcessing/ifft/r8tx.c
@@ -20,6 +20,9 @@
 /*
 ** radix 8 iteration subroutine
 */
+
+/* this function do in one turn the same computation that do radix 2 in three turns  */
+
 void r8tx ( int nxtlt,int nthpo,int lengt,
             doubleComplex* cc0,doubleComplex* cc1,doubleComplex* cc2,doubleComplex* cc3,
             doubleComplex* cc4,doubleComplex* cc5,doubleComplex* cc6,doubleComplex* cc7)
@@ -65,9 +68,11 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
 
       for(kk=j;kk<nthpo;kk+=lengt)
 	{
-	  ;/* (k-1)*2*/ /* index by twos; re & im alternate */
+	   /* (k-1)*2*/ /* index by twos; re & im alternate */
 
 
+       /*  first turn the same as calling radix 2 with the input vector */
+       /* but radix2 will have do it in three turn , radix8 do it in one */
 		Atemp0 =  zadds ( cc0[kk] , cc4[kk] ) ;
 		Atemp1 =  zadds ( cc1[kk] , cc5[kk] ) ;
 		Atemp2 =  zadds ( cc2[kk] , cc6[kk] ) ;
@@ -79,26 +84,24 @@ void r8tx ( int nxtlt,int nthpo,int lengt,
 		Atemp6 =  zdiffs ( cc2[kk] , cc6[kk] ) ;
 		Atemp7 =  zdiffs ( cc3[kk] , cc7[kk] ) ;
 
-
+        /* second turn the same as  calling radix 2 with the vector transformed by a previous call of radix2 */
+        /* the same here , three turns in one */
 		Btemp0 =  zadds  ( Atemp0 , Atemp2 ) ;
 		Btemp1 =  zadds  ( Atemp1 , Atemp3 ) ;
 		Btemp2 =  zdiffs ( Atemp0 , Atemp2 ) ;
 		Btemp3 =  zdiffs ( Atemp1 , Atemp3 ) ;
 
+        Btemp4 = DoubleComplex ( zreals ( Atemp4 ) - zimags( Atemp6 ) , zimags ( Atemp4 ) + zreals( Atemp6 ) );
+        Btemp5 = DoubleComplex ( zreals ( Atemp5 ) - zimags( Atemp7 ) , zimags ( Atemp5 ) + zreals( Atemp7 ) );
+        Btemp6 = DoubleComplex ( zreals ( Atemp4 ) + zimags( Atemp6 ) , zimags ( Atemp4 ) - zreals( Atemp6 ) );
+        Btemp7 = DoubleComplex ( zreals ( Atemp5 ) + zimags( Atemp7 ) , zimags ( Atemp5 ) - zreals( Atemp7 ) );
 
-
-
-	Btemp4 = DoubleComplex ( zreals ( Atemp4 ) - zimags( Atemp6 ) , zimags ( Atemp4 ) + zreals( Atemp6 ) );
-	Btemp5 = DoubleComplex ( zreals ( Atemp5 ) - zimags( Atemp7 ) , zimags ( Atemp5 ) + zreals( Atemp7 ) );
-	Btemp6 = DoubleComplex ( zreals ( Atemp4 ) + zimags( Atemp6 ) , zimags ( Atemp4 ) - zreals( Atemp6 ) );
-	Btemp7 = DoubleComplex ( zreals ( Atemp5 ) + zimags( Atemp7 ) , zimags ( Atemp5 ) - zreals( Atemp7 ) );
-
-
-	cc0[kk] = zadds ( Btemp0 , Btemp1 );
-
+        /*third turn the same as calling radix 2 with the vector transformed by two previous call of radix2 */
+	    cc0[kk] = zadds ( Btemp0 , Btemp1 );
 
 
 
+        /* if we are not in the first turn */
 
 	  if(j>0)
 	    {