annotate wdf/wdf_ear.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

*/

/*
    ===========================================================
    wdf_ear.c 
    ===========================================================

    Wave digital filter (WDF) implementation of the outer and 
    middle ear (EAR) electroacoustic networks.

    Author        : Christian Giguere
    First written : 01st June, 1991
    Last edited   : 18th February, 1994
    
    References:
    (1) C.Giguere and P.C.Woodland (1992). Technical Report
	CUED/F-INFENG/TR.93, University of Cambridge.
    (2) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
    (3) M.E.Lutman and A.M.Martin (1979). J.Sound.Vib. 64(1): 133-157
    (4) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327.

    Note: Ref (1) describes an implementation where the outer 
	  ear and middle ear module is _UNCOUPLED_ to the cochlea. 
	  In this case, the middle ear network is that of Lutman 
	  and Martin (1979) unmodified.

	  Ref (2) describes an implementation where the outer
	  ear and middle ear module is _COUPLED_ to the cochlea.
	  In this case, the middle ear network is a modified
	  version of Lutman and Martin (1979) network (See Fig.2)
    ===========================================================
*/


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

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

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

#define   N_OF_FREEFIELD_STATES    (           12 )
#define   N_OF_EARTUBE_STATES      (           16 )      /* per segment */
#define   N_OF_EARMIDDLE_STATES    (           48 )  

#if 0
#define   _DEBUG_   
#define   _OVERFLOW_
#endif

#if 0
#define   _IMPULSE_     	   /* bypasses input wave with a delta impulse */
#endif

#if 0
#define   _EARDRUM_		   /* outputs eardrum pressure instead of stapes velocity */
#endif

/* Middle ear circuit elements in CGS units (From Lutman and Martin (1979) */

#define   La                       (     14.0e-03 )      /* Middle ear cavities */
#define   Cp                       (      5.1e-06 )     
#define   Ra                       (         10.0 )    
#define   Rm                       (        390.0 )    
#define   Ct                       (     0.35e-06 )    
#define   Rd1                      (        200.0 )      /* Eardrum losses */
#define   Cd1                      (      0.8e-06 )    
#define   Cd2                      (      0.4e-06 )    
#define   Rd2                      (        220.0 )    
#define   Ld                       (     15.0e-03 )    
#define   Rd3                      (       5900.0 )    
#define   Cd3                      (      0.2e-06 )    
#define   Lo                       (     40.0e-03 )      /* Eardrum, mallus, incus */
#define   Co                       (      1.4e-06 )    
#define   Ro                       (         70.0 )    
#define   Cs                       (     0.25e-06 )      /* Incudo-stapedial joint */
#define   Rs                       (       3000.0 )    
#define   Lc                       (     45.0e-03 )      /* Stapes and cochlea */
#define   Cc                       (      0.6e-06 )    
#define   Rc                       (        550.0 )    

/* Additional middle ear circuit elements in CGS units (From Giguere and Woodland (1994) */

#define   Ral                      (        100.0 )      /* annular ligaments */
#define   Kst_max                  (   1./0.1e-06 )      /* max. Stapedius stiffness */
#define   Kst_min_ratio            (       0.0001 )      /* min value of Kst_ratio */


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

double Kst_ratio = Kst_min_ratio ;                       /* initial fraction of Kst_max */
							 /* updated by feedback module */
static double rn_1 ; 


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

/******************************* WDF set-up functions ************************************
* name:                 function:
*
* FreefieldWDF()        Set-up function for the WDF implementation of the external ear network
*                       (Giguere and Woodland, 1994, Figs. 1 and 8). It generates the multiplier 
*                       coefficients and initializes the states of the filter. 
*
*                       In the case of the _UNCOUPLED_ ear version, it is called by the 
*                       function ``Ear()'' in file ``ear.c''.
*                       In the case of the _COUPLED_ ear version, it is called by the 
*                       function ``TLF_GenBank()'' in file ``bank_tl.c''.
*
*                       It returns a pointer to a structure containing all the relevant 
*                       information for the computation of the whole filter.
*
* EartubeWDF()          Set-up function for the WDF implementation of the concha and auditory
*                       canal transmission lines (Giguere and Woodland, 1994, Figs. 1 and 8). 
*                       It generates the multiplier coefficients and initializes the states 
*                       of the filter.
*
*                       In the case of the _UNCOUPLED_ ear version, it is called by the 
*                       function ``Ear()'' in file ``ear.c'' once for each segment of the 
*                       transmission lines starting from the outermost segment.
*                       In the case of the _COUPLED_ ear version, it is called by the 
*                       function ``TLF_GenBank()'' in file ``bank_tl.c'' once for each segment 
*                       of the transmission lines starting from the outermost segment.
*
*                       It returns a pointer to a structure containing all the relevant 
*                       information for the current filter segment.
*
* EarmiddleWDF()        Set-up function for the WDF implementation of the middle ear network
*                       (Giguere and Woodland, 1994, Figs. 2 and 9). It generates the multiplier 
*                       coefficients and initializes the states of the filter. 
*
*                       In the case of the _UNCOUPLED_ ear version, it is called by the 
*                       function ``Ear()'' in file ``ear.c''.
*                       In the case of the _COUPLED_ ear version, it is called by the 
*                       function ``TLF_GenBank()'' in file ``bank_tl.c''.
*
*                       It returns a pointer to a structure containing all the relevant 
*                       information for the computation of the whole filter.
******************************************************************************************/

/************************** FreefieldWDF() *****************************/

WaveWDFstate *FreefieldWDF( samplerate, rho, velocity, diffraction_radius, radiation_radius )
    double  samplerate ;
    double  rho, velocity ;
    double  diffraction_radius, radiation_radius ;
{
    DeclareNew( WaveWDFstate *, filter_state ) ;
    double  fs2 ;
    double  r1, r2, r3, g1, g2, g3 ;
    double  zh, zr ;
    double  Lh, Rh, Lr, Rr ;
 
   /*** electrical circuit elements (CGS) ***/
    Rh = rho * velocity / ( Pi * diffraction_radius * diffraction_radius ) ;
    Lh = 0.5 * rho / ( Pi * diffraction_radius ) ;
    Rr = rho * velocity / ( Pi * radiation_radius * radiation_radius ) ;
    Lr = 0.7 * rho / ( Pi * radiation_radius ) ;

   /*** sampling frequency ***/
    fs2 = samplerate * 2. ;
 
   /*** WDF-port resistances and multiplier coefficients ***/
 
    /* Adaptor 0 */
    g1 = 1. / ( fs2 * Lh ) ;                                    /* (1 / Zlh) */
    g2 = 1. / Rh ;                                              /* (1 / Zrh) */
    g3 = g1 + g2 ;
    zh = 1. / g3 ;
    filter_state->gamma01 = g1 / g3 ;

    /* Adaptor 2 */
    g1 = 1. / ( fs2 * Lr ) ;                                    /* (1 / Zlr) */
    g2 = 1. / Rr ;                                              /* (1 / Zrr) */
    g3 = g1 + g2 ;
    zr = 1. / g3 ;
    filter_state->gamma21 = g1 / g3 ;

    /* Adaptor 1 */
    r1 = zh ;
    r2 = zr ;
    r3 = rn_1 = r1 + r2 ;
    filter_state->gamma11 = r1 / r3 ;

#ifdef _DEBUG_
    fprintf( stderr, "Rh=%e  Lh=%e  Rr=%e  Lr=%e\n", Rh, Lh, Rr, Lr ) ;
    fprintf( stderr, "gamma01=%f\n", filter_state->gamma01 ) ;
    fprintf( stderr, "gamma11=%f\n", filter_state->gamma11 ) ;
    fprintf( stderr, "gamma21=%f\n", filter_state->gamma21 ) ;                                            
#endif

   /*** initialise filter states ***/
    filter_state->Nstates = N_OF_FREEFIELD_STATES ;
    filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
						      "wdf_ear.c for freefield states\n" ) ;

   /*** return ***/
    return( filter_state ) ;
}

/************************** EartubeWDF() *****************************/

EartubeWDFstate *EartubeWDF( samplerate, rho, velocity, diam, seglength, attn )
    double  samplerate ;
    double  rho, velocity ;
    double  diam, seglength ;
    double  attn ;
{
    DeclareNew( EartubeWDFstate *, filter_state ) ;
    double  an, d, c, k, alpha, fs2 ;
    double  Ln, Cn, Gn, Zn ;
    double  r1, r2, r3, g1, g2, g3, r33, g43, g53 ;

   /*** acoustic tube parameters for current segment ***/
    an = Pi * diam * diam / 4. ;                /* cross-sectional area of tube (cm2) */
    d = seglength;                              /* length of tube (cm) */ 
    c = velocity ;                              /* sound velocity in air (cm/s) */
    alpha = attn ;                              /* attenuation constant of travelling waves (1/cm) */

   /*** sampling frequency ***/
    fs2 = samplerate * 2. ;

   /*** electrical circuit elements (CGS units) for current segment ***/
    Ln = ( rho * d ) / an ;                     /* series inductor (either Lch or Lcl) */
    Cn = ( an * d ) / ( rho * c * c ) ;         /* shunt capacitor (either Cch or Ccl) */
    Zn = sqrt( Ln / Cn ) ;                      /* characteristic impedance ( either Zch or Zcl) */
    Gn = 2. / Zn * alpha * d ;                  /* shunt conductance (either Gch or Gcl) */

   /*** WDF multiplier coefficients for current tube segment ***/

    /* Adaptor 3 */
    r1 = fs2 * Ln / 2. ;                                        /* 0.5 Zl */
    r2 = rn_1 ;
    r3 = r33 = r1 + r2 ;
    filter_state->gamma31 = r1 / r3 ;

    /* Adaptor 5 */
    g1 = Gn ;                                                   /* (1 / Zg ) */
    g2 = fs2 * Cn ;                                             /* (1 / Zc ) */
    g3 = g53 = g1 + g2 ;
    filter_state->gamma51 = g1 / g3 ;

    /* Adaptor 4 */
    g1 = 1 / r33 ;
    g2 = g53 ;
    g3 = g43 = g1 + g2 ;
    filter_state->gamma41 = g1 / g3 ;

    /* Adaptor 6 */
    r1 = 1 / g43 ;
    r2 = fs2 * Ln / 2. ;                                        /* 0.5 Zl */
    r3 = rn_1 = r1 + r2 ;
    filter_state->gamma61 = r1 / r3 ;

#ifdef _DEBUG_
    fprintf( stderr, "gamma31=%f ", filter_state->gamma31 ) ;
    fprintf( stderr, "gamma41=%f ", filter_state->gamma41 ) ;
    fprintf( stderr, "gamma51=%f ", filter_state->gamma51 ) ;                                            
    fprintf( stderr, "gamma61=%f\n", filter_state->gamma61 ) ;  
#endif

   /*** initialise filter states ***/
    filter_state->Nstates = N_OF_EARTUBE_STATES ;
    filter_state->states =  ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
						       "wdf_ear.c for earcanal states\n" ) ;

   /*** return ***/
    return( filter_state ) ;
}


/********************** EarmiddleWDF() ***************************/

EarmiddleWDFstate *EarmiddleWDF( samplerate, zov, output_scale, ratio )
    double  samplerate, zov, output_scale, ratio ;
{
    DeclareNew( EarmiddleWDFstate *, filter_state ) ;
    double  fs2 ;
    double  r1, r2, r3, r4, g1, g2, g3, g4 ;
    double  r73, r84, r94, r103, r114, r123, r133, r143, r153 ;
    double  r163, r174, r183, r193, r204, r214 ;
    void    update_earmiddle_init(), update_earmiddle() ;

   /*** sampling frequency ***/
    fs2 = samplerate * 2. ;

   /*** WDF-port resistances and multiplier coefficients ***/
 
    /* Adaptor 9 */
    r1 = 1. / ( fs2 * Cp ) ;                                    /* Zcp */
    r2 = Ra ;                                                   /* Zra */
    r3 = fs2 * La ;                                             /* Zla */
    r4 = r94 = r1 + r2 + r3 ;
    filter_state->gamma91 = r1 / r4 ;
    filter_state->gamma92 = r2 / r4 ;

    /* Adaptor 8 */
    g1 = 1. / r94 ;
    g2 = 1. / Rm ;                                              /* (1 / Zrm) */
    g3 = fs2 * Ct ;                                             /* (1 / Zct) */
    g4 = g1 + g2 + g3 ;
    r84 = 1. /g4 ;
    filter_state->gamma81 = g1 / g4 ;
    filter_state->gamma82 = g2 / g4 ;

    /* Adaptor 7 */
    r1 = r84 ;
    r2 = rn_1 ;
    r3 = r73 = r1 + r2 ;
    filter_state->gamma71 = r1 / r3 ;

    /* Adaptor 13 */
    r1 = 1. / ( fs2 * Cd2 ) ;                                   /* Zcd2 */
    r2 = Rd2 ;                                                  /* Zrd2 */
    r3 = r133 = r1 + r2 ;
    filter_state->gamma131 = r1 / r3 ;

    /* Adaptor 12 */
    g1 = 1. / ( fs2 * Ld ) ;                                    /* (1 / Zld) */
    g2 = 1. / r133 ;
    g3 = g1 + g2 ;
    r123 = 1. / g3 ;
    filter_state->gamma121 = g1 / g3 ;

    /* Adaptor 14 */
    g1 = 1. / Rd3 ;                                             /* (1 / Zrd3) */
    g2 = fs2 * Cd3 ;                                            /* (1 / Zcd3) */
    g3 = g1 + g2 ;
    r143 = 1. / g3 ;
    filter_state->gamma141 = g1 / g3 ;

    /* Adaptor 15 */
    r1 = 1. / ( fs2 * Cd1 ) ;                                   /* Zcd1 */
    r2 = Rd1 ;                                                  /* Zrd1 */
    r3 = r153 = r1 + r2 ;
    filter_state->gamma151 = r1 / r3 ;

    /* Adaptor 11 */
    r1 = r143 ;
    r2 = r123 ;
    r3 = r153 ;
    r4 = r114 = r1 + r2 + r3 ;
    filter_state->gamma111 = r1 / r4 ;
    filter_state->gamma112 = r2 / r4 ;

    /* Adaptor 10 */
    g1 = 1. / r114 ;
    g2 = 1. / r73 ;
    g3 = g1 + g2 ;
    r103 = 1. /g3 ;
    filter_state->gamma101 = g1 / g3 ;
 
    /* Adaptor 17 */
    r1 = 1. / ( fs2 * Co ) ;                                    /* Zco */
    r2 = Ro ;                                                   /* Zro */
    r3 = fs2 * Lo ;                                             /* Zlo */
    r4 = r174 = r1 + r2 + r3 ;
    filter_state->gamma171 = r1 / r4 ;
    filter_state->gamma172 = r2 / r4 ;

    /* Adaptor 16 */
    r1 = r103 ;
    r2 = r174 ;
    r3 = r163 = r1 + r2 ;
    filter_state->gamma161 = r1 / r3 ;

    /* Adaptor 19 */
    r1 = 1. / ( fs2 * Cs ) ;                                    /* Zcs */
    r2 = Rs ;                                                   /* Zrs */
    r3 = r193 = r1 + r2 ;
    filter_state->gamma191 = r1 / r3 ;

    /* Adaptor 18 */
    g1 = 1. / r163 ;
    g2 = 1. / r193 ;
    g3 = g1 + g2 ;
    r183 = 1. / g3 ;
    filter_state->gamma181 = g1 / g3 ;

    /* Adaptor 20 */
#ifdef _DEBUG_
    fprintf( stderr, "Zov=%e\n", zov ) ;
#endif
    if( zov != 0. ) {                                           /* coupled version ? */
      r1 = Ral ;                                                /* Zral */
      r2 = zov / ( ratio * ratio ) ;                            /* Zov / (r*r) */
    }
    else {                                                      /* uncoupled version */
      r1 = 0.0 ;
      r2 = Rc ;                                                 /* Zrc */
    }
    r3 = r183 ;
    r4 = r204 = r1 + r2 +r3 ;
    filter_state->gamma201 = r1 / r4 ;
    filter_state->gamma202 = r2 / r4 ;
    filter_state->ratio = ratio ;

    /* Adaptor 21 */
    r1 = 1. / ( fs2 * Cc ) ;                                    /* Zcc */
    r2 = r204 ;
    if( zov != 0. )                                             /* coupled version ? */
      r3 = 0.0 ;
    else                                                        /* uncoupled version */
      r3 = fs2 * Lc ;                                           /* Zlc */
    r4 = r214 = r1 + r2 + r3 ;
    filter_state->gamma211 = r1 / r4 ;
    filter_state->gamma212 = r2 / r4 ;
    
    /* Adaptor 22 */
    r1 = Kst_ratio * Kst_max / fs2 ;                            /* Zcst (initial value only) */
    r2 = r214 ;
    r3 = r1 + r2 ;
    filter_state->gamma221 = 2. * r1 / r3 ;

#ifdef _DEBUG_
    fprintf ( stderr, "gamma71 =%f            \n", filter_state->gamma71  ) ;
    fprintf ( stderr, "gamma81 =%f gamma82 =%f\n", filter_state->gamma81
						 , filter_state->gamma82  ) ;
    fprintf ( stderr, "gamma91 =%f gamma92 =%f\n", filter_state->gamma91
						 , filter_state->gamma92  ) ;
    fprintf ( stderr, "gamma101=%f            \n", filter_state->gamma101 ) ;
    fprintf ( stderr, "gamma111=%f gamma112=%f\n", filter_state->gamma111
						 , filter_state->gamma112 ) ;
    fprintf ( stderr, "gamma121=%f            \n", filter_state->gamma121 ) ;
    fprintf ( stderr, "gamma131=%f            \n", filter_state->gamma131 ) ;
    fprintf ( stderr, "gamma141=%f            \n", filter_state->gamma141 ) ;
    fprintf ( stderr, "gamma151=%f            \n", filter_state->gamma151 ) ;
    fprintf ( stderr, "gamma161=%f            \n", filter_state->gamma161 ) ;
    fprintf ( stderr, "gamma171=%f gamma172=%f\n", filter_state->gamma171
						 , filter_state->gamma172 ) ;
    fprintf ( stderr, "gamma181=%f            \n", filter_state->gamma181 ) ;
    fprintf ( stderr, "gamma191=%f            \n", filter_state->gamma191 ) ;
    fprintf ( stderr, "gamma201=%f gamma202=%f\n", filter_state->gamma201
						 , filter_state->gamma202 ) ;
    fprintf ( stderr, "gamma211=%f gamma212=%f\n", filter_state->gamma211
						 , filter_state->gamma212 ) ;
    fprintf ( stderr, "gamma221=%f            \n", filter_state->gamma221 ) ;
#endif  

   /*** specify and initialise coefficient update procedure ***/
    filter_state->update_proc = update_earmiddle ; 
    ( void ) update_earmiddle_init( fs2, r214 ) ;

   /*** output scaling ***/ 
    filter_state->out_scale = output_scale ;

   /*** initialise filter states ***/
    filter_state->Nstates = N_OF_EARMIDDLE_STATES ;
    filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates, 
						      "wdf_ear.c for middle ear states\n" ) ;

   /*** return ***/
    return( filter_state ) ;
}


/*********************************** WDF procedures ***************************************
* name:                 function:
*
* DoEarWDF()            WDF realization of the _UNCOUPLED_ version of the outer ear and 
*                       middle ear networks. It is called by function ``ear_callback()'' 
*                       in file ``ear.c'' once for a specified number of input points. 
*                       It computes the stapes velocity for each input point and keeps 
*                       track of the filter states.  
*
*                       Note: The code for the realization of the _COUPLED_ ear version         
*                             is located in file ``wdf_tl.c''.
*
*
* CloseEarWDF()         Deletes all private data structures and arrays of the outer and 
*                       middle ear filter upon completion of filtering. 
*                       In the case of the _UNCOUPLED_ ear version, it is called by the 
*                       function ``ear_callback()'' in file ``ear.c''.
*                       In the case of the _COUPLED_ ear version, it is called by the 
*                       function ``tlf_bank_callbank()'' in file ``bank_tl.c''.
******************************************************************************************/

/******************************* DoEarWDF() **************************************/

void DoEarWDF( ws, eartube_states, ms, tube_segments, inout_ptr, points )
    register WaveWDFstate       *ws ;
    register EartubeWDFstate    **eartube_states ;
    register EarmiddleWDFstate  *ms ;
    register int                tube_segments ;
    register DataType           *inout_ptr ;
    register int                points ;
{
    register EartubeWDFstate    *ts ;
    register StateType          *st, in, out ;
    register StateType          a0, an_1, b1, b2, b3, b161, b202, b212, bn ;
    register int                tube_segment = -1 ;
#ifdef _IMPULSE_
    static   int                impulse = 10000 ;
#endif


    /***** update stapedial muscle state *****/

     ( void ) ms->update_proc( ms ) ;

    /***** for all points *****/

     while( points-- > 0 ) {

     /*** input point ***/

	in = ( StateType ) *inout_ptr ;

#ifdef _IMPULSE_
	if( impulse != 0 )  {
	   in = ( StateType ) impulse ;
	   impulse = 0 ;
	}
	else 
	   in = ( StateType ) 0 ;
#endif

     /*** computation in FORWARD direction for current input point ***/

     /* freefield - external ear */

	/* get states */
	st = ( StateType * ) ws->states ;      

	/* Adaptor 0 */
	st[4] = ws->gamma01 * ( st[3] - in ) ;                /* b00 */
	st[5] = st[4] + in + in ;                             /* b03 = -a11 */

	/* Adaptor 2 */
	st[1] = -ws->gamma21 * st[0] ;                        /* b20 */
	st[2] = st[1] ;                                       /* b23 = a12 */

	/* Adaptor 1 */
	st[6] = st[5] - st[2] ;                               /* b13 = a32 */
	an_1 = st[6] ;

     /* concha and auditory canal */

	while( ++tube_segment < tube_segments ) {

	   /* get states */
	   ts = *eartube_states++ ;
	   st = ( StateType * ) ts->states ;

	   /* Adaptor 3 */
	   st[3] = st[4] - an_1 ;                             /* b33 = a41 */

	   /* Adaptor 5 */
	   st[0] = -ts->gamma51 * st[1] ;                     /* b50 */
	   st[2] = st[0] + st[1] ;                            /* b53 = a42 */
 
	   /* Adaptor 4 */
	   st[5] = st[3] - st[2] ;                            /* -a42+a41 */
	   st[6] = ts->gamma41 * st[5] ;                      /* b40 */
	   st[7] = st[6] + st[2] ;                            /* b43 = a61 */

	   /* Adaptor 6 */
	   st[8] = st[9] - st[7] ;                            /* b63 = a32 */
	   an_1 = st[8] ;
	}

     /* middle ear */

	/* get states */
	st = ( StateType * ) ms->states ;      

	/* Adaptor 9 */
	st[31] = st[30] - st[29] ;                           /* b94 = a81 */

	/* Adaptor 8 */
	st[32] = st[35] - st[31] ;                           /* a83-a81 */
	st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ;   /* b80 */ 
	st[34] = st[33] + st[35] ;                           /* b84 = a71 */

	/* Adaptor 7 */
	st[36] = -st[34] - an_1 ;                            /* b73 = a102 */

	/* Adaptor 13 */
	st[15] = -st[14] ;                                   /* b133 = a122 */

	/* Adaptor 12 */
	st[16] = st[15] + st[18] ;                           /* a122-a121 */
	st[17] = -ms->gamma121 * st[16] ;                    /* b120 */
	st[19] = st[17] + st[15] ;                           /* b123 = a112 */

	/* Adaptor 14 */
	st[21] = -ms->gamma141 * st[20] ;                    /* b140 */
	st[22] = st[21] + st[20] ;                           /* b143 = a111 */

	/* Adaptor 15 */
	st[23] = -st[24] ;                                   /* b153 = a113 */

	/* Adaptor 11 */
	st[25] = -st[22] - st[19] - st[23] ;                 /* b114 = a101 */

	/* Adaptor 10 */
	st[26] = st[25] - st[36] ;                           /* -a102+a101 */
	st[27] = ms->gamma101 * st[26] ;                     /* b100 */
	st[28] = st[27] + st[36] ;                           /* b103 = a161 */

	/* Adaptor 17 */
	st[10] = st[12] - st[11] ;                           /* b174 = a162 */
	
	/* Adaptor 16 */
	st[13] = -st[28] - st[10] ;                          /* b163 = a181 */

	/* Adaptor 19 */
	st[5] = -st[6] ;                                     /* b193 = a182 */

	/* Adaptor 18 */
	st[7] = st[13] - st[5] ;                             /* -a182+a181 */
	st[8] = ms->gamma181 * st[7] ;                       /* b180 */
	st[9] = st[8] + st[5] ;                              /* b183 = a203 */

	/* Adaptor 20 */
	st[3] = -st[9] ;                                     /* b204 = a212 */

	/* Adaptor 21 */
	st[0] = st[2] - st[1] - st[3] ;                      /* b214 = a222 */


     /*** computation in BACKWARD direction for current input point ***/

     /* middle ear */

	/* Adaptor 22 */
	a0 = st[4] + st[0] ;                                  /* a220 */
	st[4] = st[4] - ms->gamma221 * a0 ;                   /* b221 */
	b2 = -st[4] - a0 ;                                    /* b222 = a214 */

	/* Adaptor 21 */
	a0 = b2 - st[0] ;                                     /* a210 */
	st[1] = st[1] - ms->gamma211 * a0 ;                   /* b211 */
	b212 = st[3] - ms->gamma212 * a0 ;                    /* b212 = a204 */
	st[2] = -st[1] - b212 - b2 ;                          /* b213 */

	/* Adaptor 20 */
	a0 = b212 - st[3] ;                                   /* a200 */
	b202 = - ms->gamma202 * a0 ;                          /* b202 */
	b3 = ms->gamma201 * a0 - b202 - b212 ;                /* b203 = a183 */

	/* Adaptor 18 */
	b2 = st[8] + b3 ;                                     /* b182 = a193 */
	b1 = b2 - st[7] ;                                     /* b181 = a163 */

	/* Adaptor 19 */
	st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ;       /* b191 */

	/* Adaptor 16 */
	b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ;      /* b161 = a103 */
	b2 = -b161 - b1 ;                                     /* b162 = a174 */

	/* Adaptor 17 */
	a0 = b2 - st[10] ;                                    /* a170 */
	st[11] = st[11] - ms->gamma171 * a0 ;                 /* b171 */
	st[12] = ms->gamma172 * a0 - b2 - st[11] ;            /* b173 */
 
	/* Adaptor 10 */
	st[37] = st[27] + b161 ;                              /* b102 = a73 */
	st[38] = st[37] - st[26] ;                            /* b101 = a114 */

	/* Adaptor 11 */
	a0 = st[38] - st[25] ;                                /* a110 */
	b1 = st[22] - ms->gamma111 * a0 ;                     /* b111 = a143 */
	b2 = st[19] - ms->gamma112 * a0 ;                     /* b112 = a123 */
	b3 = -b1 -b2 - st[38] ;                               /* b113 = a153 */

	/* Adaptor 15 */
	st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ;    /* b151 */

	/* Adaptor 14 */
	st[20] = st[21] + b1 ;                                /* b142 */

	/* Adaptor 12 */
	b2 = st[17] + b2 ;                                    /* b122 = a133 */
	st[18] = b2 + st[16] ;                                /* b121 */

	/* Adaptor 13 */
	st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ;    /* b131 */

	/* Adaptor 7 */
	b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ;     /* b71 = a84 */
	st[39] = -b1 - st[37] ;                               /* b72 = a63 */

	/* Adaptor 8 */
	st[35] = st[33] + b1 ;                                /* b83 */
	b1 = st[35] + st[32] ;                                /* b81 = a94 */

	/* Adaptor 9 */
	a0 = b1 - st[31] ;                                    /* a90 */
	st[29] = st[29] - ms->gamma91 * a0 ;                  /* b91 */
	st[30] = -st[29] + ms->gamma92 * a0 - b1 ;            /* b93 */

     bn = st[39] ;

     /* concha and auditory canal */

	while( tube_segment-- > 0 ) {

	   /* get states */
	   ts = *--eartube_states ;
	   st = ( StateType * ) ts->states ;

	   /* Adaptor 6 */
	   b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ;        /* b61 = a43 */
	   st[9] = -b1 - bn ;                                 /* b62 */

	   /* Adaptor 4 */
	   b2 = st[6] + b1 ;                                  /* b42 = a53 */
	   b1 = b2 - st[5] ;                                  /* b41 = a33 */

	   /* Adaptor 5 */
	   st[1] = st[0] + b2 ;                               /* b52 */

	   /* Adaptor 3 */
	   st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ;     /* b31 */
	   bn = -st[4] - b1 ;                                 /* b32 = a63 */
	}

     /* freefield - external ear */

	/* get states */
	st = ( StateType * ) ws->states ;      

	/* Adaptor 1 */
	b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ;           /* b11 = -a03 */
	b2 = -b1 - bn ;                                       /* b12 = a23 */

	/* Adaptor 2 */
	st[0] = st[1] + b2 + st[0] ;                          /* b61 */

	/* Adaptor 0 */
	st[3] = b1 - st[4] + st[3] ;                          /* -b01 + p */

     /*** output ***/
#ifdef _EARDRUM_
	st = ( StateType * ) ms->states ;      
	out = 0.5 * ( st[39] + an_1 ) * ms->out_scale ; 
#else
	out = -0.5 * b202 * ms->out_scale ; 
#endif

#ifdef _OVERFLOW_
        if( out > _MaxOutput_ || out < _MinOutput_ )
	   fprintf( stderr, "Overflow error in Outer/Middle ear output=%e\n", ( double ) out ) ;
#endif     
	*inout_ptr++ = out ; 

    }
    return ;
}


/******************************* CloseEarWDF() *************************************/

void CloseEarWDF( wave_states, eartube_states, earmiddle_states, tube_segments )
  register WaveWDFstate       *wave_states ;
  register EartubeWDFstate    **eartube_states ;
  register EarmiddleWDFstate  *earmiddle_states ;
  register int                tube_segments ;
{
  Delete( earmiddle_states->states ) ;
  Delete( earmiddle_states ) ;

  for( ; --tube_segments < 0 ; ) {
     Delete( eartube_states[ tube_segments ]->states ) ;
     Delete( eartube_states[ tube_segments ] ) ;
  }
  Delete( eartube_states ) ; 

  Delete( wave_states->states ) ;
  Delete( wave_states ) ;

  return ;

}

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

static double den221 ;
static double Kst_old_ratio ; 

void update_earmiddle_init( rate, z )
  double rate, z ;
{
  Kst_old_ratio = Kst_ratio ; 
  den221 = rate * z / Kst_max ;
  return ;
}

void update_earmiddle( ms )
  EarmiddleWDFstate *ms ;
{
  if( Kst_ratio < Kst_min_ratio )
    Kst_ratio = Kst_min_ratio ;

  ms->gamma221 = 2. * Kst_ratio / ( Kst_ratio + den221 ) ;
  ms->states[4] = ms->states[4] * Kst_ratio / Kst_old_ratio ;

  Kst_old_ratio = Kst_ratio ;
 
  return ;
}