Mercurial > hg > aim92
diff model/image.c @ 0:5242703e91d3 tip
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
| author | tomwalters |
|---|---|
| date | Fri, 20 May 2011 15:19:45 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/model/image.c Fri May 20 15:19:45 2011 +0100 @@ -0,0 +1,1271 @@ +/* + 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]); + } + + +} +