 /*
  * Khoros: $Id: ldlsynth.c,v 1.2 1992/03/20 23:27:54 dkhoros Exp $
  */

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

 /*
  * $Log: ldlsynth.c,v $
 * Revision 1.2  1992/03/20  23:27:54  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, 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 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldlsynth.c
 >>>>
 >>>>      Program Name: dlsynth
 >>>>
 >>>> Date Last Updated: Mon Feb 17 11:20:40 1992 
 >>>>
 >>>>          Routines: ldlsynth - the library call for dlsynth
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldlsynth - library call for dlsynth
*
* Purpose:
*    
*    1D Lattice Filter Synthesis
*    
*    

* Input:
*    
*    lattice        pointer  to  VIFF  structure  containing   lattice
*                   weights.
*    
*    driver         pointer to VIFF structure containing driving  sig-
*                   nal.
*    
*    arith_type     type of arithmetic, scalar or vector when  dealing
*                   with complex data.
*    
*    procdir        process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

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

*
* Written By: Ramiro Jordan
*    
*    Jeremy Worley Mar 7 1990                Rewrite of ldlsynth()  to
*                                            handle  vector  and  band
*                                            oriented  signal   ensem-
*                                            bles.
*    
*    Ramiro Jordan & Jeremy Worley Feb 11 1992Fixed  bug  in  size  of
*                                            complex  data  sequences.
*                                            (Old      version      of
*                                            dload_vector     returned
*                                            dimension parameter twice
*                                            its  true  value for com-
*                                            plex.    That   bug   was
*                                            fixed,  but  was not pro-
*                                            pagated to this program.
*    
*    Jeremy Worley Feb 17 1992               1. Added explicit  return
*                                            type to ldlsynth().
*                                            2.    Made     implicitly
*                                            defined  functions expli-
*                                            cit in ldlsynth().
*                                            3.  Made  data_type   and
*                                            arith_type  arguments  to
*                                            synthesis()    explicitly
*                                            defined as int.
*                                            4. Added a return to  the
*                                            end of synthesis().
*                                            5.            Initialized
*                                            parr,pari,backr,backi  to
*                                            NULL in synthesis().
*    
*    

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


/* -library_def */
int ldlsynth ( lattice, driver, arith_type, procdir )
int arith_type,procdir;
struct xvimage *lattice, *driver;
/* -library_def_end */

/* -library_code */
{
    int  i, nc, nr, stages, data_type, synthesis();
    float *parcor, **output;
    float **temp1;
    char **dload_vector();
    int dunload_vector(),propertype(),lvconvert();
    int  size,dimension, numvects;

    char *program = "ldlsynth";

    /* initialize parameters */
    data_type = 0;

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

    /* if data storage type is integer convert to float */
    if (lattice->data_storage_type == VFF_TYP_4_BYTE )
    {
       if ( !lvconvert ( lattice, VFF_TYP_FLOAT, FALSE, TRUE, 0.0, 1.0,1) )
       {
            (void) fprintf(stderr,"lvconvert Failed\n");
            return(0);
       }
    }

    /* if data storage type is integer convert to float */
    if (driver->data_storage_type == VFF_TYP_4_BYTE )
    {
       if ( !lvconvert ( driver, VFF_TYP_FLOAT, FALSE, TRUE, 0.0, 1.0,1) )
       {
            (void) fprintf(stderr,"lvconvert Failed\n");
            return(0);
       }
    }

/*
** determine length of array to process and reorganize 
** driver data if necessary
*/

    temp1 = (float **)dload_vector(driver,&numvects,&dimension,procdir);
    if(temp1==NULL){
       fprintf(stderr,"%s:  Unable to load vector.\n",program);
       return(0);
    }

/*
** determine the number of lattice filter stages 
*/

    nc = lattice->row_size;     /* number of columns */
    nr = lattice->col_size;     /* number of rows    */
    stages = nc * nr;
    parcor = (float *) lattice->imagedata;

/*
** allocate memory for synthesized output file
*/

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

    output = (float **) malloc((unsigned int) numvects * sizeof(float));
    if(output==NULL){
       fprintf(stderr,"%s:  failure during malloc\n",program);
       return(1);
    }

/*
**  apply lattice to each signal
*/

    size=dimension;
    /* if(data_type)size/=2; */
   
     for(i=0;i<numvects;i++){
        output[i] = (float *)malloc((unsigned int)dimension*sizeof(float));
        if(output[i]==NULL){
           fprintf(stderr,"%s:  failure during malloc\n",program);
           return(0);
        }

        (void)synthesis(data_type,arith_type,temp1[i],size,parcor,stages,
                        output[i]);
    }

/*
** cleanup
*/

    (void) free ((char *)driver->imagedata);
    if(!dunload_vector((char **)output,driver,driver->data_storage_type,numvects,dimension,procdir)){
       fprintf(stderr,"%s:  Unable to unload vector.\n",program);
       return(0);
    }
    return(1);


} /* end of ldlsynth */


/**************************************************************
*
* MODULE NAME: synthesis.c
*
*     PURPOSE:  
*
*
*       INPUT MAIN: 
*
*
*      OUTPUT MAIN: 
*
*
* ROUTINES CALLED: 
*
*
**************************************************************/

int
synthesis ( data_type, arith_type, inarray, length, 
            parcor, stages, outarray )

int   length, stages, data_type, arith_type;
float *inarray, *parcor, *outarray;

{
    int   i, j, k, size, stage1, cadd(), csub(), cmul();
    char  *malloc();
    float tempr, tempi, outr, outi;
    float *parr=NULL, *pari=NULL, *backr=NULL, *backi=NULL;
  
    /* initialize parameters */
    stage1 = stages + 1;

    /* if data type is real (data_type = 0) */
    if ( !data_type )
    {
       backr = (float *) malloc((unsigned int) stage1 * sizeof(float));
       for ( k = 0; k < length; k++ )
       {
           outr = inarray[k];
           for ( i = stages; i >= 0; i-- )
           {
               outr = outr + parcor[i] * backr[i];
               backr[i+1] = backr[i] - parcor[i] * outr;
           }
           backr[0] = outr;
           outarray[k] = outr;
        }
     }
     else if ( data_type == 1 )         /* complex data */
     {
        size = length * 2;
        backr = (float *) malloc((unsigned int) stage1 * sizeof(float));
        backi = (float *) malloc((unsigned int) stage1 * sizeof(float));
        parr = (float *) malloc((unsigned int) stages * sizeof(float));
        pari = (float *) malloc((unsigned int) stages * sizeof(float));
        for ( i = 0, j = 0; i < stages*2; i += 2, j++ )
        {
            parr[j] = parcor[i];
            pari[j] = parcor[i+1];
        }
        if ( !arith_type )              /* scalar arithmetic */
        {
           for ( k = 0; k < size; k += 2  )
           {
               outr = inarray[k];
               outi = inarray[k+1];
               for ( i = stages; i >= 0; i-- )
               {
                   outr = outr + parr[i] * backr[i];
                   outi = outi + pari[i] * backi[i];
                   backr[i+1] = backr[i] - parr[i] * outr;
                   backi[i+1] = backi[i] - pari[i] * outi;
               }
               backr[0] = outr;
               backi[0] = outi;
               outarray[k] = outr;
               outarray[k+1] = outi;
           }
        }
        else if ( arith_type == 1 )     /* vector arithmetic */
        {
           for ( k = 0; k < size; k += 2  )
           {
               outr = inarray[k];
               outi = inarray[k+1];
               for ( i = stages; i >= 0; i-- )
               {
                   (void) cmul ( &tempr, &tempi, parr[i], -pari[i], 
                                 backr[i], backi[i] );
                   (void) cadd ( &outr, &outi, tempr, tempi, outr, outi );
                   (void) cmul ( &tempr, &tempi, parr[i], pari[i],
                                 outr, outi );
                   (void) csub ( &backr[i+1], &backi[i+1], backr[i], backi[i],
                                 tempr, tempi );
               } 
               backr[0] = outr;
               backi[0] = outi;
               outarray[k] = outr;
               outarray[k+1] = outi;
           }
        }
     }
 
     (void) free((char *) backr);
     if ( data_type == 1 )
     {
        (void) free((char *) backi);
        (void) free((char *) parr);
        (void) free((char *) pari);
     }

     return(1);
} /* end of synthesis */



/* -library_code_end */
