Mercurial > hg > aim92

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
author tomwalters
date Fri, 20 May 2011 15:19:45 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5242703e91d3
1 /*
2 SCCS VERSION 1.29 WITH 1.28 STRIPPED OUT.
3
4
5
6 Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
7 ===========================================================================
8
9 Permission to use, copy, modify, and distribute this software without fee
10 is hereby granted for research purposes, provided that this copyright
11 notice appears in all copies and in all supporting documentation, and that
12 the software is not redistributed for any fee (except for a nominal shipping
13 charge). Anyone wanting to incorporate all or part of this software in a
14 commercial product must obtain a license from the Medical Research Council.
15
16 The MRC makes no representations about the suitability of this
17 software for any purpose. It is provided "as is" without express or implied
18 warranty.
19 */
20
21 /***************************************************************************
22 * image.c Version with trigger at right end of sai.
23 * =======
24 * Input data derived from a cochleagram is used to construct a stabilised
25 * auditory image (sai) of size `chans' * `framewidth' points.
26 *
27 * Organisation of data.
28 * Data is 2-dimensional: in time (rows) and frequency (columns).
29 *
30 * Input data is ordered by columns. The input array consists of blocks (columns)
31 * of `chans' points, each column being one time sample of `chans' channels.
32 * The zero'th point in each column is the lowest frequency channel.
33 *
34 * Output (sai) data is ordered by rows. The output array consists of blocks
35 * (rows) of `framewidth' points.
36 * The zero'th point in each row is the "oldest" time, as the time origin is
37 * on the right edge of the sai.
38 *
39 * Each sai is called a `frame', and its dimensions are `chans' rows by
40 * `framewidth' columns. In the code, the `framewidth' is often represented
41 * as imagewidth/chans, where imagewidth is the total amount of data in an
42 * sai frame, (ie. imagewidth=chans*framewidth).
43
44 * This module has been modified to have 5 levels of strobing.
45 A Jay Datta 6/03/95
46 ***************************************************************************/
47
48 #include <math.h>
49 #include <string.h>
50
51 #include "options.h"
52 #include "stitch.h"
53 #include "source.h"
54 #include "model.h"
55 #include "units.h"
56 #include "image.h"
57 #include "calc.h"
58
59 #define MAX_BUFFER (1l<<15)
60
61 #ifdef FLOAT
62 #define INPUT_SCALE 0.25
63 #define DECAY_SCALE 0.875
64 #else
65 #define INPUT_SHIFT 2
66 #define DECAY_ROUND 7
67 #endif
68 #define DECAY_SHIFT 3
69
70 #define MAXCHAN 1000
71
72 #ifndef lint
73 static char *sccs_id = "@(#)image.c 1.27 J. Holdsworth, M. Allerhand (MRC-APU) 6/6/91" ;
74 #endif
75
76 #ifdef DSP32
77 static int segsize = 16 ;
78 #else
79 static int segsize = 50 ;
80 #endif
81
82 typedef enum { False, True } Bool ;
83 typedef struct Node{int index; short val; struct Node *next;} node;
84
85 static void doColumn();
86 static Pointer sai_callback();
87 static Pointer summary_callback();
88 static void save_callback();
89 static void addIn();
90 static void decayImage();
91
92 void output_simple_strobe_info();
93 void doList_strobe_info();
94 void output_thresh_info();
95 void decay_strobe_threshold();
96 void print_trigger_debugging_info();
97
98 void initlist();
99 node *getnode();
100 node *insertl();
101 void inss();
102
103 FILE *trigger_file;
104 FILE *infotxt_file;
105 FILE *thresho_file;
106
107 int locmax_searchtime[MAXCHAN];
108 int locmax_searchstart[MAXCHAN];
109 int val[MAXCHAN];
110 int ltime[MAXCHAN];
111 int initial_strobe_candidate_time[MAXCHAN];
112 int ttime[MAXCHAN];
113
114 static node *start, *endl, *pp, *qq;
115 static char *exposwitch;
116
117 #ifdef FLOAT
118 float *S;
119 #else
120 short *S;
121 #endif
122
123
124 /***************************************************************************
125 * struct _sai_state:
126 * The source->info->state structure which is initialised by the
127 * routine Sai, and is then the state argument for the callback
128 * function: sai_callback.
129 ***************************************************************************/
130 struct _sai_state {
131 struct _pullable_source parent ;
132 Source source;
133 void (*trigger)(); /* trigger algorithm (originally doColumn) */
134 unsigned chans; /* number of filter channels */
135 unsigned framestep; /* sai display update period */
136 DataType *image; /* sai array (DataType = short) */
137 unsigned imagewidth; /* length of sai array (framewidth*chans) */
138 unsigned imagenwidth; /* portion of sai array for transient info */
139 unsigned decay_time; /* a factor of sai decay time constant */
140 unsigned decay_count;
141 ScalarType *framedecay; /* array of factors for cochleagram frame */
142 int *cps; /* array of 1/centre-freq (in samples) */
143 int *isPulse; /* array of bools, to see if data > thres */
144
145 int Stcrit; /* int restriction */
146 int SusThresh; /* Level 2 user settable thresh */
147 ScalarType triggerdecay; /* prop/point decay for strobe threshold */
148 int *trigtime; /* time of last trigger pulse */
149 int *trigheight; /* height of last trigger pulse */
150 int *isStrobe; /* array of booleans, set when pulse peak found */
151 float *thresh; /* array of thresholds, one per channel */
152 float *tlim; /* linear decay value */
153 int time; /* keep track of sai */
154 int *previnput; /* lock till new pulse found */
155 int *def_strobe_candidate_time;
156 int *def_strobe_height;
157 char *switch_info;
158 int tlim1;
159 int stlag; /* strobe nwid first attempt to decouple nwids */
160 };
161
162 /***************************************************************************
163 * Sai:
164 * Sai is called from SaiEntry (in module model.c).
165 * Routines in this module it uses are:
166 * doColumn, (which uses: addIn)
167 * sai_callback, (which uses: decayImage)
168 *
169 * SaiEntry (in model.c) is the entry-point function to the stabilised image
170 * processing module. Its purpose is to prepare the arguments for Sai by
171 * converting parameter-option strings to integers, and then to call Sai.
172 *
173 * Sai is called with the following arguments (in model.c):
174 * Sai( source, Frameheight(), segment, width,
175 * (int) Samples( decaystr )/Framestep(), (int) Freq( samplestr ),
176 * frequencies, Samples(ttstr)/Framestep(), Samples(cgmstr)/Framestep(),
177 * (int) Freq("250Hz"), (int) Freq("20Hz") ) ;
178 *
179 * The arguments are as follows:
180 * (The option name refers to the description given by "help". The values etc
181 * can be found in table.c as described in model.docs).
182 *
183 * Argument Option name Default Comment
184 * -------- ----------- ------- -------
185 * chans (see below) Number of filter channels
186 * framestep segment_sai 16ms SAI display update period
187 * framewidth duration_sai 32ms SAI duration (ms)
188 * decay decayrate_sai 32ms SAI decay time constant (SU/ms) 15ms NOW Default
189 * samplerate samplerate 20000. Input wave sample rate (Hz)
190 * cfs (see below) An array of filter centre frequencies
191 * cgmdecay napdecay_sai 16ms NAP (ie cochleagram) decay time constant
192 *
193 * Both "chans" and "cfs" are setup from the routine "updateFrequencies" in
194 * model.c, because both arguments are derived from three basic parameters:
195 *
196 * minstr mincf_afb 220Hz Minimum center frequency (Hz)
197 * maxstr maxcf_afb 4400Hz Maximum center frequency (Hz)
198 * denstr dencf_afb 4. Filter density (filters/critical band)
199 *
200 * The "chans" argument is calculated in the routine "NumberCenterFrequencies",
201 * and the result is copied into the "frameheightstr" by a call to the routine
202 * "setFrameheight".
203 * The "cfs" argument is calculated in the routine "GenerateCenterFrequencies".
204 * (Both CenterFrequency routines are in gamma_tone.c).
205 *
206 ***************************************************************************/
207
208 Source Sai(source, chans, framestep, pwidth, nwidth, decay, samplerate, cfs, ttdecay, cgmdecay, tlim1, tlim2, suslevel, susthresh, switchinfo, expoSwitch,stlag)
209 Source source ;
210 int chans, framestep, pwidth, nwidth, decay ;
211 int samplerate;
212 double *cfs ;
213 double ttdecay, cgmdecay ;
214 int tlim1, tlim2;
215 int suslevel;
216 char *susthresh;
217 char *switchinfo;
218 char *expoSwitch;
219 int stlag;
220 {
221
222 DeclareNew( struct _sai_state *, state ) ;
223 int i, chan, startup, segment, framewidth ;
224 double pts_per_ms, lin_napdec, lin_napfac ; /* roy 11-8-92 roy */
225 char *ThreshoFile, *InfotxtFile, *TriggerFile;
226
227 exposwitch=(char *)calloc(5, sizeof(char));
228 exposwitch=expoSwitch;
229
230 ThreshoFile = (char *)calloc(256, sizeof(char));
231 InfotxtFile = (char *)calloc(256, sizeof(char));
232 TriggerFile = (char *)calloc(256, sizeof(char));
233
234 state->switch_info = (char *)calloc(256,sizeof(char));
235 state->switch_info = switchinfo;
236
237 if (!strcmp(susthresh, "on"))
238 state->SusThresh=0;
239 else if (!strcmp(susthresh, "off"))
240 state->SusThresh=1;
241 else
242 state->SusThresh=atoi(susthresh);
243
244 if (!strcmp(state->switch_info, "off"))
245 ;
246 else if (!strcmp(state->switch_info, "on"))
247 infotxt_file=stdout;
248 else{
249 strcpy(InfotxtFile, switchinfo);
250 strcat(InfotxtFile, ".info");
251 infotxt_file=fopen(InfotxtFile, "w");
252
253 strcpy(TriggerFile, switchinfo);
254 strcat(TriggerFile, ".trigger");
255 trigger_file=fopen(TriggerFile, "w");
256
257 strcpy(ThreshoFile, switchinfo);
258 strcat(ThreshoFile, ".thresh");
259 thresho_file=fopen(ThreshoFile, "w");
260 }
261
262 initlist();
263 state->tlim1=0;
264
265 framewidth = pwidth+nwidth;
266
267 pts_per_ms = (double)samplerate/1000 ;
268
269 state->stlag=stlag; /* NEW line 95April */
270 /* Initialise members of the sai state structure */
271 state->source = NewRetainingSource( source, sizeof ( DataType ) * framewidth * chans ) ;
272 state->chans = chans ;
273 state->trigger = doColumn ;
274 state->framestep = framestep * chans ;
275 state->imagewidth = framewidth * chans ; /* size of sai buffer */
276 state->imagenwidth = nwidth * chans ; /* size of buffer for transient info */
277 state->decay_time = decay * -log( 1. - 1. / ( 1 << DECAY_SHIFT ) ) ;
278 state->decay_count = state->decay_time ;
279
280 state->Stcrit = suslevel;
281 if (stlag==0 && state->Stcrit==5) /* was nwid NEW */
282 state->Stcrit = suslevel-1;
283
284 /* state->SusThresh = susthresh; */
285
286 if (!strcmp(exposwitch, "on")) fprintf(stderr, "Exponential trigger decay SWITCH is ON\n");
287 /* if (!strcmp(exposwitch, "off")) printf("SWITCH is OFF\n"); */
288
289 if (!strcmp(exposwitch, "on"))
290 state->triggerdecay = (double) 1-((ttdecay / pts_per_ms ) /100.0) ;
291 else
292 state->triggerdecay = (double) (ttdecay /pts_per_ms) / 100.0 ;
293
294 if (!strcmp(switchinfo, "off"))
295 ;
296 else
297 fprintf(infotxt_file, "\ntriggerdecay = %15e\n", state->triggerdecay);
298
299 state->isStrobe = NewZeroedArray( int, state->chans, "image.c for Strobe" ) ;
300 state->isPulse = NewZeroedArray( int, state->chans, "image.c for trigger pulse" ) ;
301 state->trigtime = NewZeroedArray( int, state->chans, "image.c for trigger time" ) ;
302 state->trigheight = NewZeroedArray( int, state->chans, "image.c for trigger height" ) ;
303 state->thresh = NewZeroedArray( float, state->chans, "image.c for thresh" ) ;
304 state->tlim = NewZeroedArray( float, state->chans, "image.c for decay" ) ;
305 state->time = 0; /* Time initialized to Zero */
306 state->previnput = NewZeroedArray( int, state->chans, "image.c for pulse lock" ) ;
307
308 state->def_strobe_candidate_time = NewZeroedArray( int, state->chans, "image.c for prevtrig time" );
309 state->def_strobe_height = NewZeroedArray( int, state->chans, "image.c for prev trig" );
310
311 for (chan=0; chan < MAXCHAN; chan++)
312 locmax_searchtime[chan]=0;
313 for (chan=0; chan < MAXCHAN; chan++)
314 locmax_searchstart[chan]=0;
315 for (chan=0; chan < MAXCHAN; chan++)
316 val[chan]=0;
317 for (chan=0; chan < MAXCHAN; chan++)
318 initial_strobe_candidate_time[chan]=0;
319 for (chan=0; chan < MAXCHAN; chan++)
320 ltime[chan]=0;
321 for (chan=0; chan < MAXCHAN; chan++)
322 ttime[chan]=0;
323
324 #ifdef FLOAT
325 S=(float *)calloc(state->chans, sizeof(float));
326 #else
327 S=(short *)calloc(state->chans, sizeof(short));
328 #endif
329
330 /* Declare new arrays for members of state structure */
331
332 state->framedecay = NewArray( ScalarType, framewidth, "image.c for decay" ) ;
333 state->image = NewZeroedArray( DataType, state->imagewidth, "image.c for image" ) ;
334 state->cps = NewZeroedArray( int, state->chans, "image.c for channel centre periods" ) ;
335
336 /* Initialise cochleagram frame decay factors */
337 if (cgmdecay > 0)
338 {
339 /* 11111 roy 11-8-92 These mods change to a linear nap decay. 1111111 roy */
340 /* It requires pts_per_ms, lin_napdec and lin_napfac declared just after DeclareNew above */
341
342 lin_napdec = ((cgmdecay/100)/pts_per_ms)/pts_per_ms ; /* roy 19-8 roy */
343
344 /* cgmdecay is div by pts_per_ms to get back to napdecay in the original units */
345 /* When napdecay is interpreted as the % to decay per ms, and since we are */
346 /* operating in points, we need napdec/pts_per_ms, as the dec for each pt within the ms */
347 /* The decay vector is set to scale the NAP by 1.0 at strobe point. So, in the sai, */
348 /* it falls to the left AND RISES to the right of 0 ms. **** roy **** */
349
350 lin_napfac = 1.0 + nwidth*lin_napdec ;
351
352 for (i=0 ; i < framewidth ; i++)
353 {
354 if ( lin_napfac > 0 )
355 state->framedecay[i] = SCALAR( lin_napfac );
356 else
357 state->framedecay[i] = SCALAR( 0 );
358 lin_napfac -= lin_napdec ;
359 }
360 }
361 else /* disable decay factor on zero argument */
362 for (i=0 ; i < framewidth ; i++)
363 state->framedecay[i] = SCALAR( 1.0 ) ;
364
365 /* Initialise centre-period for each channel */
366 for (i=0 ; i < state->chans ; i++)
367 state->cps[i] = samplerate / cfs[i];
368
369 #if defined( PC ) || defined( THINK_C )
370 if( (long) ( MAX_BUFFER - (long) state->imagewidth * sizeof ( DataType ) ) < (long) ( state->chans * sizeof ( DataType ) ) )
371 stitch_error( "Sorry, image larger than maximum buffer size\n" ) ;
372 #endif
373
374 /* for (i=0; i<state->imagewidth;i++)
375 fprintf(stderr, "Image Buffer is %f, count = %d\n", *(state->image+i), (i+1)); */
376
377
378
379 return ( SetPullableSource( state, sai_callback, "image.c stabilised image" ) ) ;
380 }
381
382 /****************************************************************************
383 * sai_callback
384 * The callback function at source->info->callback in the source
385 * returned by Sai.
386 * A "frame" is a window over the cochleagram, of size chans * points.
387 * A "segment" is a stretch of the cochleagram which is segsize (ie 50) time
388 * samples long, ie segsize columns of the cochleagram.
389 *
390 * This callback routine is executed once for each sai frame which is plotted.
391 * "*bytes" is the number of bytes in the sai, ie:
392 * *bytes = state->imagewidth * sizeof(DataType)
393 ****************************************************************************/
394 static Pointer sai_callback( state, bytes )
395 struct _sai_state *state ;
396 ByteCount *bytes ;
397 {
398 register long col, point ;
399 register long points = state->framestep; /* total num points under frame */
400 #if defined( PC )
401 int segment = ( MAX_BUFFER - state->imagewidth * sizeof ( DataType ) ) / state->chans / sizeof ( DataType ) * state->chans ;
402 #else
403 int segment = segsize * state->chans ;
404 #endif
405 register DataType *input ; /* cochleagram (input) data array */
406 static int first=1;
407 static int test=0;
408 int chan;
409 int i;
410 /* #if defined( PC )
411 test=1;
412 #else
413 if (state->SusThresh == 0)
414 segment= 30 * state->chans ;
415 #endif */
416
417 /* test++; */
418
419
420 /* If size argument is zero, then by convention free space and return null */
421 if( *bytes == 0 ) {
422 Pull( state->source, 0 ) ;
423 Delete( state->image ) ;
424 return ( DeletePullableSource( state ) ) ;
425 }
426
427 /* Initially pull enough data to buffer transient info ahead of the trigger */
428 /* Note, the pull operation is not valid when the size arg is zero. */
429 /* The result would be a segmentation fault on the next pull, below. */
430
431 if (state->Stcrit==0 || state->Stcrit==1 || state->Stcrit==2 || state->Stcrit==3)
432 if (first) {
433 if (state->imagenwidth > 0)
434 input = PullItems(state->source, state->imagenwidth, DataType);
435 first = 0;
436 }
437
438
439 /* Search frame in segment blocks for trigger-points */
440 for( point=0 ; point < points ; point += segment ) {
441 if( segment > points - point ) /* Finish at end of frame */
442 segment = points - point ;
443 /* Get one segment of the input data */
444 /* Pull 50 columns of data from the source, in buffer *input. */
445 /* Keep one sai of data, from last pull, in buffer behind *input */
446 input = PullItems(state->source, segment, DataType);
447
448 for( col=0 ; col < segment ; col += state->chans ) {
449 /* Retard the "current" data point by transient buffer width */
450 /* (This retards the trigger by the transient time, about 5ms) */
451 if (state->Stcrit==0 || state->Stcrit==1 || state->Stcrit==2 || state->Stcrit==3)
452 state->trigger(input-state->imagenwidth, state) ;
453 else{
454 state->trigger( input, state) ; /* -state->imagenwidth, state ) ; */
455 /*fprintf(stderr, " %d ", points) ; */ }
456 input += state->chans ; /* next column */
457 /* decay image periodically, once every `decay_time' columns */
458 /* Lines moved down */
459 if( --state->decay_count <= 0 ) {
460 decayImage( state ) ;
461 state->decay_count = state->decay_time ;
462 }
463 } /* for each col */
464 } /* for each seg */
465
466 /* whatever the number of bytes requested, the image is always this size */
467 *bytes = state->imagewidth * sizeof ( *state->image ) ;
468
469 /* if (state->Stcrit==1 && state->SusThresh==0){
470 if (test%5==0)
471 decayImage( state );
472 test++;}
473 else if (state->Stcrit==1 && state->SusThresh>0){ changed 2 to 3
474 if (test%10==0)
475 decayImage( state );
476 test++;}
477 else{ */
478
479 return ( (Pointer) state->image ) ;
480 }
481
482 /****************************************************************************
483 * decayImage
484 * Periodically attenuate image to give computationally efficient image decay.
485 * Attenuate all points in the sai by a constant factor.
486 ****************************************************************************/
487 static void decayImage( state )
488 struct _sai_state *state ;
489 {
490 register DataType *image_ptr, *image_end ;
491 #ifdef FLOAT
492 register FLOAT image_decay = DECAY_SCALE ;
493 #endif
494
495 image_end = state->image + state->imagewidth ; /* end of sai array */
496 #ifdef FLOAT
497 for( image_ptr=state->image ; image_ptr < image_end ; )
498 *image_ptr++ *= image_decay ;
499 #else
500 for( image_ptr=state->image ; image_ptr < image_end ; image_ptr++ )
501 *image_ptr -= *image_ptr + DECAY_ROUND >> DECAY_SHIFT ;
502 #endif
503 return ;
504 }
505
506 /****************************************************************************
507 * doColumn
508 * Trigger algorithm for Sai.
509 * Installed as source->info->state->trigger.
510 *
511 ****************************************************************************/
512 static void doColumn( input, state )
513 DataType *input ;
514 struct _sai_state *state ;
515 {
516
517 register int chan ;
518 float strobelag;
519 int s_lag;
520
521
522 short *STI;
523 static int check=0;
524 int storev;
525 char buf[2];
526 short stipts=-2000;
527 short zeroval=0;
528
529 STI=(short *)calloc(state->chans, sizeof(short));
530
531 strobelag = (double) (state->time)-(state->stlag);
532
533 s_lag=(state->time)-(state->stlag);
534
535 /* For each channel up column, add a row of data to sai if the */
536 /* current point is non-zero. Add the row from the current point */
537 /* plus state->imagenwidth, to display the transient (nwidth) info. */
538
539
540 /* 000000000000000000000000 stcrit 0000000000000000000000000000000000 */
541 /* Add on every point : Low Pass Filter */
542
543 if (state->Stcrit==0)
544 {
545 for( chan=0 ; chan < state->chans ; chan++ )
546 {
547 addIn(input+(state->imagenwidth), chan, state);
548 output_simple_strobe_info(state, stipts);
549 }
550 } /* if (state->Stcrit==0) */
551
552
553 /* 111111111111111111111111 stcrit 1111111111111111111111111111111111 */
554 /* Add on every point greater than "stthresh_ai" (default 0) */
555
556 if (state->Stcrit==1)
557 {
558 for( chan=0 ; chan < state->chans ; chan++ )
559 {
560 if (input[chan] > state->SusThresh )
561 {
562 addIn(input+(state->imagenwidth), chan, state);
563 output_simple_strobe_info(state, stipts);
564 }
565 else if (input[chan] <=state->SusThresh)
566 output_simple_strobe_info(state, zeroval);
567 }
568 } /* if (state->Stcrit==1) */
569
570
571 /* 222222222222222222222222 stcrit 2222222222222222222222222222222222 */
572 /* Add on every Nap pulse peak */
573
574 if (state->Stcrit==2)
575 {
576 for( chan=0 ; chan < state->chans ; chan++ )
577 {
578 if (state->isStrobe[chan]) /* Check to see if it is time to addIn */
579 {
580 addIn(input+((state->imagenwidth)-2*(state->chans)), chan, state);
581 state->isStrobe[chan]=0; /* Reset strobe flag */
582 doList_strobe_info(state, chan);
583 }
584
585 /* Find a pulse peak routine */
586
587 if (!state->isPulse[chan]) /* Not in pulse, looking for pulse */
588 {
589 state->thresh[chan]=input[chan];
590 if (input[chan] > 0 /*state->thresh[chan] */ ) /* Start of a Nap pulse */
591 state->isPulse[chan] =1; /* set in pulse flag */
592 }
593 else /* In pulse, looking for peak */
594 {
595 if (input[chan] < state->thresh[chan]) /* peak found */
596 {
597 state->isPulse[chan]=0;
598 state->isStrobe[chan]=1; /* set strobe pending flag */
599 }
600 else
601 state->thresh[chan]=input[chan];
602 }
603 output_thresh_info(state, chan);
604 }
605 } /* if (state->stcrit==2) */
606
607
608
609 /* 333333333333333333333333 stcrit 3333333333333333333333333333333333 */
610 /* Add on pulse peaks that exceed strobe threshold (temporal shadow ) */
611 /* Primitive Local Max algorithm */
612
613 if (state->Stcrit==3)
614 {
615 for( chan=0 ; chan < state->chans ; chan++ )
616 {
617 if (state->isStrobe[chan]) /* Check to see if it is time to addIn */
618 {
619 addIn(input+((state->imagenwidth)-2*(state->chans)), chan, state);
620 state->isStrobe[chan] = 0; /* Reset strobe flag */
621 doList_strobe_info(state, chan);
622 }
623
624 /* Find a pulse peak routine */
625
626 if (!state->isPulse[chan]) /* Not in Pulse, looking for Pulse */
627 {
628 if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan])
629 { /* Previnput check ensures threshold decay
630 does not cut into pulse */
631
632 state->isPulse[chan] = 1; /* set flag: in pulse */
633 state->thresh[chan] = input[chan]; /* start search for peak */
634 }
635 else {
636 decay_strobe_threshold(state, chan);
637 if (input[chan]==0) /* in the zeroes between nap pulses */
638 state->previnput[chan] = 0; /* Set it to Zero previnput */
639 }
640 }
641
642 if (state->isPulse[chan]) /* Else in ongoing pulse looking for peak */
643 {
644 if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan])
645 /* pulse peak found */
646 {
647 state->isPulse[chan] = 0; /* clear flag: not in pulse */
648 state->isStrobe[chan] = 1; /* set Strobe pending flag */
649 if (!strcmp(exposwitch, "off"))
650 state->tlim[chan]=state->triggerdecay * state->thresh[chan];
651 state->trigtime[chan]=(state->time)-1;
652 /* state->previnput[chan]=input[cha */
653 }
654 else /* still in pulse, so continue search for peak */
655 state->thresh[chan] = input[chan];
656
657 } /* if in Pulse */
658 state->previnput[chan]=input[chan];
659 output_thresh_info(state, chan);
660 } /* for all chans... */
661 } /* if (state->Stcrit==3) */
662
663
664
665 /* 444444444444444444444444 stcrit 4444444444444444444444444444444444 */
666 /* Add on pulse peaks that exceed strobe threshold and which are not succeeded by
667 a larger pulse in nwid ms. (Better Local Max mechanism) */
668
669 if (state->Stcrit==4)
670 {
671 for( chan=0 ; chan < state->chans ; chan++ )
672 {
673 if (state->isStrobe[chan]) /* Check to see if it is time to addIn */
674 {
675 if ( strobelag-1 >= state->trigtime[chan])
676 /* initiates strobe after nwid ms */
677 {
678 addIn(input-(state->chans), chan, state);
679 doList_strobe_info(state, chan);
680 state->isStrobe[chan] = 0; /* Reset strobe flag */
681 state->trigheight[chan] = 0; /* Reset local max value */
682
683 }
684 }
685
686 /* Find a pulse peak routine */
687
688 if (!state->isPulse[chan]) /* Not in Pulse, looking for pulse */
689 {
690 if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan])
691 { /* Previnput check ensures threshold
692 decay does not cut into pulse */
693
694 state->isPulse[chan] = 1; /* set flag: in pulse */
695 state->thresh[chan] = input[chan]; /* start search for peak */
696 /* state->previnput[chan]= 1; in pulse till the next zero value */
697 }
698 else
699 {
700 decay_strobe_threshold(state, chan);
701 if (input[chan]==0) /* in the zeroes between nap pulses */
702 state->previnput[chan] = 0; /* Free previnput */
703 }
704 }
705
706 if (state->isPulse[chan]) /* Else in ongoing pulse, looking for peak */
707 {
708 if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan])
709 /* pulse peak found */
710 {
711 state->isPulse[chan] = 0; /* clear flag: not in pulse */
712 state->isStrobe[chan] = 1; /* set Strobe pending flag */
713
714 if ( state->thresh[chan] > state->trigheight[chan] )
715 /* Check for LOCAL MAX */
716 {
717 state->trigheight[chan] = state->thresh[chan];
718 state->trigtime[chan] = (state->time)-1; /* Set strobe time */
719 if (!strcmp(exposwitch, "off"))
720 state->tlim[chan]=state->triggerdecay * state->trigheight[chan];
721 }
722 }
723 else /* still in pulse, so continue search for peak */
724 state->thresh[chan] = input[chan];
725 } /* if in Pulse */
726 state->previnput[chan]=input[chan];
727 output_thresh_info(state, chan);
728 } /* for all chans... */
729 } /* if (state->Stcrit==4) */
730
731
732
733 /* 555555555555555555555555 stcrit 5555555555555555555555555555555555 */
734 /* Add on pulse peaks that exceed storbe threshold and which are not succeeded by
735 a larger pulse in nwid ms, provided total lag is < 2*nwid ms
736 Local Max mechanism with a timeout */
737
738 if (state->Stcrit==5)
739 {
740 for (chan = 0; chan < state->chans; chan++)
741 {
742 if(state->isStrobe[chan]) /* Check to see if it is time to addIn */
743 { /* check for waiting till the first local max found */
744 if (strobelag-1 >= state->def_strobe_candidate_time[chan])
745 /* initiate strobe after nwid ms */
746 {
747 /* print_trigger_debugging_info(state, chan, s_lag, 1); */
748 addIn(input-state->chans, chan, state);
749 doList_strobe_info(state, chan);
750
751 state->def_strobe_candidate_time[chan] = 0; /* Reset strobe time */
752 state->def_strobe_height[chan] = 0; /* Reset local max value */
753 state->trigheight[chan]=0;
754 state->isStrobe[chan]=0; /* Reset strobe flag */
755
756 }
757 }
758
759 /* Find a pulse peak routine */
760
761 if (!state->isPulse[chan]) /* Not in pulse, looking for pulse */
762 {
763 if (input[chan] > state->thresh[chan] && input[chan] > state->previnput[chan])
764 { /* Previnput check ensures threshold
765 decay does not cut into pulse */
766
767 state->isPulse[chan] = 1; /* set flag: in pulse */
768 state->thresh[chan] = input[chan]; /* start search for peak */
769 /* state->previnput[chan] = 1; in pulse till the next zero value */
770 }
771 else
772 {
773 decay_strobe_threshold(state, chan);
774 if (input[chan]==0) /* in the zeroes between nap pulses */
775 state->previnput[chan]=0; /* Free previnput */
776 }
777
778 }
779
780 if (state->isPulse[chan]) /* Else in ongoing pulse, looking for peak */
781 {
782 if (input[chan] <= state->thresh[chan] && input[chan] <= state->previnput[chan])
783 /* Pulse Peak Found */
784 {
785 state->isPulse[chan]=0; /* clear flag: not in pulse */
786
787 if (locmax_searchstart[chan]==0 )
788 /* && state->time>=((ttime[chan]-initial_strobe_candidate_time[chan])+ltime[chan])) AJD 1-2-96 */
789 { /* start local max time (for timeout)
790 and shift search window */
791 locmax_searchtime[chan]=1; /* local max time started */
792 locmax_searchstart[chan]=1; /* local max flag set */
793 print_trigger_debugging_info(state, chan, s_lag, 2);
794 initial_strobe_candidate_time[chan]=(state->time)-1; /* local max start time noted */
795
796 }
797
798 /* if (locmax_searchstart[chan]==1) */ /* 8rd March, 1995 AJayDatta */
799 if (state->thresh[chan] > state->trigheight[chan])
800 /* Check for LOCAL MAX */
801 {
802 state->trigheight[chan]=state->thresh[chan];
803 state->trigtime[chan]=(state->time)-1; /* Set strobe time */
804 if (!strcmp(exposwitch, "off"))
805 state->tlim[chan]=state->triggerdecay * state->trigheight[chan];
806 if (locmax_searchstart[chan]==1)
807 print_trigger_debugging_info(state, chan, s_lag, 3);
808 } /* if greater than thresh */
809
810
811 /* Reached end of time-out period (While in the peak finding stage !)
812 End of the timeout period can occur at two stages of the algorithm,
813 at a Nap pulse peak or at any point of the pulse train. */
814
815 if (locmax_searchtime[chan]==(state->stlag)) /* imagenwidth/state->chans)) */
816 {
817 print_trigger_debugging_info(state, chan, s_lag, 4);
818
819 state->def_strobe_candidate_time[chan]=state->trigtime[chan];
820 state->def_strobe_height[chan]=state->trigheight[chan];
821 locmax_searchtime[chan]=0;
822 locmax_searchstart[chan]=0;
823 state->trigheight[chan]=0;
824 state->trigtime[chan]=0;
825 state->isStrobe[chan]=1; /* set Strobe pending flag */
826 ltime[chan]=state->time; /* End of search time */
827 ttime[chan]=state->def_strobe_candidate_time[chan]; /* Local Max time */
828 }
829 /* reset search period to zero */
830 /* as y ms search period ends, update the info for strobe and
831 clear current_ values for the next y ms */
832 } /* Pulse Peak Loop */
833 else
834 state->thresh[chan]=input[chan];
835 } /* if (state->isPulse..) */
836
837 if (locmax_searchtime[chan] == (state->stlag)) /* imagenwidth/state->chans)) */
838 /* && state->def_strobe_candidate_time[chan]==0) */
839 { /* start local max time (for timeout) */
840 print_trigger_debugging_info(state, chan, s_lag, 5);
841
842 state->def_strobe_candidate_time[chan]=state->trigtime[chan];
843 state->def_strobe_height[chan]=state->trigheight[chan];
844 locmax_searchtime[chan]=0;
845 locmax_searchstart[chan]=0;
846 state->trigheight[chan]=0;
847 state->trigtime[chan]=0;
848 state->isStrobe[chan]=1; /* set Strobe pending flag */
849 ltime[chan]=state->time;
850 ttime[chan]=state->def_strobe_candidate_time[chan];
851 }
852 /* reset search period to zero */
853 /* as y ms search period ends, update the info for strobe and
854 clear current_ values for the next y ms */
855
856
857 if (locmax_searchstart[chan]==1)
858 (locmax_searchtime[chan])++;
859 state->previnput[chan]=input[chan];
860 output_thresh_info(state, chan);
861
862 } /* for all chans */
863 } /* if (state->Stcrit== 5) */
864
865
866 qq=insertl(state->time); /* insert time to list */
867 (state->time)++;
868 return;
869
870 } /* doColumn */
871
872
873 /****************************************************************************
874 * addIn
875 * add row cochleagram into stabilised image
876 ****************************************************************************/
877 static void addIn( input, chan, state )
878 DataType *input ;
879 int chan ;
880 struct _sai_state *state ;
881 {
882 register DataType *channel_ptr, *image_ptr, *end;
883 #ifdef FLOAT
884 register FLOAT input_scale = INPUT_SCALE ;
885 #endif
886 register ScalarType *decayfactor = state->framedecay;
887
888 /* Initialize channel_ptr to the input data to be added in. */
889 /* `input' points to a particular column in the cochleagram, */
890 /* and `chan' indexes a particular channel in this column. */
891 channel_ptr = input + chan ;
892 end = channel_ptr - state->imagewidth ;
893
894 /* Initialize image_ptr to end of sai row corresponding to `chan'.*/
895 /* `state->image' points to the start of the sai, */
896 /* so increment pointer by chan*framewidth to get to row `chan'. */
897 /* To get to end of row, increment by ((chan+1)*framewidth)-1. */
898 /* But framewidth = state->imagewidth / state->chans, and hence: */
899 /* Added fix; as new spiral.c accepts nwid AJD 17-3-95 */
900 if (state->Stcrit==4 && state->stlag==0)
901 image_ptr = (state->image + state->imagewidth / state->chans * (chan+1));
902 /* 13-3-95 removed -1; ajd */
903 else
904 image_ptr = (state->image + state->imagewidth / state->chans * (chan+1))-1;
905
906 /* Decrement channel_ptr by columns, from the initial column, */
907 /* until one complete sai row has been added into, which is when */
908 /* channel_ptr has decreased by a whole imagewidth. */
909 while( channel_ptr > end ) {
910 #ifdef FLOAT
911 *image_ptr += ( ( *channel_ptr * input_scale ) * (*decayfactor++) ) ;
912 #else
913 *image_ptr += DESCALE( ( *channel_ptr >> INPUT_SHIFT ) * (*decayfactor++) ) ;
914 #endif
915 /* if (state->chans <= 40) */
916 if (*image_ptr > 32767.0)
917 *image_ptr = 32767.0; /* 32767 */
918
919 image_ptr--;
920
921 channel_ptr -= state->chans ; /* next (ie previous) input time point */
922 }
923
924 return ;
925 }
926
927
928 /***************************************************************************
929 * Summary is not called in this module.
930 * (see SummaryEntry in model.c, which is the entry point for gensas).
931 * Summary computes a row summary spectrogram and is called during the
932 * program "gensas".
933
934 * Routines in this module it uses are:
935 * summary_callback
936 ****************************************************************************/
937
938 struct _summary_state {
939 struct _fillable_source parent ;
940 Source source;
941 int framewidth; /* number of time-points or cols in an sai */
942 int frameheight; /* number of channels or rows in an sai */
943 double scale;
944 int *llim; /* integration limits (lower & upper) for sai summary data */
945 int *ulim;
946 };
947
948
949 Source Summary( source, frameheight, scale, cfs, samplerate, llimstr, ulimstr)
950 Source source ;
951 int frameheight ;
952 double scale ;
953 double *cfs, samplerate ; /* array of channel centre frequencies */
954 char *llimstr, *ulimstr; /* lower and upper limit strings */
955 {
956 int i, framewidth, max_ulim=0, min_llim=99999999;
957
958 DeclareNew( struct _summary_state *, state ) ;
959
960 /* Allocate new arrays for integration limits */
961 state->llim = NewArray(int, frameheight, "image.c for lower limit array" );
962 state->ulim = NewArray(int, frameheight, "image.c for upper limit array" );
963
964 for (i=0 ; i<frameheight ; i++) {
965
966 /* Convert strings to sample points, allowing for cycles units */
967 state->llim[i] = Cycles( llimstr, cfs[i], Samplerate() );
968 state->ulim[i] = Cycles( ulimstr, cfs[i], Samplerate() );
969
970 /* Check that llim < ulim, and quit if not so */
971 if (state->llim[i] >= state->ulim[i])
972 stitch_error("Warning: gensas integration limits badly ordered\n");
973
974 /* Find limits on required framewidth */
975 if (state->ulim[i] > max_ulim)
976 max_ulim = state->ulim[i];
977 if (state->llim[i] < min_llim)
978 min_llim = state->llim[i];
979 }
980
981 if (min_llim > 0) min_llim = 0;
982 framewidth = max_ulim - min_llim;
983
984 /* Adjust the integration limits for an sai with triggering on the */
985 /* right. (The parameters llim <= ulim, but the resulting sai indices */
986 /* state->llim[i] > state->ulim[i]). */
987 /* This is done by subtracting from (framewidth-1), where framewidth is */
988 /* the maximum required number of points, (ie ulim in points). */
989
990 for (i=0 ; i<frameheight ; i++) {
991 state->llim[i] = framewidth - state->llim[i] ;
992 state->ulim[i] = framewidth - state->ulim[i] ;
993 }
994
995 framewidth++; /* extra point as array starts from zeroth location */
996
997 /* a blockings requests up into requests of the size specifed */
998 state->source = NewBlockingSource( source, sizeof ( DataType ) * framewidth * frameheight ) ;
999 state->framewidth = framewidth ;
1000 state->frameheight = frameheight ;
1001 state->scale = scale ;
1002
1003 return ( SetFillableSource( state, summary_callback, "image.c summarising sai" ) ) ;
1004 }
1005
1006 /*********************** Routines supporting Summary ***********************/
1007
1008 static Pointer summary_callback( state, bytes, buffer )
1009 struct _summary_state *state ;
1010 ByteCount *bytes ;
1011 DataType *buffer ;
1012 {
1013 register int last = *bytes == 0 ;
1014 register int i, j, ulim, llim, point, points=ToPoints(DataType,*bytes) ;
1015 register DataType *sairow;
1016 #ifdef FLOAT
1017 register DataType sum ;
1018 #else
1019 register long sum ;
1020 #endif
1021
1022 /* Pull an sai frame */
1023
1024 for( point=0 ; point < points ; )
1025 /* For each channel (row) in the sai, sum the row between the given limits */
1026 for(i=0 ; i<state->frameheight ; i++) {
1027
1028 sairow = PullItems(state->source,state->framewidth,DataType);
1029
1030 ulim = state->ulim[i];
1031 llim = state->llim[i];
1032 sum = 0;
1033 for (j=ulim ; j<=llim ; j++)
1034 sum += sairow[j];
1035
1036 /* store the row-sum, scaled and normalized for the range of the sum */
1037 buffer[point++] = sum*state->scale / (llim-ulim);
1038 }
1039
1040 if( !last )
1041 return ( (Pointer) buffer ) ;
1042 else {
1043 Delete( state->llim ) ;
1044 Delete( state->ulim ) ;
1045
1046 return ( DeleteFillableSource( state ) ) ;
1047 }
1048 }
1049
1050 /********************************************************************************/
1051 /* */
1052 /* initlist initialises the linked list pointer before the other list functions */
1053 /* call the list. */
1054 /* */
1055 /********************************************************************************/
1056 void initlist()
1057 {
1058 static int status=0;
1059 if (status==0)
1060 start = endl = pp = qq = (node *)malloc(sizeof(node));
1061 status=1;
1062 }
1063
1064 /********************************************************************************/
1065 /* */
1066 /* Function getnode allocates storage for a link list node and returns a */
1067 /* pointer to that node */
1068 /* */
1069 /********************************************************************************/
1070 node *getnode()
1071 {
1072 node *q;
1073 q=(node *)malloc(sizeof(node));
1074 if (q==NULL) exit(66);
1075 return q;
1076 }
1077
1078 /********************************************************************************/
1079 /* */
1080 /* The function insertl inserts the value x at the end of the linked list and */
1081 /* moves the end point by another node */
1082 /* */
1083 /********************************************************************************/
1084 node *insertl(x)
1085 int x;
1086 {
1087 node *q=endl;
1088 endl=getnode();
1089 q->next=endl;
1090 q->index=x;
1091 q->val=0;
1092 /* fprintf(testfile, "insertl: index is %d, val is %d\n", q->index, q->val); */
1093 return q;
1094 }
1095
1096 /********************************************************************************/
1097 /* */
1098 /* The function inss searches the list from the start, stops 3 nodes before the */
1099 /* value x, assigns the value -2000 to the three nodes till the (x+1)th node. */
1100 /* */
1101 /********************************************************************************/
1102 void inss(x)
1103 int x;
1104 {
1105 node *q=start;
1106 int y=x-3;
1107 int z=x;
1108 while (q!=endl)
1109 {
1110 if (q->index == y)
1111 do
1112 {
1113 q->val=-2000;
1114 q=q->next;
1115 /* fprintf(stderr, "Inss(x) index=%d\n",q->index); */
1116 } while (q->index < z);
1117 q=q->next;
1118 }
1119 }
1120
1121
1122 /********************************************************************************/
1123 /* */
1124 /* output_simple_strobe_info routine creates a file pointed by trigger file */
1125 /* which writes the points in the Nap pulse which initiates strobing */
1126 /* */
1127 /********************************************************************************/
1128 void output_simple_strobe_info(state, value)
1129 struct _sai_state *state;
1130 short value;
1131 {
1132
1133 if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on")))
1134 fwrite(&(value), sizeof(short), 1 , trigger_file);
1135 }
1136
1137 /********************************************************************************/
1138 /* */
1139 /* doList_strobe_info creates a file (pointed by "trigger_file") which writes */
1140 /* in binary shorts the points in the Nap pulse which initiates a strobe */
1141 /* It uses a linked list data structure to keep track of the actual time which */
1142 /* causes a strobe, rather than the time when the strobing event occurs. */
1143 /* At stcrit 4 or 5, strobing occurs nwid ms after a nap pulse peak has been */
1144 /* located. */
1145 /* */
1146 /********************************************************************************/
1147 void doList_strobe_info(state, chan)
1148 struct _sai_state *state;
1149 int chan;
1150 {
1151 int storev;
1152
1153 if (state->Stcrit==2)
1154 storev=(state->time)-2;
1155 else if (state->Stcrit==5)
1156 storev=state->def_strobe_candidate_time[chan];
1157 /* fprintf(stderr, "Storev value of doList is = %d\n", storev);} */
1158 else
1159 storev=state->trigtime[chan];
1160
1161 inss(storev);
1162
1163 while (pp->index<=storev) /* (pp!=qq) */
1164 {
1165 if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on")))
1166 fwrite(&(pp->val), sizeof(short), 1, trigger_file);
1167 /* fprintf(testfile, "doList: index is %d, val is %d\n", pp->index, pp->val); */
1168 pp=pp->next;
1169 }
1170 /* pp=qq; */
1171 }
1172
1173 /********************************************************************************/
1174 /* */
1175 /* output_thresh_info writes threshold values (in binary shorts or floats) to */
1176 /* the file opened by the file pointer "thresho_file" */
1177 /* */
1178 /********************************************************************************/
1179 void output_thresh_info(state, chan)
1180 struct _sai_state *state;
1181 int chan;
1182 {
1183 S[chan]=state->thresh[chan];
1184
1185 if (!(!strcmp(state->switch_info, "off") || !strcmp(state->switch_info, "on")))
1186 #ifdef FLOAT
1187 fwrite((S+chan), sizeof(float), 1, thresho_file);
1188 #else
1189 fwrite((S+chan), sizeof(short), 1, thresho_file);
1190 #endif
1191
1192 }
1193
1194 /********************************************************************************/
1195 /* */
1196 /* decay_strobe_threshold decays the strobe threshold at every clock */
1197 /* tick either linearly or exponentially */
1198 /* */
1199 /********************************************************************************/
1200 void decay_strobe_threshold(state, chan)
1201 struct _sai_state *state;
1202 int chan;
1203 {
1204
1205 if (!strcmp(exposwitch, "on"))
1206 state->thresh[chan] *= state->triggerdecay;
1207 else
1208 {
1209 state->thresh[chan] -= state->tlim[chan];
1210 if (state->thresh[chan] < 0)
1211 state->thresh[chan]=0;
1212 }
1213 }
1214
1215 /********************************************************************************/
1216 /* */
1217 /* This routine prints debugging information of the local max search */
1218 /* with timeout (i.e. stcrit_ai=5) */
1219 /* */
1220 /********************************************************************************/
1221 void print_trigger_debugging_info(state, chan, strobe_lag, stage)
1222 struct _sai_state *state;
1223 int chan;
1224 int strobe_lag;
1225 int stage;
1226 {
1227 if (stage==1)
1228 {
1229 if (strobe_lag > state->def_strobe_candidate_time[chan])
1230 val[chan]=strobe_lag-state->def_strobe_candidate_time[chan];
1231 if (strobe_lag == state->def_strobe_candidate_time[chan])
1232 val[chan]=0;
1233
1234 if (!strcmp(state->switch_info,"off"));
1235 else
1236 fprintf(infotxt_file, "Chan %d: :: Actual time of trigger: s_lag=%d, trigtime=%d\n", chan, strobe_lag, state->def_strobe_candidate_time[chan]);
1237 }
1238
1239 /* 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); */
1240
1241
1242 if (stage==2)
1243 {
1244 if (!strcmp(state->switch_info,"off"));
1245 else
1246 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]));
1247 }
1248
1249 if (stage==3)
1250 {
1251 if (!strcmp(state->switch_info,"off"));
1252 else
1253 fprintf(infotxt_file, "Chan %d: Loc Max is %d Time is %d\n", chan, state->trigheight[chan], state->trigtime[chan] );
1254 }
1255
1256 if (stage==4)
1257 {
1258 if (!strcmp(state->switch_info,"off"));
1259 else
1260 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]);
1261 }
1262 if (stage==5)
1263 {
1264 if (!strcmp(state->switch_info,"off"));
1265 else
1266 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]);
1267 }
1268
1269
1270 }
1271