Mercurial > hg > aim92

comparison wdf/wdf_ear.c @ 0:5242703e91d3 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
author tomwalters
date Fri, 20 May 2011 15:19:45 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5242703e91d3
1 /*
2 Copyright (c) Applied Psychology Unit, Medical Research Council. 1988, 1989
3 ===========================================================================
4
5 Permission to use, copy, modify, and distribute this software without fee
6 is hereby granted for research purposes, provided that this copyright
7 notice appears in all copies and in all supporting documentation, and that
8 the software is not redistributed for any fee (except for a nominal shipping
9 charge). Anyone wanting to incorporate all or part of this software in a
10 commercial product must obtain a license from the Medical Research Council.
11
12 The MRC makes no representations about the suitability of this
13 software for any purpose. It is provided "as is" without express or implied
14 warranty.
15
16 THE MRC DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING
17 ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE
18 A.P.U. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY
19 DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN
20 AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
21 OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
22 */
23
24 /*
25 Acknowledgment:
26 ==============
27
28 The source code provided in this file was originally developed by
29 Christian Giguere as part of a Ph.D degree at the Department of
30 Engineering of the University of Cambridge from April 1990 to
31 November 1993. The code was subsequently adapted under a grant
32 from the Hearing Research Trust for full compatibility with
33 AIM Release 6.15.
34
35 Christian Giguere 25/03/94
36
37 */
38
39 /*
40 ===========================================================
41 wdf_ear.c
42 ===========================================================
43
44 Wave digital filter (WDF) implementation of the outer and
45 middle ear (EAR) electroacoustic networks.
46
47 Author : Christian Giguere
48 First written : 01st June, 1991
49 Last edited : 18th February, 1994
50
51 References:
52 (1) C.Giguere and P.C.Woodland (1992). Technical Report
53 CUED/F-INFENG/TR.93, University of Cambridge.
54 (2) C.Giguere and P.C.Woodland (1994). JASA 95(1): 331-342.
55 (3) M.E.Lutman and A.M.Martin (1979). J.Sound.Vib. 64(1): 133-157
56 (4) A.Fettweis (1986). Proceedings IEEE 74(2): 270-327.
57
58 Note: Ref (1) describes an implementation where the outer
59 ear and middle ear module is _UNCOUPLED_ to the cochlea.
60 In this case, the middle ear network is that of Lutman
61 and Martin (1979) unmodified.
62
63 Ref (2) describes an implementation where the outer
64 ear and middle ear module is _COUPLED_ to the cochlea.
65 In this case, the middle ear network is a modified
66 version of Lutman and Martin (1979) network (See Fig.2)
67 ===========================================================
68 */
69
70
71 /***** includes *****/
72
73 #include <math.h>
74 #include <stdio.h>
75 #include "stitch.h"
76 #include "calc.h"
77 #include "calc_tl.h"
78 #include "wdf_ear.h"
79
80 /***** defines *****/
81
82 #define N_OF_FREEFIELD_STATES ( 12 )
83 #define N_OF_EARTUBE_STATES ( 16 ) /* per segment */
84 #define N_OF_EARMIDDLE_STATES ( 48 )
85
86 #if 0
87 #define _DEBUG_
88 #define _OVERFLOW_
89 #endif
90
91 #if 0
92 #define _IMPULSE_ /* bypasses input wave with a delta impulse */
93 #endif
94
95 #if 0
96 #define _EARDRUM_ /* outputs eardrum pressure instead of stapes velocity */
97 #endif
98
99 /* Middle ear circuit elements in CGS units (From Lutman and Martin (1979) */
100
101 #define La ( 14.0e-03 ) /* Middle ear cavities */
102 #define Cp ( 5.1e-06 )
103 #define Ra ( 10.0 )
104 #define Rm ( 390.0 )
105 #define Ct ( 0.35e-06 )
106 #define Rd1 ( 200.0 ) /* Eardrum losses */
107 #define Cd1 ( 0.8e-06 )
108 #define Cd2 ( 0.4e-06 )
109 #define Rd2 ( 220.0 )
110 #define Ld ( 15.0e-03 )
111 #define Rd3 ( 5900.0 )
112 #define Cd3 ( 0.2e-06 )
113 #define Lo ( 40.0e-03 ) /* Eardrum, mallus, incus */
114 #define Co ( 1.4e-06 )
115 #define Ro ( 70.0 )
116 #define Cs ( 0.25e-06 ) /* Incudo-stapedial joint */
117 #define Rs ( 3000.0 )
118 #define Lc ( 45.0e-03 ) /* Stapes and cochlea */
119 #define Cc ( 0.6e-06 )
120 #define Rc ( 550.0 )
121
122 /* Additional middle ear circuit elements in CGS units (From Giguere and Woodland (1994) */
123
124 #define Ral ( 100.0 ) /* annular ligaments */
125 #define Kst_max ( 1./0.1e-06 ) /* max. Stapedius stiffness */
126 #define Kst_min_ratio ( 0.0001 ) /* min value of Kst_ratio */
127
128
129 /***** external variables *****/
130
131 double Kst_ratio = Kst_min_ratio ; /* initial fraction of Kst_max */
132 /* updated by feedback module */
133 static double rn_1 ;
134
135
136 /***** functions *****/
137
138 /******************************* WDF set-up functions ************************************
139 * name: function:
140 *
141 * FreefieldWDF() Set-up function for the WDF implementation of the external ear network
142 * (Giguere and Woodland, 1994, Figs. 1 and 8). It generates the multiplier
143 * coefficients and initializes the states of the filter.
144 *
145 * In the case of the _UNCOUPLED_ ear version, it is called by the
146 * function ``Ear()'' in file ``ear.c''.
147 * In the case of the _COUPLED_ ear version, it is called by the
148 * function ``TLF_GenBank()'' in file ``bank_tl.c''.
149 *
150 * It returns a pointer to a structure containing all the relevant
151 * information for the computation of the whole filter.
152 *
153 * EartubeWDF() Set-up function for the WDF implementation of the concha and auditory
154 * canal transmission lines (Giguere and Woodland, 1994, Figs. 1 and 8).
155 * It generates the multiplier coefficients and initializes the states
156 * of the filter.
157 *
158 * In the case of the _UNCOUPLED_ ear version, it is called by the
159 * function ``Ear()'' in file ``ear.c'' once for each segment of the
160 * transmission lines starting from the outermost segment.
161 * In the case of the _COUPLED_ ear version, it is called by the
162 * function ``TLF_GenBank()'' in file ``bank_tl.c'' once for each segment
163 * of the transmission lines starting from the outermost segment.
164 *
165 * It returns a pointer to a structure containing all the relevant
166 * information for the current filter segment.
167 *
168 * EarmiddleWDF() Set-up function for the WDF implementation of the middle ear network
169 * (Giguere and Woodland, 1994, Figs. 2 and 9). It generates the multiplier
170 * coefficients and initializes the states of the filter.
171 *
172 * In the case of the _UNCOUPLED_ ear version, it is called by the
173 * function ``Ear()'' in file ``ear.c''.
174 * In the case of the _COUPLED_ ear version, it is called by the
175 * function ``TLF_GenBank()'' in file ``bank_tl.c''.
176 *
177 * It returns a pointer to a structure containing all the relevant
178 * information for the computation of the whole filter.
179 ******************************************************************************************/
180
181 /************************** FreefieldWDF() *****************************/
182
183 WaveWDFstate *FreefieldWDF( samplerate, rho, velocity, diffraction_radius, radiation_radius )
184 double samplerate ;
185 double rho, velocity ;
186 double diffraction_radius, radiation_radius ;
187 {
188 DeclareNew( WaveWDFstate *, filter_state ) ;
189 double fs2 ;
190 double r1, r2, r3, g1, g2, g3 ;
191 double zh, zr ;
192 double Lh, Rh, Lr, Rr ;
193
194 /*** electrical circuit elements (CGS) ***/
195 Rh = rho * velocity / ( Pi * diffraction_radius * diffraction_radius ) ;
196 Lh = 0.5 * rho / ( Pi * diffraction_radius ) ;
197 Rr = rho * velocity / ( Pi * radiation_radius * radiation_radius ) ;
198 Lr = 0.7 * rho / ( Pi * radiation_radius ) ;
199
200 /*** sampling frequency ***/
201 fs2 = samplerate * 2. ;
202
203 /*** WDF-port resistances and multiplier coefficients ***/
204
205 /* Adaptor 0 */
206 g1 = 1. / ( fs2 * Lh ) ; /* (1 / Zlh) */
207 g2 = 1. / Rh ; /* (1 / Zrh) */
208 g3 = g1 + g2 ;
209 zh = 1. / g3 ;
210 filter_state->gamma01 = g1 / g3 ;
211
212 /* Adaptor 2 */
213 g1 = 1. / ( fs2 * Lr ) ; /* (1 / Zlr) */
214 g2 = 1. / Rr ; /* (1 / Zrr) */
215 g3 = g1 + g2 ;
216 zr = 1. / g3 ;
217 filter_state->gamma21 = g1 / g3 ;
218
219 /* Adaptor 1 */
220 r1 = zh ;
221 r2 = zr ;
222 r3 = rn_1 = r1 + r2 ;
223 filter_state->gamma11 = r1 / r3 ;
224
225 #ifdef _DEBUG_
226 fprintf( stderr, "Rh=%e Lh=%e Rr=%e Lr=%e\n", Rh, Lh, Rr, Lr ) ;
227 fprintf( stderr, "gamma01=%f\n", filter_state->gamma01 ) ;
228 fprintf( stderr, "gamma11=%f\n", filter_state->gamma11 ) ;
229 fprintf( stderr, "gamma21=%f\n", filter_state->gamma21 ) ;
230 #endif
231
232 /*** initialise filter states ***/
233 filter_state->Nstates = N_OF_FREEFIELD_STATES ;
234 filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
235 "wdf_ear.c for freefield states\n" ) ;
236
237 /*** return ***/
238 return( filter_state ) ;
239 }
240
241 /************************** EartubeWDF() *****************************/
242
243 EartubeWDFstate *EartubeWDF( samplerate, rho, velocity, diam, seglength, attn )
244 double samplerate ;
245 double rho, velocity ;
246 double diam, seglength ;
247 double attn ;
248 {
249 DeclareNew( EartubeWDFstate *, filter_state ) ;
250 double an, d, c, k, alpha, fs2 ;
251 double Ln, Cn, Gn, Zn ;
252 double r1, r2, r3, g1, g2, g3, r33, g43, g53 ;
253
254 /*** acoustic tube parameters for current segment ***/
255 an = Pi * diam * diam / 4. ; /* cross-sectional area of tube (cm2) */
256 d = seglength; /* length of tube (cm) */
257 c = velocity ; /* sound velocity in air (cm/s) */
258 alpha = attn ; /* attenuation constant of travelling waves (1/cm) */
259
260 /*** sampling frequency ***/
261 fs2 = samplerate * 2. ;
262
263 /*** electrical circuit elements (CGS units) for current segment ***/
264 Ln = ( rho * d ) / an ; /* series inductor (either Lch or Lcl) */
265 Cn = ( an * d ) / ( rho * c * c ) ; /* shunt capacitor (either Cch or Ccl) */
266 Zn = sqrt( Ln / Cn ) ; /* characteristic impedance ( either Zch or Zcl) */
267 Gn = 2. / Zn * alpha * d ; /* shunt conductance (either Gch or Gcl) */
268
269 /*** WDF multiplier coefficients for current tube segment ***/
270
271 /* Adaptor 3 */
272 r1 = fs2 * Ln / 2. ; /* 0.5 Zl */
273 r2 = rn_1 ;
274 r3 = r33 = r1 + r2 ;
275 filter_state->gamma31 = r1 / r3 ;
276
277 /* Adaptor 5 */
278 g1 = Gn ; /* (1 / Zg ) */
279 g2 = fs2 * Cn ; /* (1 / Zc ) */
280 g3 = g53 = g1 + g2 ;
281 filter_state->gamma51 = g1 / g3 ;
282
283 /* Adaptor 4 */
284 g1 = 1 / r33 ;
285 g2 = g53 ;
286 g3 = g43 = g1 + g2 ;
287 filter_state->gamma41 = g1 / g3 ;
288
289 /* Adaptor 6 */
290 r1 = 1 / g43 ;
291 r2 = fs2 * Ln / 2. ; /* 0.5 Zl */
292 r3 = rn_1 = r1 + r2 ;
293 filter_state->gamma61 = r1 / r3 ;
294
295 #ifdef _DEBUG_
296 fprintf( stderr, "gamma31=%f ", filter_state->gamma31 ) ;
297 fprintf( stderr, "gamma41=%f ", filter_state->gamma41 ) ;
298 fprintf( stderr, "gamma51=%f ", filter_state->gamma51 ) ;
299 fprintf( stderr, "gamma61=%f\n", filter_state->gamma61 ) ;
300 #endif
301
302 /*** initialise filter states ***/
303 filter_state->Nstates = N_OF_EARTUBE_STATES ;
304 filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
305 "wdf_ear.c for earcanal states\n" ) ;
306
307 /*** return ***/
308 return( filter_state ) ;
309 }
310
311
312 /********************** EarmiddleWDF() ***************************/
313
314 EarmiddleWDFstate *EarmiddleWDF( samplerate, zov, output_scale, ratio )
315 double samplerate, zov, output_scale, ratio ;
316 {
317 DeclareNew( EarmiddleWDFstate *, filter_state ) ;
318 double fs2 ;
319 double r1, r2, r3, r4, g1, g2, g3, g4 ;
320 double r73, r84, r94, r103, r114, r123, r133, r143, r153 ;
321 double r163, r174, r183, r193, r204, r214 ;
322 void update_earmiddle_init(), update_earmiddle() ;
323
324 /*** sampling frequency ***/
325 fs2 = samplerate * 2. ;
326
327 /*** WDF-port resistances and multiplier coefficients ***/
328
329 /* Adaptor 9 */
330 r1 = 1. / ( fs2 * Cp ) ; /* Zcp */
331 r2 = Ra ; /* Zra */
332 r3 = fs2 * La ; /* Zla */
333 r4 = r94 = r1 + r2 + r3 ;
334 filter_state->gamma91 = r1 / r4 ;
335 filter_state->gamma92 = r2 / r4 ;
336
337 /* Adaptor 8 */
338 g1 = 1. / r94 ;
339 g2 = 1. / Rm ; /* (1 / Zrm) */
340 g3 = fs2 * Ct ; /* (1 / Zct) */
341 g4 = g1 + g2 + g3 ;
342 r84 = 1. /g4 ;
343 filter_state->gamma81 = g1 / g4 ;
344 filter_state->gamma82 = g2 / g4 ;
345
346 /* Adaptor 7 */
347 r1 = r84 ;
348 r2 = rn_1 ;
349 r3 = r73 = r1 + r2 ;
350 filter_state->gamma71 = r1 / r3 ;
351
352 /* Adaptor 13 */
353 r1 = 1. / ( fs2 * Cd2 ) ; /* Zcd2 */
354 r2 = Rd2 ; /* Zrd2 */
355 r3 = r133 = r1 + r2 ;
356 filter_state->gamma131 = r1 / r3 ;
357
358 /* Adaptor 12 */
359 g1 = 1. / ( fs2 * Ld ) ; /* (1 / Zld) */
360 g2 = 1. / r133 ;
361 g3 = g1 + g2 ;
362 r123 = 1. / g3 ;
363 filter_state->gamma121 = g1 / g3 ;
364
365 /* Adaptor 14 */
366 g1 = 1. / Rd3 ; /* (1 / Zrd3) */
367 g2 = fs2 * Cd3 ; /* (1 / Zcd3) */
368 g3 = g1 + g2 ;
369 r143 = 1. / g3 ;
370 filter_state->gamma141 = g1 / g3 ;
371
372 /* Adaptor 15 */
373 r1 = 1. / ( fs2 * Cd1 ) ; /* Zcd1 */
374 r2 = Rd1 ; /* Zrd1 */
375 r3 = r153 = r1 + r2 ;
376 filter_state->gamma151 = r1 / r3 ;
377
378 /* Adaptor 11 */
379 r1 = r143 ;
380 r2 = r123 ;
381 r3 = r153 ;
382 r4 = r114 = r1 + r2 + r3 ;
383 filter_state->gamma111 = r1 / r4 ;
384 filter_state->gamma112 = r2 / r4 ;
385
386 /* Adaptor 10 */
387 g1 = 1. / r114 ;
388 g2 = 1. / r73 ;
389 g3 = g1 + g2 ;
390 r103 = 1. /g3 ;
391 filter_state->gamma101 = g1 / g3 ;
392
393 /* Adaptor 17 */
394 r1 = 1. / ( fs2 * Co ) ; /* Zco */
395 r2 = Ro ; /* Zro */
396 r3 = fs2 * Lo ; /* Zlo */
397 r4 = r174 = r1 + r2 + r3 ;
398 filter_state->gamma171 = r1 / r4 ;
399 filter_state->gamma172 = r2 / r4 ;
400
401 /* Adaptor 16 */
402 r1 = r103 ;
403 r2 = r174 ;
404 r3 = r163 = r1 + r2 ;
405 filter_state->gamma161 = r1 / r3 ;
406
407 /* Adaptor 19 */
408 r1 = 1. / ( fs2 * Cs ) ; /* Zcs */
409 r2 = Rs ; /* Zrs */
410 r3 = r193 = r1 + r2 ;
411 filter_state->gamma191 = r1 / r3 ;
412
413 /* Adaptor 18 */
414 g1 = 1. / r163 ;
415 g2 = 1. / r193 ;
416 g3 = g1 + g2 ;
417 r183 = 1. / g3 ;
418 filter_state->gamma181 = g1 / g3 ;
419
420 /* Adaptor 20 */
421 #ifdef _DEBUG_
422 fprintf( stderr, "Zov=%e\n", zov ) ;
423 #endif
424 if( zov != 0. ) { /* coupled version ? */
425 r1 = Ral ; /* Zral */
426 r2 = zov / ( ratio * ratio ) ; /* Zov / (r*r) */
427 }
428 else { /* uncoupled version */
429 r1 = 0.0 ;
430 r2 = Rc ; /* Zrc */
431 }
432 r3 = r183 ;
433 r4 = r204 = r1 + r2 +r3 ;
434 filter_state->gamma201 = r1 / r4 ;
435 filter_state->gamma202 = r2 / r4 ;
436 filter_state->ratio = ratio ;
437
438 /* Adaptor 21 */
439 r1 = 1. / ( fs2 * Cc ) ; /* Zcc */
440 r2 = r204 ;
441 if( zov != 0. ) /* coupled version ? */
442 r3 = 0.0 ;
443 else /* uncoupled version */
444 r3 = fs2 * Lc ; /* Zlc */
445 r4 = r214 = r1 + r2 + r3 ;
446 filter_state->gamma211 = r1 / r4 ;
447 filter_state->gamma212 = r2 / r4 ;
448
449 /* Adaptor 22 */
450 r1 = Kst_ratio * Kst_max / fs2 ; /* Zcst (initial value only) */
451 r2 = r214 ;
452 r3 = r1 + r2 ;
453 filter_state->gamma221 = 2. * r1 / r3 ;
454
455 #ifdef _DEBUG_
456 fprintf ( stderr, "gamma71 =%f \n", filter_state->gamma71 ) ;
457 fprintf ( stderr, "gamma81 =%f gamma82 =%f\n", filter_state->gamma81
458 , filter_state->gamma82 ) ;
459 fprintf ( stderr, "gamma91 =%f gamma92 =%f\n", filter_state->gamma91
460 , filter_state->gamma92 ) ;
461 fprintf ( stderr, "gamma101=%f \n", filter_state->gamma101 ) ;
462 fprintf ( stderr, "gamma111=%f gamma112=%f\n", filter_state->gamma111
463 , filter_state->gamma112 ) ;
464 fprintf ( stderr, "gamma121=%f \n", filter_state->gamma121 ) ;
465 fprintf ( stderr, "gamma131=%f \n", filter_state->gamma131 ) ;
466 fprintf ( stderr, "gamma141=%f \n", filter_state->gamma141 ) ;
467 fprintf ( stderr, "gamma151=%f \n", filter_state->gamma151 ) ;
468 fprintf ( stderr, "gamma161=%f \n", filter_state->gamma161 ) ;
469 fprintf ( stderr, "gamma171=%f gamma172=%f\n", filter_state->gamma171
470 , filter_state->gamma172 ) ;
471 fprintf ( stderr, "gamma181=%f \n", filter_state->gamma181 ) ;
472 fprintf ( stderr, "gamma191=%f \n", filter_state->gamma191 ) ;
473 fprintf ( stderr, "gamma201=%f gamma202=%f\n", filter_state->gamma201
474 , filter_state->gamma202 ) ;
475 fprintf ( stderr, "gamma211=%f gamma212=%f\n", filter_state->gamma211
476 , filter_state->gamma212 ) ;
477 fprintf ( stderr, "gamma221=%f \n", filter_state->gamma221 ) ;
478 #endif
479
480 /*** specify and initialise coefficient update procedure ***/
481 filter_state->update_proc = update_earmiddle ;
482 ( void ) update_earmiddle_init( fs2, r214 ) ;
483
484 /*** output scaling ***/
485 filter_state->out_scale = output_scale ;
486
487 /*** initialise filter states ***/
488 filter_state->Nstates = N_OF_EARMIDDLE_STATES ;
489 filter_state->states = ( char * ) NewZeroedArray( StateType, filter_state->Nstates,
490 "wdf_ear.c for middle ear states\n" ) ;
491
492 /*** return ***/
493 return( filter_state ) ;
494 }
495
496
497 /*********************************** WDF procedures ***************************************
498 * name: function:
499 *
500 * DoEarWDF() WDF realization of the _UNCOUPLED_ version of the outer ear and
501 * middle ear networks. It is called by function ``ear_callback()''
502 * in file ``ear.c'' once for a specified number of input points.
503 * It computes the stapes velocity for each input point and keeps
504 * track of the filter states.
505 *
506 * Note: The code for the realization of the _COUPLED_ ear version
507 * is located in file ``wdf_tl.c''.
508 *
509 *
510 * CloseEarWDF() Deletes all private data structures and arrays of the outer and
511 * middle ear filter upon completion of filtering.
512 * In the case of the _UNCOUPLED_ ear version, it is called by the
513 * function ``ear_callback()'' in file ``ear.c''.
514 * In the case of the _COUPLED_ ear version, it is called by the
515 * function ``tlf_bank_callbank()'' in file ``bank_tl.c''.
516 ******************************************************************************************/
517
518 /******************************* DoEarWDF() **************************************/
519
520 void DoEarWDF( ws, eartube_states, ms, tube_segments, inout_ptr, points )
521 register WaveWDFstate *ws ;
522 register EartubeWDFstate **eartube_states ;
523 register EarmiddleWDFstate *ms ;
524 register int tube_segments ;
525 register DataType *inout_ptr ;
526 register int points ;
527 {
528 register EartubeWDFstate *ts ;
529 register StateType *st, in, out ;
530 register StateType a0, an_1, b1, b2, b3, b161, b202, b212, bn ;
531 register int tube_segment = -1 ;
532 #ifdef _IMPULSE_
533 static int impulse = 10000 ;
534 #endif
535
536
537 /***** update stapedial muscle state *****/
538
539 ( void ) ms->update_proc( ms ) ;
540
541 /***** for all points *****/
542
543 while( points-- > 0 ) {
544
545 /*** input point ***/
546
547 in = ( StateType ) *inout_ptr ;
548
549 #ifdef _IMPULSE_
550 if( impulse != 0 ) {
551 in = ( StateType ) impulse ;
552 impulse = 0 ;
553 }
554 else
555 in = ( StateType ) 0 ;
556 #endif
557
558 /*** computation in FORWARD direction for current input point ***/
559
560 /* freefield - external ear */
561
562 /* get states */
563 st = ( StateType * ) ws->states ;
564
565 /* Adaptor 0 */
566 st[4] = ws->gamma01 * ( st[3] - in ) ; /* b00 */
567 st[5] = st[4] + in + in ; /* b03 = -a11 */
568
569 /* Adaptor 2 */
570 st[1] = -ws->gamma21 * st[0] ; /* b20 */
571 st[2] = st[1] ; /* b23 = a12 */
572
573 /* Adaptor 1 */
574 st[6] = st[5] - st[2] ; /* b13 = a32 */
575 an_1 = st[6] ;
576
577 /* concha and auditory canal */
578
579 while( ++tube_segment < tube_segments ) {
580
581 /* get states */
582 ts = *eartube_states++ ;
583 st = ( StateType * ) ts->states ;
584
585 /* Adaptor 3 */
586 st[3] = st[4] - an_1 ; /* b33 = a41 */
587
588 /* Adaptor 5 */
589 st[0] = -ts->gamma51 * st[1] ; /* b50 */
590 st[2] = st[0] + st[1] ; /* b53 = a42 */
591
592 /* Adaptor 4 */
593 st[5] = st[3] - st[2] ; /* -a42+a41 */
594 st[6] = ts->gamma41 * st[5] ; /* b40 */
595 st[7] = st[6] + st[2] ; /* b43 = a61 */
596
597 /* Adaptor 6 */
598 st[8] = st[9] - st[7] ; /* b63 = a32 */
599 an_1 = st[8] ;
600 }
601
602 /* middle ear */
603
604 /* get states */
605 st = ( StateType * ) ms->states ;
606
607 /* Adaptor 9 */
608 st[31] = st[30] - st[29] ; /* b94 = a81 */
609
610 /* Adaptor 8 */
611 st[32] = st[35] - st[31] ; /* a83-a81 */
612 st[33] = -ms->gamma81 * st[32] - ms->gamma82 * st[35] ; /* b80 */
613 st[34] = st[33] + st[35] ; /* b84 = a71 */
614
615 /* Adaptor 7 */
616 st[36] = -st[34] - an_1 ; /* b73 = a102 */
617
618 /* Adaptor 13 */
619 st[15] = -st[14] ; /* b133 = a122 */
620
621 /* Adaptor 12 */
622 st[16] = st[15] + st[18] ; /* a122-a121 */
623 st[17] = -ms->gamma121 * st[16] ; /* b120 */
624 st[19] = st[17] + st[15] ; /* b123 = a112 */
625
626 /* Adaptor 14 */
627 st[21] = -ms->gamma141 * st[20] ; /* b140 */
628 st[22] = st[21] + st[20] ; /* b143 = a111 */
629
630 /* Adaptor 15 */
631 st[23] = -st[24] ; /* b153 = a113 */
632
633 /* Adaptor 11 */
634 st[25] = -st[22] - st[19] - st[23] ; /* b114 = a101 */
635
636 /* Adaptor 10 */
637 st[26] = st[25] - st[36] ; /* -a102+a101 */
638 st[27] = ms->gamma101 * st[26] ; /* b100 */
639 st[28] = st[27] + st[36] ; /* b103 = a161 */
640
641 /* Adaptor 17 */
642 st[10] = st[12] - st[11] ; /* b174 = a162 */
643
644 /* Adaptor 16 */
645 st[13] = -st[28] - st[10] ; /* b163 = a181 */
646
647 /* Adaptor 19 */
648 st[5] = -st[6] ; /* b193 = a182 */
649
650 /* Adaptor 18 */
651 st[7] = st[13] - st[5] ; /* -a182+a181 */
652 st[8] = ms->gamma181 * st[7] ; /* b180 */
653 st[9] = st[8] + st[5] ; /* b183 = a203 */
654
655 /* Adaptor 20 */
656 st[3] = -st[9] ; /* b204 = a212 */
657
658 /* Adaptor 21 */
659 st[0] = st[2] - st[1] - st[3] ; /* b214 = a222 */
660
661
662 /*** computation in BACKWARD direction for current input point ***/
663
664 /* middle ear */
665
666 /* Adaptor 22 */
667 a0 = st[4] + st[0] ; /* a220 */
668 st[4] = st[4] - ms->gamma221 * a0 ; /* b221 */
669 b2 = -st[4] - a0 ; /* b222 = a214 */
670
671 /* Adaptor 21 */
672 a0 = b2 - st[0] ; /* a210 */
673 st[1] = st[1] - ms->gamma211 * a0 ; /* b211 */
674 b212 = st[3] - ms->gamma212 * a0 ; /* b212 = a204 */
675 st[2] = -st[1] - b212 - b2 ; /* b213 */
676
677 /* Adaptor 20 */
678 a0 = b212 - st[3] ; /* a200 */
679 b202 = - ms->gamma202 * a0 ; /* b202 */
680 b3 = ms->gamma201 * a0 - b202 - b212 ; /* b203 = a183 */
681
682 /* Adaptor 18 */
683 b2 = st[8] + b3 ; /* b182 = a193 */
684 b1 = b2 - st[7] ; /* b181 = a163 */
685
686 /* Adaptor 19 */
687 st[6] = st[6] + ms->gamma191 * ( st[5] - b2 ) ; /* b191 */
688
689 /* Adaptor 16 */
690 b161 = st[28] + ms->gamma161 * ( st[13] - b1 ) ; /* b161 = a103 */
691 b2 = -b161 - b1 ; /* b162 = a174 */
692
693 /* Adaptor 17 */
694 a0 = b2 - st[10] ; /* a170 */
695 st[11] = st[11] - ms->gamma171 * a0 ; /* b171 */
696 st[12] = ms->gamma172 * a0 - b2 - st[11] ; /* b173 */
697
698 /* Adaptor 10 */
699 st[37] = st[27] + b161 ; /* b102 = a73 */
700 st[38] = st[37] - st[26] ; /* b101 = a114 */
701
702 /* Adaptor 11 */
703 a0 = st[38] - st[25] ; /* a110 */
704 b1 = st[22] - ms->gamma111 * a0 ; /* b111 = a143 */
705 b2 = st[19] - ms->gamma112 * a0 ; /* b112 = a123 */
706 b3 = -b1 -b2 - st[38] ; /* b113 = a153 */
707
708 /* Adaptor 15 */
709 st[24] = st[24] + ms->gamma151 * ( st[23] - b3 ) ; /* b151 */
710
711 /* Adaptor 14 */
712 st[20] = st[21] + b1 ; /* b142 */
713
714 /* Adaptor 12 */
715 b2 = st[17] + b2 ; /* b122 = a133 */
716 st[18] = b2 + st[16] ; /* b121 */
717
718 /* Adaptor 13 */
719 st[14] = st[14] + ms->gamma131 * ( st[15] - b2 ) ; /* b131 */
720
721 /* Adaptor 7 */
722 b1 = st[34] + ms->gamma71 * ( st[36] - st[37] ) ; /* b71 = a84 */
723 st[39] = -b1 - st[37] ; /* b72 = a63 */
724
725 /* Adaptor 8 */
726 st[35] = st[33] + b1 ; /* b83 */
727 b1 = st[35] + st[32] ; /* b81 = a94 */
728
729 /* Adaptor 9 */
730 a0 = b1 - st[31] ; /* a90 */
731 st[29] = st[29] - ms->gamma91 * a0 ; /* b91 */
732 st[30] = -st[29] + ms->gamma92 * a0 - b1 ; /* b93 */
733
734 bn = st[39] ;
735
736 /* concha and auditory canal */
737
738 while( tube_segment-- > 0 ) {
739
740 /* get states */
741 ts = *--eartube_states ;
742 st = ( StateType * ) ts->states ;
743
744 /* Adaptor 6 */
745 b1 = st[7] - ts->gamma61 * ( bn - st[8] ) ; /* b61 = a43 */
746 st[9] = -b1 - bn ; /* b62 */
747
748 /* Adaptor 4 */
749 b2 = st[6] + b1 ; /* b42 = a53 */
750 b1 = b2 - st[5] ; /* b41 = a33 */
751
752 /* Adaptor 5 */
753 st[1] = st[0] + b2 ; /* b52 */
754
755 /* Adaptor 3 */
756 st[4] = ts->gamma31 * ( st[3] - b1 ) - st[4] ; /* b31 */
757 bn = -st[4] - b1 ; /* b32 = a63 */
758 }
759
760 /* freefield - external ear */
761
762 /* get states */
763 st = ( StateType * ) ws->states ;
764
765 /* Adaptor 1 */
766 b1 = ws->gamma11 * ( st[6] - bn ) - st[5] ; /* b11 = -a03 */
767 b2 = -b1 - bn ; /* b12 = a23 */
768
769 /* Adaptor 2 */
770 st[0] = st[1] + b2 + st[0] ; /* b61 */
771
772 /* Adaptor 0 */
773 st[3] = b1 - st[4] + st[3] ; /* -b01 + p */
774
775 /*** output ***/
776 #ifdef _EARDRUM_
777 st = ( StateType * ) ms->states ;
778 out = 0.5 * ( st[39] + an_1 ) * ms->out_scale ;
779 #else
780 out = -0.5 * b202 * ms->out_scale ;
781 #endif
782
783 #ifdef _OVERFLOW_
784 if( out > _MaxOutput_ || out < _MinOutput_ )
785 fprintf( stderr, "Overflow error in Outer/Middle ear output=%e\n", ( double ) out ) ;
786 #endif
787 *inout_ptr++ = out ;
788
789 }
790 return ;
791 }
792
793
794 /******************************* CloseEarWDF() *************************************/
795
796 void CloseEarWDF( wave_states, eartube_states, earmiddle_states, tube_segments )
797 register WaveWDFstate *wave_states ;
798 register EartubeWDFstate **eartube_states ;
799 register EarmiddleWDFstate *earmiddle_states ;
800 register int tube_segments ;
801 {
802 Delete( earmiddle_states->states ) ;
803 Delete( earmiddle_states ) ;
804
805 for( ; --tube_segments < 0 ; ) {
806 Delete( eartube_states[ tube_segments ]->states ) ;
807 Delete( eartube_states[ tube_segments ] ) ;
808 }
809 Delete( eartube_states ) ;
810
811 Delete( wave_states->states ) ;
812 Delete( wave_states ) ;
813
814 return ;
815
816 }
817
818 /************** low level functions **************/
819
820 static double den221 ;
821 static double Kst_old_ratio ;
822
823 void update_earmiddle_init( rate, z )
824 double rate, z ;
825 {
826 Kst_old_ratio = Kst_ratio ;
827 den221 = rate * z / Kst_max ;
828 return ;
829 }
830
831 void update_earmiddle( ms )
832 EarmiddleWDFstate *ms ;
833 {
834 if( Kst_ratio < Kst_min_ratio )
835 Kst_ratio = Kst_min_ratio ;
836
837 ms->gamma221 = 2. * Kst_ratio / ( Kst_ratio + den221 ) ;
838 ms->states[4] = ms->states[4] * Kst_ratio / Kst_old_ratio ;
839
840 Kst_old_ratio = Kst_ratio ;
841
842 return ;
843 }
844