 /*
  * Khoros: $Id: ldconvo.c,v 1.2 1992/03/20 23:25:13 dkhoros Exp $
  */

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

 /*
  * $Log: ldconvo.c,v $
 * Revision 1.2  1992/03/20  23:25:13  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: ldconvo.c
 >>>>
 >>>>      Program Name: dconvo
 >>>>
 >>>> Date Last Updated: Tue Feb 25 00:43:38 1992 
 >>>>
 >>>>          Routines: ldconvo - the library call for dconvo
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* the cleanup routines free up all allocated memory prior to a return */
#define LDA_CLEANUP {int ii; \
                     if(data1!=NULL){ \
                         for(ii=0;i<num_vects1;ii++) \
                             if(data1[ii]!=NULL)free(data1[ii]); \
                         free(data1); \
                      } \
                     if(data2!=NULL){ \
                         for(ii=0;i<num_vects2;ii++) \
                             if(data2[ii]!=NULL)free(data2[ii]); \
                         free(data2); \
                      }} 
#define LDF_CLEANUP { \
                      if(d1r!=NULL)free(d1r); \
                      if(d1i!=NULL)free(d1i); \
                      if(d2r!=NULL)free(d2r); \
                      if(d2i!=NULL)free(d2i); \
                      }
#define LDT_CLEANUP {if(convi!=NULL)free(convi);if(convr!=NULL)free(convr);}

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


/****************************************************************
*
* Routine Name: ldconvo - library call for dconvo
*
* Purpose:
*    
*    Convolution sequence of 1-D signal(s)
*    
*    

* Input:
*    
*    image1         pointer to VIFF structure containing  first  image
*                   data to be processed.
*    
*    image2         pointer to VIFF structure containing second  image
*                   data to be processed.
*    
*    arith_type     an integer specifying the type of arithmetic to be
*                   performed  if  the  data is complex.  A value of 0
*                   indicates scalar arithmetic, a value  of  1  indi-
*                   cates vector arithmetic.
*    
*    proc_option    an integer specifying what type of convolution  is
*                   to  be  performed.   If  a value of 0 is passed in
*                   (the default),  then  a  time  domain  convolution
*                   estimate  is  performed.  If a 1 passed in, then a
*                   frequency domain convolution is performed.
*    
*    procdir        process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

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

*
* Written By: Jeremy Worley, Ramiro Jordan
*    
*    Jeremy Worley 24 Feb 1992 22:15 MST
*              changed call to lfft1d in ldfconvo() to fft842_() since
*              lfft1d  is  being  phased  out.   Also  cleaned up some
*              implicitly defined functions and modularized ldfconvo()
*              which  resulted  in  the  creation of a static function
*              called _lfconvo().  Also fixed a small scaling  problem
*              with the output of ldfconvo().
*    
*    

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


/* -library_def */
int
ldconvo ( image1, image2, arith_type, proc_option, process_dir )

int proc_option, arith_type, process_dir;
struct xvimage *image1, *image2;
/* -library_def_end */

/* -library_code */
{
    int   i, data_type, dimension1, dimension2, num_vects1, num_vects2; 
    int   dunload_vector(), lscale();
    float **data1, **data2;
    char  **dload_vector();
    int ldfconvo(),ldtconvo();
    char *program = "ldconvo";

/* 
** initialize parameters 
*/

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

/*
** load data into internal working format.
*/

    if((data1 = (float **)dload_vector(image1,&num_vects1,&dimension1, 
                process_dir)) == NULL){
       fprintf(stderr,"%s: [1]dload_vector failed \n",program);
       return(0);
    }

    if((data2 = (float **)dload_vector(image2,&num_vects2,&dimension2, 
                process_dir)) == NULL){
       fprintf(stderr,"%s: [2]dload_vector failed \n",program);
       return(0);
    }

/*
** check characteristics of data and parameters to prevent errors later
*/

    if(image1->data_storage_type != image2->data_storage_type){
       fprintf(stderr,"%s: data sets must be of same data type.\n",program);
       return(0);
    }

    if(num_vects1 != num_vects2 || dimension1 != dimension2){
       fprintf(stderr,"%s:  dimensions of data sets must match.\n",program);
       return(0);
    }

/* 
** find convolution sequence 
*/

    if ( proc_option == 1 ){
       for (i=0;i<num_vects1;i++){
           if (!ldfconvo(data1[i],data2[i],dimension1,data_type,arith_type)){
              fprintf ( stderr,"ldfconvo failed\n");
              LDA_CLEANUP;
              return(0);
           }
       }
    }else{
       for(i=0;i<num_vects1;i++){
           if(!ldtconvo(data1[i],data2[i],dimension1,data_type,arith_type)){
              fprintf ( stderr,"ldtconvo failed\n");
              LDA_CLEANUP;
              return(0);
           }
       }
    }

/*
** unload data into its original format
*/

    if (!dunload_vector((char **)data1, image1, image1->data_storage_type, 
                       num_vects1, dimension1, process_dir)){
        fprintf (stderr,"%s: dunload_vector failed \n",program);
        LDA_CLEANUP;
        return(0);
    }

/*
** clean up our little mess and exit normally
*/

    LDA_CLEANUP;
    return(1);
} /* end of ldconvo */


/***********************************************************************
*
*  Routine Name: ldtconvo()
*
*          Date: 01-AUG-90
*        
*       Purpose: Performs the circular convolution of two real or 
*                complex data vectors. Complex data is assumed to
*                have real and imaginary parts interlaced.
*
*         Input: datain1    - first data vector.
*                datain2    - second data vector. 
*                dpoints    - number of data points. 
*                data_type  - 0 - real data
*                             1 - complex data
*                arith_type - 0 - scalar arithmetic on complex data
*                             1 - vector arithmetic on complex data
*
*        Output: datain1    - resultant convolution data 
*                return code: 0 - error condition
*                             1 - normal return
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int ldtconvo (datain1,datain2,dpoints,data_type,arith_type)
    int   arith_type, dpoints, data_type;
    float *datain1, *datain2;
{
    int   idx,n,h;    
    int   cadd(), cmul();
    float *convr, *convi, real, imag, tempr, tempi;
    char *program = "ldtconvo";

/*
** allocate space for working arrays 
*/

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

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

/*
** find the convolution estimate 
** first real data
*/

    if (!data_type){
       for (n=0;n<dpoints;n++){
           tempr = 0.0;
           for (h=0;h<dpoints;h++){
               idx = n-h; 
               if(idx<0)idx+=dpoints;
               tempr += datain1[h] * datain2[idx];
           }
           convr[n] = tempr;
       }
       bcopy ((char *)convr,(char *)datain1,(int)dpoints*sizeof(float));

/*
** complex data
*/

    }else{
       if (!arith_type){              /* scalar arithmetic */
          for (n=0;n<dpoints;n++){
              tempr = 0.0;
              tempi = 0.0;
              for (h=0;h<dpoints;h++){
                  idx = n-h; 
                  if(idx<0)idx += dpoints;
                  tempr += datain1[2*h]   * datain2[2*idx];
                  tempi += datain1[2*h+1] * datain2[2*idx+1];
              }
              convr[n] = tempr;
              convi[n] = tempi;
          }
          for (n=0;n<dpoints;n++) {
              datain1[2*n]   = convr[n];
              datain1[2*n+1] = convi[n];
          }
       } else {                        /* vector arithmetic */
          for (n=0;n<dpoints;n++){
              real = 0.0;
              imag = 0.0;
              tempr = 0.0;
              tempi = 0.0;
              for (h=0;h<dpoints;h++){
                  idx = n-h;
                  if(idx<0)idx += dpoints;
                  cmul (&tempr,&tempi,datain1[2*h],datain1[2*h+1],
                        datain2[2*idx],datain2[2*idx+1]);
                  cadd (&real,&imag,real,imag,tempr,tempi);
              } /* end for h */
              convr[n] = real;
              convi[n] = imag;
          } /* end for n */
          for (n=0;n<dpoints;n++) {
              datain1[2*n]   = convr[n];
              datain1[2*n+1] = convi[n];
          } /* end for n */
       } /* end else arith_type */
    } /* end else complex */

    LDT_CLEANUP;
    return(1);

} /* end of ldtconvo */


/***********************************************************************
*
*  Routine Name: ldfconvo()
*
*          Date: 01-AUG-90
*        
*       Purpose: Performs the frequency domain linear or circular 
*                convolution of two real or complex data vectors. 
*                Complex data is assumed to have real and imaginary 
*                parts interlaced.
*
*         Input: datain1     - first data vector.
*                datain2     - second data vector. 
*                dpoints     - number of data points. 
*                data_type   - 0 - real data
*                              1 - complex data
*                arith_type  - 0 - scalar arithmetic on complex data
*                              1 - vector arithmetic on complex data
*                proc_option - 0 - linear convolution
*                              1 - circular convoltion
*
*        Output: datain1     - resultant convolution data 
*                return code:  0 - error condition
*                              1 - normal return
*
*    Written By: Jeremy Worley 
*
* Modifications:
*
***********************************************************************/

int ldfconvo (datain1,datain2,dpoints,data_type,arith_type)
    int arith_type, dpoints, data_type;
    float *datain1, *datain2;
{
    int i, j, length, dpoints2;    
    int powtwo();
    int cmul();
    float *d1r, *d1i, *d2r, *d2i;

    char *program = "ldfconvo";

/*
 * initialize parameters 
 */

    dpoints2 = dpoints * 2;

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

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

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

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

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

    if((d2i = (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, multiply the two sequences together and 
 * finally determine the inverse FFT (convolution sequence) 
 */

    if ( !data_type ){
       for(i=0;i<dpoints;i++){
           d1r[i] = datain1[i];
           d2r[i] = datain2[i];
           d1i[i] = d2i[i] = 0.0;
       }

       (void)_lfconvo(d1r,d1i,d2r,d2i,length,dpoints,dpoints);

       (void)bcopy(d1r,datain1,dpoints*sizeof(float));

/*
 * complex data
 */

    }else{
       if(!arith_type){

          /*
           * do the real part first
           */

          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              d1r[j] =  datain1[i];
              d2r[j] =  datain2[i];
              d1i[j] =  d2i[j] = 0.0;
          }

          (void)_lfconvo(d1r,d1i,d2r,d2i,length,dpoints,dpoints);

         /*
          * 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++ ) {
              datain1[i] = d1r[j];
              d1r[j] =  datain1[i+1];
              d2r[j] =  datain2[i+1];
              d1i[j] =  d2i[j] = 0.0;
          }

          (void)_lfconvo(d1r,d1i,d2r,d2i,length,dpoints,dpoints);

          for(i=0,j=0;i<dpoints2;i+=2,j++) datain1[i+1] = d1r[j];

       } else {         /* vector arithmetic */
          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              d1r[j] =  datain1[i];
              d1i[j] =  datain1[i+1];
              d2r[j] =  datain2[i];
              d2i[j] =  datain2[i+1];
          }

          (void)_lfconvo(d1r,d1i,d2r,d2i,length,dpoints,dpoints);

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

          for (i=0,j=0; i < dpoints2;i+=2,j++ ) {
              datain1[i] = d1r[j];
              datain1[i+1] = d1i[j];
          }

       }
    }

    LDF_CLEANUP; 
    return(1);

} /* end of ldfconvo */

/***********************************************************************
*
*  Routine Name: static void _lfconvo()
*
*          Date: Mon Feb 17 08:51:21 MST 1992
*
*       Purpose: performs time domain convolution function in the
*                frequency domain.
*
*         Input: float *d1r   - real part of first sequence
*                float *d1i   - imaginary part of first sequence
*                float *d2r   - real part of second sequence
*                float *d2i   - imaginary part of second 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 cpoints - number of convolution points
*                              desired.
*
*        Output: returns dr as the autocorrelation sequence.
*
*    Written By:  Jeremy Worley
*
* Modifications:
*
***********************************************************************/

static void _lfconvo(d1r,d1i,d2r,d2i,points,dpoints,cpoints)
  float *d1r, *d1i,*d2r,*d2i;
  int points,dpoints,cpoints;
{
  int direction,i;
  float tempr,tempi;
  int cdiv(),cmul();
  void fft842_();

    /*
     * take FFT of data set
     */

    direction = FFT_FORWARD;
    (void) fft842_ ( &direction, &points,d1r,d1i);
    (void) fft842_ ( &direction, &points,d2r,d2i);

    for (i=0;i<points;i++) {
        cmul (&tempr,&tempi,d1r[i],d1i[i],d2r[i],d2i[i]);
        d1r[i] = tempr;
        d1i[i] = tempi;
    }

    direction = FFT_INVERSE;
    (void) fft842_ ( &direction, &points,d1r,d1i);

   /*
    * scale by number of data points
    */

    for (i = 0; i < cpoints; i++ ) {
        d1r[i] /= (float) dpoints;
        d1i[i] /= (float) dpoints;
    }

    return;
} /* end of _lfconvo() */

/* -library_code_end */
