view wdf/meddis.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.
 
    THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
    ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
    A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY 
    DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
    AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 
    OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/

/*  
    Acknowledgment:
    ==============

    The source code provided in this file was originally developed by 
    Christian Giguere as part of a Ph.D degree at the Department of
    Engineering of the University of Cambridge from April 1990 to
    November 1993. The code was subsequently adapted under a grant
    from the Hearing Research Trust for full compatibility with 
    AIM Release 6.15.

    Christian Giguere 25/03/94

*/

/*
    ===========================================================
    meddis.c
    ===========================================================

    Implementation of Ray Meddis haircell model.

    Author        : Christian Giguere
    First written : 10th September, 1991
    Last edited   : 07th March, 1994

    References:
    (1) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
    (2) R.Meddis et al. (1990). JASA 87(4): 1813-1816.
    (3) R.Meddis (1988). JASA 83(3): 1056-1063.

    Note - This is a complete re-write of the original file 
	   "haircell.c" Release 3 (20th September 1988) from  
	   J. Holdsworth and the Applied Psychology Unit. 
    ============================================================
*/


/***** includes *****/

#include <math.h>
#include <stdio.h>
#include "stitch.h"
#include "source.h"
#include "calc.h"
#include "calc_tl.h"
#include "meddis.h"


/***** defines *****/

#if 0
#define   _DEBUG_
#define   _OVERFLOW_
#endif

#if 1                                                /* define numerical implementation */
#define   _WDF_                                      /* Giguere and Woodland (1994) */
#else
#define   _FORWARD_DIFFERENCE_                       /* Meddis et al. (1990) */
#endif

#if 0  
#define   _CLAMPING_                                 /* clamping of reservoirs */
#endif
 
#define   NUMBER_OF_COEFFS     (           12 )
#define   NUMBER_OF_STATES     (           12 )


/* Feedback parameter */

#define   InputGain_max        (         24.0 )      /* max coupling gain in dB */


/* Meddis et al. (1990) model parameters */

#define   A_medium             (         10.0 )      /* medium-spontaneous rate fibre  */
#define   B_medium             (       3000.0 )
#define   g_medium             (       1000.0 )

#define   A_high               (          5.0 )      /* high-spontaneous rate fibre */
#define   B_high               (        300.0 )
#define   g_high               (       2000.0 )

#define   y_value              (         5.05 )      /* replenishment rate */
#define   l_value              (         2500 )      /* rate of loss from the cleft */
#define   r_value              (         6580 )      /* rate of return from the cleft */
#define   x_value              (        66.31 )      /* rate of release from w store */ 
#define   m_value              (          1.0 )      /* max number of quanta */
#define   h_value              (       50000. )      /* firing-rate factor */


/***** private data structures *****/

typedef struct _haircell_info HaircellInfo ;
typedef struct _haircell_bank HaircellBank ;

struct _haircell_info {
  int        number_of_states ;  
  double     output_gain ;
  StateType  *states ;
  } ;

struct _haircell_bank {
  int           channels ;
  double        *inputGain_ratio ;
  double        (*update_proc)() ;                  /* procedure to update cells input gain */
  int           number_of_coeffs ;
  CoeffType     *coeffs ;
  HaircellInfo  **cells ;
  } ;


/***** external variables *****/

double *InputGain_ratio ;


/***** functions *****/

/*******************************************************************************
* name:                 function:
*
* NewHaircells()        Set-up function for the implementation of a bank of inner
*                       haircells based on the model of Ray Meddis. It generates the 
*                       multiplier coefficients for each haircell filter and 
*                       initializes the filter states. It is called by function 
*                       ``Meddis()'' in file ``model.c''. It returns a pointer to a 
*                       structure containing all the relevant information for the 
*                       computation of the haircell filters. 
*
*                       The choice of the haircell type (medium or high-spontaneous
*                       rate fibre) is made at run time using options ``fibre_med''.
*
*                       The choice of the numerical implementation (Wave Digital 
*                       Filtering or Forward Difference algorithm) is made at 
*                       compile time using a macro substitution ``#define'' (see above).
*
*                       The choice of clamping the reservoirs' contents to non-
*                       negative values is made at compile time using a macro
*                       substitution ``#define'' (see above).
*
* DoHaircells()         Realization of the bank of haircell filters. It computes the
*                       output for each haircell for a specified number of input 
*                       points and keeps track of the filters' states. It returns
*                       the size of the output data in bytes (per point).
*
* CloseHaircells()      Deletes all private data structures and arrays of the Haircell 
*                       filters upon comletion of filtering.
********************************************************************************/


/*************************** NewHaircells() *********************************/

Pointer NewHaircells( channels, samplerate, output_gain, coupling, fibreType )
    int    channels ;
    double samplerate, output_gain ;
    double coupling ;
    char   *fibreType ;
{
    DeclareNew( HaircellBank *, bank ) ;
    HaircellInfo  *cell_info ;
    CoeffType     *cf ;
    StateType     *st ;  
    double        ydt, ldt, rdt, xdt, gdt ; 
    double        kdt_init, q_init, c_init, w_init, cell_scaling ; 
    double        A_value, B_value, g_value ;
    double        r1, r2, r3, r4, r03 ;
    double        update_inputGain() ;
    int           channel ;

   /*** scaling of input data ***/
    cell_scaling = coupling ;

   /*** allocate fibre-dependent haircell parameters ***/
    if( strncmp( fibreType, "high", 3 ) == 0 ) {
	A_value = A_high ;
	B_value = B_high ;
	g_value = g_high ;
    }
     else {
	A_value = A_medium ;
	B_value = B_medium ;
	g_value = g_medium ;
    }

#ifdef _DEBUG_
    fprintf( stderr, "Meddis Haircell fibre type=%s\n", fibreType ) ;
    fprintf( stderr, "A=%.2f B=%.2f g=%.2f\n", A_value, B_value, g_value ) ;
    fprintf( stderr, "Input scaling factor=%.4f\n", cell_scaling ) ;
#endif


   /*** scale model parameters for difference equations ***/ 
    ydt = y_value / samplerate ;
    ldt = l_value / samplerate ;
    rdt = r_value / samplerate ;
    xdt = x_value / samplerate ;
    gdt = g_value / samplerate ;


   /*** calculate initial values for infinite history of silence ***/
    kdt_init = gdt * A_value / ( A_value + B_value ) ;
    q_init   = m_value / ( 1. + kdt_init / ydt * ( 1. - rdt / ( ldt + rdt ) ) ) ;
    c_init   = q_init * kdt_init / ( ldt + rdt ) ;
    w_init   = rdt / xdt * c_init ;    


   /*** allocate and store multiplier coefficients ***/
    bank->number_of_coeffs = NUMBER_OF_COEFFS ;  
    cf = bank->coeffs = NewArray( CoeffType, bank->number_of_coeffs, 
				  "haircell_tl.c for coefficients" ) ;
#ifdef _WDF_
    output_gain = MEDDIS_SCALE * output_gain * h_value / ( 2. * r_value ) ;
    cf[0] = m_value * y_value * output_gain ;                     /* M/Zy (normalized) */
    cf[1] = A_value / cell_scaling ;
    cf[2] = ( A_value + B_value ) / cell_scaling ; 
    cf[3] = g_value ;

    /* Adaptor B0 */ 
    r1 = y_value ;                                                /* ( 1/Zy ) */
    r2 = 2. * samplerate ;                                        /* ( 1/Zcq ) */
    r3 = r03 = r1 + r2 ;
    cf[5] = r1 / r3 ;                                             /* gamma01 */
    cf[4] = r3 / cf[3] ;                                          /* needed to update gamma11 */

    /* Adaptor B1*/ 
    r1 = kdt_init * samplerate ;                                  /* ( 1/Zk ) initial value only */
    r2 = r03 ;
    r3 = r1 + r2 ;
    cf[6] = 2. * r1 / r3 ;                                        /* gamma11 updated at each sample */

    /* Adaptor B2 */ 
    r1 = l_value ;                                                /* ( 1/Zl ) */
    r2 = r_value ;                                                /* ( 1/Zr ) */
    r3 = 2. * samplerate ;                                        /* ( 1/Zcc ) */
    r4 = r1 + r2 + r3 ;
    cf[7] = r1 / r4 ;                                             /* gamma21 */
    cf[8] = r2 / r4 ;                                             /* gamma22 */

    /* Adaptor B3 */
    r1 = x_value ;                                                /* ( 1/Zx ) */
    r2 = 2. * samplerate ;                                        /* ( 1/Zcw ) */
    r3 = r1 + r2 ;
    cf[9] = 0.5 * r1 / r3 ;                                       /* 0.5 * gamma31 */

#else
    output_gain = MEDDIS_SCALE * output_gain * h_value ;
    cf[0] = m_value * output_gain ;
    cf[1] = A_value / cell_scaling ;
    cf[2] = ( A_value + B_value ) / cell_scaling ;
    cf[3] = gdt ;
    cf[4] = ydt ;
    cf[5] = ldt ;
    cf[6] = rdt ;
    cf[7] = xdt ;
#endif


   /*** loop through each haircell ***/
    bank->channels = channels ;
    bank->cells = NewArray( HaircellInfo *, bank->channels, "haircell.c for cell states" ) ;
    bank->inputGain_ratio = InputGain_ratio = NewArray( double, bank->channels, 
							"haircell_tl.c for cells input gain ratio" ) ;
    bank->update_proc = update_inputGain ;

    for( channel = 0 ; channel < bank->channels ; channel++ ) {
	 bank->inputGain_ratio[ channel ] = 0.5 ;            /* in absence of feedback, Gain = 1 */
	 cell_info = bank->cells[ channel ] = New( HaircellInfo * ) ;

	/*** allocate and initialise states ***/
	 cell_info->number_of_states = NUMBER_OF_STATES ;
	 st = cell_info->states = NewZeroedArray( StateType, cell_info->number_of_states, 
						  "haircell_tl.c for states" ) ;

#ifdef _WDF_
	st[1] =  x_value * w_init * output_gain ;               /* Ixn (normalized) */
	st[2] =  2. * q_init * samplerate * output_gain ;       /* Ixn(t-T)-b02 (normalized) */
	st[3] = -2. * c_init * samplerate * output_gain ;       /* b23 (normalized) */
	st[4] = -2. * w_init * samplerate * output_gain ;       /* b33 (normalized) */
	cell_info->output_gain = output_gain ;                  /* obsolete */

#else
	st[0] = kdt_init ;
	st[1] = q_init * output_gain ;
	st[2] = c_init * output_gain ;
	st[3] = w_init * output_gain ;
	cell_info->output_gain = output_gain ;                  /* obsolete */
#endif

    }

    return( ( Pointer ) bank ) ;
}


/************************** DoHaircells() *********************************/

int DoHaircells( bank_ptr, bytes, output_ptr, end_output_ptr, input_ptr ) 
    Pointer    *bank_ptr ;
    ByteCount  *bytes ;
    DataType   *output_ptr, *end_output_ptr, *input_ptr ;
{
    register HaircellBank *bank = ( HaircellBank * ) bank_ptr ;
    register HaircellInfo *cell_info ;
    register DataType     *input, *output ;
    register StateType    *st ;
    register StateType    in, out ;
    register CoeffType    *cf ;
    register double       inputGain ;
    register int          channel = bank->channels ;
    register int          point, npoints = ( end_output_ptr - output_ptr ) / bank->channels ;
#ifdef _WDF_
    register StateType    an, a0 ;
#else
    register StateType    replenish, eject, loss, reuptake, reprocess ;
#endif

    /*** for all channels ***/

#ifdef _DEBUG_
     fprintf( stderr, "DoHaircells() = %d points X %d channels\n", npoints, bank->channels ) ;
#endif

     cf = bank->coeffs ;

     while( channel-- ) {

	cell_info = bank->cells[ channel ] ;
	st = cell_info->states ;
	input = input_ptr + channel ;                
	output = output_ptr + channel ;
	inputGain = bank->update_proc( bank->inputGain_ratio[ channel ] ) ; 
	cf[10] = cf[1] / inputGain ;                                    /* (A/p) */
	cf[11] = cf[2] / inputGain ;                                    /* (A+B)/p */

       /*** for each channel ***/

	point = npoints ;  
	while( point-- ) { 

#ifdef _WDF_
	  /*** compute compressive permeability function k ***/
	   in = ( StateType ) *input ;                                  /* (BM vibration) */
	   if( in > -cf[10] )
	      st[0] = ( in + cf[10] ) / ( in + cf[11] ) ;               /* kn(t) / g */ 
	   else
	      st[0] = 0. ;

	  /*** update gamma11 ***/
	   cf[6] = st[0] / ( st[0] + cf[4] ) ;
	   cf[6] += cf[6] ;                                             /* gamma11 */

	  /*** compute wdf ***/

	  /* Adaptor  B0 */
	  st[5] = cf[0] + st[1] + st[2] ;                               /* -b03=a12 (normalized) */

	  /* Adaptor B1 */
	  st[6] = cf[6] * st[5] ;                                       /* -b11=2Ikn (normalized) */
	  st[7] = st[6] - st[5] ;                                       /* b12=-a03 (normalized) */

	  /* Adaptor B0 */
	  st[2] = cf[0] + st[1] - st[7] + cf[5] * ( st[7] - st[5] ) ;   /* Ixn(t-T)-b02 (normalized) */

	  /* Adaptor B2 */
	  an = st[6] - st[3] ;                                          /* a24 (normalized) */ 
	  a0 = an - st[3] ;                                             /* a20 (normalized) */
	  st[8] = cf[8] * a0 ;                                          /* -b22=2Irn (normalized) */
	  st[3] = cf[7] * a0 + st[8] - an ;                             /* b23 (normalized) */

	  /* Adaptor B3 */
	  an = st[8] - st[4] ;                                          /* a33 (normalized) */ 
	  st[1] = cf[9] * ( an - st[4] ) ;                              /* Ixn(t)=-0.5*b31 (normalized) */
	  st[4] = st[1] + st[1] - an ;                                  /* b32 (normalized) */

#ifdef _CLAMPING_
	  /*** clamping of reservoirs to non-negative quantities ***/
	  /*** not available yet ***/
#endif     

	  /*** output ***/
	   out = st[8] ;

#else

	  /*** compute change quantities ***/
	   replenish = cf[4] * ( cf[0] - st[1] ) ;                     /* y(M-q) */
	   eject     = st[0] * st[1] ;                                 /* kq */
	   loss      = cf[5] * st[2] ;                                 /* lc */ 
	   reuptake  = cf[6] * st[2] ;                                 /* rc */
	   reprocess = cf[7] * st[3] ;                                 /* xw */

	  /*** compute compressive permeability function kn_1dt ***/
	   in = ( StateType ) *input ;                                 /* (BM vibration) */
	   if( in > -cf[10] )
	      st[0] = cf[3] * ( in + cf[10] ) / ( in + cf[11] ) ;      /* kn_1(t)dt */ 
	   else
	      st[0] = 0. ;

	  /*** update reservoir quantities ***/
	   st[1] += replenish - eject + reprocess ;                    /* q */
	   st[2] += eject - loss - reuptake ;                          /* c */
	   st[3] += reuptake - reprocess ;                             /* w */
 
#ifdef _CLAMPING_
	  /*** clamping of reservoirs to non-negative quantities ***/
	   if( st[1] < 0 )
	      st[1] = 0 ;
	   if( st[2] < 0 )
	      st[2] = 0 ;
#endif

	  /*** output ***/
	   out = st[2] ; 

#endif

	  /*** output ***/
#ifdef _OVERFLOW_  
	   if( out > _MaxOutput_ ) 
	      fprintf( stderr, "Overflow error in Meddis Haircell output\%e\n", ( double ) out ) ;
#endif

	   *output = ( DataType ) out ;

	   output += bank->channels ;
	   input  += bank->channels ;
	}

    }

    return( npoints ) ;
 }


/************************** CloseHaircells() *********************************/

void CloseHaircells( bank_ptr )
    Pointer bank_ptr ;
{
    HaircellBank *bank = ( HaircellBank * ) bank_ptr ;
    HaircellInfo *cell_info ;
    int          channel ;

    for( channel = 0 ; channel < bank->channels ; channel++) {
       Delete( bank->cells[ channel ]->states ) ;
       Delete( bank->cells[ channel ] ) ;
    }
    
    Delete( bank->cells ) ;
    Delete( bank->coeffs ) ;
    Delete( bank->inputGain_ratio ) ; 
    Delete( bank ) ;

    return ;
}


/********** low level functions ************/

double update_inputGain( ratio ) 
  double ratio ;
{
  return( pow( 10., ( ratio - 0.5 ) * InputGain_max / 20. ) ) ; 
}