 /*
  * Khoros: $Id: ldltapp.c,v 1.2 1992/03/20 23:28:02 dkhoros Exp $
  */

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

 /*
  * $Log: ldltapp.c,v $
 * Revision 1.2  1992/03/20  23:28:02  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: ldltapp.c
 >>>>
 >>>>      Program Name: dltapp
 >>>>
 >>>> Date Last Updated: Sat Feb 15 16:27:11 1992 
 >>>>
 >>>>          Routines: ldltapp - the library call for dltapp
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: ldltapp - library call for dltapp
*
* Purpose:
*    
*    1D  LMS Transversal Predictor Filter
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    alpha          convergence parameter (0<alpha<1).
*    
*    delay          time delay parameter (delay >= 0).
*    
*    stages         used to determine the dimension (length or  number
*                   of filter coefficients) of the filter.  If it is a
*                   positive, non-zero integer.  Default value is  set
*                   to 1.
*    
*    arith_type     type of arithmetic to be used on complex data.   0
*                   specifies scalar arithmetic and 1 specifies vector
*                   arithmetic.
*    
*    out_option     specifies the output  file  to  be  generated.   0
*                   specifies  an error sequence, 1 specifies a signal
*                   estimate sequence, and 2 specifies a filter  coef-
*                   ficient sequence.
*    
*    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
*    
*    

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


/* -library_def */
int ldltapp ( image, alpha, delay, stages, arith_type, out_option, procdir )
int delay, stages, out_option, arith_type, procdir;
float alpha;
struct xvimage *image;

/* -library_def_end */

/* -library_code */
{
  int i, num_vects, dimension, size, stage2;
  int lltap(), lltapx(), lltapxx();
  char **dload_vector();
  int dunload_vector();
  float *error, *sigest, *tapped, **data;

  char *program = "ldltapp";

  /* 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, "ldltapp: data storage type must be of type float or complex\n\n");
       return (0);
  }

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

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

  /* allocate working arrays */
  error = (float *) malloc((unsigned int) size * sizeof(float));
  sigest = (float *) malloc((unsigned int) size * sizeof(float));
  tapped = (float *) malloc((unsigned int) stage2 * sizeof(float));

  /* set working arrays to zero */
  bzero( (char *)error, (int) (size * sizeof(float)) );
  bzero( (char *)sigest, (int) (size * sizeof(float)) );
  bzero( (char *)tapped, (int) (stage2 * sizeof(float)) );

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

  /* begin processing for all vectors */ 
  for ( i = 0; i < num_vects; i++ )
  {
      if ( image->data_storage_type == VFF_TYP_FLOAT )
      {
         lltap ( data[i], dimension, stages, alpha, delay,
                    error, sigest, tapped );
      }
      else if ( image->data_storage_type == VFF_TYP_COMPLEX )
      {
           if ( arith_type == 0 )
           {
              lltapx ( data[i], dimension, stages, alpha, delay,
                          error, sigest, tapped );
           }
           else if ( arith_type == 1 )
           {
                lltapxx ( data[i], dimension, stages, alpha, delay,
                             error, sigest, tapped );
           }
           else
           {
                (void) fprintf (stderr,"ldltapp: unknown arithmetic type, %d\n", arith_type);
           }
      }
      if ( out_option == 0 )
      {
         (void) bcopy ((char*)error, (char *) data[i], (int)size*sizeof(float)); 
      }
      else if ( out_option == 1 )
      {
           (void) bcopy ((char*)sigest, (char *) data[i], (int)size*sizeof(float
));
      }
      else if ( out_option == 2 )
      {
           (void) bcopy ((char*)tapped, (char *) data[i], (int)stage2*sizeof(float));
      }
      else
      {
           (void) fprintf (stderr,"ldltapp: unknown output file selection, %d\n", out_option);
      }

  } /* end of vector processing */

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

  return(1); 

} /* end of ldltapp */




/*********************************************************************
*
*       PROGRAM: lltap.c
*
*       PURPOSE: This routine generates a prediction error sequence,
*                a signal estimate sequence and a filter weights 
*                sequence from a transversal filter realization.
*                Uses the LMS adaptive procedure.
*
*        INPUT : - datain: inputa data record
*                - points: number of points to process
*                - stages: number of filter coefficients
*                - alpha:  convergence parameter (0<alpha<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

int
lltap ( datain, points, stages, alpha, delay, error, sigest, tapped )

int   stages, points, delay;
float alpha, *datain, *error, *sigest, *tapped;
{
    int   i, tindex, idx, limit;
    float twoalpha;

   /* initialize parameters */
    twoalpha = 2.0 * alpha;

   /* begin adaptive transversal predictor filter */
    for ( tindex = 0; tindex < points; tindex++) 
    {
        if (tindex < delay)
        {
            error[tindex] = datain[tindex];
        }
        else
        {
           /* determine how many filter weights have a signal component */
            if (tindex <= stages)
            {
                limit = tindex;
            }
            else
            {
                limit = stages;
            }
           /* evaluate an estimate for the signal */
            for (i = 0; i < limit; i++)
            {
                idx = tindex - i - delay;
                if (idx >= 0)
                {
                   sigest[tindex] = sigest[tindex] + tapped[i] * datain[idx];
                }
            }

           /* determine the error sequence */
            error[tindex] = datain[tindex] - sigest[tindex];

           /* update the filter transversal coefficients */
            for (i = 0; i < limit; i++)
            {
                idx = tindex - i - delay;
                if (idx >= 0)
                {
                   tapped[i] = tapped[i] + 
                               twoalpha * error[tindex] * datain[idx];
                }
            }
        }

    } /* end of main predictor loop */

} /* end of lltap */




/*********************************************************************
*
*       PROGRAM: lltapx.c
*
*       PURPOSE: This routine generates a prediction error sequence,
*                a signal estimate sequence and a filter weights 
*                sequence from a transversal filter realization.
*                Uses the LMS adaptive procedure.
*                Operates on complex data using real arithmetic,
*                that is, real and imaginary numbers are treated 
*                separatedly when processed.
*
*        INPUT : - datain: inputa data record
*                - points: number of points to process
*                - stages: number of filter coefficients
*                - alpha:  convergence parameter (0<alpha<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

int
lltapx ( datain, points, stages, alpha, delay, error, sigest, tapped )

int   delay, stages, points;
float alpha, *datain, *error, *sigest, *tapped;
{
    int   i, limit, twostages, idx, tindex, twodelay, size;
    float twoalpha;

   /* initialize parameters */
    size = points * 2;
    twodelay = 2 * delay;
    twoalpha = 2.0 * alpha;
    twostages = 2 * stages;

   /* begin adaptive transversal predictor filter */
    for ( tindex = 0; tindex < size; tindex += 2 ) 
    {
        if (tindex < twodelay)
        {
            error[tindex] = datain[tindex];
            error[tindex+1] = datain[tindex+1];
        }
        else
        {
           /* determine how many filter weights have a signal component */
            if (tindex <= twostages)
            {
                limit = tindex;
            }
            else
            {
                limit = twostages;
            }

           /* evaluate an estimate for the signal */
            for (i = 0; i < limit; i += 2)
            {
                idx = tindex - i - twodelay;
                if (idx >= 0)
                {
                   sigest[tindex] = sigest[tindex] + tapped[i] * datain[idx];
                   sigest[tindex+1] = sigest[tindex+1] + 
                                      tapped[i+1] * datain[idx+1];
                }
            }

           /* determine the error sequence */
            error[tindex] = datain[tindex] - sigest[tindex];
            error[tindex+1] = datain[tindex+1] - sigest[tindex+1];

           /* update the filter transversal coefficients */
            for (i = 0; i < limit; i += 2)
            {
                idx = tindex - i - twodelay;
                if (idx >= 0)
                {
                   tapped[i] = tapped[i] +
                               twoalpha * error[tindex] * datain[idx];
                   tapped[i+1] = tapped[i+1] +
                               twoalpha * error[tindex+1] * datain[idx+1];
                }
            }
        }

    } /* end of main predictor loop */

} /* end of lltapx */





/*********************************************************************
*
*       PROGRAM: lltapxx.c
*
*       PURPOSE: This routine generates a prediction error sequence,
*                a signal estimate sequence and a filter weights 
*                sequence from a transversal filter realization.
*                Uses the LMS adaptive procedure.
*                Operates on complex data using complex arithmetic,
*                that is, real and imaginary numbers are treated 
*                as one number when processed.
*
*        INPUT : - datain: inputa data record
*                - points: number of points to process
*                - stages: number of filter coefficients
*                - alpha:  convergence parameter (0<alpha<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

int
lltapxx ( datain, points, stages, alpha, delay, error, sigest, tapped )

int   stages, points, delay;
float alpha, *datain, *error, *sigest, *tapped;
{
    int     i, tindex, idx, size, limit, twostages, twodelay;
    int     cadd(), csub(), cmul(), cdiv();
    float   tempr, tempi, twoalpha;

   /* initialize parameters */
    size = points * 2;
    twoalpha = 2.0 * alpha;
    twodelay = 2.0 * delay;
    twostages = 2 * stages;

   /* begin adaptive transversal predictor filter */
    for ( tindex = 0; tindex < size; tindex += 2 ) 
    {
        if (tindex < twodelay)
        {
            error[tindex] = datain[tindex];
            error[tindex+1] = datain[tindex+1];
        }
        else
        {
           /* determine how many filter weights have a signal component */
            if (tindex <= twostages)
            {
                limit = tindex;
            }
            else
            {
                limit = twostages;
            }

           /* evaluate an estimate for the signal */
            for (i = 0; i < limit; i += 2)
            {
                idx = tindex - i - twodelay;
                if (idx >= 0)
                {
                   cmul (&tempr, &tempi, tapped[i], tapped[i+1],
                         datain[idx], datain[idx+1]);
                   cadd (&sigest[tindex], &sigest[tindex+1], 
                         sigest[tindex], sigest[tindex+1], tempr, tempi);
                }
            }

           /* determine the error sequence */
            csub (&error[tindex], &error[tindex+1], 
                  datain[tindex], datain[tindex+1],
                  sigest[tindex], sigest[tindex+1]);

           /* update the filter transversal coefficients */
            for (i = 0; i < limit; i += 2)
            {
                idx = tindex - i - twodelay;
                if (idx >= 0)
                {
                   cmul (&tempr, &tempi, error[tindex], error[tindex+1],
                         datain[idx], -datain[idx+1]);
                   cmul (&tempr, &tempi, twoalpha, 0.0,
                         tempr, tempi);
                   cadd (&tapped[i], &tapped[i+1], tapped[i], tapped[i+1],
                         tempr, tempi);
                }
            }
        }

    } /* end of main predictor loop */

} /* end of lltapxx */


/* -library_code_end */
