comparison 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. ) ) ; 
+}
+