path: root/2.3-1/src/c/signalProcessing/fft/dfftbi.c

/*
 *  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 <stdlib.h>
#include <stdio.h>
#include "max.h"
#include "fft_internal.h"

/*
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)
{

    double* rstak ;
    int*  istak ;

   int lout = 0 ;
   int lnow = 10;
   int lused= 10;

   int lbook = 10 ;


   int nfac[15] ;
   int i ;
   int in ;
   int j = 3 ;
   int j2 = 3 ;
   int j3 = 3 ;
   int jj = 9;
   int m = 0 ;
   int k ;
   int kt ;
   int kkk ;
   int nspan ;
   int nitems ;
   int ntot ;
   int maxp  = 0;
   int maxf ;
   int itype;
   int istkgt ;


   int nf = abs ( n ) ;

   ierr = 0 ;

   /*determine the factors of n */


   if ( nf == 1)
      return ;

   k = nf ;

   nspan = abs ( nf*nspn ) ;
   ntot  = abs ( nspan*nseg) ;


   if ( isn*ntot == 0 )
      {
      ierr = 1 ;
      return  ;
      }


/* we search as much 4 in the factor of vector's length as we can */

   while ( (k- (int)(k/16)*16 ) == 0 )
      {
         m++;
         nfac[m-1] = 4 ;
         k = k >> 4 ;
      }


/* we search all square factor */

   do
      {
      while ( k%jj == 0 )
         {
            m++;
            nfac[m-1] = j ;
            k /= jj ;

         }

         j+=2;
         jj= j*j ;

      }while ( jj <= k);




/* if the remaining size after all the previous division is less than 4
   then it's the last factor  */
  if ( k <= 4)
     {

      kt = m;
      nfac[m] = k;
      if ( k != 1 )
         m++;
      }
   else
     {
       if ( (k & 3) == 0 )
         {
            m++;
            nfac[m-1] = 2 ;
            k = k >> 2 ;
          }

       /*all square factor out now but k >= 5 still */
       kt = m ;
       maxp = max ( (kt+1)*2 , k-1);
       j=2;

      do
        {
         if ( k%j == 0 )
            {

               m++;
               nfac[m-1] = j ;
               k /= j ;
            }

          j = (j+1) | 1 ;

         }while ( j <= k );

      }



   if ( m <= ( kt+1) )
      maxp = m + kt + 1  ;



   if ( m + kt > 15)
    {
      ierr = 2 ;

      return ;
    }


   if ( kt != 0 )
      {
         j = kt ;

         do{
            m++;

            nfac[m-1] = nfac[j-1];
            j--;
           }while ( j != 0) ;
      }


    maxf = nfac[m-kt-1] ;

   if ( kt > 0 )
      maxf = max ( nfac[kt-1] , maxf );




   for ( kkk = 1 ; kkk <= m ; kkk++ )
      {
      maxf = max ( maxf , nfac[kkk-1]);

      }







 nitems = maxf * 4 ;
 itype = 4 ;


 istkgt = 2 + ((lnow-1)/2) ;/*lnow = 10*/
 istkgt = 6;

 /*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);

   istak[i-2] = itype ;
   istak[i-1] = lnow  ;
   lout ++ ;
   lnow = i ;
   lused = max ( lused , lnow );

   j = istkgt ;
   jj = j + maxf ;
   j2 = jj+ maxf ;
   j3 = j2+ maxf ;

   nitems = maxp ;
   itype  = 2 ;

  /*istkgt = ( lnow*isize[1] -1)/isize[1] + 2;*/
   istkgt =  lnow + 1 ;
   /*i = ( (istkgt - 1 + nitems) * isize[1] -1) / isize[1] + 3 ;*/
   i = lnow  + nitems + 2 ;
   istak = (int*) realloc ( istak ,sizeof (int) * (unsigned int) i);
   rstak = (double*) malloc ( sizeof (double) * (unsigned int) i);






   istak[i-2] = itype ;
   istak[i-1] = lnow  ;
   lout ++ ;
   lnow = i ;
   lused = max ( lused , lnow );

   k = istkgt ;

/*
c     la carte suivante est a supprimer si simple precision
c     next instruction commented by FD&MG (simulog residue?)
c    ********************************************
c      k=2*k-1
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);

   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)
         {
            ierr = 4 ;
         }

         lout-- ;
         lnow = istak[lnow-1] ;
         in-- ;
      }
   free(istak);
   free(rstak);
   
   return ;
}