annotate model/model.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

/*
    Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
    ===========================================================================

    Permission to use, copy, modify, and distribute this software without fee 
    is hereby granted for research purposes, provided that this copyright 
    notice appears in all copies and in all supporting documentation, and that 
    the software is not redistributed for any fee (except for a nominal shipping 
    charge). Anyone wanting to incorporate all or part of this software in a
    commercial product must obtain a license from the Medical Research Council.

    The MRC makes no representations about the suitability of this 
    software for any purpose.  It is provided "as is" without express or implied 
    warranty.
*/

/*
	model.c
	=====

    APU, ASP model demonstration program.


    Copyright (c), 1989  The Medical Research Council, Applied Psychology Unit.


    Author  : John Holdsworth
    Written : 22th March, 1989.

    Edited  : Mike Allerhand, 1990.

    Edited  : Roy Patterson, 1992. Changed default values from time to time.

    Edited  : MAA 10th June 1993
	      Allowed pwidth=0

    Edited  : Christian Giguere, March 1994. 
	      - Added tools for manipulating sampling rate and scaling data
	      - Added in an Entry Point for the outer/middle ear function
	      - Added in an Entry Point for the transmission line filterbank
	      - Implemented a completely new version of Meddis(). 
                Files "haircell.h, haircell.c "                             
	      - Changed option "scale_at" into "gain_at".
	          That option is now passed to newCorti() as a double instead of an integer
	          In Corti(), this meant multiplying *optr by state->scale
              - Provided facilities to print filterbank information to stderr via
                option info_afb. This supersedes option erbscale_afb.
	      - Changed some default values
              - To locate changes, search for "CG"

    Edited  : Jay Datta, November 1994.
              Added the stcrit options.

            : Jay Datta, March 1995.
              Added a function Nwidth() which returns the number
              of sample points of the nwidth section of a SAI. 
            : Silenced ulim_sas and llim_sas.  
            : Changed trise_at default to 1000 from 10000. May, 1995.
            : Changed trise_at default to 10000 from 1000. May, 1995.
            : Separated nwidth and strobe lag by creating a new
              option called stlag_ai.   May, 1995.

*/

/***************************************************************************
* This module contains:
*   Routines to access or modify string values.
*   Unit conversion routines.
*   Model entry-point functions.
****************************************************************************/

#include <string.h>
#include  <stdio.h>
#include   <math.h>

#ifdef THINK_C
#include <stdlib.h>
#endif

/* interfaces to stitch system */

#include "options.h"

#include "stitch.h"
#include "source.h"
#include  "funcs.h"
#include  "units.h"
#include   "calc.h"
#if defined( DSP32 ) || defined( PC )
#include   "oops.h"
#else
#include    "ops.h"
#endif
#include     "io.h"

/* interface to model modules */

#include "gamma_tone.h"
#include  "integrate.h"
#include   "formulae.h"              /* CG: removed haircell.h */
#include    "recurse.h"
#include     "scales.h"
#include     "spiral.h"
#include      "model.h"
#include      "corti.h"
#include      "image.h"
#include       "bank.h"

/* interface to WDF module */        /* CG */

#include "formulae_tl.h"
#include    "upsample.h"
#include     "bank_tl.h"
#include     "calc_tl.h"
#include      "meddis.h"
#include         "ear.h"

/********************* Model option parameter strings ***********************
* The strings are defined, allocated static array space, and initialized to
* default values.
****************************************************************************/

#ifndef  lint
static char *sccs_id = "@(#)model.c	1.49 John Holdsworth, Mike Allerhand, Roy Patterson, Paul Manson (MRC-APU) 12/23/92" ;
#endif

/* which model to use */

char  whichdflt[] = "none"     ,        *whichstr = whichdflt  ;
char sampledflt[] = "20000."   ,       *samplestr = sampledflt ;

/* input properties */

static char     dBdflt[] = "60."      ,           *dBstr = dBdflt     ;  /* CG */
static char   bitsdflt[] = "12"       ,          *bitstr = bitsdflt   ;
static char   swapdflt[] = "off"      ,         *swapstr = swapdflt   ;


/************************ File-format strings ******************************
* The strings are defined, allocated static array space, and initialized to
* default values.
****************************************************************************/

#define RETURN_SIZE 50l

static char      framesbuff[RETURN_SIZE] = "0",       *framesstr = framesbuff ;
static char       bytesbuff[RETURN_SIZE] = "1",   *framebytesstr = bytesbuff  ;
static char       widthbuff[RETURN_SIZE] = "1",   *framewidthstr = widthbuff  ;
static char      heightbuff[RETURN_SIZE] = "1",  *frameheightstr = heightbuff ;
static char        stepbuff[RETURN_SIZE] = "1",    *framestepstr = stepbuff   ;
static char nwidthdflt[] = "-5ms"     ,       *nwidthstr = nwidthdflt ; /* roy 23-12-92 */
static char stlagdflt[]   = "5ms"     ,       *stlagstr   = stlagdflt ; /* jay 16-05-95 */
static char   susldflt[] = "5"      ,       *suslevelstr = susldflt   ; /* jay 30-08-95 */
/************************ File-format strings *******************************
*
* The data which flows through the model is called a "file".
* This is in an array which is structured or formatted according to the way
* data is interpreted as it is pulled from each particular source.
* The array is formatted into "frames", described by the file-format strings.
* The format may be as successive points, (in which case each frame is a
* single point with unit width and height), or the format may be as successive
* 2-dimensional arrays, (in which case each frame is an array divided into
* rows and columns).
*
* The strings are defined, allocated static array space, and initialized in
* gen.c. A set of convenience routines (see below) are given to access or to
* reset the value of a string. This enables the file format to be set up on
* the fly, depending upon the format of the data passing through the source.
*
*  string:        access:     reset:          comment:
* framesstr      Frames      setFrames       num frames in file
* framebytesstr  Framebytes  setFramebytes   num bytes in frame
* framewidthstr  Framewidth  setFramewidth   num points along frame width
* frameheightst  Frameheight setFrameheight  num points along frame height
* framestepstr   Framestep   setFramestep    num points between successive frames
*
* Access routines for file-format strings:
*    Frames      returns value of "framesstr",      (initially 0)
*    Framebytes  returns value of "framebytesstr",  (initially 1)
*    Framewidth  returns value of "framewidthstr",  (initially 1)
*    Frameheight returns value of "frameheightstr", (initially 1)
*    Framestep   returns value of "framestepstr",   (initially 1)
*
****************************************************************************/

long Frames()
{
    return ( atoi( framesstr      ) ) ;
}

int Framebytes()
{
    return ( atoi( framebytesstr  ) ) ;
}

int Framewidth()
{
    return ( atoi( framewidthstr  ) ) ;
}

int Frameheight()
{
    return ( atoi( frameheightstr ) ) ;
}

int Framestep()
{
    return ( atoi( framestepstr  ) ) ;
}

int  Nwidth()
{     
      int abc = (Samples(stlagstr, Samplerate()));
      if ( (int)atoi(suslevelstr)>=4)
      /* fprintf(stderr, "Nwidth is %d\n", -abc); */
      return (-abc);
      else
      return 0;  
}


/******************* Value-to-string conversion routines. ******************/
char *ltoa( val, buffer )
long val ;
char *buffer ;
{
    (void) sprintf( buffer, "%ld", val ) ;

    return ( buffer ) ;
}

char *itoa( val, buffer )
int val ;
char *buffer ;
{
    (void) sprintf( buffer, "%d", val ) ;

    return ( buffer ) ;
}

/******************* Reset routines for file-format strings ****************
*  Convert argument (newval) to a string, and copy into the specific
*  file-format string.
*    setFrames      resets  "framesstr",      (initially 0)
*    setFramebytes  resets  "framebytesstr",  (initially 1)
*    setFramewidth  resets  "framewidthstr",  (initially 1)
*    setFrameheight resets  "frameheightstr", (initially 1)
*    setFramestep   resets  "framestepstr",   (initially 1)
*
*    updateFramebytes sets  "framebytesstr" to the total number of bytes in
*                                     the frame, (ie current width * height).
*
*    Note: updateFramebytes is automatically done at the end of each call to
*          setFramewidth or setFrameheight, to reset the "framebytesstr" to
*          account for the new frame size.
****************************************************************************/

void setFrames( newval )
long newval ;
{
    framesstr = ltoa( newval, framesbuff ) ;
    return ;
}

void setFramebytes( newval )
int newval ;
{
    framebytesstr  = itoa( newval, bytesbuff ) ;
    return ;
}

void updateFramebytes()
{
    setFramebytes( Framewidth() * Frameheight() * sizeof ( short ) ) ;
    return ;
}

void setFramewidth( newval )
int newval ;
{
    framewidthstr  = itoa( newval, widthbuff ) ;
    updateFramebytes() ;
    return ;
}

void setFrameheight( newval )
int newval ;
{
    frameheightstr = itoa( newval, heightbuff ) ;
    updateFramebytes() ;
    return ;
}

void setFramestep( newval )
int newval ;
{
    framestepstr   = itoa( newval, stepbuff ) ;
    return ;
}


static Option wavopts[] = {

 { "framebytes",   bytesbuff, &framebytesstr,   "framebytes - internal",                 SilentOption},
 { "framebytes",   bytesbuff, &framebytesstr,   "framebytes - internal",                 OutputOption},
 { "frameheight", heightbuff, &frameheightstr,  "frameheight - internal",                SilentOption},
 { "frameheight", heightbuff, &frameheightstr,  "frameheight - internal",                OutputOption},
 { "framewidth",   widthbuff, &framewidthstr,   "framewidth - internal",                 SilentOption},
 { "framewidth",   widthbuff, &framewidthstr,   "framewidth - internal",                 OutputOption},
 { "frameshift",    stepbuff, &framestepstr,    "framestep - internal",                  SilentOption},
 { "frameshift",    stepbuff, &framestepstr,    "framestep - internal",                  OutputOption},
 { "frames",      framesbuff, &framesstr,       "frames - internal",                     SilentOption},
 { "frames",      framesbuff, &framesstr,       "frames - internal",                     OutputOption},

 { "samplerate",  sampledflt, &samplestr,       "Input wave sample rate (Hz)",            InOutOption},
 { "swap_wave",     swapdflt, &swapstr,         "Swap bytes in input wave",               InOutOption},
 { "bits_wave",     bitsdflt, &bitstr,          "Significant bits in input wave",         SilentOption},
 { "dB_wave",        dBdflt,  &dBstr,           "Rel. input level (_tlf & _med only)\n",  InOutOption},  /* CG */
 { "what",         whichdflt, &whichstr,        "Type of model to use required",         SilentOption},

 (char *) 0 } ;


/***************** Tools for manipulating sampling rate *********************    * CG *
* This set of two routines allows to change and keep track of the sampling
* rate of the data which flows through the model. It is the responsability
* of each entry-point function to update the sampling rate as necessary.
*
*  routine:-       comment:-
*  Samplerate      returns the current value of the sampling rate
*  SetSamplerate   resets the sampling rate
****************************************************************************/

static double rateCache = 0.0 ;

void setSamplerate( newval )
double newval ;
{
    rateCache = newval ;
    return ;
}

double Samplerate()
{
    if( rateCache == 0 )
	rateCache = Freq( samplestr ) ;

    return( rateCache ) ;
}

/******************* Tools for checking for special strings ****************
*  OptionInt( str )
*  OptionDouble( str )
*  OptionStringsEqual( str1, str2 )     -(defined in options.c).
****************************************************************************/

double OptionDouble( str )
char *str ;
{
    if( strcmp( str, "on" ) == 0 )
	return( 1. ) ;
    else if( strcmp( str, "Not_used" ) == 0 )
	return( 0. ) ;
    else
	return( atof( str ) ) ;
}

int OptionInt( str )
char *str ;
{
    if( strcmp( str, "on" ) == 0 )
	return( 1 ) ;
    else if( strcmp( str, "Not_used" ) == 0 )
	return( 0 ) ;
    else
	return( atoi( str ) ) ;
}

int OptionLog ( str )
char *str ;
{
   if ( strcmp( str, "off" ) == 0 )
	return( 0 ) ;
    else
        return( 1 ) ;
}




/************************* Tools for scaling data **************************     * CG *
* This set of two routines allows to keep track of the scaling factor needed 
* to convert the data which flows through the model into absolute units in
* the cgs system ( 0.0002 dyne/cm2 ==> 0 dB SPL). It is the responsability of
* each entry-point function to update the correct scaling factor as necessary.
* The user must initially specify the scaling of the input wave using the
* "dB_wave" option. The default is "dB_wave=60.".
*
*  routine:-     comment:-
*  Scaling       returns a scalar to convert data to cgs units
*  SetScaling    resets the scaling factor
****************************************************************************/

static double scalingCache = 0.0 ;

void setScaling( from_rms, from_dB, to_rms, to_dB )
double from_rms, from_dB, to_rms, to_dB ;
{
    scalingCache = pow( 10., ( from_dB - to_dB ) / 20. ) * to_rms / from_rms ;
}

double Scaling()
{
    double rms = 200. ;      /* the default: 200 rms ==> 60 dB */
    double dB  =  60. ;

    if( scalingCache == 0 ) {
    
	if( strcmp( dBstr, "off" ) != 0 )
	   dB = Scalar( dBstr ) ;
	setScaling( rms, dB, 0.0002, 0. ) ;
    }

   return( scalingCache ) ;
}

/************************  entry-point functions ***************************
* Each entry-point function creates and initializes a "source object", and
* returns a pointer to this object. Each source set up by one of the routines
* below is designed to be the entry point for a specialized auditory-modelling
* process. (See also io.c for sources which read from disk).
* Source set-up functions always return a Source. They may take one or more
* sources as arguments, as a specification of the input to the process
* the source performs. The set up routine intializes the source structure with
* a pointer to a callback function, which performs the data processing when
* the source is used.
*
* Pointers to the entry-point functions are defined in the stage table,
* (see FindStage() below). The functions are called, using these pointers,
* from routine ModeledSource() (see below). The purpose of this is to call the
* functions, in the order set by the stage table, so as to initialize a chain
* of objects which can ultimately be used to execute the program.
*
* The source set-up routines are analogous to fopen(), (see stitch/source.c).
* The processing phase, (analogous to fread()), uses pull/fill/roll functions,
* and this execution is started from gen.c:main() by a call to SinkSource().
*
*   routine:-                   special subroutines:-
*   EarEntry()                  Ear()                                           * CG *
*   GenericEntry()              GenericFilterBank()
*   FilterEntry()               GenericEntry()
*   TLF_FilterEntry()           TLF_GenBank()                                   * CG *
*   FineEntry()                 EarEntry() FilterEntry() TLF_FilterEntry()      * CG *
*   EnvelopeEntry()             GenericEntry()
*   ComplexEntry()              GenericEntry()
*   PolarEntry()
*   RectifyEntry()
*   LogEntry()
*   UncompressEntry()
*   SaturateEntry()             Saturate()
*   LowpassEntry()              LowpassDataTypeSource()
*   ThresholdEntry()            Meddis()  Mfsai()
*   AdaptEntry()
*   HardEntry()                 Saturate()
*   IntegralEntry()             LowpassDataTypeSource()
*   DownSampleEntry()           DownSampleSource() BlockSampleSource()
*   SaiEntry()                  Sai()
*   SummaryEntry()              Summary()
*   NullEntry()
****************************************************************************/

/* As a means of introducing "cgm" as an alias for "fed", this NullEntry point does
   absolutely nothing to the data, but it DOES enable us to enter the queue just above
   the entry point for "fed". */

static Source NullEntry( source )
Source source ;
{
    return ( source ) ;
}


/***** upsampling parameters *****/                                     /* CG */

static char    upfdflt[] = "auto"     ,          *upfstr = upfdflt    ; /* CG */
static char cutoffdflt[] = "1.0"      ,       *cutoffstr = cutoffdflt ; /* CG */
static char  downfdflt[] = "auto"     ,        *downfstr = downfdflt  ; /* CG */

/***** ear parameters *****/                                            /* CG */
 
static char middledflt[] = "on"       ,       *middlestr = middledflt ; /* CG */
static char  lconcdflt[] = "0.90"     ,        *lconcstr =  lconcdflt ; /* CG */
static char  rconcdflt[] = "1.00"     ,        *rconcstr =  rconcdflt ; /* CG */
static char  kconcdflt[] = "0.01"     ,        *kconcstr =  kconcdflt ; /* CG */
static char  nconcdflt[] = "2"        ,        *nconcstr =  nconcdflt ; /* CG */
static char lcanaldflt[] = "2.85"     ,       *lcanalstr = lcanaldflt ; /* CG */
static char rcanaldflt[] = "0.35"     ,       *rcanalstr = rcanaldflt ; /* CG */
static char kcanaldflt[] = "0.04"     ,       *kcanalstr = kcanaldflt ; /* CG */
static char ncanaldflt[] = "4"        ,       *ncanalstr = ncanaldflt ; /* CG */

static Source EarEntry( source )                                  /* CG: new entry */
Source source ;
{
    int     upfactor, downfactor ;
    double  cutoff, gain = OptionDouble( middlestr ) ;
    DeclareNew( struct _tube_info *, concha ) ;
    DeclareNew( struct _tube_info *, canal ) ;

    if( gain ) {

      /*** set up/downsampling parameters ***/
       if( strncmp( upfstr, "auto", 4 ) == 0 ) {
	   if( strncmp( downfstr, "auto", 4 ) == 0 )
	      upfactor = MIN( AUTO, ( int ) ( ( MAXSIGNALFREQ * AUTO * 2 ) / Samplerate() + 1.0 ) ) ;
	   else 
	      upfactor = ABS( OptionInt( downfstr ) ) ;
	      downfactor = upfactor ;
       }
       else {
	   upfactor = ABS( OptionInt( upfstr ) ) ;
	   if( strncmp( downfstr, "auto", 4 ) == 0 )
	      downfactor = upfactor ; 
	   else
	      downfactor = ABS( OptionInt( downfstr ) ) ;
       }
	   
      /*** set cutoff frequency of FIR upsampling filter ***/
       cutoff = 0.5 * Samplerate() * Scalar( cutoffstr ) ;
       if( downfactor > upfactor ) {
	   upfactor = MAX( 1, upfactor ) ;
	   cutoff = cutoff * ( double ) upfactor / ( double ) downfactor ;
       }    
   
      /*** upsample input wave ***/
       if( upfactor >= 1 ) {
	   source = UpSample( source, Samplerate(), cutoff, upfactor ) ;
	   ( void ) setSamplerate( Samplerate() * upfactor ) ;
	   ( void ) setFrames( Frames() * upfactor ) ;
       }

      /*** outer/middle ear filter ***/
       concha->Nsegments = atoi( nconcstr ) ;
       concha->length = atof( lconcstr ) ;
       concha->diameter = atof( rconcstr ) * 2. ;
       concha->att_factor = atof( kconcstr ) ;
       canal->Nsegments = atoi( ncanalstr ) ;
       canal->length = atof( lcanalstr ) ;
       canal->diameter = atof( rcanalstr ) * 2. ;
       canal->att_factor = atof( kcanalstr ) ;

       source = Ear( source, Samplerate(), gain, concha, canal ) ;

       Delete( concha ) ;
       Delete( canal ) ;

      /*** downsample output data ***/
       if( downfactor >= 1 ) {
	   source = DownSampleSource( source, downfactor, 1 ) ;
	   ( void ) setSamplerate( Samplerate() / downfactor ) ;
	   ( void ) setFrames( Frames() / downfactor ) ;
       }

    }
    return ( source ) ;
}


/***** frequency scale parameters*****/

static char   qualdflt[] = "9.265"    ,         *qualstr = qualdflt   ;
static char  limitdflt[] = "24.7Hz"   ,        *limitstr = limitdflt  ;
static char     mmdflt[] = "0.89"     ,           *mmstr = mmdflt     ; /* CG */

static char interpdflt[] = "off"      ,       *interpstr = interpdflt ;

/***** filterbank parameters ****/

static char    maxdflt[] = "6000Hz"   ,          *maxstr = maxdflt    ; /* roy 23-12-92 */
static char    mindflt[] = "100Hz"    ,          *minstr = mindflt    ; /* roy 23-12-92 */
#ifdef PC
static char    dendflt[] = "off"      ,          *denstr = dendflt    ; /* roy 23-12-92 */
#else
static char    dendflt[] = "off"      ,          *denstr = dendflt    ; /* roy 23-12-92 */
#endif
static char  chansdflt[] = "75"       ,        *chansstr = chansdflt  ; /* roy 23-12-92 */
static char  audiodflt[] = "off"      ,        *audiostr = audiodflt  ; /* CG */
static char   infodflt[] = "off"      ,         *infostr = infodflt   ; /* CG */

static char filterdflt[] = "gtf"      ,       *filterstr = filterdflt ; /* CG */


/****************************************************************************
* updateFrequencies()
* Routine called from ModeledSource(), below.
* It sets the frameheight (number of filter-bank channels), and initializes
* the array of channel centre-frequencies, (called "frequencies").
* These are derived from three basic parameters:
* option-name:  string-value: default:  comment:
*  mincf_afb      minstr       220Hz    Minimum center frequency (Hz)
*  maxcf_afb      maxstr       4400Hz   Maximum center frequency (Hz)
*  dencf_afb      denstr       4.       Filter density (filters/critical band)
*
* The frameheight is calculated in the routine "NumberCenterFrequencies",
* and the result is copied into the "frameheightstr" by a call to the routine
* "setFrameheight".
* The "frequencies" are calculated in the routine "GenerateCenterFrequencies".
* (Both CenterFrequency routines are in gamma_tone.c).
*
* The frequencies are stored in a global array of center frequencies (below).
****************************************************************************/

double *frequencies ;    /* global array of centre frequencies */

static void reverseFrequencies( freqs, channels )
double *freqs ;
int channels ;
{
    int i ;
    double d ;

    for( i=0 ; i<channels/2 ; i++ ) {
	d = freqs[i] ;
	freqs[i] = freqs[channels-i] ;
	freqs[channels-i] = d ;
    }
}


static void updateFrequencies()                                     /* CG: modified */
{
    int    chans = OptionInt( chansstr ) ;
    double min, max, density = atof( denstr ) ;

    if( frequencies != (double *) 0 )
	Delete( frequencies ) ;

    SetErbParameters( Freq( limitstr ), atof( qualstr ) ) ;

    if( strncmp( filterstr, "tlf", 3 ) == 0 ) {

	if( chans == 0 && density == 0. )
	    chans = 1 ;

	min  = adjustCF( Freq( minstr ), Samplerate() ) ;
	max  = adjustCF( MAX( Freq( maxstr ), min ), Samplerate() ) ;
    }

    else {
	min  = Freq( minstr ) ;
	max  = Freq( maxstr ) ;
    }

    if( chans == 0 ) {
	frequencies = GenerateCenterFrequencies( min, max, density )   ;
	setFrameheight( NumberCenterFrequencies( min, max, density ) ) ;
    }
    else {
	frequencies = NumberedCenterFrequencies( min, max, chans )   ;
	setFrameheight( chans ) ;
    }

/*
    reverseFrequencies( Frequencies, Frameheight() ) ;
*/

    return ;
}


ChannelAxis( min, max, label )
double *min, *max ;
char **label ;
{
    *min = (int) ( ErbScale( Freq( minstr ) ) * 10. ) / 10. ;
    *max = (int) ( ErbScale( Freq( maxstr ) ) * 10. ) / 10. ;

    *label = "Center Frequency [ERBs]" ;
}

/***** GTF auditory filter parameters *****/                             /* CG */

static char  orderdflt[] = "4"        ,        *orderstr = orderdflt  ;
static char  phasedflt[] = "0"        ,        *phasestr = phasedflt  ;
static char   gaindflt[] = "4."       ,         *gainstr = gaindflt   ;
static char  floatdflt[] = "off"      ,        *floatstr = floatdflt  ;

static Source GenericEntry( source, proc )                               /* CG */
Source source ;
int (*proc)() ;
{
    int chans ;
    Source ret ;

    updateFrequencies() ;

    chans  = Frameheight() ;
    ret = GenericFilterBank( source, OptionInt( interpstr ), &chans, Samplerate(), 
                      frequencies, Erb, ErbScale, Scalar( gainstr ), OptionDouble( audiostr ),
                      atoi( orderstr ), atoi( phasestr ), atoi( bitstr ), Frames(), proc, 
                      sizeof ( DataType ), OptionInt( infostr ) ) ;      /* CG */

    setFrameheight( chans ) ;

    return ( ret ) ;
}

static Source FilterEntry( source )
Source source ;
{
    switch( OptionInt( floatstr ) ) {
#ifdef FLOAT
	case 0 :
	case 1 :
	    return ( GenericEntry( source, DoRealFilterFloatDataArray ) ) ;
#else
	case 0 :
	    return ( GenericEntry( source,     DoFilterShortDataArray ) ) ;
	case 1 :
	    return ( GenericEntry( source, DoRealFilterShortDataArray ) ) ;
#endif 
	case 2 :
	    return ( GenericEntry( source,       DoNewFilterDataArray ) ) ;
    }

    return ( source ) ;
}

/***** TLF auditory filter parameters *****/                             /* CG */

static char  tgaindflt[] = "4."       ,        *tgainstr = tgaindflt  ;  /* CG */
static char   dsatdflt[] = "5.75e-6"  ,         *dsatstr = dsatdflt   ;  /* CG */
static char   feeddflt[] = "0.99"     ,         *feedstr = feeddflt   ;  /* CG */
static char   qrefdflt[] = "2."       ,         *qrefstr = qrefdflt   ;  /* CG */ 
static char outdendflt[] = "4."       ,       *outdenstr = outdendflt ;  /* CG */
static char motiondflt[] = "vel"      ,       *motionstr = motiondflt ;  /* CG */

static Source TLF_FilterEntry( source )                                  /* CG: new entry */
Source source ;
{
    int chans ;
    int upfactor, downfactor ;
    double cutoff, dsat ;
    DeclareNew( struct _tube_info *, concha ) ;
    DeclareNew( struct _tube_info *, canal ) ;

   /*** set up/downsampling parameters ***/
    if( strncmp( upfstr, "auto", 4 ) == 0 ) {
	if( strncmp( downfstr, "auto", 4 ) == 0 )
	   upfactor = MIN( AUTO, ( int ) ( ( MAXSIGNALFREQ * AUTO * 2 ) / Samplerate() +1.0 ) ) ;
	else 
	   upfactor = ABS( OptionInt( downfstr ) ) ;
	downfactor = upfactor ;
    }
    else {
	upfactor = ABS( OptionInt( upfstr ) ) ;
	if( strncmp( downfstr, "auto", 4 ) == 0 )
	   downfactor = upfactor ; 
	else
	   downfactor = ABS( OptionInt( downfstr ) ) ;
    }
	   
   /*** set cutoff frequency of FIR upsampling filter ***/
    cutoff = 0.5 * Samplerate() * Scalar( cutoffstr ) ;
    if( downfactor > upfactor ) {
	upfactor = MAX( 1, upfactor ) ;
	cutoff = cutoff * ( double ) upfactor / ( double ) downfactor ;
    }    
  
   /*** upsample input wave ***/
    if( upfactor >= 1 ) {
	source = UpSample( source, Samplerate(), cutoff, upfactor ) ;
	( void ) setSamplerate( Samplerate() * upfactor ) ;
	( void ) setFrames( Frames() * upfactor ) ;
    }

   /*** set outer/middle ear parameters ***/
    concha->Nsegments = atoi( nconcstr ) ;
    concha->length = atof( lconcstr ) ;
    concha->diameter = atof( rconcstr ) * 2. ;
    concha->att_factor = atof( kconcstr ) ;
    canal->Nsegments = atoi( ncanalstr ) ;
    canal->length = atof( lcanalstr ) ;
    canal->diameter = atof( rcanalstr ) * 2. ;
    canal->att_factor = atof( kcanalstr ) ;

   /*** set bank parameters ***/
    updateFrequencies() ;
    SetERBscaling( atof( mmstr ) ) ;
    chans = Frameheight() ;

   /*** OHC normalization: converts from cgs units (0 dB SPL=>0.0002 dynes/cm2) to input wave units ***/
    dsat = atof( dsatstr ) / Scaling() ;

   /*** process data and downsample output ***/
    source = TLF_GenBank( source, OptionInt( interpstr ), OptionInt( infostr ), &chans, &downfactor, 
                          Samplerate(), frequencies, Scalar( tgainstr ), atof( outdenstr ),
			  atof( qrefstr ), atof( feedstr ), dsat, motionstr, concha, canal ) ;

   /*** reset sampling rate and file-format strings ***/
    ( void ) setFrameheight( chans ) ;
    ( void ) setSamplerate ( Samplerate() / downfactor ) ;  
    ( void ) setFrames( Frames() / downfactor ) ;

    Delete( concha ) ;
    Delete( canal ) ;

    return ( source ) ;
}

static Source FineEntry(source )                                 /* CG: new Entry */
Source source ;
{
    if( strncmp( filterstr, "off", 3 ) == 0 )
	source = EarEntry( source ) ;

    else {

	if( strncmp( filterstr, "tlf", 3 ) == 0 )
	    source = TLF_FilterEntry( source ) ;

	else {
       	    source = EarEntry( source ) ;
	    source = FilterEntry( source ) ;
	}
    }

    return( source ) ;
}

static Source EnvelopeEntry( source )
Source source ;
{
#ifdef FLOAT
    return ( GenericEntry( source, DoRealEnvelopeFloatDataArray ) ) ;
#else
    if( OptionInt( floatstr ) != 0 )
	return ( GenericEntry( source, DoRealEnvelopeShortDataArray ) ) ;
    else
	return ( GenericEntry( source,     DoEnvelopeShortDataArray ) ) ;
#endif
}

static Source ComplexEntry( source )
Source source ;
{
    Source ret ;

#ifdef FLOAT
    ret = GenericEntry( source, DoComplexFilterFloatDataArray ) ;
#else
    ret = GenericEntry( source, DoComplexFilterShortDataArray ) ;
#endif
    setFrameheight( Frameheight() * 2 ) ;

    return ( ret ) ;
}

/******* [fbm] [bmm] [fcp,fcr] stage: "auditory filter output" ************/

static  Option fbmopts[] = {

 {     "upfactor",    upfdflt,    &upfstr, "Upsampling factor (off, auto, int)",        SilentOption},  /* CG */
 {       "cutoff", cutoffdflt, &cutoffstr, "LP cutoff as a fraction of Nyquist freq.",  SilentOption},  /* CG */
 {   "downfactor",  downfdflt,  &downfstr, "Downsampling factor (off, auto, int)\n",    SilentOption},  /* CG */

 {   "middle_ear", middledflt, &middlestr, "Enable outer/middle ear function\n",         InOutOption},  /* CG */
 {  "nconcha_ear",  nconcdflt,  &nconcstr, "Number of segments in concha",              SilentOption},  /* CG */
 {  "lconcha_ear",  lconcdflt,  &lconcstr, "Total length of concha (cm)",               SilentOption},  /* CG */
 {  "rconcha_ear",  rconcdflt,  &rconcstr, "Radius of concha (cm)",                     SilentOption},  /* CG */
 {  "kconcha_ear",  kconcdflt,  &kconcstr, "Attn. constant of concha (1/cm)",           SilentOption},  /* CG */
 {   "ncanal_ear", ncanaldflt, &ncanalstr, "Number of segments in ear canal",           SilentOption},  /* CG */
 {   "lcanal_ear", lcanaldflt, &lcanalstr, "Length of ear canal (cm)",                  SilentOption},  /* CG */
 {   "rcanal_ear", rcanaldflt, &rcanalstr, "Radius of ear canal (cm)",                  SilentOption},  /* CG */
 {   "kcanal_ear", kcanaldflt, &kcanalstr, "Attn. constant of ear canal (1/cm)\n",      SilentOption},  /* CG */

 { "channels_afb",  chansdflt,  &chansstr, "Number of channels in filter",               InOutOption},
 {    "mincf_afb",    mindflt,    &minstr, "Minimum center frequency (Hz)",              InOutOption},
 {    "maxcf_afb",    maxdflt,    &maxstr, "Maximum center frequency (Hz)",              InOutOption},
 {    "dencf_afb",    dendflt,    &denstr, "Filter density (filters/critical band)",     InOutOption},
 {     "info_afb",  infodflt,    &infostr, "Prints filterbank information (to stderr).", InputOption},  /* CG */
 {   "interp_afb", interpdflt, &interpstr, "Levels of interpolation to apply",          OutputOption},
 {   "interp_afb", interpdflt, &interpstr, "Levels of interpolation to apply",          SilentOption},
 {"audiogram_afb",  audiodflt,  &audiostr, "Audiogram equalisation parameter",          SilentOption},  /* CG */ 
 {    "bwmin_afb",  limitdflt,  &limitstr, "Minimum filter bandwith",                    InOutOption},
 {  "quality_afb",   qualdflt,   &qualstr, "Ultimate qualtity factor of filters",        InOutOption},
 {    "mmerb_afb",     mmdflt,     &mmstr, "Length of 1 erb-rate unit along BM (mm)\n",  InOutOption},  /* CG */

 {       "filter", filterdflt, &filterstr, "Select auditory filter (gtf, tlf, off)\n",   InOutOption},  /* CG */

 {    "float_gtf",  floatdflt,  &floatstr, "Floating point filter calculations",        OutputOption},
 {    "float_gtf",  floatdflt,  &floatstr, "Floating point filter calculations",        SilentOption},
 {    "phase_gtf",  phasedflt,  &phasestr, "Phase compensation option",                  InOutOption},
 {    "order_gtf",  orderdflt,  &orderstr, "Filter order",                               InOutOption},
 {     "gain_gtf",   gaindflt,   &gainstr, "Filter output amplification\n",              InOutOption},  /* CG */

 {   "motion_tlf", motiondflt, &motionstr, "BM output motion (disp, vel)",               InOutOption},  /* CG */
 { "outdencf_tlf", outdendflt, &outdenstr, "Filter density outside display range",       InOutOption},  /* CG */
 {     "qref_tlf",   qrefdflt,   &qrefstr, "Local Q-factor of BM segment",               InOutOption},  /* CG */
 { "feedback_tlf",   feeddflt,   &feedstr, "Feedback gain of OHCs (0 to 0.999)",         InOutOption},  /* CG */
 {     "dsat_tlf",   dsatdflt,   &dsatstr, "Half-saturation displacement (cm)",          InOutOption},  /* CG */
 {     "gain_tlf",  tgaindflt,  &tgainstr, "Filter output amplification\n",              InOutOption},  /* CG */

 ( char * ) 0 } ;

static void polar_callback( state, bytes, output, end, input )
Source state ;
ByteCount *bytes ;
scomplex *output, *end, *input ;
{
    register scomplex *iptr = input ;
    register scomplex *optr = output ;
    register scomplex *eptr = end ;

    while( optr < eptr ) {
	optr->real = imB( (long) iptr->real * iptr->real + (long) iptr->imag * iptr->imag ) >> 1 ;
	optr->imag = iatan2( iptr->imag, iptr->real ) ;
	iptr++ ;
	optr++ ;
    }

    return ;
}

static Source PolarEntry( source )
Source source ;
{
    return ( NewSimpleProcessingSource( polar_callback, source, "model.c polar conversion" ) ) ;
}

static void fullwave_callback( state, bytes, buffer, end, input )
Source state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    register DataType *iptr = input ;
    register DataType *optr = buffer ;
    register DataType *eptr = end ;

    if( optr < eptr )
	do
	    if( *iptr > 0 )
		*optr++ =  *iptr++ ;
	    else
		*optr++ = -*iptr++ ;

	while( optr < eptr ) ;

    return ;
}


/* rectification */

static char   rectdflt[] =     "off",   *rectstr =   rectdflt ;

/************** [fbr] stage: "rectified filter output" ******************/

static  Option fbropts[] = {

 {      "rectify",   rectdflt,   &rectstr, "Rectify filter output",                      InOutOption},

 ( char * ) 0 } ;


static void halfwave_callback( state, bytes, buffer, end, input )
Source state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    register DataType *iptr = input ;
    register DataType *optr = buffer ;
    register DataType *eptr = end ;

    if( optr < eptr )
	do
	    if( *iptr > 0 )
		*optr++ =  *iptr++ ;
	    else
		*optr++ =0,*iptr++ ;

	while( optr < eptr ) ;

    return ;
}

static Source RectifyEntry( source )
Source source ;
{
    switch( OptionInt( rectstr ) ) {

	case 1 :
	    return ( SharingSource( NewSimpleProcessingSource( halfwave_callback, source, "model.c recification" ), source ) ) ;

	case 2:
	    return ( SharingSource( NewSimpleProcessingSource( fullwave_callback, source, "model.c recification" ), source ) ) ;

	default :
	    return ( source ) ;
    }
}


/*********** [fbc] [fec] stage: "compressed filter output" **************/

/* logaritmic compression function */

static char    logdflt[] =      "log",    *logstr =    logdflt ;
static char    powdflt[] =      "off",    *powstr =    powdflt ;

static  Option fbcopts[] = {

 {     "compress",    logdflt,    &logstr, "Compression (log, power value [0-1], off)",   InOutOption},
 ( char * ) 0 } ;


#ifdef FLOAT
static void mB_callback( state, bytes, buffer, end, input )
Pointer state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    register DataType *iptr = input  ;
    register DataType *optr = buffer ;
    register DataType *eptr = end ;
    double pw;

    if (!strcmp(logstr, "log") || !strcmp(logstr, "on"))
      ;
    else 
      {
	pw=(double)atof(logstr);
        strcpy(powstr, "on");
	if (pw > 1.0) {
          fprintf(stderr, "The value of power compression cannot be more than 1.0\n");
	  exit(-125);}
      }
    
    if( optr < eptr )
      do{
	if( *iptr  >= 1. ){
	  if (!(strcmp(powstr, "off"))){
	    *optr = log10( *iptr ) * 2000. ;} 
	  else{ 
	    *optr = (float) pow((double)*iptr, (double) pw) * 100./(pw); } 
	}
	else
	     *optr =   0. ;   /* -10.    Changed from -10,  AJD 29th August, 1995.  */
	     iptr++; optr++;  }
      while( optr < eptr ) ;

    return ;

}
#endif

static Source LogEntry( source )
Source source ;
{
    if( OptionLog( logstr ) != 0 )
#ifdef FLOAT
	return ( SharingSource( NewSimpleProcessingSource( mB_callback, source, "model.c compression" ), source ) ) ;
#else
	return ( SharingSource( milliBellShortSource( source ), source ) ) ;
#endif
    else
	return ( source ) ;
}

/*********** [fbu] [feu] stage: "uncompressed filter output" **************/

/***** uncompress data by raising to a power *****/

static char  powerdflt[] =     "off",  *powerstr =  powerdflt ;

static  Option fbuopts[] = {

 {        "power",  powerdflt,  &powerstr, "Power of Compression",                      SilentOption},

 ( char * ) 0 } ;


struct _uncompress_state { double power ; } ;

static void uncompress_callback( state, bytes, buffer, end, input )
struct _uncompress_state *state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    register DataType *iptr = input  ;
    register DataType *optr = buffer ;
    register DataType *eptr = end ;
    register double power = state->power ;

    if( optr < eptr )
	do
	    *optr++ = pow( 10., *iptr++ * power ) ;
	while( optr < eptr ) ;

    return ;
}
static Source UncompressEntry( source )
Source source ;
{
    DeclareNew( struct _uncompress_state *, state ) ;

    if( OptionDouble( powerstr ) != 0.0 ) {

	state->power  = OptionDouble( powerstr ) / 2000. ;

	return ( SharingSource (NewProcessingSource( state, uncompress_callback, stitch_free, source, "model.c compression" ), source ) ) ;
    }
    else
	return ( source ) ;
}

/*********** [fbs] [fes] stage: "saturated filter output" **************/

/***** provide for possible saturation at this point *****/

static char    satdflt[] =     "off",    *satstr =    satdflt ;

static  Option fbsopts[] = {

 {     "saturate",    satdflt,    &satstr, "Introduce Saturation of non-linearity\n",   SilentOption},

 ( char * ) 0 } ;


struct _sat_state { int saturation ; } ;

static void sat_callback( state, bytes, buffer, end, input )
struct _sat_state *state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    register DataType *iptr = input ;
    register DataType *optr = buffer ;
    register DataType *eptr = end ;

    if( optr < eptr )
	do {

	     if( *iptr > 0 )
		 *optr++ = SCALE( *iptr ) / ( SCALE( 1 ) + SCALE( *iptr ) / state->saturation ) ;
	     else
		 *optr++ = SCALE( *iptr ) / ( SCALE( 1 ) - SCALE( *iptr ) / state->saturation ) ;
	     iptr++ ;

	} while( optr < eptr ) ;

    return ;
}

static Source Saturate( source, str )
Source source ;
char *str ;
{
    DeclareNew( struct _sat_state *, state ) ;

    if( OptionInt( str ) != 0 ) {

	state->saturation = OptionInt( str ) ;

	return ( SharingSource( NewProcessingSource( (Pointer) state, sat_callback, stitch_free, source, "model.c stauration"  ), source ) ) ;
    }
    else
	return ( source ) ;
}

static Source SaturateEntry( source )
Source source ;
{
    return( Saturate( source, satstr ) ) ;
}

#if 00
/***** new processing entry *****/

static char zongdflt[] =     "off", *zongstr = zongdflt ; /* parameter default for table.c */

static Source ZongEntry( source )
Source source ;
{
    int i, chans = Frameheight() ;
    DeclareNewArray( WuState *, states, chans, "model.c for states" ) ;

    /* for all of the channels create a new wu cels object */

    for( i=0 ; i<chans ; i++ )
	states[ i ] = NewWu( zongstr ) ; /* call constructor of new wu cell */

    /* use general purpose source for processing multiplexed files with these cells */

return (source ) ;
    return ( NewMultiplexedSource( states, WuMultiplexedDataArray, DeleteWu, chans, source, "model.c wu" ) ) ;
}
#endif

/*********** [fbl] [fel] stage: "low-pass filtered filter output ***********/

/* low pass filtering */

static char lstagedflt[] =     "off", *lstagestr = lstagedflt ;
static char    lupdflt[] =   "0.5ms",    *lupstr =    lupdflt ;
static char  ldowndflt[] =       "1",  *ldownstr =  ldowndflt ;
static char  llossdflt[] =       "1",  *llossstr =  llossdflt ;
static char lvlossdflt[] =       "0", *lvlossstr = lvlossdflt ;
static char ligaindflt[] =       "1", *ligainstr = ligaindflt ;

static  Option fblopts[] = {

 {   "ligain_lpf", ligaindflt, &ligainstr, "Gain of integration stage\n",               OutputOption},
 {   "ligain_lpf", ligaindflt, &ligainstr, "Gain of integration stage\n",               SilentOption},
 {   "lvloss_lpf", lvlossdflt, &lvlossstr, "Rest level for loss",                       OutputOption},
 {   "lvloss_lpf", lvlossdflt, &lvlossstr, "Rest level for loss",                       SilentOption},
 {    "lloss_lpf",  llossdflt,  &llossstr, "Loss time constant",                        OutputOption},
 {    "lloss_lpf",  llossdflt,  &llossstr, "Loss time constant",                        SilentOption},
 {   "ltdown_lpf",  ldowndflt,  &ldownstr, "Downward time constant in ms",              OutputOption},
 {   "ltdown_lpf",  ldowndflt,  &ldownstr, "Downward time constant in ms",              SilentOption},
 {     "ltup_lpf",    lupdflt,    &lupstr, "Upward time constant in ms",                SilentOption},
 {  "lstages_lpf", lstagedflt, &lstagestr, "Stages of integration\n",                   SilentOption},

 ( char * ) 0 } ;

static Source LowpassEntry( source )
Source source ;
{
    return ( LowpassDataTypeSource( source, OptionInt( lstagestr ), Frameheight(), Samples( lupstr, Samplerate() ), Samples( ldownstr, Samplerate() ), Scalar( ligainstr ) ) ) ;
}


/*********** [nap,fbt] [fet] stage: "neural activity pattern" **************/

/***** transduction switch *****/                              /* CG */

static char  transdflt[] =    "at", *transstr = transdflt ;    /* CG */


/***** meddis inner haircell model *****/                      /* CG */

static char meddisdflt[] =    "1.", *meddisstr = meddisdflt ;
static char  fibredflt[] ="medium",  *fibrestr = fibredflt ;   /* CG */
static char thresdflt[] =     "0.",  *thresstr = thresdflt ;   /* CG */

Source Meddis( input, chans, samplerate, scalar )              /* CG: modified */
Source input ;
int chans ;
double samplerate, scalar ;
{
    char *state ;
    double coupling = 100. / pow( 10., atof( thresstr ) / 20. ) ;

   /*** convert from filterbank units to meddis units (30 dB => 1 rms ) ***/

      /* first convert from filterbank output units to cgs units (0 dB SPL => 0.0002 dynes/cm2) ***/
       ( void ) setScaling( 0.0002 / Scaling() * FILTERBANK_SCALE * 4.0, 0., 0.0002, 0. ) ;

      /* second convert from cgs units to meddis units (30 dB => 1 rms) ***/
       ( void ) setScaling( 0.0002 / Scaling(), 0., MEDDIS_RMS, MEDDIS_dB ) ;

   /*** start array of Meddis hair cells ***/
    coupling = coupling * Scaling() ;
    state = NewHaircells( chans, samplerate, scalar, coupling, fibrestr ) ;
    return ( NewProcessingSource( state, DoHaircells, CloseHaircells, input, "model.c meddis" ) ) ;
}

/***** adaptive thresholding *****/

static char    stxdflt[] =      "1.",    *stxstr =    stxdflt ; /* CG 23-03-94 */
static char   risedflt[] =  "10000.",   *risestr =   risedflt ;
static char  rapiddflt[] =     "0.6",  *rapidstr =  rapiddflt ;
static char   fastdflt[] =     "0.2",   *faststr =   fastdflt ;
static char   propdflt[] =     "0.5",   *propstr =   propdflt ;
static char    latdflt[] =     "20.",    *latstr =    latdflt ; /* roy 23-12-92 */
static char vdraindflt[] =       "5", *vdrainstr = vdraindflt ;
static char  timesdflt[] =       "2",  *timesstr =  timesdflt ;

static char   compdflt[] =     "on",    *compstr =   compdflt ; /* mha 22/1/93 */


static  Option fbtopts[] = {

 { "transduction",  transdflt,  &transstr, "Select transduction function (at, meddis, off)\n",InOutOption},  /* CG */
 
 {     "trise_at",   risedflt,   &risestr, "Threshold adaptation rate (upwards)",        InOutOption},
 {"t1recovery_at",  rapiddflt,  &rapidstr, "Initial recovery rate relative to filter",   InOutOption},
 {"t2recovery_at",   fastdflt,   &faststr, "Secondary recovery rate relative to filter", InOutOption},
 {  "propt2t1_at",   propdflt,   &propstr, "Relative height of secondary adaptation",    InOutOption},
 { "frecovery_at",    latdflt,    &latstr, "Recovery rate across frequency",             InOutOption},
 {  "reclimit_at", vdraindflt, &vdrainstr, "Limitation on recovery level",               InOutOption},
 {     "times_at",  timesdflt,  &timesstr, "Oversampling of calculation of threshold",  OutputOption},
 {     "times_at",  timesdflt,  &timesstr, "Oversampling of calculation of threshold",  SilentOption},
 {"compensate_at",   compdflt,   &compstr, "Cochlea output compensation",               SilentOption},
 {      "gain_at",    stxdflt,    &stxstr, "Adaptive thresholding output gain\n",        InOutOption},  /* CG */ 

 {    "fibre_med",  fibredflt,  &fibrestr, "Spont-rate fibre type (medium, high)",       InOutOption},  /* CG */
 {   "thresh_med", thresdflt,   &thresstr, "Haircell threshold shift (dB)",              InOutOption},  /* CG */
 {     "gain_med", meddisdflt, &meddisstr, "Meddis haircell model output gain\n",        InOutOption},  /* CG */

 ( char * ) 0 } ;

static Source ThresholdEntry( source )                           /* CG: modified */
Source source ;
{
    double density ;

    if( OptionInt( chansstr ) == 0 )
	density = OptionDouble( denstr ) ;
    else 
	density = OptionInt( chansstr ) / ( ErbScale( Freq( maxstr ) ) - ErbScale( Freq( minstr ) ) ) ;

    if( strncmp( transstr, "off", 3 ) == 0 ) 
        source = source ;

    else {

        if( strncmp( transstr, "med", 3 ) == 0 )
            source = Meddis(  source, Frameheight(), Samplerate(), OptionDouble( meddisstr ) ) ;

        else {
	   compensate = OptionDouble( compstr ) ;

	   source = NewProcessingSource(
		         newCorti(
		             Frameheight(), Samplerate(), frequencies, Erb,
			     atof(  risestr ), atof(  latstr ) * density,
			     atof( rapidstr ), atof( faststr ), atof( propstr ),
			     log10( Scalar( vdrainstr ) ) * 2000.,
			     (int) ( density * atof( latstr ) / Samplerate() + 1. ) * atoi( timesstr ),
			     OptionDouble( stxstr )
			 ), (void ( * )()) Corti, closeCorti,
			 source,
			 "model.c thresholding"
		   ) ;

         }
    }
    return ( source ) ;
}

/* modulation frequency image */

static char stepmfdflt[] =   "250Hz", *stepmfstr = stepmfdflt ;
static char  maxmfdflt[] =   "250Hz",  *maxmfstr =  maxmfdflt ;
static char  minmfdflt[] =    "50Hz",  *minmfstr =  minmfdflt ;
static char  denmfdflt[] =     "off",  *denmfstr =  denmfdflt ;

#ifndef DSP32

Source Mfsai( source, segment )                               /* CG */
Source source ;
int segment ;
{
    extern char *filestr ;
    Source mfsource = FileNameSource( filestr ) ;
    int frameheight = Frameheight() ;
    int i ;

    frequencies = GenerateCenterFrequencies( Freq( minmfstr ), Freq( maxmfstr ), OptionDouble( denmfstr ) ) ;
    setFrameheight( NumberCenterFrequencies( Freq( minmfstr ), Freq( maxmfstr ), OptionDouble( denmfstr ) ) ) ;

    mfsource = GenericEntry( mfsource, DoFilterShortDataArray ) ;
   /*  mfsource = GenericEntry( mfsource, DoNewFilterDataArray ) ; */
    mfsource = NewProcessingSource(
		newCorti(
		    Frameheight(), Samplerate(), frequencies, Erb,
		    atof( risestr ), atof( latstr ) * OptionDouble( denmfstr ),
		    atof( rapidstr ), atof( faststr ), atof( propstr ),
		    log10( Scalar( vdrainstr ) ) * 2000., /* convert to millibells */
		    (int) ( OptionDouble( denmfstr ) * atof( latstr ) / Samplerate() + 1. ) * atoi( timesstr ),
		    0
		), (void ( * )()) Corti, closeCorti,
		mfsource,
		"model.c thresholding"
	    ) ;

    return ( mfsource ) ;
}

#endif
/***** long term adaptation  - not implemented - used to sneak in modulation frequency thing *****/

static char   adapdflt[] =       "0",   *adapstr =   adapdflt ;

static Source AdaptEntry( source )
Source source ;
{
#ifndef DSP32
    if( OptionDouble( denmfstr ) > 0 )
	source = Mfsai( source ) ;
#endif
    return( source ) ;
}

/*********** [fba] [fea] stage: "adapted transduced filter output" *********/

static  Option fbaopts[] = {

 {   "mmincf_mfb",  minmfdflt,  &minmfstr, "Minimum modulation frequency (Hz)",         SilentOption},
 {   "mmaxcf_mfb",  maxmfdflt,  &maxmfstr, "Maximum modulation frequency (Hz)",         SilentOption},
 {   "mdencf_mfb",  denmfdflt,  &denmfstr, "Modulation filter density ",                SilentOption},
 {   "stepcf_mfb", stepmfdflt, &stepmfstr, "Step between images\n",                     SilentOption},

 {  "tadaptation",   adapdflt,   &adapstr, "Time constant of long term adaptation\n",   SilentOption},

 ( char * ) 0 } ;



/***** hard limiting *****/

static char   harddflt[] =     "off",   *hardstr =   harddflt ;

static Source HardEntry( source )
Source source ;
{
    return( Saturate( source, hardstr ) ) ;
}

/*********** [fbh] [feh] stage: "hard limited filter output" **************/

static  Option fbhopts[] = {

 {   "hard_limit",   harddflt,   &hardstr, "Hardlimit firing rate before integration\n",OutputOption},
 {   "hard_limit",   harddflt,   &hardstr, "Hardlimit firing rate before integration\n",SilentOption},

 ( char * ) 0 } ;


/***** [sgm,cgm,fbd,fbi] [fed,fei] stage: "integrated filter output" *******/

/***** temporal integration *****/

static char  stagedflt[] =     "off",  *stagestr =  stagedflt ;
static char     updflt[] =     "8ms",     *upstr =     updflt ;
static char   downdflt[] =       "1",   *downstr =   downdflt ;
static char   lossdflt[] =       "1",   *lossstr =   lossdflt ;
static char  vlossdflt[] =       "0",  *vlossstr =  vlossdflt ;
static char  igaindflt[] =       "1",  *igainstr =  igaindflt ;

static  Option fbiopts[] = {

 {   "stages_idt",  stagedflt,  &stagestr, "Stages of integration",                      InOutOption},
 {    "igain_idt",  igaindflt,  &igainstr, "Gain of integration stage\n",               OutputOption},
 {    "igain_idt",  igaindflt,  &igainstr, "Gain of integration stage\n",               SilentOption},

 {    "vloss_idt",  vlossdflt,  &vlossstr, "Rest level for loss",                       OutputOption},
 {    "vloss_idt",  vlossdflt,  &vlossstr, "Rest level for loss",                       SilentOption},
 {     "loss_idt",   lossdflt,   &lossstr, "Loss time constant",                        OutputOption},
 {     "loss_idt",   lossdflt,   &lossstr, "Loss time constant",                        SilentOption},
 {    "tdown_idt",   downdflt,   &downstr, "Downward time constant in ms",              OutputOption},
 {    "tdown_idt",   downdflt,   &downstr, "Downward time constant in ms",              SilentOption},
 {      "tup_idt",     updflt,     &upstr, "Low-pass filter time constant",              InOutOption},

 ( char * ) 0 } ;

static Source IntegralEntry( source )
Source source ;
{
    return ( LowpassDataTypeSource( source, OptionInt( stagestr ), Frameheight(), Samples( upstr, Samplerate() ), Samples( downstr, Samplerate() ), Scalar( igainstr ) ) ) ;
}


/* downsampling */

static char  downsdflt[] =     "off",  *downsstr =  downsdflt ;

static Option dwnopts[] = {

 { "downsample",   downsdflt, &downsstr,  "Downsampling over time",                    InOutOption},
 { "frstep_epn",   downsdflt, &downsstr,  "Downsampling over time\n",                  InOutOption},

 (char *) 0 } ;

static Source DownSampleEntry( source )
Source source ;
{
    int spacing = Samples( downsstr, Samplerate() ) ;
    int channels = Frameheight() * Framewidth() ;

    if( spacing != 0 ) {

#if defined( PC ) || defined( DSP32 )
	source = NewSegmentingSource( source, Frameheight() * Framewidth() * sizeof ( DataType ) * 4 ) ;
#endif

	if( spacing > 0 )
	    source  = DownSampleSource( source, spacing, channels ) ;
	else {
	    spacing = -spacing ;
	    source = BlockSampleSource( source, spacing, channels ) ;
	}

	setFramestep( spacing ) ;

	setFrames( Frames() / spacing ) ;
    }

    return( source ) ;
}


/* stabilised image defaults */

static char   llimdflt[] = "1.5c"     ,         *llimstr = llimdflt   ;
static char   ulimdflt[] = "24ms"     ,         *ulimstr = ulimdflt   ;

static char    saidflt[] = "1"        ,          *saistr = saidflt    ;
#ifdef PC
static char pwidthdflt[] = "20ms"     ,       *pwidthstr = pwidthdflt ;
#else
static char pwidthdflt[] = "35ms"     ,       *pwidthstr = pwidthdflt ; /* roy 23-12-92 */
#endif
/* static char nwidthdflt[] = "-5ms"     ,       *nwidthstr = nwidthdflt ; roy 23-12-92 */

static char  ltlimdflt[] = "150Hz"    ,        *ltlimstr = ltlimdflt  ;
static char  utlimdflt[] = "20Hz"      ,        *utlimstr = utlimdflt ;
static char  decaydflt[] = "30ms"     ,        *decaystr = decaydflt  ;
static char    cgmdflt[] = "2.5"     ,          *cgmstr = cgmdflt    ; /* 30ms > 1.25 roy 23-12-92 */
static char     ttdflt[] = "5"       ,           *ttstr = ttdflt     ; /* 8ms > 10 roy 23-12-92 */

static char   sustdflt[] = "0"      ,     *susthreshstr  = sustdflt   ;

static char     swdflt[] = "off"    ,            *swstr  = swdflt     ;

#ifdef PC
static char   stepdflt[] = "15.25ms"  ,         *stepstr = stepdflt    ;
#else
static char   stepdflt[] = "16ms"  ,         *stepstr = stepdflt    ;
#endif

/* spiral parameters */

static char   zerodflt[] = "off"      ,         *zerostr = zerodflt    ;
static char   formdflt[] = "arch"     ,         *formstr = formdflt    ;
static char   origdflt[] = "3.072"    ,         *origstr = origdflt    ;

static char       ArchStr[] =              "arch" ;
static char        LogStr[] =               "log" ;

static char   expodflt[] =  "off"     ,         *expostr = expodflt ;

/*********** [spl,sai] [sie] stage: "stabilized auditory image" ************/

static  Option saiopts[] = {

 {   "frstep_aid",   stepdflt,   &stepstr,  "Step between frames of the auditory image (ms)",InOutOption},
 {   "pwidth_aid", pwidthdflt, &pwidthstr,  "Positive width of auditory image  (ms)",    InOutOption},
 {   "nwidth_aid", nwidthdflt, &nwidthstr,  "Negative width of auditory image (ms)",   InOutOption},

 {     "form_spd",   formdflt,   &formstr,  "Form of spiral: Archemedian or log",        InputOption},
 { "zeroline_spd",   zerodflt,   &zerostr,  "Draw zero axis of spiral",                  InOutOption},
 {     "orig_spd",   origdflt,   &origstr,  "Spiral rotation and circuit removal\n",     InOutOption},

 {  "napdecay_ai",    cgmdflt,    &cgmstr,  "Neural activity decay rate in %/ms",        InOutOption},
 {   "stdecay_ai",     ttdflt,     &ttstr,  "Strobe-threshold decay rate in %/ms",       InOutOption},
 {"stcrit_ai",   susldflt, &suslevelstr,"Strobe criterion (1 to 5) ",  InOutOption},
 {"stlag_ai" ,  stlagdflt, &stlagstr,   "Auditory image strobe lag time (in ms)",      InOutOption},
 {   "stthresh_ai",   sustdflt, &susthreshstr, "Strobe threshold set to all points above this", SilentOption},
 {     "decay_ai",  decaydflt,  &decaystr,  "Auditory image decay time constant",      InOutOption},
 {      "stinfo_ai",     swdflt,     &swstr,  "Strobe information (off, on, filename)\n",  InOutOption},   

 {    "utrate_ai",  utlimdflt,  &utlimstr,  "Disconnected in R5.6",                      SilentOption},
 {    "ltrate_ai",  ltlimdflt,  &ltlimstr,  "Disconnected in R5.6",                      SilentOption}, 
 {     "ulim_sas",   ulimdflt,   &ulimstr,  "Upper integration limit for auditory image",   SilentOption},
 {     "llim_sas",   llimdflt,   &llimstr,  "Lower integration limit for auditory image\n", SilentOption},
 
 { "exponential_decay"  ,   expodflt,     &expostr , "Exponential trigger Decay", SilentOption},   

 ( char * ) 0 } ;

/* generate Stabilised image in image.c */

static Source SaiEntry( source )                                   /* CG */
Source source ;
{
    int framestep   = Samples( stepstr,   Samplerate() ) / Framestep() ;
    int pwidth      = Samples( pwidthstr, Samplerate() ) / Framestep() ;
    int nwidth      = Samples( nwidthstr, Samplerate() ) / Framestep() ;
    int stlag       = Samples( stlagstr,  Samplerate() ) / Framestep() ;


    /* Convert arguments for spiral drawing. See spiral.ch */

#ifndef DSP32
    if      (strncmp(formstr, ArchStr, strlen( formstr) ) == 0)    form_spl = 'A';
    else if (strncmp(formstr,  LogStr, strlen( formstr) ) == 0)    form_spl = 'L';
    else stitch_error("gensai: unknown form_spl\n");

    axis_spl      = OptionInt(    zerostr ) ;
    zero_spl      = OptionDouble( origstr ) ;
#endif

    /* Hack so that spirals never include transient info.  */ 
     if (strncmp(whichstr,"spl", strlen("spl")) == 0) 
       { nwidth=0; 
         if ( (int)atoi(suslevelstr)>=4)
            suslevelstr="3";
/*         if ( (int)atoi(suslevelstr)>=4) 
            pwidth+=stlag; 
          else
            stlag=0; */
       }
    /* Hack so that sas and sep have pwidth=max(ulim) and nwidth=min(llim) */

    if ( strcmp(whichstr, "sas") == 0 || strcmp(whichstr, "sep") == 0 ) {
	pwidth = Cycles( ulimstr, frequencies[0], Samplerate() ) ;
	nwidth = Cycles( llimstr, frequencies[0], Samplerate() ) ;
	if (nwidth > 0) nwidth=0;
    }

    /* Check ranges of pwidth and nwidth, and quit if invalid */
    if (pwidth < 0)
	stitch_error("Warning: gensai  pwidth_aid  should be positive\n");
    if (nwidth > 0)
	stitch_error("Warning: gensai  nwidth_aid  should be zero or negative\n");

    /* nwidth is given as a -ve value, but is used internally as +ve */
    nwidth = -(nwidth);

    source = Sai( source,  Frameheight(),  framestep,  pwidth,  nwidth,
		  (int)Samples(decaystr, Samplerate())/Framestep(),
		  (int)Samplerate(),  frequencies,
		  (double)atof(ttstr),
		  Samples(cgmstr,Samplerate())/Framestep(),
		  (int)Freq(utlimstr),  (int)Freq(ltlimstr), (int)atoi(suslevelstr),
                  susthreshstr, swstr, expostr, stlag);

    setFramewidth( pwidth+nwidth ) ;
    setFramestep( framestep ) ;
    setFrames( Frames()/framestep ) ;

    return( source ) ;
}

DelayAxis( min, max, label )
double *min, *max ;
char **label ;
{
    *min = -Samples( pwidthstr, Samplerate() ) / Samplerate() * 1000. ;
    *max = -Samples( nwidthstr, Samplerate() ) / Samplerate() * 1000. ;

    *label = "Integration Interval [ms]" ;
}



#if 0 /* where did this come from? */
/* table of default model options */

static char   llimdflt[] =    "3.5",    *llimstr =   llimdflt ;
static char   ulimdflt[] = "imagedurn", *ulimstr =   ulimdflt ;
static char    saidflt[] =    "0.25",    *saistr =    saidflt ;
/* static char   susldflt[] = "1",     *suslevelstr =    susldflt ;  AJD 13-3-95
   static char   sustdflt[] = "0",    *susthreshstr =    sustdflt; */
#endif



/*********** [sas] [sse] stage: "stabilized auditory spectrogram" **********/

/* summarise spectrogram by adding across rows */

static  Option sasopts[] = {

 ( char * ) 0 } ;

static Source SummaryEntry( source )                             /* CG */
Source source ;
{
    source = Summary( source, Frameheight(), Scalar( saistr ), frequencies, Samplerate(), llimstr, ulimstr ) ;

    setFramewidth( 1 ) ;

    return( source ) ;
}

#if defined(DSP32) || defined( THINK_C )
static void swab( in, out, bytes )
char *in, *out ;
int bytes ;
{
    register int i ;

    for( i=0 ; i<bytes ; i+=2 ) {
	out[i+1] = in[i] ;
	out[i] = in[i+1] ;
    }

    return ;
}
#endif

static void swab_callback( state, bytes, buffer, end, input )
Source state ;
ByteCount *bytes ;
DataType *buffer, *end, *input ;
{
    extern void swab() ;

    swab( (char *) input, (char *) buffer, *bytes ) ;

    return ;
}

static Source SwabEntry( source )
Source source ;
{
    return ( NewSimpleProcessingSource( swab_callback, source, "model.c byte swapping" ) ) ;
}

struct _test_state { struct _fillable_source parent ; Source source ; unsigned chans ; } ;

static Pointer test_callback( state, bytes, buffer )
struct _test_state *state ;
ByteCount *bytes ;
DataType *buffer ;
{
    register DataType input ;
    register DataType *optr = buffer ;
    register Pointer eptr = (Pointer) buffer + *bytes ;
    register int last = *bytes == 0 ;
    register int chan ;

    while( (Pointer) optr < eptr ) {

	input = *PullItems( state->source, 1, DataType ) ;

	for( chan=0 ; chan<state->chans ; chan++ )
		*optr++ = input ;
    }

    if( !last )
	return ( (Pointer) buffer ) ;
    else
	return ( DeleteFillableSource( state ) ) ;
}

static Source TestEntry( source )
Source source ;
{
    DeclareNew( struct _test_state *, state ) ;

    state->source = source ;
    state->chans  = Frameheight() ;

    return( SetFillableSource( state, test_callback, "test model.c" ) ) ;
}

/********************* Option defaults (application specific) ***************
* Each application may have its own set of option defaults.
* These override the generic defaults defined for each option, when the
* options table is constructed, (see constructOptions()).
* Note that options must have their full names. Abbreviations won't do here.
****************************************************************************/

/* Wave and excitation pattern modules */

#ifndef DSP32
static char zerobotstr[]   = "bottom=0" ;
static char excitestr[]    = "view=excitation" ;
static char nohiding[]     = "hiddenline=off" ;
static char downsample[]   = "downsample=10ms" ;
static char turn_on_idt[]  = "stages_idt=2" ;
static char mincfdflt[]    = "mincf_afb=50" ;
static char maxcfdflt[]    = "maxcf_afb=5000" ;
static char chandflt[]     = "channels_afb=128" ;
static char lengthdflt[]   = "length=32ms"  ;
static char greyscalestr[] = "view=greyscale" ;
#endif

char *wavdflts[] = {
#ifndef DSP32
		     "view=wave",
#endif
		      (char *)0 } ;

char *stpdflts[] = {                                             /* CG */
#ifndef DSP32
		     "view=wave",
                     "filter=off",
#endif
		      (char *)0 } ;

char *asadflts[] = {
#ifndef DSP32
		     excitestr,
		     "top=2500",               zerobotstr,
		     nohiding,
		     chandflt,                                   /* CG */
		     "transduction=off",                         /* CG */
		     turn_on_idt,                                /* CG */
		     downsample,                                 /* CG */
#endif
		      (char *)0 } ;

char *epndflts[] = {
#ifndef DSP32
		     excitestr,
		     "top=1000",                zerobotstr,      /* CG */
		     nohiding,
		     chandflt,                                   /* CG */
		     turn_on_idt,
		     downsample,
#endif
		      (char *)0 } ;

char *sepdflts[] = {
#ifndef DSP32
		     excitestr,
		     "top=200",                zerobotstr,
		     nohiding,
		     chandflt,                                   /* CG */
		     "ulim_sas=16ms",          "llim_sas=-4ms",
#endif
		      (char *)0 } ;

/* Auditory (ie landscape) modules */

char *bmmdflts[] = {
#ifndef DSP32
		     lengthdflt,
		     "top=100",                "bottom=-100",
#endif
		      (char *)0 } ;

char *napdflts[] = {
#ifndef DSP32
		     lengthdflt,
		     "top=2000",                zerobotstr,      /* CG */
		     chandflt,                                   /* CG */
		     "stages_idt=2",
		     "tup_idt=0.133ms",
		     "downsample=1p",
#endif
		      (char *)0 } ;


char *saidflts[] = {
#ifndef DSP32
		     "top=500",                zerobotstr,       /* CG */
		     "stages_idt=2",
		     "tup_idt=0.133ms",
		     "downsample=1p",
		     "frstep_epn=1p",
#endif
		      (char *)0 } ;

char *spldflts[] = {
#ifndef DSP32
		     "view=spiral",
#ifdef PC
		     "width_win=300",          "height_win=300",
#else
		     "width_win=500",          "height_win=500",
#endif
                     "box_ps=off",                 /* MAA 6/9/1995 */
                     "figurelinewidth_ps=0.25",    /* MAA 6/9/1995 */
		     "top=26",                 "bottom=25",
		     "dencf=1",                 "pen=2", 
		     "nwidth_aid=0ms",           "pwidth_aid=35ms",
                     "stlag=5ms",                "stcrit=5",  /* AJD */
#endif
		      (char *)0 } ;

/* Speech (ie greyscale) modules */

char *sgmdflts[] = {
#ifndef DSP32
		     greyscalestr,
		     "top=2500",               zerobotstr,
		     chandflt,                                   /* CG */
		     "transduction=off",                         /* CG */
		     turn_on_idt,
		     downsample,
#endif
		      (char *)0 } ;

char *cgmdflts[] = {
#ifndef DSP32
		     greyscalestr,
		     "top=1000",                zerobotstr,      /* CG */
		     chandflt,                                   /* CG */
		     turn_on_idt,
		     downsample,
#endif
		      (char *)0 } ;

char *sasdflts[] = {
#ifndef DSP32
		     greyscalestr,
		     "top=300",                zerobotstr,       /* CG */
		     chandflt,                                   /* CG */
		     "stages_idt=off",
		     "frstep_epn=off",
		     "downsample=off",
		     "ulim_sas=16ms",          "llim_sas=-4ms",
#endif
		      (char *)0 } ;



/******************************* stage table ********************************
*
* The "stages" are application names.
* A stage identifier name is the last three chars of the program name.
*    (Eg "sai" is the stage identifier for the stabilized auditory image
*     application, which is a program called "gensai". This is an alias for
*     the "gen" program, provided by one of the set of symbolic links).
* The stages are grouped into "models" (each an array of structs below):
*     envelope: calculation of the envelope of the underlying function.
*     fine:     calculation of the whole ("fine structure") of the function.
*     complex
*     nonmult
*     noncalc
*
* FindStage( which )  - find a place in the stage table for the program goal.
* Take a stage identifier, (contained in the "which" string), and return a ptr
* to the corresponding stage structure, (one line of one of the arrays below).
*
* The returned goal stage is the final processing stage of the sequence of
* stages which make up the program.
* The sequential order of processing is dictated by the stage table.
* The sequence of stages is from the first stage (texually at the bottom of
* each "model" group in the stage table), through each stage back up the table
* to the goal stage.
*
* The stage structure (defined below) is as follows:
*     struct _stage {
*         char    *ident    ;       Stage identifier name
*         Source (*entry)() ;       Pointer to function for entry to stage
*         Option  *options  ;       Model options for stage (see table.c)
*         char   **defaults ;       Option defaults for stage (see gen.c)
*         char    *help     ;       Description of stage
*     };
*
* For example, the sai stage (found when gensai is asked for) is:
*     "sai",  SaiEntry,  saiopts,  saidflts,  "stabilized auditory image",
*
* The stage struct fields are:
*
*  *ident
* Stage identifier name used by FindStage() to match a stage struct to the
* application. The stage is found by matching the ident string against those
* stored in the model arrays of stages. The default stage is "wav".
*
*  (*entry)()
* The name of the entry-point function for the stage. Most entry-point
* functions are defned in model.c. The entry-point function name is only
* accessed by a call to FindStage() from ModeledSource() in gen.c.
*
*  *options
* The options field of a stage is used by the constructOptions() routine
* to include the model options for that stage in the complete options table
* which it constructs. This is composed of the display options, (which are
* common to all applications, and are listed in array displayopts[] in gen.c),
* the model options for the stage, (which are listed in table.c), AND all the
* model options given for each of the stages back to the start stage, (at the
* bottom of the stage table, within a model group of stages).
*
*  **defaults
* The defaults field is used by the constructOptions() routine to override the
* defaults defined for the options which have been included in the complete
* options table. The reason for this is so that each application can have its
* own set of option defaults, if necessary. (If defaults are not given
* or the set is not complete, then the original defaults still apply).
* The defaults fields contain option defaults which override the generic
* defaults defined for each option.
* After the complete options table has been constructed, (constructOptions()),
* all options in the table will be:
* a) initialized to their default values (in getnopts()).
* b) re-initialized by any corresponding options in "rc" files.
* c) displayed as a help menu (if asked for).
*
*  *help
* This message is printed at the top of the help menu. The help field is only
* accessed by modelHelp(which), called from gen.c, which takes a stage
* identifier and returns the "help" string of the stage struct which has a
* matching ident field, to assign it to "helpstring".
*
*
* Its important to note that the options field refers to model options, and
* these are cumulative over the sequential stages of processing.
* The defaults field refers to the defaults to be used for the goal stage
* only, (ie the application), and these are not cumulative.
* The complete options table (constructed at run time) consists of the display
* options, and the accumulated model options (one group per processing stage).
*
* The NullEntry routine is a dummy which does no processing, and is used to
* introduce aliases into the stage sequence.
*
* Routines which call FindStage(), and access the stage table:
*
*                               main()
*                                 |
*          +----------------------+-----------------+
*          |                      |                 |
*          |                      |                 |
*  constructOptions()       ModeledSource()     modelHelp()
*          |                      |                 |
*          |                      |                 |
*          +----------------------+-----------------+
*                                 |
*                             FindStage()
*
****************************************************************************/


#define NO_OPTS     (Option *) 0
#define NO_DFLTS    (char **)  0

struct _stage *FindStage( which )
char *which ;
{
    static struct _stage envelope[] = {
	"sse",    SummaryEntry, sasopts,    sasdflts,   "stabilized auditory spectrogram",
	"sie",        SaiEntry, saiopts,    saidflts,   "stabilized auditory image",
	"sem",       NullEntry, NO_OPTS,    NO_DFLTS,   "spectrogram",
	"fed", DownSampleEntry, dwnopts,    NO_DFLTS,   "downsampled filter envelope",
	"fei",   IntegralEntry, fbiopts,    NO_DFLTS,   "integrated filter envelope",
	"feh",       HardEntry, fbhopts,    NO_DFLTS,   "hard limited filter envelope",
	"fea",      AdaptEntry, fbaopts,    NO_DFLTS,   "adapted filter envelope",
	"fet",  ThresholdEntry, fbtopts,    NO_DFLTS,   "thresholded filter envelope",
	"fel",    LowpassEntry, fblopts,    NO_DFLTS,   "low-pass filtered filter envelope",
	"fes",   SaturateEntry, fbsopts,    NO_DFLTS,   "saturated filter envelope",
	"feu", UncompressEntry, fbuopts,    NO_DFLTS,   "uncompressed filter envelope",
	"fec",   EnvelopeEntry, fbcopts,    NO_DFLTS,   "compressed filter envelope",
	"fem",       NullEntry, fbmopts,    NO_DFLTS,   "dummy entry to get options through",
	"wav",       NullEntry, wavopts,    wavdflts,   "picture of wave",
	( char * ) 0 } ;

    static struct _stage fine[] = {
	"sep",       NullEntry, NO_OPTS,    sepdflts,   "stabilized excitation pattern",
	"sas",    SummaryEntry, sasopts,    sasdflts,   "stabilized auditory spectrogram",
	"spl",       NullEntry, NO_OPTS,    spldflts,   "spiral stabilized auditory image",
	"sai",        SaiEntry, saiopts,    saidflts,   "stabilized auditory image",
	"sgm",       NullEntry, NO_OPTS,    sgmdflts,   "spectrogram",
	"cgm",       NullEntry, NO_OPTS,    cgmdflts,   "cochleogram",
	"nap",       NullEntry, NO_OPTS,    napdflts,   "neural activity pattern",
	"asa",       NullEntry, NO_OPTS,    asadflts,   "auditory spectral analysis",
	"epn",       NullEntry, NO_OPTS,    epndflts,   "excitation pattern",
	"fbd", DownSampleEntry, dwnopts,    NO_DFLTS,   "downsampled filter output",
	"fbi",   IntegralEntry, fbiopts,    NO_DFLTS,   "integrated filter output",
	"fbh",       HardEntry, fbhopts,    NO_DFLTS,   "hard limited filter output",
	"fba",      AdaptEntry, fbaopts,    NO_DFLTS,   "adapted transduced filter output",
	"fbt",  ThresholdEntry, fbtopts,    NO_DFLTS,   "thresholded filter output",
	"fbl",    LowpassEntry, fblopts,    NO_DFLTS,   "low-pass filtered filter output",
	"fbs",   SaturateEntry, fbsopts,    NO_DFLTS,   "saturated filter output",
	"fbu", UncompressEntry, fbuopts,    NO_DFLTS,   "uncompressed filter output",
	"fbc",        LogEntry, fbcopts,    NO_DFLTS,   "compressed filter output",
	"fbr",    RectifyEntry, fbropts,    NO_DFLTS,   "rectified filter output",
	"fbm",       NullEntry, NO_OPTS,    NO_DFLTS,   "auditory filter output",           /* CG */
        "stp",       FineEntry, fbmopts,    stpdflts,   "stapes vibration output",          /* CG */
	"wav",       NullEntry, wavopts,    wavdflts,   "picture of wave",
	( char * ) 0 } ;

    static struct _stage complex[] = {
	"fcp",      PolarEntry, NO_OPTS,    NO_DFLTS,   "polar co-ordinate complex output",
	"fcr",    ComplexEntry, fbmopts,    NO_DFLTS,   "rectangular co-ordinate complex output",
	"wav",       NullEntry, wavopts,    wavdflts,   "picture of wave",
	( char * ) 0 } ;

    static struct _stage nonmult[] = {
	"bmm",       NullEntry, NO_OPTS,    bmmdflts,   "basilar membrane motion",          /* CG */
        "stp",       FineEntry, fbmopts,    stpdflts,   "stapes vibration output",          /* CG */
	"wav",       NullEntry, wavopts,    wavdflts,   "draw wave",
	( char * ) 0 } ;

    static struct _stage noncalc[] = {
	"tst",       TestEntry, NO_OPTS,    NO_DFLTS,   "Test system",
	"wav",       NullEntry, wavopts,    wavdflts,   "draw wave",
	( char * ) 0 } ;

    static struct _stage *models[] = { envelope, fine, complex, nonmult, noncalc, ( struct _stage * ) 0 } ;

    struct _stage *stage ;
    int model ;

    /* for each type of model */
    for( model=0 ; models[ model ] != ( struct _stage * ) 0 ; model++ )
	/* for each stage in the model */
	for( stage=models[ model ] ; stage->ident != ( char * ) 0 ; stage++ )
	    if( strncmp( which, stage->ident, strlen( stage->ident ) ) == 0 ) {

		/* kludge to fix bug (see fixchans routine below) */
		if ( strcmp( which, "nap" ) == 0 )
		    fixchans( stage->defaults ) ;

		return ( stage ) ;
	    }

    return ( FindStage( "wav" ) ) ;
}

/* Return the "help" string for the given stage identifier */
/* This routine is called once from gen.c:main()           */
char *modelHelp( which )
char *which ;
{
    return( FindStage( which )->help ) ;
}


/*
This routine is a hack to fix a bug we caused earlier!
We changed the order of entry points in Findstage, so that nap was able
to use the lp filter. But this caused nap to inherit the defaults of all
the other entry points which use that filter. In particular, the number of
channels for nap was changed as a result of this. This hack is designed to
restore the number chans, in the case of gennap only, to their default
value.
*/

fixchans( str )
char **str ;
{
    while ( *str != (char *)0 ) {
	if ( strncmp( *str, "channels_afb", 12 ) == 0 ) {
	    strcpy( (*str)+13, chansdflt ) ;
	/*    fprintf(stderr, "fixchans: %s\n", *str ) ;        */
	    break ;
	}
	*str++ ;
    }
}


/*********************** Strings for use in DSP version ********************/

#if defined(DSPHOST) || defined(DSP32)
#if 00
char **shared[] = { &llimstr,       &ulimstr,        &saistr,       &pwidthstr,
		    &nwidthstr,     &ltlimstr,       &utlimstr,     &ttstr,
		    &cgmstr,        &decaystr,       &stepstr,      &downsstr,
		    &stagestr,      &upstr,          &downstr,      &lossstr,
		    &vlossstr,      &igainstr,       &hardstr,      &adapstr,
		    &timesstr,      &vdrainstr,      &propstr,      &faststr,
		    &latstr,        &rapidstr,       &risestr,      &stxstr,
		    &lstagestr,     &lupstr,         &ldownstr,     &llossstr,
		    &lvlossstr,     &ligainstr,      &meddisstr,    &satstr,
		    &logstr,        &rectstr,        &qualstr,      &limitstr,
		    &interpstr,     &minstr,         &maxstr,       &denstr,
		    &chansstr,      &orderstr,       &phasestr,     &gainstr,
		    &audiostr,      &floatstr,       &samplestr,    &bitstr,
		    &envstr,        &whichstr,       &framesstr,    &framebytesstr,
		    &framewidthstr, &frameheightstr, &framestepstr,
		    ( char ** ) 0 } ;
#endif
char **returned[] = { &framesstr,      &framebytesstr,  &framewidthstr,
		      &frameheightstr, &framestepstr,
		      ( char ** ) 0 } ;

#endif


/******* support for DSP32 C ******************/

#ifdef DSP32
void receiveParameters( parameters, psize )
char *parameters ;
short psize ;
{
    struct _stage *stage = FindStage( "sas" ), *sptr ;
    char *bptr ;
    int c ;

    for( sptr = stage ; sptr->ident != (char *) 0 ; sptr++ )
	;

    sptr-- ;

    for( bptr = parameters ; bptr < parameters + psize ; sptr-- )
	if( sptr->options != (Option *) 0 )
	    for( c=0 ; sptr->options[c].name != (char *) 0 ; c++ ) {
		*(sptr->options[c].value) = bptr ;
		bptr += strlen( bptr ) + 1 ;
	    }

    return ;
}
#endif

#ifdef DSPHOST
#include "utils.h"
#define  DSP_BIN_VAR     "DSPBIN"
#define  DEFAULT_DSP_BIN "c:\\bin\\dspmain"
extern char *getenv() ;

void sendParameters( which )
char *which ;
{
    register struct _stage *stage = FindStage( which ), *sptr ;
    unsigned c, bytes = 0 ;
    char *buffer, *bptr ;
    long paddr, baddr ;

    for( sptr = stage ; sptr->ident != (char *) 0 ; sptr++ )
	if( sptr->options != (Option *) 0 )
	    for( c=0 ; sptr->options[c].name != (char *) 0 ; c++ )
		bytes += strlen( *(sptr->options[c].value) ) + 1 ;

    sptr-- ;

    buffer = Allocate( bytes, "argument array model.c" ) ;

    for( bptr = buffer ; sptr >= stage ; sptr-- )
	if( sptr->options != (Option *) 0 )
	    for( c=0 ; sptr->options[c].name != (char *) 0 ; c++ ) {
		(void) strcpy( bptr, *(sptr->options[c].value) ) ;
		bptr += strlen( bptr ) + 1 ;
	    }


    paddr = find_label_name( "psize" ) ;

    dsp_dl_int( paddr, (int) ( bptr - buffer ) ) ;

    baddr = find_label_name( "params" ) ;

    while( dsp_up_long( baddr ) == 0 )
	;

    dsp_dl_bytes( deref( baddr ), (long) ( bptr - buffer ), buffer ) ;

    dsp_dl_int( paddr, 0 ) ;

    return ;
}


struct _dsp_state { unsigned long needed_addr, output_addr, errno_addr ; } ;

static void dsp_callback( state, bytes, buffer )
struct _dsp_state *state ;
ByteCount bytes ;
Pointer buffer ;
{	
    char errbuff[200] ;
	
    /* tell dsp how much you want */

    dsp_dl_int( state->needed_addr, bytes ) ;

    /* wait for it to provide it */

    while( dsp_up_int( state->needed_addr ) != 0 )
	if( dsp_up_int( state->errno_addr ) != 0 ) {
/*
	    if( dispw != ( WindowObject ) 0 ) Close( dispw ) ;
*/
	    dsp_up_bytes( find_label_name( "errstr" ), (long) sizeof ( errbuff ), errbuff ) ;
	    printf( "freemem %ld stackmax %lx\n",
		dsp_up_long( find_label_name( "freemem"  ) ),
		dsp_up_long( find_label_name( "stackmax" ) ) ) ;
	    stitch_error( errbuff, dsp_up_int( find_label_name( "errno" ) ) ) ;
     	}
    
    /* upload it into PC buffer */

    dsp_up_bytes( deref( state->output_addr ), (long) bytes, buffer ) ;

    return ;
}

Source DspSource( source )
Source source ;
{
    DeclareNew( struct _dsp_state *, state ) ;
    unsigned long addr ;
    char errbuff[200] ;
    char ***param ;
    char *binary_path = getenv( DSP_BIN_VAR ) ;

    if( binary_path == ( char * ) 0 )
	binary_path =  DEFAULT_DSP_BIN ;

    /* load dsp */

    (void) fprintf( stderr, "Downloading %s to DSP32 at i/o addres 0x%x\n", binary_path, default_addr() ) ;

    if( dsp_dl_exec( binary_path ) == 0 )
	stitch_error( "\nDSP load failed\n" ) ;

    /* start dsp */

    dsp_run() ;

    /* send parameters down to model */

    sendParameters( whichstr ) ;

    /* send down input if required */

    if( whichstr[0] != 'i' ) {
    	
	addr = find_label_name( "input" ) ;

	while( dsp_up_long( addr ) == 0 )
	    ;

	(void) fprintf( stderr, "Downloading input file of %ld points to dsp addres 0x%lx\n", Frames(), deref( addr ) ) ;

	dsp_dl_array( deref( addr ), Frames(), PullInts( source, Frames() ) ) ;
    }

    /* wait for model setup */

    state->needed_addr = find_label_name( "needed" ) ;
    state->output_addr = find_label_name( "output" ) ;
    state->errno_addr  = find_label_name( "errno"  ) ;

    addr = find_label_name( "source" ) ;

    while( dsp_up_long( addr ) == 0l && dsp_up_int( state->errno_addr ) == 0 )
	(void) fprintf( stderr, "Waiting for dsp (status:%d)\r", dsp_up_long( addr ) ) ;
      
    if( dsp_up_int( state->errno_addr ) != 0 ) {
	dsp_up_bytes( find_label_name( "errstr" ), (long) sizeof ( errbuff ), errbuff ) ;
	(void) fprintf( stderr, "dsp error: %s\n", errbuff ) ;
	stitch_error( "so there" ) ;
    }

    (void) fprintf( stderr, "DSP setup (source=0x%lx)\n", dsp_up_long( addr ) ) ;

    addr = find_label_name( "returned" ) ;

    for( param=returned ; *param != ( char ** ) 0 ; param++ ) {

	dsp_up_bytes( deref( deref( addr ) ), RETURN_SIZE, **param ) ;
	addr += 4 ;
    }
#if 0
    fprintf( stderr, "Frames %ld - %d %d %d\n\nType Return\n", Frames(), Frameheight(), Framewidth(), Framebytes() ) ;
    getchar() ;
#endif
    return ( stdAutoSource( ( Pointer ) state, dsp_callback ) ) ;
}

#endif

/***************************************************************************
* ModeledSource()
* Called from main() in gen.c to initialize a chain of objects which will
* ultimately execute the program, (see also SinkSource() in gen.c:main).
*
* The routine ModeledSource() first uses Findstage() to find the appropriate
* place in the stage table, called for by the program, (and indicated by
* "whichstr", the last three chars of the program's name).
* Then it finds the end of the table, and works back up to the stage called
* for by the program.
* At each stage, call the corresponding entry-point function, (provided it's
* not a null pointer), to create and initialize a source object.
* (Each stage in the table is a process on the way to the program as a whole.
* The order of the stage table sets the order of processing).
*
* The "check" argument is used when the "useprevious" option is on,
* (see main() and checkForFile() in gen.c), otherwise it will be a null ptr.
***************************************************************************/

Source ModeledSource( source, check )
Source source ;
Source (*check)() ;
{
    /* Find the place in the stage table called for by the program */
    register struct _stage *sptr, *stage = FindStage( whichstr ) ;
    Source stmp ;

#if !defined( DSP32 ) && !defined( THINK_C )
    if( OptionInt( swapstr ) != 0 )
	source = SwabEntry( source ) ;
#endif

#ifdef DSPHOST
    if( getenv( "DSP32" ) != ( char * ) 0 )
	return ( DspSource( source ) ) ;
    else
	(void) fprintf( stderr, "Not using DSP card\n" ) ;
#endif

    if( stage->entry == ( Source ( * )() ) 0 )
	return ( source ) ;
#if 00
    /* Set frameheight (channels) and centre-frequencies array */
    updateFrequencies() ;
#endif

#ifdef FLOAT
    source = ShortFloatSource( source ) ;
#endif

    /* Find the end of the stage table */
    for( sptr = stage ; sptr->ident != ( char * ) 0 ; sptr++ )
	;

    /* Work back up the stage table to the stage called for by the program */
    /* At each stage, call the entry-point function (if it's not null)     */
    while( sptr-- > stage ) {

	if( sptr->entry != (Source ( * )()) 0 )
	    source = sptr->entry( source ) ;

	if( check != ( Source ( * )() ) 0 ) {
	    stmp = check( sptr->ident ) ;
	    if( _SPTR( stmp ) != (struct _source *) 0 ) {

#ifdef FLOAT
		source = ShortFloatSource( stmp ) ;
#else
		source = stmp ;
#endif
		/* Set frameheight (channels) and centre-frequencies array */
		updateFrequencies() ;
	    }
	}
    }

#ifdef FLOAT
    source = FloatShortSource( source ) ;
#endif

    setFrames( Frames() ) ;

    return ( source ) ;
}