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

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

/*
    SCCS VERSION 1.29 WITH 1.28 STRIPPED OUT.



    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.
*/

/***************************************************************************
*    image.c            Version with trigger at right end of sai.
*    =======
* Input data derived from a cochleagram is used to construct a stabilised
* auditory image (sai) of size `chans' * `framewidth' points.
*
* Organisation of data.
* Data is 2-dimensional: in time (rows) and frequency (columns).
*
* Input data is ordered by columns. The input array consists of blocks (columns)
* of `chans' points, each column being one time sample of `chans' channels.
* The zero'th point in each column is the lowest frequency channel.
*
* Output (sai) data is ordered by rows. The output array consists of blocks
* (rows) of `framewidth' points.
* The zero'th point in each row is the "oldest" time, as the time origin is
* on the right edge of the sai.
*
* Each sai is called a `frame', and its dimensions are `chans' rows by
* `framewidth' columns. In the code, the `framewidth' is often represented
* as imagewidth/chans, where imagewidth is the total amount of data in an
* sai frame, (ie. imagewidth=chans*framewidth).

* This module has been modified to have 5 levels of strobing.  
                                        A Jay Datta  6/03/95 
***************************************************************************/

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

#include "options.h"
#include "stitch.h"
#include "source.h"
#include "model.h"
#include "units.h"
#include "image.h"
#include "calc.h"

#define MAX_BUFFER (1l<<15)

#ifdef FLOAT
#define INPUT_SCALE 0.25
#define DECAY_SCALE 0.875
#else
#define INPUT_SHIFT 2
#define DECAY_ROUND 7
#endif
#define DECAY_SHIFT 3

#define MAXCHAN 1000

#ifndef  lint
static char *sccs_id = "@(#)image.c	1.27     J. Holdsworth, M. Allerhand (MRC-APU)  6/6/91" ;
#endif

#ifdef DSP32
static int segsize = 16 ;
#else
static int segsize = 50 ;
#endif

typedef enum { False, True } Bool ;
typedef struct Node{int index; short val; struct Node *next;} node;

static void doColumn();
static Pointer sai_callback();
static Pointer summary_callback();
static void save_callback();
static void addIn();
static void decayImage();

void output_simple_strobe_info();
void doList_strobe_info();
void output_thresh_info();
void decay_strobe_threshold();
void print_trigger_debugging_info();

void initlist();
node *getnode();
node *insertl();
void inss();

FILE *trigger_file;
FILE *infotxt_file;
FILE *thresho_file;

int locmax_searchtime[MAXCHAN]; 
int locmax_searchstart[MAXCHAN];
int val[MAXCHAN]; 
int ltime[MAXCHAN];
int initial_strobe_candidate_time[MAXCHAN];
int ttime[MAXCHAN];

static node *start, *endl, *pp,  *qq;
static char *exposwitch;

#ifdef FLOAT
  float *S;
#else
  short *S;
#endif


/***************************************************************************
*   struct _sai_state:
*       The source->info->state structure which is initialised by the
*       routine Sai, and is then the state argument for the callback
*       function: sai_callback.
***************************************************************************/
struct _sai_state {
    struct _pullable_source parent ;
    Source     source;
    void       (*trigger)();     /* trigger algorithm (originally doColumn) */
    unsigned   chans;            /* number of filter channels               */
    unsigned   framestep;        /* sai display update period               */
    DataType   *image;           /* sai array (DataType = short)            */
    unsigned   imagewidth;       /* length of sai array (framewidth*chans)  */
    unsigned   imagenwidth;      /* portion of sai array for transient info */
    unsigned   decay_time;       /* a factor of sai decay time constant     */
    unsigned   decay_count;
    ScalarType *framedecay;      /* array of factors for cochleagram frame  */
    int        *cps;             /* array of 1/centre-freq (in samples)     */
    int        *isPulse;         /* array of bools, to see if data > thres  */

    int        Stcrit;          /* int restriction */       
    int        SusThresh;        /* Level 2 user settable thresh */
    ScalarType triggerdecay;     /* prop/point decay for strobe threshold */
    int        *trigtime;        /* time of  last trigger pulse  */
    int        *trigheight;      /* height of last trigger pulse */
    int        *isStrobe;        /* array of booleans, set when pulse peak found */
    float      *thresh;          /* array of thresholds, one per channel */
    float        *tlim;            /* linear decay value */
    int         time;            /* keep track of sai   */
    int        *previnput;       /* lock till new pulse found */
    int        *def_strobe_candidate_time;
    int        *def_strobe_height;
    char       *switch_info;
    int         tlim1;
    int         stlag; /* strobe nwid  first attempt to decouple nwids */
};

/***************************************************************************
*   Sai:
*       Sai is called from SaiEntry (in module model.c).
*       Routines in this module it uses are:
*                doColumn,       (which uses:  addIn)
*                sai_callback,   (which uses:  decayImage)
*
* SaiEntry (in model.c) is the entry-point function to the stabilised image
* processing module. Its purpose is to prepare the arguments for Sai by
* converting parameter-option strings to integers, and then to call Sai.
*
* Sai is called with the following arguments (in model.c):
*    Sai( source, Frameheight(), segment, width,
*        (int) Samples( decaystr )/Framestep(), (int) Freq( samplestr ),
*        frequencies, Samples(ttstr)/Framestep(), Samples(cgmstr)/Framestep(),
*        (int) Freq("250Hz"), (int) Freq("20Hz") ) ;
*
* The arguments are as follows:
* (The option name refers to the description given by "help". The values etc
* can be found in table.c as described in model.docs).
*
*  Argument     Option name    Default  Comment
*  --------     -----------    -------  -------
*  chans        (see below)             Number of filter channels
*  framestep    segment_sai    16ms     SAI display update period
*  framewidth   duration_sai   32ms     SAI duration (ms)
*  decay        decayrate_sai  32ms     SAI decay time constant (SU/ms)  15ms NOW Default
*  samplerate   samplerate     20000.   Input wave sample rate (Hz)
*  cfs          (see below)             An array of filter centre frequencies
*  cgmdecay     napdecay_sai   16ms     NAP (ie cochleagram) decay time constant
*
* Both "chans" and "cfs" are setup from the routine "updateFrequencies" in
* model.c, because both arguments are derived from three basic parameters:
*
*  minstr      mincf_afb      220Hz    Minimum center frequency (Hz)
*  maxstr      maxcf_afb      4400Hz   Maximum center frequency (Hz)
*  denstr      dencf_afb      4.       Filter density (filters/critical band)
*
* The "chans" argument is calculated in the routine "NumberCenterFrequencies",
* and the result is copied into the "frameheightstr" by a call to the routine
* "setFrameheight".
* The "cfs" argument is calculated in the routine "GenerateCenterFrequencies".
* (Both CenterFrequency routines are in gamma_tone.c).
*
***************************************************************************/

Source Sai(source, chans, framestep, pwidth, nwidth, decay, samplerate, cfs, ttdecay, cgmdecay, tlim1, tlim2, suslevel, susthresh, switchinfo, expoSwitch,stlag)
Source source ;
int    chans, framestep, pwidth, nwidth, decay ;
int    samplerate;
double *cfs ;
double ttdecay, cgmdecay ; 
int    tlim1, tlim2;
int    suslevel;
char   *susthresh;
char   *switchinfo; 
char   *expoSwitch;
int stlag;
{
    
    DeclareNew( struct _sai_state *, state ) ;
    int i, chan, startup, segment, framewidth ;
    double pts_per_ms, lin_napdec, lin_napfac ;                 /* roy 11-8-92 roy */
    char *ThreshoFile, *InfotxtFile, *TriggerFile;

    exposwitch=(char *)calloc(5, sizeof(char));
    exposwitch=expoSwitch;

    ThreshoFile = (char *)calloc(256, sizeof(char));
    InfotxtFile = (char *)calloc(256, sizeof(char));
    TriggerFile = (char *)calloc(256, sizeof(char));

    state->switch_info = (char *)calloc(256,sizeof(char));
    state->switch_info = switchinfo;
    
    if (!strcmp(susthresh, "on"))
      state->SusThresh=0;
    else if (!strcmp(susthresh, "off"))
      state->SusThresh=1;
    else
      state->SusThresh=atoi(susthresh);

     if (!strcmp(state->switch_info, "off"))
         ;
    else if (!strcmp(state->switch_info, "on"))
        infotxt_file=stdout;
    else{
        strcpy(InfotxtFile, switchinfo);
        strcat(InfotxtFile, ".info");
        infotxt_file=fopen(InfotxtFile, "w");

        strcpy(TriggerFile, switchinfo);
        strcat(TriggerFile, ".trigger");
        trigger_file=fopen(TriggerFile, "w");

        strcpy(ThreshoFile, switchinfo);
        strcat(ThreshoFile, ".thresh");
        thresho_file=fopen(ThreshoFile, "w");
      }

    initlist(); 
    state->tlim1=0;

    framewidth = pwidth+nwidth;

    pts_per_ms = (double)samplerate/1000 ;

    state->stlag=stlag; /* NEW line 95April */
    /* Initialise members of the sai state structure */
    state->source = NewRetainingSource( source, sizeof ( DataType ) * framewidth * chans ) ;
    state->chans  = chans ;
    state->trigger = doColumn ;
    state->framestep = framestep  * chans ;
    state->imagewidth = framewidth * chans ;    /* size of sai buffer */
    state->imagenwidth = nwidth * chans ; /* size of buffer for transient info */
    state->decay_time  = decay * -log( 1. - 1. / ( 1 << DECAY_SHIFT ) ) ;
    state->decay_count = state->decay_time ;

    state->Stcrit   = suslevel;
    if (stlag==0 && state->Stcrit==5) /* was nwid  NEW */
    state->Stcrit = suslevel-1;
 
    /* state->SusThresh = susthresh; */

    if (!strcmp(exposwitch, "on")) fprintf(stderr, "Exponential trigger decay SWITCH is ON\n");
/*    if (!strcmp(exposwitch, "off")) printf("SWITCH is OFF\n"); */

    if (!strcmp(exposwitch, "on"))
      state->triggerdecay = (double) 1-((ttdecay / pts_per_ms ) /100.0) ;
    else 
      state->triggerdecay = (double) (ttdecay /pts_per_ms) / 100.0 ;

    if (!strcmp(switchinfo, "off"))
       ;
    else 
      fprintf(infotxt_file, "\ntriggerdecay = %15e\n", state->triggerdecay);  

    state->isStrobe   = NewZeroedArray( int,           state->chans, "image.c for Strobe"  ) ;
    state->isPulse    = NewZeroedArray( int,           state->chans, "image.c for trigger pulse"  ) ;
    state->trigtime   = NewZeroedArray( int,           state->chans, "image.c for trigger time" ) ;
    state->trigheight = NewZeroedArray( int,           state->chans, "image.c for trigger height" ) ;
    state->thresh     = NewZeroedArray( float,      state->chans, "image.c for thresh" ) ;
    state->tlim       = NewZeroedArray( float,           state->chans, "image.c for decay"  ) ;
    state->time       = 0;  /* Time initialized to Zero */
    state->previnput  = NewZeroedArray( int,           state->chans, "image.c for pulse lock" ) ;

 state->def_strobe_candidate_time      = NewZeroedArray( int,    state->chans, "image.c for prevtrig time" );
 state->def_strobe_height    = NewZeroedArray( int,    state->chans, "image.c for prev trig" );

    for (chan=0; chan < MAXCHAN; chan++)
        locmax_searchtime[chan]=0; 
    for (chan=0; chan < MAXCHAN; chan++)
        locmax_searchstart[chan]=0;
    for (chan=0; chan < MAXCHAN; chan++)
        val[chan]=0;
    for (chan=0; chan < MAXCHAN; chan++)
        initial_strobe_candidate_time[chan]=0;
    for (chan=0; chan < MAXCHAN; chan++)
        ltime[chan]=0;
    for (chan=0; chan < MAXCHAN; chan++)
        ttime[chan]=0;

#ifdef FLOAT
  S=(float *)calloc(state->chans, sizeof(float));
#else
  S=(short *)calloc(state->chans, sizeof(short));
#endif

    /* Declare new arrays for members of state structure */

    state->framedecay = NewArray(       ScalarType,      framewidth, "image.c for decay"   ) ;
    state->image      = NewZeroedArray( DataType, state->imagewidth, "image.c for image"   ) ;
    state->cps        = NewZeroedArray( int,           state->chans, "image.c for channel centre periods" ) ;

    /* Initialise cochleagram frame decay factors */
    if (cgmdecay > 0)    
      {
      /* 11111 roy 11-8-92  These mods change to a linear nap decay.         1111111 roy */
      /* It requires pts_per_ms, lin_napdec and lin_napfac declared just after DeclareNew above */

         lin_napdec = ((cgmdecay/100)/pts_per_ms)/pts_per_ms ; /* roy 19-8 roy */

	 /* cgmdecay is div by pts_per_ms to get back to napdecay in the original units */
	 /* When napdecay is interpreted as the % to decay per ms, and since we are  */
	 /* operating in points, we need napdec/pts_per_ms, as the dec for each pt within the ms */
	 /* The decay vector is set to scale the NAP by 1.0 at strobe point. So, in the sai, */
	 /* it falls to the left AND RISES to the right of 0 ms. **** roy **** */

         lin_napfac = 1.0 + nwidth*lin_napdec ;

	 for (i=0 ; i < framewidth ; i++)
	    {
              if ( lin_napfac > 0 ) 
                 state->framedecay[i] = SCALAR( lin_napfac ); 
              else
        	 state->framedecay[i] = SCALAR( 0 ); 
              lin_napfac -= lin_napdec ; 
	    }
       }
    else                        /* disable decay factor on zero argument */
	for (i=0 ; i < framewidth ; i++)
	    state->framedecay[i] = SCALAR( 1.0 ) ;

    /* Initialise centre-period for each channel */
    for (i=0 ; i < state->chans ; i++)
	state->cps[i] = samplerate / cfs[i];

#if defined( PC ) || defined( THINK_C )
    if( (long) ( MAX_BUFFER - (long) state->imagewidth * sizeof ( DataType ) ) < (long) ( state->chans * sizeof ( DataType ) ) )
	stitch_error( "Sorry, image larger than maximum buffer size\n" ) ;
#endif

/* for (i=0; i<state->imagewidth;i++)
 fprintf(stderr, "Image Buffer is %f, count = %d\n", *(state->image+i), (i+1)); */



    return ( SetPullableSource( state, sai_callback, "image.c stabilised image" ) ) ;
}

/****************************************************************************
*    sai_callback
*        The callback function at source->info->callback in the source
*        returned by Sai.
* A "frame" is a window over the cochleagram, of size chans * points.
* A "segment" is a stretch of the cochleagram which is segsize (ie 50) time
*  samples long, ie segsize columns of the cochleagram.
*
* This callback routine is executed once for each sai frame which is plotted.
* "*bytes" is the number of bytes in the sai, ie:
*               *bytes = state->imagewidth * sizeof(DataType)
****************************************************************************/
static Pointer sai_callback( state, bytes )
struct _sai_state *state ;
ByteCount *bytes ;
{
    register long col, point ;
    register long points = state->framestep;  /* total num points under frame   */
#if defined( PC )
    int segment = ( MAX_BUFFER - state->imagewidth * sizeof ( DataType ) ) / state->chans / sizeof ( DataType ) * state->chans ;
#else
    int segment = segsize * state->chans ;
#endif
    register DataType *input ;              /* cochleagram (input) data array */
    static   int  first=1;
    static   int  test=0;
             int chan;
             int i;
/* #if defined( PC )
    test=1;
#else
    if (state->SusThresh == 0)
       segment= 30 * state->chans ;
#endif  */

/* test++; */


    /* If size argument is zero, then by convention free space and return null */
    if( *bytes == 0 ) {
	Pull( state->source, 0 ) ;
	Delete( state->image   ) ;
	return ( DeletePullableSource( state ) ) ;
    }

    /* Initially pull enough data to buffer transient info ahead of the trigger */
    /* Note, the pull operation is not valid when the size arg is zero.  */
    /* The result would be a segmentation fault on the next pull, below. */

if (state->Stcrit==0 || state->Stcrit==1 || state->Stcrit==2 || state->Stcrit==3)
    if (first) {
	if (state->imagenwidth > 0)
	    input = PullItems(state->source, state->imagenwidth, DataType); 
	first = 0;
    } 


    /* Search frame in segment blocks for trigger-points */
    for( point=0 ; point < points ; point += segment ) {
	if( segment > points - point )  /* Finish at end of frame */
	    segment = points - point ;
	/* Get one segment of the input data */
	/* Pull 50 columns of data from the source, in buffer *input. */
	/* Keep one sai of data, from last pull, in buffer behind *input */
	input = PullItems(state->source, segment, DataType);
        
	for( col=0 ; col < segment ; col += state->chans ) {
	    /* Retard the "current" data point by transient buffer width */
	    /* (This retards the trigger by the transient time, about 5ms) */
            if (state->Stcrit==0 || state->Stcrit==1 || state->Stcrit==2 || state->Stcrit==3)
               state->trigger(input-state->imagenwidth, state) ;
            else{
	        state->trigger( input, state) ;     /* -state->imagenwidth, state ) ; */
                /*fprintf(stderr, " %d ", points)  ; */ }
	    input += state->chans ;     /* next column */
	    /* decay image  periodically, once every `decay_time' columns */
            /* Lines moved down */  
	    if( --state->decay_count <= 0 ) {
		decayImage( state ) ;
		state->decay_count = state->decay_time ;
	      }
	 } /* for each col */
      }    /* for each seg */

    /* whatever the number of bytes requested, the image is always this size */
     *bytes = state->imagewidth * sizeof ( *state->image ) ;

/* if (state->Stcrit==1 && state->SusThresh==0){
                if (test%5==0)
                   decayImage( state );
                test++;}
            else if (state->Stcrit==1 && state->SusThresh>0){  changed 2 to 3 
                 if (test%10==0)
                    decayImage( state );
                 test++;} 
            else{ */

    return ( (Pointer) state->image ) ; 
 }

/****************************************************************************
*   decayImage
* Periodically attenuate image to give computationally efficient image decay.
* Attenuate all points in the sai by a constant factor.
****************************************************************************/
static void decayImage( state )
struct _sai_state *state ;
{
    register DataType *image_ptr, *image_end ;
#ifdef FLOAT
    register FLOAT image_decay = DECAY_SCALE ;
#endif

    image_end = state->image + state->imagewidth ;  /* end of sai array */
#ifdef FLOAT
    for( image_ptr=state->image ; image_ptr < image_end ;   )
	*image_ptr++ *= image_decay ;
#else
    for( image_ptr=state->image ; image_ptr < image_end ; image_ptr++ )
	*image_ptr   -= *image_ptr + DECAY_ROUND >> DECAY_SHIFT ;
#endif
    return ;
}

/****************************************************************************
*   doColumn
*       Trigger algorithm for Sai.
*       Installed as source->info->state->trigger.
*
****************************************************************************/
static void doColumn( input, state )
DataType *input ;
struct _sai_state *state ;
{
  
  register int chan ;
  float strobelag;
  int s_lag;

  
  short *STI;
  static int   check=0;
  int storev;
  char buf[2];
  short stipts=-2000;
  short zeroval=0;

  STI=(short *)calloc(state->chans, sizeof(short));

  strobelag = (double) (state->time)-(state->stlag);
  
  s_lag=(state->time)-(state->stlag);
    
  /* For each channel up column, add a row of data to sai if the      */
  /* current point is non-zero. Add the row from the current point    */
  /* plus state->imagenwidth, to display the transient (nwidth) info. */
  

/* 000000000000000000000000 stcrit 0000000000000000000000000000000000 */
/* Add on every point : Low Pass Filter */

if (state->Stcrit==0)
  {          
    for( chan=0 ; chan < state->chans ; chan++ )
      {
	addIn(input+(state->imagenwidth), chan, state); 
        output_simple_strobe_info(state, stipts);
      }
  } /* if (state->Stcrit==0)  */ 


/* 111111111111111111111111 stcrit 1111111111111111111111111111111111 */
/* Add on every point greater than "stthresh_ai" (default 0) */
       
if (state->Stcrit==1)
  {                          
    for( chan=0 ; chan < state->chans ; chan++ )
      {   
	if (input[chan] > state->SusThresh )
	  {   
	    addIn(input+(state->imagenwidth), chan, state);
	    output_simple_strobe_info(state, stipts);
	  }
	else if (input[chan] <=state->SusThresh)
	  output_simple_strobe_info(state, zeroval);
      }
  } /* if (state->Stcrit==1)  */
 

/* 222222222222222222222222 stcrit 2222222222222222222222222222222222 */
/* Add on every Nap pulse peak  */ 

if (state->Stcrit==2)
  {
    for( chan=0 ; chan < state->chans ; chan++ )
      {
	if (state->isStrobe[chan])  /* Check to see if it is time to addIn */ 
	  {
	    addIn(input+((state->imagenwidth)-2*(state->chans)), chan, state); 
	    state->isStrobe[chan]=0;    /* Reset strobe flag */
            doList_strobe_info(state,  chan);
	  }
	
/* Find a pulse peak routine */

	if (!state->isPulse[chan])   /* Not in pulse, looking for pulse */
	  {
	    state->thresh[chan]=input[chan];
	    if (input[chan] > 0 /*state->thresh[chan] */ ) /* Start of a Nap pulse */
	      state->isPulse[chan] =1;  /* set in pulse flag */ 
	  }
	else                        /* In pulse, looking for peak */ 
	  {
	    if (input[chan] < state->thresh[chan]) /* peak found */
	      {
		state->isPulse[chan]=0;  
		state->isStrobe[chan]=1; /* set strobe pending flag */
	      }
	    else
	      state->thresh[chan]=input[chan];
	  }
	output_thresh_info(state, chan);
      }
  }        /* if (state->stcrit==2) */



/* 333333333333333333333333 stcrit 3333333333333333333333333333333333 */
/* Add on pulse peaks that exceed strobe threshold  (temporal shadow ) */
/* Primitive Local Max algorithm */

if (state->Stcrit==3)  
  {
    for( chan=0 ; chan < state->chans ; chan++ )
      {
	if (state->isStrobe[chan])  /* Check to see if it is time to addIn  */
	  {
	    addIn(input+((state->imagenwidth)-2*(state->chans)), chan, state); 
	    state->isStrobe[chan]   = 0;  /* Reset strobe flag */ 
            doList_strobe_info(state, chan);
	  }

/* Find a pulse peak routine */

	if (!state->isPulse[chan])    /* Not in Pulse, looking for Pulse */
	  {
	    if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan])     
	      {                             /* Previnput check ensures threshold decay 
                                               does not cut into pulse */     

		state->isPulse[chan]  = 1;            /* set flag: in pulse */
		state->thresh[chan]   = input[chan];  /* start search for peak */
	      }
            else {
	      decay_strobe_threshold(state, chan);
	      if (input[chan]==0)    /* in the zeroes between nap pulses */
		state->previnput[chan] = 0;  /* Set it to Zero previnput */  
	    }
	  }

	if (state->isPulse[chan])  /* Else in ongoing pulse looking for peak */
	  {
	    if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan]) 
                                                       /*  pulse peak found */
	      {	      
		state->isPulse[chan]  = 0;       /* clear flag: not in pulse */
		state->isStrobe[chan] = 1;       /* set Strobe pending flag  */
                if (!strcmp(exposwitch, "off")) 
		  state->tlim[chan]=state->triggerdecay * state->thresh[chan];
                state->trigtime[chan]=(state->time)-1;
                /* state->previnput[chan]=input[cha */
	      }
	    else		    /* still in pulse, so continue search for peak */
	        state->thresh[chan] = input[chan]; 

	  }        /* if in Pulse */     
	state->previnput[chan]=input[chan];
	output_thresh_info(state, chan);
      }                /* for all chans... */
  }                 /* if (state->Stcrit==3) */



/* 444444444444444444444444 stcrit 4444444444444444444444444444444444 */
/* Add on pulse peaks that exceed strobe threshold and which are not succeeded by 
   a larger pulse in nwid ms.  (Better Local Max mechanism) */

if (state->Stcrit==4) 
  {
    for( chan=0 ; chan < state->chans ; chan++ )
      {
	if (state->isStrobe[chan])  /* Check to see if it is time to addIn  */
	  {
	    if ( strobelag-1 >= state->trigtime[chan])  
                 /* initiates strobe after nwid ms */
	      {
		addIn(input-(state->chans), chan, state);
                doList_strobe_info(state, chan);
		state->isStrobe[chan]   = 0;    /* Reset strobe flag */
		state->trigheight[chan] = 0;    /* Reset local max value */
		
	      }                   
	  }

/* Find a pulse peak routine */

	if (!state->isPulse[chan])    /* Not in Pulse, looking for pulse */
	  {
	    if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan])   
	      {                            /* Previnput check ensures threshold
                                              decay does not cut into pulse    */

		state->isPulse[chan]  = 1;            /* set flag: in pulse */
		state->thresh[chan]   = input[chan];  /* start search for peak */
          /* state->previnput[chan]= 1;     in pulse till the next zero value */
	      }
	    else 
	      {
		decay_strobe_threshold(state, chan);
                if (input[chan]==0) /* in the zeroes between nap pulses */ 
		  state->previnput[chan] = 0;  /* Free previnput */
 	      }                             
	  }

	if (state->isPulse[chan])  /* Else in ongoing pulse, looking for peak */
	  {
	    if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan])    
                                                      /* pulse peak found */
	      {
		state->isPulse[chan]  = 0;       /* clear flag: not in pulse */
		state->isStrobe[chan] = 1;       /* set Strobe pending flag  */
		
		if ( state->thresh[chan] > state->trigheight[chan] ) 
                            /* Check for LOCAL MAX */
		  {
		    state->trigheight[chan] = state->thresh[chan];
		    state->trigtime[chan] = (state->time)-1;  /* Set strobe time */
                    if (!strcmp(exposwitch, "off"))
		      state->tlim[chan]=state->triggerdecay * state->trigheight[chan];
		  }
	      }
	    else	     /* still in pulse, so continue search for peak */
	      state->thresh[chan] = input[chan]; 
	  }   /* if in Pulse */     
        state->previnput[chan]=input[chan];
	output_thresh_info(state, chan);	
  }             /* for all chans... */
}            /* if (state->Stcrit==4) */



/* 555555555555555555555555 stcrit 5555555555555555555555555555555555 */
/* Add on pulse peaks that exceed storbe threshold and which are not succeeded by
   a larger pulse in nwid ms, provided total lag is < 2*nwid ms  
   Local Max mechanism with a timeout                            */

if (state->Stcrit==5) 
 { 
  for (chan = 0;  chan < state->chans;  chan++) 
    {
     if(state->isStrobe[chan])   /* Check to see if it is time to addIn */
       {    /* check for waiting till the first local max found */
	 if (strobelag-1 >= state->def_strobe_candidate_time[chan])
	   /* initiate strobe after nwid ms  */
	   {  
	     /* print_trigger_debugging_info(state, chan, s_lag, 1); */
	     addIn(input-state->chans, chan, state);
	     doList_strobe_info(state, chan);	       
	     
	     state->def_strobe_candidate_time[chan] = 0; /* Reset strobe time */
	     state->def_strobe_height[chan] = 0; /* Reset local max value */
	     state->trigheight[chan]=0; 
	     state->isStrobe[chan]=0;       /* Reset strobe flag */
	     
	   }
       }
     
/* Find a pulse peak routine */

    if (!state->isPulse[chan])  /* Not in pulse, looking for pulse */
      { 
	if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan]) 
	  {                            /* Previnput check ensures threshold
					  decay does not cut into pulse */
	    
	    state->isPulse[chan]   = 1;           /* set flag: in pulse */
	    state->thresh[chan]    = input[chan]; /* start search for peak */
       /* state->previnput[chan] = 1;     in pulse till the next zero value */
	  }
        else
          {
	    decay_strobe_threshold(state, chan);
	    if (input[chan]==0)  /* in the zeroes between nap pulses */
	      state->previnput[chan]=0; /* Free previnput */
	  }

      }

    if (state->isPulse[chan]) /* Else in ongoing pulse, looking for peak */ 
     {  
       if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan])  
	 /* Pulse Peak Found */
	 {
	   state->isPulse[chan]=0;  /* clear flag: not in pulse */
	   
	   if (locmax_searchstart[chan]==0  )
           /* && state->time>=((ttime[chan]-initial_strobe_candidate_time[chan])+ltime[chan]))  AJD 1-2-96 */
	     {                              /* start local max time (for timeout) 
					       and shift search window */
	       locmax_searchtime[chan]=1;   /* local max time started */
	       locmax_searchstart[chan]=1;  /* local max flag set */
	       print_trigger_debugging_info(state, chan, s_lag, 2);
	       initial_strobe_candidate_time[chan]=(state->time)-1;   /* local max start time noted */
	       
	     }
	   
	   /* if (locmax_searchstart[chan]==1) */   /* 8rd March, 1995  AJayDatta */
	   if (state->thresh[chan] > state->trigheight[chan])
	     /* Check for LOCAL MAX */
	     {
	       state->trigheight[chan]=state->thresh[chan];
	       state->trigtime[chan]=(state->time)-1; /* Set strobe time */
	       if (!strcmp(exposwitch, "off"))
		 state->tlim[chan]=state->triggerdecay * state->trigheight[chan];
	       if (locmax_searchstart[chan]==1)
		 print_trigger_debugging_info(state, chan, s_lag, 3);	  
	     }       /* if greater than thresh */    
	    

    /* Reached end of time-out period  (While in the peak finding stage !)  
       End of the timeout period can occur at two stages of the algorithm,
       at a Nap pulse peak or at any point of the pulse train.            */ 
        
 	  if (locmax_searchtime[chan]==(state->stlag)) /* imagenwidth/state->chans)) */
	    {
	      print_trigger_debugging_info(state, chan, s_lag, 4);

	      state->def_strobe_candidate_time[chan]=state->trigtime[chan];
	      state->def_strobe_height[chan]=state->trigheight[chan];
	      locmax_searchtime[chan]=0;
	      locmax_searchstart[chan]=0; 
	      state->trigheight[chan]=0; 
	      state->trigtime[chan]=0;   
	      state->isStrobe[chan]=1;  /* set Strobe pending flag */
	      ltime[chan]=state->time;  /* End of search time */
	      ttime[chan]=state->def_strobe_candidate_time[chan]; /* Local Max time */ 
	    }
	  /* reset search period to zero */
	  /* as y ms search period ends, update the info for strobe and 
	     clear current_ values for the next y ms */
	}   /* Pulse Peak Loop  */ 
      else
	state->thresh[chan]=input[chan];
    } /* if (state->isPulse..) */
    
    if (locmax_searchtime[chan] == (state->stlag)) /* imagenwidth/state->chans)) */ 
    /* && state->def_strobe_candidate_time[chan]==0) */
      {                              /* start local max time (for timeout) */  
	print_trigger_debugging_info(state, chan, s_lag, 5);

	state->def_strobe_candidate_time[chan]=state->trigtime[chan];
	state->def_strobe_height[chan]=state->trigheight[chan];
	locmax_searchtime[chan]=0;
	locmax_searchstart[chan]=0;
        state->trigheight[chan]=0; 
        state->trigtime[chan]=0;   
	state->isStrobe[chan]=1;  /* set Strobe pending flag */
	ltime[chan]=state->time;
	ttime[chan]=state->def_strobe_candidate_time[chan];
      }
    /* reset search period to zero */
    /* as y ms search period ends, update the info for strobe and 
       clear current_ values for the next y ms */
    
 
    if (locmax_searchstart[chan]==1)
      (locmax_searchtime[chan])++;
  state->previnput[chan]=input[chan];
  output_thresh_info(state, chan);

  } /* for all chans */
} /* if (state->Stcrit== 5) */

  
  qq=insertl(state->time); /* insert time to list */ 
  (state->time)++;  
  return;
  
}  /* doColumn */ 


/****************************************************************************
*  addIn
*    add row cochleagram into stabilised image
****************************************************************************/
static void addIn( input, chan, state )
     DataType *input ;
     int chan ;
     struct _sai_state *state ;
{
  register DataType *channel_ptr, *image_ptr, *end;
#ifdef FLOAT
  register FLOAT input_scale = INPUT_SCALE ;
#endif
  register ScalarType *decayfactor = state->framedecay;
  
  /* Initialize channel_ptr to the input data to be added in.  */
  /* `input' points to a particular column in the cochleagram, */
  /* and `chan' indexes a particular channel in this column.   */
  channel_ptr = input + chan ;
  end = channel_ptr - state->imagewidth ;
  
  /* Initialize image_ptr to end of sai row corresponding to `chan'.*/
  /* `state->image' points to the start of the sai,                 */
  /* so increment pointer by chan*framewidth to get to row `chan'.  */
  /* To get to end of row, increment by ((chan+1)*framewidth)-1.    */
  /* But  framewidth = state->imagewidth / state->chans, and hence: */
  /* Added fix; as new spiral.c accepts nwid  AJD 17-3-95 */
  if (state->Stcrit==4 && state->stlag==0)
  image_ptr = (state->image + state->imagewidth / state->chans * (chan+1));
                         /* 13-3-95 removed -1; ajd */
  else
  image_ptr = (state->image + state->imagewidth / state->chans * (chan+1))-1; 
  
  /* Decrement channel_ptr by columns, from the initial column,    */
  /* until one complete sai row has been added into, which is when */
  /* channel_ptr has decreased by a whole imagewidth.              */
  while( channel_ptr > end ) {
#ifdef FLOAT
    *image_ptr +=        ( ( *channel_ptr *  input_scale ) * (*decayfactor++) ) ;
#else
    *image_ptr += DESCALE( ( *channel_ptr >> INPUT_SHIFT ) * (*decayfactor++) ) ;
#endif
    /* if (state->chans <= 40) */
    if (*image_ptr > 32767.0) 
      *image_ptr = 32767.0;    /* 32767 */  
    
    image_ptr--; 
    
    channel_ptr -= state->chans ;   /* next (ie previous) input time point */
  }
  
  return ;
}


/***************************************************************************
 *        Summary is not called in this module.
 *        (see SummaryEntry in model.c, which is the entry point for gensas).
 *        Summary computes a row summary spectrogram and is called during the
 *        program "gensas".
 
 *        Routines in this module it uses are:
 *                summary_callback
 ****************************************************************************/

struct  _summary_state {
  struct _fillable_source parent ;
  Source source;
  int    framewidth;          /* number of time-points or cols in an sai */
  int    frameheight;         /* number of channels or rows in an sai */
  double scale;
  int   *llim;  /* integration limits (lower & upper) for sai summary data */
  int   *ulim;
};


Source Summary( source, frameheight, scale, cfs, samplerate, llimstr, ulimstr)
Source  source ;
int     frameheight ;
double  scale ;
double *cfs, samplerate ;       /* array of channel centre frequencies */
char   *llimstr, *ulimstr;      /* lower and upper limit strings */
{
    int    i, framewidth, max_ulim=0, min_llim=99999999;

    DeclareNew( struct _summary_state *, state ) ;

    /* Allocate new arrays for integration limits */
    state->llim = NewArray(int, frameheight, "image.c for lower limit array" );
    state->ulim = NewArray(int, frameheight, "image.c for upper limit array" );

    for (i=0 ; i<frameheight ; i++) {

	/* Convert strings to sample points, allowing for cycles units */
	state->llim[i] = Cycles( llimstr, cfs[i], Samplerate() );
	state->ulim[i] = Cycles( ulimstr, cfs[i], Samplerate() );

	/* Check that llim < ulim, and quit if not so */
	if (state->llim[i] >= state->ulim[i])
	    stitch_error("Warning: gensas integration limits badly ordered\n");

	/* Find limits on required framewidth */
	if (state->ulim[i] > max_ulim)
	    max_ulim = state->ulim[i];
	if (state->llim[i] < min_llim)
	    min_llim = state->llim[i];
    }

    if (min_llim > 0) min_llim = 0;
    framewidth = max_ulim - min_llim;

    /* Adjust the integration limits for an sai with triggering on the      */
    /* right. (The parameters llim <= ulim, but the resulting sai indices   */
    /* state->llim[i] > state->ulim[i]).                                    */
    /* This is done by subtracting from (framewidth-1), where framewidth is */
    /* the maximum required number of points, (ie ulim in points).          */

    for (i=0 ; i<frameheight ; i++) {
	state->llim[i] = framewidth - state->llim[i] ;
	state->ulim[i] = framewidth - state->ulim[i] ;
    }

    framewidth++;   /* extra point as array starts from zeroth location */

    /* a blockings requests up into requests of the size specifed */
    state->source      = NewBlockingSource( source, sizeof ( DataType ) * framewidth * frameheight ) ;
    state->framewidth  = framewidth ;
    state->frameheight = frameheight ;
    state->scale       = scale ;

    return ( SetFillableSource( state, summary_callback, "image.c summarising sai" ) ) ;
}

/*********************** Routines supporting Summary ***********************/

static Pointer summary_callback( state, bytes, buffer )
struct _summary_state *state ;
ByteCount *bytes ;
DataType *buffer ;
{
    register int last = *bytes == 0 ;
    register int i, j, ulim, llim, point, points=ToPoints(DataType,*bytes) ;
    register DataType *sairow;
#ifdef FLOAT
    register DataType sum ;
#else
    register long sum ;
#endif

    /* Pull an sai frame */

    for( point=0 ; point < points ; )
	/* For each channel (row) in the sai, sum the row between the given limits */
	for(i=0 ; i<state->frameheight ; i++) {

	    sairow = PullItems(state->source,state->framewidth,DataType);

	    ulim = state->ulim[i];
	    llim = state->llim[i];
	    sum = 0;
	    for (j=ulim ; j<=llim ; j++)
		sum += sairow[j];

	    /* store the row-sum, scaled and normalized for the range of the sum */
	    buffer[point++] = sum*state->scale / (llim-ulim);
	}

    if( !last )
	return ( (Pointer) buffer ) ;
    else {
	Delete( state->llim ) ;
	Delete( state->ulim ) ;

	return ( DeleteFillableSource( state ) ) ;
    }
}

/********************************************************************************/
/*                                                                              */
/* initlist initialises the linked list pointer before the other list functions */
/* call the list.                                                               */
/*                                                                              */
/********************************************************************************/
void initlist()
{ 
  static int status=0;
  if (status==0)
    start = endl = pp = qq = (node *)malloc(sizeof(node));
  status=1;
}

/********************************************************************************/
/*                                                                              */
/* Function getnode allocates storage for a link list node and returns a        */
/* pointer to that node                                                         */
/*                                                                              */
/********************************************************************************/
node *getnode()
{
  node *q;
  q=(node *)malloc(sizeof(node));
  if (q==NULL) exit(66);
  return q;
}

/********************************************************************************/
/*                                                                              */
/* The function insertl inserts the value x at the end of the linked list and   */
/* moves the end point by another node                                          */
/*                                                                              */
/********************************************************************************/
node *insertl(x)
int x;
{
  node *q=endl; 
  endl=getnode();
  q->next=endl;
  q->index=x;
  q->val=0;
  /* fprintf(testfile, "insertl: index is %d, val is %d\n", q->index, q->val); */
  return q;
}

/********************************************************************************/
/*                                                                              */
/* The function inss searches the list from the start, stops 3 nodes before the */
/* value x, assigns the value -2000 to the three nodes till the (x+1)th node.   */
/*                                                                              */
/********************************************************************************/
void inss(x)
int x;
{
  node *q=start;
  int   y=x-3;
  int   z=x;
  while (q!=endl)
    {
      if (q->index == y)
	do
	  {
	    q->val=-2000;
	    q=q->next;
            /* fprintf(stderr, "Inss(x) index=%d\n",q->index); */
	  } while (q->index < z);
      q=q->next;
    }
}


/********************************************************************************/
/*                                                                              */
/* output_simple_strobe_info routine creates a file pointed by trigger file     */
/* which writes the points in the Nap pulse which initiates strobing            */
/*                                                                              */
/********************************************************************************/
void output_simple_strobe_info(state, value)
struct _sai_state *state;
short value;
{

if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on"))) 
     fwrite(&(value), sizeof(short), 1 , trigger_file);
}

/********************************************************************************/
/*                                                                              */
/* doList_strobe_info creates a file (pointed by "trigger_file") which writes   */
/* in binary shorts the points in the Nap pulse which initiates a strobe        */
/* It uses a linked list data structure to keep track of the actual time which  */
/* causes a strobe, rather than the time when the strobing event occurs.        */
/* At stcrit 4 or 5, strobing occurs nwid ms after a nap pulse peak has been    */
/* located.                                                                     */ 
/*                                                                              */
/********************************************************************************/
void doList_strobe_info(state, chan)
struct _sai_state *state;
int chan;
{
  int storev;
	    
  if (state->Stcrit==2)
    storev=(state->time)-2;
  else if (state->Stcrit==5)
    storev=state->def_strobe_candidate_time[chan];
 /*   fprintf(stderr, "Storev value of doList is = %d\n", storev);} */
  else
    storev=state->trigtime[chan];

  inss(storev);
 
  while (pp->index<=storev)  /* (pp!=qq) */
    {
      if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on"))) 
	fwrite(&(pp->val), sizeof(short), 1, trigger_file);
        /* fprintf(testfile, "doList: index is %d, val is %d\n", pp->index, pp->val); */
      pp=pp->next;
    }
  /* pp=qq; */
}

/********************************************************************************/
/*                                                                              */
/* output_thresh_info writes threshold values (in binary shorts or floats) to   */
/* the file opened by the file pointer "thresho_file"                           */
/*                                                                              */
/********************************************************************************/
void output_thresh_info(state, chan)
struct _sai_state *state;
int chan;
{
  S[chan]=state->thresh[chan];

if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on")))
#ifdef FLOAT
	fwrite((S+chan), sizeof(float), 1, thresho_file); 
#else
	fwrite((S+chan), sizeof(short), 1, thresho_file);
#endif 

}

/********************************************************************************/
/*                                                                              */
/* decay_strobe_threshold decays the strobe threshold at every clock            */
/* tick either linearly or exponentially                                        */
/*                                                                              */
/********************************************************************************/
void decay_strobe_threshold(state, chan)
struct _sai_state *state;
int chan;
{

  if (!strcmp(exposwitch, "on"))
    state->thresh[chan] *= state->triggerdecay;
  else
    {
      state->thresh[chan]  -= state->tlim[chan];
      if (state->thresh[chan] < 0) 
	state->thresh[chan]=0;
    }
}

/********************************************************************************/
/*                                                                              */
/* This routine prints debugging information of the local max search            */
/* with timeout (i.e. stcrit_ai=5)                                              */
/*                                                                              */
/********************************************************************************/
void print_trigger_debugging_info(state, chan, strobe_lag, stage)
struct _sai_state *state;
int chan;
int strobe_lag;
int stage;
{
  if (stage==1)
    {
      if (strobe_lag > state->def_strobe_candidate_time[chan])
	val[chan]=strobe_lag-state->def_strobe_candidate_time[chan];
      if (strobe_lag == state->def_strobe_candidate_time[chan])
	val[chan]=0; 
      
      if (!strcmp(state->switch_info,"off"));
      else
	fprintf(infotxt_file, "Chan %d: :: Actual time of trigger: s_lag=%d, trigtime=%d\n", chan, strobe_lag, state->def_strobe_candidate_time[chan]);
    }

/*   fprintf(infotxt_file, "Chan %d: Actual time of trigger: s_lag=%d, trigtime=%d, shift=%d\n", chan, strobe_lag, state->def_strobe_candidate_time[chan], val[chan]*state->chans); */
       
  
  if (stage==2)
    {
      if (!strcmp(state->switch_info,"off"));
      else
	fprintf(infotxt_file, "Chan %d: Loc Max Start Time = %d, trigcurrent %d, Threshold=%f, shift=%d\n", chan, state->time, state->trigheight[chan], state->thresh[chan], (ttime[chan]-initial_strobe_candidate_time[chan]));
    }

  if (stage==3)
    {
      if (!strcmp(state->switch_info,"off"));
      else
	fprintf(infotxt_file, "Chan %d:    Loc Max is %d Time is %d\n", chan, state->trigheight[chan], state->trigtime[chan] );
    }

  if (stage==4)
    {
      if (!strcmp(state->switch_info,"off"));
      else
	fprintf(infotxt_file, "Chan %d: At the End of Search (At peak), trigtime = %d time = %d height = %d\n", chan, state->trigtime[chan], state->time, state->trigheight[chan]);
    }
  if (stage==5)
    {
      if (!strcmp(state->switch_info,"off"));
      else
	fprintf(infotxt_file, "Chan %d: At the End of Search, trigtime = %d time = %d height = %d\n", chan, state->trigtime[chan], state->time, state->trigheight[chan]);  
    }
  
  
}