tomwalters@0: /*
tomwalters@0:     Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
tomwalters@0:     ===========================================================================
tomwalters@0: 
tomwalters@0:     Permission to use, copy, modify, and distribute this software without fee 
tomwalters@0:     is hereby granted for research purposes, provided that this copyright 
tomwalters@0:     notice appears in all copies and in all supporting documentation, and that 
tomwalters@0:     the software is not redistributed for any fee (except for a nominal shipping 
tomwalters@0:     charge). Anyone wanting to incorporate all or part of this software in a
tomwalters@0:     commercial product must obtain a license from the Medical Research Council.
tomwalters@0: 
tomwalters@0:     The MRC makes no representations about the suitability of this 
tomwalters@0:     software for any purpose.  It is provided "as is" without express or implied 
tomwalters@0:     warranty.
tomwalters@0:  
tomwalters@0:     THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
tomwalters@0:     ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
tomwalters@0:     A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY 
tomwalters@0:     DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
tomwalters@0:     AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 
tomwalters@0:     OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*  
tomwalters@0:     Acknowledgment:
tomwalters@0:     ==============
tomwalters@0: 
tomwalters@0:     The source code provided in this file was originally developed by 
tomwalters@0:     Christian Giguere as part of a Ph.D degree at the Department of
tomwalters@0:     Engineering of the University of Cambridge from April 1990 to
tomwalters@0:     November 1993. The code was subsequently adapted under a grant
tomwalters@0:     from the Hearing Research Trust for full compatibility with 
tomwalters@0:     AIM Release 6.15.
tomwalters@0: 
tomwalters@0:     Christian Giguere 25/03/94
tomwalters@0: 
tomwalters@0: */
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:     ===========================================================
tomwalters@0:     wdf_ear.c 
tomwalters@0:     ===========================================================
tomwalters@0: 
tomwalters@0:     Wave digital filter (WDF) implementation of the outer and 
tomwalters@0:     middle ear (EAR) electroacoustic networks.
tomwalters@0: 
tomwalters@0:     Author        : Christian Giguere
tomwalters@0:     First written : 01st June, 1991
tomwalters@0:     Last edited   : 18th February, 1994
tomwalters@0:     
tomwalters@0:     References:
tomwalters@0:     (1) C.Giguere and P.C.Woodland (1992). Technical Report
tomwalters@0: 	CUED/F-INFENG/TR.93, University of Cambridge.
tomwalters@0:     (2) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
tomwalters@0:     (3) M.E.Lutman and A.M.Martin (1979). J.Sound.Vib. 64(1): 133-157
tomwalters@0:     (4) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327.
tomwalters@0: 
tomwalters@0:     Note: Ref (1) describes an implementation where the outer 
tomwalters@0: 	  ear and middle ear module is _UNCOUPLED_ to the cochlea. 
tomwalters@0: 	  In this case, the middle ear network is that of Lutman 
tomwalters@0: 	  and Martin (1979) unmodified.
tomwalters@0: 
tomwalters@0: 	  Ref (2) describes an implementation where the outer
tomwalters@0: 	  ear and middle ear module is _COUPLED_ to the cochlea.
tomwalters@0: 	  In this case, the middle ear network is a modified
tomwalters@0: 	  version of Lutman and Martin (1979) network (See Fig.2)
tomwalters@0:     ===========================================================
tomwalters@0: */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** includes *****/
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: #include <stdio.h>
tomwalters@0: #include "stitch.h"
tomwalters@0: #include "calc.h"
tomwalters@0: #include "calc_tl.h"
tomwalters@0: #include "wdf_ear.h"
tomwalters@0: 
tomwalters@0: /***** defines *****/
tomwalters@0: 
tomwalters@0: #define   N_OF_FREEFIELD_STATES    (           12 )
tomwalters@0: #define   N_OF_EARTUBE_STATES      (           16 )      /* per segment */
tomwalters@0: #define   N_OF_EARMIDDLE_STATES    (           48 )  
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _DEBUG_   
tomwalters@0: #define   _OVERFLOW_
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _IMPULSE_     	   /* bypasses input wave with a delta impulse */
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #if 0
tomwalters@0: #define   _EARDRUM_		   /* outputs eardrum pressure instead of stapes velocity */
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: /* Middle ear circuit elements in CGS units (From Lutman and Martin (1979) */
tomwalters@0: 
tomwalters@0: #define   La                       (     14.0e-03 )      /* Middle ear cavities */
tomwalters@0: #define   Cp                       (      5.1e-06 )     
tomwalters@0: #define   Ra                       (         10.0 )    
tomwalters@0: #define   Rm                       (        390.0 )    
tomwalters@0: #define   Ct                       (     0.35e-06 )    
tomwalters@0: #define   Rd1                      (        200.0 )      /* Eardrum losses */
tomwalters@0: #define   Cd1                      (      0.8e-06 )    
tomwalters@0: #define   Cd2                      (      0.4e-06 )    
tomwalters@0: #define   Rd2                      (        220.0 )    
tomwalters@0: #define   Ld                       (     15.0e-03 )    
tomwalters@0: #define   Rd3                      (       5900.0 )    
tomwalters@0: #define   Cd3                      (      0.2e-06 )    
tomwalters@0: #define   Lo                       (     40.0e-03 )      /* Eardrum, mallus, incus */
tomwalters@0: #define   Co                       (      1.4e-06 )    
tomwalters@0: #define   Ro                       (         70.0 )    
tomwalters@0: #define   Cs                       (     0.25e-06 )      /* Incudo-stapedial joint */
tomwalters@0: #define   Rs                       (       3000.0 )    
tomwalters@0: #define   Lc                       (     45.0e-03 )      /* Stapes and cochlea */
tomwalters@0: #define   Cc                       (      0.6e-06 )    
tomwalters@0: #define   Rc                       (        550.0 )    
tomwalters@0: 
tomwalters@0: /* Additional middle ear circuit elements in CGS units (From Giguere and Woodland (1994) */
tomwalters@0: 
tomwalters@0: #define   Ral                      (        100.0 )      /* annular ligaments */
tomwalters@0: #define   Kst_max                  (   1./0.1e-06 )      /* max. Stapedius stiffness */
tomwalters@0: #define   Kst_min_ratio            (       0.0001 )      /* min value of Kst_ratio */
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** external variables *****/
tomwalters@0: 
tomwalters@0: double Kst_ratio = Kst_min_ratio ;                       /* initial fraction of Kst_max */
tomwalters@0: 							 /* updated by feedback module */
tomwalters@0: static double rn_1 ; 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***** functions *****/
tomwalters@0: 
tomwalters@0: /******************************* WDF set-up functions ************************************
tomwalters@0: * name:                 function:
tomwalters@0: *
tomwalters@0: * FreefieldWDF()        Set-up function for the WDF implementation of the external ear network
tomwalters@0: *                       (Giguere and Woodland, 1994, Figs. 1 and 8). It generates the multiplier 
tomwalters@0: *                       coefficients and initializes the states of the filter. 
tomwalters@0: *
tomwalters@0: *                       In the case of the _UNCOUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``Ear()'' in file ``ear.c''.
tomwalters@0: *                       In the case of the _COUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``TLF_GenBank()'' in file ``bank_tl.c''.
tomwalters@0: *
tomwalters@0: *                       It returns a pointer to a structure containing all the relevant 
tomwalters@0: *                       information for the computation of the whole filter.
tomwalters@0: *
tomwalters@0: * EartubeWDF()          Set-up function for the WDF implementation of the concha and auditory
tomwalters@0: *                       canal transmission lines (Giguere and Woodland, 1994, Figs. 1 and 8). 
tomwalters@0: *                       It generates the multiplier coefficients and initializes the states 
tomwalters@0: *                       of the filter.
tomwalters@0: *
tomwalters@0: *                       In the case of the _UNCOUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``Ear()'' in file ``ear.c'' once for each segment of the 
tomwalters@0: *                       transmission lines starting from the outermost segment.
tomwalters@0: *                       In the case of the _COUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``TLF_GenBank()'' in file ``bank_tl.c'' once for each segment 
tomwalters@0: *                       of the transmission lines starting from the outermost segment.
tomwalters@0: *
tomwalters@0: *                       It returns a pointer to a structure containing all the relevant 
tomwalters@0: *                       information for the current filter segment.
tomwalters@0: *
tomwalters@0: * EarmiddleWDF()        Set-up function for the WDF implementation of the middle ear network
tomwalters@0: *                       (Giguere and Woodland, 1994, Figs. 2 and 9). It generates the multiplier 
tomwalters@0: *                       coefficients and initializes the states of the filter. 
tomwalters@0: *
tomwalters@0: *                       In the case of the _UNCOUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``Ear()'' in file ``ear.c''.
tomwalters@0: *                       In the case of the _COUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``TLF_GenBank()'' in file ``bank_tl.c''.
tomwalters@0: *
tomwalters@0: *                       It returns a pointer to a structure containing all the relevant 
tomwalters@0: *                       information for the computation of the whole filter.
tomwalters@0: ******************************************************************************************/
tomwalters@0: 
tomwalters@0: /************************** FreefieldWDF() *****************************/
tomwalters@0: 
tomwalters@0: WaveWDFstate *FreefieldWDF( samplerate, rho, velocity, diffraction_radius, radiation_radius )
tomwalters@0:     double  samplerate ;
tomwalters@0:     double  rho, velocity ;
tomwalters@0:     double  diffraction_radius, radiation_radius ;
tomwalters@0: {
tomwalters@0:     DeclareNew( WaveWDFstate *, filter_state ) ;
tomwalters@0:     double  fs2 ;
tomwalters@0:     double  r1, r2, r3, g1, g2, g3 ;
tomwalters@0:     double  zh, zr ;
tomwalters@0:     double  Lh, Rh, Lr, Rr ;
tomwalters@0:  
tomwalters@0:    /*** electrical circuit elements (CGS) ***/
tomwalters@0:     Rh = rho * velocity / ( Pi * diffraction_radius * diffraction_radius ) ;
tomwalters@0:     Lh = 0.5 * rho / ( Pi * diffraction_radius ) ;
tomwalters@0:     Rr = rho * velocity / ( Pi * radiation_radius * radiation_radius ) ;
tomwalters@0:     Lr = 0.7 * rho / ( Pi * radiation_radius ) ;
tomwalters@0: 
tomwalters@0:    /*** sampling frequency ***/
tomwalters@0:     fs2 = samplerate * 2. ;
tomwalters@0:  
tomwalters@0:    /*** WDF-port resistances and multiplier coefficients ***/
tomwalters@0:  
tomwalters@0:     /* Adaptor 0 */
tomwalters@0:     g1 = 1. / ( fs2 * Lh ) ;                                    /* (1 / Zlh) */
tomwalters@0:     g2 = 1. / Rh ;                                              /* (1 / Zrh) */
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     zh = 1. / g3 ;
tomwalters@0:     filter_state->gamma01 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 2 */
tomwalters@0:     g1 = 1. / ( fs2 * Lr ) ;                                    /* (1 / Zlr) */
tomwalters@0:     g2 = 1. / Rr ;                                              /* (1 / Zrr) */
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     zr = 1. / g3 ;
tomwalters@0:     filter_state->gamma21 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 1 */
tomwalters@0:     r1 = zh ;
tomwalters@0:     r2 = zr ;
tomwalters@0:     r3 = rn_1 = r1 + r2 ;
tomwalters@0:     filter_state->gamma11 = r1 / r3 ;
tomwalters@0: 
tomwalters@0: #ifdef _DEBUG_
tomwalters@0:     fprintf( stderr, "Rh=%e  Lh=%e  Rr=%e  Lr=%e\n", Rh, Lh, Rr, Lr ) ;
tomwalters@0:     fprintf( stderr, "gamma01=%f\n", filter_state->gamma01 ) ;
tomwalters@0:     fprintf( stderr, "gamma11=%f\n", filter_state->gamma11 ) ;
tomwalters@0:     fprintf( stderr, "gamma21=%f\n", filter_state->gamma21 ) ;                                            
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:    /*** initialise filter states ***/
tomwalters@0:     filter_state->Nstates = N_OF_FREEFIELD_STATES ;
tomwalters@0:     filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
tomwalters@0: 						      "wdf_ear.c for freefield states\n" ) ;
tomwalters@0: 
tomwalters@0:    /*** return ***/
tomwalters@0:     return( filter_state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /************************** EartubeWDF() *****************************/
tomwalters@0: 
tomwalters@0: EartubeWDFstate *EartubeWDF( samplerate, rho, velocity, diam, seglength, attn )
tomwalters@0:     double  samplerate ;
tomwalters@0:     double  rho, velocity ;
tomwalters@0:     double  diam, seglength ;
tomwalters@0:     double  attn ;
tomwalters@0: {
tomwalters@0:     DeclareNew( EartubeWDFstate *, filter_state ) ;
tomwalters@0:     double  an, d, c, k, alpha, fs2 ;
tomwalters@0:     double  Ln, Cn, Gn, Zn ;
tomwalters@0:     double  r1, r2, r3, g1, g2, g3, r33, g43, g53 ;
tomwalters@0: 
tomwalters@0:    /*** acoustic tube parameters for current segment ***/
tomwalters@0:     an = Pi * diam * diam / 4. ;                /* cross-sectional area of tube (cm2) */
tomwalters@0:     d = seglength;                              /* length of tube (cm) */ 
tomwalters@0:     c = velocity ;                              /* sound velocity in air (cm/s) */
tomwalters@0:     alpha = attn ;                              /* attenuation constant of travelling waves (1/cm) */
tomwalters@0: 
tomwalters@0:    /*** sampling frequency ***/
tomwalters@0:     fs2 = samplerate * 2. ;
tomwalters@0: 
tomwalters@0:    /*** electrical circuit elements (CGS units) for current segment ***/
tomwalters@0:     Ln = ( rho * d ) / an ;                     /* series inductor (either Lch or Lcl) */
tomwalters@0:     Cn = ( an * d ) / ( rho * c * c ) ;         /* shunt capacitor (either Cch or Ccl) */
tomwalters@0:     Zn = sqrt( Ln / Cn ) ;                      /* characteristic impedance ( either Zch or Zcl) */
tomwalters@0:     Gn = 2. / Zn * alpha * d ;                  /* shunt conductance (either Gch or Gcl) */
tomwalters@0: 
tomwalters@0:    /*** WDF multiplier coefficients for current tube segment ***/
tomwalters@0: 
tomwalters@0:     /* Adaptor 3 */
tomwalters@0:     r1 = fs2 * Ln / 2. ;                                        /* 0.5 Zl */
tomwalters@0:     r2 = rn_1 ;
tomwalters@0:     r3 = r33 = r1 + r2 ;
tomwalters@0:     filter_state->gamma31 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 5 */
tomwalters@0:     g1 = Gn ;                                                   /* (1 / Zg ) */
tomwalters@0:     g2 = fs2 * Cn ;                                             /* (1 / Zc ) */
tomwalters@0:     g3 = g53 = g1 + g2 ;
tomwalters@0:     filter_state->gamma51 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 4 */
tomwalters@0:     g1 = 1 / r33 ;
tomwalters@0:     g2 = g53 ;
tomwalters@0:     g3 = g43 = g1 + g2 ;
tomwalters@0:     filter_state->gamma41 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 6 */
tomwalters@0:     r1 = 1 / g43 ;
tomwalters@0:     r2 = fs2 * Ln / 2. ;                                        /* 0.5 Zl */
tomwalters@0:     r3 = rn_1 = r1 + r2 ;
tomwalters@0:     filter_state->gamma61 = r1 / r3 ;
tomwalters@0: 
tomwalters@0: #ifdef _DEBUG_
tomwalters@0:     fprintf( stderr, "gamma31=%f ", filter_state->gamma31 ) ;
tomwalters@0:     fprintf( stderr, "gamma41=%f ", filter_state->gamma41 ) ;
tomwalters@0:     fprintf( stderr, "gamma51=%f ", filter_state->gamma51 ) ;                                            
tomwalters@0:     fprintf( stderr, "gamma61=%f\n", filter_state->gamma61 ) ;  
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:    /*** initialise filter states ***/
tomwalters@0:     filter_state->Nstates = N_OF_EARTUBE_STATES ;
tomwalters@0:     filter_state->states =  ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
tomwalters@0: 						       "wdf_ear.c for earcanal states\n" ) ;
tomwalters@0: 
tomwalters@0:    /*** return ***/
tomwalters@0:     return( filter_state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /********************** EarmiddleWDF() ***************************/
tomwalters@0: 
tomwalters@0: EarmiddleWDFstate *EarmiddleWDF( samplerate, zov, output_scale, ratio )
tomwalters@0:     double  samplerate, zov, output_scale, ratio ;
tomwalters@0: {
tomwalters@0:     DeclareNew( EarmiddleWDFstate *, filter_state ) ;
tomwalters@0:     double  fs2 ;
tomwalters@0:     double  r1, r2, r3, r4, g1, g2, g3, g4 ;
tomwalters@0:     double  r73, r84, r94, r103, r114, r123, r133, r143, r153 ;
tomwalters@0:     double  r163, r174, r183, r193, r204, r214 ;
tomwalters@0:     void    update_earmiddle_init(), update_earmiddle() ;
tomwalters@0: 
tomwalters@0:    /*** sampling frequency ***/
tomwalters@0:     fs2 = samplerate * 2. ;
tomwalters@0: 
tomwalters@0:    /*** WDF-port resistances and multiplier coefficients ***/
tomwalters@0:  
tomwalters@0:     /* Adaptor 9 */
tomwalters@0:     r1 = 1. / ( fs2 * Cp ) ;                                    /* Zcp */
tomwalters@0:     r2 = Ra ;                                                   /* Zra */
tomwalters@0:     r3 = fs2 * La ;                                             /* Zla */
tomwalters@0:     r4 = r94 = r1 + r2 + r3 ;
tomwalters@0:     filter_state->gamma91 = r1 / r4 ;
tomwalters@0:     filter_state->gamma92 = r2 / r4 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 8 */
tomwalters@0:     g1 = 1. / r94 ;
tomwalters@0:     g2 = 1. / Rm ;                                              /* (1 / Zrm) */
tomwalters@0:     g3 = fs2 * Ct ;                                             /* (1 / Zct) */
tomwalters@0:     g4 = g1 + g2 + g3 ;
tomwalters@0:     r84 = 1. /g4 ;
tomwalters@0:     filter_state->gamma81 = g1 / g4 ;
tomwalters@0:     filter_state->gamma82 = g2 / g4 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 7 */
tomwalters@0:     r1 = r84 ;
tomwalters@0:     r2 = rn_1 ;
tomwalters@0:     r3 = r73 = r1 + r2 ;
tomwalters@0:     filter_state->gamma71 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 13 */
tomwalters@0:     r1 = 1. / ( fs2 * Cd2 ) ;                                   /* Zcd2 */
tomwalters@0:     r2 = Rd2 ;                                                  /* Zrd2 */
tomwalters@0:     r3 = r133 = r1 + r2 ;
tomwalters@0:     filter_state->gamma131 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 12 */
tomwalters@0:     g1 = 1. / ( fs2 * Ld ) ;                                    /* (1 / Zld) */
tomwalters@0:     g2 = 1. / r133 ;
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     r123 = 1. / g3 ;
tomwalters@0:     filter_state->gamma121 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 14 */
tomwalters@0:     g1 = 1. / Rd3 ;                                             /* (1 / Zrd3) */
tomwalters@0:     g2 = fs2 * Cd3 ;                                            /* (1 / Zcd3) */
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     r143 = 1. / g3 ;
tomwalters@0:     filter_state->gamma141 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 15 */
tomwalters@0:     r1 = 1. / ( fs2 * Cd1 ) ;                                   /* Zcd1 */
tomwalters@0:     r2 = Rd1 ;                                                  /* Zrd1 */
tomwalters@0:     r3 = r153 = r1 + r2 ;
tomwalters@0:     filter_state->gamma151 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 11 */
tomwalters@0:     r1 = r143 ;
tomwalters@0:     r2 = r123 ;
tomwalters@0:     r3 = r153 ;
tomwalters@0:     r4 = r114 = r1 + r2 + r3 ;
tomwalters@0:     filter_state->gamma111 = r1 / r4 ;
tomwalters@0:     filter_state->gamma112 = r2 / r4 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 10 */
tomwalters@0:     g1 = 1. / r114 ;
tomwalters@0:     g2 = 1. / r73 ;
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     r103 = 1. /g3 ;
tomwalters@0:     filter_state->gamma101 = g1 / g3 ;
tomwalters@0:  
tomwalters@0:     /* Adaptor 17 */
tomwalters@0:     r1 = 1. / ( fs2 * Co ) ;                                    /* Zco */
tomwalters@0:     r2 = Ro ;                                                   /* Zro */
tomwalters@0:     r3 = fs2 * Lo ;                                             /* Zlo */
tomwalters@0:     r4 = r174 = r1 + r2 + r3 ;
tomwalters@0:     filter_state->gamma171 = r1 / r4 ;
tomwalters@0:     filter_state->gamma172 = r2 / r4 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 16 */
tomwalters@0:     r1 = r103 ;
tomwalters@0:     r2 = r174 ;
tomwalters@0:     r3 = r163 = r1 + r2 ;
tomwalters@0:     filter_state->gamma161 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 19 */
tomwalters@0:     r1 = 1. / ( fs2 * Cs ) ;                                    /* Zcs */
tomwalters@0:     r2 = Rs ;                                                   /* Zrs */
tomwalters@0:     r3 = r193 = r1 + r2 ;
tomwalters@0:     filter_state->gamma191 = r1 / r3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 18 */
tomwalters@0:     g1 = 1. / r163 ;
tomwalters@0:     g2 = 1. / r193 ;
tomwalters@0:     g3 = g1 + g2 ;
tomwalters@0:     r183 = 1. / g3 ;
tomwalters@0:     filter_state->gamma181 = g1 / g3 ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 20 */
tomwalters@0: #ifdef _DEBUG_
tomwalters@0:     fprintf( stderr, "Zov=%e\n", zov ) ;
tomwalters@0: #endif
tomwalters@0:     if( zov != 0. ) {                                           /* coupled version ? */
tomwalters@0:       r1 = Ral ;                                                /* Zral */
tomwalters@0:       r2 = zov / ( ratio * ratio ) ;                            /* Zov / (r*r) */
tomwalters@0:     }
tomwalters@0:     else {                                                      /* uncoupled version */
tomwalters@0:       r1 = 0.0 ;
tomwalters@0:       r2 = Rc ;                                                 /* Zrc */
tomwalters@0:     }
tomwalters@0:     r3 = r183 ;
tomwalters@0:     r4 = r204 = r1 + r2 +r3 ;
tomwalters@0:     filter_state->gamma201 = r1 / r4 ;
tomwalters@0:     filter_state->gamma202 = r2 / r4 ;
tomwalters@0:     filter_state->ratio = ratio ;
tomwalters@0: 
tomwalters@0:     /* Adaptor 21 */
tomwalters@0:     r1 = 1. / ( fs2 * Cc ) ;                                    /* Zcc */
tomwalters@0:     r2 = r204 ;
tomwalters@0:     if( zov != 0. )                                             /* coupled version ? */
tomwalters@0:       r3 = 0.0 ;
tomwalters@0:     else                                                        /* uncoupled version */
tomwalters@0:       r3 = fs2 * Lc ;                                           /* Zlc */
tomwalters@0:     r4 = r214 = r1 + r2 + r3 ;
tomwalters@0:     filter_state->gamma211 = r1 / r4 ;
tomwalters@0:     filter_state->gamma212 = r2 / r4 ;
tomwalters@0:     
tomwalters@0:     /* Adaptor 22 */
tomwalters@0:     r1 = Kst_ratio * Kst_max / fs2 ;                            /* Zcst (initial value only) */
tomwalters@0:     r2 = r214 ;
tomwalters@0:     r3 = r1 + r2 ;
tomwalters@0:     filter_state->gamma221 = 2. * r1 / r3 ;
tomwalters@0: 
tomwalters@0: #ifdef _DEBUG_
tomwalters@0:     fprintf ( stderr, "gamma71 =%f            \n", filter_state->gamma71  ) ;
tomwalters@0:     fprintf ( stderr, "gamma81 =%f gamma82 =%f\n", filter_state->gamma81
tomwalters@0: 						 , filter_state->gamma82  ) ;
tomwalters@0:     fprintf ( stderr, "gamma91 =%f gamma92 =%f\n", filter_state->gamma91
tomwalters@0: 						 , filter_state->gamma92  ) ;
tomwalters@0:     fprintf ( stderr, "gamma101=%f            \n", filter_state->gamma101 ) ;
tomwalters@0:     fprintf ( stderr, "gamma111=%f gamma112=%f\n", filter_state->gamma111
tomwalters@0: 						 , filter_state->gamma112 ) ;
tomwalters@0:     fprintf ( stderr, "gamma121=%f            \n", filter_state->gamma121 ) ;
tomwalters@0:     fprintf ( stderr, "gamma131=%f            \n", filter_state->gamma131 ) ;
tomwalters@0:     fprintf ( stderr, "gamma141=%f            \n", filter_state->gamma141 ) ;
tomwalters@0:     fprintf ( stderr, "gamma151=%f            \n", filter_state->gamma151 ) ;
tomwalters@0:     fprintf ( stderr, "gamma161=%f            \n", filter_state->gamma161 ) ;
tomwalters@0:     fprintf ( stderr, "gamma171=%f gamma172=%f\n", filter_state->gamma171
tomwalters@0: 						 , filter_state->gamma172 ) ;
tomwalters@0:     fprintf ( stderr, "gamma181=%f            \n", filter_state->gamma181 ) ;
tomwalters@0:     fprintf ( stderr, "gamma191=%f            \n", filter_state->gamma191 ) ;
tomwalters@0:     fprintf ( stderr, "gamma201=%f gamma202=%f\n", filter_state->gamma201
tomwalters@0: 						 , filter_state->gamma202 ) ;
tomwalters@0:     fprintf ( stderr, "gamma211=%f gamma212=%f\n", filter_state->gamma211
tomwalters@0: 						 , filter_state->gamma212 ) ;
tomwalters@0:     fprintf ( stderr, "gamma221=%f            \n", filter_state->gamma221 ) ;
tomwalters@0: #endif  
tomwalters@0: 
tomwalters@0:    /*** specify and initialise coefficient update procedure ***/
tomwalters@0:     filter_state->update_proc = update_earmiddle ; 
tomwalters@0:     ( void ) update_earmiddle_init( fs2, r214 ) ;
tomwalters@0: 
tomwalters@0:    /*** output scaling ***/ 
tomwalters@0:     filter_state->out_scale = output_scale ;
tomwalters@0: 
tomwalters@0:    /*** initialise filter states ***/
tomwalters@0:     filter_state->Nstates = N_OF_EARMIDDLE_STATES ;
tomwalters@0:     filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates, 
tomwalters@0: 						      "wdf_ear.c for middle ear states\n" ) ;
tomwalters@0: 
tomwalters@0:    /*** return ***/
tomwalters@0:     return( filter_state ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*********************************** WDF procedures ***************************************
tomwalters@0: * name:                 function:
tomwalters@0: *
tomwalters@0: * DoEarWDF()            WDF realization of the _UNCOUPLED_ version of the outer ear and 
tomwalters@0: *                       middle ear networks. It is called by function ``ear_callback()'' 
tomwalters@0: *                       in file ``ear.c'' once for a specified number of input points. 
tomwalters@0: *                       It computes the stapes velocity for each input point and keeps 
tomwalters@0: *                       track of the filter states.  
tomwalters@0: *
tomwalters@0: *                       Note: The code for the realization of the _COUPLED_ ear version         
tomwalters@0: *                             is located in file ``wdf_tl.c''.
tomwalters@0: *
tomwalters@0: *
tomwalters@0: * CloseEarWDF()         Deletes all private data structures and arrays of the outer and 
tomwalters@0: *                       middle ear filter upon completion of filtering. 
tomwalters@0: *                       In the case of the _UNCOUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``ear_callback()'' in file ``ear.c''.
tomwalters@0: *                       In the case of the _COUPLED_ ear version, it is called by the 
tomwalters@0: *                       function ``tlf_bank_callbank()'' in file ``bank_tl.c''.
tomwalters@0: ******************************************************************************************/
tomwalters@0: 
tomwalters@0: /******************************* DoEarWDF() **************************************/
tomwalters@0: 
tomwalters@0: void DoEarWDF( ws, eartube_states, ms, tube_segments, inout_ptr, points )
tomwalters@0:     register WaveWDFstate       *ws ;
tomwalters@0:     register EartubeWDFstate    **eartube_states ;
tomwalters@0:     register EarmiddleWDFstate  *ms ;
tomwalters@0:     register int                tube_segments ;
tomwalters@0:     register DataType           *inout_ptr ;
tomwalters@0:     register int                points ;
tomwalters@0: {
tomwalters@0:     register EartubeWDFstate    *ts ;
tomwalters@0:     register StateType          *st, in, out ;
tomwalters@0:     register StateType          a0, an_1, b1, b2, b3, b161, b202, b212, bn ;
tomwalters@0:     register int                tube_segment = -1 ;
tomwalters@0: #ifdef _IMPULSE_
tomwalters@0:     static   int                impulse = 10000 ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /***** update stapedial muscle state *****/
tomwalters@0: 
tomwalters@0:      ( void ) ms->update_proc( ms ) ;
tomwalters@0: 
tomwalters@0:     /***** for all points *****/
tomwalters@0: 
tomwalters@0:      while( points-- > 0 ) {
tomwalters@0: 
tomwalters@0:      /*** input point ***/
tomwalters@0: 
tomwalters@0: 	in = ( StateType ) *inout_ptr ;
tomwalters@0: 
tomwalters@0: #ifdef _IMPULSE_
tomwalters@0: 	if( impulse != 0 )  {
tomwalters@0: 	   in = ( StateType ) impulse ;
tomwalters@0: 	   impulse = 0 ;
tomwalters@0: 	}
tomwalters@0: 	else 
tomwalters@0: 	   in = ( StateType ) 0 ;
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0:      /*** computation in FORWARD direction for current input point ***/
tomwalters@0: 
tomwalters@0:      /* freefield - external ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ws->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 0 */
tomwalters@0: 	st[4] = ws->gamma01 * ( st[3] - in ) ;                /* b00 */
tomwalters@0: 	st[5] = st[4] + in + in ;                             /* b03 = -a11 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 2 */
tomwalters@0: 	st[1] = -ws->gamma21 * st[0] ;                        /* b20 */
tomwalters@0: 	st[2] = st[1] ;                                       /* b23 = a12 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 1 */
tomwalters@0: 	st[6] = st[5] - st[2] ;                               /* b13 = a32 */
tomwalters@0: 	an_1 = st[6] ;
tomwalters@0: 
tomwalters@0:      /* concha and auditory canal */
tomwalters@0: 
tomwalters@0: 	while( ++tube_segment < tube_segments ) {
tomwalters@0: 
tomwalters@0: 	   /* get states */
tomwalters@0: 	   ts = *eartube_states++ ;
tomwalters@0: 	   st = ( StateType * ) ts->states ;
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 3 */
tomwalters@0: 	   st[3] = st[4] - an_1 ;                             /* b33 = a41 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 5 */
tomwalters@0: 	   st[0] = -ts->gamma51 * st[1] ;                     /* b50 */
tomwalters@0: 	   st[2] = st[0] + st[1] ;                            /* b53 = a42 */
tomwalters@0:  
tomwalters@0: 	   /* Adaptor 4 */
tomwalters@0: 	   st[5] = st[3] - st[2] ;                            /* -a42+a41 */
tomwalters@0: 	   st[6] = ts->gamma41 * st[5] ;                      /* b40 */
tomwalters@0: 	   st[7] = st[6] + st[2] ;                            /* b43 = a61 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 6 */
tomwalters@0: 	   st[8] = st[9] - st[7] ;                            /* b63 = a32 */
tomwalters@0: 	   an_1 = st[8] ;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:      /* middle ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ms->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 9 */
tomwalters@0: 	st[31] = st[30] - st[29] ;                           /* b94 = a81 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 8 */
tomwalters@0: 	st[32] = st[35] - st[31] ;                           /* a83-a81 */
tomwalters@0: 	st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ;   /* b80 */ 
tomwalters@0: 	st[34] = st[33] + st[35] ;                           /* b84 = a71 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 7 */
tomwalters@0: 	st[36] = -st[34] - an_1 ;                            /* b73 = a102 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 13 */
tomwalters@0: 	st[15] = -st[14] ;                                   /* b133 = a122 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 12 */
tomwalters@0: 	st[16] = st[15] + st[18] ;                           /* a122-a121 */
tomwalters@0: 	st[17] = -ms->gamma121 * st[16] ;                    /* b120 */
tomwalters@0: 	st[19] = st[17] + st[15] ;                           /* b123 = a112 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 14 */
tomwalters@0: 	st[21] = -ms->gamma141 * st[20] ;                    /* b140 */
tomwalters@0: 	st[22] = st[21] + st[20] ;                           /* b143 = a111 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 15 */
tomwalters@0: 	st[23] = -st[24] ;                                   /* b153 = a113 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 11 */
tomwalters@0: 	st[25] = -st[22] - st[19] - st[23] ;                 /* b114 = a101 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 10 */
tomwalters@0: 	st[26] = st[25] - st[36] ;                           /* -a102+a101 */
tomwalters@0: 	st[27] = ms->gamma101 * st[26] ;                     /* b100 */
tomwalters@0: 	st[28] = st[27] + st[36] ;                           /* b103 = a161 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 17 */
tomwalters@0: 	st[10] = st[12] - st[11] ;                           /* b174 = a162 */
tomwalters@0: 	
tomwalters@0: 	/* Adaptor 16 */
tomwalters@0: 	st[13] = -st[28] - st[10] ;                          /* b163 = a181 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 19 */
tomwalters@0: 	st[5] = -st[6] ;                                     /* b193 = a182 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 18 */
tomwalters@0: 	st[7] = st[13] - st[5] ;                             /* -a182+a181 */
tomwalters@0: 	st[8] = ms->gamma181 * st[7] ;                       /* b180 */
tomwalters@0: 	st[9] = st[8] + st[5] ;                              /* b183 = a203 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 20 */
tomwalters@0: 	st[3] = -st[9] ;                                     /* b204 = a212 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 21 */
tomwalters@0: 	st[0] = st[2] - st[1] - st[3] ;                      /* b214 = a222 */
tomwalters@0: 
tomwalters@0: 
tomwalters@0:      /*** computation in BACKWARD direction for current input point ***/
tomwalters@0: 
tomwalters@0:      /* middle ear */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 22 */
tomwalters@0: 	a0 = st[4] + st[0] ;                                  /* a220 */
tomwalters@0: 	st[4] = st[4] - ms->gamma221 * a0 ;                   /* b221 */
tomwalters@0: 	b2 = -st[4] - a0 ;                                    /* b222 = a214 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 21 */
tomwalters@0: 	a0 = b2 - st[0] ;                                     /* a210 */
tomwalters@0: 	st[1] = st[1] - ms->gamma211 * a0 ;                   /* b211 */
tomwalters@0: 	b212 = st[3] - ms->gamma212 * a0 ;                    /* b212 = a204 */
tomwalters@0: 	st[2] = -st[1] - b212 - b2 ;                          /* b213 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 20 */
tomwalters@0: 	a0 = b212 - st[3] ;                                   /* a200 */
tomwalters@0: 	b202 = - ms->gamma202 * a0 ;                          /* b202 */
tomwalters@0: 	b3 = ms->gamma201 * a0 - b202 - b212 ;                /* b203 = a183 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 18 */
tomwalters@0: 	b2 = st[8] + b3 ;                                     /* b182 = a193 */
tomwalters@0: 	b1 = b2 - st[7] ;                                     /* b181 = a163 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 19 */
tomwalters@0: 	st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ;       /* b191 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 16 */
tomwalters@0: 	b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ;      /* b161 = a103 */
tomwalters@0: 	b2 = -b161 - b1 ;                                     /* b162 = a174 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 17 */
tomwalters@0: 	a0 = b2 - st[10] ;                                    /* a170 */
tomwalters@0: 	st[11] = st[11] - ms->gamma171 * a0 ;                 /* b171 */
tomwalters@0: 	st[12] = ms->gamma172 * a0 - b2 - st[11] ;            /* b173 */
tomwalters@0:  
tomwalters@0: 	/* Adaptor 10 */
tomwalters@0: 	st[37] = st[27] + b161 ;                              /* b102 = a73 */
tomwalters@0: 	st[38] = st[37] - st[26] ;                            /* b101 = a114 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 11 */
tomwalters@0: 	a0 = st[38] - st[25] ;                                /* a110 */
tomwalters@0: 	b1 = st[22] - ms->gamma111 * a0 ;                     /* b111 = a143 */
tomwalters@0: 	b2 = st[19] - ms->gamma112 * a0 ;                     /* b112 = a123 */
tomwalters@0: 	b3 = -b1 -b2 - st[38] ;                               /* b113 = a153 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 15 */
tomwalters@0: 	st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ;    /* b151 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 14 */
tomwalters@0: 	st[20] = st[21] + b1 ;                                /* b142 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 12 */
tomwalters@0: 	b2 = st[17] + b2 ;                                    /* b122 = a133 */
tomwalters@0: 	st[18] = b2 + st[16] ;                                /* b121 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 13 */
tomwalters@0: 	st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ;    /* b131 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 7 */
tomwalters@0: 	b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ;     /* b71 = a84 */
tomwalters@0: 	st[39] = -b1 - st[37] ;                               /* b72 = a63 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 8 */
tomwalters@0: 	st[35] = st[33] + b1 ;                                /* b83 */
tomwalters@0: 	b1 = st[35] + st[32] ;                                /* b81 = a94 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 9 */
tomwalters@0: 	a0 = b1 - st[31] ;                                    /* a90 */
tomwalters@0: 	st[29] = st[29] - ms->gamma91 * a0 ;                  /* b91 */
tomwalters@0: 	st[30] = -st[29] + ms->gamma92 * a0 - b1 ;            /* b93 */
tomwalters@0: 
tomwalters@0:      bn = st[39] ;
tomwalters@0: 
tomwalters@0:      /* concha and auditory canal */
tomwalters@0: 
tomwalters@0: 	while( tube_segment-- > 0 ) {
tomwalters@0: 
tomwalters@0: 	   /* get states */
tomwalters@0: 	   ts = *--eartube_states ;
tomwalters@0: 	   st = ( StateType * ) ts->states ;
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 6 */
tomwalters@0: 	   b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ;        /* b61 = a43 */
tomwalters@0: 	   st[9] = -b1 - bn ;                                 /* b62 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 4 */
tomwalters@0: 	   b2 = st[6] + b1 ;                                  /* b42 = a53 */
tomwalters@0: 	   b1 = b2 - st[5] ;                                  /* b41 = a33 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 5 */
tomwalters@0: 	   st[1] = st[0] + b2 ;                               /* b52 */
tomwalters@0: 
tomwalters@0: 	   /* Adaptor 3 */
tomwalters@0: 	   st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ;     /* b31 */
tomwalters@0: 	   bn = -st[4] - b1 ;                                 /* b32 = a63 */
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:      /* freefield - external ear */
tomwalters@0: 
tomwalters@0: 	/* get states */
tomwalters@0: 	st = ( StateType * ) ws->states ;      
tomwalters@0: 
tomwalters@0: 	/* Adaptor 1 */
tomwalters@0: 	b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ;           /* b11 = -a03 */
tomwalters@0: 	b2 = -b1 - bn ;                                       /* b12 = a23 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 2 */
tomwalters@0: 	st[0] = st[1] + b2 + st[0] ;                          /* b61 */
tomwalters@0: 
tomwalters@0: 	/* Adaptor 0 */
tomwalters@0: 	st[3] = b1 - st[4] + st[3] ;                          /* -b01 + p */
tomwalters@0: 
tomwalters@0:      /*** output ***/
tomwalters@0: #ifdef _EARDRUM_
tomwalters@0: 	st = ( StateType * ) ms->states ;      
tomwalters@0: 	out = 0.5 * ( st[39] + an_1 ) * ms->out_scale ; 
tomwalters@0: #else
tomwalters@0: 	out = -0.5 * b202 * ms->out_scale ; 
tomwalters@0: #endif
tomwalters@0: 
tomwalters@0: #ifdef _OVERFLOW_
tomwalters@0:         if( out > _MaxOutput_ || out < _MinOutput_ )
tomwalters@0: 	   fprintf( stderr, "Overflow error in Outer/Middle ear output=%e\n", ( double ) out ) ;
tomwalters@0: #endif     
tomwalters@0: 	*inout_ptr++ = out ; 
tomwalters@0: 
tomwalters@0:     }
tomwalters@0:     return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /******************************* CloseEarWDF() *************************************/
tomwalters@0: 
tomwalters@0: void CloseEarWDF( wave_states, eartube_states, earmiddle_states, tube_segments )
tomwalters@0:   register WaveWDFstate       *wave_states ;
tomwalters@0:   register EartubeWDFstate    **eartube_states ;
tomwalters@0:   register EarmiddleWDFstate  *earmiddle_states ;
tomwalters@0:   register int                tube_segments ;
tomwalters@0: {
tomwalters@0:   Delete( earmiddle_states->states ) ;
tomwalters@0:   Delete( earmiddle_states ) ;
tomwalters@0: 
tomwalters@0:   for( ; --tube_segments < 0 ; ) {
tomwalters@0:      Delete( eartube_states[ tube_segments ]->states ) ;
tomwalters@0:      Delete( eartube_states[ tube_segments ] ) ;
tomwalters@0:   }
tomwalters@0:   Delete( eartube_states ) ; 
tomwalters@0: 
tomwalters@0:   Delete( wave_states->states ) ;
tomwalters@0:   Delete( wave_states ) ;
tomwalters@0: 
tomwalters@0:   return ;
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /************** low level functions **************/
tomwalters@0: 
tomwalters@0: static double den221 ;
tomwalters@0: static double Kst_old_ratio ; 
tomwalters@0: 
tomwalters@0: void update_earmiddle_init( rate, z )
tomwalters@0:   double rate, z ;
tomwalters@0: {
tomwalters@0:   Kst_old_ratio = Kst_ratio ; 
tomwalters@0:   den221 = rate * z / Kst_max ;
tomwalters@0:   return ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: void update_earmiddle( ms )
tomwalters@0:   EarmiddleWDFstate *ms ;
tomwalters@0: {
tomwalters@0:   if( Kst_ratio < Kst_min_ratio )
tomwalters@0:     Kst_ratio = Kst_min_ratio ;
tomwalters@0: 
tomwalters@0:   ms->gamma221 = 2. * Kst_ratio / ( Kst_ratio + den221 ) ;
tomwalters@0:   ms->states[4] = ms->states[4] * Kst_ratio / Kst_old_ratio ;
tomwalters@0: 
tomwalters@0:   Kst_old_ratio = Kst_ratio ;
tomwalters@0:  
tomwalters@0:   return ;
tomwalters@0: }
tomwalters@0: