diff model/corti.c @ 0:5242703e91d3 tip

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

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