Mercurial > hg > aim92
diff wdf/wdf_ear.c @ 0:5242703e91d3 tip
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
| author | tomwalters |
|---|---|
| date | Fri, 20 May 2011 15:19:45 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/wdf/wdf_ear.c Fri May 20 15:19:45 2011 +0100 @@ -0,0 +1,844 @@ +/* + 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 ; +} +