 /*
  * Khoros: $Id: signal.c,v 1.1 1991/05/10 15:41:26 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: signal.c,v 1.1 1991/05/10 15:41:26 khoros Exp $";
#endif

 /*
  * $Log: signal.c,v $
 * Revision 1.1  1991/05/10  15:41:26  khoros
 * Initial revision
 *
  */ 


/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, 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 too 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.
 *----------------------------------------------------------------------
 */

/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>	    file name: signal.c                               <<<<
   >>>>                                                       <<<<
   >>>>   description: contains routines that generate 2d     <<<<
   >>>>                periodic "images"  If either the row   <<<<
   >>>>                or column arguments are set to 1, then <<<<
   >>>>                a 1d  data set is generated.           <<<<
   >>>>                In all cases, the data is stored in    <<<<
   >>>>                float vector that must be recast as    <<<<
   >>>>                necessary.                             <<<<
   >>>>                                                       <<<<
   >>>>                This is part of vmath.a                <<<<
   >>>>                                                       <<<<
   >>>>      routines: 					      <<<<
   >>>>                                                       <<<<
   >>>> modifications: See module headers                     <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */


#include "unmcopyright.h"	 /* Copyright 1990 by UNM */
#include "vinclude.h"	

/***********************************************************************
*
*  Routine Name: FMODULO (macro)
*
*          Date: February 1990
*
*       Purpose: Implements a "floating point" mod function that can 
*                handle remainders that are non integer.
*
*         Input: x - first argument of mod
*                y - second argument of mod
*
*        Output: (float)(x mod y)
*
*    Written By: Jeremy Worley
*
* Modifications: March 1990 - Casted x and y so that they are "type
*                             independant"
*
***********************************************************************/

#define FMODULO(x,y) ((float)((double)x-(double)y*(int)((double)x/(double)y)))


/***********************************************************************
*
*  Routine Name: gsin()
*
*          Date: April 23, 1990
*
*       Purpose: Generates a 2d sinusoid data set
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                sample   - (float) sampling frequency
*                amp      - (float) amplitude of data set
*                rfreq    - (float) frequency of data set in row direction 
*                cfreq    - (float) frequency of data set in column direction 
*                rphase   - (float) phase of data set in row direction
*                cphase   - (float) phase of data set in column direction
*
*        Output: 1 - error condition
*                0 - no error
*
*    Written By: Jeremy Worley
*
* Modifications:
*   6-Dec-90 Scott Wilson - Fix code so that the sine generated is
*                           amp*sin(cfreq*x+rphase+rfreq*y+cphase) as is
*                           described by the Fourier integral.
*                           Also note that rphase and cphase are
*                           redundant.
*
***********************************************************************/

int gsin(image,rows,cols,sample,amp,rfreq,cfreq,rphase,cphase)
    float *image,sample,amp,rfreq,cfreq,rphase,cphase;
    int rows,cols;
{
    int i,j;
    float r,c;
    double a;

    if(image==NULL)return(0);
    if(sample==0.0)return(0);

    for(i=0;i<rows;i++)
      {
        r = (float)i*rfreq/sample;
        for(j=0;j<cols;j++)
          {
            c = (float)j*cfreq/sample;
            a = r+rphase+c+cphase;
            *image++ = amp*sin(a);
          }
      }

    return(1);
}


/***********************************************************************
*
*  Routine Name: gpwl()
*
*          Date: April 23, 1990
*
*       Purpose: Generates a 2d "piecewise linear" data set.  By properly
*                setting arguments to this routine, you can generate
*                impulse data (spikes), triangular waves, sawtooth waves,
*                reverse sawtooth waves, and square waves.
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                sample   - (float) sampling frequency
*                minimum  - (float) minimum value of data set.
*                maximum  - (float) maximum value of the data set.
*                rperiod  - (float) period of function in the row direction
*                cperiod  - (float) period of function in the column direction
*                rrise    - (float) rise time of function in the row direction
*                crise    - (float) rise time of function in the column 
*                           direction
*                rwidth   - (float) width (when high) in the row direction
*                cwidth   - (float) width (when high) in the column direction
*
*        Output: 1 - error condition
*                0 - no error
*
*    Written By: Jeremy Worley
*
* Modifications: Thu Oct 18 15:16:34 MDT 1990 Jeremy Worley
*			Fixed a bug that made the fall part of the pulse
*			start at a value less than maximum and drop into
*			negative values:
*				(section-rwidth) --> (section-rwidth-rrise)
*				(section-cwidth) --> (section-cwidth-crise)
*
***********************************************************************/

int gpwl(vect,rows,cols,sample,minimum,maximum,rperiod,cperiod,rrise,crise,
         rfall,cfall,rwidth,cwidth)
    float *vect,sample,minimum,maximum,rperiod,cperiod,rrise,crise,rfall,cfall,
          rwidth,cwidth;
    int rows,cols;

{
    int r,c;
    float ramp,camp,amp,section,rpart,cpart;

  amp = maximum - minimum;

  if(rows==1){
     ramp=0.0;
  }else
     ramp=amp;

  if(cols==1){
     camp=0.0;
     ramp=amp;
  }else
     camp=amp;


    for(r=0;r<rows;r++){
        for(c=0;c<cols;c++){
/*
** rows
*/
            section = (FMODULO(r,(rperiod*sample))/sample);
            if(section<rrise && rrise!=0.0){
               rpart = section/rrise;
            }else if(section<(rrise + rwidth) && rwidth!=0.0){
               rpart = 1;
            }else if(section<(rrise+rfall+rwidth) && rfall!=0.0){
               rpart = 1.0 - (section-rwidth-rrise)/rfall; 
            }else{
               rpart = 0.0;
            }
/*
** columns
*/
            section = (FMODULO(c,(cperiod*sample))/sample);
            if(section<crise && crise!=0.0){
               cpart = section/crise;
            }else if(section<(crise + cwidth) && cwidth!=0.0){
               cpart = 1;
            }else if(section<(crise+cfall+cwidth) && cfall!=0.0){
               cpart = 1.0 - (section-cwidth-crise)/cfall;
            }else{
               cpart = 0.0;
            }
            
            *vect++ = ramp*rpart + camp*cpart + minimum;

        }/* end c loop */
    } /* end r loop */

    return(1);
}


/***********************************************************************
*
*  Routine Name: gsinc()
*
*          Date: April 23, 1990
*
*       Purpose: Generates a 2d Sinc data set
*
*         Input: vect     - (*float) where the data set goes
*                rows     - (int) number of rows
*                cols     - (int) number of columns
*                sample   - (float) sampling frequency
*                amp      - (float) amplitude of signal
*                rfreq    - (float) frequency of signal in the row direction
*                cfreq    - (float) frequency of signal in the column direction
*                center   - (int) centering option: 
*                                 0 - center it on the data set (i.e. the
*                                     zero point of the data is at the midpoint
*                                     of vect)
*                                 1 - left justify it (the center of the data
*                                     is at vect[0], and you only get the 
*                                     right hand side of the sinc)
*
*        Output: 1 - error condition
*                0 - no error
*
*    Written By: Jeremy Worley
*
* Modifications:
*
***********************************************************************/

#define CENTER 0
#define LEFTJST 1

int gsinc(image,rows,cols,sample,amp,rfreq,cfreq,center)
    float *image,sample,amp,rfreq,cfreq;
    int rows,cols,center;
{
    int i,j,rshift,cshift;
    float ramp,camp,r,c;

    if(image==NULL)return(0);
    if(sample==0.0)return(0);

    if(rows==1)
       ramp=0.0;
    else
       ramp=amp/2;

    if(cols==1)
       camp=0.0;
    else
       camp=amp/2;

    if(rows==1 || cols==1){
       camp*=2;
       ramp*=2;
    }

    switch(center){
       case CENTER  : rshift = rows/2;
                      cshift = cols/2;
                      break;
       case LEFTJST : rshift = 0;
                      cshift = 0;
                      break;
    }

    for(i=0;i<rows;i++){
        r = (float)(i-rshift)*rfreq/sample;
        for(j=0;j<cols;j++){
            c = (float)(j-cshift)*cfreq/sample;
            if((r+c)!=0.0){
               *image++ = (ramp*sin(r) + camp*sin(c))/(r + c);
            }else{
               *image++ = amp;
            }
        }
    }

    return(1);
}

