comparison 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 ) ;
+}
+
+