 /*
  * Khoros: $Id: ldwindow.c,v 1.2 1992/03/20 23:29:23 dkhoros Exp $
  */

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

 /*
  * $Log: ldwindow.c,v $
 * Revision 1.2  1992/03/20  23:29:23  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: ldwindow.c
 >>>>
 >>>>      Program Name: dwindow
 >>>>
 >>>> Date Last Updated: Sat Feb 15 16:37:06 1992 
 >>>>
 >>>>          Routines: ldwindow - the library call for dwindow
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: ldwindow - library call for dwindow
*
* Purpose:
*    
*    1D Window (Apodization) Functions.
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    window_option  specifies the window function to be applied to the
*                   input data.
*    
*    theta          theta parameter to be supplied to the kaiser func-
*                   tion.
*    
*    damp           damping factor for exponential functions.
*    
*    width          desired window width for user defined window.
*    
*    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, Glen Abousleman
*    
*    28 Feb 90 Jeremy Worley  Added process direction code and  multi-
*                             signal capability
*    
*    23 Jul 90  Jeremy  Worley  Corrected  incorrect  data  types  and
*                             function/variable  conflicts  in  bess()
*                             and fact().
*    
*    04 Aug 90 Jeremy Worley  Changed  dload_vector()  declaration  to
*                             char**;  added explicit return values in
*                             each window function; corrected indexing
*                             problems  in  each  of  the window func-
*                             tions; fixed indexing  inconsistancy  in
*                             driver  code;  moved  low level routines
*                             into vmath.a;
*    
*    

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


/* -library_def */
int
ldwindow ( image, window_option, theta, damp, width, procdir )

int    window_option, width, procdir;
float  theta, damp;
struct xvimage *image;
/* -library_def_end */

/* -library_code */
{
    int dimension, numvects, data_type,i;
    float *array,**temp1;
    char **dload_vector();
    int dunload_vector();
    char *program = "ldwindow";

/*
** initialize parameters
*/

    data_type = 0;

/* 
** perform error checking
*/

    if ( !(propertype(program, image, VFF_TYP_FLOAT, FALSE)) &&
         !(propertype(program, image, VFF_TYP_COMPLEX, FALSE)) &&
         !(propertype(program, image, VFF_TYP_4_BYTE, FALSE)) )
    {
          (void) fprintf (stderr, "\n\n%s:   ", program);
          (void) fprintf (stderr, "ldwindow: 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 (image->data_storage_type == VFF_TYP_4_BYTE )
    {
       if ( !lvconvert ( image, VFF_TYP_FLOAT, FALSE, TRUE, 0.0, 1.0,0 ) )
       {
            (void) fprintf(stderr,"lvconvert Failed\n");
            return(0);
       }
    }

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

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

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

/*
**  loop through all signals and apply window to each
*/

    for(i=0;i<numvects;i++){
        array = temp1[i];

        switch ( window_option ) {
        case 0:
             break;
        case 1:
             (void) bartlett ( array, dimension, data_type );
             break;
        case 2:
             (void) hamming ( array, dimension, data_type );
             break;
        case 3:
             (void) hanning ( array, dimension, data_type );
             break;
        case 4:
             (void) parzen ( array, dimension, data_type );
             break;
        case 5:
             (void) kaiser ( array, dimension, data_type, theta );
             break;
        case 6:
             (void) exponential ( array, dimension, data_type, damp );
             break;
        case 7:
             (void) blackman ( array, dimension, data_type );
             break;
        case 8:
             (void) blarris ( array, dimension, data_type );
             break;
        case 9:
             (void) rtriangular ( array, dimension, data_type );
             break;
        case 10:
             (void) ltriangular ( array, dimension, data_type );
             break;
        case 11:
             (void) rexponential ( array, dimension, data_type, damp );
             break;
        case 12:
             (void) lexponential ( array, dimension, data_type, damp );
             break;
        case 13:
             (void) utriangular ( array, dimension, data_type, width );
             break;
        }/* end switch */
    } /* end loop */

/*
** call dunload to return vector to its "original" format
*/

    if(!dunload_vector((char **)temp1,image,image->data_storage_type,
        numvects,dimension,procdir)){
       fprintf(stderr,"ldwindow:  Unable to unload vector.\n");
       return(0);
    }

/*
** normal return
*/

    return(1);

} /* end of ldwindow */



/******************************************************************************
*
*
*       PURPOSE: Window routines. All window functions assume an even
*                number of samples from 0 to N-1.
*
*                Window functions available:
*
*                Rectangular
*                Bartlett
*                Hamming
*                Hanning
*                Kaiser
*                Parzen
*                Blackman
*                Exponential
*                Blackman-Harris
*                Left-Triangular
*                Right- Triangular
*                Left-Exponential
*                Right-Exponential
*                User-defined Triangular
*
*
******************************************************************************/

#define PI 3.14159265

/* bartlett window */
int
bartlett ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double temp;
 
   for ( n = 0; n < length/2; n++ ) {
       temp = (double) ( (2.0*n) / (double)length );
       if(data_type) {
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp*window[n];
       }
   }
   for ( n = length/2; n < length; n++ ) {
       temp = (double) ( (2.0*(length-n)) / (double)length );
       if (data_type){
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp * window[n];
       }
   }
   return(1);
}



/* hamming window */
int
hamming ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double temp;

   for ( n = 0; n < length; n++ ) {
       temp = (double)( 0.54-0.46*cos(2.0*PI*(double)n/(double)length) );
       if(data_type){
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp * window[n];
       }
    }
   return(1);
}


/* hanning window */
int
hanning ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double temp;

   for ( n = 0; n < length; n++ ) {
       temp = (double)( 0.5-0.50*cos(2.0*PI*(double)n/(double)length) );
       if (data_type) {
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp * window[n];
       }
   }
   return(1);
}


/* blackman window */
int
blackman ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double temp1, temp2;

   for ( n = 0; n < length; n++ ) {
       temp1 = (double)( 2.0*PI*(double)n/(double)length );
       temp2 = (double)( 0.42 - 0.50*cos(temp1) + 0.08*cos(2.0*temp1) );
       if (data_type) {
          window[n*2] = (float)temp2 * window[n*2];
          window[n*2+1] = (float)temp2 * window[n*2+1];
       }else{
          window[n] = (float)temp2 * window[n];
       }
   }
   return(1);
}


/* blackman-harris window */
int
blarris ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double term1, term2, term3, term4;

   for ( n = 0; n < length; n++ ) {
       term1 = 0.35875;
       term2 = 0.48829*cos(2.0*PI*1*(double)n/((double)length));
       term3 = 0.14128*cos(2.0*PI*2*(double)n/((double)length));
       term4 = 0.01168*cos(2.0*PI*3*(double)n/((double)length));
       if (data_type) {
          window[n*2] = (term1 - term2 + term3 - term4) * window[n*2];
          window[n*2+1] = (term1 - term2 + term3 - term4) * window[n*2+1];
       }else{
          window[n] = (term1 - term2 + term3 - term4) * window[n];
       }
   }
   return(1);
}


/* parzen window */
int
parzen ( window, length, data_type )
int length, data_type;
float *window;
{
   int n;
   double arg1, absolute();

   for ( n = 0; n < length; n++ ) {
       arg1 = ((double)n - 0.5*((double)length))/(0.5*((double)length+2));
       if(data_type){
          window[n*2] = (1.0 - absolute(arg1)) * window[n*2];
          window[n*2+1] = (1.0 - absolute(arg1)) * window[n*2+1];
       }else{
          window[n] = (1.0 - absolute(arg1)) * window[n];
       }
   }
   return(1);
}


/* kaiser window */
int
kaiser ( window, length, data_type, theta )
int length, data_type;
float *window, theta;
{
   int n;
   double fact(), bess(), arg1, arg2, arg3;

   arg1 = (double)length / 2.0;
   for ( n = 0; n < length; n++ ) {
       arg2 = (double)n - arg1;
       arg3 = sqrt(arg1*arg1-arg2*arg2);
       if (data_type) {
          window[n*2] = (bess(theta*arg3)/bess(theta*arg1))*window[n*2];
          window[n*2+1] = (bess(theta*arg3)/bess(theta*arg1))*window[n*2+1];
       }else{
          window[n] = (bess(theta*arg3)/bess(theta*arg1))*window[n];
       }
   }
   return(1);
}


/* rexponential window */
int
rexponential ( window, length, data_type, tau )
int length, data_type;
float *window, tau;
{
   int n;
   double temp;

   for (n=0;n<length;n++){
       temp = exp(-PI * tau * n );
       if (data_type){
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp*window[n];
       }
   }
   return(1);
}





/* lexponential window */
int
lexponential ( window, length, data_type, tau )
int length, data_type;
float *window, tau;
{
   int n, offset;
   double temp;

   offset = length - 1;
   for ( n = offset; n >= 0; n-- ) {
       temp = exp( -PI * tau * (offset-n) );
       if(data_type){
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp *window[n];
       }
   }
   return(1);
}



/* exponential window */
int
exponential ( window, length, data_type, tau )
int length, data_type;
float *window, tau;
{
   int n, offset;
   double temp;

   offset = length/2 - 1;
   for ( n = offset; n >= 0; n-- ) {
       temp = exp( -PI * tau * (offset-n) );
       if(data_type){
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp * window[n];
       }
   }
   for ( n = length/2; n < length; n++ ) {
       temp = exp(-PI * tau * n );
       if(data_type) {
          window[n*2] = (float)temp * window[n*2];
          window[n*2+1] = (float)temp * window[n*2+1];
       }else{
          window[n] = (float)temp * window[n];
       }
   }
   return(1);
}




/* left-triangular window */
int
ltriangular ( window, length, data_type )
int length, data_type;
float *window;
{
  int n;

  for ( n = 0; n < length; n++ ){
      if(data_type){ 
         window[n*2] = (1.0 - (double)n/(double)(length-1))*window[n*2];
         window[n*2+1] = (1.0 - (double)n/(double)(length-1))*window[n*2+1];
      }else{
         window[n] = (1.0 - (double)n/(double)(length-1))*window[n];
      }
  }
  return(1);
}



/* right-triangular window */
int
rtriangular ( window, length, data_type )
int length, data_type;
float *window;
{
  int n;
        
  for ( n = 0; n < length; n++ ){
      if(data_type) {
         window[n*2] = ((double)n/(double)(length-1))*window[n*2];
         window[n*2+1] = ((double)n/(double)(length-1))*window[n*2+1];
      }else{
         window[n] = ((double)n/(double)(length-1))*window[n];
      }
  }
  return(1);
}



/* user-defined triangular window */
int
utriangular ( window, length, data_type, center )
int length, data_type, center;
float *window;
{
  int n, i;
  double arg1;

  arg1 = (double)(length-center)/2.0;
  for ( n = 0; n < arg1; n++ ){
      if(data_type){
         window[n*2] = (double)n/arg1*window[n*2];
         window[n*2+1] = (double)n/arg1*window[n*2+1];
      }else{
         window[n] = (double)n/arg1*window[n];
      }
  }
  i = (int)arg1;
  for(n=((int)arg1+center);n<length;n++){
      if(data_type){
         window[n*2] = (double)i/arg1*window[n*2];
         window[n*2+1] = (double)i/arg1*window[n*2+1];
      }else{
         window[n] = (double)i/arg1*window[n];
      }
      i = i-1; 
  }
  return(1);
}

/* -library_code_end */
