 /*
  * Khoros: $Id: ldfilconv.c,v 1.3 1992/03/25 17:30:34 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldfilconv.c,v 1.3 1992/03/25 17:30:34 dkhoros Exp $";
#endif

 /*
  * $Log: ldfilconv.c,v $
 * Revision 1.3  1992/03/25  17:30:34  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: ldfilconv.c
 >>>>
 >>>>      Program Name: dfilconv
 >>>>
 >>>> Date Last Updated: Tue Feb 25 00:01:31 1992 
 >>>>
 >>>>          Routines: ldfilconv - the library call for dfilconv
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */

static int partran(), tranpar();

/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldfilconv - library call for dfilconv
*
* Purpose:
*    
*    1D  Filter Coefficient Conversion
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    arith_type     type of arithmetic to be used on complex data.   0
*                   specifies  scalar  arithmetic and 1 specifies com-
*                   plex arithmetic.
*    
*    out_option     specifies the output  file  to  be  generated.   0
*                   specifies   lattice   to  transversal  coefficient
*                   conversion, and 1 specifies transversal to lattice
*                   conversion.
*    
*    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
*    
*    Mon Jan 27 09:50:26 MST 1992 Worley     changed  declarations  of
*                                            partran  and  tranpar  in
*                                            calling    routines    to
*                                            static.
*    
*    

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


/* -library_def */
int ldfilconv ( image, out_option, arith_type, procdir )
int out_option, arith_type, procdir;
struct xvimage *image;

/* -library_def_end */

/* -library_code */
{
  float **datain,*dataout;
  char **dload_vector();
  int   stages, size, data_type, i;
  int dunload_vector();
  int  dimension, num_vects;

  char *program = "ldfilconv";

  /* perform error checking of data types */
  if ( !(propertype(program, image, VFF_TYP_FLOAT, FALSE)) &&
       !(propertype(program, image, VFF_TYP_COMPLEX, FALSE)) )
  {
       (void) fprintf (stderr, "\n\n%s:   ", program);
       (void) fprintf (stderr, "ldfilconv: data storage type must be of type float or complex\n\n");
       return (0);
  }

  /* reorganize the data into our internal representation */
  if ((datain = (float **)dload_vector(image, &num_vects, &dimension, procdir))
            == NULL)
  {
     (void) fprintf(stderr,"%s: dload_vector failed \n",program);
     return(0);
  }

  /* determine length of array to process */
  stages = dimension;
  size = dimension;

  if ( image->data_storage_type == VFF_TYP_COMPLEX )
  {
     size *= 2;
  }

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

  /* allocate working arrays */
  dataout = (float *) malloc((unsigned int) size * sizeof(float));

  /* set working arrays to zero */
  bzero( (char *)dataout, (int) (size * sizeof(float)) );

  free ((char *) image->imagedata);

  /* begin processing for all vectors */
  for ( i = 0; i < num_vects; i++ )
  {
      if ( !out_option )
      {
         (void) partran ( data_type, arith_type, stages, datain[i], dataout );
      }
      else if ( out_option == 1 )
      {
         (void) tranpar ( data_type, arith_type, stages, dataout, datain[i] );
      }
      else
      {
         (void) fprintf (stderr,"ldfilconv: unknown filter conversion type, %d\n", out_option);
      }
      (void) bcopy ((char *)dataout, (char *)datain[i], (int)size*sizeof(float));
  } /* end of vector processing */

  /* return data to its original format */
  if (!dunload_vector((char **)datain, image, image->data_storage_type, num_vects,
 dimension, procdir))
  {
     (void) fprintf (stderr,"%s: dunload_vector failed \n",program);
     return(0);
  }

  (void) free (dataout);
  return(1);

} /* end of ldfilconv */



/************************************************************************
*
*       MODULE NAME: partran.c
*
*           PURPOSE: converts lattice filter coefficients (PARCOR)
*                    to transversal (tapped-delay or AR) filter
*                    coefficients.  Can operate on real or complex
*                    data. Also, it can perform either scalar
*                    or complex arithmetic on the complex data. 
*
**************************************************************************/

static int
partran ( datatype, arith_type, stages, parcor, tapped )

int stages, datatype, arith_type;
float *parcor, *tapped;
{
    int   i, j, row, col, stage1, cmul(), csub();
    float real, imag;
    float *fptr, *fpti, *prevr, *previ, *parr, *pari;

   /* initialize parameters */
    stage1 = stages + 1;
    real = imag = 0.0;

   /* allocate memory for the needed pointers */
    fptr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    prevr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    if ( datatype == 1)
    {
       fpti = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       previ = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       parr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       pari = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    }

   /* if data is real (datatype=0) */
    if ( !datatype )
    {
       prevr[1] = 1.0;
       prevr[0] = -1.0 * parcor[0];
      /* perform the conversion */
       for ( row = 1; row <= stages; row++ )
       {
           fptr[0] = -1.0 * parcor[row-1];
           fptr[row] = 1.0;
           for ( col = 1; col < row; col++ )
           {
               fptr[col] = prevr[col-1] - parcor[row-1] * prevr[row-1-col];
           }
           bcopy((char *) fptr, (char *) prevr, (int) stages * sizeof(float));
       }
       for ( i = 0; i < stages; i++ )
       {
           tapped[i] = fptr[stages-i-1];
       }
    }
    else if ( datatype == 1 )
    {
       for ( i = 0, j = 0; i < stages*2; i += 2, j ++ )
       {
           parr[j] = parcor[i];
           pari[j] = parcor[i+1];
       }
       prevr[1] = 1.0;
       prevr[0] = -1.0 * parr[0];
       if(!arith_type)
          previ[1] = 1.0;
       else
          previ[1] = 0.0 ;
       previ[0] = -1.0 * pari[0];
      /* perform the conversion */
       for ( row = 1; row <= stages; row++ )
       {
           fptr[0] = -1.0 * parr[row-1];
           fptr[row] = 1.0;
           fpti[0] = -1.0 * pari[row-1];
           if(!arith_type)
              fpti[row] = 1.0;
           else
              fpti[row] = 0.0 ;
           for ( col = 1; col < row; col++ )
           {
              /* scalar type of arithmetic */
               if ( !arith_type )
               {
                  fptr[col] = prevr[col-1] - parr[row-1] * prevr[row-1-col];
                  fpti[col] = previ[col-1] - pari[row-1] * previ[row-1-col];
               }
              /* complex type of arithmetic */
               else if ( arith_type == 1 )
               {
                  cmul ( &real, &imag, parr[row-1], pari[row-1],
                         prevr[row-1-col], previ[row-1-col] );
                  csub ( &fptr[col], &fpti[col],  prevr[col-1], previ[col-1],
                        real,imag );
               }
           }
           bcopy((char *) fptr, (char *) prevr, (int) stages * sizeof(float));
           bcopy((char *) fpti, (char *) previ, (int) stages * sizeof(float));
       }
       for ( i = 0, j = 0; i < stages*2; i += 2, j++ )
       {
           tapped[i] = fptr[stages-j-1];
           tapped[i+1] = fpti[stages-j-1];
       }
    }

   /* release memory */
    free((char *)fptr);
    free((char *)prevr);
    if ( datatype == 1 )
    {
       free((char *)fpti);
       free((char *)previ);
       free((char *)parr);
       free((char *)pari);
    }

} /* end of partran */


/*************************************************************************
*
*       MODULE NAME: tranpar.c
*
*           PURPOSE: converts transversal filter coefficients (AR) 
*                    to lattice (PARCOR) filter coefficients.  
*                    Can operate on real or complex data.
*                    Also, it can perform either scalar or
*                    complex arithmetic on the complex data.
*
**************************************************************************/

static int
tranpar ( datatype, arith_type, stages, parcor, tapped )

int stages, datatype, arith_type;
float *parcor, *tapped;
{
    int     i, j, row, col, stage1, cmul(), cdiv(), cadd();
    char    *malloc();
    float   doofusr,doofusi, numer, numei, denomr, denomi;
    float   *fptr, *fpti, *prevr, *previ, *parr, *pari;
    double  cmagsq();
    float tmp2r,tmp2i;

   /* initialize parameters */
    stage1 = stages + 1;
    doofusr = 0.0;
    doofusi = 0.0;

   /* allocate memory for the needed pointers */
    fptr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    prevr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    if ( datatype == 1)
    {
       fpti = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       previ = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       parr = (float *) malloc((unsigned int) (stage1)*sizeof(float));
       pari = (float *) malloc((unsigned int) (stage1)*sizeof(float));
    }

   /* if data is real (datatype=0) */
    if ( !datatype )
    {
       for ( i = 0; i < stages; i++ )
       {
           prevr[i] = tapped[stages-i-1];
       }
      /* perform the conversion */
       for ( row = stages-1; row >= 0; row-- ) {
           numer = -1.0 * prevr[0];
           denomr = 1.0 - numer * numer;
           parcor[row] = numer;
           doofusr = numer / denomr;
           for ( col = 0; col < row; col++ )
           {
               fptr[col] = prevr[col+1] / denomr + doofusr * prevr[row-col];
           }
           bcopy((char *) fptr, (char *) prevr, (int) stages * sizeof(float));
       }
    } 
    else if ( datatype == 1 )
    {
       for ( i = 0, j = 0; i < stages*2; i +=2, j++ )
       {
           prevr[stages-j-1] = tapped[i];
           previ[stages-j-1] = tapped[i+1];
       }
      /* perform the conversion */
       for ( row = stages-1; row >= 0; row-- )
       {
          /* scalar type of arithmetic */
           if ( !arith_type )
           {
              numer = -1.0 * prevr[0];
              numei = -1.0 * previ[0];
              parr[row] = numer;
              pari[row] = numei;
              denomr = 1.0 - numer * numer;
              denomi = 1.0 - numei * numei;
              doofusr = numer / denomr;
              doofusi = numei / denomi;
              for ( col = 0; col < row; col++ )
              {
                  fptr[col] = prevr[col+1]/denomr + doofusr * prevr[row-col];
                  fpti[col] = previ[col+1]/denomi + doofusi * previ[row-col];
              }
           }
           else if ( arith_type == 1 )
           {
              numer = -1.0 * prevr[0];
              numei = -1.0 * previ[0];
              doofusr = (float) cmagsq ( numer, numei );
              doofusi = 0.0;
              denomr = 1.0 - doofusr;
              denomi = 0.0 - doofusi;
              parr[row] = numer;
              pari[row] = numei;
              cdiv ( &doofusr, &doofusi, numer, numei, denomr, denomi );
              for ( col = 0; col < row; col++ )
              {
                  tmp2r = doofusr; 
                  tmp2i = doofusi;
                  cmul ( &tmp2r, &tmp2i, prevr[row-col], previ[row-col],
                          tmp2r, tmp2i );

                  cdiv ( &fptr[col],&fpti[col],prevr[col+1],previ[col+1],
                        denomr,denomi);

                  cadd ( &fptr[col], &fpti[col], tmp2r, tmp2i,
                         fptr[col],fpti[col] );
              }
           }
           bcopy((char *) fptr, (char *) prevr, (int) stages * sizeof(float));
           bcopy((char *) fpti, (char *) previ, (int) stages * sizeof(float));
       }
       for ( i = 0, j = 0; i < stages*2; i +=2, j++ )
       {
           parcor[i] = parr[j];
           parcor[i+1] = pari[j];
       }
    }

   /* release memory */
    free((char *)fptr);
    free((char *)prevr);
    if ( datatype == 1 )
    {
       free((char *)fpti);
       free((char *)previ);
       free((char *)parr);
       free((char *)pari);
    }

} /* end of tranpar */


/* -library_code_end */
