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

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

 /*
  * $Log: ldmtapp.c,v $
 * Revision 1.2  1992/03/20  23:28:33  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: ldmtapp.c
 >>>>
 >>>>      Program Name: dmtapp
 >>>>
 >>>> Date Last Updated: Thu Mar  5 10:57:36 1992 
 >>>>
 >>>>          Routines: ldmtapp - the library call for dmtapp
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: ldmtapp - library call for dmtapp
*
* Purpose:
*    
*    1D MLMS Transversal Predictor Filter
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    alpha          convergence parameter (0<alpha<1).
*    
*    beta           smoothing parameter (0<beta<1), common practice is
*                   to  have beta = 1-alpha, the optimum value is beta
*                   = 2.0/(filter order).
*    
*    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
*    
*    Jeremy Worley 05 Mar 1992 10:55 MST
*              Added  an  explicit  declaration  of  propertype()   to
*              ldmtapp().   Also  added an explicit declaration of the
*              argument delay to lmtapxx().
*    
*    

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


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

/* -library_def_end */

/* -library_code */
{
  int i, num_vects, dimension, size, stage2;
  int lmtap(), lmtapx(), lmtapxx();
  char **dload_vector();
  int  dunload_vector();
  float *error, *sigest, *tapped, **data;
  int propertype();

  char *program = "ldmtapp";

  /* 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, "ldmtapp: 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 )
      {
         lmtap ( data[i], dimension, stages, alpha, beta,
                  delay, error, sigest, tapped );
      }
      else if ( image->data_storage_type == VFF_TYP_COMPLEX )
      {
           if ( arith_type == 0 )
           {
              lmtapx ( data[i], dimension, stages, alpha, beta,
                        delay, error, sigest, tapped );
           }
           else if ( arith_type == 1 )
           {
                lmtapxx ( data[i], dimension, stages, alpha, beta,
                           delay, error, sigest, tapped );
           }
           else
           {
                (void) fprintf (stderr,"ldmtapp: 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,"ldmtapp: 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 ldmtapp */




/*********************************************************************
*
*       PROGRAM: lmtap.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 MLMS adaptive procedure.
*
*        INPUT : - datain: inputa data record
*                - points: number of points to process
*                - stages: number of filter coefficients
*                - alpha:  convergence parameter (0<alpha<1)
*                - beta:   smoothing parameter (0<beta<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

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

int   stages, points, delay;
float alpha, beta;
float *datain, *error, *sigest, *tapped;
{
    int   i, tidx, idx, limit;
    float sigma;

   /* initialize parameters */
    sigma = 1.0;

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

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

           /* evaluate the instantaneous power of the sample */
            sigma = beta*sigma + (1.0-beta)*datain[tidx]*datain[tidx];
            if (sigma <= 0.0)
            {
                sigma = 0.0000001;
            }

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

    } /* end of main predictor loop */

} /* end of lmtap */




/*********************************************************************
*
*       PROGRAM: lmtapx.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 MLMS 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)
*                - beta:   smoothing parameter (0<beta<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

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

int   delay, stages, points;
float alpha, beta;
float *datain, *error, *sigest, *tapped;
{
    int   i, limit, twostages, idx, tidx, twodelay, size;
    float sigmar, sigmai;

   /* initialize parameters */
    size = points * 2;
    sigmar = sigmai = 1.0;
    twodelay = 2 * delay;
    twostages = 2 * stages;

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

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

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

           /* evaluate the instantaneous power of the sample */
            sigmar = beta*sigmar + (1.0-beta)*datain[tidx]*datain[tidx];
            sigmai = beta*sigmai + (1.0-beta)*datain[tidx+1]*datain[tidx+1];
            if (sigmar <= 0.0)
            {
                sigmar = 0.0000001;
            }
            if (sigmai <= 0.0)
            {
                sigmai = 0.0000001;
            }

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

    } /* end of main predictor loop */

} /* end of lmtapx */





/*********************************************************************
*
*       PROGRAM: lmtapxx.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 MLMS 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)
*                - beta:   smoothing parameter (0<beta<1)
*                - delay:  time-delay parameter
*
*       OUTPUT : - error:  error sequence
*                - sigest: signal estimate
*                - tapped: filter coefficients 
*
*********************************************************************/

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

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

   /* initialize parameters */
    newbeta = 1.0 - beta;
    size = points * 2;
    real = 1.0;
    imag = tempr = tempi = 0.0;
    twodelay = 2.0 * delay;
    twostages = 2 * stages;

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

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

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

           /* evaluate the instantaneous power of the sample */
            cmul (&real, &imag, real, imag, beta, 0.0);
            cmul (&tempr, &tempi, datain[tidx], datain[tidx+1],
                  datain[tidx], -datain[tidx+1]);
            cmul (&tempr, &tempi, tempr, tempi, newbeta, 0.0);
            cadd (&real, &imag, real, imag, tempr, tempi);
            if (real <= 0.0)
            {
                real = 0.0000001;
            }
            if (imag <= 0.0)
            {
                imag = 0.0000001;
            }

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

    } /* end of main predictor loop */

} /* end of lmtapxx */
/* -library_code_end */
