tomwalters@0: /*
tomwalters@0:     Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
tomwalters@0:     ===========================================================================
tomwalters@0: 
tomwalters@0:     Permission to use, copy, modify, and distribute this software without fee 
tomwalters@0:     is hereby granted for research purposes, provided that this copyright 
tomwalters@0:     notice appears in all copies and in all supporting documentation, and that 
tomwalters@0:     the software is not redistributed for any fee (except for a nominal shipping 
tomwalters@0:     charge). Anyone wanting to incorporate all or part of this software in a
tomwalters@0:     commercial product must obtain a license from the Medical Research Council.
tomwalters@0: 
tomwalters@0:     The MRC makes no representations about the suitability of this 
tomwalters@0:     software for any purpose.  It is provided "as is" without express or implied 
tomwalters@0:     warranty.
tomwalters@0:  
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:     corti.c
tomwalters@0:     =======
tomwalters@0: 
tomwalters@0:     attempt at model of organ of corti
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     Copyright (c), 1989  The Medical Research Council, Applied Psychology Unit.
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     Author  : John Holdsworth
tomwalters@0:     Written : 30th November, 1988.
tomwalters@0: 
tomwalters@0:     Edited  : 
tomwalters@0: 
tomwalters@0:               Christian Giguere, March 1994. 
tomwalters@0: 	      - Changed option "scale_at" into "gain_at".
tomwalters@0: 	      - That option is now passed to newCorti() as a double instead of an integer
tomwalters@0: 	      - In Corti(), this meant multiplying *optr by state->scale
tomwalters@0:               - To locate changes, search for "CG"
tomwalters@0: 
tomwalters@0: */
tomwalters@0: #include <stdio.h>
tomwalters@0: #include  <math.h>
tomwalters@0: 
tomwalters@0: #include "stitch.h"
tomwalters@0: #include "source.h"
tomwalters@0: #include  "corti.h"
tomwalters@0: #include   "calc.h"
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define DEBUG
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #ifndef  lint
tomwalters@0: static char *sccs_id = "@(#)corti.c	1.15     J. Holdsworth (MRC-APU)  5/31/91" ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: /* calmB is the mB of the magnitude where the slope of temporal recovery is calculated */
tomwalters@0: /* 5497 corresponds to a magnitude of 500 units */
tomwalters@0: 
tomwalters@0: static double calmB = 5497. ;
tomwalters@0: 
tomwalters@0: struct _cell_state {
tomwalters@0: StoreType microphonic, rapid_limit, absolute_limit ;
tomwalters@0: ScalarType rapid_rise, fast_rise, rapid_decay, fast_decay ;
tomwalters@0: double compensate ;             /* Changed to double from StoreType (mha:22/1/93) */
tomwalters@0: } ;
tomwalters@0: 
tomwalters@0: struct _corti_state {
tomwalters@0: struct _cell_state *cell_states ;
tomwalters@0: StoreType scale ;
tomwalters@0: ScalarType t_lat ;
tomwalters@0: unsigned cells, times ;
tomwalters@0: } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* global variable for spectral tilt control to compensate for differences */
tomwalters@0: /* between channels in terms of mean firing rate                           */
tomwalters@0: 
tomwalters@0: double compensate = 1. ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: char *newCorti( cells, samplerate, center_frequencies, bandwidth_function, t_rise, t_lat, t_rapid, t_fast, prop, absthresh, times, scale )
tomwalters@0: int cells ;
tomwalters@0: double samplerate, *center_frequencies, (*bandwidth_function)(), t_rise, t_lat, t_rapid, t_fast, prop, absthresh ;
tomwalters@0: int times ;
tomwalters@0: double scale ;   /* CG: 22/03/94 */
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _corti_state *, state ) ;
tomwalters@0:     struct _cell_state *cell_ptr ;
tomwalters@0:     double ratio, k_rise ;
tomwalters@0:     double expected_rate, expected_reference ;
tomwalters@0: 
tomwalters@0:     /* parameters that are constant across cells */
tomwalters@0: 
tomwalters@0:     state->t_lat = SCALAR( t_lat  / ( samplerate * times ) ) ;
tomwalters@0: 
tomwalters@0:     state->cells  = cells ;
tomwalters@0:     state->times  = times ;
tomwalters@0:     state->scale  = scale ; 
tomwalters@0: 
tomwalters@0:     state->cell_states = NewArray( struct _cell_state, state->cells, "corti.c for cell states" ) ;
tomwalters@0: 
tomwalters@0:     for( cell_ptr = state->cell_states ; cell_ptr < state->cell_states + cells ; cell_ptr++ ) {
tomwalters@0: 
tomwalters@0: 	/* corti states */
tomwalters@0: 
tomwalters@0: 	cell_ptr->microphonic    = SCALAR( absthresh ) ;
tomwalters@0: 	cell_ptr->rapid_limit    = SCALAR( absthresh ) ;
tomwalters@0: 	cell_ptr->absolute_limit = SCALAR( absthresh ) ;
tomwalters@0: 
tomwalters@0: 	/* recovery time constants */
tomwalters@0: 
tomwalters@0: 	cell_ptr->rapid_decay = SCALAR( bandwidth_function( center_frequencies[cell_ptr-state->cell_states] ) *
tomwalters@0: 				 t_rapid * calmB / ( calmB - absthresh ) / samplerate /* / ( 1. - prop ) */ ) ;
tomwalters@0: 
tomwalters@0: 	cell_ptr->fast_decay  = SCALAR( 1. / t_fast / samplerate ) ;
tomwalters@0: 
tomwalters@0: 	 /* adaptation time constants */
tomwalters@0: 
tomwalters@0: 	ratio = prop / ( 1. - prop ) / cell_ptr->rapid_decay * cell_ptr->fast_decay ;
tomwalters@0: 
tomwalters@0: 	if( t_rise > 0 )
tomwalters@0: 	    k_rise = t_rise / ( samplerate * times ) ;
tomwalters@0: 	else
tomwalters@0: 	    k_rise = 1. ;
tomwalters@0: 
tomwalters@0: 	cell_ptr->rapid_rise = SCALAR( k_rise * 1.    / ( 1. + ratio ) ) ;
tomwalters@0: 	cell_ptr->fast_rise  = SCALAR( k_rise * ratio / ( 1. + ratio ) ) ;
tomwalters@0: 
tomwalters@0: 	/* introduce variable amount of compensation w.r.t expected firing rate variation */
tomwalters@0: 
tomwalters@0: 	expected_rate = UNSCALE( cell_ptr->rapid_decay ) * ( calmB - absthresh ) ;
tomwalters@0: 
tomwalters@0: 	if( cell_ptr == state->cell_states )
tomwalters@0: 	    expected_reference = expected_rate ;
tomwalters@0: 
tomwalters@0: 	/* New formula for computing cell_ptr->compensate  (mha: 20/1/93)  */
tomwalters@0: 	/* Parameters chosen to flatten pulse-train spectrum.              */
tomwalters@0: 	/* _cell_state member "compensate" changed to double.              */
tomwalters@0: 	/* The exponent used is 1.128*compensate so that the given         */
tomwalters@0: 	/* parameter (compensate_at) can be set as "on" (ie 1.0) for       */
tomwalters@0: 	/* optimal compensation (instead of a magic number)                */
tomwalters@0: 	/* The value for "compensate_at=off", (calmB - absthresh) was      */
tomwalters@0: 	/* chosen as similar to the original, and because the results      */
tomwalters@0: 	/* were appropriate for an un-compensated spectrum.                */
tomwalters@0: 
tomwalters@0: 	if (compensate)
tomwalters@0: 	    cell_ptr->compensate =  pow(center_frequencies[cell_ptr - state->cell_states], 1.128*compensate) + 800. ;
tomwalters@0: 	else
tomwalters@0: 	    cell_ptr->compensate =  calmB - absthresh ;
tomwalters@0: 
tomwalters@0: 	/* Originally, compensate used the SCALAR macro, and so it failed  */
tomwalters@0: 	/* for a floating-point version (defined in calc.h).               */
tomwalters@0: 	/* To fix this, the new compensate formula is scaled up using the  */
tomwalters@0: 	/* original integer version of the UNSCALE macro (see calc.h).     */
tomwalters@0: 
tomwalters@0: #ifdef FLOAT
tomwalters@0: 	cell_ptr->compensate /=  (float) ( 1l << 16l ) ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #ifdef DEBUG
tomwalters@0: 	(void) fprintf( stderr, "rapid_decay:%f fast_decay:%f ratio:%f k_rise:%f rapid_rise:%f fast_rise:%f compensation:%f\n",
tomwalters@0: 			  UNSCALE(cell_ptr->rapid_decay), UNSCALE(cell_ptr->fast_decay),
tomwalters@0: 	    ratio, k_rise, UNSCALE(cell_ptr->rapid_rise), UNSCALE(cell_ptr->fast_rise),
tomwalters@0: 			1./UNSCALE(cell_ptr->compensate) ) ;
tomwalters@0: #endif
tomwalters@0:     }
tomwalters@0: #ifdef DEBUG
tomwalters@0:     (void) fprintf( stderr, "times:%d t_lat:%f\n", times, UNSCALE( state->t_lat ) ) ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:     return( (char *) state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: int Corti( state_ptr, bytes, out, end, in )
tomwalters@0: char *state_ptr ;
tomwalters@0: ByteCount *bytes ;
tomwalters@0: DataType *out, *end, *in ;
tomwalters@0: {
tomwalters@0:     register struct _corti_state *state = (struct _corti_state *) state_ptr ;
tomwalters@0:     register struct _cell_state *cell_ptr, *cell_end = state->cell_states + state->cells ;
tomwalters@0:     register DataType *iptr, *optr, *eptr = end ;
tomwalters@0:     register StoreType delta, old_last_mic ;
tomwalters@0:     register int time ;
tomwalters@0: 
tomwalters@0:     while( out < eptr ) {
tomwalters@0: 
tomwalters@0: 	for( time=0 ; time < state->times ; time++ ) {
tomwalters@0: 
tomwalters@0: 	    cell_ptr = state->cell_states ;
tomwalters@0: 
tomwalters@0: 	    iptr = in  ;
tomwalters@0: 	    optr = out ;
tomwalters@0: 
tomwalters@0: 	    /* raise threhsolds */
tomwalters@0: 
tomwalters@0: 	    for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ ) {
tomwalters@0: 
tomwalters@0: 		if( time == 0 )
tomwalters@0: 		    *optr = 0 ;
tomwalters@0: 
tomwalters@0: 		delta = *iptr++ - DESCALE( cell_ptr->microphonic ) ;
tomwalters@0: 
tomwalters@0: 		if( delta > 0 ) {
tomwalters@0: 
tomwalters@0: 		    cell_ptr->microphonic += delta * cell_ptr->rapid_rise ;
tomwalters@0: 		    cell_ptr->rapid_limit += delta * cell_ptr->fast_rise  ;
tomwalters@0: 
tomwalters@0: 		    *optr                 += delta ;
tomwalters@0: 		}
tomwalters@0: 
tomwalters@0: 		optr++ ;
tomwalters@0: 	    }
tomwalters@0: 
tomwalters@0: 	    /* leak thresholds symetrically across channels */
tomwalters@0: 
tomwalters@0: 	    old_last_mic = state->cell_states->microphonic ;
tomwalters@0: 
tomwalters@0: 	    for( cell_ptr = state->cell_states+1 ; cell_ptr < cell_end ; cell_ptr++ ) {
tomwalters@0: 
tomwalters@0: 		delta = DESCALE( old_last_mic - cell_ptr->microphonic ) * state->t_lat ;
tomwalters@0: 		old_last_mic = cell_ptr->microphonic ;
tomwalters@0: 
tomwalters@0: 		if( cell_ptr-1 != state->cell_states )
tomwalters@0: 		    (cell_ptr-1)->microphonic -= delta ;
tomwalters@0: 
tomwalters@0: 		if( cell_ptr+1 != cell_end )
tomwalters@0: 		    (cell_ptr  )->microphonic += delta ;
tomwalters@0: 	    }
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0: 	/* decay potentails */
tomwalters@0: 
tomwalters@0: 	for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ )
tomwalters@0: 	    cell_ptr->microphonic -= DESCALE( cell_ptr->microphonic - cell_ptr->rapid_limit ) * cell_ptr->rapid_decay ;
tomwalters@0: 
tomwalters@0: 	for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ )
tomwalters@0: 	    cell_ptr->rapid_limit -= DESCALE( cell_ptr->rapid_limit - cell_ptr->absolute_limit ) * cell_ptr->fast_decay ;
tomwalters@0: 
tomwalters@0: 	/* generate output */
tomwalters@0: 
tomwalters@0: 	iptr = in ;
tomwalters@0: 	optr = out ;
tomwalters@0: 
tomwalters@0: 	if( state->scale >= 0 )                                               /* CG: 22/3/94 */
tomwalters@0: 	    for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ ) {
tomwalters@0: 		delta = SCALE( *iptr++ ) - cell_ptr->microphonic ;
tomwalters@0: 		if( delta > 0 )
tomwalters@0: 		    *optr++ = state->scale * delta / cell_ptr->compensate ;   /* CG: 22/3/94 */
tomwalters@0: 		else
tomwalters@0: 		    *optr++ = 0 ;
tomwalters@0: 	    }
tomwalters@0: 	else if( state->scale < 0 )
tomwalters@0: 	    switch( (int) state->scale ) {
tomwalters@0: 
tomwalters@0: 		/* output internal states */
tomwalters@0: 
tomwalters@0: 		case -1 :
tomwalters@0: 		    for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ )
tomwalters@0: 			*optr++ = DESCALE( cell_ptr->microphonic ) ;
tomwalters@0: 
tomwalters@0: 		    break ;
tomwalters@0: 
tomwalters@0: 		case -2 :
tomwalters@0: 		    for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ )
tomwalters@0: 			*optr++ = DESCALE( cell_ptr->rapid_limit ) ;
tomwalters@0: 
tomwalters@0: 		    break ;
tomwalters@0: 
tomwalters@0: 		default :
tomwalters@0: 		    for( cell_ptr = state->cell_states ; cell_ptr < cell_end ; cell_ptr++ )
tomwalters@0: 			*optr++ = ( *iptr++ - DESCALE( cell_ptr->microphonic ) ) * -state->scale ;
tomwalters@0: 
tomwalters@0: 		    break ;
tomwalters@0: 
tomwalters@0: 	    }
tomwalters@0: 	/* else leave value as is */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 	/* next frame */
tomwalters@0: 
tomwalters@0: 	in  += state->cells ;
tomwalters@0: 	out += state->cells ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* return amount processed */
tomwalters@0: 
tomwalters@0:     return ( *bytes ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /* for compatability */
tomwalters@0: 
tomwalters@0: int OldCorti( state_ptr, cells, in, out )
tomwalters@0: char *state_ptr ;
tomwalters@0: unsigned cells ;
tomwalters@0: DataType *in, *out ;
tomwalters@0: {
tomwalters@0:     ByteCount bytes = ToBytes( DataType, cells ) ;
tomwalters@0: 
tomwalters@0:     return ( ToPoints( DataType, Corti( state_ptr, &bytes, out, out+cells, in ) ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: void closeCorti( state )
tomwalters@0: struct _corti_state *state ;
tomwalters@0: {
tomwalters@0:     Delete( state->cell_states ) ;
tomwalters@0:     Delete( state ) ;
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* experimental corti  simulation */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: struct _corti2_state { int *working, *falls, slope, scale ; short *forward, *backward ; } ;
tomwalters@0: 
tomwalters@0: char *newCorti2( cells, samplerate, center_frequencies, bandwidth_function, bandwidth_scalar, t_slope, scale, forward, backward )
tomwalters@0: int cells ;
tomwalters@0: double samplerate, *center_frequencies, (*bandwidth_function)(), bandwidth_scalar, t_slope, scale ;
tomwalters@0: short *forward, *backward ;
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _corti2_state *, state ) ;
tomwalters@0:     int i ;
tomwalters@0: 
tomwalters@0:     state->working =       NewArray( int, cells, "corti.c for works" ) ;
tomwalters@0:     state->falls   = NewZeroedArray( int, cells, "corti.c for falls" ) ;
tomwalters@0: 
tomwalters@0:     for( i=0 ; i < cells ; i++ )
tomwalters@0: 	state->falls[i] = bandwidth_function( center_frequencies[i] ) * bandwidth_scalar /
tomwalters@0: 				 samplerate * t_slope ;
tomwalters@0: 
tomwalters@0:     state->forward  = forward  ;
tomwalters@0:     state->backward = backward ;
tomwalters@0: 
tomwalters@0:     state->scale = scale ;
tomwalters@0: 
tomwalters@0:     return( ( char * ) state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: int Corti2( state_ptr, cells, in, out )
tomwalters@0: char *state_ptr ;
tomwalters@0: unsigned cells ;
tomwalters@0: short *in, *out ;
tomwalters@0: {
tomwalters@0:     struct _corti2_state *state ;
tomwalters@0:     register int i, when ;
tomwalters@0: 
tomwalters@0:     state = ( struct _corti2_state * ) state_ptr ;
tomwalters@0: 
tomwalters@0:     for( i=0 ; i<cells ; i++ ) {
tomwalters@0: 	state->working[i] -= state->falls[i] ;
tomwalters@0: 	out[i] = in[i] - state->working[i] ;
tomwalters@0: 	if( out[i] > 0 )
tomwalters@0: 	    state->working[i] += out[i] ;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     when = 0 ;
tomwalters@0: 
tomwalters@0:     for( i=1 ; i<cells ; i++ )
tomwalters@0: 	if( state->working[i] < state->working[ when ] + state->forward[  i - when ] ) {
tomwalters@0: 	    state->working[i] = state->working[ when ] + state->forward[  i - when ] ;
tomwalters@0: 	    out[i] = 0 ;
tomwalters@0: 	}
tomwalters@0: 	else
tomwalters@0: 	    when = i ;
tomwalters@0: 
tomwalters@0:     when = cells-1 ;
tomwalters@0: 
tomwalters@0:     for( i=cells-2 ; i>=0 ; i-- )
tomwalters@0: 	if( state->working[i] < state->working[ when ] + state->backward[ when - i ] ) {
tomwalters@0: 	    state->working[i] = state->working[ when ] + state->backward[ when - i ] ;
tomwalters@0: 	    out[i] = 0 ;
tomwalters@0: 	}
tomwalters@0: 	else
tomwalters@0: 	    when = i ;
tomwalters@0: 
tomwalters@0:     for( i=0 ; i<cells ; i++ )
tomwalters@0: 	if( state->scale != 0 )
tomwalters@0: 	    out[i] = out[i] * state->scale / state->falls[i] ;
tomwalters@0: 	else
tomwalters@0: 	    out[i] = state->working[i] ;
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: void closeCorti2( state )
tomwalters@0: struct _corti2_state *state ;
tomwalters@0: {
tomwalters@0:     stitch_free( ( char * ) state->working ) ;
tomwalters@0:     stitch_free( ( char * ) state->falls ) ;
tomwalters@0: 
tomwalters@0:     Delete( state ) ;
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: 
tomwalters@0: /* obsolete code stored here for now.. stitch entry points */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: struct _stx_state { Source source ; unsigned cells ; char *corti ; int (*model)() ; } ;
tomwalters@0: 
tomwalters@0: static void stx_callback( state, bytes, buffer )
tomwalters@0: struct _stx_state *state ;
tomwalters@0: ByteCount bytes ;
tomwalters@0: DataType *buffer ;
tomwalters@0: {
tomwalters@0:     int points      = ToPoints( DataType, bytes ) ;
tomwalters@0:     DataType *input = PullItems( state->source, points, DataType ) ;
tomwalters@0:     int pointer ;
tomwalters@0: 
tomwalters@0:     for( pointer=0 ; pointer < points ; pointer += state->cells )
tomwalters@0: 	(void) state->model( state->corti, state->cells, input+pointer, buffer+pointer ) ;
tomwalters@0: 
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: Source Stx( source, chans, samplerate, center_frequencies, bandwidth_function, bandwidth_scalar, t_rise, t_lat, t_rapid, absthresh, times, scale )
tomwalters@0: Source source ;
tomwalters@0: int chans ;
tomwalters@0: double samplerate, *center_frequencies, (*bandwidth_function)(), bandwidth_scalar, t_rise, t_lat, t_rapid, absthresh ;
tomwalters@0: int times, scale ;
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _stx_state *, state ) ;
tomwalters@0: 
tomwalters@0:     state->source = source ;
tomwalters@0:     state->cells  = chans ;
tomwalters@0: 
tomwalters@0:     state->corti = newCorti( chans, samplerate, center_frequencies, bandwidth_function, t_rise, t_lat, t_rapid, absthresh, times, scale ) ;
tomwalters@0:     state->model = OldCorti ;
tomwalters@0: 
tomwalters@0:     return ( stdAutoSource( ( Pointer ) state, stx_callback ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: Source Stx2( source, chans, samplerate, center_frequencies, bandwidth_function, bandwidth_scalar, rapid_rises, scale )
tomwalters@0: Source source ;
tomwalters@0: int chans ;
tomwalters@0: double samplerate, *center_frequencies, (*bandwidth_function)(), bandwidth_scalar, rapid_rises, scale ;
tomwalters@0: {
tomwalters@0:     DeclareNew( struct _stx_state *, state ) ;
tomwalters@0:     short *forward = ( short * ) 0, *backward = ( short * ) 0 ;
tomwalters@0: 
tomwalters@0:     state->source = source ;
tomwalters@0:     state->cells  = chans ;
tomwalters@0: 
tomwalters@0:     state->corti = newCorti2( chans, samplerate, center_frequencies, bandwidth_function, bandwidth_scalar, rapid_rises, scale, forward, backward ) ;
tomwalters@0:     state->model = Corti2 ;
tomwalters@0: 
tomwalters@0:     return ( stdAutoSource( ( Pointer ) state, stx_callback ) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: #endif