 /*
  * Khoros: $Id: ldacorr.c,v 1.2 1992/03/20 23:24:59 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldacorr.c,v 1.2 1992/03/20 23:24:59 dkhoros Exp $";
#endif

 /*
  * $Log: ldacorr.c,v $
 * Revision 1.2  1992/03/20  23:24:59  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldacorr.c
 >>>>
 >>>>      Program Name: dacorr
 >>>>
 >>>> Date Last Updated: Mon Mar  9 20:14:34 1992 
 >>>>
 >>>>          Routines: ldacorr - the library call for dacorr
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* the cleanup routines free up all allocated memory prior to a return */
#define LDA_CLEANUP { \
                     if(data!=NULL){ \
                         for(i=0;i<num_vects;i++) \
                             if(data[i]!=NULL)free(data[i]); \
                         free(data); \
                      }} 
#define LDF_CLEANUP {if(dr!=NULL)free(dr);if(di!=NULL)free(di);}
#define LDT_CLEANUP {if(autoi!=NULL)free(autoi);if(autor!=NULL)free(autor);}

/* the following two defines implement a log2() function */ 
#define MAGIC 3.32192810
#define LOG2(x) (log10((double)x)*(double)MAGIC)

#define FFT_FORWARD 0
#define FFT_INVERSE 1

static void _lfacorr();

/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldacorr - library call for dacorr
*
* Purpose:
*    
*    Autocorrelation sequence of 1-D signal(s)
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    procdir        process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

* Output:
*    
*    image          pointer to VIFF structure  containing  image  data
*                   after processing.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Ramiro Jordan, Jeremy Worley
*    
*    Jeremy Worley 7 Jun 1990 14:31 MST
*              added error checking to all memory allocations
*    
*    Jeremy Worley 19 Jul 19 1990 16:27 MST
*              changed apoints to dpoints in minor  loop  of  ldtacorr
*              for  both  real  and complex parts to conform to proper
*              definition of autocorrelation.
*    
*    Jeremy Worley 17 Feb 1992 10:16 MST
*              changed call to lfft1d in ldfacorr() to fft842_() since
*              lfft1d  is  being  phased  out.   Also  cleaned up some
*              implicitly defined functions and modularized ldfacorr()
*              which  resulted  in  the  creation of a static function
*              called _lfacorr(). Also fixed a small  scaling  problem
*              with the output of ldfacorr().
*    
*    Jeremy Worley 09 Mar 1992 19:54 MST
*              Deleted a declaration to cdiv() since it was not used.
*    
*    

****************************************************************/


/* -library_def */
int
ldacorr ( image, arith_type, out_option, points, process_dir )

int out_option, arith_type, points, process_dir;
struct xvimage *image;
/* -library_def_end */

/* -library_code */
{
    int   i, data_type, dimension, num_vects, apoints;
    int   dunload_vector(), lscale(), ldfacorr(), ldtacorr();
    float **data;
    char  **dload_vector();

    char *program = "ldacorr";

    /* initialize parameters */
    data_type = 0;
    apoints = points;

    if ( image->data_storage_type == VFF_TYP_COMPLEX )
    {
       data_type = 1;
    }

    if((data = (float **)dload_vector(image, &num_vects, &dimension, process_dir)) == NULL)
    {
      (void) fprintf(stderr,"%s: dload_vector failed \n",program);
      return(0);
    }

    if ( apoints > dimension )
    {
          (void) fprintf (stderr, "\n\n%s:   ", program);
          (void) fprintf (stderr, "%s: number of autocorrelation points ",
                          program);
          (void) fprintf (stderr, "cannot be greater than number of data points\n\n");
          LDA_CLEANUP;
          return (0);
    }

    /* find autocorrelation sequence */
    if ( out_option == 2 )              
    {
       for ( i = 0; i < num_vects; i++ )
       {
           if ( !ldfacorr ( data[i], dimension, apoints, data_type, arith_type ) )
           {
              (void) fprintf ( stderr,"ldfacorr failed\n");
              LDA_CLEANUP;
              return(0);
           }
       }
    }
    else 
    {
       for ( i = 0; i < num_vects; i++ )
       {
           if ( !ldtacorr ( data[i], dimension, apoints, data_type, out_option, arith_type ) )
           {
              (void) fprintf ( stderr,"ldtacorr failed\n");
              LDA_CLEANUP;
              return(0);
           }
       }
    }

    if (!dunload_vector((char **)data, image, 
        (unsigned long)image->data_storage_type, num_vects, apoints, 
        process_dir))
    {
       (void) fprintf (stderr,"%s: dunload_vector failed \n",program);
       LDA_CLEANUP;
       return(0);
    }

    LDA_CLEANUP;
    return(1);

} /* end of ldacorr */


/*
*****************************************************************************
*
*       datain:         input - input data points
*                       output - autocorrelation sequence
*
*       apoints:        number of points of desired autocorrelation
*                       sequence
*
*       dpoints:        number of points of input data sequence 
*
*       data_type:      0 - real data
*                       1 - complex data
*
*       out_option:     0 - biased autocorrelation estimate
*                       1 - unbiased autocorrelation estimate
*
*       arith_type:     0 - scalar arithmetic 
*                       1 - vector arithmetic
*
*****************************************************************************
*/

int
ldtacorr ( datain, dpoints, apoints, data_type, out_option, arith_type )

int out_option, arith_type, apoints, dpoints, data_type;
float *datain;
{
    int i, j, lag, length, size;    
    int   cadd(), cmul();
    float *autor, *autoi, real, imag, tempr, tempi;
    char  *malloc();
    char *program = "ldtacorr";

    /* initialize parameters */
    size = dpoints;

    if ( apoints > dpoints )
    {
          (void) fprintf (stderr, "\n\n%s:   ", program);
          (void) fprintf (stderr, "ldtacorr: number of autocorrelation points ");
          (void) fprintf (stderr, "cannot be greater than number of data points\n\n");
          return (0);
    }


   /* allocate space for working arrays */
    autor = (float *) malloc((unsigned)(apoints * sizeof(float)));
    if(autor==NULL){
       fprintf(stderr,"%s: [1]memory allocation failed.\n",program);
       return(0);
    }

       autoi = (float *) malloc((unsigned)(apoints * sizeof(float)));
       if(autoi==NULL){
          fprintf(stderr,"%s: [2]memory allocation failed.\n",program);
          LDT_CLEANUP;
          return(0);
       }

    /* find the autocorrelation estimate */
    if ( !data_type )                   /* real data */
    {
       for ( lag = 0; lag < apoints; lag++ )
       {
           tempr = 0.0;
           for ( j = 0; j < (dpoints-lag); j++ )
           {
               tempr = tempr + datain[lag+j] * datain[j];
           }
           if ( !out_option )
           {
              autor[lag] = tempr / (float) (size);
           }
           else if ( out_option == 1 )
           {
              autor[lag] = tempr / (float) (size-lag);
           }
       }
       bcopy ( (char *)autor, (char *)datain, (int)apoints*sizeof(float) );
    }
    else if ( data_type == 1 )          /* complex data */
    {
       if ( !arith_type )               /* scalar arithmetic */
       {
          length = apoints * 2;
          for ( i = 0, lag = 0; lag < length; lag += 2, i++ )
          {
              tempr = 0.0;
              tempi = 0.0;
              for ( j = 0; j < (dpoints*2-lag); j += 2 )
              {
                  tempr = tempr + datain[lag+j] * datain[j];
                  tempi = tempi + datain[lag+j+1] * datain[j+1];
              }
              if ( !out_option )
              {
                 autor[i] = tempr / (float) (size);
                 autoi[i] = tempi / (float) (size);
              }
              else if ( out_option == 1 )
              {
                 autor[i] = tempr / (float) (size-i);
                 autoi[i] = tempi / (float) (size-i);
              }
          }
          for ( i = 0, lag = 0; lag < length; lag += 2, i++ )
          {
              datain[lag] = autor[i];
              datain[lag+1] = autoi[i];
          }
       }
       else if ( arith_type == 1 )      /* vector arithmetic */
       {
          length = apoints * 2;
          for ( i = 0, lag = 0; lag < length; lag += 2, i++ )
          {
              real = 0.0;
              imag = 0.0;
              tempr = 0.0;
              tempi = 0.0;
              for ( j = 0; j < (length-lag); j += 2 )
              {
                  cmul ( &tempr, &tempi, datain[j], datain[j+1],
                         datain[lag+j], -datain[lag+j+1] );
                  cadd ( &real, &imag, real, imag, tempr, tempi );

              }
              if ( !out_option )
              {
                 autor[i] = real / (float) (size);
                 autoi[i] = imag / (float) (size);
              }
              else if ( out_option == 1 )
              {
                 autor[i] = real / (float) (size-i);
                 autoi[i] = imag / (float) (size-i);
              }
          }
          for ( i = 0, lag = 0; lag < length; lag += 2, i++ )
          {
              datain[lag] = autor[i];
              datain[lag+1] = autoi[i];
          }
       }
    }

    LDT_CLEANUP;
    return(1);

} /* end of ldtacorr */



/*
*****************************************************************************
*
*       datain:         input - input data points
*                       output - autocorrelation sequence
*
*       apoints:        number of points of desired autocorrelation
*                       sequence
*
*       dpoints:        number of points of input data sequence 
*
*       data_type:      0 - real data
*                       1 - complex data
*
*       out_option:     0 - biased autocorrelation estimate
*                       1 - unbiased autocorrelation estimate
*
*       arith_type:     0 - scalar arithmetic 
*                       1 - vector arithmetic
*
*****************************************************************************
*/

int ldfacorr ( datain, dpoints, apoints, data_type, arith_type )
   int arith_type, apoints, dpoints, data_type;
   float *datain;
{
   int i, j, length, dpoints2;    
   int powtwo();
   int cmul();
   float *dr, *di;

   char *program = "ldfacorr";

/* 
 * initialize parameters 
 */

    length = dpoints + apoints; /* necessary size of data into fft */
    dpoints2 = dpoints * 2;

/* 
 * force length to be a power of two 
 */

    if (!powtwo(length)){
      int temp;
      temp = (int)(LOG2(length)+1);
      length = (int) (pow((double)2,(double)temp));
    }

    if((dr = (float *)malloc((unsigned)length * sizeof(float)))==NULL){
       fprintf(stderr,"%s: [1]memory allocation failed.\n",program);
       return(0);
    }
      
    if((di = (float *)malloc((unsigned)length * sizeof(float)))==NULL){
       fprintf(stderr,"%s: [2]memory allocation failed.\n",program);
       return(0);
    }

/* 
 * convert a real data array into a complex data array make sure the 
 * imaginary part of each number is zero take the forward FFT, find 
 * the squared magnitude and finally determine the inverse FFT 
 * (autocorrelation sequence) 
 */

    if(!data_type){
       for(i=0;i<dpoints;i++){
           dr[i] = datain[i];
           di[i] = 0.0;
       }

       (void)_lfacorr(dr,di,length,dpoints,apoints);

       (void)bcopy(dr,datain,apoints*sizeof(float));

    } else if ( data_type == 1 ) {
       if ( !arith_type ) {

         /*
          * do the real part first
          */
          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              dr[j] =  datain[i];
              di[j] =  0.0;
          }

          (void)_lfacorr(dr,di,length,dpoints,apoints);

         /*
          * copy dr back into datain real part and prepare dr,di for
          * imaginary part of calculation
          */

          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              datain[i] = dr[j];
              dr[j] =  datain[i+1];
              di[j] =  0.0;
          }

          (void)_lfacorr(dr,di,length,dpoints,apoints);
 
          for(i=0,j=0;i<dpoints2;i+=2,j++) datain[i+1] = dr[j];

       } else if ( arith_type == 1) {
          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              dr[j] =  datain[i];
              di[j] =  datain[i+1];
          }

          (void)_lfacorr(dr,di,length,dpoints,apoints);

         /*
          * copy dr back into datain real part and zero imaginary 
          * part
          */

          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              datain[i] = dr[j];
              datain[i+1] = 0.0;
          }

       }
    }

    LDF_CLEANUP; 
    return(1);

} /* end of ldfacorr */


/***********************************************************************
*
*  Routine Name: static void _lfacorr()
*
*          Date: Mon Feb 17 08:51:21 MST 1992
*        
*       Purpose: performs time domain autocorrelation function in the
*                frequency domain. 
*
*         Input: float *dr   - real part of sequence
*                float *di   - imaginary part of sequence
*                int points  - number of points in input sequence (must
*                              be a power of two) 
*                int dpoints - number of points in original data sequence
*                              (need not be a power of two)
*                int apoints - number of autocorrelation points
*                              desired.
*
*        Output: returns dr as the autocorrelation sequence. 
*
*    Written By:  Jeremy Worley
*
* Modifications:
*
***********************************************************************/

static void _lfacorr(dr,di,points,dpoints,apoints)
  float *dr, *di;
  int points,dpoints,apoints;
{
  int direction,i;
  float tempr,tempi;
  int cmul();
  void fft842_();

    /* 
     * take FFT of data set 
     */

    direction = FFT_FORWARD;
    (void) fft842_ ( &direction, &points,dr,di);

    for (i=0;i<points;i++) {
        cmul (&tempr,&tempi,dr[i],di[i],dr[i],-di[i]);
        dr[i] = tempr;
        di[i] = tempi;
    }

    direction = FFT_INVERSE;
    (void) fft842_ ( &direction, &points,dr,di);

   /* 
    * scale by number of data points 
    */

    for (i = 0; i < apoints; i++ ) {
        dr[i] /= (float) dpoints;
    }

    return;
} /* end of _lacorr() */
/* -library_code_end */
