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

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

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